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 387 : init_MR_Jaeschke(MR_Jaeschke_t *S, GEN n)
32 : {
33 387 : S->n = n = absi_shallow(n);
34 387 : S->t = subiu(n,1);
35 387 : S->r1 = vali(S->t);
36 387 : S->t1 = shifti(S->t, -S->r1);
37 387 : S->sqrt1 = cgeti(lg(n)); S->sqrt1[1] = evalsigne(0)|evallgefint(2);
38 387 : S->sqrt2 = cgeti(lg(n)); S->sqrt2[1] = evalsigne(0)|evallgefint(2);
39 387 : }
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 690 : ispsp(MR_Jaeschke_t *S, ulong a)
47 : {
48 690 : pari_sp av = avma;
49 690 : GEN c = Fp_pow(utoipos(a), S->t1, S->n);
50 : long r;
51 :
52 690 : if (is_pm1(c) || equalii(S->t, c)) return 1;
53 : /* go fishing for -1, not for 1 (saves one squaring) */
54 708 : for (r = S->r1 - 1; r; r--) /* r1 - 1 squarings */
55 : {
56 624 : GEN c2 = c;
57 624 : c = remii(sqri(c), S->n);
58 624 : if (equalii(S->t, c))
59 : {
60 271 : if (!signe(S->sqrt1))
61 : {
62 166 : affii(subii(S->n, c2), S->sqrt2);
63 166 : affii(c2, S->sqrt1); return 1;
64 : }
65 : /* saw one earlier: too many sqrt(-1)s mod n ? */
66 105 : return equalii(c2, S->sqrt1) || equalii(c2, S->sqrt2);
67 : }
68 353 : 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 84 : 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 137995 : is2psp(GEN n)
81 : {
82 137995 : GEN c, t = subiu(n, 1);
83 137995 : pari_sp av = avma;
84 137995 : long e = vali(t);
85 :
86 137996 : c = Fp_pow(gen_2, shifti(t, -e), n);
87 137996 : if (is_pm1(c) || equalii(t, c)) return 1;
88 247284 : while (--e)
89 : { /* go fishing for -1, not for 1 (e - 1 squaring) */
90 133829 : c = remii(sqri(c), n);
91 133831 : if (equalii(t, c)) return 1;
92 : /* can return 0 if (c == 1) but very infrequent */
93 122464 : 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 113455 : return 0;
100 : }
101 : static int
102 7999412 : uispsp_pre(ulong a, ulong n, ulong ni)
103 : {
104 7999412 : ulong c, t = n - 1;
105 7999412 : long e = vals(t);
106 :
107 8000435 : c = Fl_powu_pre(a, t >> e, n, ni);
108 7989866 : if (c == 1 || c == t) return 1;
109 14481663 : while (--e)
110 : { /* go fishing for -1, not for 1 (saves one squaring) */
111 7718297 : c = Fl_sqr_pre(c, n, ni);
112 7724927 : if (c == t) return 1;
113 : /* can return 0 if (c == 1) but very infrequent */
114 : }
115 6763366 : return 0;
116 : }
117 : int
118 7813921 : uispsp(ulong a, ulong n)
119 : {
120 : ulong c, t;
121 : long e;
122 :
123 7813921 : if (n & HIGHMASK) return uispsp_pre(a, n, get_Fl_red(n));
124 7538930 : t = n - 1;
125 7538930 : e = vals(t);
126 7550590 : c = Fl_powu(a, t >> e, n);
127 7588793 : if (c == 1 || c == t) return 1;
128 8354241 : while (--e)
129 : { /* go fishing for -1, not for 1 (e - 1 squaring) */
130 5813191 : c = Fl_sqr(c, n);
131 5813866 : if (c == t) return 1;
132 : /* can return 0 if (c == 1) but very infrequent */
133 : }
134 2541050 : 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 3673 : MR_Jaeschke(GEN n)
175 : {
176 : MR_Jaeschke_t S;
177 : pari_sp av;
178 :
179 3673 : if (lgefint(n) == 3) return uisprime(uel(n,2));
180 380 : if (!mod2(n)) return 0;
181 380 : av = avma; init_MR_Jaeschke(&S, n);
182 380 : 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 24541 : LucasMod(GEN n, ulong P, GEN N)
195 : {
196 24541 : pari_sp av = avma;
197 24541 : GEN nd = int_MSW(n);
198 24541 : ulong m = *nd;
199 : long i, j;
200 24541 : GEN v = utoipos(P), v1 = utoipos(P*P - 2);
201 :
202 24541 : if (m == 1)
203 1843 : j = 0;
204 : else
205 : {
206 22698 : j = 1+bfffo(m); /* < BIL */
207 22698 : m <<= j; j = BITS_IN_LONG - j;
208 : }
209 24541 : for (i=lgefint(n)-2;;) /* cf. leftright_pow */
210 : {
211 3288621 : for (; j; m<<=1,j--)
212 : { /* v = v_k, v1 = v_{k+1} */
213 3216507 : if (m&HIGHBIT)
214 : { /* set v = v_{2k+1}, v1 = v_{2k+2} */
215 1139945 : v = subiu(mulii(v,v1), P);
216 1139769 : v1= subiu(sqri(v1), 2);
217 : }
218 : else
219 : {/* set v = v_{2k}, v1 = v_{2k+1} */
220 2076562 : v1= subiu(mulii(v,v1), P);
221 2076396 : v = subiu(sqri(v), 2);
222 : }
223 3216180 : v = modii(v, N);
224 3216161 : v1= modii(v1,N);
225 3216159 : 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 72114 : if (--i == 0) return v;
232 47924 : j = BITS_IN_LONG;
233 47924 : nd=int_precW(nd);
234 47924 : 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 1046883 : u_LucasMod_pre(ulong n, ulong P, ulong N, ulong NI)
242 : {
243 : ulong v, v1, m;
244 : long j;
245 :
246 1046883 : if (n == 1) return P;
247 1046871 : j = 1 + bfffo(n); /* < BIL */
248 1046871 : v = P; v1 = P*P - 2;
249 1046871 : m = n<<j; j = BITS_IN_LONG - j;
250 64257384 : for (; j; m<<=1,j--)
251 : { /* v = v_k, v1 = v_{k+1} */
252 63228329 : if (m & HIGHBIT)
253 : { /* set v = v_{2k+1}, v1 = v_{2k+2} */
254 6399410 : v = Fl_sub(Fl_mul_pre(v,v1,N,NI), P, N);
255 6399298 : v1= Fl_sub(Fl_sqr_pre(v1,N,NI), 2UL, N);
256 : }
257 : else
258 : {/* set v = v_{2k}, v1 = v_{2k+1} */
259 56828919 : v1= Fl_sub(Fl_mul_pre(v,v1,N,NI),P, N);
260 56828443 : v = Fl_sub(Fl_sqr_pre(v,N,NI), 2UL, N);
261 : }
262 : }
263 1029055 : return v;
264 : }
265 :
266 : static ulong
267 1046703 : get_disc(ulong n)
268 : {
269 : ulong b;
270 : long i;
271 1046703 : for (b = 3, i = 0;; b += 2, i++)
272 1272967 : {
273 2319670 : ulong c = b*b - 4; /* = 1 mod 4 */
274 2319670 : if (krouu(n % c, c) < 0) break;
275 1272967 : if (i == 64 && uissquareall(n, &c)) return 0; /* oo loop if N = m^2 */
276 : }
277 1046884 : return b;
278 : }
279 : static int
280 1046710 : uislucaspsp_pre(ulong n, ulong ni)
281 : {
282 : long i, v;
283 1046710 : ulong b, z, m = n + 1;
284 :
285 1046710 : if (!m) return 0; /* neither 2^32-1 nor 2^64-1 are Lucas-pp */
286 1046710 : b = get_disc(n); if (!b) return 0;
287 1046881 : v = vals(m); m >>= v;
288 1046872 : z = u_LucasMod_pre(m, b, n, ni);
289 1047026 : if (z == 2 || z == n-2) return 1;
290 907836 : for (i=1; i<v; i++)
291 : {
292 907831 : if (!z) return 1;
293 478111 : z = Fl_sub(Fl_sqr_pre(z,n,ni), 2UL, n);
294 478112 : if (z == 2) return 0;
295 : }
296 5 : return 0;
297 : }
298 : /* never called; no need to optimize */
299 : int
300 0 : uislucaspsp(ulong n)
301 0 : { return uislucaspsp_pre(n, get_Fl_red(n)); }
302 :
303 : /* N > 3. Caller should check that N is not a square first (taken care of here,
304 : * but inefficient) */
305 : static int
306 24541 : islucaspsp(GEN N)
307 : {
308 24541 : pari_sp av = avma;
309 : GEN m, z;
310 : long i, v;
311 : ulong b;
312 :
313 24541 : for (b=3;; b+=2)
314 27911 : {
315 52452 : ulong c = b*b - 4; /* = 1 mod 4 */
316 52452 : if (b == 129 && Z_issquare(N)) return 0; /* avoid oo loop if N = m^2 */
317 52452 : if (krouu(umodiu(N,c), c) < 0) break;
318 : }
319 24541 : m = addiu(N,1); v = vali(m); m = shifti(m,-v);
320 24541 : z = LucasMod(m, b, N);
321 24541 : if (absequaliu(z, 2)) return 1;
322 19120 : if (equalii(z, subiu(N,2))) return 1;
323 36398 : for (i=1; i<v; i++)
324 : {
325 36171 : if (!signe(z)) return 1;
326 27352 : z = modii(subiu(sqri(z), 2), N);
327 27352 : if (absequaliu(z, 2)) return 0;
328 27352 : if (gc_needed(av,1))
329 : {
330 0 : if(DEBUGMEM>1) pari_warn(warnmem,"islucaspsp");
331 0 : z = gerepileupto(av, z);
332 : }
333 : }
334 227 : return 0;
335 : }
336 :
337 : /* composite strong 2-pseudoprime < 1016801 whose prime divisors are > 101.
338 : * All have a prime divisor <= 661 */
339 : static int
340 3766 : is_2_prp_101(ulong n)
341 : {
342 3766 : switch(n) {
343 0 : case 42799:
344 : case 49141:
345 : case 88357:
346 : case 90751:
347 : case 104653:
348 : case 130561:
349 : case 196093:
350 : case 220729:
351 : case 253241:
352 : case 256999:
353 : case 271951:
354 : case 280601:
355 : case 357761:
356 : case 390937:
357 : case 458989:
358 : case 486737:
359 : case 489997:
360 : case 514447:
361 : case 580337:
362 : case 741751:
363 : case 838861:
364 : case 873181:
365 : case 877099:
366 : case 916327:
367 : case 976873:
368 0 : case 983401: return 1;
369 3766 : } return 0;
370 : }
371 :
372 : static int
373 7683279 : _uispsp(ulong a, long n) { a %= n; return !a || uispsp(a, n); }
374 : static int
375 14759183 : _uisprime(ulong n)
376 : {
377 : #ifdef LONG_IS_64BIT
378 : ulong ni;
379 14505035 : if (n < 341531)
380 5698670 : return _uispsp(9345883071009581737UL, n);
381 8806365 : if (n < 1050535501)
382 1034364 : return _uispsp(336781006125UL, n)
383 1034364 : && _uispsp(9639812373923155UL, n);
384 7772001 : if (n < 350269456337)
385 50254 : return _uispsp(4230279247111683200UL, n)
386 14237 : && _uispsp(14694767155120705706UL, n)
387 64491 : && _uispsp(16641139526367750375UL, n);
388 : /* n & HIGHMASK */
389 7721747 : ni = get_Fl_red(n);
390 7724056 : return (uispsp_pre(2, n, ni) && uislucaspsp_pre(n,ni));
391 : #else
392 254148 : if (n < 360018361) return _uispsp(1143370UL, n) && _uispsp(2350307676UL, n);
393 73441 : return uispsp(15, n) && uispsp(176006322UL, n) && _uispsp(4221622697UL, n);
394 : #endif
395 : }
396 :
397 : int
398 75146155 : uisprime(ulong n)
399 : {
400 75146155 : if (n < 103)
401 7339886 : switch(n)
402 : {
403 5420938 : case 2:
404 : case 3:
405 : case 5:
406 : case 7:
407 : case 11:
408 : case 13:
409 : case 17:
410 : case 19:
411 : case 23:
412 : case 29:
413 : case 31:
414 : case 37:
415 : case 41:
416 : case 43:
417 : case 47:
418 : case 53:
419 : case 59:
420 : case 61:
421 : case 67:
422 : case 71:
423 : case 73:
424 : case 79:
425 : case 83:
426 : case 89:
427 : case 97:
428 5420938 : case 101: return 1;
429 1918948 : default: return 0;
430 : }
431 : /* gcd-extraction is much slower */
432 107386592 : return odd(n) && n % 3 && n % 5 && n % 7 && n % 11 && n % 13 && n % 17
433 20350749 : && n % 19 && n % 23 && n % 29 && n % 31 && n % 37 && n % 41
434 107390267 : && (n < 1849 || _uisprime(n));
435 : }
436 :
437 : /* assume no prime divisor <= 101 */
438 : int
439 16487 : uisprime_101(ulong n)
440 : {
441 16487 : if (n < 1016801) return n < 10609? 1: (uispsp(2, n) && !is_2_prp_101(n));
442 12708 : return _uisprime(n);
443 : }
444 :
445 : /* assume no prime divisor <= 661 */
446 : int
447 177016 : uisprime_661(ulong n)
448 : {
449 177016 : if (n < 1016801) return n < 452929? 1: uispsp(2, n);
450 120026 : return _uisprime(n);
451 : }
452 :
453 : static int
454 373640 : iu_coprime(GEN N, ulong u) { return ugcd(u, umodiu(N,u)) == 1; }
455 : long
456 50329073 : BPSW_psp(GEN N)
457 : {
458 : pari_sp av;
459 50329073 : if (typ(N) != t_INT) pari_err_TYPE("BPSW_psp",N);
460 50329993 : if (signe(N) <= 0) return 0;
461 50268239 : if (lgefint(N) == 3) return uisprime(uel(N,2));
462 196643 : if (!mod2(N)) return 0;
463 : #ifdef LONG_IS_64BIT
464 : /* 16294579238595022365 = 3*5*7*11*13*17*19*23*29*31*37*41*43*47*53
465 : * 7145393598349078859 = 59*61*67*71*73*79*83*89*97*101 */
466 176917 : if (!iu_coprime(N, 16294579238595022365UL) ||
467 109616 : !iu_coprime(N, 7145393598349078859UL)) return 0;
468 : #else
469 : /* 4127218095 = 3*5*7*11*13*17*19*23*37
470 : * 3948078067 = 29*31*41*43*47*53
471 : * 4269855901 = 59*83*89*97*101
472 : * 1673450759 = 61*67*71*73*79 */
473 114222 : if (!iu_coprime(N, 4127218095UL) ||
474 90294 : !iu_coprime(N, 3948078067UL) ||
475 82455 : !iu_coprime(N, 1673450759UL) ||
476 68224 : !iu_coprime(N, 4269855901UL)) return 0;
477 : #endif
478 : /* no prime divisor < 103 */
479 105510 : av = avma;
480 105510 : return gc_long(av, is2psp(N) && islucaspsp(N));
481 : }
482 :
483 : /* can we write n = x^k ? Assume N has no prime divisor <= 2^14, else may
484 : * miss som powers. Not memory clean */
485 : long
486 46910 : isanypower_nosmalldiv(GEN N, GEN *px)
487 : {
488 46910 : GEN x = N, y;
489 46910 : ulong mask = 7;
490 46910 : long ex, k = 1;
491 : forprime_t T;
492 59330 : while (Z_issquareall(x, &y)) { k <<= 1; x = y; }
493 47154 : while ( (ex = is_357_power(x, &y, &mask)) ) { k *= ex; x = y; }
494 46910 : (void)u_forprime_init(&T, 11, ULONG_MAX);
495 : /* stop when x^(1/k) < 2^14 */
496 46910 : while ( (ex = is_pth_power(x, &y, &T, 15)) ) { k *= ex; x = y; }
497 46908 : *px = x; return k;
498 : }
499 :
500 : /* no prime divisor <= 2^14 (> 661) */
501 : long
502 51313 : BPSW_psp_nosmalldiv(GEN N)
503 : {
504 : pari_sp av;
505 51313 : long l = lgefint(N);
506 :
507 51313 : if (l == 3) return uisprime_661(uel(N,2));
508 32498 : av = avma;
509 : /* N large: test for pure power, rarely succeeds, but requires < 1% of
510 : * compositeness test times */
511 32498 : if (bit_accuracy(l) > 512 && isanypower_nosmalldiv(N,&N) != 1)
512 14 : return gc_long(av,0);
513 32482 : N = absi_shallow(N);
514 32482 : return gc_long(av, is2psp(N) && islucaspsp(N));
515 : }
516 :
517 : /***********************************************************************/
518 : /** **/
519 : /** Pocklington-Lehmer **/
520 : /** P-1 primality test **/
521 : /** **/
522 : /***********************************************************************/
523 : /* Assume x BPSW pseudoprime. Check whether it's small enough to be certified
524 : * prime (< 2^64). Reference for strong 2-pseudoprimes:
525 : * http://www.cecm.sfu.ca/Pseudoprimes/index-2-to-64.html */
526 : static int
527 4131113 : BPSW_isprime_small(GEN x)
528 : {
529 4131113 : long l = lgefint(x);
530 : #ifdef LONG_IS_64BIT
531 4065206 : return (l == 3);
532 : #else
533 65907 : return (l <= 4);
534 : #endif
535 : }
536 :
537 : /* Assume N > 1, p^e || N-1, p prime. Find a witness a(p) such that
538 : * a^(N-1) = 1 (mod N)
539 : * a^(N-1)/p - 1 invertible mod N.
540 : * Proves that any divisor of N is 1 mod p^e. Return 0 if N is composite */
541 : static ulong
542 13640 : pl831(GEN N, GEN p)
543 : {
544 13640 : GEN b, c, g, Nmunp = diviiexact(subiu(N,1), p);
545 13640 : pari_sp av = avma;
546 : ulong a;
547 20028 : for(a = 2;; a++, set_avma(av))
548 : {
549 20028 : b = Fp_pow(utoipos(a), Nmunp, N);
550 20028 : if (!equali1(b)) break;
551 : }
552 13640 : c = Fp_pow(b,p,N);
553 13640 : g = gcdii(subiu(b,1), N); /* 0 < g < N */
554 13640 : return (equali1(c) && equali1(g))? a: 0;
555 : }
556 :
557 : /* Brillhart, Lehmer, Selfridge test (Crandall & Pomerance, Th 4.1.5)
558 : * N^(1/3) <= f fully factored, f | N-1. If pl831(p) is true for
559 : * any prime divisor p of f, then any divisor of N is 1 mod f.
560 : * In that case return 1 iff N is prime */
561 : static int
562 49 : BLS_test(GEN N, GEN f)
563 : {
564 : GEN c1, c2, r, q;
565 49 : q = dvmdii(N, f, &r);
566 49 : if (!is_pm1(r)) return 0;
567 49 : c2 = dvmdii(q, f, &c1);
568 : /* N = 1 + f c1 + f^2 c2, 0 <= c_i < f; check whether it is of the form
569 : * (1 + fa)(1 + fb) */
570 49 : return ! Z_issquare(subii(sqri(c1), shifti(c2,2)));
571 : }
572 :
573 : /* BPSW_psp(N) && !BPSW_isprime_small(N). Decide between Pocklington-Lehmer
574 : * and APRCL/ECPP. Return a vector of (small) primes such that PL-witnesses
575 : * guarantee the primality of N. Return NULL if PL is likely too expensive.
576 : * Return gen_0 if BLS test finds N to be composite */
577 : static GEN
578 4429 : BPSW_try_PL(GEN N)
579 : {
580 4429 : ulong B = minuu(1UL<<19, maxprime());
581 4429 : GEN E, p, U, F, N_1 = subiu(N,1);
582 4429 : GEN fa = absZ_factor_limit_strict(N_1, B, &U), P = gel(fa,1);
583 :
584 4429 : if (!U) return P; /* N-1 fully factored */
585 1707 : p = gel(U,1);
586 1707 : E = gel(fa,2);
587 1707 : U = powii(p, gel(U,2)); /* unfactored part of N-1 */
588 1707 : F = (lg(P) == 2)? powii(gel(P,1), gel(E,1)): diviiexact(N_1, U);
589 :
590 : /* N-1 = F U, F factored, U composite, (U,F) = 1 */
591 1707 : if (cmpii(F, U) > 0) return P; /* 1/2-smooth */
592 1700 : if (cmpii(sqri(F), U) > 0) return BLS_test(N,F)? P: gen_0; /* 1/3-smooth */
593 1658 : return NULL; /* not smooth enough */
594 : }
595 :
596 : static GEN isprimePL(GEN N);
597 : /* F is a vector whose entries are primes. For each of them, find a PL
598 : * witness. Return 0 if caller lied and F contains a composite */
599 : static long
600 2771 : PL_certify(GEN N, GEN F)
601 : {
602 2771 : long i, l = lg(F);
603 16334 : for(i = 1; i < l; i++)
604 13563 : if (! pl831(N, gel(F,i))) return 0;
605 2771 : return 1;
606 : }
607 : /* F is a vector whose entries are *believed* to be primes (BPSW_psp).
608 : * For each of them, recording a witness and recursive primality certificate */
609 : static GEN
610 98 : PL_certificate(GEN N, GEN F)
611 : {
612 98 : long i, l = lg(F);
613 : GEN C;
614 98 : if (BPSW_isprime_small(N)) return N;
615 98 : C = cgetg(l, t_VEC);
616 574 : for (i = 1; i < l; i++)
617 : {
618 476 : GEN p = gel(F,i), C0;
619 : ulong w;
620 :
621 476 : if (BPSW_isprime_small(p)) { gel(C,i) = p; continue; }
622 77 : w = pl831(N,p); if (!w) return gen_0;
623 77 : C0 = isprimePL(p);
624 77 : if (isintzero(C0))
625 : { /* composite in prime factorisation ! */
626 0 : err_printf("Not a prime: %Ps", p);
627 0 : pari_err_BUG("PL_certificate [false prime number]");
628 : }
629 77 : gel(C,i) = mkvec3(p,utoipos(w), C0);
630 : }
631 98 : return mkvec2(N, C);
632 : }
633 : /* M a t_MAT */
634 : static int
635 84 : PL_isvalid(GEN v)
636 : {
637 : GEN C, F, F2, N, N1, U;
638 : long i, l;
639 84 : switch(typ(v))
640 : {
641 0 : case t_INT: return BPSW_isprime_small(v) && BPSW_psp(v);
642 84 : case t_VEC: if (lg(v) == 3) break;/*fall through */
643 0 : default: return 0;
644 : }
645 84 : N = gel(v,1);
646 84 : C = gel(v,2);
647 84 : if (typ(N) != t_INT || signe(N) <= 0 || typ(C) != t_VEC) return 0;
648 84 : U = N1 = subiu(N,1);
649 84 : l = lg(C);
650 427 : for (i = 1; i < l; i++)
651 : {
652 350 : GEN p = gel(C,i), a = NULL, C0 = NULL, ap;
653 : long vp;
654 350 : if (typ(p) != t_INT)
655 : {
656 70 : if (lg(p) != 4) return 0;
657 70 : a = gel(p,2); C0 = gel(p,3); p = gel(p,1);
658 70 : if (typ(p) != t_INT || typ(a) != t_INT || !PL_isvalid(C0)) return 0;
659 : }
660 350 : vp = Z_pvalrem(U, p, &U); if (!vp) return 0;
661 343 : if (!a)
662 : {
663 280 : if (!BPSW_isprime_small(p)) return 0;
664 280 : continue;
665 : }
666 63 : if (!equalii(gel(C0,1), p)) return 0;
667 63 : ap = Fp_pow(a, diviiexact(N1,p), N);
668 63 : if (!equali1(gcdii(subiu(ap,1), N)) || !equali1(Fp_pow(ap, p, N))) return 0;
669 : }
670 77 : F = diviiexact(N1, U); /* factored part of N-1 */
671 77 : F2= sqri(F);
672 77 : if (cmpii(F2, N) > 0) return 1;
673 0 : if (cmpii(mulii(F2,F), N) <= 0) return 0;
674 0 : return BLS_test(N,F);
675 : }
676 :
677 : /* Assume N is a strong BPSW pseudoprime, Pocklington-Lehmer primality proof.
678 : * Return gen_0 (nonprime), N (small prime), matrix (large prime)
679 : *
680 : * The matrix has 3 columns, [a,b,c] with
681 : * a[i] prime factor of N-1,
682 : * b[i] witness for a[i] as in pl831
683 : * c[i] check_prime(a[i]) */
684 : static GEN
685 119 : isprimePL(GEN N)
686 : {
687 : GEN cbrtN, N_1, F, f;
688 119 : if (BPSW_isprime_small(N)) return N;
689 98 : cbrtN = sqrtnint(N,3);
690 98 : N_1 = subiu(N,1);
691 98 : F = Z_factor_until(N_1, sqri(cbrtN));
692 98 : f = factorback(F); /* factored part of N-1, f^3 > N */
693 98 : if (DEBUGLEVEL>3)
694 : {
695 0 : GEN r = divri(itor(f,LOWDEFAULTPREC), N);
696 0 : err_printf("Pocklington-Lehmer: proving primality of N = %Ps\n", N);
697 0 : err_printf("Pocklington-Lehmer: N-1 factored up to %Ps! (%.3Ps%%)\n", f, r);
698 : }
699 : /* if N-1 is only N^(1/3)-smooth, BLS test */
700 98 : if (!equalii(f,N_1) && cmpii(sqri(f),N) <= 0 && !BLS_test(N,f))
701 0 : return gen_0; /* Failed, N is composite */
702 98 : F = gel(F,1); settyp(F, t_VEC);
703 98 : return PL_certificate(N, F);
704 : }
705 :
706 : /* assume N a BPSW pseudoprime, in particular, it is odd > 2. Prove N prime */
707 : long
708 4130783 : BPSW_isprime(GEN N)
709 : {
710 : pari_sp av;
711 : long t;
712 : GEN P;
713 4130783 : if (BPSW_isprime_small(N)) return 1;
714 4392 : av = avma; P = BPSW_try_PL(N);
715 4429 : if (!P) /* not smooth enough */
716 1658 : t = expi(N) < 768? isprimeAPRCL(N): isprimeECPP(N);
717 : else
718 2771 : t = (typ(P) == t_INT)? 0: PL_certify(N,P);
719 4429 : return gc_long(av,t);
720 : }
721 :
722 : static long
723 42 : _isprimePL(GEN x)
724 : {
725 42 : pari_sp av = avma;
726 42 : if (!BPSW_psp(x)) return 0;
727 35 : return gc_long(av, !isintzero(isprimePL(x)));
728 : }
729 : GEN
730 39367765 : gisprime(GEN x, long flag)
731 : {
732 39367765 : switch (flag)
733 : {
734 39368557 : case 0: return map_proto_lG(isprime,x);
735 28 : case 1: return map_proto_lG(_isprimePL,x);
736 14 : case 2: return map_proto_lG(isprimeAPRCL,x);
737 21 : case 3: return map_proto_lG(isprimeECPP,x);
738 : }
739 0 : pari_err_FLAG("gisprime");
740 : return NULL;/*LCOV_EXCL_LINE*/
741 : }
742 :
743 : long
744 40070654 : isprime(GEN x) { return BPSW_psp(x) && BPSW_isprime(x); }
745 :
746 : GEN
747 84 : primecert0(GEN x, long flag, long stopat)
748 : {
749 84 : if ((flag || typ(x) == t_INT) && !BPSW_psp(x)) return gen_0;
750 77 : switch(flag)
751 : {
752 70 : case 0: return ecpp0(x, stopat);
753 7 : case 1: { pari_sp av = avma; return gerepilecopy(av, isprimePL(x)); }
754 : }
755 0 : pari_err_FLAG("primecert");
756 : return NULL;/*LCOV_EXCL_LINE*/
757 : }
758 :
759 : GEN
760 0 : primecert(GEN x, long flag)
761 0 : { return primecert0(x, flag, 0); }
762 :
763 : enum { c_VOID = 0, c_ECPP, c_N1 };
764 : static long
765 98 : cert_type(GEN c)
766 : {
767 98 : switch(typ(c))
768 : {
769 0 : case t_INT: return c_ECPP;
770 98 : case t_VEC:
771 98 : if (lg(c) == 3 && typ(gel(c,1)) == t_INT) return c_N1;
772 84 : return c_ECPP;
773 : }
774 0 : return c_VOID;
775 : }
776 :
777 : long
778 56 : primecertisvalid(GEN c)
779 : {
780 56 : switch(typ(c))
781 : {
782 7 : case t_INT: return BPSW_isprime_small(c) && BPSW_psp(c);
783 49 : case t_VEC:
784 49 : return cert_type(c) == c_ECPP? ecppisvalid(c): PL_isvalid(c);
785 : }
786 0 : return 0;
787 : }
788 :
789 : static long
790 392 : check_ecppcertentry(GEN c)
791 : {
792 : GEN v;
793 392 : long i,l = lg(c);
794 392 : if (typ(c)!=t_VEC || l!=6) return 0;
795 1925 : for(i=1; i<=4; i++)
796 1540 : if (typ(gel(c,i))!=t_INT) return 0;
797 385 : v = gel(c,5);
798 385 : if(typ(v)!=t_VEC) return 0;
799 1155 : for(i=1; i<=2; i++)
800 770 : if (typ(gel(v,i))!=t_INT) return 0;
801 385 : return 1;
802 : }
803 :
804 : long
805 56 : check_ecppcert(GEN c)
806 : {
807 : long i, l;
808 56 : switch(typ(c))
809 : {
810 0 : case t_INT: return signe(c) >= 0;
811 56 : case t_VEC: break;
812 0 : default: return 0;
813 : }
814 56 : l = lg(c); if (l == 1) return 0;
815 434 : for (i = 1; i < l; i++)
816 392 : if (check_ecppcertentry(gel(c,i))==0) return 0;
817 42 : return 1;
818 : }
819 :
820 : GEN
821 49 : primecertexport(GEN c, long flag)
822 : {
823 49 : if (cert_type(c) != c_ECPP) pari_err_IMPL("N-1 certificate");
824 49 : if (!check_ecppcert(c))
825 14 : pari_err_TYPE("primecertexport - invalid certificate", c);
826 35 : return ecppexport(c, flag);
827 : }
828 :
829 : /***********************************************************************/
830 : /** **/
831 : /** PRIME NUMBERS **/
832 : /** **/
833 : /***********************************************************************/
834 :
835 : static struct {
836 : ulong p;
837 : long n;
838 : } prime_table[] = {
839 : { 0, 0},
840 : { 7919, 1000},
841 : { 17389, 2000},
842 : { 27449, 3000},
843 : { 37813, 4000},
844 : { 48611, 5000},
845 : { 59359, 6000},
846 : { 70657, 7000},
847 : { 81799, 8000},
848 : { 93179, 9000},
849 : { 104729, 10000},
850 : { 224737, 20000},
851 : { 350377, 30000},
852 : { 479909, 40000},
853 : { 611953, 50000},
854 : { 746773, 60000},
855 : { 882377, 70000},
856 : { 1020379, 80000},
857 : { 1159523, 90000},
858 : { 1299709, 100000},
859 : { 2750159, 200000},
860 : { 7368787, 500000},
861 : { 15485863, 1000000},
862 : { 32452843, 2000000},
863 : { 86028121, 5000000},
864 : { 179424673, 10000000},
865 : { 373587883, 20000000},
866 : { 982451653, 50000000},
867 : { 2038074743, 100000000},
868 : { 4000000483UL,189961831},
869 : { 4222234741UL,200000000},
870 : #if BITS_IN_LONG == 64
871 : { 5336500537UL, 250000000L},
872 : { 6461335109UL, 300000000L},
873 : { 7594955549UL, 350000000L},
874 : { 8736028057UL, 400000000L},
875 : { 9883692017UL, 450000000L},
876 : { 11037271757UL, 500000000L},
877 : { 13359555403UL, 600000000L},
878 : { 15699342107UL, 700000000L},
879 : { 18054236957UL, 800000000L},
880 : { 20422213579UL, 900000000L},
881 : { 22801763489UL, 1000000000L},
882 : { 47055833459UL, 2000000000L},
883 : { 71856445751UL, 3000000000L},
884 : { 97011687217UL, 4000000000L},
885 : {122430513841UL, 5000000000L},
886 : {148059109201UL, 6000000000L},
887 : {173862636221UL, 7000000000L},
888 : {200000000507UL, 8007105083L},
889 : {225898512559UL, 9000000000L},
890 : {252097800623UL,10000000000L},
891 : {384489816343UL,15000000000L},
892 : {518649879439UL,20000000000L},
893 : {654124187867UL,25000000000L},
894 : {790645490053UL,30000000000L},
895 : {928037044463UL,35000000000L},
896 : {1066173339601UL,40000000000L},
897 : {1344326694119UL,50000000000L},
898 : {1624571841097UL,60000000000L},
899 : {1906555030411UL,70000000000L},
900 : {2190026988349UL,80000000000L},
901 : {2474799787573UL,90000000000L},
902 : {2760727302517UL,100000000000L}
903 : #endif
904 : };
905 : static const int prime_table_len = numberof(prime_table);
906 :
907 : /* find prime closest to n in prime_table. */
908 : static long
909 17835222 : prime_table_closest_p(ulong n)
910 : {
911 : long i;
912 20357954 : for (i = 1; i < prime_table_len; i++)
913 : {
914 20358149 : ulong p = prime_table[i].p;
915 20358149 : if (p > n)
916 : {
917 17835417 : ulong u = n - prime_table[i-1].p;
918 17835417 : if (p - n > u) i--;
919 17835417 : break;
920 : }
921 : }
922 17835222 : if (i == prime_table_len) i = prime_table_len - 1;
923 17835222 : return i;
924 : }
925 :
926 : /* return the n-th successor of prime p > 2; n > 0 */
927 : static GEN
928 56 : prime_successor(ulong p, ulong n)
929 : {
930 56 : GEN Q = utoipos(p), P = NULL;
931 : ulong i;
932 56 : RESET:
933 : for(;;)
934 : {
935 : forprime_t S;
936 56 : GEN L = dbltor(4*n * log(gtodouble(Q) + 1)), R = addii(Q, ceil_safe(L));
937 :
938 56 : forprime_init(&S, addiu(Q, 1), R);
939 56 : Q = R;
940 2325176 : for (i = 1; i <= n; i++)
941 : {
942 2325120 : P = forprime_next(&S);
943 2325120 : if (!P) { n -= i-1; goto RESET; }
944 : }
945 56 : return P;
946 : }
947 : }
948 : /* find the N-th prime */
949 : static GEN
950 273 : prime_table_find_n(ulong N)
951 : {
952 : byteptr d;
953 273 : ulong n, p, maxp = maxprime();
954 : long i;
955 2240 : for (i = 1; i < prime_table_len; i++)
956 : {
957 2240 : n = prime_table[i].n;
958 2240 : if (n > N)
959 : {
960 273 : ulong u = N - prime_table[i-1].n;
961 273 : if (n - N > u) i--;
962 273 : break;
963 : }
964 : }
965 273 : if (i == prime_table_len) i = prime_table_len - 1;
966 273 : p = prime_table[i].p;
967 273 : n = prime_table[i].n;
968 273 : if (n > N && p > maxp)
969 : {
970 0 : i--;
971 0 : p = prime_table[i].p;
972 0 : n = prime_table[i].n;
973 : }
974 : /* if beyond prime table, then n <= N */
975 273 : d = diffptr + n;
976 273 : if (n > N)
977 : {
978 28 : n -= N;
979 113120 : do { n--; PREC_PRIME_VIADIFF(p,d); } while (n) ;
980 : }
981 245 : else if (n < N)
982 : {
983 245 : ulong maxpN = maxprimeN();
984 245 : if (N >= maxpN)
985 : {
986 56 : if (N == maxpN) return utoipos(maxp);
987 56 : if (p < maxp) { p = maxp; n = maxpN; }
988 56 : return prime_successor(p, N - n);
989 : }
990 : /* can find prime(N) in table */
991 189 : n = N - n;
992 56007 : do { n--; NEXT_PRIME_VIADIFF(p,d); } while (n) ;
993 : }
994 217 : return utoipos(p);
995 : }
996 :
997 : ulong
998 0 : uprime(long N)
999 : {
1000 0 : pari_sp av = avma;
1001 : GEN p;
1002 0 : if (N <= 0) pari_err_DOMAIN("prime", "n", "<=",gen_0, stoi(N));
1003 0 : p = prime_table_find_n(N);
1004 0 : if (lgefint(p) != 3) pari_err_OVERFLOW("uprime");
1005 0 : return gc_ulong(av, p[2]);
1006 : }
1007 : GEN
1008 280 : prime(long N)
1009 : {
1010 280 : pari_sp av = avma;
1011 : GEN p;
1012 280 : if (N <= 0) pari_err_DOMAIN("prime", "n", "<=",gen_0, stoi(N));
1013 273 : new_chunk(4); /*HACK*/
1014 273 : p = prime_table_find_n(N);
1015 273 : set_avma(av); return icopy(p);
1016 : }
1017 :
1018 : static void
1019 133 : prime_interval(GEN N, GEN *pa, GEN *pb, GEN *pd)
1020 : {
1021 : GEN a, b, d;
1022 133 : switch(typ(N))
1023 : {
1024 77 : case t_INT:
1025 77 : a = gen_2;
1026 77 : b = subiu(N,1); /* between 2 and N-1 */
1027 77 : d = subiu(N,2);
1028 77 : if (signe(d) <= 0)
1029 7 : pari_err_DOMAIN("randomprime","N", "<=", gen_2, N);
1030 70 : break;
1031 56 : case t_VEC:
1032 56 : if (lg(N) != 3) pari_err_TYPE("randomprime",N);
1033 56 : a = gel(N,1);
1034 56 : b = gel(N,2);
1035 56 : if (gcmp(b, a) < 0)
1036 7 : pari_err_DOMAIN("randomprime","b-a", "<", gen_0, mkvec2(a,b));
1037 49 : if (typ(a) != t_INT)
1038 : {
1039 7 : a = gceil(a);
1040 7 : if (typ(a) != t_INT) pari_err_TYPE("randomprime",a);
1041 : }
1042 49 : if (typ(b) != t_INT)
1043 : {
1044 7 : b = gfloor(b);
1045 7 : if (typ(b) != t_INT) pari_err_TYPE("randomprime",b);
1046 : }
1047 49 : if (cmpiu(a, 2) < 0)
1048 : {
1049 7 : a = gen_2;
1050 7 : d = subiu(b,1);
1051 : }
1052 : else
1053 42 : d = addiu(subii(b,a), 1);
1054 49 : if (signe(d) <= 0)
1055 14 : pari_err_DOMAIN("randomprime","floor(b) - max(ceil(a),2)", "<",
1056 : gen_0, mkvec2(a,b));
1057 35 : break;
1058 0 : default:
1059 0 : pari_err_TYPE("randomprime", N);
1060 : a = b = d = NULL; /*LCOV_EXCL_LINE*/
1061 : }
1062 105 : *pa = a; *pb = b; *pd = d;
1063 105 : }
1064 :
1065 : /* random b-bit prime */
1066 : GEN
1067 91 : randomprime(GEN N)
1068 : {
1069 91 : pari_sp av = avma, av2;
1070 : GEN a, b, d;
1071 91 : if (!N)
1072 : for(;;)
1073 56 : {
1074 63 : ulong p = random_bits(31);
1075 63 : if (uisprime(p)) return utoipos(p);
1076 : }
1077 84 : prime_interval(N, &a, &b, &d); av2 = avma;
1078 : for (;;)
1079 1685 : {
1080 1741 : GEN p = addii(a, randomi(d));
1081 1741 : if (BPSW_psp(p)) return gerepileuptoint(av, p);
1082 1685 : set_avma(av2);
1083 : }
1084 : }
1085 : GEN
1086 140 : randomprime0(GEN N, GEN q)
1087 : {
1088 140 : pari_sp av = avma, av2;
1089 : GEN C, D, a, b, d, r;
1090 140 : if (!q) return randomprime(N);
1091 49 : switch(typ(q))
1092 : {
1093 28 : case t_INT: C = gen_1; D = q; break;
1094 21 : case t_INTMOD: C = gel(q,2); D = gel(q,1); break;
1095 0 : default:
1096 0 : pari_err_TYPE("randomprime", q);
1097 : return NULL;/*LCOV_EXCL_LINE*/
1098 : }
1099 49 : if (!N) N = int2n(31);
1100 49 : prime_interval(N, &a, &b, &d);
1101 49 : r = modii(subii(C, a), D);
1102 49 : if (signe(r)) { a = addii(a, r); d = subii(d, r); }
1103 49 : if (!equali1(gcdii(C,D)))
1104 : {
1105 14 : if (isprime(a)) return gerepilecopy(av, a);
1106 7 : pari_err_COPRIME("randomprime", C, D);
1107 : }
1108 35 : d = divii(d, D); if (!signe(d)) d = gen_1;
1109 35 : av2 = avma;
1110 : for (;;)
1111 3584 : {
1112 3619 : GEN p = addii(a, mulii(D, randomi(d)));
1113 3619 : if (BPSW_psp(p)) return gerepileuptoint(av, p);
1114 3584 : set_avma(av2);
1115 : }
1116 : return NULL;
1117 : }
1118 :
1119 : /* set *pp = nextprime(a) = p
1120 : * *pd so that NEXT_PRIME_VIADIFF(d, p) = nextprime(p+1)
1121 : * *pn so that p = the n-th prime
1122 : * error if nextprime(a) is out of primetable bounds */
1123 : void
1124 17837292 : prime_table_next_p(ulong a, byteptr *pd, ulong *pp, ulong *pn)
1125 : {
1126 : byteptr d;
1127 17837292 : ulong p, n, maxp = maxprime();
1128 17835589 : long i = prime_table_closest_p(a);
1129 17835303 : p = prime_table[i].p;
1130 17835303 : if (p > a && p > maxp)
1131 : {
1132 0 : i--;
1133 0 : p = prime_table[i].p;
1134 : }
1135 : /* if beyond prime table, then p <= a */
1136 17835303 : n = prime_table[i].n;
1137 17835303 : d = diffptr + n;
1138 17835303 : if (p < a)
1139 : {
1140 17126933 : if (a > maxp) pari_err_MAXPRIME(a);
1141 46372504 : do { n++; NEXT_PRIME_VIADIFF(p,d); } while (p < a);
1142 : }
1143 708370 : else if (p != a)
1144 : {
1145 338653212 : do { n--; PREC_PRIME_VIADIFF(p,d); } while (p > a) ;
1146 708370 : if (p < a) { NEXT_PRIME_VIADIFF(p,d); n++; }
1147 : }
1148 17835303 : *pn = n;
1149 17835303 : *pp = p;
1150 17835303 : *pd = d;
1151 17835303 : }
1152 :
1153 : ulong
1154 18423 : uprimepi(ulong a)
1155 : {
1156 18423 : ulong p, n, maxp = maxprime();
1157 18423 : if (a <= maxp)
1158 : {
1159 : byteptr d;
1160 18312 : prime_table_next_p(a, &d, &p, &n);
1161 18312 : return p == a? n: n-1;
1162 : }
1163 : else
1164 : {
1165 111 : long i = prime_table_closest_p(a);
1166 : forprime_t S;
1167 111 : p = prime_table[i].p;
1168 111 : if (p > a)
1169 : {
1170 28 : i--;
1171 28 : p = prime_table[i].p;
1172 : }
1173 : /* p = largest prime in table <= a */
1174 111 : n = prime_table[i].n;
1175 111 : (void)u_forprime_init(&S, p+1, a);
1176 52121158 : for (; p; n++) p = u_forprime_next(&S);
1177 111 : return n-1;
1178 : }
1179 : }
1180 :
1181 : GEN
1182 252 : primepi(GEN x)
1183 : {
1184 252 : pari_sp av = avma;
1185 252 : GEN pp, nn, N = typ(x) == t_INT? x: gfloor(x);
1186 : forprime_t S;
1187 : ulong n, p;
1188 : long i;
1189 252 : if (typ(N) != t_INT) pari_err_TYPE("primepi",N);
1190 252 : if (signe(N) <= 0) return gen_0;
1191 252 : if (lgefint(N) == 3) return gc_utoi(av, uprimepi(N[2]));
1192 1 : i = prime_table_len-1;
1193 1 : p = prime_table[i].p;
1194 1 : n = prime_table[i].n;
1195 1 : (void)forprime_init(&S, utoipos(p+1), N);
1196 1 : nn = setloop(utoipos(n));
1197 1 : pp = gen_0;
1198 3280223 : for (; pp; incloop(nn)) pp = forprime_next(&S);
1199 1 : return gerepileuptoint(av, subiu(nn,1));
1200 : }
1201 :
1202 : /* pi(x) < x/log x * (1 + 1/log x + 2.51/log^2 x)), x>=355991 [ Dusart ]
1203 : * pi(x) < x/(log x - 1.1), x >= 60184 [ Dusart ]
1204 : * ? \p9
1205 : * ? M = 0; for(x = 4, 60184, M = max(M, log(x) - x/primepi(x))); M
1206 : * %1 = 1.11196252 */
1207 : double
1208 402582 : primepi_upper_bound(double x)
1209 : {
1210 402582 : if (x >= 355991)
1211 : {
1212 1904 : double L = 1/log(x);
1213 1904 : return x * L * (1 + L + 2.51*L*L);
1214 : }
1215 400678 : if (x >= 60184) return x / (log(x) - 1.1);
1216 400678 : if (x < 5) return 2; /* don't bother */
1217 273462 : return x / (log(x) - 1.111963);
1218 : }
1219 : /* pi(x) > x/log x (1 + 1/log x), x >= 599 [ Dusart ]
1220 : * pi(x) > x / (log x + 2), x >= 55 [ Rosser ] */
1221 : double
1222 1827 : primepi_lower_bound(double x)
1223 : {
1224 1827 : if (x >= 599)
1225 : {
1226 14 : double L = 1/log(x);
1227 14 : return x * L * (1 + L);
1228 : }
1229 1813 : if (x < 55) return 0; /* don't bother */
1230 0 : return x / (log(x) + 2.);
1231 : }
1232 : GEN
1233 1 : gprimepi_upper_bound(GEN x)
1234 : {
1235 1 : pari_sp av = avma;
1236 : double L;
1237 1 : if (typ(x) != t_INT) x = gfloor(x);
1238 1 : if (expi(x) <= 1022)
1239 : {
1240 1 : set_avma(av);
1241 1 : return dbltor(primepi_upper_bound(gtodouble(x)));
1242 : }
1243 0 : x = itor(x, LOWDEFAULTPREC);
1244 0 : L = 1 / rtodbl(logr_abs(x));
1245 0 : x = mulrr(x, dbltor(L * (1 + L + 2.51*L*L)));
1246 0 : return gerepileuptoleaf(av, x);
1247 : }
1248 : GEN
1249 1 : gprimepi_lower_bound(GEN x)
1250 : {
1251 1 : pari_sp av = avma;
1252 : double L;
1253 1 : if (typ(x) != t_INT) x = gfloor(x);
1254 1 : if (abscmpiu(x, 55) <= 0) return gen_0;
1255 1 : if (expi(x) <= 1022)
1256 : {
1257 1 : set_avma(av);
1258 1 : return dbltor(primepi_lower_bound(gtodouble(x)));
1259 : }
1260 0 : x = itor(x, LOWDEFAULTPREC);
1261 0 : L = 1 / rtodbl(logr_abs(x));
1262 0 : x = mulrr(x, dbltor(L * (1 + L)));
1263 0 : return gerepileuptoleaf(av, x);
1264 : }
1265 :
1266 : GEN
1267 98 : primes(long n)
1268 : {
1269 : forprime_t S;
1270 : long i;
1271 : GEN y;
1272 98 : if (n <= 0) return cgetg(1, t_VEC);
1273 98 : y = cgetg(n+1, t_VEC);
1274 98 : (void)new_chunk(3*n); /*HACK*/
1275 98 : u_forprime_init(&S, 2, (ulong)n > maxprimeN()? ULONG_MAX: maxprime());
1276 98 : set_avma((pari_sp)y);
1277 5992 : for (i = 1; i <= n; i++) gel(y, i) = utoipos( u_forprime_next(&S) );
1278 98 : return y;
1279 : }
1280 : GEN
1281 91 : primes_zv(long n)
1282 : {
1283 : forprime_t S;
1284 : long i;
1285 : GEN y;
1286 91 : if (n <= 0) return cgetg(1, t_VECSMALL);
1287 91 : y = cgetg(n+1,t_VECSMALL);
1288 : /* don't initialize sieve if all fits in primetable */
1289 91 : u_forprime_init(&S, 2, (ulong)n > maxprimeN()? ULONG_MAX: maxprime());
1290 3731 : for (i = 1; i <= n; i++) y[i] = u_forprime_next(&S);
1291 91 : set_avma((pari_sp)y); return y;
1292 : }
1293 : GEN
1294 224 : primes0(GEN N)
1295 : {
1296 224 : switch(typ(N))
1297 : {
1298 98 : case t_INT: return primes(itos(N));
1299 126 : case t_VEC:
1300 126 : if (lg(N) == 3) return primes_interval(gel(N,1),gel(N,2));
1301 : }
1302 0 : pari_err_TYPE("primes", N);
1303 : return NULL;/*LCOV_EXCL_LINE*/
1304 : }
1305 :
1306 : GEN
1307 189 : primes_interval(GEN a, GEN b)
1308 : {
1309 189 : pari_sp av = avma;
1310 : forprime_t S;
1311 : long i, n;
1312 : GEN y, d, p;
1313 189 : if (typ(a) != t_INT)
1314 : {
1315 7 : a = gceil(a);
1316 7 : if (typ(a) != t_INT) pari_err_TYPE("primes_interval",a);
1317 : }
1318 189 : if (typ(b) != t_INT)
1319 : {
1320 14 : b = gfloor(b);
1321 14 : if (typ(b) != t_INT) pari_err_TYPE("primes_interval",b);
1322 : }
1323 182 : if (signe(a) < 0) a = gen_2;
1324 182 : d = subii(b, a);
1325 182 : if (signe(d) < 0 || signe(b) <= 0) { set_avma(av); return cgetg(1, t_VEC); }
1326 182 : if (lgefint(b) == 3)
1327 : {
1328 156 : set_avma(av);
1329 156 : y = primes_interval_zv(itou(a), itou(b));
1330 156 : n = lg(y); settyp(y, t_VEC);
1331 470082 : for (i = 1; i < n; i++) gel(y,i) = utoipos(y[i]);
1332 156 : return y;
1333 : }
1334 : /* at most d+1 primes in [a,b]. If d large, try better bound to lower
1335 : * memory use */
1336 26 : if (abscmpiu(d,100000) > 0)
1337 : {
1338 1 : GEN D = gsub(gprimepi_upper_bound(b), gprimepi_lower_bound(a));
1339 1 : D = ceil_safe(D);
1340 1 : if (cmpii(D, d) < 0) d = D;
1341 : }
1342 26 : n = itos(d)+1;
1343 26 : forprime_init(&S, a, b);
1344 26 : y = cgetg(n+1, t_VEC); i = 1;
1345 5995 : while ((p = forprime_next(&S))) gel(y, i++) = icopy(p);
1346 26 : setlg(y, i); return gerepileupto(av, y);
1347 : }
1348 :
1349 : /* a <= b, at most d primes in [a,b]. Return them */
1350 : static GEN
1351 14516 : primes_interval_i(ulong a, ulong b, ulong d)
1352 : {
1353 14516 : ulong p, i = 1, n = d + 1;
1354 : forprime_t S;
1355 14516 : GEN y = cgetg(n+1, t_VECSMALL);
1356 14516 : pari_sp av = avma;
1357 14516 : u_forprime_init(&S, a, b);
1358 163296643 : while ((p = u_forprime_next(&S))) y[i++] = p;
1359 14516 : set_avma(av); setlg(y, i); stackdummy((pari_sp)(y + i), (pari_sp)(y + n+1));
1360 14516 : return y;
1361 : }
1362 : GEN
1363 14327 : primes_interval_zv(ulong a, ulong b)
1364 : {
1365 : ulong d;
1366 14327 : if (a < 3) return primes_upto_zv(b);
1367 14222 : if (b < a) return cgetg(1, t_VECSMALL);
1368 14222 : d = b - a;
1369 14222 : if (d > 100000UL)
1370 : {
1371 1826 : ulong D = (ulong)ceil(primepi_upper_bound(b)-primepi_lower_bound(a));
1372 1826 : if (D < d) d = D;
1373 : }
1374 14222 : return primes_interval_i(a, b, d);
1375 : }
1376 : GEN
1377 294 : primes_upto_zv(ulong b)
1378 : {
1379 : ulong d;
1380 294 : if (b < 2) return cgetg(1, t_VECSMALL);
1381 294 : d = (b > 100000UL)? (ulong)primepi_upper_bound(b): b;
1382 294 : return primes_interval_i(2, b, d);
1383 : }
1384 :
1385 : /***********************************************************************/
1386 : /** **/
1387 : /** PRIVATE PRIME TABLE **/
1388 : /** **/
1389 : /***********************************************************************/
1390 :
1391 : void
1392 349693 : pari_set_primetab(GEN global_primetab)
1393 : {
1394 349693 : if (global_primetab)
1395 : {
1396 347879 : long i, l = lg(global_primetab);
1397 347879 : primetab = cgetg_block(l, t_VEC);
1398 386213 : for (i = 1; i < l; i++)
1399 38502 : gel(primetab,i) = gclone(gel(global_primetab,i));
1400 1814 : } else primetab = cgetg_block(1, t_VEC);
1401 349534 : }
1402 :
1403 : /* delete dummy NULL entries */
1404 : static void
1405 21 : cleanprimetab(GEN T)
1406 : {
1407 21 : long i,j, l = lg(T);
1408 70 : for (i = j = 1; i < l; i++)
1409 49 : if (T[i]) T[j++] = T[i];
1410 21 : setlg(T,j);
1411 21 : }
1412 : /* remove p from T */
1413 : static void
1414 28 : rmprime(GEN T, GEN p)
1415 : {
1416 : long i;
1417 28 : if (typ(p) != t_INT) pari_err_TYPE("removeprimes",p);
1418 28 : i = ZV_search(T, p);
1419 28 : if (!i)
1420 7 : pari_err_DOMAIN("removeprimes","prime","not in",
1421 : strtoGENstr("primetable"), p);
1422 21 : gunclone(gel(T,i)); gel(T,i) = NULL;
1423 21 : cleanprimetab(T);
1424 21 : }
1425 :
1426 : /*stolen from ZV_union_shallow() : clone entries from y */
1427 : static GEN
1428 56 : addp_union(GEN x, GEN y)
1429 : {
1430 56 : long i, j, k, lx = lg(x), ly = lg(y);
1431 56 : GEN z = cgetg(lx + ly - 1, t_VEC);
1432 56 : i = j = k = 1;
1433 84 : while (i<lx && j<ly)
1434 : {
1435 28 : int s = cmpii(gel(x,i), gel(y,j));
1436 28 : if (s < 0)
1437 21 : gel(z,k++) = gel(x,i++);
1438 7 : else if (s > 0)
1439 0 : gel(z,k++) = gclone(gel(y,j++));
1440 : else {
1441 7 : gel(z,k++) = gel(x,i++);
1442 7 : j++;
1443 : }
1444 : }
1445 56 : while (i<lx) gel(z,k++) = gel(x,i++);
1446 147 : while (j<ly) gel(z,k++) = gclone(gel(y,j++));
1447 56 : setlg(z, k); return z;
1448 : }
1449 :
1450 : /* p is NULL, or a single element or a row vector with "primes" to add to
1451 : * prime table. */
1452 : static GEN
1453 189 : addp(GEN *T, GEN p)
1454 : {
1455 189 : pari_sp av = avma;
1456 : long i, l;
1457 : GEN v;
1458 :
1459 189 : if (!p || lg(p) == 1) return *T;
1460 70 : if (!is_vec_t(typ(p))) p = mkvec(p);
1461 :
1462 70 : RgV_check_ZV(p, "addprimes");
1463 63 : v = gen_indexsort_uniq(p, (void*)&cmpii, &cmp_nodata);
1464 63 : p = vecpermute(p, v);
1465 63 : if (abscmpiu(gel(p,1), 2) < 0) pari_err_DOMAIN("addprimes", "p", "<", gen_2,p);
1466 56 : p = addp_union(*T, p);
1467 56 : l = lg(p);
1468 56 : if (l != lg(*T))
1469 : {
1470 56 : GEN old = *T, t = cgetg_block(l, t_VEC);
1471 175 : for (i = 1; i < l; i++) gel(t,i) = gel(p,i);
1472 56 : *T = t; gunclone(old);
1473 : }
1474 56 : set_avma(av); return *T;
1475 : }
1476 : GEN
1477 189 : addprimes(GEN p) { return addp(&primetab, p); }
1478 :
1479 : static GEN
1480 28 : rmprimes(GEN T, GEN prime)
1481 : {
1482 : long i,tx;
1483 :
1484 28 : if (!prime) return T;
1485 28 : tx = typ(prime);
1486 28 : if (is_vec_t(tx))
1487 : {
1488 14 : if (prime == T)
1489 : {
1490 14 : for (i=1; i < lg(prime); i++) gunclone(gel(prime,i));
1491 7 : setlg(prime, 1);
1492 : }
1493 : else
1494 : {
1495 21 : for (i=1; i < lg(prime); i++) rmprime(T, gel(prime,i));
1496 : }
1497 14 : return T;
1498 : }
1499 14 : rmprime(T, prime); return T;
1500 : }
1501 : GEN
1502 28 : removeprimes(GEN prime) { return rmprimes(primetab, prime); }
|
