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

Generated by: LCOV version 1.13