File Coverage

primality.c
Criterion Covered Total %
statement 593 673 88.1
branch 660 992 66.5
condition n/a
subroutine n/a
pod n/a
total 1253 1665 75.2


line stmt bran cond sub pod time code
1             #include
2             #include
3             #include
4             #include
5              
6             #include "ptypes.h"
7             #include "primality.h"
8             #include "mulmod.h"
9             #define FUNC_gcd_ui 1
10             #define FUNC_is_perfect_square
11             #include "util.h"
12             #include "montmath.h" /* Fast Montgomery math */
13              
14             /* Primality related functions */
15              
16             #if !USE_MONTMATH
17             static const UV mr_bases_const2[1] = {2};
18             #endif
19             /******************************************************************************/
20              
21              
22 13123           static int jacobi_iu(IV in, UV m) {
23 13123           int j = 1;
24 13123           UV n = (in < 0) ? -in : in;
25              
26 13123 50         if (m <= 0 || (m%2) == 0) return 0;
    50          
27 13123 100         if (in < 0 && (m%4) == 3) j = -j;
    100          
28 43812 100         while (n != 0) {
29 56988 100         while ((n % 2) == 0) {
30 26299           n >>= 1;
31 26299 100         if ( (m % 8) == 3 || (m % 8) == 5 ) j = -j;
    100          
32             }
33 30689           { UV t = n; n = m; m = t; }
34 30689 100         if ( (n % 4) == 3 && (m % 4) == 3 ) j = -j;
    100          
35 30689           n = n % m;
36             }
37 13123 100         return (m == 1) ? j : 0;
38             }
39              
40 5731           static UV select_extra_strong_parameters(UV n, UV increment) {
41             int j;
42 5731           UV D, P = 3;
43             while (1) {
44 12858           D = P*P - 4;
45 12858           j = jacobi_iu(D, n);
46 12858 50         if (j == 0) return 0;
47 12858 100         if (j == -1) break;
48 7127 100         if (P == (3+20*increment) && is_perfect_square(n)) return 0;
    50          
49 7127           P += increment;
50 7127 50         if (P > 65535)
51 0           croak("lucas_extrastrong_params: P exceeded 65535");
52 7127           }
53 5731 50         if (P >= n) P %= n; /* Never happens with increment < 4 */
54 5731           return P;
55             }
56              
57             /* Fermat pseudoprime */
58 86           int is_pseudoprime(UV const n, UV a)
59             {
60 86 50         if (n < 4) return (n == 2 || n == 3);
    0          
    0          
61 86 100         if (!(n&1) && !(a&1)) return 0;
    50          
62 86 50         if (a < 2) croak("Base %"UVuf" is invalid", a);
63 86 50         if (a >= n) {
64 0           a %= n;
65 0 0         if (a <= 1) return (a == 1);
66 0 0         if (a == n-1) return !(a & 1);
67             }
68              
69             #if USE_MONTMATH
70 86 100         if (n & 1) { /* The Montgomery code only works for odd n */
71 84           const uint64_t npi = mont_inverse(n), mont1 = mont_get1(n);
72 84 100         const uint64_t monta = (a == 2) ? mont_get2(n) : mont_geta(a, n);
73 84           return mont_powmod(monta, n-1, n) == mont1;
74             }
75             #endif
76 2           return powmod(a, n-1, n) == 1; /* a^(n-1) = 1 mod n */
77             }
78              
79             /* Euler (aka Euler-Jacobi) pseudoprime: a^((n-1)/2) = (a|n) mod n */
80 109           int is_euler_pseudoprime(UV const n, UV a)
81             {
82 109 50         if (n < 5) return (n == 2 || n == 3);
    0          
    0          
83 109 50         if (!(n&1)) return 0;
84 109 50         if (a < 2) croak("Base %"UVuf" is invalid", a);
85 109 100         if (a > 2) {
86 73 100         if (a >= n) {
87 1           a %= n;
88 1 50         if (a <= 1) return (a == 1);
89 1 50         if (a == n-1) return !(a & 1);
90             }
91 72 50         if ((n % a) == 0) return 0;
92             }
93             {
94             #if USE_MONTMATH
95 108           const uint64_t npi = mont_inverse(n), mont1 = mont_get1(n);
96 108           const uint64_t monta = mont_geta(a, n);
97 108           UV ap = mont_powmod(monta, (n-1)>>1, n);
98 108 100         if (ap != mont1 && ap != n-mont1) return 0;
    50          
99 108 100         if (a == 2) {
100 36           uint32_t nmod8 = n & 0x7;
101 36 100         return (nmod8 == 1 || nmod8 == 7) ? (ap == mont1) : (ap == n-mont1);
    100          
102             } else {
103 72 100         return (kronecker_uu(a,n) >= 0) ? (ap == mont1) : (ap == n-mont1);
104             }
105             #else
106             UV ap = powmod(a, (n-1)>>1, n);
107             if (ap != 1 && ap != n-1) return 0;
108             if (a == 2) {
109             uint32_t nmod8 = n & 0x7;
110             return (nmod8 == 1 || nmod8 == 7) ? (ap == 1) : (ap == n-1);
111             } else {
112             return (kronecker_uu(a,n) >= 0) ? (ap == 1) : (ap == n-1);
113             }
114             #endif
115             }
116             }
117              
118             /* Colin Plumb's extended Euler Criterion test.
119             * A tiny bit (~1 percent) faster than base 2 Fermat or M-R.
120             * More stringent than base 2 Fermat, but a subset of base 2 M-R.
121             */
122 1195           int is_euler_plumb_pseudoprime(UV const n)
123             {
124             UV ap;
125 1195           uint32_t nmod8 = n & 0x7;
126 1195 50         if (n < 5) return (n == 2 || n == 3);
    0          
    0          
127 1195 50         if (!(n&1)) return 0;
128             #if USE_MONTMATH
129             {
130 1195           const uint64_t npi = mont_inverse(n), mont1 = mont_get1(n);
131 1195           const uint64_t mont2 = mont_get2(n);
132 1195 100         ap = mont_powmod(mont2, (n-1) >> (1 + (nmod8 == 1)), n);
133 1195 100         if (ap == mont1) return (nmod8 == 1 || nmod8 == 7);
    100          
    50          
134 770 100         if (ap == n-mont1) return (nmod8 == 1 || nmod8 == 3 || nmod8 == 5);
    100          
    100          
    50          
135             }
136             #else
137             ap = powmod(2, (n-1) >> (1 + (nmod8 == 1)), n);
138             if (ap == 1) return (nmod8 == 1 || nmod8 == 7);
139             if (ap == n-1) return (nmod8 == 1 || nmod8 == 3 || nmod8 == 5);
140             #endif
141 9           return 0;
142             }
143              
144             /* Miller-Rabin probabilistic primality test
145             * Returns 1 if probably prime relative to the bases, 0 if composite.
146             * Bases must be between 2 and n-2
147             */
148 88337           int miller_rabin(UV const n, const UV *bases, int nbases)
149             {
150             #if USE_MONTMATH
151 88337 50         MPUassert(n > 3, "MR called with n <= 3");
152 88337 100         if ((n & 1) == 0) return 0;
153             {
154 88336           const uint64_t npi = mont_inverse(n), mont1 = mont_get1(n);
155 88336           uint64_t a, ma, md, u = n-1;
156 88336           int i, j, t = 0;
157              
158 263346 100         while (!(u&1)) { t++; u >>= 1; }
159 122521 100         for (j = 0; j < nbases; j++) {
160 92079           a = bases[j];
161 92079 100         if (a < 2) croak("Base %"UVuf" is invalid", (UV)a);
162 92077 100         if (a >= n) {
163 76           a %= n;
164 76 100         if (a == 0 || (a == n-1 && a&1)) return 0;
    100          
    50          
165             }
166 92071           ma = mont_geta(a,n);
167 92071 100         if (a == 1 || a == n-1 || !ma) continue;
    100          
    50          
168 92059           md = mont_powmod(ma, u, n);
169 92059 100         if (md != mont1 && md != n-mont1) {
    100          
170 138079 100         for (i=1; i
171 80232 50         md = mont_sqrmod(md, n);
172 80232 100         if (md == mont1) return 0;
173 80193 100         if (md == n-mont1) break;
174             }
175 67796 100         if (i == t)
176 57847           return 0;
177             }
178             }
179             }
180             #else
181             UV d = n-1;
182             int b, r, s = 0;
183              
184             MPUassert(n > 3, "MR called with n <= 3");
185             if ((n & 1) == 0) return 0;
186              
187             while (!(d&1)) { s++; d >>= 1; }
188             for (b = 0; b < nbases; b++) {
189             UV x, a = bases[b];
190             if (a < 2) croak("Base %"UVuf" is invalid", a);
191             if (a >= n) {
192             a %= n;
193             if (a == 0 || (a == n-1 && a&1)) return 0;
194             }
195             if (a == 1 || a == n-1) continue;
196             /* n is a strong pseudoprime to base a if either:
197             * a^d = 1 mod n
198             * a^(d2^r) = -1 mod n for some r: 0 <= r <= s-1
199             */
200             x = powmod(a, d, n);
201             if ( (x == 1) || (x == n-1) ) continue;
202             for (r = 1; r < s; r++) { /* r=0 was just done, test r = 1 to s-1 */
203             x = sqrmod(x, n);
204             if ( x == n-1 ) break;
205             if ( x == 1 ) return 0;
206             }
207             if (r >= s) return 0;
208             }
209             #endif
210 30442           return 1;
211             }
212              
213 8015           int BPSW(UV const n)
214             {
215 8015 50         if (n < 7) return (n == 2 || n == 3 || n == 5);
    0          
    0          
    0          
216 8015 50         if ((n % 2) == 0 || n == UV_MAX) return 0;
    50          
217              
218             #if !USE_MONTMATH
219             return miller_rabin(n, mr_bases_const2, 1)
220             && is_almost_extra_strong_lucas_pseudoprime(n,1);
221             #else
222             {
223 8015           const uint64_t npi = mont_inverse(n), mont1 = mont_get1(n);
224 8015           const uint64_t mont2 = mont_get2(n);
225 8015           uint64_t md, u = n-1;
226 8015           int i, t = 0;
227             UV P, V, d, s;
228              
229             /* M-R with base 2 */
230 24175 100         while (!(u&1)) { t++; u >>= 1; }
231 8015           md = mont_powmod(mont2, u, n);
232 8015 100         if (md != mont1 && md != n-mont1) {
    100          
233 9685 100         for (i=1; i
234 6195 100         md = mont_sqrmod(md, n);
235 6195 50         if (md == mont1) return 0;
236 6195 100         if (md == n-mont1) break;
237             }
238 5337 100         if (i == t)
239 3490           return 0;
240             }
241             /* AES Lucas test */
242 4525           P = select_extra_strong_parameters(n, 1);
243 4525 50         if (P == 0) return 0;
244              
245 4525           d = n+1;
246 4525           s = 0;
247 13755 100         while ( (d & 1) == 0 ) { s++; d >>= 1; }
248              
249             {
250 4525           const uint64_t montP = mont_geta(P, n);
251             UV W, b;
252 4525 100         W = submod( mont_mulmod( montP, montP, n), mont2, n);
253 4525           V = montP;
254 218916 100         { UV v = d; b = 1; while (v >>= 1) b++; }
255 218916 100         while (b-- > 1) {
256 214391 100         UV T = submod( mont_mulmod(V, W, n), montP, n);
257 214391 100         if ( (d >> (b-1)) & UVCONST(1) ) {
258 114665           V = T;
259 114665 100         W = submod( mont_mulmod(W, W, n), mont2, n);
260             } else {
261 99726           W = T;
262 99726 100         V = submod( mont_mulmod(V, V, n), mont2, n);
263             }
264             }
265             }
266              
267 4525 100         if (V == mont2 || V == (n-mont2))
    100          
268 2506           return 1;
269 4452 100         while (s-- > 1) {
270 4432 100         if (V == 0)
271 1999           return 1;
272 2433 100         V = submod( mont_mulmod(V, V, n), mont2, n);
273 2433 50         if (V == mont2)
274 0           return 0;
275             }
276             }
277 20           return 0;
278             #endif
279             }
280              
281             /* Alternate modular lucas sequence code.
282             * A bit slower than the normal one, but works with even valued n. */
283 1           static void alt_lucas_seq(UV* Uret, UV* Vret, UV* Qkret, UV n, UV Pmod, UV Qmod, UV k)
284             {
285             UV Uh, Vl, Vh, Ql, Qh;
286             int j, s, m;
287              
288 1           Uh = 1; Vl = 2; Vh = Pmod; Ql = 1; Qh = 1;
289 1           s = 0; m = 0;
290 1 50         { UV v = k; while (!(v & 1)) { v >>= 1; s++; } }
291 24 100         { UV v = k; while (v >>= 1) m++; }
292              
293 1 50         if (Pmod == 1 && Qmod == (n-1)) {
    50          
294 1           int Sl = Ql, Sh = Qh;
295 24 100         for (j = m; j > s; j--) {
296 23           Sl *= Sh;
297 23 100         Ql = (Sl==1) ? 1 : n-1;
298 23 100         if ( (k >> j) & UVCONST(1) ) {
299 8           Sh = -Sl;
300 8           Uh = mulmod(Uh, Vh, n);
301 8           Vl = submod(mulmod(Vh, Vl, n), Ql, n);
302 8 100         Vh = submod(sqrmod(Vh, n), (Sh==1) ? 2 : n-2, n);
303             } else {
304 15           Sh = Sl;
305 15           Uh = submod(mulmod(Uh, Vl, n), Ql, n);
306 15           Vh = submod(mulmod(Vh, Vl, n), Ql, n);
307 15 100         Vl = submod(sqrmod(Vl, n), (Sl==1) ? 2 : n-2, n);
308             }
309             }
310 1           Sl *= Sh;
311 1 50         Ql = (Sl==1) ? 1 : n-1;
312 1           Uh = submod(mulmod(Uh, Vl, n), Ql, n);
313 1           Vl = submod(mulmod(Vh, Vl, n), Ql, n);
314 1 50         for (j = 0; j < s; j++) {
315 0           Uh = mulmod(Uh, Vl, n);
316 0 0         Vl = submod(sqrmod(Vl, n), (j>0) ? 2 : n-2, n);
317             }
318 1           *Uret = Uh;
319 1           *Vret = Vl;
320 1 50         *Qkret = (s>0)?1:n-1;
321 1           return;
322             }
323              
324 0 0         for (j = m; j > s; j--) {
325 0           Ql = mulmod(Ql, Qh, n);
326 0 0         if ( (k >> j) & UVCONST(1) ) {
327 0           Qh = mulmod(Ql, Qmod, n);
328 0           Uh = mulmod(Uh, Vh, n);
329 0           Vl = submod(mulmod(Vh, Vl, n), mulmod(Pmod, Ql, n), n);
330 0           Vh = submod(sqrmod(Vh, n), mulmod(2, Qh, n), n);
331             } else {
332 0           Qh = Ql;
333 0           Uh = submod(mulmod(Uh, Vl, n), Ql, n);
334 0           Vh = submod(mulmod(Vh, Vl, n), mulmod(Pmod, Ql, n), n);
335 0           Vl = submod(sqrmod(Vl, n), mulmod(2, Ql, n), n);
336             }
337             }
338 0           Ql = mulmod(Ql, Qh, n);
339 0           Qh = mulmod(Ql, Qmod, n);
340 0           Uh = submod(mulmod(Uh, Vl, n), Ql, n);
341 0           Vl = submod(mulmod(Vh, Vl, n), mulmod(Pmod, Ql, n), n);
342 0           Ql = mulmod(Ql, Qh, n);
343 0 0         for (j = 0; j < s; j++) {
344 0           Uh = mulmod(Uh, Vl, n);
345 0           Vl = submod(sqrmod(Vl, n), mulmod(2, Ql, n), n);
346 0           Ql = sqrmod(Ql, n);
347             }
348 0           *Uret = Uh;
349 0           *Vret = Vl;
350 0           *Qkret = Ql;
351             }
352              
353             /* Generic Lucas sequence for any appropriate P and Q */
354 26319           void lucas_seq(UV* Uret, UV* Vret, UV* Qkret, UV n, IV P, IV Q, UV k)
355             {
356             UV U, V, b, Dmod, Qmod, Pmod, Qk;
357              
358 26319 50         MPUassert(n > 1, "lucas_sequence: modulus n must be > 1");
359 26319 100         if (k == 0) {
360 25           *Uret = 0;
361 25           *Vret = 2;
362 25           *Qkret = Q;
363 25           return;
364             }
365              
366 26294 100         Qmod = (Q < 0) ? (UV) (Q + (IV)(((-Q/n)+1)*n)) : (UV)Q % n;
367 26294 50         Pmod = (P < 0) ? (UV) (P + (IV)(((-P/n)+1)*n)) : (UV)P % n;
368 26294           Dmod = submod( mulmod(Pmod, Pmod, n), mulmod(4, Qmod, n), n);
369 26294 100         if (Dmod == 0) {
370 13           b = Pmod >> 1;
371 13           *Uret = mulmod(k, powmod(b, k-1, n), n);
372 13           *Vret = mulmod(2, powmod(b, k, n), n);
373 13           *Qkret = powmod(Qmod, k, n);
374 13           return;
375             }
376 26281 100         if ((n % 2) == 0) {
377 1           alt_lucas_seq(Uret, Vret, Qkret, n, Pmod, Qmod, k);
378 1           return;
379             }
380 26280           U = 1;
381 26280           V = Pmod;
382 26280           Qk = Qmod;
383 391336 100         { UV v = k; b = 0; while (v >>= 1) b++; }
384              
385 26280 100         if (Q == 1) {
386 186 100         while (b--) {
387 133           U = mulmod(U, V, n);
388 133           V = mulsubmod(V, V, 2, n);
389 133 100         if ( (k >> b) & UVCONST(1) ) {
390 51           UV t2 = mulmod(U, Dmod, n);
391 51           U = muladdmod(U, Pmod, V, n);
392 51 100         if (U & 1) { U = (n>>1) + (U>>1) + 1; } else { U >>= 1; }
393 51           V = muladdmod(V, Pmod, t2, n);
394 51 100         if (V & 1) { V = (n>>1) + (V>>1) + 1; } else { V >>= 1; }
395             }
396             }
397 52171 100         } else if (P == 1 && Q == -1) {
    100          
398             /* This is about 30% faster than the generic code below. Since 50% of
399             * Lucas and strong Lucas tests come here, I think it's worth doing. */
400 25944           int sign = Q;
401 389691 100         while (b--) {
402 363747           U = mulmod(U, V, n);
403 363747 100         if (sign == 1) V = mulsubmod(V, V, 2, n);
404 193981           else V = muladdmod(V, V, 2, n);
405 363747           sign = 1; /* Qk *= Qk */
406 363747 100         if ( (k >> b) & UVCONST(1) ) {
407 168054           UV t2 = mulmod(U, Dmod, n);
408 168054           U = addmod(U, V, n);
409 168054 100         if (U & 1) { U = (n>>1) + (U>>1) + 1; } else { U >>= 1; }
410 168054           V = addmod(V, t2, n);
411 168054 100         if (V & 1) { V = (n>>1) + (V>>1) + 1; } else { V >>= 1; }
412 168054           sign = -1; /* Qk *= Q */
413             }
414             }
415 25944 100         if (sign == 1) Qk = 1;
416             } else {
417 1459 100         while (b--) {
418 1176           U = mulmod(U, V, n);
419 1176           V = mulsubmod(V, V, addmod(Qk,Qk,n), n);
420 1176           Qk = sqrmod(Qk, n);
421 1176 100         if ( (k >> b) & UVCONST(1) ) {
422 536           UV t2 = mulmod(U, Dmod, n);
423 536           U = muladdmod(U, Pmod, V, n);
424 536 100         if (U & 1) { U = (n>>1) + (U>>1) + 1; } else { U >>= 1; }
425 536           V = muladdmod(V, Pmod, t2, n);
426 536 100         if (V & 1) { V = (n>>1) + (V>>1) + 1; } else { V >>= 1; }
427 536           Qk = mulmod(Qk, Qmod, n);
428             }
429             }
430             }
431 26280           *Uret = U;
432 26280           *Vret = V;
433 26280           *Qkret = Qk;
434             }
435              
436             #define OVERHALF(v) ( (UV)((v>=0)?v:-v) > (UVCONST(1) << (BITS_PER_WORD/2-1)) )
437 199           int lucasu(IV* U, IV P, IV Q, UV k)
438             {
439             IV Uh, Vl, Vh, Ql, Qh;
440             int j, s, n;
441              
442 199 50         if (U == 0) return 0;
443 199 100         if (k == 0) { *U = 0; return 1; }
444              
445 185           Uh = 1; Vl = 2; Vh = P; Ql = 1; Qh = 1;
446 185           s = 0; n = 0;
447 332 100         { UV v = k; while (!(v & 1)) { v >>= 1; s++; } }
448 588 100         { UV v = k; while (v >>= 1) n++; }
449              
450 441 100         for (j = n; j > s; j--) {
451 256 50         if (OVERHALF(Uh) || OVERHALF(Vh) || OVERHALF(Vl) || OVERHALF(Ql) || OVERHALF(Qh)) return 0;
    50          
    50          
    50          
    50          
452 256           Ql *= Qh;
453 256 100         if ( (k >> j) & UVCONST(1) ) {
454 175           Qh = Ql * Q;
455 175           Uh = Uh * Vh;
456 175           Vl = Vh * Vl - P * Ql;
457 175           Vh = Vh * Vh - 2 * Qh;
458             } else {
459 81           Qh = Ql;
460 81           Uh = Uh * Vl - Ql;
461 81           Vh = Vh * Vl - P * Ql;
462 81           Vl = Vl * Vl - 2 * Ql;
463             }
464             }
465 185 50         if (OVERHALF(Ql) || OVERHALF(Qh)) return 0;
    50          
466 185           Ql = Ql * Qh;
467 185           Qh = Ql * Q;
468 185 50         if (OVERHALF(Uh) || OVERHALF(Vh) || OVERHALF(Vl) || OVERHALF(Ql) || OVERHALF(Qh)) return 0;
    50          
    50          
    50          
    50          
469 185           Uh = Uh * Vl - Ql;
470 185           Vl = Vh * Vl - P * Ql;
471 185           Ql = Ql * Qh;
472 332 100         for (j = 0; j < s; j++) {
473 147 50         if (OVERHALF(Uh) || OVERHALF(Vl) || OVERHALF(Ql)) return 0;
    50          
    50          
474 147           Uh *= Vl;
475 147           Vl = Vl * Vl - 2 * Ql;
476 147           Ql *= Ql;
477             }
478 185           *U = Uh;
479 185           return 1;
480             }
481 152           int lucasv(IV* V, IV P, IV Q, UV k)
482             {
483             IV Vl, Vh, Ql, Qh;
484             int j, s, n;
485              
486 152 50         if (V == 0) return 0;
487 152 100         if (k == 0) { *V = 2; return 1; }
488              
489 141           Vl = 2; Vh = P; Ql = 1; Qh = 1;
490 141           s = 0; n = 0;
491 251 100         { UV v = k; while (!(v & 1)) { v >>= 1; s++; } }
492 445 100         { UV v = k; while (v >>= 1) n++; }
493              
494 335 100         for (j = n; j > s; j--) {
495 194 50         if (OVERHALF(Vh) || OVERHALF(Vl) || OVERHALF(Ql) || OVERHALF(Qh)) return 0;
    50          
    50          
    50          
496 194           Ql *= Qh;
497 194 100         if ( (k >> j) & UVCONST(1) ) {
498 130           Qh = Ql * Q;
499 130           Vl = Vh * Vl - P * Ql;
500 130           Vh = Vh * Vh - 2 * Qh;
501             } else {
502 64           Qh = Ql;
503 64           Vh = Vh * Vl - P * Ql;
504 64           Vl = Vl * Vl - 2 * Ql;
505             }
506             }
507 141 50         if (OVERHALF(Ql) || OVERHALF(Qh)) return 0;
    50          
508 141           Ql = Ql * Qh;
509 141           Qh = Ql * Q;
510 141 50         if (OVERHALF(Vh) || OVERHALF(Vl) || OVERHALF(Ql) || OVERHALF(Qh)) return 0;
    50          
    50          
    50          
511 141           Vl = Vh * Vl - P * Ql;
512 141           Ql = Ql * Qh;
513 251 100         for (j = 0; j < s; j++) {
514 110 50         if (OVERHALF(Vl) || OVERHALF(Ql)) return 0;
    50          
515 110           Vl = Vl * Vl - 2 * Ql;
516 110           Ql *= Ql;
517             }
518 141           *V = Vl;
519 141           return 1;
520             }
521              
522             /* Lucas tests:
523             * 0: Standard
524             * 1: Strong
525             * 2: Stronger (Strong + page 1401 extra tests)
526             * 3: Extra Strong (Mo/Jones/Grantham)
527             *
528             * None of them have any false positives for the BPSW test. Also see the
529             * "almost extra strong" test.
530             */
531 71           int is_lucas_pseudoprime(UV n, int strength)
532             {
533             IV P, Q, D;
534             UV U, V, Qk, d, s;
535              
536 71 100         if (n < 7) return (n == 2 || n == 3 || n == 5);
    50          
    0          
    0          
537 70 100         if ((n % 2) == 0 || n == UV_MAX) return 0;
    50          
538              
539 65 100         if (strength < 3) {
540 43           UV Du = 5;
541 43           IV sign = 1;
542             int j;
543             while (1) {
544 108           D = Du * sign;
545 108           j = jacobi_iu(D, n);
546 108 100         if (j != 1 && Du != n) break;
    100          
547 65 50         if (Du == 21 && is_perfect_square(n)) return 0;
    0          
548 65           Du += 2;
549 65           sign = -sign;
550 65           }
551 43 100         if (j != -1) return 0;
552 42           P = 1;
553 42           Q = (1 - D) / 4;
554 42 50         if (strength == 2 && Q == -1) P=Q=D=5; /* Method A* */
    0          
555             /* Check gcd(n,2QD). gcd(n,2D) already done. */
556 42 100         Qk = (Q >= 0) ? Q % n : n-(((UV)(-Q)) % n);
557 42 50         if (gcd_ui(Qk,n) != 1) return 0;
558             } else {
559 22           P = select_extra_strong_parameters(n, 1);
560 22 50         if (P == 0) return 0;
561 22           Q = 1;
562 22           D = P*P - 4;
563             }
564 64 50         MPUassert( D == (P*P - 4*Q) , "is_lucas_pseudoprime: incorrect DPQ");
565              
566             #if 0 /* Condition 2, V_n+1 = 2Q mod n */
567             { UV us, vs, qs; lucas_seq(&us, &vs, &qs, n, P, Q, n+1); return (vs == addmod(Q,Q,n)); }
568             #endif
569             #if 0 /* Condition 3, n is a epsp(Q) */
570             return is_euler_pseudoprime(n,Qk);
571             #endif
572              
573 64           d = n+1;
574 64           s = 0;
575 64 100         if (strength > 0)
576 142 100         while ( (d & 1) == 0 ) { s++; d >>= 1; }
577              
578             #if USE_MONTMATH
579             {
580 64           const uint64_t npi = mont_inverse(n), mont1 = mont_get1(n);
581 64           const uint64_t mont2 = mont_get2(n);
582 64           const uint64_t montP = (P == 1) ? mont1
583 86 100         : (P >= 0) ? mont_geta(P, n)
584 22 50         : n - mont_geta(-P, n);
585 64           const uint64_t montQ = (Q == 1) ? mont1
586 106 100         : (Q >= 0) ? mont_geta(Q, n)
587 42 100         : n - mont_geta(-Q, n);
588 106           const uint64_t montD = (D >= 0) ? mont_geta(D, n)
589 64 100         : n - mont_geta(-D, n);
590             UV b;
591 999 100         { UV v = d; b = 0; while (v >>= 1) b++; }
592              
593             /* U, V, Qk, and mont* are in Montgomery space */
594 64           U = mont1;
595 64           V = montP;
596              
597 102 100         if (Q == 1 || Q == -1) { /* Faster code for |Q|=1, also opt for P=1 */
    100          
598 38           int sign = Q;
599 610 100         while (b--) {
600 572 50         U = mont_mulmod(U, V, n);
601 572 100         if (sign == 1) V = submod( mont_sqrmod(V,n), mont2, n);
    50          
602 101 50         else V = addmod( mont_sqrmod(V,n), mont2, n);
603 572           sign = 1;
604 572 100         if ( (d >> b) & UVCONST(1) ) {
605 238 50         UV t2 = mont_mulmod(U, montD, n);
606 238 100         if (P == 1) {
607 92           U = addmod(U, V, n);
608 92           V = addmod(V, t2, n);
609             } else {
610 146 50         U = addmod( mont_mulmod(U, montP, n), V, n);
611 146 50         V = addmod( mont_mulmod(V, montP, n), t2, n);
612             }
613 238 100         if (U & 1) { U = (n>>1) + (U>>1) + 1; } else { U >>= 1; }
614 238 100         if (V & 1) { V = (n>>1) + (V>>1) + 1; } else { V >>= 1; }
615 238           sign = Q;
616             }
617             }
618 38 100         Qk = (sign == 1) ? mont1 : n-mont1;
619             } else {
620 26           Qk = montQ;
621 389 100         while (b--) {
622 363 50         U = mont_mulmod(U, V, n);
623 363 50         V = submod( mont_sqrmod(V,n), addmod(Qk,Qk,n), n);
624 363 50         Qk = mont_sqrmod(Qk,n);
625 363 100         if ( (d >> b) & UVCONST(1) ) {
626 186 50         UV t2 = mont_mulmod(U, montD, n);
627 186 50         U = addmod( mont_mulmod(U, montP, n), V, n);
628 186 100         if (U & 1) { U = (n>>1) + (U>>1) + 1; } else { U >>= 1; }
629 186 50         V = addmod( mont_mulmod(V, montP, n), t2, n);
630 186 100         if (V & 1) { V = (n>>1) + (V>>1) + 1; } else { V >>= 1; }
631 186 50         Qk = mont_mulmod(Qk, montQ, n);
632             }
633             }
634             }
635              
636 64 100         if (strength == 0) {
637 20 50         if (U == 0)
638 20           return 1;
639 44 100         } else if (strength == 1) {
640 22 100         if (U == 0)
641 8           return 1;
642 39 100         while (s--) {
643 36 100         if (V == 0)
644 11           return 1;
645 25 100         if (s) {
646 22 50         V = submod( mont_sqrmod(V,n), addmod(Qk,Qk,n), n);
647 22 50         Qk = mont_sqrmod(Qk,n);
648             }
649             }
650 22 50         } else if (strength == 2) {
651 0           UV Ql = 0, Qj = 0;
652 0           int qjacobi, is_slpsp = 0;
653 0 0         if (U == 0)
654 0           is_slpsp = 1;
655 0 0         while (s--) {
656 0 0         if (V == 0)
657 0           is_slpsp = 1;
658 0           Ql = Qk;
659 0 0         V = submod( mont_sqrmod(V,n), addmod(Qk,Qk,n), n);
660 0 0         Qk = mont_sqrmod(Qk,n);
661             }
662 0 0         if (!is_slpsp) return 0; /* slpsp */
663 0 0         if (V != addmod(montQ,montQ,n)) return 0; /* V_{n+1} != 2Q mod n */
664 0           qjacobi = jacobi_iu(Q,n);
665 0 0         Qj = (qjacobi == 0) ? 0 : (qjacobi == 1) ? montQ : n-montQ;
    0          
666 0 0         if (Ql != Qj) return 0; /* n is epsp base Q */
667 0           return 1;
668             } else {
669 22 100         if ( U == 0 && (V == mont2 || V == (n-mont2)) )
    100          
    50          
670 14           return 1;
671 8           s--;
672 20 50         while (s--) {
673 20 100         if (V == 0)
674 8           return 1;
675 12 50         if (s)
676 12 50         V = submod( mont_sqrmod(V,n), mont2, n);
677             }
678             }
679 3           return 0;
680             }
681             #else
682             lucas_seq(&U, &V, &Qk, n, P, Q, d);
683              
684             if (strength == 0) {
685             if (U == 0)
686             return 1;
687             } else if (strength == 1) {
688             if (U == 0)
689             return 1;
690             /* Now check to see if V_{d*2^r} == 0 for any 0 <= r < s */
691             while (s--) {
692             if (V == 0)
693             return 1;
694             if (s) {
695             V = mulsubmod(V, V, addmod(Qk,Qk,n), n);
696             Qk = sqrmod(Qk, n);
697             }
698             }
699             } else if (strength == 2) {
700             UV Ql, Qj = 0;
701             UV Qu = (Q >= 0) ? Q % n : n-(((UV)(-Q)) % n);
702             int qjacobi, is_slpsp = 0;
703             if (U == 0)
704             is_slpsp = 1;
705             while (s--) {
706             if (V == 0)
707             is_slpsp = 1;
708             Ql = Qk;
709             V = mulsubmod(V, V, addmod(Qk,Qk,n), n);
710             Qk = sqrmod(Qk, n);
711             }
712             if (!is_slpsp) return 0; /* slpsp */
713             if (V != addmod(Qu,Qu,n)) return 0; /* V_{n+1} != 2Q mod n */
714             qjacobi = jacobi_iu(Q,n);
715             Qj = (qjacobi == 0) ? 0 : (qjacobi == 1) ? Qu : n-Qu;
716             if (Ql != Qj) return 0; /* n is epsp base Q */
717             return 1;
718             } else {
719             if ( U == 0 && (V == 2 || V == (n-2)) )
720             return 1;
721             /* Now check to see if V_{d*2^r} == 0 for any 0 <= r < s-1 */
722             s--;
723             while (s--) {
724             if (V == 0)
725             return 1;
726             if (s)
727             V = mulsubmod(V, V, 2, n);
728             }
729             }
730             return 0;
731             #endif
732             }
733              
734             /* A generalization of Pari's shortcut to the extra-strong Lucas test.
735             *
736             * This only calculates and tests V, which means less work, but it does result
737             * in a few more pseudoprimes than the full extra-strong test.
738             *
739             * I've added a gcd check at the top, which needs to be done and also results
740             * in fewer pseudoprimes. Pari always does trial division to 100 first so
741             * is unlikely to come up there.
742             *
743             * increment: 1 for Baillie OEIS, 2 for Pari.
744             *
745             * With increment = 1, these results will be a subset of the extra-strong
746             * Lucas pseudoprimes. With increment = 2, we produce Pari's results.
747             */
748 1184           int is_almost_extra_strong_lucas_pseudoprime(UV n, UV increment)
749             {
750             UV P, V, W, d, s, b;
751              
752 1184 50         if (n < 13) return (n == 2 || n == 3 || n == 5 || n == 7 || n == 11);
    0          
    0          
    0          
    0          
    0          
753 1184 50         if ((n % 2) == 0 || n == UV_MAX) return 0;
    50          
754 1184 50         if (increment < 1 || increment > 256)
    50          
755 0           croak("Invalid lucas parameter increment: %"UVuf"\n", increment);
756              
757             /* Ensure small primes work with large increments. */
758 1184 50         if ( (increment >= 16 && n <= 331) || (increment > 148 && n <= 631) )
    0          
    50          
    0          
759 0           return !!is_prob_prime(n);
760              
761 1184           P = select_extra_strong_parameters(n, increment);
762 1184 50         if (P == 0) return 0;
763              
764 1184           d = n+1;
765 1184           s = 0;
766 3535 100         while ( (d & 1) == 0 ) { s++; d >>= 1; }
767 19410 100         { UV v = d; b = 0; while (v >>= 1) b++; }
768              
769             #if USE_MONTMATH
770             {
771 1184           const uint64_t npi = mont_inverse(n), mont1 = mont_get1(n);
772 1184           const uint64_t mont2 = mont_get2(n);
773 1184           const uint64_t montP = mont_geta(P, n);
774 1184 50         W = submod( mont_mulmod( montP, montP, n), mont2, n);
775 1184           V = montP;
776 19410 100         while (b--) {
777 18226 50         UV T = submod( mont_mulmod(V, W, n), montP, n);
778 18226 100         if ( (d >> b) & UVCONST(1) ) {
779 9373           V = T;
780 9373 50         W = submod( mont_mulmod(W, W, n), mont2, n);
781             } else {
782 8853           W = T;
783 8853 50         V = submod( mont_mulmod(V, V, n), mont2, n);
784             }
785             }
786              
787 1184 100         if (V == mont2 || V == (n-mont2))
    100          
788 684           return 1;
789 500           s--;
790 1077 50         while (s--) {
791 1077 100         if (V == 0)
792 500           return 1;
793 577 50         if (s)
794 577 50         V = submod( mont_mulmod(V, V, n), mont2, n);
795             }
796 0           return 0;
797             }
798             #else
799             W = mulsubmod(P, P, 2, n);
800             V = P;
801             while (b--) {
802             UV T = mulsubmod(V, W, P, n);
803             if ( (d >> b) & UVCONST(1) ) {
804             V = T;
805             W = mulsubmod(W, W, 2, n);
806             } else {
807             W = T;
808             V = mulsubmod(V, V, 2, n);
809             }
810             }
811             if (V == 2 || V == (n-2))
812             return 1;
813             while (s-- > 1) {
814             if (V == 0)
815             return 1;
816             V = mulsubmod(V, V, 2, n);
817             if (V == 2)
818             return 0;
819             }
820             return 0;
821             #endif
822             }
823              
824              
825             typedef struct {
826             unsigned short div;
827             unsigned short period;
828             unsigned short offset;
829             } _perrin;
830             #define NPERRINDIV 19
831             /* 1112 mask bytes */
832             static const uint32_t _perrinmask[] = {22,523,514,65890,8519810,130,4259842,0,526338,2147483904U,1644233728,1,8194,1073774592,1024,134221824,128,512,181250,2048,0,1,134217736,1049600,524545,2147500288U,0,524290,536870912,32768,33554432,2048,0,2,2,256,65536,64,536875010,32768,256,64,0,32,1073741824,0,1048576,1048832,371200000,0,0,536887552,32,2147487744U,2097152,32768,1024,0,1024,536870912,128,512,0,0,512,0,2147483650U,45312,128,0,8388640,0,8388608,8388608,0,2048,4096,92800000,262144,0,65536,4,0,4,4,4194304,8388608,1075838976,536870956,0,134217728,8192,0,8192,8192,0,2,0,268435458,134223392,1073741824,268435968,2097152,67108864,0,8192,1073741840,0,0,128,0,0,512,1450000,8,131136,536870928,0,4,2097152,4096,64,0,32768,0,0,131072,371200000,2048,33570816,4096,32,1024,536870912,1048576,16384,0,8388608,0,0,0,2,512,0,128,0,134217728,2,32,0,0,0,0,8192,0,1073742080,536870912,0,4096,16777216,526336,32,0,65536,33554448,708,67108864,2048,0,0,536870912,0,536870912,33554432,33554432,2147483648U,512,64,0,1074003968,512,0,524288,0,0,0,67108864,524288,1048576,0,131076,0,33554432,131072,0,2,8390656,16384,16777216,134217744,0,131104,0,2,32768,0,0,0,1450000,32768,0,0,0,0,0,16,0,1024,16400,1048576,32,1024,0,260,536870912,269484032,0,16384,0,524290,0,0,512,65536,0,0,0,134217732,0,67108880,536887296,0,0,32,0,65568,0,524288,2147483648U,0,4096,4096,134217984,268500992,0,33554432,131072,0,0,0,16777216,0,0,0,0,0,524288,0,0,67108864,0,0,2,0,2,32,1024,0};
833             static _perrin _perrindata[NPERRINDIV] = {
834             {2, 7, 0},
835             {3, 13, 1},
836             {4, 14, 2},
837             {5, 24, 3},
838             {7, 48, 4},
839             {9, 39, 6},
840             {11, 120, 8},
841             {13, 183, 12},
842             {17, 288, 18},
843             {19, 180, 27},
844             {23, 22, 33},
845             {25, 120, 34},
846             {29, 871, 38},
847             {31, 993, 66},
848             {37, 1368, 98},
849             {41, 1723, 141},
850             {43, 231, 195},
851             {47, 2257, 203},
852             {223, 111, 274}
853             };
854              
855             /* Calculate signature using the doubling rule from Adams and Shanks 1982 */
856 20           static void calc_perrin_sig(UV* S, UV n) {
857             #if USE_MONTMATH
858 20           uint64_t npi = 0, mont1;
859             int i;
860             #endif
861             UV T[6], T01, T34, T45;
862             int b;
863              
864             /* Signature for n = 1 */
865 20           S[0] = 1; S[1] = n-1; S[2] = 3; S[3] = 3; S[4] = 0; S[5] = 2;
866 20 50         if (n <= 1) return;
867              
868             #if USE_MONTMATH
869 20 100         if ( (n&1) ) {
870 18           npi = mont_inverse(n);
871 18           mont1 = mont_get1(n);
872 18           S[0] = mont1; S[1] = n-mont1; S[5] = addmod(mont1,mont1,n);
873 18           S[2] = addmod(S[5],mont1,n); S[3] = S[2];
874             }
875             #endif
876              
877             /* Bits in n */
878 625 100         { UV v = n; b = 1; while (v >>= 1) b++; }
879              
880 625 100         while (b-- > 1) {
881             /* Double */
882             #if USE_MONTMATH
883 605 100         if (n&1) {
884 558 50         T[0] = submod(submod(mont_sqrmod(S[0],n), S[5],n), S[5],n);
885 558 50         T[1] = submod(submod(mont_sqrmod(S[1],n), S[4],n), S[4],n);
886 558 50         T[2] = submod(submod(mont_sqrmod(S[2],n), S[3],n), S[3],n);
887 558 50         T[3] = submod(submod(mont_sqrmod(S[3],n), S[2],n), S[2],n);
888 558 50         T[4] = submod(submod(mont_sqrmod(S[4],n), S[1],n), S[1],n);
889 558 50         T[5] = submod(submod(mont_sqrmod(S[5],n), S[0],n), S[0],n);
890             } else
891             #endif
892             {
893 47           T[0] = submod(submod(sqrmod(S[0],n), S[5],n), S[5],n);
894 47           T[1] = submod(submod(sqrmod(S[1],n), S[4],n), S[4],n);
895 47           T[2] = submod(submod(sqrmod(S[2],n), S[3],n), S[3],n);
896 47           T[3] = submod(submod(sqrmod(S[3],n), S[2],n), S[2],n);
897 47           T[4] = submod(submod(sqrmod(S[4],n), S[1],n), S[1],n);
898 47           T[5] = submod(submod(sqrmod(S[5],n), S[0],n), S[0],n);
899             }
900             /* Move to S, filling in */
901 605           T01 = submod(T[2], T[1], n);
902 605           T34 = submod(T[5], T[4], n);
903 605           T45 = addmod(T34, T[3], n);
904 605 100         if ( (n >> (b-1)) & 1U ) {
905 281           S[0] = T[0]; S[1] = T01; S[2] = T[1];
906 281           S[3] = T[4]; S[4] = T45; S[5] = T[5];
907             } else {
908 324           S[0] = T01; S[1] = T[1]; S[2] = addmod(T01,T[0],n);
909 324           S[3] = T34; S[4] = T[4]; S[5] = T45;
910             }
911             }
912             #if USE_MONTMATH
913 20 100         if (n&1) { /* Recover result from Montgomery form */
914 128 100         for (i = 0; i < 6; i++)
915 108 50         S[i] = mont_recover(S[i],n);
916             }
917             #endif
918             }
919              
920 20           int is_perrin_pseudoprime(UV n, int restricted)
921             {
922             int jacobi, i;
923             UV S[6];
924              
925 20 50         if (n < 3) return (n >= 2);
926 20 100         if (!(n&1) && restricted > 2) return 0; /* Odds only for restrict > 2 */
    50          
927              
928             /* Hard code the initial tests. 60% of composites caught by 4 tests. */
929             {
930 20           uint32_t n32 = n % 10920;
931 20 100         if (!(n32&1) && !(( 22 >> (n32% 7)) & 1)) return 0;
    50          
932 20 50         if (!(n32%3) && !(( 523 >> (n32%13)) & 1)) return 0;
    0          
933 20 100         if (!(n32%5) && !((65890 >> (n32%24)) & 1)) return 0;
    50          
934 20 50         if (!(n32%4) && !(( 514 >> (n32%14)) & 1)) return 0;
    0          
935             }
936 320 100         for (i = 4; i < NPERRINDIV; i++) {
937 300 100         if ((n % _perrindata[i].div) == 0) {
938 12           const uint32_t *mask = _perrinmask + _perrindata[i].offset;
939 12           unsigned short mod = n % _perrindata[i].period;
940 12 50         if (!((mask[mod/32] >> (mod%32)) & 1))
941 0           return 0;
942             }
943             }
944             /* Depending on which filters are used, 10-20% of composites are left. */
945              
946 20           calc_perrin_sig(S, n);
947              
948 20 50         if (S[4] != 0) return 0; /* P(n) = 0 mod n */
949 20 100         if (restricted == 0) return 1;
950              
951 5 100         if (S[1] != n-1) return 0; /* P(-n) = -1 mod n */
952 4 100         if (restricted == 1) return 1;
953              
954             /* Full restricted test looks for an acceptable signature.
955             *
956             * restrict = 2 is Adams/Shanks without quadratic form test
957             *
958             * restrict = 3 is Arno or Grantham: No qform, also reject mults of 2 and 23
959             *
960             * See:
961             * Adams/Shanks 1982 pages 257-261
962             * Arno 1991 pages 371-372
963             * Grantham 2000 pages 5-6
964             */
965              
966 3           jacobi = kronecker_su(-23,n);
967              
968 3 100         if (jacobi == -1) { /* Q-type */
969              
970 1           UV B = S[2], B2 = sqrmod(B,n);
971 1           UV A = submod(addmod(1,mulmod(B,3,n),n),B2,n);
972 1           UV C = submod(mulmod(B2,3,n),2,n);
973 1 50         if (S[0] == A && S[2] == B && S[3] == B && S[5] == C &&
    50          
    50          
    50          
    0          
974 0 0         B != 3 && submod(mulmod(B2,B,n),B,n) == 1) {
975 0 0         if (_XS_get_verbose()>1) printf("%"UVuf" Q-Type %"UVuf" -1 %"UVuf" %"UVuf" 0 %"UVuf"\n", n, A, B, B, C);
976 1           return 1;
977             }
978              
979             } else { /* S-Type or I-Type */
980              
981 2 50         if (jacobi == 0 && n != 23 && restricted > 2) {
    50          
    100          
982 1 50         if (_XS_get_verbose()>1) printf("%"UVuf" Jacobi %d\n",n,jacobi);
983 1           return 0; /* Adams/Shanks allows (-23|n) = 0 for S-Type */
984             }
985              
986 1 50         if (S[0] == 1 && S[2] == 3 && S[3] == 3 && S[5] == 2) {
    50          
    50          
    50          
987 1 50         if (_XS_get_verbose()>1) printf("%"UVuf" S-Type 1 -1 3 3 0 2\n",n);
988 1           return 1;
989 0 0         } else if (S[0] == 0 && S[5] == n-1 && S[2] != S[3] &&
    0          
990 0 0         addmod(S[2],S[3],n) == n-3 && sqrmod(submod(S[2],S[3],n),n) == n-(23%n)) {
991 0 0         if (_XS_get_verbose()>1) printf("%"UVuf" I-Type 0 -1 %"UVuf" %"UVuf" 0 -1\n",n, S[2], S[3]);
992 0           return 1;
993             }
994              
995             }
996 1 50         if (_XS_get_verbose()>1) printf("%"UVuf" ? %2d ? %"UVuf" -1 %"UVuf" %"UVuf" 0 %"UVuf"\n", n, jacobi, S[0],S[2],S[3],S[5]);
997 20           return 0;
998             }
999              
1000 28           int is_frobenius_pseudoprime(UV n, IV P, IV Q)
1001             {
1002             UV U, V, Qk, Vcomp;
1003 28           int k = 0;
1004             IV D;
1005             UV Du, Pu, Qu;
1006              
1007 28 50         if (n < 7) return (n == 2 || n == 3 || n == 5);
    0          
    0          
    0          
1008 28 50         if ((n % 2) == 0 || n == UV_MAX) return 0;
    50          
1009              
1010 28 50         if (P == 0 && Q == 0) {
    0          
1011 0           P = -1; Q = 2;
1012 0 0         if (n == 7) P = 1; /* So we don't test kronecker(-7,7) */
1013             do {
1014 0           P += 2;
1015 0 0         if (P == 3) P = 5; /* P=3,Q=2 -> D=9-8=1 => k=1, so skip */
1016 0           D = P*P-4*Q;
1017 0           Du = D >= 0 ? D : -D;
1018 0           k = kronecker_su(D, n);
1019 0 0         if (P == 10001 && is_perfect_square(n)) return 0;
    0          
1020 0 0         } while (k == 1);
1021 0 0         if (k == 0) return 0;
1022             /* D=P^2-8 will not be a perfect square */
1023 0 0         if (_XS_get_verbose()) printf("%"UVuf" Frobenius (%"IVdf",%"IVdf") : x^2 - %"IVdf"x + %"IVdf"\n", n, P, Q, P, Q);
1024 0           Vcomp = 4;
1025             } else {
1026 28           D = P*P-4*Q;
1027 28           Du = D >= 0 ? D : -D;
1028 28 100         if (D != 5 && is_perfect_square(Du))
    50          
1029 0           croak("Frobenius invalid P,Q: (%"IVdf",%"IVdf")", P, Q);
1030             }
1031 28           Pu = (P >= 0 ? P : -P) % n;
1032 28           Qu = (Q >= 0 ? Q : -Q) % n;
1033              
1034 28           Qk = gcd_ui(n, Pu*Qu*Du);
1035 28 50         if (Qk != 1) {
1036 0 0         if (Qk == n) return !!is_prob_prime(n);
1037 0           return 0;
1038             }
1039 28 50         if (k == 0) {
1040 28           k = kronecker_su(D, n);
1041 28 50         if (k == 0) return 0;
1042 28 100         if (k == 1) {
1043 24           Vcomp = 2;
1044             } else {
1045 4           Qu = addmod(Qu,Qu,n);
1046 4 50         Vcomp = (Q >= 0) ? Qu : n-Qu;
1047             }
1048             }
1049              
1050 28           lucas_seq(&U, &V, &Qk, n, P, Q, n-k);
1051             /* if (_XS_get_verbose()) printf("%"UVuf" Frobenius U = %"UVuf" V = %"UVuf"\n", n, U, V); */
1052 28 50         if (U == 0 && V == Vcomp) return 1;
    50          
1053 28           return 0;
1054             }
1055              
1056             /*
1057             * Khashin, 2013, "Counterexamples for Frobenius primality test"
1058             * http://arxiv.org/abs/1307.7920
1059             * 1. Select c as first odd prime where (c,n)=-1.
1060             * 2. Check (1 + sqrt(c))^n mod n equiv (1 - sqrt(c) mod n
1061             *
1062             * His Sep 2016 talk starts with c = -1,2 using checks:
1063             * (2+sqrt(c)^n = 2-sqrt(c) mod n for c = -1,2
1064             * (1+sqrt(c)^n = 1-sqrt(c) mod n for c = odd prime
1065             * There doesn't seem to be a big advantage for this change.
1066             */
1067 102           int is_frobenius_khashin_pseudoprime(UV n)
1068             {
1069             int k;
1070 102           UV ra, rb, a, b, d = n-1, c = 1;
1071              
1072 102 50         if (n < 7) return (n == 2 || n == 3 || n == 5);
    0          
    0          
    0          
1073 102 50         if ((n % 2) == 0 || n == UV_MAX) return 0;
    50          
1074 102 50         if (is_perfect_square(n)) return 0;
1075              
1076             /* c = first odd prime where (c|n)=-1 */
1077             do {
1078 188           c += 2;
1079 188 100         if (c==9 || (c>=15 && (!(c%3) || !(c%5) || !(c%7) || !(c%11) || !(c%13))))
    100          
    100          
    100          
    50          
    50          
    50          
1080 12           continue;
1081 176           k = kronecker_uu(c, n);
1082 188 100         } while (k == 1);
1083 102 100         if (k == 0) return 0;
1084              
1085             #if USE_MONTMATH
1086             {
1087 61           const uint64_t npi = mont_inverse(n);
1088 61           const uint64_t mont1 = mont_get1(n);
1089 61           const uint64_t montc = mont_geta(c, n);
1090 61           ra = rb = a = b = mont1;
1091 1915 100         while (d) {
1092 1854 100         if (d & 1) {
1093 907           UV ta=ra, tb=rb;
1094 907 50         ra = addmod( mont_mulmod(ta,a,n), mont_mulmod(mont_mulmod(tb,b,n),montc,n), n );
    0          
    50          
    50          
1095 907 50         rb = addmod( mont_mulmod(tb,a,n), mont_mulmod(ta,b,n), n);
    50          
1096             }
1097 1854           d >>= 1;
1098 1854 100         if (d) {
1099 1793 50         UV t = mont_mulmod(mont_mulmod(b,b,n),montc,n);
    0          
    50          
1100 1793 50         b = mont_mulmod(b,a,n);
1101 1793           b = addmod(b,b,n);
1102 1793 50         a = addmod(mont_mulmod(a,a,n),t,n);
1103             }
1104             }
1105 61 100         return (ra == mont1 && rb == n-mont1);
    50          
1106             }
1107             #else
1108             ra = rb = a = b = 1;
1109             while (d) {
1110             if (d & 1) {
1111             /* This is faster than the 3-mulmod 5-addmod version */
1112             UV ta=ra, tb=rb;
1113             ra = addmod( mulmod(ta,a,n), mulmod(mulmod(tb,b,n),c,n), n );
1114             rb = addmod( mulmod(tb,a,n), mulmod(ta,b,n), n);
1115             }
1116             d >>= 1;
1117             if (d) {
1118             UV t = mulmod(sqrmod(b,n),c,n);
1119             b = mulmod(b,a,n);
1120             b = addmod(b,b,n);
1121             a = addmod(sqrmod(a,n),t,n);
1122             }
1123             }
1124             return (ra == 1 && rb == n-1);
1125             #endif
1126             }
1127              
1128             /*
1129             * The Frobenius-Underwood test has no known counterexamples below 2^50, but
1130             * has not been extensively tested above that. This is the Minimal Lambda+2
1131             * test from section 9 of "Quadratic Composite Tests" by Paul Underwood.
1132             *
1133             * It is generally slower than the AES Lucas test, but for large values is
1134             * competitive with the BPSW test. Since our BPSW is known to have no
1135             * counterexamples under 2^64, while the results of this test are unknown,
1136             * it is mainly useful for numbers larger than 2^64 as an additional
1137             * non-correlated test.
1138             */
1139 102           int is_frobenius_underwood_pseudoprime(UV n)
1140             {
1141             int j, bit;
1142             UV x, result, a, b, np1, len, t1;
1143             IV t;
1144              
1145 102 50         if (n < 7) return (n == 2 || n == 3 || n == 5);
    0          
    0          
    0          
1146 102 50         if ((n % 2) == 0 || n == UV_MAX) return 0;
    50          
1147              
1148 171 50         for (x = 0; x < 1000000; x++) {
1149 171 100         if (x==2 || x==4 || x==7 || x==8 || x==10 || x==14 || x==16 || x==18)
    100          
    50          
    50          
    50          
    50          
    50          
    50          
1150 14           continue;
1151 157           t = (IV)(x*x) - 4;
1152 157           j = jacobi_iu(t, n);
1153 157 100         if (j == -1) break;
1154 67 100         if (j == 0 || (x == 20 && is_perfect_square(n)))
    50          
    0          
1155 12           return 0;
1156             }
1157 90 50         if (x >= 1000000) croak("FU test failure, unable to find suitable a");
1158 90           t1 = gcd_ui(n, (x+4)*(2*x+5));
1159 90 100         if (t1 != 1 && t1 != n)
    50          
1160 29           return 0;
1161 61           np1 = n+1;
1162 1858 100         { UV v = np1; len = 1; while (v >>= 1) len++; }
1163              
1164             #if USE_MONTMATH
1165             {
1166 61           const uint64_t npi = mont_inverse(n), mont1 = mont_get1(n);
1167 61           const uint64_t mont2 = mont_get2(n);
1168 61           const uint64_t mont5 = mont_geta(5, n);
1169              
1170 61           x = mont_geta(x, n);
1171 61           a = mont1;
1172 61           b = mont2;
1173              
1174 61 100         if (x == 0) {
1175 44           result = mont5;
1176 1336 100         for (bit = len-2; bit >= 0; bit--) {
1177 1292           t1 = addmod(b, b, n);
1178 1292 50         b = mont_mulmod(submod(b, a, n), addmod(b, a, n), n);
1179 1292 50         a = mont_mulmod(a, t1, n);
1180 1292 100         if ( (np1 >> bit) & UVCONST(1) ) {
1181 568           t1 = b;
1182 568           b = submod( addmod(b, b, n), a, n);
1183 568           a = addmod( addmod(a, a, n), t1, n);
1184             }
1185             }
1186             } else {
1187 17           UV multiplier = addmod(x, mont2, n);
1188 17           result = addmod( addmod(x, x, n), mont5, n);
1189 522 100         for (bit = len-2; bit >= 0; bit--) {
1190 505 50         t1 = addmod( mont_mulmod(a, x, n), addmod(b, b, n), n);
1191 505 50         b = mont_mulmod(submod(b, a, n), addmod(b, a, n), n);
1192 505 50         a = mont_mulmod(a, t1, n);
1193 505 100         if ( (np1 >> bit) & UVCONST(1) ) {
1194 262           t1 = b;
1195 262           b = submod( addmod(b, b, n), a, n);
1196 262 50         a = addmod( mont_mulmod(a, multiplier, n), t1, n);
1197             }
1198             }
1199             }
1200 61 100         return (a == 0 && b == result);
    50          
1201             }
1202             #else
1203             a = 1;
1204             b = 2;
1205              
1206             if (x == 0) {
1207             result = 5;
1208             for (bit = len-2; bit >= 0; bit--) {
1209             t1 = addmod(b, b, n);
1210             b = mulmod( submod(b, a, n), addmod(b, a, n), n);
1211             a = mulmod(a, t1, n);
1212             if ( (np1 >> bit) & UVCONST(1) ) {
1213             t1 = b;
1214             b = submod( addmod(b, b, n), a, n);
1215             a = addmod( addmod(a, a, n), t1, n);
1216             }
1217             }
1218             } else {
1219             UV multiplier = addmod(x, 2, n);
1220             result = addmod( addmod(x, x, n), 5, n);
1221             for (bit = len-2; bit >= 0; bit--) {
1222             t1 = addmod( mulmod(a, x, n), addmod(b, b, n), n);
1223             b = mulmod(submod(b, a, n), addmod(b, a, n), n);
1224             a = mulmod(a, t1, n);
1225             if ( (np1 >> bit) & UVCONST(1) ) {
1226             t1 = b;
1227             b = submod( addmod(b, b, n), a, n);
1228             a = addmod( mulmod(a, multiplier, n), t1, n);
1229             }
1230             }
1231             }
1232              
1233             if (_XS_get_verbose()>1) printf("%"UVuf" is %s with x = %"UVuf"\n", n, (a == 0 && b == result) ? "probably prime" : "composite", x);
1234             if (a == 0 && b == result)
1235             return 1;
1236             return 0;
1237             #endif
1238             }
1239              
1240             /* We have a native-UV Lucas-Lehmer test with simple pretest. If 2^p-1 is
1241             * prime but larger than a UV, we'll have to bail, and they'll run the nice
1242             * GMP version. However, they're just asking if this is a Mersenne prime, and
1243             * there are millions of CPU years that have gone into enumerating them, so
1244             * instead we'll use a table. */
1245             #define NUM_KNOWN_MERSENNE_PRIMES 49
1246             static const uint32_t _mersenne_primes[NUM_KNOWN_MERSENNE_PRIMES] = {2,3,5,7,13,17,19,31,61,89,107,127,521,607,1279,2203,2281,3217,4253,4423,9689,9941,11213,19937,21701,23209,44497,86243,110503,132049,216091,756839,859433,1257787,1398269,2976221,3021377,6972593,13466917,20996011,24036583,25964951,30402457,32582657,37156667,42643801,43112609,57885161,74207281};
1247             #define LAST_CHECKED_MERSENNE 40364833
1248 2282           int is_mersenne_prime(UV p)
1249             {
1250             int i;
1251 113403 100         for (i = 0; i < NUM_KNOWN_MERSENNE_PRIMES; i++)
1252 111138 100         if (p == _mersenne_primes[i])
1253 17           return 1;
1254 2265 50         return (p < LAST_CHECKED_MERSENNE) ? 0 : -1;
1255             }
1256 0           int lucas_lehmer(UV p)
1257             {
1258             UV k, V, mp;
1259              
1260 0 0         if (p == 2) return 1;
1261 0 0         if (!is_prob_prime(p)) return 0;
1262 0 0         if (p > BITS_PER_WORD) croak("lucas_lehmer with p > BITS_PER_WORD");
1263 0           V = 4;
1264 0           mp = UV_MAX >> (BITS_PER_WORD - p);
1265 0 0         for (k = 3; k <= p; k++) {
1266 0           V = mulsubmod(V, V, 2, mp);
1267             }
1268 0           return (V == 0);
1269             }
1270              
1271             /******************************************************************************/
1272              
1273             /* Hashing similar to Forišek and Jančina 2015 */
1274             static const uint16_t mr_bases_hash32[256] = {
1275             157,1150,304,8758,362,15524,1743,212,1056,1607,140,3063,160,913,5842,2013,598,1929,696,1474,3006,524,155,705,694,1238,1851,1053,585,626,603,222,1109,1105,604,646,606,1249,1553,5609,515,548,1371,152,2824,532,3556,831,88,185,1355,501,1556,317,582,4739,4710,145,1045,2976,2674,318,1293,10934,1434,1178,3159,26,3526,1859,6467,602,699,5113,3152,2002,2361,101,464,68,813,446,1368,4637,368,1068,307,2820,6189,10457,569,1690,551,237,226,3235,405,3179,1101,610,56,14647,1687,247,8109,5172,1725,1248,536,2869,1047,899,12285,1026,250,1867,1432,336,5175,1632,5169,39,362,290,1372,11988,1329,2168,34,8781,495,399,34,29,4333,1669,166,6405,7357,694,579,746,1278,6347,7751,179,1085,11734,1615,3575,4253,7894,3097,591,1354,1676,151,702,7,5607,2565,440,566,112,3622,1241,1193,2324,1530,1423,548,3341,2012,6305,2410,39,106,3046,1507,1325,1807,2323,5645,1524,1301,1522,238,1226,2476,2126,1677,3288,1981,18481,287,1011,2877,563,7654,1231,776,3907,117,174,1124,199,16838,164,41,313,1692,1574,1021,2804,1093,1263,956,8508,1221,3743,1318,1304,1344,7628,10739,228,30,520,103,1621,6278,847,4537,272,2213,1989,1826,915,318,401,924,227,911,15505,1670,212,1391,700,3254,4931,3637,2822,1726,137,1843,1300
1276             };
1277              
1278              
1279 177009           int is_prob_prime(UV n)
1280             {
1281             int ret;
1282              
1283 177009 100         if (n < 11) {
1284 6295 50         if (n == 2 || n == 3 || n == 5 || n == 7) return 2;
    100          
    100          
    100          
1285 207           else return 0;
1286             }
1287              
1288             #if BITS_PER_WORD == 64
1289 170714 100         if (n <= UVCONST(4294967295)) {
1290             #else
1291             {
1292             #endif
1293 167339           uint32_t x = n;
1294             UV base;
1295 205515 100         if (!(x%2) || !(x%3) || !(x%5) || !(x%7)) return 0;
    100          
    100          
    100          
1296 98933 100         if (x < 121) /* 11*11 */ return 2;
1297 96338 100         if (!(x%11) || !(x%13) || !(x%17) || !(x%19) ||
    100          
    100          
    100          
    100          
1298 73892 100         !(x%23) || !(x%29) || !(x%31) || !(x%37) ||
    100          
    100          
    100          
1299 96338 100         !(x%41) || !(x%43) || !(x%47) || !(x%53)) return 0;
    100          
    100          
1300 61741 100         if (x < 3481) /* 59*59 */ return 2;
1301             /* Trial division crossover point depends on platform */
1302             if (!USE_MONTMATH && n < 500000) {
1303             uint32_t f = 59;
1304             uint32_t limit = isqrt(n);
1305             while (f <= limit) {
1306             if ((x%f) == 0) return 0; f += 2;
1307             if ((x%f) == 0) return 0; f += 6;
1308             if ((x%f) == 0) return 0; f += 4;
1309             if ((x%f) == 0) return 0; f += 2;
1310             if ((x%f) == 0) return 0; f += 4;
1311             if ((x%f) == 0) return 0; f += 2;
1312             if ((x%f) == 0) return 0; f += 4;
1313             if ((x%f) == 0) return 0; f += 6;
1314             }
1315             return 2;
1316             }
1317             /* Use Mueller-like 32-bit hash to find single M-R base to use. */
1318 60757           x = (((x >> 16) ^ x) * 0x45d9f3b) & 0xFFFFFFFFUL;
1319 60757           x = ((x >> 16) ^ x) & 255;
1320 60757           base = mr_bases_hash32[x];
1321 60757           ret = miller_rabin(n, &base, 1);
1322             #if BITS_PER_WORD == 64
1323             } else { /* input is >= 2^32, UV is 64-bit*/
1324 3375 50         if (!(n%2) || !(n%3) || !(n%5) || !(n%7)) return 0;
    100          
    100          
    100          
1325 1841 100         if (!(n%11) || !(n%13) || !(n%17) || !(n%19) ||
    100          
    100          
    100          
    100          
1326 1438 100         !(n%23) || !(n%29) || !(n%31) || !(n%37) ||
    100          
    100          
    100          
1327 1841 100         !(n%41) || !(n%43) || !(n%47) || !(n%53)) return 0;
    100          
    100          
1328 1234 100         if (!(n%59) || !(n%61) || !(n%67) || !(n%71)) return 0;
    100          
    100          
    100          
1329 1180 100         if (!(n%73) || !(n%79) || !(n%83) || !(n%89)) return 0;
    100          
    100          
    100          
1330             /*
1331             * AESLSP test costs about 1.5 Selfridges, vs. ~2.2 for strong Lucas.
1332             * This makes the full BPSW test cost about 2.5x M-R tests for a prime.
1333             */
1334 1144           ret = BPSW(n);
1335             #endif
1336             }
1337 61901           return 2*ret;
1338             }