1 #include <gmp.h>
2 #include "ptypes.h"
3
4 #include "primality.h"
5 #include "gmp_main.h" /* primality_pretest */
6 #include "bls75.h"
7 #include "ecpp.h"
8 #include "factor.h"
9
10 #define FUNC_is_perfect_square 1
11 #define FUNC_mpz_logn
12 #include "utility.h"
13
14 #define NSMALLPRIMES 54
15 static const unsigned char sprimes[NSMALLPRIMES] = {2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79,83,89,97,101,103,107,109,113,127,131,137,139,149,151,157,163,167,173,179,181,191,193,197,199,211,223,227,229,233,239,241,251};
16
17
18 /* 0 fail, 1 pass, -1 nothing. Modifies base a */
_preprocess_base(mpz_t n,mpz_t a)19 static int _preprocess_base(mpz_t n, mpz_t a)
20 {
21 if (mpz_cmp_ui(a, 1) <= 0)
22 croak("Base %ld is invalid", mpz_get_si(a));
23 if (mpz_cmp_ui(n, 3) <= 0)
24 return (mpz_cmp_ui(n, 2) >= 0);
25
26 if (mpz_cmp_ui(a, 2) > 0) {
27 if (mpz_cmp(a, n) >= 0) {
28 mpz_mod(a, a, n);
29 if (mpz_cmp_ui(a, 1) <= 0)
30 return mpz_sgn(a);
31 }
32 }
33 return -1;
34 }
35
is_pseudoprime(mpz_t n,mpz_t a)36 int is_pseudoprime(mpz_t n, mpz_t a)
37 {
38 mpz_t nm1;
39 int res;
40
41 if ((res = _preprocess_base(n, a)) >= 0)
42 return res;
43
44 mpz_init(nm1);
45 mpz_sub_ui(nm1, n, 1);
46 mpz_powm(nm1, a, nm1, n);
47 res = (mpz_cmp_ui(nm1, 1) == 0);
48 mpz_clear(nm1);
49 return res;
50 }
51
is_euler_pseudoprime(mpz_t n,mpz_t a)52 int is_euler_pseudoprime(mpz_t n, mpz_t a)
53 {
54 mpz_t nm1, ap;
55 int res;
56
57 if (mpz_even_p(n))
58 return (mpz_cmp_ui(n,2) == 0);
59 if ((res = _preprocess_base(n, a)) >= 0)
60 return res;
61
62 mpz_init(ap);
63 if (mpz_gcd(ap, a, n), mpz_cmp_ui(ap, 1) != 0) {
64 mpz_clear(ap);
65 return 0;
66 }
67 mpz_init(nm1);
68 mpz_sub_ui(nm1, n, 1);
69 mpz_tdiv_q_2exp(ap, nm1, 1);
70 mpz_powm(ap, a, ap, n);
71
72 if (mpz_cmp_ui(ap, 1) && mpz_cmp(ap, nm1))
73 res = 0;
74 else if (mpz_kronecker(a, n) >= 0)
75 res = (mpz_cmp_ui(ap, 1) == 0);
76 else
77 res = (mpz_cmp(ap, nm1) == 0);
78
79 mpz_clear(nm1);
80 mpz_clear(ap);
81 return res;
82 }
83
mrx(mpz_t x,mpz_t d,mpz_t n,UV s)84 static int mrx(/*destroyed*/mpz_t x, /*destroyed*/ mpz_t d, mpz_t n, UV s)
85 {
86 UV r;
87 mpz_powm(x, x, d, n);
88 mpz_sub_ui(d, n, 1);
89 if (!mpz_cmp_ui(x, 1) || !mpz_cmp(x, d))
90 return 1;
91 for (r = 1; r < s; r++) {
92 mpz_powm_ui(x, x, 2, n);
93 if (!mpz_cmp_ui(x, 1))
94 break;
95 if (!mpz_cmp(x, d))
96 return 1;
97 }
98 return 0;
99 }
100
101
miller_rabin(mpz_t n,mpz_t a)102 int miller_rabin(mpz_t n, mpz_t a)
103 {
104 mpz_t d, x;
105 int cmpr, rval = 1;
106
107 cmpr = mpz_cmp_ui(n, 2);
108 if (cmpr == 0) return 1; /* 2 is prime */
109 if (cmpr < 0) return 0; /* below 2 is composite */
110 if (mpz_even_p(n)) return 0; /* multiple of 2 is composite */
111 if (mpz_cmp_ui(a, 1) <= 0) croak("Base %ld is invalid", mpz_get_si(a));
112
113 mpz_init_set(x, a);
114 mpz_init_set(d, n);
115 mpz_sub_ui(d, d, 1);
116
117 /* Handle large and small bases. Use x so we don't modify their input a. */
118 if (mpz_cmp(x, n) >= 0)
119 mpz_mod(x, x, n);
120 if ( (mpz_cmp_ui(x, 1) > 0) && (mpz_cmp(x, d) < 0) ) {
121 UV s = mpz_scan1(d, 0);
122 mpz_tdiv_q_2exp(d, d, s);
123 rval = mrx(x, d, n, s);
124 }
125 mpz_clear(d);
126 mpz_clear(x);
127 return rval;
128 }
miller_rabin_ui(mpz_t n,unsigned long a)129 int miller_rabin_ui(mpz_t n, unsigned long a)
130 {
131 mpz_t d, x;
132 int cmpr, rval = 1;
133
134 cmpr = mpz_cmp_ui(n, 2);
135 if (cmpr == 0) return 1; /* 2 is prime */
136 if (cmpr < 0) return 0; /* below 2 is composite */
137 if (mpz_even_p(n)) return 0; /* multiple of 2 is composite */
138 if (a <= 1) croak("Base %lu is invalid", a);
139
140 mpz_init_set_ui(x, a);
141 mpz_init_set(d, n);
142 mpz_sub_ui(d, d, 1);
143
144 if (mpz_cmp(x, n) >= 0)
145 mpz_mod(x, x, n);
146 if ( (mpz_cmp_ui(x, 1) > 0) && (mpz_cmp(x, d) < 0) ) {
147 UV s = mpz_scan1(d, 0);
148 mpz_tdiv_q_2exp(d, d, s);
149 rval = mrx(x, d, n, s);
150 }
151 mpz_clear(d);
152 mpz_clear(x);
153 return rval;
154 }
155
is_miller_prime(mpz_t n,int assume_grh)156 int is_miller_prime(mpz_t n, int assume_grh)
157 {
158 mpz_t d, x, D;
159 UV s, maxa, a;
160 int rval;
161
162 {
163 int cmpr = mpz_cmp_ui(n, 2);
164 if (cmpr == 0) return 1; /* 2 is prime */
165 if (cmpr < 0) return 0; /* below 2 is composite */
166 if (mpz_even_p(n)) return 0; /* multiple of 2 is composite */
167 }
168
169 if (mpz_cmp_ui(n, 1373653) < 0) {
170 maxa = 3;
171 } else if (assume_grh) {
172 double logn = mpz_logn(n);
173 double dmaxa = 2 * logn * logn; /* Bach (1990) */
174 /* Wedeniwski (2001) claims the following, but it might be wrong
175 * double dmaxa = 1.5L * logn * logn - 44.0L/5.0L * logn + 13; */
176 if (dmaxa >= (double)ULONG_MAX)
177 croak("is_miller_prime: n is too large for GRH DMR");
178 maxa = ceil(dmaxa);
179 } else { /* Bober and Goldmakher 2015 (http://arxiv.org/abs/1311.7556) */
180 /* n_p < p^(1/(4*sqrt(e))+epsilon). Do it with logs */
181 double dmaxa = exp( (1.0L/6.5948850828L) * mpz_logn(n) );
182 if (dmaxa >= (double)ULONG_MAX)
183 croak("is_miller_prime: n is too large for unconditional DMR");
184 maxa = ceil(dmaxa);
185 }
186
187 if (mpz_cmp_ui(n, maxa) <= 0)
188 maxa = mpz_get_ui(n) - 1;
189 if (get_verbose_level() > 1)
190 printf("Deterministic Miller-Rabin testing bases from 2 to %"UVuf"\n", maxa);
191
192 mpz_init_set(d, n);
193 mpz_sub_ui(d, d, 1);
194 s = mpz_scan1(d, 0);
195 mpz_tdiv_q_2exp(d, d, s);
196 mpz_init(D);
197 mpz_init(x);
198
199 for (a = 2, rval = 1; rval && a <= maxa; a++) {
200 mpz_set_ui(x, a);
201 mpz_set(D, d);
202 rval = mrx(x, D, n, s);
203 }
204 mpz_clear(x);
205 mpz_clear(D);
206 mpz_clear(d);
207 return rval;
208 }
209
miller_rabin_random(mpz_t n,UV numbases,char * seedstr)210 int miller_rabin_random(mpz_t n, UV numbases, char* seedstr)
211 {
212 gmp_randstate_t randstate;
213 mpz_t t, base;
214 UV i;
215
216 if (numbases == 0) return 1;
217 if (mpz_cmp_ui(n, 100) < 0) /* tiny n */
218 return (_GMP_is_prob_prime(n) > 0);
219
220 /* See if they've asked for a ludicrous number of bases */
221 mpz_init(t);
222 mpz_mul_ui(t, n, 3);
223 mpz_div_ui(t, t, 4);
224 if (mpz_cmp_ui(t, numbases) <= 0) {
225 int res = is_bpsw_dmr_prime(n);
226 if (res != 1) {
227 mpz_clear(t);
228 return !!res;
229 }
230 numbases = mpz_get_ui(t);
231 }
232
233 mpz_init(base);
234 mpz_sub_ui(t, n, 3);
235
236 if (seedstr == 0) { /* Normal case with no seed string. Use our CSPRNG. */
237 for (i = 0; i < numbases; i++) {
238 mpz_isaac_urandomm(base, t); /* base 0 .. (n-3)-1 */
239 mpz_add_ui(base, base, 2); /* base 2 .. n-2 */
240 if (miller_rabin(n, base) == 0)
241 break;
242 }
243 } else {
244 gmp_randinit_default(randstate);
245 mpz_set_str(base, seedstr, 0);
246 gmp_randseed(randstate, base);
247 for (i = 0; i < numbases; i++) {
248 mpz_urandomm(base, randstate, t); /* base 0 .. (n-3)-1 */
249 mpz_add_ui(base, base, 2); /* base 2 .. n-2 */
250 if (miller_rabin(n, base) == 0)
251 break;
252 }
253 gmp_randclear(randstate);
254 }
255 mpz_clear(base); mpz_clear(t);
256 return (i >= numbases);
257 }
258
259
is_euler_plumb_pseudoprime(mpz_t n)260 int is_euler_plumb_pseudoprime(mpz_t n)
261 {
262 unsigned int nmod8;
263 mpz_t x, two;
264 int result = 0;
265 if (mpz_cmp_ui(n,5) < 0)
266 return (mpz_cmp_ui(n,2) == 0 || mpz_cmp_ui(n,3) == 0);
267 if (mpz_even_p(n))
268 return 0;
269 nmod8 = mpz_fdiv_ui(n, 8);
270 mpz_init(x); mpz_init_set_ui(two,2);
271 mpz_sub_ui(x, n, 1);
272 mpz_fdiv_q_2exp(x, x, 1 + (nmod8 == 1));
273 mpz_powm(x, two, x, n);
274 if (mpz_cmp_ui(x,1) == 0) {
275 result = (nmod8 == 1 || nmod8 == 7);
276 } else {
277 mpz_add_ui(x,x,1);
278 result = (mpz_cmp(x,n) == 0 && (nmod8 == 1 || nmod8 == 3 || nmod8 == 5));
279 }
280 mpz_clear(two); mpz_clear(x);
281 return result;
282 }
283
284 /* Returns Lucas sequence U_k mod n and V_k mod n defined by P,Q */
lucas_seq(mpz_t U,mpz_t V,mpz_t n,IV P,IV Q,mpz_t k,mpz_t Qk,mpz_t t)285 void lucas_seq(mpz_t U, mpz_t V, mpz_t n, IV P, IV Q, mpz_t k,
286 mpz_t Qk, mpz_t t)
287 {
288 UV b = mpz_sizeinbase(k, 2);
289 IV D = P*P - 4*Q;
290
291 if (mpz_cmp_ui(n, 2) < 0) croak("Lucas sequence modulus n must be > 1");
292 MPUassert( mpz_cmp_ui(k, 0) >= 0, "lucas_seq: k is negative" );
293 MPUassert( mpz_cmp_si(n,(P>=0) ? P : -P) > 0, "lucas_seq: P is out of range");
294 MPUassert( mpz_cmp_si(n,(Q>=0) ? Q : -Q) > 0, "lucas_seq: Q is out of range");
295 MPUassert( D != 0, "lucas_seq: D is zero" );
296
297 if (mpz_cmp_ui(k, 0) <= 0) {
298 mpz_set_ui(U, 0);
299 mpz_set_ui(V, 2);
300 return;
301 }
302 if (mpz_even_p(n)) {
303 /* If n is even, then we can't divide by 2. Do it differently. */
304 alt_lucas_seq(U, V, n, P, Q, k, Qk, t);
305 return;
306 }
307
308 mpz_set_ui(U, 1);
309 mpz_set_si(V, P);
310 mpz_set_si(Qk, Q);
311
312 if (Q == 1) {
313 /* Use the fast V method if possible. Much faster with small n. */
314 mpz_set_si(t, P*P-4);
315 if (P > 2 && mpz_invert(t, t, n)) {
316 /* Compute V_k and V_{k+1}, then computer U_k from them. */
317 mpz_set_si(V, P);
318 mpz_set_si(U, P*P-2);
319 while (b > 1) {
320 b--;
321 if (mpz_tstbit(k, b-1)) {
322 mpz_mul(V, V, U); mpz_sub_ui(V, V, P); mpz_mod(V, V, n);
323 mpz_mul(U, U, U); mpz_sub_ui(U, U, 2); mpz_mod(U, U, n);
324 } else {
325 mpz_mul(U, V, U); mpz_sub_ui(U, U, P); mpz_mod(U, U, n);
326 mpz_mul(V, V, V); mpz_sub_ui(V, V, 2); mpz_mod(V, V, n);
327 }
328 }
329 mpz_mul_ui(U, U, 2);
330 mpz_submul_ui(U, V, P);
331 mpz_mul(U, U, t);
332 } else {
333 /* Fast computation of U_k and V_k, specific to Q = 1 */
334 while (b > 1) {
335 mpz_mulmod(U, U, V, n, t); /* U2k = Uk * Vk */
336 mpz_mul(V, V, V);
337 mpz_sub_ui(V, V, 2);
338 mpz_mod(V, V, n); /* V2k = Vk^2 - 2 Q^k */
339 b--;
340 if (mpz_tstbit(k, b-1)) {
341 mpz_mul_si(t, U, D);
342 /* U: U2k+1 = (P*U2k + V2k)/2 */
343 mpz_mul_si(U, U, P);
344 mpz_add(U, U, V);
345 if (mpz_odd_p(U)) mpz_add(U, U, n);
346 mpz_fdiv_q_2exp(U, U, 1);
347 /* V: V2k+1 = (D*U2k + P*V2k)/2 */
348 mpz_mul_si(V, V, P);
349 mpz_add(V, V, t);
350 if (mpz_odd_p(V)) mpz_add(V, V, n);
351 mpz_fdiv_q_2exp(V, V, 1);
352 }
353 }
354 }
355 } else {
356 while (b > 1) {
357 mpz_mulmod(U, U, V, n, t); /* U2k = Uk * Vk */
358 mpz_mul(V, V, V);
359 mpz_submul_ui(V, Qk, 2);
360 mpz_mod(V, V, n); /* V2k = Vk^2 - 2 Q^k */
361 mpz_mul(Qk, Qk, Qk); /* Q2k = Qk^2 */
362 b--;
363 if (mpz_tstbit(k, b-1)) {
364 mpz_mul_si(t, U, D);
365 /* U: U2k+1 = (P*U2k + V2k)/2 */
366 mpz_mul_si(U, U, P);
367 mpz_add(U, U, V);
368 if (mpz_odd_p(U)) mpz_add(U, U, n);
369 mpz_fdiv_q_2exp(U, U, 1);
370 /* V: V2k+1 = (D*U2k + P*V2k)/2 */
371 mpz_mul_si(V, V, P);
372 mpz_add(V, V, t);
373 if (mpz_odd_p(V)) mpz_add(V, V, n);
374 mpz_fdiv_q_2exp(V, V, 1);
375
376 mpz_mul_si(Qk, Qk, Q);
377 }
378 mpz_mod(Qk, Qk, n);
379 }
380 }
381 mpz_mod(U, U, n);
382 mpz_mod(V, V, n);
383 }
384
alt_lucas_seq(mpz_t Uh,mpz_t Vl,mpz_t n,IV P,IV Q,mpz_t k,mpz_t Ql,mpz_t t)385 void alt_lucas_seq(mpz_t Uh, mpz_t Vl, mpz_t n, IV P, IV Q, mpz_t k,
386 mpz_t Ql, mpz_t t)
387 {
388 mpz_t Vh, Qh;
389 int j, s = mpz_scan1(k,0), b = mpz_sizeinbase(k,2);
390
391 if (mpz_sgn(k) <= 0) {
392 mpz_set_ui(Uh, 0);
393 mpz_set_ui(Vl, 2);
394 return;
395 }
396
397 mpz_set_ui(Uh, 1);
398 mpz_set_ui(Vl, 2);
399 mpz_set_ui(Ql, 1);
400 mpz_init_set_si(Vh,P);
401 mpz_init_set_ui(Qh,1);
402
403 for (j = b; j > s; j--) {
404 mpz_mul(Ql, Ql, Qh);
405 mpz_mod(Ql, Ql, n);
406 if (mpz_tstbit(k, j)) {
407 mpz_mul_si(Qh, Ql, Q);
408 mpz_mul(Uh, Uh, Vh);
409 mpz_mul_si(t, Ql, P); mpz_mul(Vl, Vl, Vh); mpz_sub(Vl, Vl, t);
410 mpz_mul(Vh, Vh, Vh); mpz_sub(Vh, Vh, Qh); mpz_sub(Vh, Vh, Qh);
411 } else {
412 mpz_set(Qh, Ql);
413 mpz_mul(Uh, Uh, Vl); mpz_sub(Uh, Uh, Ql);
414 mpz_mul_si(t, Ql, P); mpz_mul(Vh, Vh, Vl); mpz_sub(Vh, Vh, t);
415 mpz_mul(Vl, Vl, Vl); mpz_sub(Vl, Vl, Ql); mpz_sub(Vl, Vl, Ql);
416 }
417 mpz_mod(Qh, Qh, n);
418 mpz_mod(Uh, Uh, n);
419 mpz_mod(Vh, Vh, n);
420 mpz_mod(Vl, Vl, n);
421 }
422 mpz_mul(Ql, Ql, Qh);
423 mpz_mul_si(Qh, Ql, Q);
424 mpz_mul(Uh, Uh, Vl); mpz_sub(Uh, Uh, Ql);
425 mpz_mul_si(t, Ql, P); mpz_mul(Vl, Vl, Vh); mpz_sub(Vl, Vl, t);
426 mpz_mul(Ql, Ql, Qh);
427 mpz_clear(Qh); mpz_clear(Vh);
428 mpz_mod(Ql, Ql, n);
429 mpz_mod(Uh, Uh, n);
430 mpz_mod(Vl, Vl, n);
431 for (j = 0; j < s; j++) {
432 mpz_mul(Uh, Uh, Vl);
433 mpz_mul(Vl, Vl, Vl); mpz_sub(Vl, Vl, Ql); mpz_sub(Vl, Vl, Ql);
434 mpz_mul(Ql, Ql, Ql);
435 mpz_mod(Ql, Ql, n);
436 mpz_mod(Uh, Uh, n);
437 mpz_mod(Vl, Vl, n);
438 }
439 }
440
lucasuv(mpz_t Uh,mpz_t Vl,IV P,IV Q,mpz_t k)441 void lucasuv(mpz_t Uh, mpz_t Vl, IV P, IV Q, mpz_t k)
442 {
443 mpz_t Vh, Ql, Qh, t;
444 int j, s, n;
445
446 if (mpz_sgn(k) <= 0) {
447 mpz_set_ui(Uh, 0);
448 mpz_set_ui(Vl, 2);
449 return;
450 }
451
452 mpz_set_ui(Uh, 1);
453 mpz_set_ui(Vl, 2);
454 mpz_init_set_si(Vh,P);
455 mpz_init(t);
456
457 s = mpz_scan1(k, 0); /* number of zero bits at the end */
458 n = mpz_sizeinbase(k,2);
459
460 /* It is tempting to try to pull out the various Q operations when Q=1 or
461 * Q=-1. This doesn't lead to any immediate savings. Don't bother unless
462 * there is a way to reduce the actual operations involving U and V. */
463 mpz_init_set_ui(Ql,1);
464 mpz_init_set_ui(Qh,1);
465
466 for (j = n; j > s; j--) {
467 mpz_mul(Ql, Ql, Qh);
468 if (mpz_tstbit(k, j)) {
469 mpz_mul_si(Qh, Ql, Q);
470 mpz_mul(Uh, Uh, Vh);
471 mpz_mul_si(t, Ql, P); mpz_mul(Vl, Vl, Vh); mpz_sub(Vl, Vl, t);
472 mpz_mul(Vh, Vh, Vh); mpz_sub(Vh, Vh, Qh); mpz_sub(Vh, Vh, Qh);
473 } else {
474 mpz_set(Qh, Ql);
475 mpz_mul(Uh, Uh, Vl); mpz_sub(Uh, Uh, Ql);
476 mpz_mul_si(t, Ql, P); mpz_mul(Vh, Vh, Vl); mpz_sub(Vh, Vh, t);
477 mpz_mul(Vl, Vl, Vl); mpz_sub(Vl, Vl, Ql); mpz_sub(Vl, Vl, Ql);
478 }
479 }
480 mpz_mul(Ql, Ql, Qh);
481 mpz_mul_si(Qh, Ql, Q);
482 mpz_mul(Uh, Uh, Vl); mpz_sub(Uh, Uh, Ql);
483 mpz_mul_si(t, Ql, P); mpz_mul(Vl, Vl, Vh); mpz_sub(Vl, Vl, t);
484 mpz_mul(Ql, Ql, Qh);
485 mpz_clear(Qh); mpz_clear(t); mpz_clear(Vh);
486 for (j = 0; j < s; j++) {
487 mpz_mul(Uh, Uh, Vl);
488 mpz_mul(Vl, Vl, Vl); mpz_sub(Vl, Vl, Ql); mpz_sub(Vl, Vl, Ql);
489 mpz_mul(Ql, Ql, Ql);
490 }
491 mpz_clear(Ql);
492 }
493
lucas_lehmer(UV p)494 int lucas_lehmer(UV p)
495 {
496 UV k, tlim;
497 int res, pbits;
498 mpz_t V, mp, t;
499
500 if (p == 2) return 1;
501 if (!(p&1)) return 0;
502
503 mpz_init_set_ui(t, p);
504 if (!_GMP_is_prob_prime(t)) /* p must be prime */
505 { mpz_clear(t); return 0; }
506 if (p < 23)
507 { mpz_clear(t); return (p != 11); }
508
509 pbits = mpz_sizeinbase(t,2);
510 mpz_init(mp);
511 mpz_setbit(mp, p);
512 mpz_sub_ui(mp, mp, 1);
513
514 /* If p=3 mod 4 and p,2p+1 both prime, then 2p+1 | 2^p-1. Cheap test. */
515 if (p > 3 && p % 4 == 3) {
516 mpz_mul_ui(t, t, 2);
517 mpz_add_ui(t, t, 1);
518 if (_GMP_is_prob_prime(t) && mpz_divisible_p(mp, t))
519 { mpz_clear(mp); mpz_clear(t); return 0; }
520 }
521
522 /* Do a little trial division first. Saves quite a bit of time. */
523 tlim = (p < 1500) ? p/2 : (p < 5000) ? p : 2*p;
524 if (tlim > UV_MAX/(2*p)) tlim = UV_MAX/(2*p);
525 for (k = 1; k < tlim; k++) {
526 UV q = 2*p*k+1;
527 if ( (q%8==1 || q%8==7) && /* factor must be 1 or 7 mod 8 */
528 q % 3 && q % 5 && q % 7 && q % 11 && q % 13) { /* and must be prime */
529 if (1 && q < (UVCONST(1) << (BITS_PER_WORD/2)) ) {
530 UV b = 1, j = pbits;
531 while (j--) {
532 b = (b*b) % q;
533 if (p & (UVCONST(1) << j)) { b *= 2; if (b >= q) b -= q; }
534 }
535 if (b == 1)
536 { mpz_clear(mp); mpz_clear(t); return 0; }
537 } else {
538 if( mpz_divisible_ui_p(mp, q) )
539 { mpz_clear(mp); mpz_clear(t); return 0; }
540 }
541 }
542 }
543 /* We could do some specialized p+1 factoring here. */
544
545 mpz_init_set_ui(V, 4);
546
547 for (k = 3; k <= p; k++) {
548 mpz_mul(V, V, V);
549 mpz_sub_ui(V, V, 2);
550 /* mpz_mod(V, V, mp) but more efficiently done given mod 2^p-1 */
551 if (mpz_sgn(V) < 0) mpz_add(V, V, mp);
552 /* while (n > mp) { n = (n >> p) + (n & mp) } if (n==mp) n=0 */
553 /* but in this case we can have at most one loop plus a carry */
554 mpz_tdiv_r_2exp(t, V, p);
555 mpz_tdiv_q_2exp(V, V, p);
556 mpz_add(V, V, t);
557 while (mpz_cmp(V, mp) >= 0) mpz_sub(V, V, mp);
558 }
559 res = !mpz_sgn(V);
560 mpz_clear(t); mpz_clear(mp); mpz_clear(V);
561 return res;
562 }
563
564 /* Returns: -1 unknown, 0 composite, 2 definitely prime */
llr(mpz_t N)565 int llr(mpz_t N)
566 {
567 mpz_t v, k, V, U, Qk, t;
568 UV i, n, P;
569 int res = -1;
570
571 if (mpz_cmp_ui(N,100) <= 0) return (_GMP_is_prob_prime(N) ? 2 : 0);
572 if (mpz_even_p(N) || mpz_divisible_ui_p(N, 3)) return 0;
573 mpz_init(v); mpz_init(k);
574 mpz_add_ui(v, N, 1);
575 n = mpz_scan1(v, 0);
576 mpz_tdiv_q_2exp(k, v, n);
577 /* N = k * 2^n - 1 */
578 if (mpz_cmp_ui(k,1) == 0) {
579 res = lucas_lehmer(n) ? 2 : 0;
580 goto DONE_LLR;
581 }
582 if (mpz_sizeinbase(k,2) > n)
583 goto DONE_LLR;
584
585 mpz_init(V);
586 mpz_init(U); mpz_init(Qk); mpz_init(t);
587 if (!mpz_divisible_ui_p(k, 3)) { /* Select V for 3 not divis k */
588 lucas_seq(U, V, N, 4, 1, k, Qk, t);
589 } else if ((n % 4 == 0 || n % 4 == 3) && mpz_cmp_ui(k,3)==0) {
590 mpz_set_ui(V, 5778);
591 } else {
592 /* Öystein J. Rödseth: http://www.uib.no/People/nmaoy/papers/luc.pdf */
593 for (P=5; P < 1000; P++) {
594 mpz_set_ui(t, P-2);
595 if (mpz_jacobi(t, N) == 1) {
596 mpz_set_ui(t, P+2);
597 if (mpz_jacobi(t, N) == -1) {
598 break;
599 }
600 }
601 }
602 if (P >= 1000) {
603 mpz_clear(t); mpz_clear(Qk); mpz_clear(U);
604 mpz_clear(V);
605 goto DONE_LLR;
606 }
607 lucas_seq(U, V, N, P, 1, k, Qk, t);
608 }
609 mpz_clear(t); mpz_clear(Qk); mpz_clear(U);
610
611 for (i = 3; i <= n; i++) {
612 mpz_mul(V, V, V);
613 mpz_sub_ui(V, V, 2);
614 mpz_mod(V, V, N);
615 }
616 res = mpz_sgn(V) ? 0 : 2;
617 mpz_clear(V);
618
619 DONE_LLR:
620 if (res != -1 && get_verbose_level() > 1) printf("N shown %s with LLR\n", res ? "prime" : "composite");
621 mpz_clear(k); mpz_clear(v);
622 return res;
623 }
624 /* Returns: -1 unknown, 0 composite, 2 definitely prime */
proth(mpz_t N)625 int proth(mpz_t N)
626 {
627 mpz_t v, k, a;
628 UV n;
629 int i, res = -1;
630 /* TODO: Should have a flag to skip pretests */
631 if (mpz_cmp_ui(N,100) <= 0) return (_GMP_is_prob_prime(N) ? 2 : 0);
632 if (mpz_even_p(N) || mpz_divisible_ui_p(N, 3)) return 0;
633 mpz_init(v); mpz_init(k);
634 mpz_sub_ui(v, N, 1);
635 n = mpz_scan1(v, 0);
636 mpz_tdiv_q_2exp(k, v, n);
637 /* N = k * 2^n + 1 */
638 if (mpz_sizeinbase(k,2) <= n) {
639 mpz_init(a);
640 for (i = 0; i < 25 && res == -1; i++) {
641 mpz_set_ui(a, sprimes[i]);
642 if (mpz_jacobi(a, N) == -1)
643 res = 0;
644 }
645 /* (a,N) = -1 if res=0. Do Proth primality test. */
646 if (res == 0) {
647 mpz_tdiv_q_2exp(k, v, 1);
648 mpz_powm(a, a, k, N);
649 if (mpz_cmp(a, v) == 0)
650 res = 2;
651 }
652 mpz_clear(a);
653 }
654 /* TODO: look into Rao (2018): k*2^n+1 for n>1, k prime */
655 /* if (n > 1 && res == -1 && _GMP_BPSW(k)) { ... } */
656 mpz_clear(k); mpz_clear(v);
657 if (res != -1 && get_verbose_level() > 1) printf("N shown %s with Proth\n", res ? "prime" : "composite");
658 fflush(stdout);
659 return res;
660 }
is_proth_form(mpz_t N)661 int is_proth_form(mpz_t N)
662 {
663 mpz_t v, k;
664 UV n;
665 int res = 0;
666 if (mpz_cmp_ui(N,100) <= 0) return (_GMP_is_prob_prime(N) ? 2 : 0);
667 if (mpz_even_p(N) || mpz_divisible_ui_p(N, 3)) return 0;
668 mpz_init(v); mpz_init(k);
669 mpz_sub_ui(v, N, 1);
670 n = mpz_scan1(v, 0);
671 mpz_tdiv_q_2exp(k, v, n);
672 /* N = k * 2^n + 1 */
673 if (mpz_sizeinbase(k,2) <= n) res = 1;
674 mpz_clear(k); mpz_clear(v);
675 return res;
676 }
677
lucas_selfridge_params(IV * P,IV * Q,mpz_t n,mpz_t t)678 static int lucas_selfridge_params(IV* P, IV* Q, mpz_t n, mpz_t t)
679 {
680 IV D = 5;
681 UV Dui = (UV) D;
682 while (1) {
683 UV gcd = mpz_gcd_ui(NULL, n, Dui);
684 if ((gcd > 1) && mpz_cmp_ui(n, gcd) != 0)
685 return 0;
686 mpz_set_si(t, D);
687 if (mpz_jacobi(t, n) == -1)
688 break;
689 if (Dui == 21 && mpz_perfect_square_p(n))
690 return 0;
691 Dui += 2;
692 D = (D > 0) ? -Dui : Dui;
693 if (Dui > 1000000)
694 croak("lucas_selfridge_params: D exceeded 1e6");
695 }
696 if (P) *P = 1;
697 if (Q) *Q = (1 - D) / 4;
698 return 1;
699 }
700
lucas_extrastrong_params(IV * P,IV * Q,mpz_t n,mpz_t t,UV inc)701 static int lucas_extrastrong_params(IV* P, IV* Q, mpz_t n, mpz_t t, UV inc)
702 {
703 UV tP = 3;
704 if (inc < 1 || inc > 256)
705 croak("Invalid lucas parameter increment: %"UVuf"\n", inc);
706 while (1) {
707 UV D = tP*tP - 4;
708 UV gcd = mpz_gcd_ui(NULL, n, D);
709 if (gcd > 1 && mpz_cmp_ui(n, gcd) != 0)
710 return 0;
711 mpz_set_ui(t, D);
712 if (mpz_jacobi(t, n) == -1)
713 break;
714 if (tP == (3+20*inc) && mpz_perfect_square_p(n))
715 return 0;
716 tP += inc;
717 if (tP > 65535)
718 croak("lucas_extrastrong_params: P exceeded 65535");
719 }
720 if (P) *P = (IV)tP;
721 if (Q) *Q = 1;
722 return 1;
723 }
724
725
726
727 /* This code was verified against Feitsma's psps-below-2-to-64.txt file.
728 * is_strong_pseudoprime reduced it from 118,968,378 to 31,894,014.
729 * all three variations of the Lucas test reduce it to 0.
730 * The test suite should check that they generate the correct pseudoprimes.
731 *
732 * The standard and strong versions use the method A (Selfridge) parameters,
733 * while the extra strong version uses Baillie's parameters from OEIS A217719.
734 *
735 * Using the strong version, we can implement the strong BPSW test as
736 * specified by Baillie and Wagstaff, 1980, page 1401.
737 *
738 * Testing on my x86_64 machine, the strong Lucas code is over 35% faster than
739 * T.R. Nicely's implementation, and over 40% faster than David Cleaver's.
740 */
_GMP_is_lucas_pseudoprime(mpz_t n,int strength)741 int _GMP_is_lucas_pseudoprime(mpz_t n, int strength)
742 {
743 mpz_t d, U, V, Qk, t;
744 IV P, Q;
745 UV s = 0;
746 int rval;
747 int _verbose = get_verbose_level();
748
749 {
750 int cmpr = mpz_cmp_ui(n, 2);
751 if (cmpr == 0) return 1; /* 2 is prime */
752 if (cmpr < 0) return 0; /* below 2 is composite */
753 if (mpz_even_p(n)) return 0; /* multiple of 2 is composite */
754 }
755
756 mpz_init(t);
757 rval = (strength < 2) ? lucas_selfridge_params(&P, &Q, n, t)
758 : lucas_extrastrong_params(&P, &Q, n, t, 1);
759 if (!rval) {
760 mpz_clear(t);
761 return 0;
762 }
763 if (_verbose>3) gmp_printf("N: %Zd D: %"IVdf" P: %"UVuf" Q: %"IVdf"\n", n, P*P-4*Q, P, Q);
764
765 mpz_init(U); mpz_init(V); mpz_init(Qk);
766 mpz_init_set(d, n);
767 mpz_add_ui(d, d, 1);
768
769 if (strength > 0) {
770 s = mpz_scan1(d, 0);
771 mpz_tdiv_q_2exp(d, d, s);
772 }
773
774 lucas_seq(U, V, n, P, Q, d, Qk, t);
775 mpz_clear(d);
776
777 rval = 0;
778 if (strength == 0) {
779 /* Standard checks U_{n+1} = 0 mod n. */
780 rval = (mpz_sgn(U) == 0);
781 } else if (strength == 1) {
782 if (mpz_sgn(U) == 0) {
783 rval = 1;
784 } else {
785 while (s--) {
786 if (mpz_sgn(V) == 0) {
787 rval = 1;
788 break;
789 }
790 if (s) {
791 mpz_mul(V, V, V);
792 mpz_submul_ui(V, Qk, 2);
793 mpz_mod(V, V, n);
794 mpz_mulmod(Qk, Qk, Qk, n, t);
795 }
796 }
797 }
798 } else {
799 mpz_sub_ui(t, n, 2);
800 if ( mpz_sgn(U) == 0 && (mpz_cmp_ui(V, 2) == 0 || mpz_cmp(V, t) == 0) ) {
801 rval = 1;
802 } else {
803 s--; /* The extra strong test tests r < s-1 instead of r < s */
804 while (s--) {
805 if (mpz_sgn(V) == 0) {
806 rval = 1;
807 break;
808 }
809 if (s) {
810 mpz_mul(V, V, V);
811 mpz_sub_ui(V, V, 2);
812 mpz_mod(V, V, n);
813 }
814 }
815 }
816 }
817 mpz_clear(Qk); mpz_clear(V); mpz_clear(U); mpz_clear(t);
818 return rval;
819 }
820
821 /* Pari's clever method. It's an extra-strong Lucas test, but without
822 * computing U_d. This makes it faster, but yields more pseudoprimes.
823 *
824 * increment: 1 for Baillie OEIS, 2 for Pari.
825 *
826 * With increment = 1, these results will be a subset of the extra-strong
827 * Lucas pseudoprimes. With increment = 2, we produce Pari's results (we've
828 * added the necessary GCD with D so we produce somewhat fewer).
829 */
_GMP_is_almost_extra_strong_lucas_pseudoprime(mpz_t n,UV increment)830 int _GMP_is_almost_extra_strong_lucas_pseudoprime(mpz_t n, UV increment)
831 {
832 mpz_t d, V, W, t;
833 UV P, s;
834 int rval;
835
836 {
837 int cmpr = mpz_cmp_ui(n, 2);
838 if (cmpr == 0) return 1; /* 2 is prime */
839 if (cmpr < 0) return 0; /* below 2 is composite */
840 if (mpz_even_p(n)) return 0; /* multiple of 2 is composite */
841 }
842
843 mpz_init(t);
844 {
845 IV PP;
846 if (! lucas_extrastrong_params(&PP, 0, n, t, increment) ) {
847 mpz_clear(t);
848 return 0;
849 }
850 P = (UV) PP;
851 }
852
853 mpz_init(d);
854 mpz_add_ui(d, n, 1);
855
856 s = mpz_scan1(d, 0);
857 mpz_tdiv_q_2exp(d, d, s);
858
859 /* Calculate V_d */
860 {
861 UV b = mpz_sizeinbase(d, 2);
862 mpz_init_set_ui(V, P);
863 mpz_init_set_ui(W, P*P-2); /* V = V_{k}, W = V_{k+1} */
864
865 while (b > 1) {
866 b--;
867 if (mpz_tstbit(d, b-1)) {
868 mpz_mul(V, V, W);
869 mpz_sub_ui(V, V, P);
870
871 mpz_mul(W, W, W);
872 mpz_sub_ui(W, W, 2);
873 } else {
874 mpz_mul(W, V, W);
875 mpz_sub_ui(W, W, P);
876
877 mpz_mul(V, V, V);
878 mpz_sub_ui(V, V, 2);
879 }
880 mpz_mod(V, V, n);
881 mpz_mod(W, W, n);
882 }
883 }
884 mpz_clear(d);
885
886 rval = 0;
887 mpz_sub_ui(t, n, 2);
888
889 #if 0
890 /* According to Stan Wagon's 1999 "Mathematica in Action", this method,
891 * equivalent to our Extra Strong test, is used, though with a different
892 * (unspecified) parameter selection method. I'm not convinced Mathematica
893 * actually uses this method, without any confirmation from Wolfram.
894 */
895 {
896 /* V must be all 2 or x,x,-2,2,2,2. First V=+/-2 must have a W=P. */
897 int must_have_2 = 0, bad = 0;
898
899 if (mpz_cmp_ui(V,2) == 0) {
900 must_have_2 = 1;
901 if (mpz_cmp_ui(W, P) != 0)
902 bad = 1;
903 }
904 while (s-- && !bad) {
905 if (must_have_2) {
906 if (mpz_cmp_ui(V,2) != 0)
907 bad = 1;
908 } else {
909 if (mpz_cmp(V,t) == 0) { /* Found -2. Check W=-P. Look for 2's. */
910 mpz_sub(W, n, W);
911 if (mpz_cmp_ui(W, P) != 0)
912 bad = 1;
913 must_have_2 = 1;
914 } else { /* Keep looking for a -2. */
915 mpz_mul(W, V, W);
916 mpz_sub_ui(W, W, P);
917 mpz_mod(W, W, n);
918 }
919 }
920 mpz_mul(V, V, V);
921 mpz_sub_ui(V, V, 2);
922 mpz_mod(V, V, n);
923 }
924 /* Fail if bad val found or didn't have a 2 at the end */
925 rval = !bad && must_have_2 && !mpz_cmp_ui(V,2);
926 }
927 #else
928 if ( mpz_cmp_ui(V, 2) == 0 || mpz_cmp(V, t) == 0 ) {
929 rval = 1;
930 } else {
931 s--; /* The extra strong test tests r < s-1 instead of r < s */
932 while (s--) {
933 if (mpz_sgn(V) == 0) {
934 rval = 1;
935 break;
936 }
937 if (s) {
938 mpz_mul(V, V, V);
939 mpz_sub_ui(V, V, 2);
940 mpz_mod(V, V, n);
941 }
942 }
943 }
944 #endif
945 mpz_clear(W); mpz_clear(V); mpz_clear(t);
946 return rval;
947 }
948
is_perrin_pseudoprime(mpz_t n,int restricted)949 int is_perrin_pseudoprime(mpz_t n, int restricted)
950 {
951 mpz_t S[6], T[6], T01, T34, T45, t;
952 int cmpr, i, j, rval;
953
954 cmpr = mpz_cmp_ui(n, 2);
955 if (cmpr == 0) return 1; /* 2 is prime */
956 if (cmpr < 0) return 0; /* below 2 is composite */
957 if (restricted > 2 && mpz_even_p(n)) return 0;
958
959 { /* Simple filter for composites */
960 uint32_t n32 = mpz_fdiv_ui(n, 2762760);
961 if (!(n32& 1) && !(( 22 >> (n32% 7)) & 1)) return 0;
962 if (!(n32% 3) && !(( 523 >> (n32%13)) & 1)) return 0;
963 if (!(n32% 5) && !((65890 >> (n32%24)) & 1)) return 0;
964 if (!(n32% 4) && !(( 514 >> (n32%14)) & 1)) return 0;
965 if (!(n32%23) && !(( 2 >> (n32%22)) & 1)) return 0;
966 }
967
968 /* Calculate signature using Adams/Shanks doubling rule. */
969 mpz_init(t);
970 mpz_init_set_ui(S[0], 1);
971 mpz_init(S[1]); mpz_sub_ui(S[1], n, 1);
972 mpz_init_set_ui(S[2], 3);
973 mpz_init_set_ui(S[3], 3);
974 mpz_init_set_ui(S[4], 0);
975 mpz_init_set_ui(S[5], 2);
976
977 for (i=0; i < 6; i++)
978 mpz_init(T[i]);
979 mpz_init(T01); mpz_init(T34); mpz_init(T45);
980
981 for (i = mpz_sizeinbase(n,2)-2; i >= 0; i--) {
982 mpz_mul(t,S[0],S[0]); mpz_sub(t,t,S[5]); mpz_sub(t,t,S[5]); mpz_mod(T[0],t,n);
983 mpz_mul(t,S[1],S[1]); mpz_sub(t,t,S[4]); mpz_sub(t,t,S[4]); mpz_mod(T[1],t,n);
984 mpz_mul(t,S[2],S[2]); mpz_sub(t,t,S[3]); mpz_sub(t,t,S[3]); mpz_mod(T[2],t,n);
985 mpz_mul(t,S[3],S[3]); mpz_sub(t,t,S[2]); mpz_sub(t,t,S[2]); mpz_mod(T[3],t,n);
986 mpz_mul(t,S[4],S[4]); mpz_sub(t,t,S[1]); mpz_sub(t,t,S[1]); mpz_mod(T[4],t,n);
987 mpz_mul(t,S[5],S[5]); mpz_sub(t,t,S[0]); mpz_sub(t,t,S[0]); mpz_mod(T[5],t,n);
988 mpz_sub(t,T[2],T[1]); mpz_mod(T01,t,n);
989 mpz_sub(t,T[5],T[4]); mpz_mod(T34,t,n);
990 mpz_add(t,T34, T[3]); mpz_mod(T45,t,n);
991 if (mpz_tstbit(n, i)) {
992 mpz_set(S[0],T[0]); mpz_set(S[1],T01); mpz_set(S[2],T[1]);
993 mpz_set(S[3],T[4]); mpz_set(S[4],T45); mpz_set(S[5],T[5]);
994 } else {
995 mpz_add(t,T01,T[0]);
996 mpz_set(S[0],T01); mpz_set(S[1],T[1]); mpz_mod(S[2],t,n);
997 mpz_set(S[3],T34); mpz_set(S[4],T[4]); mpz_set(S[5],T45);
998 }
999 }
1000
1001 for (i=0; i < 6; i++)
1002 mpz_clear(T[i]);
1003 mpz_clear(T01); mpz_clear(T34); mpz_clear(T45);
1004
1005 rval = !mpz_sgn(S[4]);
1006 if (rval == 0 || restricted == 0)
1007 goto DONE_PERRIN;
1008
1009 mpz_sub_ui(t,n,1);
1010 rval = !mpz_cmp(S[1],t);
1011 if (rval == 0 || restricted == 1)
1012 goto DONE_PERRIN;
1013
1014 /* Adams/Shanks or Arno,Grantham full signature test */
1015 rval = 0;
1016 j = mpz_si_kronecker(-23, n);
1017
1018 if (j == -1) {
1019 mpz_t A, B, C;
1020 mpz_init_set(B, S[2]); mpz_init(A); mpz_init(C);
1021
1022 mpz_mul(t,B,B); mpz_mod(t,t,n);
1023 mpz_mul_ui(A,B,3); mpz_add_ui(A,A,1); mpz_sub(A,A,t); mpz_mod(A,A,n);
1024 mpz_mul_ui(C,t,3); mpz_sub_ui(C,C,2); mpz_mod(C,C,n);
1025
1026 mpz_mul(t,t,B); mpz_sub(t,t,B); mpz_mod(t,t,n);
1027 rval = !mpz_cmp(S[0],A) && !mpz_cmp(S[2],B) && !mpz_cmp(S[3],B) && !mpz_cmp(S[5],C) && mpz_cmp_ui(B,3) && !mpz_cmp_ui(t,1);
1028 } else if (restricted > 2 && j == 0 && mpz_cmp_ui(n,23)) {
1029 /* Jacobi symbol says 23|n. n is composite if != 23 */
1030 rval = 0;
1031 } else {
1032 if (!mpz_cmp(S[2],S[3])) {
1033 rval = !mpz_cmp_ui(S[0],1) && !mpz_cmp_ui(S[2],3) && !mpz_cmp_ui(S[3],3) && !mpz_cmp_ui(S[5],2);
1034 } else {
1035 mpz_sub_ui(t, n, 1);
1036 rval = !mpz_cmp_ui(S[0],0) && !mpz_cmp(S[5],t) &&
1037 (mpz_add(t,S[2],S[3]),mpz_add_ui(t,t,3),mpz_mod(t,t,n),!mpz_sgn(t)) &&
1038 (mpz_sub(t,S[2],S[3]),mpz_mul(t,t,t),mpz_add_ui(t,t,23),mpz_mod(t,t,n),!mpz_sgn(t));
1039 }
1040 }
1041
1042 DONE_PERRIN:
1043 for (i=0; i < 6; i++)
1044 mpz_clear(S[i]);
1045 mpz_clear(t);
1046 return rval;
1047 }
1048
is_frobenius_pseudoprime(mpz_t n,IV P,IV Q)1049 int is_frobenius_pseudoprime(mpz_t n, IV P, IV Q)
1050 {
1051 mpz_t t, Vcomp, d, U, V, Qk;
1052 IV D;
1053 int k = 0;
1054 int rval;
1055
1056 {
1057 int cmpr = mpz_cmp_ui(n, 2);
1058 if (cmpr == 0) return 1; /* 2 is prime */
1059 if (cmpr < 0) return 0; /* below 2 is composite */
1060 if (mpz_even_p(n)) return 0; /* multiple of 2 is composite */
1061 }
1062 mpz_init(t);
1063 if (P == 0 && Q == 0) {
1064 P = 1; Q = 2;
1065 do {
1066 P += 2;
1067 if (P == 3) P = 5; /* P=3,Q=2 -> D=9-8=1 => k=1, so skip */
1068 if (P == 21 && mpz_perfect_square_p(n))
1069 { mpz_clear(t); return 0; }
1070 D = P*P-4*Q;
1071 if (mpz_cmp_ui(n, P >= 0 ? P : -P) <= 0) break;
1072 if (mpz_cmp_ui(n, D >= 0 ? D : -D) <= 0) break;
1073 mpz_set_si(t, D);
1074 k = mpz_jacobi(t, n);
1075 } while (k == 1);
1076 } else {
1077 D = P*P-4*Q;
1078 if (is_perfect_square( D >= 0 ? D : -D, 0 ))
1079 croak("Frobenius invalid P,Q: (%"IVdf",%"IVdf")", P, Q);
1080 mpz_set_si(t, D);
1081 k = mpz_jacobi(t, n);
1082 }
1083
1084 /* Check initial conditions */
1085 {
1086 UV Pu = P >= 0 ? P : -P;
1087 UV Qu = Q >= 0 ? Q : -Q;
1088 UV Du = D >= 0 ? D : -D;
1089
1090 /* If abs(P) or abs(Q) or abs(D) >= n, exit early. */
1091 if (mpz_cmp_ui(n, Pu) <= 0 || mpz_cmp_ui(n, Qu) <= 0 || mpz_cmp_ui(n, Du) <= 0) {
1092 mpz_clear(t);
1093 return _GMP_trial_factor(n, 2, Du+Pu+Qu) ? 0 : 1;
1094 }
1095 /* If k = 0, then D divides n */
1096 if (k == 0) {
1097 mpz_clear(t);
1098 return 0;
1099 }
1100 /* If n is not coprime to P*Q*D then we found a factor */
1101 if (mpz_gcd_ui(NULL, n, Du*Pu*Qu) > 1) {
1102 mpz_clear(t);
1103 return 0;
1104 }
1105 }
1106
1107 mpz_init(Vcomp);
1108 if (k == 1) {
1109 mpz_set_si(Vcomp, 2);
1110 } else {
1111 mpz_set_si(Vcomp, Q);
1112 mpz_mul_ui(Vcomp, Vcomp, 2);
1113 mpz_mod(Vcomp, Vcomp, n);
1114 }
1115
1116 mpz_init(U); mpz_init(V); mpz_init(Qk); mpz_init(d);
1117 if (k == 1) mpz_sub_ui(d, n, 1);
1118 else mpz_add_ui(d, n, 1);
1119
1120 lucas_seq(U, V, n, P, Q, d, Qk, t);
1121 rval = ( mpz_sgn(U) == 0 && mpz_cmp(V, Vcomp) == 0 );
1122
1123 mpz_clear(d); mpz_clear(Qk); mpz_clear(V); mpz_clear(U);
1124 mpz_clear(Vcomp); mpz_clear(t);
1125
1126 return rval;
1127 }
1128
1129 /* Use Crandall/Pomerance, steps from Loebenberger 2008 */
is_frobenius_cp_pseudoprime(mpz_t n,UV ntests)1130 int is_frobenius_cp_pseudoprime(mpz_t n, UV ntests)
1131 {
1132 mpz_t t, a, b, d, w1, wm, wm1, m;
1133 UV tn;
1134 int j;
1135 int result = 1;
1136
1137 if (mpz_cmp_ui(n, 100) < 0) /* tiny n */
1138 return (_GMP_is_prob_prime(n) > 0);
1139 if (mpz_even_p(n))
1140 return 0;
1141
1142 mpz_init(t); mpz_init(a); mpz_init(b); mpz_init(d);
1143 mpz_init(w1); mpz_init(wm); mpz_init(wm1); mpz_init(m);
1144 for (tn = 0; tn < ntests; tn++) {
1145 /* Step 1: choose a and b in 1..n-1 and d=a^2-4b not square and coprime */
1146 do {
1147 mpz_sub_ui(t, n, 1);
1148 mpz_isaac_urandomm(a, t);
1149 mpz_add_ui(a, a, 1);
1150 mpz_isaac_urandomm(b, t);
1151 mpz_add_ui(b, b, 1);
1152 /* Check d and gcd */
1153 mpz_mul(d, a, a);
1154 mpz_mul_ui(t, b, 4);
1155 mpz_sub(d, d, t);
1156 } while (mpz_perfect_square_p(d));
1157 mpz_mul(t, a, b);
1158 mpz_mul(t, t, d);
1159 mpz_gcd(t, t, n);
1160 if (mpz_cmp_ui(t, 1) != 0 && mpz_cmp(t, n) != 0)
1161 { result = 0; break; }
1162 /* Step 2: W1 = a^2b^-1 - 2 mod n */
1163 if (!mpz_invert(t, b, n))
1164 { result = 0; break; }
1165 mpz_mul(w1, a, a);
1166 mpz_mul(w1, w1, t);
1167 mpz_sub_ui(w1, w1, 2);
1168 mpz_mod(w1, w1, n);
1169 /* Step 3: m = (n-(d|n))/2 */
1170 j = mpz_jacobi(d, n);
1171 if (j == -1) mpz_add_ui(m, n, 1);
1172 else if (j == 0) mpz_set(m, n);
1173 else if (j == 1) mpz_sub_ui(m, n, 1);
1174 mpz_tdiv_q_2exp(m, m, 1);
1175 /* Step 8 here: B = b^(n-1)/2 mod n (stored in d) */
1176 mpz_sub_ui(t, n, 1);
1177 mpz_tdiv_q_2exp(t, t, 1);
1178 mpz_powm(d, b, t, n);
1179 /* Quick Euler test */
1180 mpz_sub_ui(t, n, 1);
1181 if (mpz_cmp_ui(d, 1) != 0 && mpz_cmp(d, t) != 0)
1182 { result = 0; break; }
1183 /* Step 4: calculate Wm,Wm+1 */
1184 mpz_set_ui(wm, 2);
1185 mpz_set(wm1, w1);
1186 {
1187 UV bit = mpz_sizeinbase(m, 2);
1188 while (bit-- > 0) {
1189 if (mpz_tstbit(m, bit)) {
1190 mpz_mul(t, wm, wm1);
1191 mpz_sub(wm, t, w1);
1192 mpz_mul(t, wm1, wm1);
1193 mpz_sub_ui(wm1, t, 2);
1194 } else {
1195 mpz_mul(t, wm, wm1);
1196 mpz_sub(wm1, t, w1);
1197 mpz_mul(t, wm, wm);
1198 mpz_sub_ui(wm, t, 2);
1199 }
1200 mpz_mod(wm, wm, n);
1201 mpz_mod(wm1, wm1, n);
1202 }
1203 }
1204 /* Step 5-7: compare w1/wm */
1205 mpz_mul(t, w1, wm);
1206 mpz_mod(t, t, n);
1207 mpz_mul_ui(wm1, wm1, 2);
1208 mpz_mod(wm1, wm1, n);
1209 if (mpz_cmp(t, wm1) != 0)
1210 { result = 0; break; }
1211 /* Step 8 was earlier */
1212 /* Step 9: See if Bwm = 2 mod n */
1213 mpz_mul(wm, wm, d);
1214 mpz_mod(wm, wm, n);
1215 if (mpz_cmp_ui(wm, 2) != 0)
1216 { result = 0; break; }
1217 }
1218 mpz_clear(w1); mpz_clear(wm); mpz_clear(wm1); mpz_clear(m);
1219 mpz_clear(t); mpz_clear(a); mpz_clear(b); mpz_clear(d);
1220 return result;
1221 }
1222
1223 /* New code based on draft paper */
_GMP_is_frobenius_underwood_pseudoprime(mpz_t n)1224 int _GMP_is_frobenius_underwood_pseudoprime(mpz_t n)
1225 {
1226 mpz_t temp1, temp2, n_plus_1, s, t;
1227 unsigned long a, ap2, len;
1228 int bit, j, rval = 0;
1229 int _verbose = get_verbose_level();
1230
1231 {
1232 int cmpr = mpz_cmp_ui(n, 2);
1233 if (cmpr == 0) return 1; /* 2 is prime */
1234 if (cmpr < 0) return 0; /* below 2 is composite */
1235 if (mpz_even_p(n)) return 0; /* multiple of 2 is composite */
1236 }
1237
1238 mpz_init(temp1);
1239
1240 for (a = 0; a < 1000000; a++) {
1241 if (a==2 || a==4 || a==7 || a==8 || a==10 || a==14 || a==16 || a==18)
1242 continue;
1243 mpz_set_si(temp1, (long)(a*a) - 4);
1244 j = mpz_jacobi(temp1, n);
1245 if (j == -1) break;
1246 if (j == 0 || (a == 20 && mpz_perfect_square_p(n)))
1247 { mpz_clear(temp1); return 0; }
1248 }
1249 if (a >= 1000000)
1250 { mpz_clear(temp1); croak("FU test failure, unable to find suitable a"); }
1251 if (mpz_gcd_ui(NULL, n, (a+4)*(2*a+5)) != 1)
1252 { mpz_clear(temp1); return 0; }
1253
1254 mpz_init(temp2); mpz_init(n_plus_1),
1255
1256 ap2 = a+2;
1257 mpz_add_ui(n_plus_1, n, 1);
1258 len = mpz_sizeinbase(n_plus_1, 2);
1259 mpz_init_set_ui(s, 1);
1260 mpz_init_set_ui(t, 2);
1261
1262 for (bit = len-2; bit >= 0; bit--) {
1263 mpz_add(temp2, t, t);
1264 if (a != 0) {
1265 mpz_mul_ui(temp1, s, a);
1266 mpz_add(temp2, temp1, temp2);
1267 }
1268 mpz_mul(temp1, temp2, s);
1269 mpz_sub(temp2, t, s);
1270 mpz_add(s, s, t);
1271 mpz_mul(t, s, temp2);
1272 mpz_mod(t, t, n);
1273 mpz_mod(s, temp1, n);
1274 if ( mpz_tstbit(n_plus_1, bit) ) {
1275 if (a == 0) mpz_add(temp1, s, s);
1276 else mpz_mul_ui(temp1, s, ap2);
1277 mpz_add(temp1, temp1, t);
1278 mpz_add(temp2, t, t);
1279 mpz_sub(t, temp2, s);
1280 mpz_set(s, temp1);
1281 }
1282 }
1283 /* n+1 always has an even last bit, so s and t always modded */
1284 mpz_set_ui(temp1, 2*a+5);
1285 mpz_mod(temp1, temp1, n);
1286 if (mpz_cmp_ui(s, 0) == 0 && mpz_cmp(t, temp1) == 0)
1287 rval = 1;
1288 if (_verbose>1) { gmp_printf("%Zd is %s with a = %"UVuf"\n", n, (rval) ? "probably prime" : "composite", a); fflush(stdout); }
1289
1290 mpz_clear(temp1); mpz_clear(temp2); mpz_clear(n_plus_1);
1291 mpz_clear(s); mpz_clear(t);
1292 return rval;
1293 }
1294
_GMP_is_frobenius_khashin_pseudoprime(mpz_t n)1295 int _GMP_is_frobenius_khashin_pseudoprime(mpz_t n)
1296 {
1297 mpz_t t, ta, tb, ra, rb, a, b, n_minus_1;
1298 unsigned long c = 1;
1299 int bit, len, k, rval = 0;
1300
1301 {
1302 int cmpr = mpz_cmp_ui(n, 2);
1303 if (cmpr == 0) return 1; /* 2 is prime */
1304 if (cmpr < 0) return 0; /* below 2 is composite */
1305 if (mpz_even_p(n)) return 0; /* multiple of 2 is composite */
1306 }
1307 if (mpz_perfect_square_p(n)) return 0;
1308
1309 mpz_init(t);
1310 do {
1311 c += 2;
1312 mpz_set_ui(t, c);
1313 k = mpz_jacobi(t, n);
1314 } while (k == 1);
1315 if (k == 0) {
1316 mpz_clear(t);
1317 return 0;
1318 }
1319
1320 mpz_init_set_ui(ra, 1); mpz_init_set_ui(rb, 1);
1321 mpz_init_set_ui(a, 1); mpz_init_set_ui(b, 1);
1322 mpz_init(ta); mpz_init(tb);
1323 mpz_init(n_minus_1);
1324 mpz_sub_ui(n_minus_1, n, 1);
1325
1326 len = mpz_sizeinbase(n_minus_1, 2);
1327 for (bit = 0; bit < len; bit++) {
1328 if ( mpz_tstbit(n_minus_1, bit) ) {
1329 mpz_mul(ta, ra, a);
1330 mpz_mul(tb, rb, b);
1331 mpz_add(t, ra, rb);
1332 mpz_add(rb, a, b);
1333 mpz_mul(rb, rb, t);
1334 mpz_sub(rb, rb, ta);
1335 mpz_sub(rb, rb, tb);
1336 mpz_mod(rb, rb, n);
1337 mpz_mul_ui(tb, tb, c);
1338 mpz_add(ra, ta, tb);
1339 mpz_mod(ra, ra, n);
1340 }
1341 if (bit < len-1) {
1342 mpz_mul(t, b, b);
1343 mpz_mul_ui(t, t, c);
1344 mpz_mul(b, b, a);
1345 mpz_add(b, b, b);
1346 mpz_mod(b, b, n);
1347 mpz_mul(a, a, a);
1348 mpz_add(a, a, t);
1349 mpz_mod(a, a, n);
1350 }
1351 }
1352 if ( (mpz_cmp_ui(ra,1) == 0) && (mpz_cmp(rb, n_minus_1) == 0) )
1353 rval = 1;
1354
1355 mpz_clear(n_minus_1);
1356 mpz_clear(ta); mpz_clear(tb);
1357 mpz_clear(a); mpz_clear(b);
1358 mpz_clear(ra); mpz_clear(rb);
1359 mpz_clear(t);
1360 return rval;
1361 }
1362
1363
1364
1365
_GMP_BPSW(mpz_t n)1366 int _GMP_BPSW(mpz_t n)
1367 {
1368 if (mpz_cmp_ui(n, 4) < 0)
1369 return (mpz_cmp_ui(n, 1) <= 0) ? 0 : 2;
1370
1371 if (miller_rabin_ui(n, 2) == 0) /* Miller Rabin with base 2 */
1372 return 0;
1373
1374 if (_GMP_is_lucas_pseudoprime(n, 2 /*extra strong*/) == 0)
1375 return 0;
1376
1377 if (mpz_sizeinbase(n, 2) <= 64) /* BPSW is deterministic below 2^64 */
1378 return 2;
1379
1380 return 1;
1381 }
1382
1383 /* Assume n is a BPSW PRP, return 1 (no result), 0 (composite), 2 (prime) */
is_deterministic_miller_rabin_prime(mpz_t n)1384 int is_deterministic_miller_rabin_prime(mpz_t n)
1385 {
1386 mpz_t t;
1387 int i, res = 1, maxp = 0;
1388
1389 if (mpz_sizeinbase(n, 2) <= 82) {
1390 mpz_init(t);
1391 /* n < 3825123056546413051 => maxp=9, but BPSW should have handled */
1392 if (mpz_set_str(t, "318665857834031151167461",10), mpz_cmp(n,t) < 0)
1393 maxp = 12;
1394 else if (mpz_set_str(t,"3317044064679887385961981",10), mpz_cmp(n,t) < 0)
1395 maxp = 13;
1396 if (maxp > 0) {
1397 for (i = 1; i < maxp && res; i++) {
1398 res = miller_rabin_ui(n, sprimes[i]);
1399 }
1400 if (res == 1) res = 2;
1401 }
1402 mpz_clear(t);
1403 }
1404 return res;
1405 }
1406
1407
1408 /*
1409 * is_prob_prime BPSW -- fast, no known counterexamples
1410 * is_prime is_prob_prime + a little extra
1411 * is_provable_prime really prove it, which could take a very long time
1412 *
1413 * They're all identical for numbers <= 2^64.
1414 *
1415 * The extra M-R test in is_prime start actually costing something after
1416 * 1000 bits or so. Primality proving will get quite a bit slower as the
1417 * number of bits increases.
1418 *
1419 * Both is_prime and is_provable prime start with some trial division
1420 * followed by an extra strong BPSW test, just like is_prob_prime. For
1421 * 65+ bit inputs, they do more. A single extra random-base M-R test is
1422 * next done. In the worst case this has a 1/4 chance of rejecting the
1423 * composite, but an average change of < 1e-12. is_prime will return at
1424 * this point. is_provable_prime tries various special proofs, followed
1425 * by ECPP. ECPP returns either:
1426 * 2 "definitely prime",
1427 * 0 "definitely composite, or
1428 * 1 "probably prime"
1429 * where the latter doesn't really tell us anything. It's unusual, but
1430 * possible. If this even happs, we do a couple Frobenius-type tests, so
1431 * even an answer of 1 (not proven) will have less chance of false results.
1432 * A composite would have to have no small factors, be the first known
1433 * counterexample to each of three different and distinct tests, pass a
1434 * random-base M-R test, and manage to have the ECPP implementation not
1435 * find it a composite.
1436 */
1437
_GMP_is_prob_prime(mpz_t n)1438 int _GMP_is_prob_prime(mpz_t n)
1439 {
1440 /* Step 1: Look for small divisors. This is done purely for performance.
1441 * It is *not* a requirement for the BPSW test. */
1442 int res = primality_pretest(n);
1443 if (res != 1) return res;
1444
1445 /* We'd like to do the LLR test here, but it screws with certificates. */
1446
1447 /* Step 2: The BPSW test. spsp base 2 and slpsp. */
1448 return _GMP_BPSW(n);
1449 }
1450
is_bpsw_dmr_prime(mpz_t n)1451 int is_bpsw_dmr_prime(mpz_t n)
1452 {
1453 int prob_prime = _GMP_BPSW(n);
1454 if (prob_prime == 1) {
1455 prob_prime = is_deterministic_miller_rabin_prime(n);
1456 if (prob_prime == 0) gmp_printf("\n\n**** BPSW counter-example found? ****\n**** N = %Zd ****\n\n", n);
1457 }
1458 return prob_prime;
1459 }
1460
_GMP_is_prime(mpz_t n)1461 int _GMP_is_prime(mpz_t n)
1462 {
1463 UV nbits;
1464 /* Similar to is_prob_prime, but put LLR before BPSW, then do more testing */
1465
1466 /* First, simple pretesting */
1467 int prob_prime = primality_pretest(n);
1468 if (prob_prime != 1) return prob_prime;
1469
1470 /* If the number is of form N=k*2^n-1 and we have a fast proof, do it. */
1471 prob_prime = llr(n);
1472 if (prob_prime == 0 || prob_prime == 2) return prob_prime;
1473
1474 /* If the number is of form N=k*2^n+1 and we have a fast proof, do it. */
1475 prob_prime = proth(n);
1476 if (prob_prime == 0 || prob_prime == 2) return prob_prime;
1477
1478 /* Start with BPSW */
1479 prob_prime = _GMP_BPSW(n);
1480 nbits = mpz_sizeinbase(n, 2);
1481
1482 /* Use Sorenson/Webster 2015 deterministic M-R if possible */
1483 if (prob_prime == 1) {
1484 prob_prime = is_deterministic_miller_rabin_prime(n);
1485 if (prob_prime == 0) gmp_printf("\n\n**** BPSW counter-example found? ****\n**** N = %Zd ****\n\n", n);
1486 }
1487
1488 /* n has passed the ES BPSW test, making it quite unlikely it is a
1489 * composite (and it cannot be if n < 2^64). */
1490
1491 /* For small numbers, try a quick BLS75 proof. */
1492 if (prob_prime == 1) {
1493 if (is_proth_form(n))
1494 prob_prime = _GMP_primality_bls_nm1(n, 2 /* effort */, 0 /* cert */);
1495 else if (nbits <= 150)
1496 prob_prime = _GMP_primality_bls_nm1(n, 0 /* effort */, 0 /* cert */);
1497 /* We could do far better by calling bls75_hybrid, especially with a
1498 * larger effort. But that takes more time. I've decided it isn't worth
1499 * the extra time here. We'll still try a very quick N-1 proof. */
1500 }
1501
1502 /* If prob_prime is still 1, let's run some extra tests.
1503 * Argument against: Nobody has yet found a BPSW counterexample, so
1504 * this is wasted time. They can make a call to one of
1505 * the other tests themselves if they care.
1506 * Argument for: is_prime() should be as correct as reasonably possible.
1507 * They can call is_prob_prime() to skip extra tests.
1508 *
1509 * Choices include:
1510 *
1511 * - A number of fixed-base Miller-Rabin tests.
1512 * - A number of random-base Miller-Rabin tests.
1513 * - A Frobenius-Underwood test.
1514 * - A Frobenius-Khashin test.
1515 *
1516 * The Miller-Rabin tests have better performance.
1517 *
1518 * We will use random bases to make it more difficult to get a false
1519 * result even if someone has a way to generate BPSW pseudoprimes.
1520 *
1521 * We can estimate probabilities of random odd 'k'-bit composites passing
1522 * 't' random-base Miller-Rabin tests using papers such as Kim and
1523 * Pomerance; Damg??rd, Landrock, and Pomerance; Lichtman and Pomerance.
1524 * We can also perform random sampling to get empirical data, which shows
1525 * a probability of less than 1e-12 at 65 bits, and lowering as the number
1526 * of bits increases. One extra test is all that is necessary.
1527 */
1528
1529 if (prob_prime == 1) {
1530 UV ntests = 1;
1531 prob_prime = miller_rabin_random(n, ntests, 0);
1532 /* prob_prime = _GMP_is_frobenius_underwood_pseudoprime(n); */
1533 /* prob_prime = _GMP_is_frobenius_khashin_pseudoprime(n); */
1534 }
1535
1536 /* We have now done trial division, an SPSP-2 test, an extra-strong Lucas
1537 * test, and a random-base Miller-Rabin test. We have exceeded the testing
1538 * done in most math packages. Any composites will have no small factors,
1539 * be the first known BPSW counterexample, and gotten past the random-base
1540 * Miller-Rabin test.
1541 */
1542
1543 return prob_prime;
1544 }
1545
1546
_GMP_is_provable_prime(mpz_t n,char ** prooftext)1547 int _GMP_is_provable_prime(mpz_t n, char** prooftext)
1548 {
1549 int prob_prime = primality_pretest(n);
1550 if (prob_prime != 1) return prob_prime;
1551
1552 /* Try LLR and Proth if they don't need a proof certificate. */
1553 if (prooftext == 0) {
1554 prob_prime = llr(n);
1555 if (prob_prime == 0 || prob_prime == 2) return prob_prime;
1556 prob_prime = proth(n);
1557 if (prob_prime == 0 || prob_prime == 2) return prob_prime;
1558 }
1559
1560 /* Start with BPSW */
1561 prob_prime = _GMP_BPSW(n);
1562 if (prob_prime != 1) return prob_prime;
1563
1564 /* Use Sorenson/Webster 2015 deterministic M-R if possible */
1565 if (prooftext == 0) {
1566 prob_prime = is_deterministic_miller_rabin_prime(n);
1567 if (prob_prime != 1) return prob_prime;
1568 }
1569
1570 /* Run one more M-R test, just in case. */
1571 prob_prime = miller_rabin_random(n, 1, 0);
1572 if (prob_prime != 1) return prob_prime;
1573
1574 /* We can choose a primality proving algorithm:
1575 * AKS _GMP_is_aks_prime really slow, don't bother
1576 * N-1 _GMP_primality_bls_nm1 small or special numbers
1577 * ECPP _GMP_ecpp fastest in general
1578 */
1579
1580 /* Give n-1 a small go */
1581 prob_prime = _GMP_primality_bls_nm1(n, is_proth_form(n) ? 3 : 1, prooftext);
1582 if (prob_prime != 1) return prob_prime;
1583
1584 /* ECPP */
1585 prob_prime = _GMP_ecpp(n, prooftext);
1586
1587 /* While extremely unusual, it is not impossible for our ECPP implementation
1588 * to give up. If this happens, we won't return 2 with a proof, but let's
1589 * at least run some more tests. */
1590 if (prob_prime == 1)
1591 prob_prime = _GMP_is_frobenius_underwood_pseudoprime(n);
1592 if (prob_prime == 1)
1593 prob_prime = _GMP_is_frobenius_khashin_pseudoprime(n);
1594
1595 return prob_prime;
1596 }
1597