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

Generated by: LCOV version 1.13