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