1 /*
2 Copyright (C) 2014 William Hart
3 Copyright (C) 2020 Fredrik Johansson
4
5 This file is part of FLINT.
6
7 FLINT is free software: you can redistribute it and/or modify it under
8 the terms of the GNU Lesser General Public License (LGPL) as published
9 by the Free Software Foundation; either version 2.1 of the License, or
10 (at your option) any later version. See <http://www.gnu.org/licenses/>.
11 */
12
13 #define ulong ulongxx /* interferes with system includes */
14 #include <stdlib.h>
15 #include <stdio.h>
16 #include <math.h>
17 #undef ulong
18 #define ulong mp_limb_t
19 #include <gmp.h>
20 #include "flint.h"
21 #include "ulong_extras.h"
22 #include "fmpz.h"
23 #include "fmpz_vec.h"
24 #include "aprcl.h"
25 #include "mpn_extras.h"
26
fmpz_is_prime(const fmpz_t n)27 int fmpz_is_prime(const fmpz_t n)
28 {
29 double logd;
30 ulong p, ppi, limit;
31 ulong * pp1, * pm1;
32 slong i, l, num, num_pp1, num_pm1;
33 const ulong * primes;
34 const double * pinv;
35
36 fmpz_t F1, Fsqr, Fcub, R, t;
37 int res = -1;
38
39 if (fmpz_cmp_ui(n, 1) <= 0)
40 return 0;
41
42 if (fmpz_abs_fits_ui(n))
43 return n_is_prime(fmpz_get_ui(n));
44
45 if (fmpz_is_even(n))
46 return 0;
47
48 if (flint_mpn_factor_trial(COEFF_TO_PTR(*n)->_mp_d, COEFF_TO_PTR(*n)->_mp_size, 1, fmpz_bits(n)))
49 return 0;
50
51 /* todo: use fmpz_is_perfect_power? */
52 if (fmpz_is_square(n))
53 return 0;
54
55 /* Fast deterministic Miller-Rabin test up to about 81 bits. This choice of
56 bases certifies primality for n < 3317044064679887385961981;
57 see https://doi.org/10.1090/mcom/3134 */
58 fmpz_init(t);
59 fmpz_tdiv_q_2exp(t, n, 64);
60 if (fmpz_cmp_ui(t, 179817) < 0)
61 {
62 static const char bases[] = { 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 0 };
63
64 for (i = 0; bases[i] != 0; i++)
65 {
66 fmpz_set_ui(t, bases[i]);
67 if (!fmpz_is_strong_probabprime(n, t))
68 return 0; /* no need to clear t since it is small */
69 }
70
71 return 1;
72 }
73
74 /* Do a single base-2 test to rule out most composites */
75 fmpz_set_ui(t, 2);
76 if (!fmpz_is_strong_probabprime(n, t))
77 return 0;
78
79 fmpz_clear(t);
80
81 logd = fmpz_dlog(n);
82 limit = (ulong) (logd*logd*logd/100.0) + 20;
83
84 fmpz_init(F1);
85 fmpz_init(R);
86 fmpz_init(Fsqr);
87 fmpz_init(Fcub);
88
89 for (l = 0; l < 4 && res == -1; l++, limit *= 10)
90 {
91 num_pm1 = num_pp1 = 0;
92
93 /* number of primes multiplied that will fit in a word */
94 num = FLINT_BITS/FLINT_BIT_COUNT(limit);
95
96 /* compute remainders of n mod p for primes p up to limit (approx.) */
97
98 n_prime_pi_bounds(&ppi, &ppi, limit); /* precompute primes */
99 primes = n_primes_arr_readonly(ppi + FLINT_BITS);
100 pinv = n_prime_inverses_arr_readonly(ppi + FLINT_BITS);
101
102 pm1 = _nmod_vec_init(2 + (ulong) logd); /* space for primes dividing n - 1 */
103 pp1 = _nmod_vec_init(2 + (ulong) logd); /* space for primes dividing n + 1 */
104
105 while (primes[0] < limit)
106 {
107 /* multiply batch of primes */
108
109 p = primes[0];
110 for (i = 1; i < num; i++)
111 p *= primes[i];
112
113 /* multi-modular reduction */
114
115 p = fmpz_tdiv_ui(n, p);
116
117 /* check for factors */
118 for (i = 0; i < num; i++)
119 {
120 ulong r = n_mod2_precomp(p, primes[i], pinv[i]);
121
122 if (r == 1) /* n - 1 = 0 mod p */
123 pm1[num_pm1++] = primes[i];
124
125 if (r == primes[i] - 1) /* n + 1 = 0 mod p */
126 pp1[num_pp1++] = primes[i];
127 }
128
129 /* get next batch of primes */
130 primes += num;
131 pinv += num;
132 }
133
134 /* p - 1 test */
135 res = fmpz_is_prime_pocklington(F1, R, n, pm1, num_pm1);
136
137 if (res == 1)
138 {
139 fmpz_mul(Fsqr, F1, F1);
140 if (fmpz_cmp(Fsqr, n) < 0)
141 {
142 fmpz_mul(Fcub, Fsqr, F1);
143 if (fmpz_cmp(Fcub, n) >= 0) /* Brillhart, Lehmer, Selfridge test */
144 {
145 fmpz_t n1, c2, c1;
146
147 fmpz_init(n1);
148 fmpz_init(c2);
149 fmpz_init(c1);
150
151 fmpz_sub_ui(n1, n, 1); /* n is 1 mod F1 */
152 fmpz_tdiv_q(n1, n1, F1);
153
154 fmpz_tdiv_qr(c2, c1, n1, F1); /* Let n = c2*F^2 + c1*F + 1 */
155
156 fmpz_mul(c1, c1, c1); /* check if c1^2 - 4*c2 is a square */
157 fmpz_submul_ui(c1, c2, 4);
158
159 if (fmpz_is_square(c1))
160 res = 0;
161 /* else n is prime (res == 1) */
162
163 fmpz_clear(n1);
164 fmpz_clear(c2);
165 fmpz_clear(c1);
166 } else /* p + 1 test */
167 {
168 fmpz_t F2, Fm1;
169
170 fmpz_init(F2);
171 fmpz_init(Fm1);
172
173 res = fmpz_is_prime_morrison(F2, R, n, pp1, num_pp1);
174
175 if (res == 1)
176 {
177 fmpz_sub_ui(Fm1, F2, 1); /* need F2 - 1 > sqrt(n) */
178 fmpz_mul(Fsqr, Fm1, Fm1);
179
180 if (fmpz_cmp(Fsqr, n) <= 0)
181 {
182 fmpz_mul(Fcub, Fsqr, Fm1);
183
184 if (fmpz_cmp(Fcub, n) > 0) /* Improved n + 1 test */
185 {
186 fmpz_t r1, r0, b, r, t;
187
188 fmpz_init(r1);
189 fmpz_init(r0);
190 fmpz_init(b);
191 fmpz_init(r);
192 fmpz_init(t);
193
194 fmpz_tdiv_qr(r1, r0, R, F2); /* R = r1*F2 + r0 */
195
196 /* check if x^2 + r0*x - r1 has positive integral root */
197 fmpz_mul(t, r0, r0); /* b = sqrt(r0^2 - 4(-r1)) */
198 fmpz_addmul_ui(t, r1, 4);
199 fmpz_sqrtrem(b, r, t);
200
201 if (fmpz_is_zero(r) && fmpz_cmp(b, r0) > 0) /* if so, composite */
202 res = 0;
203
204 /* check if x^2 + (r0 - F2)*x - r1 - 1 has positive integral root */
205 fmpz_sub(r0, r0, F2);
206 fmpz_add_ui(r1, r1, 1);
207
208 fmpz_mul(t, r0, r0); /* b = sqrt((r0 - F2)^2 - 4(-r1 - 1)) */
209 fmpz_addmul_ui(t, r1, 4);
210 fmpz_sqrtrem(b, r, t);
211
212 if (fmpz_is_zero(r) && fmpz_cmp(b, r0) > 0) /* if so, composite */
213 res = 0;
214
215 fmpz_clear(t);
216 fmpz_clear(b);
217 fmpz_clear(r);
218 fmpz_clear(r1);
219 fmpz_clear(r0);
220 } else /* Brillhart, Lehmer, Selfridge combined p-1, p+1 test */
221 {
222 fmpz_t F, nmodF;
223
224 fmpz_init(F);
225
226 fmpz_mul(F, F1, F2); /* F = lcm(F1, F2), F1 | n - 1, F2 | n + 1 */
227 if (fmpz_is_even(F1) && fmpz_is_even(F2))
228 fmpz_tdiv_q_2exp(F, F, 1);
229
230 fmpz_mul(Fsqr, F, F);
231
232 if (fmpz_cmp(Fsqr, n) > 0) /* lcm(F1, F2) > sqrt(n) */
233 {
234 fmpz_init(nmodF);
235
236 fmpz_mod(nmodF, n, F); /* check n mod F not factor of n */
237
238 if (!fmpz_equal(nmodF, n) && !fmpz_is_one(nmodF)
239 && fmpz_divisible(n, nmodF))
240 res = 0;
241
242 fmpz_clear(nmodF);
243 } else
244 {
245 fmpz_t d;
246
247 fmpz_init(d);
248
249 fmpz_mul(Fcub, Fsqr, F);
250
251 if (fmpz_cmp(Fcub, n) > 0) /* Lenstra's divisors in residue class */
252 {
253 fmpz_t r;
254
255 fmpz_init(r);
256
257 fmpz_set_ui(r, 1);
258 if (fmpz_divisor_in_residue_class_lenstra(d, n, r, F))
259 res = 0;
260
261 fmpz_mod(r, n, F);
262 if (fmpz_divisor_in_residue_class_lenstra(d, n, r, F))
263 res = 0;
264
265 fmpz_clear(r);
266 } else /* apr-cl primality test */
267 {
268 res = aprcl_is_prime(n);
269 }
270
271 fmpz_clear(d);
272 }
273
274 fmpz_clear(F);
275 }
276 }
277 /* else n is prime, i.e. res = 1 */
278 }
279
280 fmpz_clear(F2);
281 fmpz_clear(Fm1);
282 }
283 }
284 }
285
286 _nmod_vec_clear(pm1);
287 _nmod_vec_clear(pp1);
288
289 }
290
291 /* aprcl_is_prime() actually throws, but it does not hurt to have
292 this fallback here */
293 if (res < 0)
294 {
295 flint_printf("Exception in fmpz_is_prime: failed to prove ");
296 fmpz_print(n);
297 flint_printf(" prime or composite\n");
298 flint_abort();
299 }
300
301 fmpz_clear(F1);
302 fmpz_clear(R);
303 fmpz_clear(Fsqr);
304 fmpz_clear(Fcub);
305
306 return res;
307 }
308