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