Line data Source code
1 : /* Copyright (C) 2000 The PARI group.
2 :
3 : This file is part of the PARI/GP package.
4 :
5 : PARI/GP is free software; you can redistribute it and/or modify it under the
6 : terms of the GNU General Public License as published by the Free Software
7 : Foundation; either version 2 of the License, or (at your option) any later
8 : version. It is distributed in the hope that it will be useful, but WITHOUT
9 : ANY WARRANTY WHATSOEVER.
10 :
11 : Check the License for details. You should have received a copy of it, along
12 : with the package; see the file 'COPYING'. If not, write to the Free Software
13 : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
14 :
15 : #include "pari.h"
16 : #include "paripriv.h"
17 :
18 : #define DEBUGLEVEL DEBUGLEVEL_isprime
19 :
20 : /*********************************************************************/
21 : /** **/
22 : /** PSEUDO PRIMALITY (MILLER-RABIN) **/
23 : /** **/
24 : /*********************************************************************/
25 : typedef struct {
26 : GEN n, sqrt1, sqrt2, t1, t;
27 : long r1;
28 : } MR_Jaeschke_t;
29 :
30 : static void
31 411 : init_MR_Jaeschke(MR_Jaeschke_t *S, GEN n)
32 : {
33 411 : S->n = n = absi_shallow(n);
34 411 : S->t = subiu(n,1);
35 411 : S->r1 = vali(S->t);
36 411 : S->t1 = shifti(S->t, -S->r1);
37 411 : S->sqrt1 = cgeti(lg(n)); S->sqrt1[1] = evalsigne(0)|evallgefint(2);
38 411 : S->sqrt2 = cgeti(lg(n)); S->sqrt2[1] = evalsigne(0)|evallgefint(2);
39 411 : }
40 :
41 : /* is n strong pseudo-prime for base a ? 'End matching' (check for square
42 : * roots of -1): if ends do mismatch, then we have factored n, and this
43 : * information should be made available to the factoring machinery. But so
44 : * exceedingly rare... besides we use BSPW now. */
45 : static int
46 721 : ispsp(MR_Jaeschke_t *S, ulong a)
47 : {
48 721 : pari_sp av = avma;
49 721 : GEN c = Fp_pow(utoipos(a), S->t1, S->n);
50 : long r;
51 :
52 721 : if (is_pm1(c) || equalii(S->t, c)) return 1;
53 : /* go fishing for -1, not for 1 (saves one squaring) */
54 768 : for (r = S->r1 - 1; r; r--) /* r1 - 1 squarings */
55 : {
56 667 : GEN c2 = c;
57 667 : c = remii(sqri(c), S->n);
58 667 : if (equalii(S->t, c))
59 : {
60 278 : if (!signe(S->sqrt1))
61 : {
62 177 : affii(subii(S->n, c2), S->sqrt2);
63 177 : affii(c2, S->sqrt1); return 1;
64 : }
65 : /* saw one earlier: too many sqrt(-1)s mod n ? */
66 101 : return equalii(c2, S->sqrt1) || equalii(c2, S->sqrt2);
67 : }
68 389 : if (gc_needed(av,1))
69 : {
70 0 : if(DEBUGMEM>1) pari_warn(warnmem,"ispsp, r = %ld", r);
71 0 : c = gerepileuptoint(av, c);
72 : }
73 : }
74 101 : return 0;
75 : }
76 :
77 : /* is n > 0 strong pseudo-prime for base 2 ? Only used when lgefint(n) > 3,
78 : * so don't test */
79 : static int
80 131705 : is2psp(GEN n)
81 : {
82 131705 : GEN c, t = subiu(n, 1);
83 131702 : pari_sp av = avma;
84 131702 : long e = vali(t);
85 :
86 131702 : c = Fp_pow(gen_2, shifti(t, -e), n);
87 131703 : if (is_pm1(c) || equalii(t, c)) return 1;
88 249396 : while (--e)
89 : { /* go fishing for -1, not for 1 (e - 1 squaring) */
90 141021 : c = remii(sqri(c), n);
91 141022 : if (equalii(t, c)) return 1;
92 : /* can return 0 if (c == 1) but very infrequent */
93 130065 : if (gc_needed(av,1))
94 : {
95 0 : if(DEBUGMEM>1) pari_warn(warnmem,"is2psp, r = %ld", e);
96 0 : c = gerepileuptoint(av, c);
97 : }
98 : }
99 108375 : return 0;
100 : }
101 : static int
102 5788309 : uispsp_pre(ulong a, ulong n, ulong ni)
103 : {
104 5788309 : ulong c, t = n - 1;
105 5788309 : long e = vals(t);
106 :
107 5789001 : c = Fl_powu_pre(a, t >> e, n, ni);
108 5776881 : if (c == 1 || c == t) return 1;
109 10419835 : while (--e)
110 : { /* go fishing for -1, not for 1 (saves one squaring) */
111 5523288 : c = Fl_sqr_pre(c, n, ni);
112 5531100 : if (c == t) return 1;
113 : /* can return 0 if (c == 1) but very infrequent */
114 : }
115 4896547 : return 0;
116 : }
117 : int
118 7409744 : uispsp(ulong a, ulong n)
119 : {
120 : ulong c, t;
121 : long e;
122 :
123 7409744 : if (n & HIGHMASK) return uispsp_pre(a, n, get_Fl_red(n));
124 7216392 : t = n - 1;
125 7216392 : e = vals(t);
126 7232733 : c = Fl_powu(a, t >> e, n);
127 7278840 : if (c == 1 || c == t) return 1;
128 8069361 : while (--e)
129 : { /* go fishing for -1, not for 1 (e - 1 squaring) */
130 5623791 : c = Fl_sqr(c, n);
131 5625262 : if (c == t) return 1;
132 : /* can return 0 if (c == 1) but very infrequent */
133 : }
134 2445570 : return 0;
135 : }
136 : int
137 0 : uis2psp(ulong n) { return uispsp(2, n); }
138 :
139 : /* Miller-Rabin test for k random bases */
140 : long
141 28 : millerrabin(GEN n, long k)
142 : {
143 28 : pari_sp av2, av = avma;
144 : ulong r;
145 : long i;
146 : MR_Jaeschke_t S;
147 :
148 28 : if (typ(n) != t_INT) pari_err_TYPE("millerrabin",n);
149 28 : if (signe(n) <= 0) return 0;
150 : /* If |n| <= 3, check if n = +- 1 */
151 28 : if (lgefint(n) == 3 && uel(n,2) <= 3) return uel(n,2) != 1;
152 :
153 14 : if (!mod2(n)) return 0;
154 7 : init_MR_Jaeschke(&S, n); av2 = avma;
155 21 : for (i = 1; i <= k; i++)
156 : {
157 20 : do r = umodui(pari_rand(), n); while (!r);
158 14 : if (DEBUGLEVEL > 4) err_printf("Miller-Rabin: testing base %ld\n", r);
159 14 : if (!ispsp(&S, r)) return gc_long(av,0);
160 14 : set_avma(av2);
161 : }
162 7 : return gc_long(av,1);
163 : }
164 :
165 : GEN
166 14 : gispseudoprime(GEN x, long flag)
167 14 : { return flag? map_proto_lGL(millerrabin, x, flag): map_proto_lG(BPSW_psp,x); }
168 :
169 : long
170 0 : ispseudoprime(GEN x, long flag)
171 0 : { return flag? millerrabin(x, flag): BPSW_psp(x); }
172 :
173 : int
174 4432 : MR_Jaeschke(GEN n)
175 : {
176 : MR_Jaeschke_t S;
177 : pari_sp av;
178 :
179 4432 : if (lgefint(n) == 3) return uisprime(uel(n,2));
180 404 : if (!mod2(n)) return 0;
181 404 : av = avma; init_MR_Jaeschke(&S, n);
182 404 : return gc_int(av, ispsp(&S, 31) && ispsp(&S, 73));
183 : }
184 :
185 : /*********************************************************************/
186 : /** **/
187 : /** PSEUDO PRIMALITY (LUCAS) **/
188 : /** **/
189 : /*********************************************************************/
190 : /* compute n-th term of Lucas sequence modulo N.
191 : * v_{k+2} = P v_{k+1} - v_k, v_0 = 2, v_1 = P.
192 : * Assume n > 0 */
193 : static GEN
194 23330 : LucasMod(GEN n, ulong P, GEN N)
195 : {
196 23330 : pari_sp av = avma;
197 23330 : GEN nd = int_MSW(n);
198 23330 : ulong m = *nd;
199 : long i, j;
200 23330 : GEN v = utoipos(P), v1 = utoipos(P*P - 2);
201 :
202 23330 : if (m == 1)
203 1845 : j = 0;
204 : else
205 : {
206 21485 : j = 1+bfffo(m); /* < BIL */
207 21485 : m <<= j; j = BITS_IN_LONG - j;
208 : }
209 23330 : for (i=lgefint(n)-2;;) /* cf. leftright_pow */
210 : {
211 2457080 : for (; j; m<<=1,j--)
212 : { /* v = v_k, v1 = v_{k+1} */
213 2400344 : if (m&HIGHBIT)
214 : { /* set v = v_{2k+1}, v1 = v_{2k+2} */
215 737084 : v = subiu(mulii(v,v1), P);
216 736755 : v1= subiu(sqri(v1), 2);
217 : }
218 : else
219 : {/* set v = v_{2k}, v1 = v_{2k+1} */
220 1663260 : v1= subiu(mulii(v,v1), P);
221 1662880 : v = subiu(sqri(v), 2);
222 : }
223 2399723 : v = modii(v, N);
224 2399883 : v1= modii(v1,N);
225 2399906 : if (gc_needed(av,1))
226 : {
227 0 : if(DEBUGMEM>1) pari_warn(warnmem,"LucasMod");
228 0 : gerepileall(av, 2, &v,&v1);
229 : }
230 : }
231 56736 : if (--i == 0) return v;
232 33848 : j = BITS_IN_LONG;
233 33848 : nd=int_precW(nd);
234 33848 : m = *nd;
235 : }
236 : }
237 : /* compute n-th term of Lucas sequence modulo N.
238 : * v_{k+2} = P v_{k+1} - v_k, v_0 = 2, v_1 = P.
239 : * Assume n > 0 */
240 : static ulong
241 752645 : u_LucasMod_pre(ulong n, ulong P, ulong N, ulong NI)
242 : {
243 : ulong v, v1, m;
244 : long j;
245 :
246 752645 : if (n == 1) return P;
247 752633 : j = 1 + bfffo(n); /* < BIL */
248 752633 : v = P; v1 = P*P - 2;
249 752633 : m = n<<j; j = BITS_IN_LONG - j;
250 45993966 : for (; j; m<<=1,j--)
251 : { /* v = v_k, v1 = v_{k+1} */
252 45258455 : if (m & HIGHBIT)
253 : { /* set v = v_{2k+1}, v1 = v_{2k+2} */
254 4622419 : v = Fl_sub(Fl_mul_pre(v,v1,N,NI), P, N);
255 4622463 : v1= Fl_sub(Fl_sqr_pre(v1,N,NI), 2UL, N);
256 : }
257 : else
258 : {/* set v = v_{2k}, v1 = v_{2k+1} */
259 40636036 : v1= Fl_sub(Fl_mul_pre(v,v1,N,NI),P, N);
260 40636149 : v = Fl_sub(Fl_sqr_pre(v,N,NI), 2UL, N);
261 : }
262 : }
263 735511 : return v;
264 : }
265 :
266 : /* !(N & HIGHMASK) */
267 : static ulong
268 0 : u_LucasMod(ulong n, ulong P, ulong N)
269 : {
270 : ulong v, v1, m;
271 : long j;
272 :
273 0 : if (n == 1) return P;
274 0 : j = 1 + bfffo(n); /* < BIL */
275 0 : v = P; v1 = P*P - 2;
276 0 : m = n<<j; j = BITS_IN_LONG - j;
277 0 : for (; j; m<<=1,j--)
278 : { /* v = v_k, v1 = v_{k+1} */
279 0 : if (m & HIGHBIT)
280 : { /* set v = v_{2k+1}, v1 = v_{2k+2} */
281 0 : v = Fl_sub((v*v1) % N, P, N);
282 0 : v1= Fl_sub((v1*v1)% N, 2UL, N);
283 : }
284 : else
285 : {/* set v = v_{2k}, v1 = v_{2k+1} */
286 0 : v1= Fl_sub((v*v1) % N ,P, N);
287 0 : v = Fl_sub((v*v) % N, 2UL, N);
288 : }
289 : }
290 0 : return v;
291 : }
292 :
293 : static ulong
294 752407 : get_disc(ulong n)
295 : {
296 : ulong b;
297 : long i;
298 752407 : for (b = 3, i = 0;; b += 2, i++)
299 871110 : {
300 1623517 : ulong c = b*b - 4; /* = 1 mod 4 */
301 1623517 : if (krouu(n % c, c) < 0) break;
302 871110 : if (i == 64 && uissquareall(n, &c)) return 0; /* oo loop if N = m^2 */
303 : }
304 752654 : return b;
305 : }
306 : static int
307 752412 : uislucaspsp_pre(ulong n, ulong ni)
308 : {
309 : long i, v;
310 752412 : ulong b, z, m = n + 1;
311 :
312 752412 : if (!m) return 0; /* neither 2^32-1 nor 2^64-1 are Lucas-pp */
313 752412 : b = get_disc(n); if (!b) return 0;
314 752651 : v = vals(m); m >>= v;
315 752644 : z = u_LucasMod_pre(m, b, n, ni);
316 752775 : if (z == 2 || z == n-2) return 1;
317 658905 : for (i=1; i<v; i++)
318 : {
319 658904 : if (!z) return 1;
320 347963 : z = Fl_sub(Fl_sqr_pre(z,n,ni), 2UL, n);
321 347964 : if (z == 2) return 0;
322 : }
323 5 : return 0;
324 : }
325 : int
326 0 : uislucaspsp(ulong n)
327 : {
328 : long i, v;
329 : ulong b, z, m;
330 :
331 0 : if (n & HIGHMASK) return uislucaspsp_pre(n, get_Fl_red(n));
332 0 : m = n + 1;
333 0 : if (!m) return 0; /* neither 2^32-1 nor 2^64-1 are Lucas-pp */
334 0 : b = get_disc(n); if (!b) return 0;
335 0 : v = vals(m); m >>= v;
336 0 : z = u_LucasMod(m, b, n);
337 0 : if (z == 2 || z == n-2) return 1;
338 0 : for (i=1; i<v; i++)
339 : {
340 0 : if (!z) return 1;
341 0 : z = Fl_sub((z*z) % n, 2UL, n);
342 0 : if (z == 2) return 0;
343 : }
344 0 : return 0;
345 : }
346 : /* N > 3. Caller should check that N is not a square first (taken care of here,
347 : * but inefficient) */
348 : static int
349 23330 : islucaspsp(GEN N)
350 : {
351 23330 : pari_sp av = avma;
352 : GEN m, z;
353 : long i, v;
354 : ulong b;
355 :
356 23330 : for (b=3;; b+=2)
357 27795 : {
358 51125 : ulong c = b*b - 4; /* = 1 mod 4 */
359 51125 : if (b == 129 && Z_issquare(N)) return 0; /* avoid oo loop if N = m^2 */
360 51125 : if (krouu(umodiu(N,c), c) < 0) break;
361 : }
362 23330 : m = addiu(N,1); v = vali(m); m = shifti(m,-v);
363 23330 : z = LucasMod(m, b, N);
364 23330 : if (absequaliu(z, 2)) return 1;
365 18109 : if (equalii(z, subiu(N,2))) return 1;
366 19737 : for (i=1; i<v; i++)
367 : {
368 19533 : if (!signe(z)) return 1;
369 11345 : z = modii(subiu(sqri(z), 2), N);
370 11345 : if (absequaliu(z, 2)) return 0;
371 11345 : if (gc_needed(av,1))
372 : {
373 0 : if(DEBUGMEM>1) pari_warn(warnmem,"islucaspsp");
374 0 : z = gerepileupto(av, z);
375 : }
376 : }
377 204 : return 0;
378 : }
379 :
380 : /* composite strong 2-pseudoprime < 1016801 whose prime divisors are > 101.
381 : * All have a prime divisor <= 661 */
382 : static int
383 3766 : is_2_prp_101(ulong n)
384 : {
385 3766 : switch(n) {
386 0 : case 42799:
387 : case 49141:
388 : case 88357:
389 : case 90751:
390 : case 104653:
391 : case 130561:
392 : case 196093:
393 : case 220729:
394 : case 253241:
395 : case 256999:
396 : case 271951:
397 : case 280601:
398 : case 357761:
399 : case 390937:
400 : case 458989:
401 : case 486737:
402 : case 489997:
403 : case 514447:
404 : case 580337:
405 : case 741751:
406 : case 838861:
407 : case 873181:
408 : case 877099:
409 : case 916327:
410 : case 976873:
411 0 : case 983401: return 1;
412 3766 : } return 0;
413 : }
414 :
415 : static int
416 7376334 : _uispsp(ulong a, long n) { a %= n; return !a || uispsp(a, n); }
417 : static int
418 12256554 : _uisprime(ulong n)
419 : {
420 : #ifdef LONG_IS_64BIT
421 12091550 : if (n < 341531)
422 5365379 : return _uispsp(9345883071009581737UL, n);
423 6726171 : if (n < 1050535501)
424 1098179 : return _uispsp(336781006125UL, n)
425 1098179 : && _uispsp(9639812373923155UL, n);
426 5627992 : if (n < 350269456337)
427 36325 : return _uispsp(4230279247111683200UL, n)
428 14878 : && _uispsp(14694767155120705706UL, n)
429 51203 : && _uispsp(16641139526367750375UL, n);
430 5591667 : if (n & HIGHMASK)
431 : {
432 5591691 : ulong ni = get_Fl_red(n);
433 5594091 : return (uispsp_pre(2, n, ni) && uislucaspsp_pre(n,ni));
434 : }
435 0 : return uispsp(2, n) && uislucaspsp(n);
436 : #else
437 165004 : if (n < 360018361) return _uispsp(1143370UL, n) && _uispsp(2350307676UL, n);
438 14183 : return uispsp(15, n) && uispsp(176006322UL, n) && _uispsp(4221622697UL, n);
439 : #endif
440 : }
441 :
442 : int
443 65123873 : uisprime(ulong n)
444 : {
445 65123873 : if (n < 103)
446 7105607 : switch(n)
447 : {
448 5395085 : case 2:
449 : case 3:
450 : case 5:
451 : case 7:
452 : case 11:
453 : case 13:
454 : case 17:
455 : case 19:
456 : case 23:
457 : case 29:
458 : case 31:
459 : case 37:
460 : case 41:
461 : case 43:
462 : case 47:
463 : case 53:
464 : case 59:
465 : case 61:
466 : case 67:
467 : case 71:
468 : case 73:
469 : case 79:
470 : case 83:
471 : case 89:
472 : case 97:
473 5395085 : case 101: return 1;
474 1710522 : default: return 0;
475 : }
476 : /* gcd-extraction is much slower */
477 91773605 : return odd(n) && n % 3 && n % 5 && n % 7 && n % 11 && n % 13 && n % 17
478 17349824 : && n % 19 && n % 23 && n % 29 && n % 31 && n % 37 && n % 41
479 91775808 : && (n < 1849 || _uisprime(n));
480 : }
481 :
482 : /* assume no prime divisor <= 101 */
483 : int
484 16487 : uisprime_101(ulong n)
485 : {
486 16487 : if (n < 1016801) return n < 10609? 1: (uispsp(2, n) && !is_2_prp_101(n));
487 12708 : return _uisprime(n);
488 : }
489 :
490 : /* assume no prime divisor <= 661 */
491 : int
492 64144 : uisprime_661(ulong n)
493 : {
494 64144 : if (n < 1016801) return n < 452929? 1: uispsp(2, n);
495 53221 : return _uisprime(n);
496 : }
497 :
498 : static int
499 343201 : iu_coprime(GEN N, ulong u) { return ugcd(u, umodiu(N,u)) == 1; }
500 : long
501 48360487 : BPSW_psp(GEN N)
502 : {
503 : pari_sp av;
504 48360487 : if (typ(N) != t_INT) pari_err_TYPE("BPSW_psp",N);
505 48360955 : if (signe(N) <= 0) return 0;
506 48360955 : if (lgefint(N) == 3) return uisprime(uel(N,2));
507 186907 : if (!mod2(N)) return 0;
508 : #ifdef LONG_IS_64BIT
509 : /* 16294579238595022365 = 3*5*7*11*13*17*19*23*29*31*37*41*43*47*53
510 : * 7145393598349078859 = 59*61*67*71*73*79*83*89*97*101 */
511 154773 : if (!iu_coprime(N, 16294579238595022365UL) ||
512 98950 : !iu_coprime(N, 7145393598349078859UL)) return 0;
513 : #else
514 : /* 4127218095 = 3*5*7*11*13*17*19*23*37
515 : * 3948078067 = 29*31*41*43*47*53
516 : * 4269855901 = 59*83*89*97*101
517 : * 1673450759 = 61*67*71*73*79 */
518 109958 : if (!iu_coprime(N, 4127218095UL) ||
519 86177 : !iu_coprime(N, 3948078067UL) ||
520 78389 : !iu_coprime(N, 1673450759UL) ||
521 66031 : !iu_coprime(N, 4269855901UL)) return 0;
522 : #endif
523 : /* no prime divisor < 103 */
524 92047 : av = avma;
525 92047 : return gc_long(av, is2psp(N) && islucaspsp(N));
526 : }
527 :
528 : /* can we write n = x^k ? Assume N has no prime divisor <= 2^14.
529 : * Not memory clean */
530 : long
531 31830 : isanypower_nosmalldiv(GEN N, GEN *px)
532 : {
533 31830 : GEN x = N, y;
534 31830 : ulong mask = 7;
535 31830 : long ex, k = 1;
536 : forprime_t T;
537 46867 : while (Z_issquareall(x, &y)) { k <<= 1; x = y; }
538 32131 : while ( (ex = is_357_power(x, &y, &mask)) ) { k *= ex; x = y; }
539 31829 : (void)u_forprime_init(&T, 11, ULONG_MAX);
540 : /* stop when x^(1/k) < 2^14 */
541 31828 : while ( (ex = is_pth_power(x, &y, &T, 15)) ) { k *= ex; x = y; }
542 31830 : *px = x; return k;
543 : }
544 :
545 : /* no prime divisor <= 2^14 (> 661) */
546 : long
547 66706 : BPSW_psp_nosmalldiv(GEN N)
548 : {
549 : pari_sp av;
550 66706 : long l = lgefint(N);
551 :
552 66706 : if (l == 3) return uisprime_661(uel(N,2));
553 41669 : av = avma;
554 : /* N large: test for pure power, rarely succeeds, but requires < 1% of
555 : * compositeness test times */
556 41669 : if (bit_accuracy(l) > 512 && isanypower_nosmalldiv(N,&N) != 1)
557 2028 : return gc_long(av,0);
558 39641 : N = absi_shallow(N);
559 39641 : return gc_long(av, is2psp(N) && islucaspsp(N));
560 : }
561 :
562 : /***********************************************************************/
563 : /** **/
564 : /** Pocklington-Lehmer **/
565 : /** P-1 primality test **/
566 : /** **/
567 : /***********************************************************************/
568 : /* Assume x BPSW pseudoprime. Check whether it's small enough to be certified
569 : * prime (< 2^64). Reference for strong 2-pseudoprimes:
570 : * http://www.cecm.sfu.ca/Pseudoprimes/index-2-to-64.html */
571 : static int
572 4116385 : BPSW_isprime_small(GEN x)
573 : {
574 4116385 : long l = lgefint(x);
575 : #ifdef LONG_IS_64BIT
576 4050473 : return (l == 3);
577 : #else
578 65912 : return (l <= 4);
579 : #endif
580 : }
581 :
582 : /* Assume N > 1, p^e || N-1, p prime. Find a witness a(p) such that
583 : * a^(N-1) = 1 (mod N)
584 : * a^(N-1)/p - 1 invertible mod N.
585 : * Proves that any divisor of N is 1 mod p^e. Return 0 if N is composite */
586 : static ulong
587 13542 : pl831(GEN N, GEN p)
588 : {
589 13542 : GEN b, c, g, Nmunp = diviiexact(subiu(N,1), p);
590 13542 : pari_sp av = avma;
591 : ulong a;
592 19916 : for(a = 2;; a++, set_avma(av))
593 : {
594 19916 : b = Fp_pow(utoipos(a), Nmunp, N);
595 19916 : if (!equali1(b)) break;
596 : }
597 13542 : c = Fp_pow(b,p,N);
598 13542 : g = gcdii(subiu(b,1), N); /* 0 < g < N */
599 13542 : return (equali1(c) && equali1(g))? a: 0;
600 : }
601 :
602 : /* Brillhart, Lehmer, Selfridge test (Crandall & Pomerance, Th 4.1.5)
603 : * N^(1/3) <= f fully factored, f | N-1. If pl831(p) is true for
604 : * any prime divisor p of f, then any divisor of N is 1 mod f.
605 : * In that case return 1 iff N is prime */
606 : static int
607 49 : BLS_test(GEN N, GEN f)
608 : {
609 : GEN c1, c2, r, q;
610 49 : q = dvmdii(N, f, &r);
611 49 : if (!is_pm1(r)) return 0;
612 49 : c2 = dvmdii(q, f, &c1);
613 : /* N = 1 + f c1 + f^2 c2, 0 <= c_i < f; check whether it is of the form
614 : * (1 + fa)(1 + fb) */
615 49 : return ! Z_issquare(subii(sqri(c1), shifti(c2,2)));
616 : }
617 :
618 : /* BPSW_psp(N) && !BPSW_isprime_small(N). Decide between Pocklington-Lehmer
619 : * and APRCL/ECPP. Return a vector of (small) primes such that PL-witnesses
620 : * guarantee the primality of N. Return NULL if PL is likely too expensive.
621 : * Return gen_0 if BLS test finds N to be composite */
622 : static GEN
623 4429 : BPSW_try_PL(GEN N)
624 : {
625 4429 : ulong B = minuu(1UL<<19, maxprime());
626 4429 : GEN E, p, U, F, N_1 = subiu(N,1);
627 4429 : GEN fa = absZ_factor_limit_strict(N_1, B, &U), P = gel(fa,1);
628 :
629 4429 : if (!U) return P; /* N-1 fully factored */
630 1721 : p = gel(U,1);
631 1721 : E = gel(fa,2);
632 1721 : U = powii(p, gel(U,2)); /* unfactored part of N-1 */
633 1721 : F = (lg(P) == 2)? powii(gel(P,1), gel(E,1)): diviiexact(N_1, U);
634 :
635 : /* N-1 = F U, F factored, U composite, (U,F) = 1 */
636 1721 : if (cmpii(F, U) > 0) return P; /* 1/2-smooth */
637 1714 : if (cmpii(sqri(F), U) > 0) return BLS_test(N,F)? P: gen_0; /* 1/3-smooth */
638 1672 : return NULL; /* not smooth enough */
639 : }
640 :
641 : static GEN isprimePL(GEN N);
642 : /* F is a vector whose entries are primes. For each of them, find a PL
643 : * witness. Return 0 if caller lied and F contains a composite */
644 : static long
645 2757 : PL_certify(GEN N, GEN F)
646 : {
647 2757 : long i, l = lg(F);
648 16222 : for(i = 1; i < l; i++)
649 13465 : if (! pl831(N, gel(F,i))) return 0;
650 2757 : return 1;
651 : }
652 : /* F is a vector whose entries are *believed* to be primes (BPSW_psp).
653 : * For each of them, recording a witness and recursive primality certificate */
654 : static GEN
655 98 : PL_certificate(GEN N, GEN F)
656 : {
657 98 : long i, l = lg(F);
658 : GEN C;
659 98 : if (BPSW_isprime_small(N)) return N;
660 98 : C = cgetg(l, t_VEC);
661 574 : for (i = 1; i < l; i++)
662 : {
663 476 : GEN p = gel(F,i), C0;
664 : ulong w;
665 :
666 476 : if (BPSW_isprime_small(p)) { gel(C,i) = p; continue; }
667 77 : w = pl831(N,p); if (!w) return gen_0;
668 77 : C0 = isprimePL(p);
669 77 : if (isintzero(C0))
670 : { /* composite in prime factorisation ! */
671 0 : err_printf("Not a prime: %Ps", p);
672 0 : pari_err_BUG("PL_certificate [false prime number]");
673 : }
674 77 : gel(C,i) = mkvec3(p,utoipos(w), C0);
675 : }
676 98 : return mkvec2(N, C);
677 : }
678 : /* M a t_MAT */
679 : static int
680 84 : PL_isvalid(GEN v)
681 : {
682 : GEN C, F, F2, N, N1, U;
683 : long i, l;
684 84 : switch(typ(v))
685 : {
686 0 : case t_INT: return BPSW_isprime_small(v) && BPSW_psp(v);
687 84 : case t_VEC: if (lg(v) == 3) break;/*fall through */
688 0 : default: return 0;
689 : }
690 84 : N = gel(v,1);
691 84 : C = gel(v,2);
692 84 : if (typ(N) != t_INT || signe(N) <= 0 || typ(C) != t_VEC) return 0;
693 84 : U = N1 = subiu(N,1);
694 84 : l = lg(C);
695 427 : for (i = 1; i < l; i++)
696 : {
697 350 : GEN p = gel(C,i), a = NULL, C0 = NULL, ap;
698 : long vp;
699 350 : if (typ(p) != t_INT)
700 : {
701 70 : if (lg(p) != 4) return 0;
702 70 : a = gel(p,2); C0 = gel(p,3); p = gel(p,1);
703 70 : if (typ(p) != t_INT || typ(a) != t_INT || !PL_isvalid(C0)) return 0;
704 : }
705 350 : vp = Z_pvalrem(U, p, &U); if (!vp) return 0;
706 343 : if (!a)
707 : {
708 280 : if (!BPSW_isprime_small(p)) return 0;
709 280 : continue;
710 : }
711 63 : if (!equalii(gel(C0,1), p)) return 0;
712 63 : ap = Fp_pow(a, diviiexact(N1,p), N);
713 63 : if (!equali1(gcdii(subiu(ap,1), N)) || !equali1(Fp_pow(ap, p, N))) return 0;
714 : }
715 77 : F = diviiexact(N1, U); /* factored part of N-1 */
716 77 : F2= sqri(F);
717 77 : if (cmpii(F2, N) > 0) return 1;
718 0 : if (cmpii(mulii(F2,F), N) <= 0) return 0;
719 0 : return BLS_test(N,F);
720 : }
721 :
722 : /* Assume N is a strong BPSW pseudoprime, Pocklington-Lehmer primality proof.
723 : * Return gen_0 (nonprime), N (small prime), matrix (large prime)
724 : *
725 : * The matrix has 3 columns, [a,b,c] with
726 : * a[i] prime factor of N-1,
727 : * b[i] witness for a[i] as in pl831
728 : * c[i] check_prime(a[i]) */
729 : static GEN
730 119 : isprimePL(GEN N)
731 : {
732 : GEN cbrtN, N_1, F, f;
733 119 : if (BPSW_isprime_small(N)) return N;
734 98 : cbrtN = sqrtnint(N,3);
735 98 : N_1 = subiu(N,1);
736 98 : F = Z_factor_until(N_1, sqri(cbrtN));
737 98 : f = factorback(F); /* factored part of N-1, f^3 > N */
738 98 : if (DEBUGLEVEL>3)
739 : {
740 0 : GEN r = divri(itor(f,LOWDEFAULTPREC), N);
741 0 : err_printf("Pocklington-Lehmer: proving primality of N = %Ps\n", N);
742 0 : err_printf("Pocklington-Lehmer: N-1 factored up to %Ps! (%.3Ps%%)\n", f, r);
743 : }
744 : /* if N-1 is only N^(1/3)-smooth, BLS test */
745 98 : if (!equalii(f,N_1) && cmpii(sqri(f),N) <= 0 && !BLS_test(N,f))
746 0 : return gen_0; /* Failed, N is composite */
747 98 : F = gel(F,1); settyp(F, t_VEC);
748 98 : return PL_certificate(N, F);
749 : }
750 :
751 : /* assume N a BPSW pseudoprime, in particular, it is odd > 2. Prove N prime */
752 : long
753 4116377 : BPSW_isprime(GEN N)
754 : {
755 : pari_sp av;
756 : long t;
757 : GEN P;
758 4116377 : if (BPSW_isprime_small(N)) return 1;
759 4379 : av = avma; P = BPSW_try_PL(N);
760 4429 : if (!P) /* not smooth enough */
761 1672 : t = expi(N) < 768? isprimeAPRCL(N): isprimeECPP(N);
762 : else
763 2757 : t = (typ(P) == t_INT)? 0: PL_certify(N,P);
764 4429 : return gc_long(av,t);
765 : }
766 :
767 : static long
768 42 : _isprimePL(GEN x)
769 : {
770 42 : pari_sp av = avma;
771 42 : if (!BPSW_psp(x)) return 0;
772 35 : return gc_long(av, !isintzero(isprimePL(x)));
773 : }
774 : GEN
775 39086513 : gisprime(GEN x, long flag)
776 : {
777 39086513 : switch (flag)
778 : {
779 39087029 : case 0: return map_proto_lG(isprime,x);
780 28 : case 1: return map_proto_lG(_isprimePL,x);
781 14 : case 2: return map_proto_lG(isprimeAPRCL,x);
782 21 : case 3: return map_proto_lG(isprimeECPP,x);
783 : }
784 0 : pari_err_FLAG("gisprime");
785 : return NULL;/*LCOV_EXCL_LINE*/
786 : }
787 :
788 : long
789 39786712 : isprime(GEN x) { return BPSW_psp(x) && BPSW_isprime(x); }
790 :
791 : GEN
792 84 : primecert0(GEN x, long flag, long stopat)
793 : {
794 84 : if ((flag || typ(x) == t_INT) && !BPSW_psp(x)) return gen_0;
795 77 : switch(flag)
796 : {
797 70 : case 0: return ecpp0(x, stopat);
798 7 : case 1: { pari_sp av = avma; return gerepilecopy(av, isprimePL(x)); }
799 : }
800 0 : pari_err_FLAG("primecert");
801 : return NULL;/*LCOV_EXCL_LINE*/
802 : }
803 :
804 : GEN
805 0 : primecert(GEN x, long flag)
806 0 : { return primecert0(x, flag, 0); }
807 :
808 : enum { c_VOID = 0, c_ECPP, c_N1 };
809 : static long
810 98 : cert_type(GEN c)
811 : {
812 98 : switch(typ(c))
813 : {
814 0 : case t_INT: return c_ECPP;
815 98 : case t_VEC:
816 98 : if (lg(c) == 3 && typ(gel(c,1)) == t_INT) return c_N1;
817 84 : return c_ECPP;
818 : }
819 0 : return c_VOID;
820 : }
821 :
822 : long
823 56 : primecertisvalid(GEN c)
824 : {
825 56 : switch(typ(c))
826 : {
827 7 : case t_INT: return BPSW_isprime_small(c) && BPSW_psp(c);
828 49 : case t_VEC:
829 49 : return cert_type(c) == c_ECPP? ecppisvalid(c): PL_isvalid(c);
830 : }
831 0 : return 0;
832 : }
833 :
834 : static long
835 392 : check_ecppcertentry(GEN c)
836 : {
837 : GEN v;
838 392 : long i,l = lg(c);
839 392 : if (typ(c)!=t_VEC || l!=6) return 0;
840 1925 : for(i=1; i<=4; i++)
841 1540 : if (typ(gel(c,i))!=t_INT) return 0;
842 385 : v = gel(c,5);
843 385 : if(typ(v)!=t_VEC) return 0;
844 1155 : for(i=1; i<=2; i++)
845 770 : if (typ(gel(v,i))!=t_INT) return 0;
846 385 : return 1;
847 : }
848 :
849 : long
850 56 : check_ecppcert(GEN c)
851 : {
852 : long i, l;
853 56 : switch(typ(c))
854 : {
855 0 : case t_INT: return signe(c) >= 0;
856 56 : case t_VEC: break;
857 0 : default: return 0;
858 : }
859 56 : l = lg(c); if (l == 1) return 0;
860 434 : for (i = 1; i < l; i++)
861 392 : if (check_ecppcertentry(gel(c,i))==0) return 0;
862 42 : return 1;
863 : }
864 :
865 : GEN
866 49 : primecertexport(GEN c, long flag)
867 : {
868 49 : if (cert_type(c) != c_ECPP) pari_err_IMPL("N-1 certificate");
869 49 : if (!check_ecppcert(c))
870 14 : pari_err_TYPE("primecertexport - invalid certificate", c);
871 35 : return ecppexport(c, flag);
872 : }
873 :
874 : /***********************************************************************/
875 : /** **/
876 : /** PRIME NUMBERS **/
877 : /** **/
878 : /***********************************************************************/
879 :
880 : static struct {
881 : ulong p;
882 : long n;
883 : } prime_table[] = {
884 : { 0, 0},
885 : { 7919, 1000},
886 : { 17389, 2000},
887 : { 27449, 3000},
888 : { 37813, 4000},
889 : { 48611, 5000},
890 : { 59359, 6000},
891 : { 70657, 7000},
892 : { 81799, 8000},
893 : { 93179, 9000},
894 : { 104729, 10000},
895 : { 224737, 20000},
896 : { 350377, 30000},
897 : { 479909, 40000},
898 : { 611953, 50000},
899 : { 746773, 60000},
900 : { 882377, 70000},
901 : { 1020379, 80000},
902 : { 1159523, 90000},
903 : { 1299709, 100000},
904 : { 2750159, 200000},
905 : { 7368787, 500000},
906 : { 15485863, 1000000},
907 : { 32452843, 2000000},
908 : { 86028121, 5000000},
909 : { 179424673, 10000000},
910 : { 373587883, 20000000},
911 : { 982451653, 50000000},
912 : { 2038074743, 100000000},
913 : { 4000000483UL,189961831},
914 : { 4222234741UL,200000000},
915 : #if BITS_IN_LONG == 64
916 : { 5336500537UL, 250000000L},
917 : { 6461335109UL, 300000000L},
918 : { 7594955549UL, 350000000L},
919 : { 8736028057UL, 400000000L},
920 : { 9883692017UL, 450000000L},
921 : { 11037271757UL, 500000000L},
922 : { 13359555403UL, 600000000L},
923 : { 15699342107UL, 700000000L},
924 : { 18054236957UL, 800000000L},
925 : { 20422213579UL, 900000000L},
926 : { 22801763489UL, 1000000000L},
927 : { 47055833459UL, 2000000000L},
928 : { 71856445751UL, 3000000000L},
929 : { 97011687217UL, 4000000000L},
930 : {122430513841UL, 5000000000L},
931 : {148059109201UL, 6000000000L},
932 : {173862636221UL, 7000000000L},
933 : {200000000507UL, 8007105083L},
934 : {225898512559UL, 9000000000L},
935 : {252097800623UL,10000000000L},
936 : {384489816343UL,15000000000L},
937 : {518649879439UL,20000000000L},
938 : {654124187867UL,25000000000L},
939 : {790645490053UL,30000000000L},
940 : {928037044463UL,35000000000L},
941 : {1066173339601UL,40000000000L},
942 : {1344326694119UL,50000000000L},
943 : {1624571841097UL,60000000000L},
944 : {1906555030411UL,70000000000L},
945 : {2190026988349UL,80000000000L},
946 : {2474799787573UL,90000000000L},
947 : {2760727302517UL,100000000000L}
948 : #endif
949 : };
950 : static const int prime_table_len = numberof(prime_table);
951 :
952 : /* find prime closest to n in prime_table. */
953 : static long
954 15893805 : prime_table_closest_p(ulong n)
955 : {
956 : long i;
957 16206113 : for (i = 1; i < prime_table_len; i++)
958 : {
959 16206258 : ulong p = prime_table[i].p;
960 16206258 : if (p > n)
961 : {
962 15893950 : ulong u = n - prime_table[i-1].p;
963 15893950 : if (p - n > u) i--;
964 15893950 : break;
965 : }
966 : }
967 15893805 : if (i == prime_table_len) i = prime_table_len - 1;
968 15893805 : return i;
969 : }
970 :
971 : /* return the n-th successor of prime p > 2; n > 0 */
972 : static GEN
973 77 : prime_successor(ulong p, ulong n)
974 : {
975 77 : GEN Q = utoipos(p), P = NULL;
976 : ulong i;
977 77 : RESET:
978 : for(;;)
979 : {
980 : forprime_t S;
981 77 : GEN L = dbltor(4*n * log(gtodouble(Q) + 1)), R = addii(Q, ceil_safe(L));
982 :
983 77 : forprime_init(&S, addiu(Q, 1), R);
984 77 : Q = R;
985 2402708 : for (i = 1; i <= n; i++)
986 : {
987 2402631 : P = forprime_next(&S);
988 2402631 : if (!P) { n -= i-1; goto RESET; }
989 : }
990 77 : return P;
991 : }
992 : }
993 : /* find the N-th prime */
994 : static GEN
995 273 : prime_table_find_n(ulong N)
996 : {
997 : byteptr d;
998 273 : ulong n, p, maxp = maxprime();
999 : long i;
1000 2240 : for (i = 1; i < prime_table_len; i++)
1001 : {
1002 2240 : n = prime_table[i].n;
1003 2240 : if (n > N)
1004 : {
1005 273 : ulong u = N - prime_table[i-1].n;
1006 273 : if (n - N > u) i--;
1007 273 : break;
1008 : }
1009 : }
1010 273 : if (i == prime_table_len) i = prime_table_len - 1;
1011 273 : p = prime_table[i].p;
1012 273 : n = prime_table[i].n;
1013 273 : if (n > N && p > maxp)
1014 : {
1015 14 : i--;
1016 14 : p = prime_table[i].p;
1017 14 : n = prime_table[i].n;
1018 : }
1019 : /* if beyond prime table, then n <= N */
1020 273 : d = diffptr + n;
1021 273 : if (n > N)
1022 : {
1023 14 : n -= N;
1024 50624 : do { n--; PREC_PRIME_VIADIFF(p,d); } while (n) ;
1025 : }
1026 259 : else if (n < N)
1027 : {
1028 259 : ulong maxpN = maxprimeN();
1029 259 : if (N >= maxpN)
1030 : {
1031 77 : if (N == maxpN) return utoipos(maxp);
1032 77 : if (p < maxp) { p = maxp; n = maxpN; }
1033 77 : return prime_successor(p, N - n);
1034 : }
1035 : /* can find prime(N) in table */
1036 182 : n = N - n;
1037 45234 : do { n--; NEXT_PRIME_VIADIFF(p,d); } while (n) ;
1038 : }
1039 196 : return utoipos(p);
1040 : }
1041 :
1042 : ulong
1043 0 : uprime(long N)
1044 : {
1045 0 : pari_sp av = avma;
1046 : GEN p;
1047 0 : if (N <= 0) pari_err_DOMAIN("prime", "n", "<=",gen_0, stoi(N));
1048 0 : p = prime_table_find_n(N);
1049 0 : if (lgefint(p) != 3) pari_err_OVERFLOW("uprime");
1050 0 : return gc_ulong(av, p[2]);
1051 : }
1052 : GEN
1053 280 : prime(long N)
1054 : {
1055 280 : pari_sp av = avma;
1056 : GEN p;
1057 280 : if (N <= 0) pari_err_DOMAIN("prime", "n", "<=",gen_0, stoi(N));
1058 273 : new_chunk(4); /*HACK*/
1059 273 : p = prime_table_find_n(N);
1060 273 : set_avma(av); return icopy(p);
1061 : }
1062 :
1063 : static void
1064 98 : prime_interval(GEN N, GEN *pa, GEN *pb, GEN *pd)
1065 : {
1066 : GEN a, b, d;
1067 98 : switch(typ(N))
1068 : {
1069 56 : case t_INT:
1070 56 : a = gen_2;
1071 56 : b = subiu(N,1); /* between 2 and N-1 */
1072 56 : d = subiu(N,2);
1073 56 : if (signe(d) <= 0)
1074 7 : pari_err_DOMAIN("randomprime","N", "<=", gen_2, N);
1075 49 : break;
1076 42 : case t_VEC:
1077 42 : if (lg(N) != 3) pari_err_TYPE("randomprime",N);
1078 42 : a = gel(N,1);
1079 42 : b = gel(N,2);
1080 42 : if (gcmp(b, a) < 0)
1081 7 : pari_err_DOMAIN("randomprime","b-a", "<", gen_0, mkvec2(a,b));
1082 35 : if (typ(a) != t_INT)
1083 : {
1084 7 : a = gceil(a);
1085 7 : if (typ(a) != t_INT) pari_err_TYPE("randomprime",a);
1086 : }
1087 35 : if (typ(b) != t_INT)
1088 : {
1089 7 : b = gfloor(b);
1090 7 : if (typ(b) != t_INT) pari_err_TYPE("randomprime",b);
1091 : }
1092 35 : if (cmpiu(a, 2) < 0)
1093 : {
1094 7 : a = gen_2;
1095 7 : d = subiu(b,1);
1096 : }
1097 : else
1098 28 : d = addiu(subii(b,a), 1);
1099 35 : if (signe(d) <= 0)
1100 14 : pari_err_DOMAIN("randomprime","floor(b) - max(ceil(a),2)", "<",
1101 : gen_0, mkvec2(a,b));
1102 21 : break;
1103 0 : default:
1104 0 : pari_err_TYPE("randomprime", N);
1105 : a = b = d = NULL; /*LCOV_EXCL_LINE*/
1106 : }
1107 70 : *pa = a; *pb = b; *pd = d;
1108 70 : }
1109 :
1110 : /* random b-bit prime */
1111 : GEN
1112 56 : randomprime(GEN N)
1113 : {
1114 56 : pari_sp av = avma, av2;
1115 : GEN a, b, d;
1116 56 : if (!N)
1117 : for(;;)
1118 56 : {
1119 63 : ulong p = random_bits(31);
1120 63 : if (uisprime(p)) return utoipos(p);
1121 : }
1122 49 : prime_interval(N, &a, &b, &d); av2 = avma;
1123 : for (;;)
1124 4997 : {
1125 5018 : GEN p = addii(a, randomi(d));
1126 5018 : if (BPSW_psp(p)) return gerepileuptoint(av, p);
1127 4997 : set_avma(av2);
1128 : }
1129 : }
1130 : GEN
1131 105 : randomprime0(GEN N, GEN q)
1132 : {
1133 105 : pari_sp av = avma, av2;
1134 : GEN C, D, a, b, d, r;
1135 105 : if (!q) return randomprime(N);
1136 49 : switch(typ(q))
1137 : {
1138 28 : case t_INT: C = gen_1; D = q; break;
1139 21 : case t_INTMOD: C = gel(q,2); D = gel(q,1); break;
1140 0 : default:
1141 0 : pari_err_TYPE("randomprime", q);
1142 : return NULL;/*LCOV_EXCL_LINE*/
1143 : }
1144 49 : if (!N) N = int2n(31);
1145 49 : prime_interval(N, &a, &b, &d);
1146 49 : r = modii(subii(C, a), D);
1147 49 : if (signe(r)) { a = addii(a, r); d = subii(d, r); }
1148 49 : if (!equali1(gcdii(C,D)))
1149 : {
1150 14 : if (isprime(a)) return gerepilecopy(av, a);
1151 7 : pari_err_COPRIME("randomprime", C, D);
1152 : }
1153 35 : d = divii(d, D); if (!signe(d)) d = gen_1;
1154 35 : av2 = avma;
1155 : for (;;)
1156 3584 : {
1157 3619 : GEN p = addii(a, mulii(D, randomi(d)));
1158 3619 : if (BPSW_psp(p)) return gerepileuptoint(av, p);
1159 3584 : set_avma(av2);
1160 : }
1161 : return NULL;
1162 : }
1163 :
1164 : /* set *pp = nextprime(a) = p
1165 : * *pd so that NEXT_PRIME_VIADIFF(d, p) = nextprime(p+1)
1166 : * *pn so that p = the n-th prime
1167 : * error if nextprime(a) is out of primetable bounds */
1168 : void
1169 15895322 : prime_table_next_p(ulong a, byteptr *pd, ulong *pp, ulong *pn)
1170 : {
1171 : byteptr d;
1172 15895322 : ulong p, n, maxp = maxprime();
1173 15894258 : long i = prime_table_closest_p(a);
1174 15893671 : p = prime_table[i].p;
1175 15893671 : if (p > a && p > maxp)
1176 : {
1177 0 : i--;
1178 0 : p = prime_table[i].p;
1179 : }
1180 : /* if beyond prime table, then p <= a */
1181 15893671 : n = prime_table[i].n;
1182 15893671 : d = diffptr + n;
1183 15893671 : if (p < a)
1184 : {
1185 15808162 : if (a > maxp) pari_err_MAXPRIME(a);
1186 41196250 : do { n++; NEXT_PRIME_VIADIFF(p,d); } while (p < a);
1187 : }
1188 85509 : else if (p != a)
1189 : {
1190 39300624 : do { n--; PREC_PRIME_VIADIFF(p,d); } while (p > a) ;
1191 85509 : if (p < a) { NEXT_PRIME_VIADIFF(p,d); n++; }
1192 : }
1193 15893671 : *pn = n;
1194 15893671 : *pp = p;
1195 15893671 : *pd = d;
1196 15893671 : }
1197 :
1198 : ulong
1199 18101 : uprimepi(ulong a)
1200 : {
1201 18101 : ulong p, n, maxp = maxprime();
1202 18101 : if (a <= maxp)
1203 : {
1204 : byteptr d;
1205 17976 : prime_table_next_p(a, &d, &p, &n);
1206 17976 : return p == a? n: n-1;
1207 : }
1208 : else
1209 : {
1210 125 : long i = prime_table_closest_p(a);
1211 : forprime_t S;
1212 125 : p = prime_table[i].p;
1213 125 : if (p > a)
1214 : {
1215 28 : i--;
1216 28 : p = prime_table[i].p;
1217 : }
1218 : /* p = largest prime in table <= a */
1219 125 : n = prime_table[i].n;
1220 125 : (void)u_forprime_init(&S, p+1, a);
1221 52168646 : for (; p; n++) p = u_forprime_next(&S);
1222 125 : return n-1;
1223 : }
1224 : }
1225 :
1226 : GEN
1227 252 : primepi(GEN x)
1228 : {
1229 252 : pari_sp av = avma;
1230 252 : GEN pp, nn, N = typ(x) == t_INT? x: gfloor(x);
1231 : forprime_t S;
1232 : ulong n, p;
1233 : long i;
1234 252 : if (typ(N) != t_INT) pari_err_TYPE("primepi",N);
1235 252 : if (signe(N) <= 0) return gen_0;
1236 252 : if (lgefint(N) == 3) { n = N[2]; set_avma(av); return utoi(uprimepi(n)); }
1237 1 : i = prime_table_len-1;
1238 1 : p = prime_table[i].p;
1239 1 : n = prime_table[i].n;
1240 1 : (void)forprime_init(&S, utoipos(p+1), N);
1241 1 : nn = setloop(utoipos(n));
1242 1 : pp = gen_0;
1243 3280223 : for (; pp; incloop(nn)) pp = forprime_next(&S);
1244 1 : return gerepileuptoint(av, subiu(nn,1));
1245 : }
1246 :
1247 : /* pi(x) < x/log x * (1 + 1/log x + 2.51/log^2 x)), x>=355991 [ Dusart ]
1248 : * pi(x) < x/(log x - 1.1), x >= 60184 [ Dusart ]
1249 : * ? \p9
1250 : * ? M = 0; for(x = 4, 60184, M = max(M, log(x) - x/primepi(x))); M
1251 : * %1 = 1.11196252 */
1252 : double
1253 399029 : primepi_upper_bound(double x)
1254 : {
1255 399029 : if (x >= 355991)
1256 : {
1257 91 : double L = 1/log(x);
1258 91 : return x * L * (1 + L + 2.51*L*L);
1259 : }
1260 398938 : if (x >= 60184) return x / (log(x) - 1.1);
1261 398938 : if (x < 5) return 2; /* don't bother */
1262 272030 : return x / (log(x) - 1.111963);
1263 : }
1264 : /* pi(x) > x/log x (1 + 1/log x), x >= 599 [ Dusart ]
1265 : * pi(x) > x / (log x + 2), x >= 55 [ Rosser ] */
1266 : double
1267 14 : primepi_lower_bound(double x)
1268 : {
1269 14 : if (x >= 599)
1270 : {
1271 14 : double L = 1/log(x);
1272 14 : return x * L * (1 + L);
1273 : }
1274 0 : if (x < 55) return 0; /* don't bother */
1275 0 : return x / (log(x) + 2.);
1276 : }
1277 : GEN
1278 1 : gprimepi_upper_bound(GEN x)
1279 : {
1280 1 : pari_sp av = avma;
1281 : double L;
1282 1 : if (typ(x) != t_INT) x = gfloor(x);
1283 1 : if (expi(x) <= 1022)
1284 : {
1285 1 : set_avma(av);
1286 1 : return dbltor(primepi_upper_bound(gtodouble(x)));
1287 : }
1288 0 : x = itor(x, LOWDEFAULTPREC);
1289 0 : L = 1 / rtodbl(logr_abs(x));
1290 0 : x = mulrr(x, dbltor(L * (1 + L + 2.51*L*L)));
1291 0 : return gerepileuptoleaf(av, x);
1292 : }
1293 : GEN
1294 1 : gprimepi_lower_bound(GEN x)
1295 : {
1296 1 : pari_sp av = avma;
1297 : double L;
1298 1 : if (typ(x) != t_INT) x = gfloor(x);
1299 1 : if (abscmpiu(x, 55) <= 0) return gen_0;
1300 1 : if (expi(x) <= 1022)
1301 : {
1302 1 : set_avma(av);
1303 1 : return dbltor(primepi_lower_bound(gtodouble(x)));
1304 : }
1305 0 : x = itor(x, LOWDEFAULTPREC);
1306 0 : L = 1 / rtodbl(logr_abs(x));
1307 0 : x = mulrr(x, dbltor(L * (1 + L)));
1308 0 : return gerepileuptoleaf(av, x);
1309 : }
1310 :
1311 : GEN
1312 77 : primes(long n)
1313 : {
1314 : forprime_t S;
1315 : long i;
1316 : GEN y;
1317 77 : if (n <= 0) return cgetg(1, t_VEC);
1318 77 : y = cgetg(n+1, t_VEC);
1319 77 : (void)new_chunk(3*n); /*HACK*/
1320 77 : u_forprime_init(&S, 2, ULONG_MAX);
1321 77 : set_avma((pari_sp)y);
1322 3871 : for (i = 1; i <= n; i++) gel(y, i) = utoipos( u_forprime_next(&S) );
1323 77 : return y;
1324 : }
1325 : GEN
1326 91 : primes_zv(long n)
1327 : {
1328 : forprime_t S;
1329 : long i;
1330 : GEN y;
1331 91 : if (n <= 0) return cgetg(1, t_VECSMALL);
1332 91 : y = cgetg(n+1,t_VECSMALL);
1333 91 : u_forprime_init(&S, 2, ULONG_MAX);
1334 3731 : for (i = 1; i <= n; i++) y[i] = u_forprime_next(&S);
1335 91 : set_avma((pari_sp)y); return y;
1336 : }
1337 : GEN
1338 182 : primes0(GEN N)
1339 : {
1340 182 : switch(typ(N))
1341 : {
1342 77 : case t_INT: return primes(itos(N));
1343 105 : case t_VEC:
1344 105 : if (lg(N) == 3) return primes_interval(gel(N,1),gel(N,2));
1345 : }
1346 0 : pari_err_TYPE("primes", N);
1347 : return NULL;/*LCOV_EXCL_LINE*/
1348 : }
1349 :
1350 : GEN
1351 168 : primes_interval(GEN a, GEN b)
1352 : {
1353 168 : pari_sp av = avma;
1354 : forprime_t S;
1355 : long i, n;
1356 : GEN y, d, p;
1357 168 : if (typ(a) != t_INT)
1358 : {
1359 0 : a = gceil(a);
1360 0 : if (typ(a) != t_INT) pari_err_TYPE("primes_interval",a);
1361 : }
1362 168 : if (typ(b) != t_INT)
1363 : {
1364 7 : b = gfloor(b);
1365 7 : if (typ(b) != t_INT) pari_err_TYPE("primes_interval",b);
1366 : }
1367 161 : if (signe(a) < 0) a = gen_2;
1368 161 : d = subii(b, a);
1369 161 : if (signe(d) < 0 || signe(b) <= 0) { set_avma(av); return cgetg(1, t_VEC); }
1370 161 : if (lgefint(b) == 3)
1371 : {
1372 143 : set_avma(av);
1373 143 : y = primes_interval_zv(itou(a), itou(b));
1374 143 : n = lg(y); settyp(y, t_VEC);
1375 469988 : for (i = 1; i < n; i++) gel(y,i) = utoipos(y[i]);
1376 143 : return y;
1377 : }
1378 : /* at most d+1 primes in [a,b]. If d large, try better bound to lower
1379 : * memory use */
1380 18 : if (abscmpiu(d,100000) > 0)
1381 : {
1382 1 : GEN D = gsub(gprimepi_upper_bound(b), gprimepi_lower_bound(a));
1383 1 : D = ceil_safe(D);
1384 1 : if (cmpii(D, d) < 0) d = D;
1385 : }
1386 18 : n = itos(d)+1;
1387 18 : forprime_init(&S, a, b);
1388 18 : y = cgetg(n+1, t_VEC); i = 1;
1389 5900 : while ((p = forprime_next(&S))) gel(y, i++) = icopy(p);
1390 18 : setlg(y, i); return gerepileupto(av, y);
1391 : }
1392 :
1393 : /* a <= b, at most d primes in [a,b]. Return them */
1394 : static GEN
1395 12554 : primes_interval_i(ulong a, ulong b, ulong d)
1396 : {
1397 12554 : ulong p, i = 1, n = d + 1;
1398 : forprime_t S;
1399 12554 : GEN y = cgetg(n+1, t_VECSMALL);
1400 12554 : pari_sp av = avma;
1401 12554 : u_forprime_init(&S, a, b);
1402 14582261 : while ((p = u_forprime_next(&S))) y[i++] = p;
1403 12554 : set_avma(av); setlg(y, i); stackdummy((pari_sp)(y + i), (pari_sp)(y + n+1));
1404 12554 : return y;
1405 : }
1406 : GEN
1407 12337 : primes_interval_zv(ulong a, ulong b)
1408 : {
1409 : ulong d;
1410 12337 : if (!a) return primes_upto_zv(b);
1411 12337 : if (b < a) return cgetg(1, t_VECSMALL);
1412 12337 : d = b - a;
1413 12337 : if (d > 100000UL)
1414 : {
1415 13 : ulong D = (ulong)ceil(primepi_upper_bound(b)-primepi_lower_bound(a));
1416 13 : if (D < d) d = D;
1417 : }
1418 12337 : return primes_interval_i(a, b, d);
1419 : }
1420 : GEN
1421 217 : primes_upto_zv(ulong b)
1422 : {
1423 : ulong d;
1424 217 : if (b < 2) return cgetg(1, t_VECSMALL);
1425 217 : d = (b > 100000UL)? (ulong)primepi_upper_bound(b): b;
1426 217 : return primes_interval_i(2, b, d);
1427 : }
1428 :
1429 : /***********************************************************************/
1430 : /** **/
1431 : /** PRIVATE PRIME TABLE **/
1432 : /** **/
1433 : /***********************************************************************/
1434 :
1435 : void
1436 228970 : pari_set_primetab(GEN global_primetab)
1437 : {
1438 228970 : if (global_primetab)
1439 : {
1440 227230 : long i, l = lg(global_primetab);
1441 227230 : primetab = cgetg_block(l, t_VEC);
1442 257679 : for (i = 1; i < l; i++)
1443 30673 : gel(primetab,i) = gclone(gel(global_primetab,i));
1444 1740 : } else primetab = cgetg_block(1, t_VEC);
1445 228808 : }
1446 :
1447 : /* delete dummy NULL entries */
1448 : static void
1449 21 : cleanprimetab(GEN T)
1450 : {
1451 21 : long i,j, l = lg(T);
1452 70 : for (i = j = 1; i < l; i++)
1453 49 : if (T[i]) T[j++] = T[i];
1454 21 : setlg(T,j);
1455 21 : }
1456 : /* remove p from T */
1457 : static void
1458 28 : rmprime(GEN T, GEN p)
1459 : {
1460 : long i;
1461 28 : if (typ(p) != t_INT) pari_err_TYPE("removeprimes",p);
1462 28 : i = ZV_search(T, p);
1463 28 : if (!i)
1464 7 : pari_err_DOMAIN("removeprimes","prime","not in",
1465 : strtoGENstr("primetable"), p);
1466 21 : gunclone(gel(T,i)); gel(T,i) = NULL;
1467 21 : cleanprimetab(T);
1468 21 : }
1469 :
1470 : /*stolen from ZV_union_shallow() : clone entries from y */
1471 : static GEN
1472 56 : addp_union(GEN x, GEN y)
1473 : {
1474 56 : long i, j, k, lx = lg(x), ly = lg(y);
1475 56 : GEN z = cgetg(lx + ly - 1, t_VEC);
1476 56 : i = j = k = 1;
1477 84 : while (i<lx && j<ly)
1478 : {
1479 28 : int s = cmpii(gel(x,i), gel(y,j));
1480 28 : if (s < 0)
1481 21 : gel(z,k++) = gel(x,i++);
1482 7 : else if (s > 0)
1483 0 : gel(z,k++) = gclone(gel(y,j++));
1484 : else {
1485 7 : gel(z,k++) = gel(x,i++);
1486 7 : j++;
1487 : }
1488 : }
1489 56 : while (i<lx) gel(z,k++) = gel(x,i++);
1490 147 : while (j<ly) gel(z,k++) = gclone(gel(y,j++));
1491 56 : setlg(z, k); return z;
1492 : }
1493 :
1494 : /* p is NULL, or a single element or a row vector with "primes" to add to
1495 : * prime table. */
1496 : static GEN
1497 189 : addp(GEN *T, GEN p)
1498 : {
1499 189 : pari_sp av = avma;
1500 : long i, l;
1501 : GEN v;
1502 :
1503 189 : if (!p || lg(p) == 1) return *T;
1504 70 : if (!is_vec_t(typ(p))) p = mkvec(p);
1505 :
1506 70 : RgV_check_ZV(p, "addprimes");
1507 63 : v = gen_indexsort_uniq(p, (void*)&cmpii, &cmp_nodata);
1508 63 : p = vecpermute(p, v);
1509 63 : if (abscmpiu(gel(p,1), 2) < 0) pari_err_DOMAIN("addprimes", "p", "<", gen_2,p);
1510 56 : p = addp_union(*T, p);
1511 56 : l = lg(p);
1512 56 : if (l != lg(*T))
1513 : {
1514 56 : GEN old = *T, t = cgetg_block(l, t_VEC);
1515 175 : for (i = 1; i < l; i++) gel(t,i) = gel(p,i);
1516 56 : *T = t; gunclone(old);
1517 : }
1518 56 : set_avma(av); return *T;
1519 : }
1520 : GEN
1521 189 : addprimes(GEN p) { return addp(&primetab, p); }
1522 :
1523 : static GEN
1524 28 : rmprimes(GEN T, GEN prime)
1525 : {
1526 : long i,tx;
1527 :
1528 28 : if (!prime) return T;
1529 28 : tx = typ(prime);
1530 28 : if (is_vec_t(tx))
1531 : {
1532 14 : if (prime == T)
1533 : {
1534 14 : for (i=1; i < lg(prime); i++) gunclone(gel(prime,i));
1535 7 : setlg(prime, 1);
1536 : }
1537 : else
1538 : {
1539 21 : for (i=1; i < lg(prime); i++) rmprime(T, gel(prime,i));
1540 : }
1541 14 : return T;
1542 : }
1543 14 : rmprime(T, prime); return T;
1544 : }
1545 : GEN
1546 28 : removeprimes(GEN prime) { return rmprimes(primetab, prime); }
|
