Code coverage tests

This page documents the degree to which the PARI/GP source code is tested by our public test suite, distributed with the source distribution in directory src/test/. This is measured by the gcov utility; we then process gcov output using the lcov frond-end.

We test a few variants depending on Configure flags on the pari.math.u-bordeaux.fr machine (x86_64 architecture), and agregate them in the final report:

The target is to exceed 90% coverage for all mathematical modules (given that branches depending on DEBUGLEVEL or DEBUGMEM are not covered). This script is run to produce the results below.

LCOV - code coverage report
Current view: top level - basemath - prime.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.16.2 lcov report (development 29454-e7b6e21fb9) Lines: 640 703 91.0 %
Date: 2024-07-18 09:03:47 Functions: 70 76 92.1 %
Legend: Lines: hit not hit

          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      137840 : is2psp(GEN n)
      81             : {
      82      137840 :   GEN c, t = subiu(n, 1);
      83      137836 :   pari_sp av = avma;
      84      137836 :   long e = vali(t);
      85             : 
      86      137836 :   c = Fp_pow(gen_2, shifti(t, -e), n);
      87      137847 :   if (is_pm1(c) || equalii(t, c)) return 1;
      88      247054 :   while (--e)
      89             :   { /* go fishing for -1, not for 1 (e - 1 squaring) */
      90      133740 :     c = remii(sqri(c), n);
      91      133740 :     if (equalii(t, c)) return 1;
      92             :     /* can return 0 if (c == 1) but very infrequent */
      93      122370 :     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      113314 :   return 0;
     100             : }
     101             : static int
     102     7838608 : uispsp_pre(ulong a, ulong n, ulong ni)
     103             : {
     104     7838608 :   ulong c, t = n - 1;
     105     7838608 :   long e = vals(t);
     106             : 
     107     7839212 :   c = Fl_powu_pre(a, t >> e, n, ni);
     108     7827998 :   if (c == 1 || c == t) return 1;
     109    14257405 :   while (--e)
     110             :   { /* go fishing for -1, not for 1 (saves one squaring) */
     111     7591794 :     c = Fl_sqr_pre(c, n, ni);
     112     7598526 :     if (c == t) return 1;
     113             :     /* can return 0 if (c == 1) but very infrequent */
     114             :   }
     115     6665611 :   return 0;
     116             : }
     117             : int
     118     1030398 : uispsp(ulong a, ulong n)
     119             : {
     120             :   ulong c, t;
     121             :   long e;
     122             : 
     123     1030398 :   if (n & HIGHMASK) return uispsp_pre(a, n, get_Fl_red(n));
     124      820474 :   t = n - 1;
     125      820474 :   e = vals(t);
     126      820251 :   c = Fl_powu(a, t >> e, n);
     127      824180 :   if (c == 1 || c == t) return 1;
     128      923148 :   while (--e)
     129             :   { /* go fishing for -1, not for 1 (e - 1 squaring) */
     130      617171 :     c = Fl_sqr(c, n);
     131      617059 :     if (c == t) return 1;
     132             :     /* can return 0 if (c == 1) but very infrequent */
     133             :   }
     134      305977 :   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     3307079 :     for (; j; m<<=1,j--)
     212             :     { /* v = v_k, v1 = v_{k+1} */
     213     3234584 :       if (m&HIGHBIT)
     214             :       { /* set v = v_{2k+1}, v1 = v_{2k+2} */
     215     1149275 :         v = subiu(mulii(v,v1), P);
     216     1149161 :         v1= subiu(sqri(v1), 2);
     217             :       }
     218             :       else
     219             :       {/* set v = v_{2k}, v1 = v_{2k+1} */
     220     2085309 :         v1= subiu(mulii(v,v1), P);
     221     2085148 :         v = subiu(sqri(v), 2);
     222             :       }
     223     3234345 :       v = modii(v, N);
     224     3234253 :       v1= modii(v1,N);
     225     3234302 :       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       72495 :     if (--i == 0) return v;
     232       48265 :     j = BITS_IN_LONG;
     233       48265 :     nd=int_precW(nd);
     234       48265 :     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     1032375 : u_LucasMod_pre(ulong n, ulong P, ulong N, ulong NI)
     242             : {
     243             :   ulong v, v1, m;
     244             :   long j;
     245             : 
     246     1032375 :   if (n == 1) return P;
     247     1032363 :   j = 1 + bfffo(n); /* < BIL */
     248     1032363 :   v = P; v1 = P*P - 2;
     249     1032363 :   m = n<<j; j = BITS_IN_LONG - j;
     250    63341664 :   for (; j; m<<=1,j--)
     251             :   { /* v = v_k, v1 = v_{k+1} */
     252    62328927 :     if (m & HIGHBIT)
     253             :     { /* set v = v_{2k+1}, v1 = v_{2k+2} */
     254     6304343 :       v = Fl_sub(Fl_mul_pre(v,v1,N,NI), P, N);
     255     6304310 :       v1= Fl_sub(Fl_sqr_pre(v1,N,NI), 2UL, N);
     256             :     }
     257             :     else
     258             :     {/* set v = v_{2k}, v1 = v_{2k+1} */
     259    56024584 :       v1= Fl_sub(Fl_mul_pre(v,v1,N,NI),P, N);
     260    56022134 :       v = Fl_sub(Fl_sqr_pre(v,N,NI), 2UL, N);
     261             :     }
     262             :   }
     263     1012737 :   return v;
     264             : }
     265             : 
     266             : static ulong
     267     1032159 : get_disc(ulong n)
     268             : {
     269             :   ulong b;
     270             :   long i;
     271     1032159 :   for (b = 3, i = 0;; b += 2, i++)
     272     1253822 :   {
     273     2285981 :     ulong c = b*b - 4; /* = 1 mod 4 */
     274     2285981 :     if (krouu(n % c, c) < 0) break;
     275     1253822 :     if (i == 64 && uissquareall(n, &c)) return 0; /* oo loop if N = m^2 */
     276             :   }
     277     1032380 :   return b;
     278             : }
     279             : static int
     280     1032165 : uislucaspsp_pre(ulong n, ulong ni)
     281             : {
     282             :   long i, v;
     283     1032165 :   ulong b, z, m = n + 1;
     284             : 
     285     1032165 :   if (!m) return 0; /* neither 2^32-1 nor 2^64-1 are Lucas-pp */
     286     1032165 :   b = get_disc(n); if (!b) return 0;
     287     1032380 :   v = vals(m); m >>= v;
     288     1032370 :   z = u_LucasMod_pre(m, b, n, ni);
     289     1032489 :   if (z == 2 || z == n-2) return 1;
     290      895805 :   for (i=1; i<v; i++)
     291             :   {
     292      895798 :     if (!z) return 1;
     293      472244 :     z = Fl_sub(Fl_sqr_pre(z,n,ni), 2UL, n);
     294      472244 :     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       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      928907 : _uispsp(ulong a, long n) { a %= n; return !a || uispsp(a, n); }
     374             : static int
     375     8308109 : _uisprime(ulong n)
     376             : {
     377             : #ifdef LONG_IS_64BIT
     378             :   ulong ni;
     379     8206619 :   if (n < 341531)
     380           0 :     return _uispsp(9345883071009581737UL, n);
     381     8206619 :   if (n < 1050535501)
     382      544004 :     return _uispsp(336781006125UL, n)
     383      546127 :        && _uispsp(9639812373923155UL, n);
     384     7662615 :   if (n < 350269456337)
     385       37409 :     return _uispsp(4230279247111683200UL, n)
     386       14417 :         && _uispsp(14694767155120705706UL, n)
     387       51826 :         && _uispsp(16641139526367750375UL, n);
     388             :   /* n & HIGHMASK */
     389     7625206 :   ni = get_Fl_red(n);
     390     7627906 :   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    73480257 : uisprime(ulong n)
     399             : {
     400    73480257 :   if (!odd(n)) return n == 2;
     401    43883595 :   if (n <= maxprimelim()) return PRIMES_search(n) > 0;
     402             :   /* gcd-extraction is much slower */
     403    13560959 :   return n % 3 && n % 5 && n % 7 && n % 11 && n % 13 && n % 17
     404    10221676 :       && n % 19 && n % 23 && n % 29 && n % 31 && n % 37 && n % 41
     405    28107079 :       && _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       93397 : uisprime_661(ulong n)
     420             : {
     421       93397 :   if (n <= maxprimelim()) return PRIMES_search(n) > 0;
     422       77260 :   if (n < 1016801) return n < 452929? 1: uispsp(2, n);
     423       77260 :   return _uisprime(n);
     424             : }
     425             : 
     426             : static int
     427      373680 : iu_coprime(GEN N, ulong u) { return ugcd(u, umodiu(N,u)) == 1; }
     428             : long
     429    48341267 : BPSW_psp(GEN N)
     430             : {
     431             :   pari_sp av;
     432    48341267 :   if (typ(N) != t_INT) pari_err_TYPE("BPSW_psp",N);
     433    48341765 :   if (signe(N) <= 0) return 0;
     434    48280011 :   if (lgefint(N) == 3) return uisprime(uel(N,2));
     435      198963 :   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      176882 :   if (!iu_coprime(N, 16294579238595022365UL) ||
     440      109588 :       !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      105509 :   av = avma;
     453      105509 :   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       45827 : Z_isanypower_nosmalldiv(GEN N, ulong cutoff, GEN *px)
     460             : {
     461       45827 :   GEN x = N, y;
     462       45827 :   ulong mask = 7;
     463       45827 :   long ex, k = 1;
     464             :   forprime_t T;
     465       57516 :   while (Z_issquareall(x, &y)) { k <<= 1; x = y; }
     466       45892 :   while ( (ex = is_357_power(x, &y, &mask)) ) { k *= ex; x = y; }
     467       45830 :   (void)u_forprime_init(&T, 11, ULONG_MAX);
     468             :   /* stop when x^(1/k) < 2^cutoff */
     469       45829 :   while ( (ex = is_pth_power(x, &y, &T, cutoff)) ) { k *= ex; x = y; }
     470       45837 :   *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       49980 : BPSW_psp_nosmalldiv(GEN N)
     477             : {
     478             :   pari_sp av;
     479       49980 :   long l = lgefint(N);
     480             : 
     481       49980 :   if (l == 3) return uisprime_661(uel(N,2));
     482       32333 :   av = avma;
     483             :   /* N large: test for pure power, rarely succeeds, but requires < 1% of
     484             :    * compositeness test times */
     485       32333 :   if (bit_accuracy(l) > 512 && Z_isanypower_nosmalldiv(N, 15, &N) != 1)
     486          14 :     return gc_long(av,0);
     487       32326 :   N = absi_shallow(N);
     488       32326 :   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     3950993 : BPSW_isprime_small(GEN x)
     502             : {
     503     3950993 :   long l = lgefint(x);
     504             : #ifdef LONG_IS_64BIT
     505     3883958 :   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     3949574 : BPSW_isprime(GEN N)
     683             : {
     684             :   pari_sp av;
     685             :   long t;
     686             :   GEN P;
     687     3949574 :   if (BPSW_isprime_small(N)) return 1;
     688        4503 :   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    37363064 : gisprime(GEN x, long flag)
     705             : {
     706    37363064 :   switch (flag)
     707             :   {
     708    37364040 :     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    38075205 : 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    57410797 : PRIMES_search(ulong x)
    1091             : {
    1092    57410797 :   pari_prime *T = pari_PRIMES;
    1093    57410797 :   ulong i, u = minuu(T[0], (x + 2) >> ((x < 122UL)? 1: 2)), l = 1;
    1094             :   do
    1095             :   {
    1096   426779797 :     i = (l+u) >> 1;
    1097   426779797 :     if (x < T[i]) u = i-1; else if (x > T[i]) l = i+1; else return i;
    1098   407833405 :   } while (u > l);
    1099    38455030 :   if (u == l) { i = l; if (x == T[i]) return i; }
    1100    29175551 :   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      402989 : primepi_upper_bound(double x)
    1154             : {
    1155      402989 :   if (x >= 355991)
    1156             :   {
    1157        1890 :     double L = 1/log(x);
    1158        1890 :     return x * L * (1 + L + 2.51*L*L);
    1159             :   }
    1160      401099 :   if (x >= 60184) return x / (log(x) - 1.1);
    1161      401099 :   if (x < 5) return 2; /* don't bother */
    1162      273770 :   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      324441 : pari_set_primetab(GEN global_primetab)
    1357             : {
    1358      324441 :   if (global_primetab)
    1359             :   {
    1360      322589 :     long i, l = lg(global_primetab);
    1361      322589 :     primetab = cgetg_block(l, t_VEC);
    1362      352969 :     for (i = 1; i < l; i++)
    1363       30499 :       gel(primetab,i) = gclone(gel(global_primetab,i));
    1364        1852 :   } else primetab = cgetg_block(1, t_VEC);
    1365      324321 : }
    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); }

Generated by: LCOV version 1.16