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.10.0 lcov report (development 21059-cbe0d6a) Lines: 583 641 91.0 %
Date: 2017-09-22 06:24:58 Functions: 66 68 97.1 %
Legend: Lines: hit not hit

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

Generated by: LCOV version 1.11