1 #include "cado.h" // IWYU pragma: keep
2 #include <stdio.h>
3 /* Routines that are the same for modredc_15ul.c and modredc_2ul2.c */
4
5
6 /* Divide residue by 3. Returns 1 if division is possible, 0 otherwise.
7 Assumes that a+3m does not overflow */
8
9 #include "mod_common.c"
10 #include "macros.h"
11
12 int
mod_div3(residue_t r,const residue_t a,const modulus_t m)13 mod_div3 (residue_t r, const residue_t a, const modulus_t m)
14 {
15 residue_t t;
16 unsigned long a3 = (a[1] % 256UL + a[1] / 256UL +
17 a[0] % 256UL + a[0] / 256UL) % 3UL;
18 const unsigned long m3 = (m[0].m[0] % 256UL + m[0].m[0] / 256UL +
19 m[0].m[1] % 256UL + m[0].m[1] / 256UL) % 3UL;
20 #ifdef WANT_ASSERT_EXPENSIVE
21 residue_t a_backup;
22 #endif
23
24 if (m3 == 0)
25 return 0;
26
27 #ifdef WANT_ASSERT_EXPENSIVE
28 mod_init_noset0 (a_backup, m);
29 mod_set (a_backup, a, m);
30 #endif
31
32 mod_init_noset0 (t, m);
33 t[1] = a[1];
34 t[0] = a[0];
35
36 /* Make t[1]:t[0] divisible by 3 */
37 if (a3 != 0UL)
38 {
39 if (a3 + m3 == 3UL) /* Hence a3 == 1, m3 == 2 or a3 == 2, m3 == 1 */
40 {
41 ularith_add_2ul_2ul (&(t[0]), &(t[1]), m[0].m[0], m[0].m[1]);
42 }
43 else /* a3 == 1, m3 == 1 or a3 == 2, m3 == 2 */
44 {
45 ularith_add_2ul_2ul (&(t[0]), &(t[1]), m[0].m[0], m[0].m[1]);
46 ularith_add_2ul_2ul (&(t[0]), &(t[1]), m[0].m[0], m[0].m[1]);
47 }
48
49 /* Now t[1]:t[0] is divisible by 3 */
50 ASSERT_EXPENSIVE ((t[0] % 3UL + t[1] % 3UL) % 3UL == 0UL);
51 }
52
53 /* a = a1 * 2^w + a0, 3|a
54 Let a = a' * 3 * 2^w + a'', a'' < 3 * 2^w.
55 3 | a'', a'' / 3 < 2^w
56 So a / 3 = a' * w + a'' / 3
57 a' = trunc(a1 / 3)
58 a'' = a0 * 3^{-1} (mod 2^w)
59 Hence we get the correct result with one one-word multiplication
60 and one one-word truncating division by a small constant.
61 */
62
63 r[1] = t[1] / 3UL;
64 #if LONG_BIT == 32
65 r[0] = t[0] * 0xaaaaaaabUL; /* 1/3 (mod 2^32) */
66 #elif LONG_BIT == 64
67 r[0] = t[0] * 0xaaaaaaaaaaaaaaabUL; /* 1/3 (mod 2^64) */
68 #else
69 #error Unknown word size
70 #endif
71
72 #ifdef WANT_ASSERT_EXPENSIVE
73 mod_add (t, r, r, m);
74 mod_add (t, t, r, m);
75 ASSERT_EXPENSIVE (mod_equal (a_backup, t, m));
76 mod_clear (a_backup, m);
77 #endif
78 mod_clear (t, m);
79
80 return 1;
81 }
82
83
84 /* Divide residue by 5. Returns 1 if division is possible, 0 otherwise */
85
86 int
mod_div5(residue_t r,const residue_t a,const modulus_t m)87 mod_div5 (residue_t r, const residue_t a, const modulus_t m)
88 {
89 /* inv_5[i] = -1/i (mod 5) */
90 const unsigned long inv_5[5] = {0,4,2,3,1};
91 #if LONG_BIT == 32
92 const unsigned long c = 0xcccccccdUL; /* 1/5 (mod 2^32) */
93 #elif LONG_BIT == 64
94 const unsigned long c = 0xcccccccccccccccdUL; /* 1/5 (mod 2^64) */
95 #else
96 #error Unknown word size
97 #endif
98
99 return mod_divn (r, a, 5UL, 1UL, inv_5, c, m);
100 }
101
102
103 /* Divide residue by 7. Returns 1 if division is possible, 0 otherwise */
104
105 int
mod_div7(residue_t r,const residue_t a,const modulus_t m)106 mod_div7 (residue_t r, const residue_t a, const modulus_t m)
107 {
108 const unsigned long w_mod_7 = (sizeof (unsigned long) == 4) ? 4UL : 2UL;
109 /* inv_7[i] = -1/i (mod 7) */
110 const unsigned long inv_7[7] = {0,6,3,2,5,4,1};
111 #if LONG_BIT == 32
112 const unsigned long c = 0xb6db6db7UL; /* 1/7 (mod 2^32) */
113 #elif LONG_BIT == 64
114 const unsigned long c = 0x6db6db6db6db6db7UL; /* 1/7 (mod 2^64) */
115 #else
116 #error Unknown word size
117 #endif
118
119 return mod_divn (r, a, 7UL, w_mod_7, inv_7, c, m);
120 }
121
122
123 /* Divide residue by 11. Returns 1 if division is possible, 0 otherwise */
124
125 int
mod_div11(residue_t r,const residue_t a,const modulus_t m)126 mod_div11 (residue_t r, const residue_t a, const modulus_t m)
127 {
128 const unsigned long w_mod_11 = (sizeof (unsigned long) == 4) ? 4UL : 5UL;
129 /* inv_11[i] = -1/i (mod 11) */
130 const unsigned long inv_11[11] = {0, 10, 5, 7, 8, 2, 9, 3, 4, 6, 1};
131 #if LONG_BIT == 32
132 const unsigned long c = 0xba2e8ba3UL; /* 1/11 (mod 2^32) */
133 #elif LONG_BIT == 64
134 const unsigned long c = 0x2e8ba2e8ba2e8ba3UL; /* 1/11 (mod 2^64) */
135 #else
136 #error Unknown word size
137 #endif
138
139 return mod_divn (r, a, 11UL, w_mod_11, inv_11, c, m);
140 }
141
142
143 /* Divide residue by 13. Returns 1 if division is possible, 0 otherwise */
144
145 int
mod_div13(residue_t r,const residue_t a,const modulus_t m)146 mod_div13 (residue_t r, const residue_t a, const modulus_t m)
147 {
148 const unsigned long w_mod_13 = (sizeof (unsigned long) == 4) ? 9UL : 3UL;
149 /* inv_13[i] = -1/i (mod 13) */
150 const unsigned long inv_13[13] = {0, 12, 6, 4, 3, 5, 2, 11, 8, 10, 9, 7, 1};
151 #if LONG_BIT == 32
152 const unsigned long c = 0xc4ec4ec5UL; /* 1/13 (mod 2^32) */
153 #elif LONG_BIT == 64
154 const unsigned long c = 0x4ec4ec4ec4ec4ec5UL; /* 1/13 (mod 2^64) */
155 #else
156 #error Unknown word size
157 #endif
158
159 return mod_divn (r, a, 13UL, w_mod_13, inv_13, c, m);
160 }
161
162
163 void
mod_gcd(modint_t r,const residue_t A,const modulus_t m)164 mod_gcd (modint_t r, const residue_t A, const modulus_t m)
165 {
166 unsigned long a[2], b[2];
167 int sh;
168
169 /* Since we do REDC arithmetic, we must have m odd */
170 ASSERT_EXPENSIVE (m[0].m[0] & 1UL);
171
172 if (A[0] == 0UL && A[1] == 0UL)
173 {
174 r[0] = m[0].m[0];
175 r[1] = m[0].m[1];
176 return;
177 }
178
179 a[0] = A[0];
180 a[1] = A[1];
181 b[0] = m[0].m[0];
182 b[1] = m[0].m[1];
183
184 while (a[1] != 0UL || a[0] != 0UL)
185 {
186 /* Make a odd */
187 #if LOOKUP_TRAILING_ZEROS
188 do {
189 sh = trailing_zeros [(unsigned char) a[0]];
190 ularith_shrd (&(a[0]), a[1], a[0], sh);
191 *(long *) &(a[1]) >>= sh;
192 } while (sh == 8);
193 #else
194 if (a[0] == 0UL) /* ctzl does not like zero input */
195 {
196 a[0] = a[1];
197 a[1] = ((long)a[1] < 0L) ? (unsigned long) (-1L) : 0UL;
198 }
199 sh = ularith_ctz (a[0]);
200 ularith_shrd (&(a[0]), a[1], a[0], sh);
201 *(long *) &(a[1]) >>= sh;
202 #endif
203
204 /* Try to make the low two bits of b[0] zero */
205 ASSERT_EXPENSIVE (a[0] % 2UL == 1UL);
206 ASSERT_EXPENSIVE (b[0] % 2UL == 1UL);
207 if ((a[0] ^ b[0]) & 2UL)
208 ularith_add_2ul_2ul (&(b[0]), &(b[1]), a[0], a[1]);
209 else
210 ularith_sub_2ul_2ul (&(b[0]), &(b[1]), a[0], a[1]);
211
212 if (b[0] == 0UL && b[1] == 0UL)
213 {
214 if ((long) a[1] < 0L)
215 {
216 a[1] = -a[1];
217 if (a[0] != 0UL)
218 a[1]--;
219 a[0] = -a[0];
220 }
221 r[0] = a[0];
222 r[1] = a[1];
223 return;
224 }
225
226 /* Make b odd */
227 #if LOOKUP_TRAILING_ZEROS
228 do {
229 sh = trailing_zeros [(unsigned char) b[0]];
230 ularith_shrd (&(b[0]), b[1], b[0], sh);
231 *(long *) &(b[1]) >>= sh;
232 } while (sh == 8);
233 #else
234 if (b[0] == 0UL) /* ctzl does not like zero input */
235 {
236 b[0] = b[1];
237 b[1] = ((long)b[1] < 0) ? (unsigned long) (-1L) : 0UL;
238 }
239 sh = ularith_ctz (b[0]);
240 ularith_shrd (&(b[0]), b[1], b[0], sh);
241 *(long *) &(b[1]) >>= sh;
242 #endif
243 ASSERT_EXPENSIVE (a[0] % 2UL == 1UL);
244 ASSERT_EXPENSIVE (b[0] % 2UL == 1UL);
245
246 if ((a[0] ^ b[0]) & 2)
247 ularith_add_2ul_2ul (&(a[0]), &(a[1]), b[0], b[1]);
248 else
249 ularith_sub_2ul_2ul (&(a[0]), &(a[1]), b[0], b[1]);
250 }
251
252 if ((long) b[1] < 0)
253 {
254 b[1] = -b[1];
255 if (b[0] != 0UL)
256 b[1]--;
257 b[0] = -b[0];
258 }
259 r[0] = b[0];
260 r[1] = b[1];
261
262 return;
263 }
264
265
266 /* Simple addition chains for small multipliers, used in the powering
267 functions with small bases */
268 static inline void
simple_mul(residue_t r,const residue_t a,const unsigned long b,const modulus_t m)269 simple_mul (residue_t r, const residue_t a, const unsigned long b,
270 const modulus_t m)
271 {
272 if (b == 2UL) {
273 mod_add (r, a, a, m);
274 } else if (b == 3UL) {
275 mod_add (r, a, a, m);
276 mod_add (r, r, a, m);
277 } else if (b == 5UL) {
278 mod_add (r, a, a, m);
279 mod_add (r, r, r, m);
280 mod_add (r, r, a, m);
281 } else {
282 ASSERT (b == 7UL);
283 mod_add (r, a, a, m);
284 mod_add (r, r, r, m);
285 mod_add (r, r, r, m);
286 mod_sub (r, r, a, m);
287 }
288 }
289
290 /* Compute r = b^e, where b is a small integer, currently b=2,3,5,7 are
291 implemented. Here, e is an unsigned long */
292 static inline void
mod_npow_ul(residue_t r,const unsigned long b,const unsigned long e,const modulus_t m)293 mod_npow_ul (residue_t r, const unsigned long b, const unsigned long e,
294 const modulus_t m)
295 {
296 unsigned long mask;
297 residue_t t, u;
298
299 if (e == 0UL)
300 {
301 mod_set1 (r, m);
302 return;
303 }
304
305 mod_init_noset0 (t, m);
306
307 if (b == 2UL) {
308 mod_set1 (t, m);
309 mod_add (t, t, t, m); /* t = 2 */
310 } else {
311 mod_init_noset0 (u, m);
312 mod_set1 (u, m);
313 simple_mul (t, u, b, m); /* t = b */
314 }
315
316 mask = (1UL << (LONG_BIT - 1)) >> ularith_clz (e);
317 ASSERT (e & mask);
318 mask >>= 1;
319
320 while (mask > 0UL)
321 {
322 mod_sqr (t, t, m);
323 if (b == 2UL) {
324 mod_intshl (t, t, (e & mask) ? 1 : 0);
325 ularith_sub_2ul_2ul_ge (&(t[0]), &(t[1]), m[0].m[0], m[0].m[1]);
326 } else {
327 simple_mul (u, t, b, m);
328 if (e & mask)
329 mod_set (t, u, m);
330 }
331 mask >>= 1;
332 }
333 mod_set (r, t, m);
334 mod_clear (t, m);
335 if (b != 2UL)
336 mod_clear (u, m);
337 }
338
339
340 /* Compute r = 2^e mod m. Here, e is an unsigned long */
341 void
mod_2pow_ul(residue_t r,const unsigned long e,const modulus_t m)342 mod_2pow_ul (residue_t r, const unsigned long e, const modulus_t m)
343 {
344 mod_npow_ul (r, 2UL, e, m);
345 }
346
347
348 /* Compute r = b^e, where b is a small integer, currently b=2,3,5,7 are
349 implemented. Here e is a multiple precision integer
350 sum_{i=0}^{e_nrwords-1} e[i] * (machine word base)^i */
351 static inline void
mod_npow_mp(residue_t r,const unsigned long b,const unsigned long * e,const int e_nrwords,const modulus_t m)352 mod_npow_mp (residue_t r, const unsigned long b, const unsigned long *e,
353 const int e_nrwords, const modulus_t m)
354 {
355 residue_t t, u;
356 int i = e_nrwords - 1;
357 unsigned long mask, ei;
358
359 if (e_nrwords == 0 || e[i] == 0UL)
360 {
361 mod_set1 (r, m);
362 return;
363 }
364
365 mod_init_noset0 (t, m);
366 if (b == 2UL) {
367 mod_set1 (t, m);
368 mod_add (t, t, t, m); /* t = 2 */
369 } else {
370 mod_init_noset0 (u, m);
371 mod_set1 (u, m);
372 simple_mul (t, u, b, m); /* t = b */
373 }
374
375 mask = (1UL << (LONG_BIT - 1)) >> ularith_clz (e[i]);
376 mask >>= 1;
377
378 for ( ; i >= 0; i--)
379 {
380 ei = e[i];
381 while (mask > 0UL)
382 {
383 mod_sqr (t, t, m);
384 if (b == 2UL) {
385 mod_intshl (t, t, (ei & mask) ? 1 : 0);
386 ularith_sub_2ul_2ul_ge (&(t[0]), &(t[1]), m[0].m[0], m[0].m[1]);
387 } else {
388 simple_mul (u, t, b, m);
389 if (ei & mask)
390 mod_set (t, u, m);
391 }
392 mask >>= 1; /* (r^2)^(mask/2) * b^e = r^mask * b^e */
393 }
394 mask = ~0UL - (~0UL >> 1);
395 }
396 mod_set (r, t, m);
397 mod_clear (t, m);
398 if (b != 2UL)
399 mod_clear (u, m);
400 }
401
402
403 /* Compute r = 2^e mod m. Here e is a multiple precision integer
404 sum_{i=0}^{e_nrwords-1} e[i] * (machine word base)^i */
405 void
mod_2pow_mp(residue_t r,const unsigned long * e,const int e_nrwords,const modulus_t m)406 mod_2pow_mp (residue_t r, const unsigned long *e, const int e_nrwords,
407 const modulus_t m)
408 {
409 mod_npow_mp (r, 2UL, e, e_nrwords, m);
410 }
411
412
413 /* Returns 1 if m is a strong probable prime wrt base b, 0 otherwise.
414 We assume m is odd. */
415 int
mod_sprp(const residue_t b,const modulus_t m)416 mod_sprp (const residue_t b, const modulus_t m)
417 {
418 residue_t r, minusone;
419 int i = 0, po2 = 0;
420 modint_t mm1;
421
422 mod_getmod_int (mm1, m);
423
424 if (mod_intequal_ul (mm1, 1UL))
425 return 0;
426
427 /* Let m-1 = 2^l * k, k odd. Set mm1 = k, po2 = l */
428 mm1[0]--; /* No borrow since m is odd */
429 if (mm1[0] == 0UL)
430 {
431 mm1[0] = mm1[1];
432 mm1[1] = 0UL;
433 po2 += LONG_BIT;
434 }
435 ASSERT (mm1[0] != 0UL);
436 i = ularith_ctz (mm1[0]);
437 mod_intshr (mm1, mm1, i);
438 po2 += i;
439
440 mod_init_noset0 (r, m);
441 mod_init_noset0 (minusone, m);
442 mod_set1 (minusone, m);
443 mod_neg (minusone, minusone, m);
444
445 /* Exponentiate */
446 if (mm1[1] > 0UL)
447 mod_pow_mp (r, b, mm1, 2UL, m);
448 else
449 mod_pow_ul (r, b, mm1[0], m);
450
451 i = find_minus1 (r, minusone, po2, m);
452
453 mod_clear (r, m);
454 mod_clear (minusone, m);
455 return i;
456 }
457
458
459 /* Returns 1 if m is a strong probable prime wrt base 2, 0 otherwise.
460 We assume m is odd. */
461 int
mod_sprp2(const modulus_t m)462 mod_sprp2 (const modulus_t m)
463 {
464 residue_t r, minusone;
465 int i = 0, po2 = 0;
466 modint_t mm1;
467
468 mod_getmod_int (mm1, m);
469
470 /* Let m-1 = 2^l * k, k odd. Set mm1 = k, po2 = l */
471 mm1[0]--; /* No borrow since m is odd */
472
473 if (mm1[0] == 0UL)
474 {
475 mm1[0] = mm1[1];
476 mm1[1] = 0UL;
477 po2 += LONG_BIT;
478 }
479 ASSERT (mm1[0] != 0UL);
480 i = ularith_ctz (mm1[0]);
481 mod_intshr (mm1, mm1, i);
482 po2 += i;
483
484 mod_init_noset0 (r, m);
485 mod_init_noset0 (minusone, m);
486 mod_set1 (minusone, m);
487 mod_neg (minusone, minusone, m);
488
489 /* Exponentiate */
490 if (mm1[1] > 0UL)
491 mod_2pow_mp (r, mm1, 2UL, m);
492 else
493 mod_2pow_ul (r, mm1[0], m);
494
495 i = find_minus1 (r, minusone, po2, m);
496
497 mod_clear (r, m);
498 mod_clear (minusone, m);
499 return i;
500 }
501
502
503 int
mod_isprime(const modulus_t m)504 mod_isprime (const modulus_t m)
505 {
506 residue_t b, minusone, r1;
507 modint_t n, mm1;
508 int r = 0, po2 = 0, i;
509
510 mod_getmod_int (n, m);
511
512 if (mod_intcmp_ul (n, 1UL) == 0)
513 return 0;
514
515 if (n[0] % 2UL == 0UL)
516 {
517 r = (mod_intcmp_ul(n, 2UL) == 0);
518 #if defined(PARI)
519 printf ("isprime(%lu) == %d /* PARI */\n", n[0], r);
520 #endif
521 return r;
522 }
523
524 /* Set mm1 to the odd part of m-1 */
525 mod_intset (mm1, n);
526 mm1[0]--;
527 if (mm1[0] == 0UL)
528 {
529 mm1[0] = mm1[1];
530 mm1[1] = 0UL;
531 po2 += LONG_BIT;
532 }
533 ASSERT (mm1[0] != 0UL);
534 i = ularith_ctz (mm1[0]);
535 mod_intshr (mm1, mm1, i);
536 po2 += i;
537
538 mod_init_noset0 (b, m);
539 mod_init_noset0 (minusone, m);
540 mod_init_noset0 (r1, m);
541 mod_set1 (minusone, m);
542 mod_neg (minusone, minusone, m);
543
544 /* Do base 2 SPRP test */
545 if (mm1[1] != 0UL)
546 mod_2pow_mp (r1, mm1, 2UL, m);
547 else
548 mod_2pow_ul (r1, mm1[0], m);
549 /* If n is prime and 1 or 7 (mod 8), then 2 is a square (mod n)
550 and one less squaring must suffice. This does not strengthen the
551 test but saves one squaring for composite input */
552 if (n[0] % 8 == 7)
553 {
554 if (!mod_is1 (r1, m))
555 goto end;
556 }
557 else if (!find_minus1 (r1, minusone, po2 - ((n[0] % 8 == 7) ? 1 : 0), m))
558 goto end; /* Not prime */
559
560 /* Base 3 is poor at identifying composites == 1 (mod 3), but good at
561 identifying composites == 2 (mod 3). Thus we use it only for 2 (mod 3) */
562 i = n[0] % 3UL + n[1] % 3;
563 if (i == 1 || i == 4)
564 {
565 if (mm1[1] != 0UL)
566 mod_npow_mp (r1, 7UL, mm1, 2UL, m); /* r = 7^mm1 mod m */
567 else
568 mod_npow_ul (r1, 7UL, mm1[0], m);
569 if (!find_minus1 (r1, minusone, po2, m))
570 goto end; /* Not prime */
571
572 mod_set_ul_reduced (b, 61UL, m); /* Use addition chain? */
573 if (mm1[1] != 0UL)
574 mod_pow_mp (r1, b, mm1, 2UL, m); /* r = 61^mm1 mod m */
575 else
576 mod_pow_ul (r1, b, mm1[0], m);
577 if (!find_minus1 (r1, minusone, po2, m))
578 goto end; /* Not prime */
579
580 #if LONG_BIT == 32
581 {
582 /* These are the two base 2,7,61 SPSP below 11207066041 */
583 const modint_t
584 c4759123141 = {464155845UL, 1UL},
585 c8411807377 = {4116840081UL, 1UL},
586 c11207066041 = {2617131449UL, 2UL};
587 if (mod_intcmp (n, c11207066041) < 0)
588 {
589 r = mod_intcmp (n, c4759123141) != 0 &&
590 mod_intcmp (n, c8411807377) != 0;
591 goto end;
592 }
593 }
594 #endif
595
596 if (mm1[1] != 0UL)
597 mod_npow_mp (r1, 5UL, mm1, 2UL, m); /* r = 5^mm1 mod m */
598 else
599 mod_npow_ul (r1, 5UL, mm1[0], m);
600 if (!find_minus1 (r1, minusone, po2, m))
601 goto end; /* Not prime */
602
603 #if LONG_BIT == 32
604 {
605 /* These are the base 2,5,7,61 SPSP < 10^13 and n == 1 (mod 3) */
606 const modint_t
607 c30926647201 = {861876129UL,7UL},
608 c45821738881 = {2872065921UL,10UL},
609 c74359744201 = {1345300169UL,17UL},
610 c90528271681 = {333958465UL,21UL},
611 c110330267041 = {2956084641UL,25UL},
612 c373303331521 = {3936144065UL,86UL},
613 c440478111067 = {2391446875UL,102UL},
614 c1436309367751 = {1790290887UL,334UL},
615 c1437328758421 = {2809681557UL,334UL},
616 c1858903385041 = {3477513169UL,432UL},
617 c4897239482521 = {976765081UL,1140UL},
618 c5026103290981 = {991554661UL,1170UL},
619 c5219055617887 = {670353247UL,1215UL},
620 c5660137043641 = {3665114809UL,1317UL},
621 c6385803726241 = {3482324385UL,1486UL};
622
623 r = mod_intcmp (n, c30926647201) != 0 &&
624 mod_intcmp (n, c45821738881) != 0 &&
625 mod_intcmp (n, c74359744201) != 0 &&
626 mod_intcmp (n, c90528271681) != 0 &&
627 mod_intcmp (n, c110330267041) != 0 &&
628 mod_intcmp (n, c373303331521) != 0 &&
629 mod_intcmp (n, c440478111067) != 0 &&
630 mod_intcmp (n, c1436309367751) != 0 &&
631 mod_intcmp (n, c1437328758421) != 0 &&
632 mod_intcmp (n, c1858903385041) != 0 &&
633 mod_intcmp (n, c4897239482521) != 0 &&
634 mod_intcmp (n, c5026103290981) != 0 &&
635 mod_intcmp (n, c5219055617887) != 0 &&
636 mod_intcmp (n, c5660137043641) != 0 &&
637 mod_intcmp (n, c6385803726241) != 0;
638 }
639 #else
640 /* For 64 bit arithmetic, a two-word modulus is neccessarily too
641 large for any deterministic test (with the lists of SPSP
642 currently available). A modulus >2^64 and == 1 (mod 3) that
643 survived base 2,5,7,61 is assumed to be prime. */
644 r = 1;
645 #endif
646 }
647 else
648 {
649 /* Case n % 3 == 0, 2 */
650
651 if (mm1[1] != 0UL)
652 mod_npow_mp (r1, 3UL, mm1, 2UL, m); /* r = 3^mm1 mod m */
653 else
654 mod_npow_ul (r1, 3UL, mm1[0], m);
655 if (!find_minus1 (r1, minusone, po2, m))
656 goto end; /* Not prime */
657
658 if (mm1[1] != 0UL)
659 mod_npow_mp (r1, 5UL, mm1, 2UL, m); /* r = 5^mm1 mod m */
660 else
661 mod_npow_ul (r1, 5UL, mm1[0], m);
662 if (!find_minus1 (r1, minusone, po2, m))
663 goto end; /* Not prime */
664
665 #if LONG_BIT == 32
666 {
667 /* These are the base 2,3,5 SPSP < 10^13 and n == 2 (mod 3) */
668 const modint_t
669 c244970876021 = {157740149UL,57UL},
670 c405439595861 = {1712670037UL,94UL},
671 c1566655993781 = {3287898037UL,364UL},
672 c3857382025841 = {501394033UL,898UL},
673 c4074652846961 = {3023850353UL,948UL},
674 c5783688565841 = {2662585425UL,1346UL};
675
676 r = mod_intcmp (n, c244970876021) != 0 &&
677 mod_intcmp (n, c405439595861) != 0 &&
678 mod_intcmp (n, c1566655993781) != 0 &&
679 mod_intcmp (n, c3857382025841) != 0 &&
680 mod_intcmp (n, c4074652846961) != 0 &&
681 mod_intcmp (n, c5783688565841) != 0;
682 }
683 #else
684 r = 1;
685 #endif
686 }
687
688 end:
689 #if defined(PARI)
690 printf ("isprime(%lu) == %d /* PARI */\n", n, r);
691 #endif
692 mod_clear (b, m);
693 mod_clear (minusone, m);
694 mod_clear (r1, m);
695 return r;
696 }
697
698
699
700 int
mod_jacobi(const residue_t a_par,const modulus_t m_par)701 mod_jacobi (const residue_t a_par, const modulus_t m_par)
702 {
703 modint_t a, m, s;
704 int t = 1;
705
706 mod_get_int (a, a_par, m_par);
707 mod_getmod_int (m, m_par);
708
709 while (!mod_intequal_ul (a, 0UL))
710 {
711 while (a[0] % 2UL == 0UL) /* TODO speedup */
712 {
713 mod_intshr (a, a, 1);
714 if (m[0] % 8UL == 3UL || m[0] % 8UL == 5UL)
715 t = -t;
716 }
717 mod_intset (s, a); /* swap a and m */
718 mod_intset (a, m);
719 mod_intset (m, s);
720 if (a[0] % 4UL == 3UL && m[0] % 4UL == 3UL)
721 t = -t;
722
723 /* m is odd here */
724 if (mod_intcmp (a, m) >= 0)
725 mod_intmod (a, a, m);
726 }
727 if (m[1] != 0UL || m[0] != 1UL)
728 t = 0;
729
730 return t;
731 }
732