CoCalc Logo Icon
StoreFeaturesDocsShareSupportNewsAboutSign UpSign In

Real-time collaboration for Jupyter Notebooks, Linux Terminals, LaTeX, VS Code, R IDE, and more,
all in one place.

| Download

GAP 4.8.9 installation with standard packages -- copy to your CoCalc project to get it

Views: 418346
1
/* Functions needed for bootstrapping the gmp build, based on mini-gmp.
2
3
Copyright 2001, 2002, 2004, 2011, 2012 Free Software Foundation, Inc.
4
5
This file is part of the GNU MP Library.
6
7
The GNU MP Library is free software; you can redistribute it and/or modify
8
it under the terms of either:
9
10
* the GNU Lesser General Public License as published by the Free
11
Software Foundation; either version 3 of the License, or (at your
12
option) any later version.
13
14
or
15
16
* the GNU General Public License as published by the Free Software
17
Foundation; either version 2 of the License, or (at your option) any
18
later version.
19
20
or both in parallel, as here.
21
22
The GNU MP Library is distributed in the hope that it will be useful, but
23
WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
24
or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
25
for more details.
26
27
You should have received copies of the GNU General Public License and the
28
GNU Lesser General Public License along with the GNU MP Library. If not,
29
see https://www.gnu.org/licenses/. */
30
31
32
#include "mini-gmp/mini-gmp.c"
33
34
#define MIN(l,o) ((l) < (o) ? (l) : (o))
35
#define PTR(x) ((x)->_mp_d)
36
#define SIZ(x) ((x)->_mp_size)
37
38
#define xmalloc gmp_default_alloc
39
40
int
41
isprime (unsigned long int t)
42
{
43
unsigned long int q, r, d;
44
45
if (t < 32)
46
return (0xa08a28acUL >> t) & 1;
47
if ((t & 1) == 0)
48
return 0;
49
50
if (t % 3 == 0)
51
return 0;
52
if (t % 5 == 0)
53
return 0;
54
if (t % 7 == 0)
55
return 0;
56
57
for (d = 11;;)
58
{
59
q = t / d;
60
r = t - q * d;
61
if (q < d)
62
return 1;
63
if (r == 0)
64
break;
65
d += 2;
66
q = t / d;
67
r = t - q * d;
68
if (q < d)
69
return 1;
70
if (r == 0)
71
break;
72
d += 4;
73
}
74
return 0;
75
}
76
77
int
78
log2_ceil (int n)
79
{
80
int e;
81
assert (n >= 1);
82
for (e = 0; ; e++)
83
if ((1 << e) >= n)
84
break;
85
return e;
86
}
87
88
/* Set inv to the inverse of d, in the style of invert_limb, ie. for
89
udiv_qrnnd_preinv. */
90
void
91
mpz_preinv_invert (mpz_t inv, mpz_t d, int numb_bits)
92
{
93
mpz_t t;
94
int norm;
95
assert (SIZ(d) > 0);
96
97
norm = numb_bits - mpz_sizeinbase (d, 2);
98
assert (norm >= 0);
99
mpz_init_set_ui (t, 1L);
100
mpz_mul_2exp (t, t, 2*numb_bits - norm);
101
mpz_tdiv_q (inv, t, d);
102
mpz_set_ui (t, 1L);
103
mpz_mul_2exp (t, t, numb_bits);
104
mpz_sub (inv, inv, t);
105
106
mpz_clear (t);
107
}
108
109
/* Calculate r satisfying r*d == 1 mod 2^n. */
110
void
111
mpz_invert_2exp (mpz_t r, mpz_t a, unsigned long n)
112
{
113
unsigned long i;
114
mpz_t inv, prod;
115
116
assert (mpz_odd_p (a));
117
118
mpz_init_set_ui (inv, 1L);
119
mpz_init (prod);
120
121
for (i = 1; i < n; i++)
122
{
123
mpz_mul (prod, inv, a);
124
if (mpz_tstbit (prod, i) != 0)
125
mpz_setbit (inv, i);
126
}
127
128
mpz_mul (prod, inv, a);
129
mpz_tdiv_r_2exp (prod, prod, n);
130
assert (mpz_cmp_ui (prod, 1L) == 0);
131
132
mpz_set (r, inv);
133
134
mpz_clear (inv);
135
mpz_clear (prod);
136
}
137
138
/* Calculate inv satisfying r*a == 1 mod 2^n. */
139
void
140
mpz_invert_ui_2exp (mpz_t r, unsigned long a, unsigned long n)
141
{
142
mpz_t az;
143
mpz_init_set_ui (az, a);
144
mpz_invert_2exp (r, az, n);
145
mpz_clear (az);
146
}
147
148