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 523 : init_MR_Jaeschke(MR_Jaeschke_t *S, GEN n)
32 : {
33 523 : S->n = n = absi_shallow(n);
34 523 : S->t = subiu(n,1);
35 523 : S->r1 = vali(S->t);
36 523 : S->t1 = shifti(S->t, -S->r1);
37 523 : S->sqrt1 = cgeti(lg(n)); S->sqrt1[1] = evalsigne(0)|evallgefint(2);
38 523 : S->sqrt2 = cgeti(lg(n)); S->sqrt2[1] = evalsigne(0)|evallgefint(2);
39 523 : }
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 917 : ispsp(MR_Jaeschke_t *S, ulong a)
47 : {
48 917 : pari_sp av = avma;
49 917 : GEN c = Fp_pow(utoipos(a), S->t1, S->n);
50 : long r;
51 :
52 917 : if (is_pm1(c) || equalii(S->t, c)) return 1;
53 : /* go fishing for -1, not for 1 (saves one squaring) */
54 1009 : for (r = S->r1 - 1; r; r--) /* r1 - 1 squarings */
55 : {
56 880 : GEN c2 = c;
57 880 : c = remii(sqri(c), S->n);
58 880 : if (equalii(S->t, c))
59 : {
60 344 : if (!signe(S->sqrt1))
61 : {
62 213 : affii(subii(S->n, c2), S->sqrt2);
63 213 : affii(c2, S->sqrt1); return 1;
64 : }
65 : /* saw one earlier: too many sqrt(-1)s mod n ? */
66 131 : return equalii(c2, S->sqrt1) || equalii(c2, S->sqrt2);
67 : }
68 536 : 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 129 : 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 137828 : is2psp(GEN n)
81 : {
82 137828 : GEN c, t = subiu(n, 1);
83 137821 : pari_sp av = avma;
84 137821 : long e = vali(t);
85 :
86 137821 : c = Fp_pow(gen_2, shifti(t, -e), n);
87 137838 : if (is_pm1(c) || equalii(t, c)) return 1;
88 247034 : while (--e)
89 : { /* go fishing for -1, not for 1 (e - 1 squaring) */
90 133728 : c = remii(sqri(c), n);
91 133728 : if (equalii(t, c)) return 1;
92 : /* can return 0 if (c == 1) but very infrequent */
93 122358 : 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 113306 : return 0;
100 : }
101 : static int
102 7836593 : uispsp_pre(ulong a, ulong n, ulong ni)
103 : {
104 7836593 : ulong c, t = n - 1;
105 7836593 : long e = vals(t);
106 :
107 7837226 : c = Fl_powu_pre(a, t >> e, n, ni);
108 7826005 : if (c == 1 || c == t) return 1;
109 14253774 : while (--e)
110 : { /* go fishing for -1, not for 1 (saves one squaring) */
111 7589712 : c = Fl_sqr_pre(c, n, ni);
112 7596676 : if (c == t) return 1;
113 : /* can return 0 if (c == 1) but very infrequent */
114 : }
115 6664062 : return 0;
116 : }
117 : int
118 1030349 : uispsp(ulong a, ulong n)
119 : {
120 : ulong c, t;
121 : long e;
122 :
123 1030349 : if (n & HIGHMASK) return uispsp_pre(a, n, get_Fl_red(n));
124 820425 : t = n - 1;
125 820425 : e = vals(t);
126 820218 : c = Fl_powu(a, t >> e, n);
127 824250 : if (c == 1 || c == t) return 1;
128 923317 : while (--e)
129 : { /* go fishing for -1, not for 1 (e - 1 squaring) */
130 617322 : c = Fl_sqr(c, n);
131 617215 : if (c == t) return 1;
132 : /* can return 0 if (c == 1) but very infrequent */
133 : }
134 305995 : 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 5669 : MR_Jaeschke(GEN n)
175 : {
176 : MR_Jaeschke_t S;
177 : pari_sp av;
178 :
179 5669 : if (lgefint(n) == 3) return uisprime(uel(n,2));
180 516 : if (!mod2(n)) return 0;
181 516 : av = avma; init_MR_Jaeschke(&S, n);
182 516 : 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 24532 : LucasMod(GEN n, ulong P, GEN N)
195 : {
196 24532 : pari_sp av = avma;
197 24532 : GEN nd = int_MSW(n);
198 24532 : ulong m = *nd;
199 : long i, j;
200 24532 : GEN v = utoipos(P), v1 = utoipos(P*P - 2);
201 :
202 24532 : if (m == 1)
203 1845 : j = 0;
204 : else
205 : {
206 22687 : j = 1+bfffo(m); /* < BIL */
207 22687 : m <<= j; j = BITS_IN_LONG - j;
208 : }
209 24532 : for (i=lgefint(n)-2;;) /* cf. leftright_pow */
210 : {
211 3307157 : for (; j; m<<=1,j--)
212 : { /* v = v_k, v1 = v_{k+1} */
213 3234640 : if (m&HIGHBIT)
214 : { /* set v = v_{2k+1}, v1 = v_{2k+2} */
215 1149299 : v = subiu(mulii(v,v1), P);
216 1149154 : v1= subiu(sqri(v1), 2);
217 : }
218 : else
219 : {/* set v = v_{2k}, v1 = v_{2k+1} */
220 2085341 : v1= subiu(mulii(v,v1), P);
221 2085186 : v = subiu(sqri(v), 2);
222 : }
223 3234332 : v = modii(v, N);
224 3234306 : v1= modii(v1,N);
225 3234362 : 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 72517 : if (--i == 0) return v;
232 48264 : j = BITS_IN_LONG;
233 48264 : nd=int_precW(nd);
234 48264 : 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 1032146 : u_LucasMod_pre(ulong n, ulong P, ulong N, ulong NI)
242 : {
243 : ulong v, v1, m;
244 : long j;
245 :
246 1032146 : if (n == 1) return P;
247 1032134 : j = 1 + bfffo(n); /* < BIL */
248 1032134 : v = P; v1 = P*P - 2;
249 1032134 : m = n<<j; j = BITS_IN_LONG - j;
250 63345362 : for (; j; m<<=1,j--)
251 : { /* v = v_k, v1 = v_{k+1} */
252 62330539 : if (m & HIGHBIT)
253 : { /* set v = v_{2k+1}, v1 = v_{2k+2} */
254 6303527 : v = Fl_sub(Fl_mul_pre(v,v1,N,NI), P, N);
255 6303542 : v1= Fl_sub(Fl_sqr_pre(v1,N,NI), 2UL, N);
256 : }
257 : else
258 : {/* set v = v_{2k}, v1 = v_{2k+1} */
259 56027012 : v1= Fl_sub(Fl_mul_pre(v,v1,N,NI),P, N);
260 56028207 : v = Fl_sub(Fl_sqr_pre(v,N,NI), 2UL, N);
261 : }
262 : }
263 1014823 : return v;
264 : }
265 :
266 : static ulong
267 1031950 : get_disc(ulong n)
268 : {
269 : ulong b;
270 : long i;
271 1031950 : for (b = 3, i = 0;; b += 2, i++)
272 1253320 : {
273 2285270 : ulong c = b*b - 4; /* = 1 mod 4 */
274 2285270 : if (krouu(n % c, c) < 0) break;
275 1253320 : if (i == 64 && uissquareall(n, &c)) return 0; /* oo loop if N = m^2 */
276 : }
277 1032152 : return b;
278 : }
279 : static int
280 1031952 : uislucaspsp_pre(ulong n, ulong ni)
281 : {
282 : long i, v;
283 1031952 : ulong b, z, m = n + 1;
284 :
285 1031952 : if (!m) return 0; /* neither 2^32-1 nor 2^64-1 are Lucas-pp */
286 1031952 : b = get_disc(n); if (!b) return 0;
287 1032148 : v = vals(m); m >>= v;
288 1032148 : z = u_LucasMod_pre(m, b, n, ni);
289 1032270 : if (z == 2 || z == n-2) return 1;
290 895405 : for (i=1; i<v; i++)
291 : {
292 895400 : if (!z) return 1;
293 471955 : z = Fl_sub(Fl_sqr_pre(z,n,ni), 2UL, n);
294 471954 : 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 24532 : islucaspsp(GEN N)
307 : {
308 24532 : pari_sp av = avma;
309 : GEN m, z;
310 : long i, v;
311 : ulong b;
312 :
313 24532 : for (b=3;; b+=2)
314 27890 : {
315 52422 : ulong c = b*b - 4; /* = 1 mod 4 */
316 52422 : if (b == 129 && Z_issquare(N)) return 0; /* avoid oo loop if N = m^2 */
317 52422 : if (krouu(umodiu(N,c), c) < 0) break;
318 : }
319 24532 : m = addiu(N,1); v = vali(m); m = shifti(m,-v);
320 24532 : z = LucasMod(m, b, N);
321 24532 : if (absequaliu(z, 2)) return 1;
322 19095 : if (equalii(z, subiu(N,2))) return 1;
323 36382 : for (i=1; i<v; i++)
324 : {
325 36161 : if (!signe(z)) return 1;
326 27354 : z = modii(subiu(sqri(z), 2), N);
327 27354 : if (absequaliu(z, 2)) return 0;
328 27354 : 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 928925 : _uispsp(ulong a, long n) { a %= n; return !a || uispsp(a, n); }
374 : static int
375 8306142 : _uisprime(ulong n)
376 : {
377 : #ifdef LONG_IS_64BIT
378 : ulong ni;
379 8204652 : if (n < 341531)
380 0 : return _uispsp(9345883071009581737UL, n);
381 8204652 : if (n < 1050535501)
382 543883 : return _uispsp(336781006125UL, n)
383 546305 : && _uispsp(9639812373923155UL, n);
384 7660769 : if (n < 350269456337)
385 37409 : return _uispsp(4230279247111683200UL, n)
386 14417 : && _uispsp(14694767155120705706UL, n)
387 51826 : && _uispsp(16641139526367750375UL, n);
388 : /* n & HIGHMASK */
389 7623360 : ni = get_Fl_red(n);
390 7625857 : return (uispsp_pre(2, n, ni) && uislucaspsp_pre(n,ni));
391 : #else
392 101490 : if (n < 360018361) return _uispsp(1143370UL, n) && _uispsp(2350307676UL, n);
393 73446 : return uispsp(15, n) && uispsp(176006322UL, n) && _uispsp(4221622697UL, n);
394 : #endif
395 : }
396 :
397 : int
398 73413647 : uisprime(ulong n)
399 : {
400 73413647 : if (!odd(n)) return n == 2;
401 43857325 : if (n <= maxprimelim()) return PRIMES_search(n) > 0;
402 : /* gcd-extraction is much slower */
403 13558528 : return n % 3 && n % 5 && n % 7 && n % 11 && n % 13 && n % 17
404 10219662 : && n % 19 && n % 23 && n % 29 && n % 31 && n % 37 && n % 41
405 28101878 : && _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 93396 : uisprime_661(ulong n)
420 : {
421 93396 : if (n <= maxprimelim()) return PRIMES_search(n) > 0;
422 77259 : if (n < 1016801) return n < 452929? 1: uispsp(2, n);
423 77259 : return _uisprime(n);
424 : }
425 :
426 : static int
427 373681 : iu_coprime(GEN N, ulong u) { return ugcd(u, umodiu(N,u)) == 1; }
428 : long
429 48281348 : BPSW_psp(GEN N)
430 : {
431 : pari_sp av;
432 48281348 : if (typ(N) != t_INT) pari_err_TYPE("BPSW_psp",N);
433 48281826 : if (signe(N) <= 0) return 0;
434 48220072 : if (lgefint(N) == 3) return uisprime(uel(N,2));
435 198672 : 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 176883 : if (!iu_coprime(N, 16294579238595022365UL) ||
440 109582 : !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 105510 : av = avma;
453 105510 : 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 45826 : Z_isanypower_nosmalldiv(GEN N, ulong cutoff, GEN *px)
460 : {
461 45826 : GEN x = N, y;
462 45826 : ulong mask = 7;
463 45826 : long ex, k = 1;
464 : forprime_t T;
465 57517 : while (Z_issquareall(x, &y)) { k <<= 1; x = y; }
466 45887 : while ( (ex = is_357_power(x, &y, &mask)) ) { k *= ex; x = y; }
467 45827 : (void)u_forprime_init(&T, 11, ULONG_MAX);
468 : /* stop when x^(1/k) < 2^cutoff */
469 45821 : while ( (ex = is_pth_power(x, &y, &T, cutoff)) ) { k *= ex; x = y; }
470 45826 : *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 49975 : BPSW_psp_nosmalldiv(GEN N)
477 : {
478 : pari_sp av;
479 49975 : long l = lgefint(N);
480 :
481 49975 : if (l == 3) return uisprime_661(uel(N,2));
482 32328 : av = avma;
483 : /* N large: test for pure power, rarely succeeds, but requires < 1% of
484 : * compositeness test times */
485 32328 : if (bit_accuracy(l) > 512 && Z_isanypower_nosmalldiv(N, 15, &N) != 1)
486 14 : return gc_long(av,0);
487 32313 : N = absi_shallow(N);
488 32313 : 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 3947208 : BPSW_isprime_small(GEN x)
502 : {
503 3947208 : long l = lgefint(x);
504 : #ifdef LONG_IS_64BIT
505 3880173 : return (l == 3);
506 : #else
507 67035 : 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 13640 : pl831(GEN N, GEN p)
517 : {
518 13640 : GEN b, c, g, Nmunp = diviiexact(subiu(N,1), p);
519 13640 : pari_sp av = avma;
520 : ulong a;
521 20028 : for(a = 2;; a++, set_avma(av))
522 : {
523 20028 : b = Fp_pow(utoipos(a), Nmunp, N);
524 20028 : if (!equali1(b)) break;
525 : }
526 13640 : c = Fp_pow(b,p,N);
527 13640 : g = gcdii(subiu(b,1), N); /* 0 < g < N */
528 13640 : 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 49 : BLS_test(GEN N, GEN f)
537 : {
538 : GEN c1, c2, r, q;
539 49 : q = dvmdii(N, f, &r);
540 49 : if (!is_pm1(r)) return 0;
541 49 : 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 49 : 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 4429 : BPSW_try_PL(GEN N)
553 : {
554 4429 : ulong B = minuu(1UL<<19, maxprime());
555 4429 : GEN E, p, U, F, N_1 = subiu(N,1);
556 4429 : GEN fa = absZ_factor_limit_strict(N_1, B, &U), P = gel(fa,1);
557 :
558 4429 : if (!U) return P; /* N-1 fully factored */
559 1707 : p = gel(U,1);
560 1707 : E = gel(fa,2);
561 1707 : U = powii(p, gel(U,2)); /* unfactored part of N-1 */
562 1707 : 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 1707 : if (cmpii(F, U) > 0) return P; /* 1/2-smooth */
566 1700 : if (cmpii(sqri(F), U) > 0) return BLS_test(N,F)? P: gen_0; /* 1/3-smooth */
567 1658 : return NULL; /* not smooth enough */
568 : }
569 :
570 : static GEN isprimePL(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 2771 : PL_certify(GEN N, GEN F)
575 : {
576 2771 : long i, l = lg(F);
577 16334 : for(i = 1; i < l; i++)
578 13563 : if (! pl831(N, gel(F,i))) return 0;
579 2771 : 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 98 : PL_certificate(GEN N, GEN F)
585 : {
586 98 : long i, l = lg(F);
587 : GEN C;
588 98 : if (BPSW_isprime_small(N)) return N;
589 98 : C = cgetg(l, t_VEC);
590 574 : for (i = 1; i < l; i++)
591 : {
592 476 : GEN p = gel(F,i), C0;
593 : ulong w;
594 :
595 476 : if (BPSW_isprime_small(p)) { gel(C,i) = p; continue; }
596 77 : w = pl831(N,p); if (!w) return gen_0;
597 77 : C0 = isprimePL(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 98 : 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 119 : isprimePL(GEN N)
660 : {
661 : GEN cbrtN, N_1, F, f;
662 119 : if (BPSW_isprime_small(N)) return N;
663 98 : cbrtN = sqrtnint(N,3);
664 98 : N_1 = subiu(N,1);
665 98 : F = Z_factor_until(N_1, sqri(cbrtN));
666 98 : f = factorback(F); /* factored part of N-1, f^3 > N */
667 98 : 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 98 : if (!equalii(f,N_1) && cmpii(sqri(f),N) <= 0 && !BLS_test(N,f))
675 0 : return gen_0; /* Failed, N is composite */
676 98 : F = gel(F,1); settyp(F, t_VEC);
677 98 : return PL_certificate(N, F);
678 : }
679 :
680 : /* assume N a BPSW pseudoprime, in particular, it is odd > 2. Prove N prime */
681 : long
682 3945868 : BPSW_isprime(GEN N)
683 : {
684 : pari_sp av;
685 : long t;
686 : GEN P;
687 3945868 : if (BPSW_isprime_small(N)) return 1;
688 4174 : av = avma; P = BPSW_try_PL(N);
689 4429 : if (!P) /* not smooth enough */
690 1658 : t = expi(N) < 768? isprimeAPRCL(N): isprimeECPP(N);
691 : else
692 2771 : t = (typ(P) == t_INT)? 0: PL_certify(N,P);
693 4429 : return gc_long(av,t);
694 : }
695 :
696 : static long
697 42 : _isprimePL(GEN x)
698 : {
699 42 : pari_sp av = avma;
700 42 : if (!BPSW_psp(x)) return 0;
701 35 : return gc_long(av, !isintzero(isprimePL(x)));
702 : }
703 : GEN
704 37289895 : gisprime(GEN x, long flag)
705 : {
706 37289895 : switch (flag)
707 : {
708 37301748 : case 0: return map_proto_lG(isprime,x);
709 28 : case 1: return map_proto_lG(_isprimePL,x);
710 14 : case 2: return map_proto_lG(isprimeAPRCL,x);
711 21 : case 3: return map_proto_lG(isprimeECPP,x);
712 : }
713 0 : pari_err_FLAG("gisprime");
714 : return NULL;/*LCOV_EXCL_LINE*/
715 : }
716 :
717 : long
718 38010850 : isprime(GEN x) { return BPSW_psp(x) && BPSW_isprime(x); }
719 :
720 : GEN
721 84 : primecert0(GEN x, long flag, long stopat)
722 : {
723 84 : if ((flag || typ(x) == t_INT) && !BPSW_psp(x)) return gen_0;
724 77 : switch(flag)
725 : {
726 70 : case 0: return ecpp0(x, stopat);
727 7 : case 1: { pari_sp av = avma; return gerepilecopy(av, isprimePL(x)); }
728 : }
729 0 : pari_err_FLAG("primecert");
730 : return NULL;/*LCOV_EXCL_LINE*/
731 : }
732 :
733 : GEN
734 0 : primecert(GEN x, long flag)
735 0 : { return primecert0(x, flag, 0); }
736 :
737 : enum { c_VOID = 0, c_ECPP, c_N1 };
738 : static long
739 98 : cert_type(GEN c)
740 : {
741 98 : switch(typ(c))
742 : {
743 0 : case t_INT: return c_ECPP;
744 98 : case t_VEC:
745 98 : if (lg(c) == 3 && typ(gel(c,1)) == t_INT) return c_N1;
746 84 : return c_ECPP;
747 : }
748 0 : return c_VOID;
749 : }
750 :
751 : long
752 56 : primecertisvalid(GEN c)
753 : {
754 56 : switch(typ(c))
755 : {
756 7 : case t_INT: return BPSW_isprime_small(c) && BPSW_psp(c);
757 49 : case t_VEC:
758 49 : return cert_type(c) == c_ECPP? ecppisvalid(c): PL_isvalid(c);
759 : }
760 0 : return 0;
761 : }
762 :
763 : static long
764 392 : check_ecppcertentry(GEN c)
765 : {
766 : GEN v;
767 392 : long i,l = lg(c);
768 392 : if (typ(c)!=t_VEC || l!=6) return 0;
769 1925 : for(i=1; i<=4; i++)
770 1540 : if (typ(gel(c,i))!=t_INT) return 0;
771 385 : v = gel(c,5);
772 385 : if(typ(v)!=t_VEC) return 0;
773 1155 : for(i=1; i<=2; i++)
774 770 : if (typ(gel(v,i))!=t_INT) return 0;
775 385 : return 1;
776 : }
777 :
778 : long
779 56 : check_ecppcert(GEN c)
780 : {
781 : long i, l;
782 56 : switch(typ(c))
783 : {
784 0 : case t_INT: return signe(c) >= 0;
785 56 : case t_VEC: break;
786 0 : default: return 0;
787 : }
788 56 : l = lg(c); if (l == 1) return 0;
789 434 : for (i = 1; i < l; i++)
790 392 : if (check_ecppcertentry(gel(c,i))==0) return 0;
791 42 : return 1;
792 : }
793 :
794 : GEN
795 49 : primecertexport(GEN c, long flag)
796 : {
797 49 : if (cert_type(c) != c_ECPP) pari_err_IMPL("N-1 certificate");
798 49 : if (!check_ecppcert(c))
799 14 : pari_err_TYPE("primecertexport - invalid certificate", c);
800 35 : return ecppexport(c, flag);
801 : }
802 :
803 : /***********************************************************************/
804 : /** **/
805 : /** PRIME NUMBERS **/
806 : /** **/
807 : /***********************************************************************/
808 :
809 : static struct {
810 : ulong p;
811 : long n;
812 : } prime_table[] = {
813 : { 0, 0},
814 : { 7919, 1000},
815 : { 17389, 2000},
816 : { 27449, 3000},
817 : { 37813, 4000},
818 : { 48611, 5000},
819 : { 59359, 6000},
820 : { 70657, 7000},
821 : { 81799, 8000},
822 : { 93179, 9000},
823 : { 104729, 10000},
824 : { 224737, 20000},
825 : { 350377, 30000},
826 : { 479909, 40000},
827 : { 611953, 50000},
828 : { 746773, 60000},
829 : { 882377, 70000},
830 : { 1020379, 80000},
831 : { 1159523, 90000},
832 : { 1299709, 100000},
833 : { 2750159, 200000},
834 : { 7368787, 500000},
835 : { 15485863, 1000000},
836 : { 32452843, 2000000},
837 : { 86028121, 5000000},
838 : { 179424673, 10000000},
839 : { 373587883, 20000000},
840 : { 982451653, 50000000},
841 : { 2038074743, 100000000},
842 : { 4000000483UL,189961831},
843 : { 4222234741UL,200000000},
844 : #if BITS_IN_LONG == 64
845 : { 5336500537UL, 250000000L},
846 : { 6461335109UL, 300000000L},
847 : { 7594955549UL, 350000000L},
848 : { 8736028057UL, 400000000L},
849 : { 9883692017UL, 450000000L},
850 : { 11037271757UL, 500000000L},
851 : { 13359555403UL, 600000000L},
852 : { 15699342107UL, 700000000L},
853 : { 18054236957UL, 800000000L},
854 : { 20422213579UL, 900000000L},
855 : { 22801763489UL, 1000000000L},
856 : { 47055833459UL, 2000000000L},
857 : { 71856445751UL, 3000000000L},
858 : { 97011687217UL, 4000000000L},
859 : {122430513841UL, 5000000000L},
860 : {148059109201UL, 6000000000L},
861 : {173862636221UL, 7000000000L},
862 : {200000000507UL, 8007105083L},
863 : {225898512559UL, 9000000000L},
864 : {252097800623UL,10000000000L},
865 : {384489816343UL,15000000000L},
866 : {518649879439UL,20000000000L},
867 : {654124187867UL,25000000000L},
868 : {790645490053UL,30000000000L},
869 : {928037044463UL,35000000000L},
870 : {1066173339601UL,40000000000L},
871 : {1344326694119UL,50000000000L},
872 : {1624571841097UL,60000000000L},
873 : {1906555030411UL,70000000000L},
874 : {2190026988349UL,80000000000L},
875 : {2474799787573UL,90000000000L},
876 : {2760727302517UL,100000000000L}
877 : #endif
878 : };
879 : static const int prime_table_len = numberof(prime_table);
880 :
881 : /* find prime closest to n in prime_table. */
882 : static long
883 97 : prime_table_closest_p(ulong n)
884 : {
885 : long i;
886 2209 : for (i = 1; i < prime_table_len; i++)
887 : {
888 2209 : ulong p = prime_table[i].p;
889 2209 : if (p > n)
890 : {
891 97 : ulong u = n - prime_table[i-1].p;
892 97 : if (p - n > u) i--;
893 97 : break;
894 : }
895 : }
896 97 : if (i == prime_table_len) i = prime_table_len - 1;
897 97 : return i;
898 : }
899 :
900 : /* return the n-th successor of prime p > 2; n > 0 */
901 : static GEN
902 56 : prime_successor(ulong p, ulong n)
903 : {
904 56 : GEN Q = utoipos(p), P = NULL;
905 : ulong i;
906 56 : RESET:
907 : for(;;)
908 : {
909 : forprime_t S;
910 56 : GEN L = dbltor(4*n * log(gtodouble(Q) + 1)), R = addii(Q, ceil_safe(L));
911 :
912 56 : forprime_init(&S, addiu(Q, 1), R);
913 56 : Q = R;
914 2325176 : for (i = 1; i <= n; i++)
915 : {
916 2325120 : P = forprime_next(&S);
917 2325120 : if (!P) { n -= i-1; goto RESET; }
918 : }
919 56 : return P;
920 : }
921 : }
922 : /* find the N-th prime */
923 : static GEN
924 273 : prime_table_find_n(ulong N)
925 : {
926 : ulong n, p, maxp;
927 : long i;
928 273 : if (N <= (ulong)pari_PRIMES[0]) return utoipos((ulong)pari_PRIMES[N]);
929 : /* not in table */
930 1204 : for (i = 1; i < prime_table_len; i++)
931 : {
932 1204 : n = prime_table[i].n;
933 1204 : if (n > N)
934 : {
935 56 : ulong u = N - prime_table[i-1].n;
936 56 : if (n - N > u) i--;
937 56 : break;
938 : }
939 : }
940 56 : if (i == prime_table_len) i = prime_table_len - 1;
941 56 : p = prime_table[i].p;
942 56 : n = prime_table[i].n;
943 56 : if (n > N)
944 : {
945 0 : i--;
946 0 : p = prime_table[i].p;
947 0 : n = prime_table[i].n;
948 : }
949 : /* n <= N */
950 56 : if (n == N) return utoipos(p);
951 56 : maxp = maxprime();
952 56 : if (p < maxp) { p = maxp; n = pari_PRIMES[0]; }
953 56 : return prime_successor(p, N - n);
954 : }
955 :
956 : ulong
957 0 : uprime(long N)
958 : {
959 0 : pari_sp av = avma;
960 : GEN p;
961 0 : if (N <= 0) pari_err_DOMAIN("prime", "n", "<=",gen_0, stoi(N));
962 0 : p = prime_table_find_n(N);
963 0 : if (lgefint(p) != 3) pari_err_OVERFLOW("uprime");
964 0 : return gc_ulong(av, p[2]);
965 : }
966 : GEN
967 280 : prime(long N)
968 : {
969 280 : pari_sp av = avma;
970 : GEN p;
971 280 : if (N <= 0) pari_err_DOMAIN("prime", "n", "<=",gen_0, stoi(N));
972 273 : new_chunk(4); /*HACK*/
973 273 : p = prime_table_find_n(N);
974 273 : set_avma(av); return icopy(p);
975 : }
976 :
977 : /* N encodes an interval in which to draw a random uniform prime;
978 : * decoded as [a, b], d = b-a+1 */
979 : static void
980 147 : prime_interval(GEN N, GEN *pa, GEN *pb, GEN *pd)
981 : {
982 : GEN a, b, d;
983 147 : switch(typ(N))
984 : {
985 91 : case t_INT:
986 91 : a = gen_2;
987 91 : b = subiu(N,1); /* between 2 and N-1 */
988 91 : d = subiu(N,2);
989 91 : if (signe(d) <= 0)
990 14 : pari_err_DOMAIN("randomprime","N", "<=", gen_2, N);
991 77 : break;
992 56 : case t_VEC:
993 56 : if (lg(N) != 3) pari_err_TYPE("randomprime",N);
994 56 : a = gel(N,1);
995 56 : b = gel(N,2);
996 56 : if (gcmp(b, a) < 0)
997 7 : pari_err_DOMAIN("randomprime","b-a", "<", gen_0, mkvec2(a,b));
998 49 : if (typ(a) != t_INT)
999 : {
1000 7 : a = gceil(a);
1001 7 : if (typ(a) != t_INT) pari_err_TYPE("randomprime",a);
1002 : }
1003 49 : if (typ(b) != t_INT)
1004 : {
1005 7 : b = gfloor(b);
1006 7 : if (typ(b) != t_INT) pari_err_TYPE("randomprime",b);
1007 : }
1008 49 : if (cmpiu(a, 2) < 0)
1009 : {
1010 7 : a = gen_2;
1011 7 : d = subiu(b,1);
1012 : }
1013 : else
1014 42 : d = addiu(subii(b,a), 1);
1015 49 : if (signe(d) <= 0)
1016 14 : pari_err_DOMAIN("randomprime","floor(b) - max(ceil(a),2)", "<",
1017 : gen_0, mkvec2(a,b));
1018 35 : break;
1019 0 : default:
1020 0 : pari_err_TYPE("randomprime", N);
1021 : a = b = d = NULL; /*LCOV_EXCL_LINE*/
1022 : }
1023 112 : *pa = a; *pb = b; *pd = d;
1024 112 : }
1025 :
1026 : /* random prime */
1027 : GEN
1028 91 : randomprime(GEN N)
1029 : {
1030 91 : pari_sp av = avma, av2;
1031 : GEN a, b, d;
1032 91 : if (!N)
1033 : for(;;)
1034 56 : {
1035 63 : ulong p = random_bits(31);
1036 63 : if (uisprime(p)) return utoipos(p);
1037 : }
1038 84 : prime_interval(N, &a, &b, &d); av2 = avma;
1039 : for (;;)
1040 1685 : {
1041 1741 : GEN p = addii(a, randomi(d));
1042 1741 : if (BPSW_psp(p)) return gerepileuptoint(av, p);
1043 1685 : set_avma(av2);
1044 : }
1045 : }
1046 : GEN
1047 154 : randomprime0(GEN N, GEN q)
1048 : {
1049 154 : pari_sp av = avma, av2;
1050 : GEN C, D, a, b, d, r;
1051 154 : if (!q) return randomprime(N);
1052 63 : switch(typ(q))
1053 : {
1054 42 : case t_INT: C = gen_1; D = q; break;
1055 21 : case t_INTMOD: C = gel(q,2); D = gel(q,1); break;
1056 0 : default:
1057 0 : pari_err_TYPE("randomprime", q);
1058 : return NULL;/*LCOV_EXCL_LINE*/
1059 : }
1060 63 : if (!N) N = int2n(31);
1061 63 : prime_interval(N, &a, &b, &d);
1062 56 : r = modii(subii(C, a), D);
1063 56 : if (signe(r))
1064 : {
1065 49 : a = addii(a, r);
1066 49 : if (cmpii(a,b) > 0)
1067 7 : pari_err(e_MISC, "no prime satisfies congruence in interval");
1068 42 : d = subii(d, r);
1069 : }
1070 49 : if (!equali1(gcdii(C,D)))
1071 : {
1072 14 : if (isprime(a)) return gerepilecopy(av, a);
1073 7 : pari_err_COPRIME("randomprime", C, D);
1074 : }
1075 35 : d = divii(d, D); if (!signe(d)) d = gen_1;
1076 35 : av2 = avma;
1077 : for (;;)
1078 3584 : {
1079 3619 : GEN p = addii(a, mulii(D, randomi(d)));
1080 3619 : if (BPSW_psp(p)) return gerepileuptoint(av, p);
1081 3584 : set_avma(av2);
1082 : }
1083 : return NULL;
1084 : }
1085 :
1086 : /* cf gen_search; assume 0 <= x <= maxprimelim(). Return i > 0 such that
1087 : * x = T[i] or -i < 0 such that T[i-1] < x < T[i]; return -1 for the special
1088 : * case 0 <= x < 2 */
1089 : long
1090 57368859 : PRIMES_search(ulong x)
1091 : {
1092 57368859 : pari_prime *T = pari_PRIMES;
1093 57368859 : ulong i, u = minuu(T[0], (x + 2) >> ((x < 122UL)? 1: 2)), l = 1;
1094 : do
1095 : {
1096 426233046 : i = (l+u) >> 1;
1097 426233046 : if (x < T[i]) u = i-1; else if (x > T[i]) l = i+1; else return i;
1098 407292919 : } while (u > l);
1099 38428132 : if (u == l) { i = l; if (x == T[i]) return i; }
1100 29149650 : return -(long)(T[i] > x? i: i+1);
1101 : }
1102 :
1103 : ulong
1104 18423 : uprimepi(ulong a)
1105 : {
1106 : ulong maxp, p, n;
1107 : forprime_t S;
1108 : long i;
1109 18423 : if (a <= maxprimelim())
1110 : {
1111 18326 : i = PRIMES_search(a);
1112 18326 : return (ulong)(i > 0? i: -i-1);
1113 : }
1114 : /* a > maxprimelim >= maxprime */
1115 97 : i = prime_table_closest_p(a);
1116 97 : p = prime_table[i].p;
1117 97 : if (p > a) p = prime_table[--i].p;
1118 : /* p = largest prime in prime_table <= a and a > maxprime */
1119 97 : maxp = maxprime(); /* replace p by maxprime if larger */
1120 97 : if (p > maxp) n = prime_table[i].n; else { p = maxp; n = maxprimeN(); }
1121 97 : (void)u_forprime_init(&S, p+1, a);
1122 52092780 : for (; p; n++) p = u_forprime_next(&S);
1123 97 : return n-1;
1124 : }
1125 :
1126 : GEN
1127 252 : primepi(GEN x)
1128 : {
1129 252 : pari_sp av = avma;
1130 252 : GEN pp, nn, N = typ(x) == t_INT? x: gfloor(x);
1131 : forprime_t S;
1132 : ulong n, p;
1133 : long i;
1134 252 : if (typ(N) != t_INT) pari_err_TYPE("primepi",N);
1135 252 : if (signe(N) <= 0) return gen_0;
1136 252 : if (lgefint(N) == 3) return gc_utoi(av, uprimepi(N[2]));
1137 1 : i = prime_table_len-1;
1138 1 : p = prime_table[i].p;
1139 1 : n = prime_table[i].n;
1140 1 : (void)forprime_init(&S, utoipos(p+1), N);
1141 1 : nn = setloop(utoipos(n));
1142 1 : pp = gen_0;
1143 3280223 : for (; pp; incloop(nn)) pp = forprime_next(&S);
1144 1 : return gerepileuptoint(av, subiu(nn,1));
1145 : }
1146 :
1147 : /* pi(x) < x/log x * (1 + 1/log x + 2.51/log^2 x)), x>=355991 [ Dusart ]
1148 : * pi(x) < x/(log x - 1.1), x >= 60184 [ Dusart ]
1149 : * ? \p9
1150 : * ? M = 0; for(x = 4, 60184, M = max(M, log(x) - x/primepi(x))); M
1151 : * %1 = 1.11196252 */
1152 : double
1153 402995 : primepi_upper_bound(double x)
1154 : {
1155 402995 : if (x >= 355991)
1156 : {
1157 1890 : double L = 1/log(x);
1158 1890 : return x * L * (1 + L + 2.51*L*L);
1159 : }
1160 401105 : if (x >= 60184) return x / (log(x) - 1.1);
1161 401105 : if (x < 5) return 2; /* don't bother */
1162 273778 : return x / (log(x) - 1.111963);
1163 : }
1164 : /* pi(x) > x/log x (1 + 1/log x), x >= 599 [ Dusart ]
1165 : * pi(x) > x / (log x + 2), x >= 55 [ Rosser ] */
1166 : double
1167 14 : primepi_lower_bound(double x)
1168 : {
1169 14 : if (x >= 599)
1170 : {
1171 14 : double L = 1/log(x);
1172 14 : return x * L * (1 + L);
1173 : }
1174 0 : if (x < 55) return 0; /* don't bother */
1175 0 : return x / (log(x) + 2.);
1176 : }
1177 : GEN
1178 1 : gprimepi_upper_bound(GEN x)
1179 : {
1180 1 : pari_sp av = avma;
1181 : double L;
1182 1 : if (typ(x) != t_INT) x = gfloor(x);
1183 1 : if (expi(x) <= 1022)
1184 : {
1185 1 : set_avma(av);
1186 1 : return dbltor(primepi_upper_bound(gtodouble(x)));
1187 : }
1188 0 : x = itor(x, LOWDEFAULTPREC);
1189 0 : L = 1 / rtodbl(logr_abs(x));
1190 0 : x = mulrr(x, dbltor(L * (1 + L + 2.51*L*L)));
1191 0 : return gerepileuptoleaf(av, x);
1192 : }
1193 : GEN
1194 1 : gprimepi_lower_bound(GEN x)
1195 : {
1196 1 : pari_sp av = avma;
1197 : double L;
1198 1 : if (typ(x) != t_INT) x = gfloor(x);
1199 1 : if (abscmpiu(x, 55) <= 0) return gen_0;
1200 1 : if (expi(x) <= 1022)
1201 : {
1202 1 : set_avma(av);
1203 1 : return dbltor(primepi_lower_bound(gtodouble(x)));
1204 : }
1205 0 : x = itor(x, LOWDEFAULTPREC);
1206 0 : L = 1 / rtodbl(logr_abs(x));
1207 0 : x = mulrr(x, dbltor(L * (1 + L)));
1208 0 : return gerepileuptoleaf(av, x);
1209 : }
1210 :
1211 : GEN
1212 105 : primes(long n)
1213 : {
1214 : forprime_t S;
1215 : long i;
1216 : GEN y;
1217 105 : if (n <= 0) return cgetg(1, t_VEC);
1218 105 : y = cgetg(n+1, t_VEC);
1219 105 : (void)new_chunk(3*n); /*HACK*/
1220 105 : u_forprime_init(&S, 2, (ulong)n > maxprimeN()? ULONG_MAX: maxprime());
1221 105 : set_avma((pari_sp)y);
1222 6069 : for (i = 1; i <= n; i++) gel(y, i) = utoipos( u_forprime_next(&S) );
1223 105 : return y;
1224 : }
1225 : GEN
1226 91 : primes_zv(long n)
1227 : {
1228 : forprime_t S;
1229 : long i;
1230 : GEN y;
1231 91 : if (n <= 0) return cgetg(1, t_VECSMALL);
1232 91 : y = cgetg(n+1,t_VECSMALL);
1233 : /* don't initialize sieve if all fits in primetable */
1234 91 : u_forprime_init(&S, 2, (ulong)n > maxprimeN()? ULONG_MAX: maxprime());
1235 3731 : for (i = 1; i <= n; i++) y[i] = u_forprime_next(&S);
1236 91 : set_avma((pari_sp)y); return y;
1237 : }
1238 : GEN
1239 231 : primes0(GEN N)
1240 : {
1241 231 : switch(typ(N))
1242 : {
1243 105 : case t_INT: return primes(itos(N));
1244 126 : case t_VEC:
1245 126 : if (lg(N) == 3) return primes_interval(gel(N,1),gel(N,2));
1246 : }
1247 0 : pari_err_TYPE("primes", N);
1248 : return NULL;/*LCOV_EXCL_LINE*/
1249 : }
1250 :
1251 : GEN
1252 189 : primes_interval(GEN a, GEN b)
1253 : {
1254 189 : pari_sp av = avma;
1255 : forprime_t S;
1256 : long i, n;
1257 : GEN y, d, p;
1258 189 : if (typ(a) != t_INT)
1259 : {
1260 7 : a = gceil(a);
1261 7 : if (typ(a) != t_INT) pari_err_TYPE("primes_interval",a);
1262 : }
1263 189 : if (typ(b) != t_INT)
1264 : {
1265 14 : b = gfloor(b);
1266 14 : if (typ(b) != t_INT) pari_err_TYPE("primes_interval",b);
1267 : }
1268 182 : if (signe(a) < 0) a = gen_2;
1269 182 : d = subii(b, a);
1270 182 : if (signe(d) < 0 || signe(b) <= 0) { set_avma(av); return cgetg(1, t_VEC); }
1271 182 : if (lgefint(b) == 3)
1272 : {
1273 156 : set_avma(av);
1274 156 : y = primes_interval_zv(itou(a), itou(b));
1275 156 : n = lg(y); settyp(y, t_VEC);
1276 470082 : for (i = 1; i < n; i++) gel(y,i) = utoipos(y[i]);
1277 156 : return y;
1278 : }
1279 : /* at most d+1 primes in [a,b]. If d large, try better bound to lower
1280 : * memory use */
1281 26 : if (abscmpiu(d,100000) > 0)
1282 : {
1283 1 : GEN D = gsub(gprimepi_upper_bound(b), gprimepi_lower_bound(a));
1284 1 : D = ceil_safe(D);
1285 1 : if (cmpii(D, d) < 0) d = D;
1286 : }
1287 26 : n = itos(d)+1;
1288 26 : forprime_init(&S, a, b);
1289 26 : y = cgetg(n+1, t_VEC); i = 1;
1290 5995 : while ((p = forprime_next(&S))) gel(y, i++) = icopy(p);
1291 26 : setlg(y, i); return gerepileupto(av, y);
1292 : }
1293 :
1294 : /* simpler case a <= b <= maxprimelim() */
1295 : static GEN
1296 14465 : PRIMES_interval(ulong a, ulong b)
1297 : {
1298 14465 : long k, i = PRIMES_search(a), j = PRIMES_search(b);
1299 : GEN y;
1300 14465 : if (i < 0) i = -i;
1301 14465 : if (j < 0) j = -j - 1;
1302 : /* T[i] = smallest p >= a, T[j] = largest p <= b */
1303 14465 : y = cgetg(j - i + 2, t_VECSMALL);
1304 14465 : if (i == 1)
1305 3462809 : for (; i <= j; i++) y[i] = pari_PRIMES[i];
1306 : else
1307 151391158 : for (k = 1; i <= j; i++) y[k++] = pari_PRIMES[i];
1308 14465 : return y;
1309 : }
1310 : /* a <= b, at most d primes in [a,b]. Return them.
1311 : * PRIMES_interval is more efficient if b <= maxprimelim() */
1312 : static GEN
1313 79 : primes_interval_i(ulong a, ulong b, ulong d)
1314 : {
1315 79 : ulong p, i = 1, n = d+1;
1316 79 : GEN y = cgetg(n+1, t_VECSMALL);
1317 79 : pari_sp av = avma;
1318 : forprime_t S;
1319 79 : u_forprime_init(&S, a, b);
1320 10739901 : while ((p = u_forprime_next(&S))) y[i++] = p;
1321 79 : set_avma(av); setlg(y, i); stackdummy((pari_sp)(y + i), (pari_sp)(y + n+1));
1322 79 : return y;
1323 : }
1324 : GEN
1325 14355 : primes_interval_zv(ulong a, ulong b)
1326 : {
1327 : ulong d;
1328 14355 : if (a < 3) return primes_upto_zv(b);
1329 14250 : if (b < a) return cgetg(1, t_VECSMALL);
1330 14250 : if (b <= maxprimelim()) return PRIMES_interval(a, b);
1331 44 : d = b - a;
1332 44 : if (d > 100000UL)
1333 : {
1334 13 : ulong D = (ulong)ceil(primepi_upper_bound(b)-primepi_lower_bound(a));
1335 13 : if (D < d) d = D;
1336 : }
1337 44 : return primes_interval_i(a, b, d);
1338 : }
1339 : GEN
1340 294 : primes_upto_zv(ulong b)
1341 : {
1342 : ulong d;
1343 294 : if (b < 2) return cgetg(1, t_VECSMALL);
1344 294 : if (b <= maxprimelim()) return PRIMES_interval(2, b);
1345 35 : d = (b > 100000UL)? (ulong)primepi_upper_bound(b): b;
1346 35 : return primes_interval_i(2, b, d);
1347 : }
1348 :
1349 : /***********************************************************************/
1350 : /** **/
1351 : /** PRIVATE PRIME TABLE **/
1352 : /** **/
1353 : /***********************************************************************/
1354 :
1355 : void
1356 324287 : pari_set_primetab(GEN global_primetab)
1357 : {
1358 324287 : if (global_primetab)
1359 : {
1360 322438 : long i, l = lg(global_primetab);
1361 322438 : primetab = cgetg_block(l, t_VEC);
1362 352879 : for (i = 1; i < l; i++)
1363 30468 : gel(primetab,i) = gclone(gel(global_primetab,i));
1364 1849 : } else primetab = cgetg_block(1, t_VEC);
1365 324262 : }
1366 :
1367 : /* delete dummy NULL entries */
1368 : static void
1369 21 : cleanprimetab(GEN T)
1370 : {
1371 21 : long i,j, l = lg(T);
1372 70 : for (i = j = 1; i < l; i++)
1373 49 : if (T[i]) T[j++] = T[i];
1374 21 : setlg(T,j);
1375 21 : }
1376 : /* remove p from T */
1377 : static void
1378 28 : rmprime(GEN T, GEN p)
1379 : {
1380 : long i;
1381 28 : if (typ(p) != t_INT) pari_err_TYPE("removeprimes",p);
1382 28 : i = ZV_search(T, p);
1383 28 : if (!i)
1384 7 : pari_err_DOMAIN("removeprimes","prime","not in",
1385 : strtoGENstr("primetable"), p);
1386 21 : gunclone(gel(T,i)); gel(T,i) = NULL;
1387 21 : cleanprimetab(T);
1388 21 : }
1389 :
1390 : /*stolen from ZV_union_shallow() : clone entries from y */
1391 : static GEN
1392 56 : addp_union(GEN x, GEN y)
1393 : {
1394 56 : long i, j, k, lx = lg(x), ly = lg(y);
1395 56 : GEN z = cgetg(lx + ly - 1, t_VEC);
1396 56 : i = j = k = 1;
1397 84 : while (i<lx && j<ly)
1398 : {
1399 28 : int s = cmpii(gel(x,i), gel(y,j));
1400 28 : if (s < 0)
1401 21 : gel(z,k++) = gel(x,i++);
1402 7 : else if (s > 0)
1403 0 : gel(z,k++) = gclone(gel(y,j++));
1404 : else {
1405 7 : gel(z,k++) = gel(x,i++);
1406 7 : j++;
1407 : }
1408 : }
1409 56 : while (i<lx) gel(z,k++) = gel(x,i++);
1410 147 : while (j<ly) gel(z,k++) = gclone(gel(y,j++));
1411 56 : setlg(z, k); return z;
1412 : }
1413 :
1414 : /* p is NULL, or a single element or a row vector with "primes" to add to
1415 : * prime table. */
1416 : static GEN
1417 189 : addp(GEN *T, GEN p)
1418 : {
1419 189 : pari_sp av = avma;
1420 : long i, l;
1421 : GEN v;
1422 :
1423 189 : if (!p || lg(p) == 1) return *T;
1424 70 : if (!is_vec_t(typ(p))) p = mkvec(p);
1425 :
1426 70 : RgV_check_ZV(p, "addprimes");
1427 63 : v = gen_indexsort_uniq(p, (void*)&cmpii, &cmp_nodata);
1428 63 : p = vecpermute(p, v);
1429 63 : if (abscmpiu(gel(p,1), 2) < 0) pari_err_DOMAIN("addprimes", "p", "<", gen_2,p);
1430 56 : p = addp_union(*T, p);
1431 56 : l = lg(p);
1432 56 : if (l != lg(*T))
1433 : {
1434 56 : GEN old = *T, t = cgetg_block(l, t_VEC);
1435 175 : for (i = 1; i < l; i++) gel(t,i) = gel(p,i);
1436 56 : *T = t; gunclone(old);
1437 : }
1438 56 : set_avma(av); return *T;
1439 : }
1440 : GEN
1441 189 : addprimes(GEN p) { return addp(&primetab, p); }
1442 :
1443 : static GEN
1444 28 : rmprimes(GEN T, GEN prime)
1445 : {
1446 : long i,tx;
1447 :
1448 28 : if (!prime) return T;
1449 28 : tx = typ(prime);
1450 28 : if (is_vec_t(tx))
1451 : {
1452 14 : if (prime == T)
1453 : {
1454 14 : for (i=1; i < lg(prime); i++) gunclone(gel(prime,i));
1455 7 : setlg(prime, 1);
1456 : }
1457 : else
1458 : {
1459 21 : for (i=1; i < lg(prime); i++) rmprime(T, gel(prime,i));
1460 : }
1461 14 : return T;
1462 : }
1463 14 : rmprime(T, prime); return T;
1464 : }
1465 : GEN
1466 28 : removeprimes(GEN prime) { return rmprimes(primetab, prime); }
|
