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-bordeaux1.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.8.0 lcov report (development 18590-b5f7c1c) Lines: 564 619 91.1 %
Date: 2016-02-09 Functions: 61 64 95.3 %
Legend: Lines: hit not hit | Branches: + taken - not taken # not executed Branches: 407 523 77.8 %

           Branch data     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                 :      67819 : init_MR_Jaeschke(MR_Jaeschke_t *S, GEN n)
      33                 :            : {
      34         [ -  + ]:      67819 :   if (signe(n) < 0) n = absi(n);
      35                 :      67819 :   S->n = n;
      36                 :      67819 :   S->t = addsi(-1,n);
      37                 :      67819 :   S->r1 = vali(S->t);
      38                 :      67819 :   S->t1 = shifti(S->t, -S->r1);
      39                 :      67819 :   S->sqrt1 = cgeti(lg(n)); S->sqrt1[1] = evalsigne(0)|evallgefint(2);
      40                 :      67819 :   S->sqrt2 = cgeti(lg(n)); S->sqrt2[1] = evalsigne(0)|evallgefint(2);
      41                 :      67819 : }
      42                 :            : static void
      43                 :    5273995 : Fl_init_MR_Jaeschke(Fl_MR_Jaeschke_t *S, ulong n)
      44                 :            : {
      45                 :    5273995 :   S->n = n;
      46                 :    5273995 :   S->t = n-1;
      47                 :    5273995 :   S->r1 = vals(S->t);
      48                 :    5278001 :   S->t1 = S->t >> S->r1;
      49                 :    5278001 :   S->sqrt1 = 0;
      50                 :    5278001 :   S->sqrt2 = 0;
      51                 :    5278001 : }
      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                 :       6039 : MR_Jaeschke_ok(MR_Jaeschke_t *S, GEN c)
      59                 :            : {
      60         [ +  + ]:       6039 :   if (signe(S->sqrt1))
      61                 :            :   { /* saw one earlier: compare */
      62 [ +  + ][ -  + ]:         61 :     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                 :       5978 :     affii(c, S->sqrt1);
      72                 :       5978 :     affii(subii(S->n, c), S->sqrt2);
      73                 :            :   }
      74                 :       6039 :   return 0;
      75                 :            : }
      76                 :            : static int
      77                 :     996552 : Fl_MR_Jaeschke_ok(Fl_MR_Jaeschke_t *S, ulong c)
      78                 :            : {
      79         [ +  + ]:     996552 :   if (S->sqrt1)
      80                 :            :   { /* saw one earlier: compare */
      81 [ +  + ][ -  + ]:        535 :     if (c != S->sqrt1 && c != S->sqrt2) return 1;
      82                 :            :   } else { /* remember */
      83                 :     996017 :     S->sqrt1 = c;
      84                 :     996017 :     S->sqrt2 = S->n - c;
      85                 :            :   }
      86                 :     996552 :   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                 :      67933 : bad_for_base(MR_Jaeschke_t *S, GEN a)
      93                 :            : {
      94                 :      67933 :   pari_sp av = avma;
      95                 :            :   long r;
      96                 :      67933 :   GEN c2, c = Fp_pow(a, S->t1, S->n);
      97                 :            : 
      98 [ +  + ][ +  + ]:      67933 :   if (is_pm1(c) || equalii(S->t, c)) return 0;
      99                 :            : 
     100                 :            :   /* go fishing for -1, not for 1 (saves one squaring) */
     101         [ +  + ]:     117203 :   for (r = S->r1 - 1; r; r--) /* r1 - 1 squarings */
     102                 :            :   {
     103                 :      63088 :     c2 = c; c = remii(sqri(c), S->n);
     104         [ +  + ]:      63088 :     if (equalii(S->t, c)) return MR_Jaeschke_ok(S, c2);
     105         [ -  + ]:      57049 :     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                 :      67933 :   return 1;
     112                 :            : }
     113                 :            : static int
     114                 :    5277412 : Fl_bad_for_base(Fl_MR_Jaeschke_t *S, ulong a)
     115                 :            : {
     116                 :            :   long r;
     117                 :    5277412 :   ulong c2, c = Fl_powu(a, S->t1, S->n);
     118                 :            : 
     119 [ +  + ][ +  + ]:    5286497 :   if (c == 1 || c == S->t) return 0;
     120                 :            : 
     121                 :            :   /* go fishing for -1, not for 1 (saves one squaring) */
     122         [ +  + ]:    5651254 :   for (r = S->r1 - 1; r; r--) /* r1 - 1 squarings */
     123                 :            :   {
     124                 :    3119045 :     c2 = c; c = Fl_sqr(c, S->n);
     125         [ +  + ]:    3119478 :     if (c == S->t) return Fl_MR_Jaeschke_ok(S, c2);
     126                 :            :   }
     127                 :    5286853 :   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                 :         28 :   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                 :          0 : ispseudoprime(GEN x, long flag)
     162         [ #  # ]:          0 : { 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                 :       1744 : Fl_MR_Jaeschke(ulong n, long k)
     195                 :            : {
     196                 :       1744 :   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         [ -  + ]:       1744 :   if (!(n & 1)) return 0;
     204         [ -  + ]:       1744 :   if (k == 16)
     205                 :            :   { /* use smaller (faster) bases if possible */
     206         [ #  # ]:          0 :     p = (n < 3215031751UL)? pr: pr+13;
     207                 :          0 :     k = 4;
     208                 :            :   }
     209         [ +  - ]:       1744 :   else if (k == 17)
     210                 :            :   {
     211         [ +  + ]:       1744 :     p = (n < 1373653UL)? pr: pr+11;
     212                 :       1744 :     k = 2;
     213                 :            :   }
     214                 :          0 :   else p = pr; /* 2,3,5,... */
     215                 :       1744 :   Fl_init_MR_Jaeschke(&S, n);
     216         [ +  + ]:       5136 :   for (i=1; i<=k; i++)
     217                 :            :   {
     218         [ -  + ]:       3440 :     r = p[i] % n; if (!r) break;
     219         [ +  + ]:       3440 :     if (Fl_bad_for_base(&S, r)) return 0;
     220                 :            :   }
     221                 :       1744 :   return 1;
     222                 :            : }
     223                 :            : 
     224                 :            : int
     225                 :       1873 : MR_Jaeschke(GEN n, long k)
     226                 :            : {
     227                 :       1873 :   pari_sp av2, av = avma;
     228                 :       1873 :   const ulong pr[] =
     229                 :            :     { 0, 2,3,5,7,11,13,17,19,23,29, 31,73, 2,13,23,1662803UL, };
     230                 :            :   const ulong *p;
     231                 :            :   long i;
     232                 :            :   MR_Jaeschke_t S;
     233                 :            : 
     234         [ +  + ]:       1873 :   if (lgefint(n) == 3) return Fl_MR_Jaeschke(uel(n,2), k);
     235                 :            : 
     236         [ -  + ]:        129 :   if (!mod2(n)) return 0;
     237         [ -  + ]:        129 :   if      (k == 16) { p = pr+13; k = 4; } /* 2,13,23,1662803 */
     238         [ +  - ]:        129 :   else if (k == 17) { p = pr+11; k = 2; } /* 31,73 */
     239                 :          0 :   else p = pr; /* 2,3,5,... */
     240                 :        129 :   init_MR_Jaeschke(&S, n); av2 = avma;
     241         [ +  + ]:        343 :   for (i=1; i<=k; i++)
     242                 :            :   {
     243         [ +  + ]:        236 :     if (bad_for_base(&S, utoipos(p[i]))) { avma = av; return 0; }
     244                 :        214 :     avma = av2;
     245                 :            :   }
     246                 :       1873 :   avma = av; return 1;
     247                 :            : }
     248                 :            : 
     249                 :            : /*********************************************************************/
     250                 :            : /**                                                                 **/
     251                 :            : /**                      PSEUDO PRIMALITY (LUCAS)                   **/
     252                 :            : /**                                                                 **/
     253                 :            : /*********************************************************************/
     254                 :            : /* compute n-th term of Lucas sequence modulo N.
     255                 :            :  * v_{k+2} = P v_{k+1} - v_k, v_0 = 2, v_1 = P.
     256                 :            :  * Assume n > 0 */
     257                 :            : static GEN
     258                 :      13590 : LucasMod(GEN n, ulong P, GEN N)
     259                 :            : {
     260                 :      13590 :   pari_sp av = avma;
     261                 :      13590 :   GEN nd = int_MSW(n);
     262                 :      13590 :   ulong m = *nd;
     263 [ +  + ][ +  + ]:      13590 :   long i, j = 1+bfffo(m);
         [ +  + ][ +  + ]
     264                 :      13590 :   GEN v = utoipos(P), v1 = utoipos(P*P - 2);
     265                 :            : 
     266                 :      13590 :   m <<= j; j = BITS_IN_LONG - j;
     267                 :      13590 :   for (i=lgefint(n)-2;;) /* cf. leftright_pow */
     268                 :            :   {
     269         [ +  + ]:     900973 :     for (; j; m<<=1,j--)
     270                 :            :     { /* v = v_k, v1 = v_{k+1} */
     271         [ +  + ]:     876442 :       if (m&HIGHBIT)
     272                 :            :       { /* set v = v_{2k+1}, v1 = v_{2k+2} */
     273                 :     149076 :         v = subis(mulii(v,v1), (long)P);
     274                 :     149073 :         v1= subis(sqri(v1), 2);
     275                 :            :       }
     276                 :            :       else
     277                 :            :       {/* set v = v_{2k}, v1 = v_{2k+1} */
     278                 :     727366 :         v1= subis(mulii(v,v1), (long)P);
     279                 :     727370 :         v = subis(sqri(v), 2);
     280                 :            :       }
     281                 :     876441 :       v = modii(v, N);
     282                 :     876440 :       v1= modii(v1,N);
     283         [ -  + ]:     876443 :       if (gc_needed(av,1))
     284                 :            :       {
     285         [ #  # ]:          0 :         if(DEBUGMEM>1) pari_warn(warnmem,"LucasMod");
     286                 :          0 :         gerepileall(av, 2, &v,&v1);
     287                 :            :       }
     288                 :            :     }
     289         [ +  + ]:      24531 :     if (--i == 0) return v;
     290                 :      10941 :     j = BITS_IN_LONG;
     291                 :      10941 :     nd=int_precW(nd);
     292                 :      10941 :     m = *nd;
     293                 :      10941 :   }
     294                 :            : }
     295                 :            : /* compute n-th term of Lucas sequence modulo N.
     296                 :            :  * v_{k+2} = P v_{k+1} - v_k, v_0 = 2, v_1 = P.
     297                 :            :  * Assume n > 0 */
     298                 :            : static ulong
     299                 :     942122 : u_LucasMod(ulong n, ulong P, ulong N)
     300                 :            : {
     301 [ +  + ][ +  + ]:     942122 :   long j = 1 + bfffo(n);
         [ +  + ][ +  + ]
     302                 :     942122 :   ulong v = P, v1 = P*P - 2, mP = N - P, m2 = N - 2, m = n << j;
     303                 :            : 
     304                 :     942122 :   j = BITS_IN_LONG - j;
     305 [ +  + ][ +  + ]:   50128546 :   for (; j; m<<=1,j--)
     306                 :            :   { /* v = v_k, v1 = v_{k+1} */
     307         [ +  + ]:   49186424 :     if (m & HIGHBIT)
     308                 :            :     { /* set v = v_{2k+1}, v1 = v_{2k+2} */
     309                 :    4017128 :       v = Fl_add(Fl_mul(v,v1,N), mP, N);
     310                 :    4017128 :       v1= Fl_add(Fl_mul(v1,v1,N),m2, N);
     311                 :            :     }
     312                 :            :     else
     313                 :            :     {/* set v = v_{2k}, v1 = v_{2k+1} */
     314                 :   45169296 :       v1= Fl_add(Fl_mul(v,v1,N),mP, N);
     315                 :   45169296 :       v = Fl_add(Fl_mul(v,v,N), m2, N);
     316                 :            :     }
     317                 :            :   }
     318                 :     942122 :   return v;
     319                 :            : }
     320                 :            : 
     321                 :            : int
     322                 :     942129 : uislucaspsp(ulong n)
     323                 :            : {
     324                 :            :   long i, v;
     325                 :     942129 :   ulong b, z, m2, m = n + 1;
     326                 :            : 
     327                 :     942129 :   for (b=3, i=0;; b+=2, i++)
     328                 :            :   {
     329                 :    1555779 :     ulong c = b*b - 4; /* = 1 mod 4 */
     330         [ +  + ]:    1555779 :     if (krouu(n % c, c) < 0) break;
     331 [ +  + ][ +  - ]:     613657 :     if (i == 64 && uissquareall(n, &c)) return 0; /* oo loop if N = m^2 */
     332                 :     613650 :   }
     333         [ -  + ]:     942122 :   if (!m) return 0; /* neither 2^32-1 nor 2^64-1 are Lucas-pp */
     334                 :     942122 :   v = vals(m); m >>= v;
     335                 :     942122 :   z = u_LucasMod(m, b, n);
     336         [ +  + ]:     942122 :   if (z == 2) return 1;
     337                 :     894314 :   m2 = n - 2;
     338         [ +  + ]:     894314 :   if (z == m2) return 1;
     339         [ +  + ]:     466307 :   for (i=1; i<v; i++)
     340                 :            :   {
     341         [ +  + ]:     466265 :     if (!z) return 1;
     342                 :     217396 :     z = Fl_add(Fl_mul(z,z, n), m2, n);
     343         [ -  + ]:     217396 :     if (z == 2) return 0;
     344                 :            :   }
     345                 :     942129 :   return 0;
     346                 :            : }
     347                 :            : /* N > 3. Caller should check that N is not a square first (taken care of here,
     348                 :            :  * but inefficient) */
     349                 :            : static int
     350                 :      13590 : IsLucasPsP(GEN N)
     351                 :            : {
     352                 :      13590 :   pari_sp av = avma;
     353                 :            :   GEN N_2, m, z;
     354                 :            :   long i, v;
     355                 :            :   ulong b;
     356                 :            : 
     357                 :      13590 :   for (b=3;; b+=2)
     358                 :            :   {
     359                 :      32563 :     ulong c = b*b - 4; /* = 1 mod 4 */
     360 [ -  + ][ #  # ]:      32563 :     if (b == 129 && Z_issquare(N)) return 0; /* avoid oo loop if N = m^2 */
     361         [ +  + ]:      32563 :     if (krouu(umodiu(N,c), c) < 0) break;
     362                 :      18973 :   }
     363                 :      13590 :   m = addis(N,1); v = vali(m); m = shifti(m,-v);
     364                 :      13590 :   z = LucasMod(m, b, N);
     365         [ +  + ]:      13590 :   if (equaliu(z, 2)) return 1;
     366                 :      11739 :   N_2 = subis(N,2);
     367         [ +  + ]:      11739 :   if (equalii(z, N_2)) return 1;
     368         [ +  + ]:      12768 :   for (i=1; i<v; i++)
     369                 :            :   {
     370         [ +  + ]:      12665 :     if (!signe(z)) return 1;
     371                 :       7005 :     z = modii(subis(sqri(z), 2), N);
     372         [ -  + ]:       7005 :     if (equaliu(z, 2)) return 0;
     373         [ -  + ]:       7005 :     if (gc_needed(av,1))
     374                 :            :     {
     375         [ #  # ]:          0 :       if(DEBUGMEM>1) pari_warn(warnmem,"IsLucasPsP");
     376                 :          0 :       z = gerepileupto(av, z);
     377                 :            :     }
     378                 :            :   }
     379                 :      13590 :   return 0;
     380                 :            : }
     381                 :            : 
     382                 :            : /* assume u odd, u > 1 */
     383                 :            : static int
     384                 :     230698 : iu_coprime(GEN N, ulong u)
     385                 :            : {
     386                 :     230698 :   const ulong n = umodiu(N, u);
     387 [ +  + ][ +  + ]:     230701 :   return (n == 1 || gcduodd(n, u) == 1);
     388                 :            : }
     389                 :            : /* assume u odd, u > 1 */
     390                 :            : static int
     391                 :   20821825 : uu_coprime(ulong n, ulong u)
     392                 :            : {
     393                 :   20821825 :   return gcduodd(n, u) == 1;
     394                 :            : }
     395                 :            : 
     396                 :            : /* composite strong 2-pseudoprime < 1016801 whose prime divisors are > 101 */
     397                 :            : static int
     398                 :    2107176 : is_2_prp_101(ulong n)
     399                 :            : {
     400         [ +  + ]:    2107176 :   switch(n) {
     401                 :            :   case 42799:
     402                 :            :   case 49141:
     403                 :            :   case 88357:
     404                 :            :   case 90751:
     405                 :            :   case 104653:
     406                 :            :   case 130561:
     407                 :            :   case 196093:
     408                 :            :   case 220729:
     409                 :            :   case 253241:
     410                 :            :   case 256999:
     411                 :            :   case 271951:
     412                 :            :   case 280601:
     413                 :            :   case 357761:
     414                 :            :   case 390937:
     415                 :            :   case 458989:
     416                 :            :   case 486737:
     417                 :            :   case 489997:
     418                 :            :   case 514447:
     419                 :            :   case 580337:
     420                 :            :   case 741751:
     421                 :            :   case 838861:
     422                 :            :   case 873181:
     423                 :            :   case 877099:
     424                 :            :   case 916327:
     425                 :            :   case 976873:
     426                 :        230 :   case 983401: return 1;
     427                 :    2107176 :   } return 0;
     428                 :            : }
     429                 :            : 
     430                 :            : static int
     431                 :    5272459 : u_2_prp(ulong n)
     432                 :            : {
     433                 :            :   Fl_MR_Jaeschke_t S;
     434                 :    5272459 :   Fl_init_MR_Jaeschke(&S, n);
     435                 :    5276396 :   return Fl_bad_for_base(&S, 2) == 0;
     436                 :            : }
     437                 :            : static int
     438 [ +  + ][ +  + ]:    3164584 : uBPSW_psp(ulong n) { return (u_2_prp(n) && uislucaspsp(n)); }
     439                 :            : 
     440                 :            : int
     441                 :   15609520 : uisprime(ulong n)
     442                 :            : {
     443         [ +  + ]:   15609520 :   if (n < 103)
     444         [ +  + ]:     849875 :     switch(n)
     445                 :            :     {
     446                 :            :       case 2:
     447                 :            :       case 3:
     448                 :            :       case 5:
     449                 :            :       case 7:
     450                 :            :       case 11:
     451                 :            :       case 13:
     452                 :            :       case 17:
     453                 :            :       case 19:
     454                 :            :       case 23:
     455                 :            :       case 29:
     456                 :            :       case 31:
     457                 :            :       case 37:
     458                 :            :       case 41:
     459                 :            :       case 43:
     460                 :            :       case 47:
     461                 :            :       case 53:
     462                 :            :       case 59:
     463                 :            :       case 61:
     464                 :            :       case 67:
     465                 :            :       case 71:
     466                 :            :       case 73:
     467                 :            :       case 79:
     468                 :            :       case 83:
     469                 :            :       case 89:
     470                 :            :       case 97:
     471                 :     648147 :       case 101: return 1;
     472                 :     201728 :       default: return 0;
     473                 :            :     }
     474         [ +  + ]:   14759645 :   if (!odd(n)) return 0;
     475                 :            : #ifdef LONG_IS_64BIT
     476                 :            :   /* 16294579238595022365 = 3*5*7*11*13*17*19*23*29*31*37*41*43*47*53
     477                 :            :    *  7145393598349078859 = 59*61*67*71*73*79*83*89*97*101 */
     478   [ +  +  +  + ]:   17362026 :   if (!uu_coprime(n, 16294579238595022365UL) ||
     479                 :   11990777 :       !uu_coprime(n,  7145393598349078859UL)) return 0;
     480                 :            : #else
     481                 :            :   /* 4127218095 = 3*5*7*11*13*17*19*23*37
     482                 :            :    * 3948078067 = 29*31*41*43*47*53
     483                 :            :    * 4269855901 = 59*83*89*97*101
     484                 :            :    * 1673450759 = 61*67*71*73*79 */
     485   [ +  +  +  + ]:    2101995 :   if (!uu_coprime(n, 4127218095UL) ||
     486         [ +  + ]:    1624697 :       !uu_coprime(n, 3948078067UL) ||
     487         [ +  + ]:    1470313 :       !uu_coprime(n, 1673450759UL) ||
     488                 :    1269922 :       !uu_coprime(n, 4269855901UL)) return 0;
     489                 :            : #endif
     490         [ +  + ]:    6242666 :   if (n < 10427) return 1;
     491 [ +  + ][ +  + ]:    5228912 :   if (n < 1016801) return !is_2_prp_101(n) && u_2_prp(n);
                 [ +  + ]
     492                 :   16236156 :   return uBPSW_psp(n);
     493                 :            : }
     494                 :            : 
     495                 :            : /* assume no prime divisor <= 101 */
     496                 :            : int
     497                 :      16481 : uisprime_101(ulong n)
     498                 :            : {
     499         [ +  + ]:      16481 :   if (n < 10427) return 1;
     500 [ +  + ][ +  - ]:      16474 :   if (n < 1016801) return !is_2_prp_101(n) && u_2_prp(n);
                 [ +  - ]
     501                 :      16481 :   return uBPSW_psp(n);
     502                 :            : }
     503                 :            : 
     504                 :            : /* assume no prime divisor <= 661 */
     505                 :            : int
     506                 :      26346 : uisprime_661(ulong n) { return uBPSW_psp(n); }
     507                 :            : 
     508                 :            : long
     509                 :    4647323 : BPSW_psp(GEN N)
     510                 :            : {
     511                 :            :   pari_sp av;
     512                 :            :   MR_Jaeschke_t S;
     513                 :            :   int k;
     514                 :            : 
     515         [ -  + ]:    4647323 :   if (typ(N) != t_INT) pari_err_TYPE("BPSW_psp",N);
     516         [ +  + ]:    4835874 :   if (signe(N) <= 0) return 0;
     517         [ +  + ]:    4841357 :   if (lgefint(N) == 3) return uisprime(uel(N,2));
     518         [ +  + ]:      93730 :   if (!mod2(N)) return 0;
     519                 :            : #ifdef LONG_IS_64BIT
     520                 :            :   /* 16294579238595022365 = 3*5*7*11*13*17*19*23*29*31*37*41*43*47*53
     521                 :            :    *  7145393598349078859 = 59*61*67*71*73*79*83*89*97*101 */
     522   [ +  +  +  + ]:      68543 :   if (!iu_coprime(N, 16294579238595022365UL) ||
     523                 :      36722 :       !iu_coprime(N,  7145393598349078859UL)) return 0;
     524                 :            : #else
     525                 :            :   /* 4127218095 = 3*5*7*11*13*17*19*23*37
     526                 :            :    * 3948078067 = 29*31*41*43*47*53
     527                 :            :    * 4269855901 = 59*83*89*97*101
     528                 :            :    * 1673450759 = 61*67*71*73*79 */
     529   [ +  +  +  + ]:      93909 :   if (!iu_coprime(N, 4127218095UL) ||
     530         [ +  + ]:      74932 :       !iu_coprime(N, 3948078067UL) ||
     531         [ +  + ]:      68248 :       !iu_coprime(N, 1673450759UL) ||
     532                 :      55620 :       !iu_coprime(N, 4269855901UL)) return 0;
     533                 :            : #endif
     534                 :            :   /* no prime divisor < 103 */
     535                 :      63427 :   av = avma;
     536                 :      63427 :   init_MR_Jaeschke(&S, N);
     537 [ +  + ][ +  + ]:      63427 :   k = (!bad_for_base(&S, gen_2) && IsLucasPsP(N));
     538                 :    4971562 :   avma = av; return k;
     539                 :            : }
     540                 :            : 
     541                 :            : /* can we write n = x^k ? Assume N has no prime divisor <= 2^14.
     542                 :            :  * Not memory clean */
     543                 :            : long
     544                 :       7683 : isanypower_nosmalldiv(GEN N, GEN *px)
     545                 :            : {
     546                 :       7683 :   GEN x = N, y;
     547                 :       7683 :   ulong mask = 7;
     548                 :       7683 :   long ex, k = 1;
     549                 :            :   forprime_t T;
     550         [ +  + ]:      11205 :   while (Z_issquareall(x, &y)) { k <<= 1; x = y; }
     551         [ +  + ]:       7711 :   while ( (ex = is_357_power(x, &y, &mask)) ) { k *= ex; x = y; }
     552                 :       7683 :   (void)u_forprime_init(&T, 11, ULONG_MAX);
     553                 :            :   /* stop when x^(1/k) < 2^14 */
     554         [ -  + ]:       7683 :   while ( (ex = is_pth_power(x, &y, &T, 15)) ) { k *= ex; x = y; }
     555                 :       7683 :   *px = x; return k;
     556                 :            : }
     557                 :            : 
     558                 :            : /* no prime divisor <= 2^14 (> 661) */
     559                 :            : long
     560                 :      14356 : BPSW_psp_nosmalldiv(GEN N)
     561                 :            : {
     562                 :            :   pari_sp av;
     563                 :            :   MR_Jaeschke_t S;
     564                 :      14356 :   long l = lgefint(N);
     565                 :            :   int k;
     566                 :            : 
     567         [ +  + ]:      14356 :   if (l == 3) return uisprime_661(uel(N,2));
     568                 :       4277 :   av = avma;
     569                 :            :   /* N large: test for pure power, rarely succeeds, but requires < 1% of
     570                 :            :    * compositeness test times */
     571 [ +  + ][ +  + ]:       4277 :   if (bit_accuracy(l) > 512 && isanypower_nosmalldiv(N, &N) != 1)
     572                 :            :   {
     573                 :         21 :     avma = av; return 0;
     574                 :            :   }
     575                 :       4256 :   init_MR_Jaeschke(&S, N);
     576 [ +  + ][ +  + ]:       4256 :   k = (!bad_for_base(&S, gen_2) && IsLucasPsP(N));
     577                 :      14356 :   avma = av; return k;
     578                 :            : }
     579                 :            : 
     580                 :            : /***********************************************************************/
     581                 :            : /**                                                                   **/
     582                 :            : /**                       Pocklington-Lehmer                          **/
     583                 :            : /**                        P-1 primality test                         **/
     584                 :            : /**                                                                   **/
     585                 :            : /***********************************************************************/
     586                 :            : /* Assume x BPSW pseudoprime. Check whether it's small enough to be certified
     587                 :            :  * prime (< 2^64). Reference for strong 2-pseudoprimes:
     588                 :            :  *   http://www.cecm.sfu.ca/Pseudoprimes/index-2-to-64.html */
     589                 :            : static int
     590                 :     951174 : BPSW_isprime_small(GEN x)
     591                 :            : {
     592                 :     951174 :   long l = lgefint(x);
     593                 :            : #ifdef LONG_IS_64BIT
     594                 :     885982 :   return (l == 3);
     595                 :            : #else
     596                 :      65192 :   return (l <= 4);
     597                 :            : #endif
     598                 :            : }
     599                 :            : 
     600                 :            : /* Assume N > 1, p^e || N-1, p prime. Find a witness a(p) such that
     601                 :            :  *   a^(N-1) = 1 (mod N)
     602                 :            :  *   a^(N-1)/p - 1 invertible mod N.
     603                 :            :  * Proves that any divisor of N is 1 mod p^e. Return 0 if N is composite */
     604                 :            : static ulong
     605                 :      15231 : pl831(GEN N, GEN p)
     606                 :            : {
     607                 :      15231 :   GEN b, c, g, Nmunp = diviiexact(addis(N,-1), p);
     608                 :      15231 :   pari_sp av = avma;
     609                 :            :   ulong a;
     610                 :      15231 :   for(a = 2;; a++, avma = av)
     611                 :            :   {
     612                 :      22328 :     b = Fp_pow(utoipos(a), Nmunp, N);
     613         [ +  + ]:      22328 :     if (!equali1(b)) break;
     614                 :       7097 :   }
     615                 :      15231 :   c = Fp_pow(b,p,N);
     616                 :      15231 :   g = gcdii(addis(b,-1), N); /* 0 < g < N */
     617 [ +  + ][ +  - ]:      15231 :   return (equali1(c) && equali1(g))? a: 0;
     618                 :            : }
     619                 :            : 
     620                 :            : /* Brillhart, Lehmer, Selfridge test (Crandall & Pomerance, Th 4.1.5)
     621                 :            :  * N^(1/3) <= f fully factored, f | N-1. If pl831(p) is true for
     622                 :            :  * any prime divisor p of f, then any divisor of N is 1 mod f.
     623                 :            :  * In that case return 1 iff N is prime */
     624                 :            : static int
     625                 :         63 : BLS_test(GEN N, GEN f)
     626                 :            : {
     627                 :            :   GEN c1, c2, r, q;
     628                 :         63 :   q = dvmdii(N, f, &r);
     629         [ -  + ]:         63 :   if (!is_pm1(r)) return 0;
     630                 :         63 :   c2 = dvmdii(q, f, &c1);
     631                 :            :   /* N = 1 + f c1 + f^2 c2, 0 <= c_i < f; check whether it is of the form
     632                 :            :    * (1 + fa)(1 + fb) */
     633                 :         63 :   return ! Z_issquare(subii(sqri(c1), shifti(c2,2)));
     634                 :            : }
     635                 :            : 
     636                 :            : /* BPSW_psp(N) && !BPSW_isprime_small(N). Decide between Pocklington-Lehmer
     637                 :            :  * and APRCL. Return a vector of (small) primes such that PL-witnesses
     638                 :            :  * guarantee the primality of N. Return NULL if PL is likely too expensive.
     639                 :            :  * Return gen_0 if BLS test finds N to be composite */
     640                 :            : static GEN
     641                 :       5008 : BPSW_try_PL(GEN N)
     642                 :            : {
     643                 :       5008 :   ulong B = minuu(1UL<<19, maxprime());
     644                 :       5008 :   GEN E, p, U, F, N_1 = subiu(N,1);
     645                 :       5008 :   GEN fa = Z_factor_limit(N_1, B), P = gel(fa,1);
     646                 :       5008 :   long n = lg(P)-1;
     647                 :            : 
     648                 :       5008 :   p = gel(P,n);
     649                 :            :   /* if p prime, then N-1 is fully factored */
     650 [ +  + ][ +  + ]:       5008 :   if (cmpii(p,sqru(B)) <= 0 || (BPSW_psp_nosmalldiv(p) && BPSW_isprime(p)))
                 [ +  - ]
     651                 :       3012 :     return P;
     652                 :            : 
     653                 :       1996 :   E = gel(fa,2);
     654                 :       1996 :   U = powii(p, gel(E,n)); /* unfactored part of N-1 */
     655                 :            :   /* factored part of N-1; n >= 2 since 2p | N-1 */
     656         [ +  + ]:       1996 :   F = (n == 2)? powii(gel(P,1), gel(E,1)): diviiexact(N_1,  U);
     657                 :       1996 :   setlg(P, n); /* remove last (composite) entry */
     658                 :            : 
     659                 :            :   /* N-1 = F U, F factored, U possibly composite, (U,F) = 1 */
     660         [ +  + ]:       1996 :   if (cmpii(F, U) > 0) return P; /* 1/2-smooth */
     661 [ +  + ][ +  - ]:       1989 :   if (cmpii(sqri(F), U) > 0) return BLS_test(N,F)? P: gen_0; /* 1/3-smooth */
     662                 :       5008 :   return NULL; /* not smooth enough */
     663                 :            : }
     664                 :            : 
     665                 :            : static GEN isprimePL(GEN N);
     666                 :            : static GEN PL_certificate(GEN N, GEN F);
     667                 :            : /* Assume N a BPSW pseudoprime. Return 0 if not prime, and a primality label
     668                 :            :  * otherwise: 1 (small), 2 (APRCL), or PL-certificate  */
     669                 :            : static GEN
     670                 :         49 : check_prime(GEN N)
     671                 :            : {
     672                 :            :   GEN P;
     673         [ +  - ]:         49 :   if (BPSW_isprime_small(N)) return gen_1;
     674                 :            :   /* PL for small N: APRCL is faster but we prefer a certificate */
     675         [ #  # ]:          0 :   if (expi(N) <= 250) return isprimePL(N);
     676                 :          0 :   P = BPSW_try_PL(N);
     677                 :            :   /* if PL likely too expensive: give up certificate and use APRCL */
     678 [ #  # ][ #  # ]:          0 :   if (!P) return isprimeAPRCL(N)? gen_2: gen_0;
     679         [ #  # ]:         49 :   return typ(P) != t_INT? PL_certificate(N,P): gen_0;
     680                 :            : }
     681                 :            : 
     682                 :            : /* F is a vector whose entries are primes. For each of them, find a PL
     683                 :            :  * witness. Return 0 if caller lied and F contains a composite */
     684                 :            : static long
     685                 :       3075 : PL_certify(GEN N, GEN F)
     686                 :            : {
     687                 :       3075 :   long i, l = lg(F);
     688         [ +  + ]:      18250 :   for(i = 1; i < l; i++)
     689         [ -  + ]:      15175 :     if (! pl831(N, gel(F,i))) return 0;
     690                 :       3075 :   return 1;
     691                 :            : }
     692                 :            : /* F is a vector whose entries are *believed* to be primes. For each of them,
     693                 :            :  * recording a witness and recursive primality certificate */
     694                 :            : static GEN
     695                 :         28 : PL_certificate(GEN N, GEN F)
     696                 :            : {
     697                 :         28 :   long i, l = lg(F);
     698                 :         28 :   GEN W = cgetg(l,t_COL);
     699                 :         28 :   GEN R = cgetg(l,t_COL);
     700         [ +  + ]:         77 :   for(i=1; i<l; i++)
     701                 :            :   {
     702                 :         56 :     GEN p = gel(F,i);
     703                 :         56 :     ulong witness = pl831(N,p);
     704         [ +  + ]:         56 :     if (!witness) return gen_0;
     705                 :         49 :     gel(W,i) = utoipos(witness);
     706                 :         49 :     gel(R,i) = check_prime(p);
     707         [ -  + ]:         49 :     if (isintzero(gel(R,i)))
     708                 :            :     { /* composite in prime factorisation ! */
     709                 :          0 :       err_printf("Not a prime: %Ps", p);
     710                 :          0 :       pari_err_BUG("PL_certificate [false prime number]");
     711                 :            :     }
     712                 :            :   }
     713                 :         28 :   return mkmat3(F, W, R);
     714                 :            : }
     715                 :            : /* Assume N is a strong BPSW pseudoprime, Pocklington-Lehmer primality proof.
     716                 :            :  * Return gen_0 (non-prime), gen_1 (small prime), matrix (large prime)
     717                 :            :  *
     718                 :            :  * The matrix has 3 columns, [a,b,c] with
     719                 :            :  * a[i] prime factor of N-1,
     720                 :            :  * b[i] witness for a[i] as in pl831
     721                 :            :  * c[i] check_prime(a[i]) */
     722                 :            : static GEN
     723                 :         35 : isprimePL(GEN N)
     724                 :            : {
     725                 :         35 :   pari_sp ltop = avma;
     726                 :            :   GEN cbrtN, N_1, F, f;
     727                 :            : 
     728         [ -  + ]:         35 :   if (typ(N) != t_INT) pari_err_TYPE("isprimePL",N);
     729      [ -  +  + ]:         35 :   switch(cmpis(N,2))
     730                 :            :   {
     731                 :          0 :     case -1:return gen_0;
     732                 :          7 :     case 0: return gen_1;
     733                 :            :   }
     734                 :            :   /* N > 2 */
     735                 :         28 :   cbrtN = sqrtnint(N, 3);
     736                 :         28 :   N_1 = addis(N,-1);
     737                 :         28 :   F = Z_factor_until(N_1, sqri(cbrtN));
     738                 :         28 :   f = factorback(F); /* factored part of N-1, f^3 > N */
     739         [ -  + ]:         28 :   if (DEBUGLEVEL>3)
     740                 :            :   {
     741                 :          0 :     GEN r = divri(itor(f,LOWDEFAULTPREC), N);
     742                 :          0 :     err_printf("Pocklington-Lehmer: proving primality of N = %Ps\n", N);
     743                 :          0 :     err_printf("Pocklington-Lehmer: N-1 factored up to %Ps! (%.3Ps%%)\n", f, r);
     744                 :            :   }
     745                 :            :   /* if N-1 is only N^(1/3)-smooth, BLS test */
     746 [ +  + ][ +  + ]:         28 :   if (!equalii(f,N_1) && cmpii(sqri(f),N) <= 0 && !BLS_test(N,f))
                 [ -  + ]
     747                 :          0 :   { avma = ltop; return gen_0; } /* Failed, N is composite */
     748                 :         35 :   return gerepilecopy(ltop, PL_certificate(N, gel(F,1)));
     749                 :            : }
     750                 :            : 
     751                 :            : /* assume N a BPSW pseudoprime, in particular, it is odd > 2. Prove N prime */
     752                 :            : long
     753                 :     951106 : BPSW_isprime(GEN N)
     754                 :            : {
     755                 :            :   pari_sp av;
     756                 :            :   long t;
     757                 :            :   GEN P;
     758         [ +  + ]:     951106 :   if (BPSW_isprime_small(N)) return 1;
     759                 :       5008 :   av = avma; P = BPSW_try_PL(N);
     760         [ +  + ]:       5008 :   if (!P)
     761                 :       1933 :     t = isprimeAPRCL(N); /* not smooth enough */
     762                 :            :   else
     763         [ +  - ]:       3075 :     t = (typ(P) == t_INT)? 0: PL_certify(N,P);
     764                 :     951119 :   avma = av; return t;
     765                 :            : }
     766                 :            : 
     767                 :            : GEN
     768                 :    3664422 : gisprime(GEN x, long flag)
     769                 :            : {
     770   [ +  +  +  - ]:    3664422 :   switch (flag)
     771                 :            :   {
     772                 :    3664387 :     case 0: return map_proto_lG(isprime,x);
     773                 :         21 :     case 1: return map_proto_G(isprimePL,x);
     774                 :         14 :     case 2: return map_proto_lG(isprimeAPRCL,x);
     775                 :            :   }
     776                 :          0 :   pari_err_FLAG("gisprime");
     777                 :    3909875 :   return NULL;
     778                 :            : }
     779                 :            : 
     780                 :            : long
     781 [ +  + ][ +  + ]:    4414316 : isprime(GEN x) { return BPSW_psp(x) && BPSW_isprime(x); }
     782                 :            : 
     783                 :            : /***********************************************************************/
     784                 :            : /**                                                                   **/
     785                 :            : /**                          PRIME NUMBERS                            **/
     786                 :            : /**                                                                   **/
     787                 :            : /***********************************************************************/
     788                 :            : 
     789                 :            : static struct {
     790                 :            :   ulong p;
     791                 :            :   long n;
     792                 :            : } prime_table[] = {
     793                 :            :   {           0,          0},
     794                 :            :   {        7919,       1000},
     795                 :            :   {       17389,       2000},
     796                 :            :   {       27449,       3000},
     797                 :            :   {       37813,       4000},
     798                 :            :   {       48611,       5000},
     799                 :            :   {       59359,       6000},
     800                 :            :   {       70657,       7000},
     801                 :            :   {       81799,       8000},
     802                 :            :   {       93179,       9000},
     803                 :            :   {      104729,      10000},
     804                 :            :   {      224737,      20000},
     805                 :            :   {      350377,      30000},
     806                 :            :   {      479909,      40000},
     807                 :            :   {      611953,      50000},
     808                 :            :   {      746773,      60000},
     809                 :            :   {      882377,      70000},
     810                 :            :   {     1020379,      80000},
     811                 :            :   {     1159523,      90000},
     812                 :            :   {     1299709,     100000},
     813                 :            :   {     2750159,     200000},
     814                 :            :   {     7368787,     500000},
     815                 :            :   {    15485863,    1000000},
     816                 :            :   {    32452843,    2000000},
     817                 :            :   {    86028121,    5000000},
     818                 :            :   {   179424673,   10000000},
     819                 :            :   {   373587883,   20000000},
     820                 :            :   {   982451653,   50000000},
     821                 :            :   {  2038074743,  100000000},
     822                 :            :   {  4000000483UL,189961831},
     823                 :            :   {  4222234741UL,200000000},
     824                 :            : #if BITS_IN_LONG == 64
     825                 :            :   { 11037271757,  500000000},
     826                 :            :   { 22801763489, 1000000000},
     827                 :            :   { 47055833459, 2000000000},
     828                 :            :   {122430513841, 5000000000},
     829                 :            :   {200000000507, 8007105083},
     830                 :            : #endif
     831                 :            : };
     832                 :            : static const int prime_table_len = numberof(prime_table);
     833                 :            : 
     834                 :            : /* find prime closest to n in prime_table. */
     835                 :            : static long
     836                 :    8667138 : prime_table_closest_p(ulong n)
     837                 :            : {
     838                 :            :   long i;
     839         [ +  - ]:   21686079 :   for (i = 1; i < prime_table_len; i++)
     840                 :            :   {
     841                 :   21686079 :     ulong p = prime_table[i].p;
     842         [ +  + ]:   21686079 :     if (p > n)
     843                 :            :     {
     844                 :    8667138 :       ulong u = n - prime_table[i-1].p;
     845         [ +  + ]:    8667138 :       if (p - n > u) i--;
     846                 :    8667138 :       break;
     847                 :            :     }
     848                 :            :   }
     849         [ -  + ]:    8667138 :   if (i == prime_table_len) i = prime_table_len - 1;
     850                 :    8667138 :   return i;
     851                 :            : }
     852                 :            : 
     853                 :            : /* return the n-th successor of prime p > 2 */
     854                 :            : static GEN
     855                 :         70 : prime_successor(ulong p, ulong n)
     856                 :            : {
     857                 :            :   forprime_t S;
     858                 :            :   ulong i;
     859                 :         70 :   forprime_init(&S, utoipos(p+1), NULL);
     860         [ +  + ]:    2402624 :   for (i = 1; i < n; i++) (void)forprime_next(&S);
     861                 :         70 :   return forprime_next(&S);
     862                 :            : }
     863                 :            : /* find the N-th prime */
     864                 :            : static GEN
     865                 :        161 : prime_table_find_n(ulong N)
     866                 :            : {
     867                 :            :   byteptr d;
     868                 :        161 :   ulong n, p, maxp = maxprime();
     869                 :            :   long i;
     870         [ +  - ]:       2037 :   for (i = 1; i < prime_table_len; i++)
     871                 :            :   {
     872                 :       2037 :     n = prime_table[i].n;
     873         [ +  + ]:       2037 :     if (n > N)
     874                 :            :     {
     875                 :        161 :       ulong u = N - prime_table[i-1].n;
     876         [ +  + ]:        161 :       if (n - N > u) i--;
     877                 :        161 :       break;
     878                 :            :     }
     879                 :            :   }
     880         [ -  + ]:        161 :   if (i == prime_table_len) i = prime_table_len - 1;
     881                 :        161 :   p = prime_table[i].p;
     882                 :        161 :   n = prime_table[i].n;
     883 [ +  + ][ +  + ]:        161 :   if (n > N && p > maxp)
     884                 :            :   {
     885                 :         14 :     i--;
     886                 :         14 :     p = prime_table[i].p;
     887                 :         14 :     n = prime_table[i].n;
     888                 :            :   }
     889                 :            :   /* if beyond prime table, then n <= N */
     890                 :        161 :   d = diffptr + n;
     891         [ +  + ]:        161 :   if (n > N)
     892                 :            :   {
     893                 :         14 :     n -= N;
     894         [ +  + ]:      50624 :     do { n--; PREC_PRIME_VIADIFF(p,d); } while (n) ;
     895                 :            :   }
     896         [ +  - ]:        147 :   else if (n < N)
     897                 :            :   {
     898                 :        147 :     n = N-n;
     899         [ +  + ]:        147 :     if (p > maxp) return prime_successor(p, n);
     900                 :            :     do {
     901         [ -  + ]:      44492 :       if (!*d) return prime_successor(p, n);
     902                 :      44492 :       n--; NEXT_PRIME_VIADIFF(p,d);
     903         [ +  + ]:      44492 :     } while (n) ;
     904                 :            :   }
     905                 :        161 :   return utoipos(p);
     906                 :            : }
     907                 :            : 
     908                 :            : ulong
     909                 :          0 : uprime(long N)
     910                 :            : {
     911                 :          0 :   pari_sp av = avma;
     912                 :            :   GEN p;
     913         [ #  # ]:          0 :   if (N <= 0) pari_err_DOMAIN("prime", "n", "<=",gen_0, stoi(N));
     914                 :          0 :   p = prime_table_find_n(N);
     915         [ #  # ]:          0 :   if (lgefint(p) != 3) pari_err_OVERFLOW("uprime");
     916                 :          0 :   avma = av; return p[2];
     917                 :            : }
     918                 :            : GEN
     919                 :        168 : prime(long N)
     920                 :            : {
     921                 :        168 :   pari_sp av = avma;
     922                 :            :   GEN p;
     923         [ +  + ]:        168 :   if (N <= 0) pari_err_DOMAIN("prime", "n", "<=",gen_0, stoi(N));
     924                 :        161 :   new_chunk(4); /*HACK*/
     925                 :        161 :   p = prime_table_find_n(N);
     926                 :        161 :   avma = av; return icopy(p);
     927                 :            : }
     928                 :            : 
     929                 :            : /* random b-bit prime */
     930                 :            : GEN
     931                 :         49 : randomprime(GEN N)
     932                 :            : {
     933                 :         49 :   pari_sp av = avma, av2;
     934                 :            :   GEN a, b, d;
     935         [ +  + ]:         49 :   if (!N)
     936                 :            :     for(;;)
     937                 :            :     {
     938                 :         63 :       ulong p = random_bits(31);
     939         [ +  + ]:         63 :       if (uisprime(p)) return utoipos(p);
     940                 :         56 :     }
     941      [ +  +  - ]:         42 :   switch(typ(N))
     942                 :            :   {
     943                 :            :     case t_INT:
     944                 :         14 :       a = gen_2;
     945                 :         14 :       b = subiu(N,1); /* between 2 and N-1 */
     946                 :         14 :       d = subiu(N,2);
     947         [ +  + ]:         14 :       if (signe(d) <= 0)
     948                 :          7 :         pari_err_DOMAIN("randomprime","N", "<=", gen_2, N);
     949                 :          7 :       break;
     950                 :            :     case t_VEC:
     951         [ -  + ]:         28 :       if (lg(N) != 3) pari_err_TYPE("randomprime",N);
     952                 :         28 :       a = gel(N,1);
     953                 :         28 :       b = gel(N,2);
     954         [ +  + ]:         28 :       if (gcmp(b, a) < 0)
     955                 :          7 :         pari_err_DOMAIN("randomprime","b-a", "<", gen_0, mkvec2(a,b));
     956         [ +  + ]:         21 :       if (typ(a) != t_INT)
     957                 :            :       {
     958                 :          7 :         a = gceil(a);
     959         [ -  + ]:          7 :         if (typ(a) != t_INT) pari_err_TYPE("randomprime",a);
     960                 :            :       }
     961         [ +  + ]:         21 :       if (typ(b) != t_INT)
     962                 :            :       {
     963                 :          7 :         b = gfloor(b);
     964         [ -  + ]:          7 :         if (typ(b) != t_INT) pari_err_TYPE("randomprime",b);
     965                 :            :       }
     966         [ +  + ]:         21 :       if (cmpis(a, 2) < 0)
     967                 :            :       {
     968                 :          7 :         a = gen_2;
     969                 :          7 :         d = subiu(b,1);
     970                 :            :       }
     971                 :            :       else
     972                 :         14 :         d = addiu(subii(b,a), 1);
     973         [ +  + ]:         21 :       if (signe(d) <= 0)
     974                 :         14 :         pari_err_DOMAIN("randomprime","floor(b) - max(ceil(a),2)", "<",
     975                 :            :                         gen_0, mkvec2(a,b));
     976                 :          7 :       break;
     977                 :            :     default:
     978                 :          0 :       pari_err_TYPE("randomprime", N);
     979                 :          0 :       return NULL; /*notreached*/
     980                 :            :   }
     981                 :         14 :   av2 = avma;
     982                 :            :   for (;;)
     983                 :            :   {
     984                 :        210 :     GEN p = addii(a, randomi(d));
     985         [ +  + ]:        210 :     if (BPSW_psp(p)) return gerepileuptoint(av, p);
     986                 :        196 :     avma = av2;
     987                 :        217 :   }
     988                 :            : }
     989                 :            : 
     990                 :            : /* set *pp = nextprime(a) = p
     991                 :            :  *     *pd so that NEXT_PRIME_VIADIFF(d, p) = nextprime(p+1)
     992                 :            :  *     *pn so that p = the n-th prime
     993                 :            :  * error if nextprime(a) is out of primetable bounds */
     994                 :            : void
     995                 :    8667019 : prime_table_next_p(ulong a, byteptr *pd, ulong *pp, ulong *pn)
     996                 :            : {
     997                 :            :   byteptr d;
     998                 :    8667019 :   ulong p, n, maxp = maxprime();
     999                 :    8667021 :   long i = prime_table_closest_p(a);
    1000                 :    8667021 :   p = prime_table[i].p;
    1001 [ +  + ][ -  + ]:    8667021 :   if (p > a && p > maxp)
    1002                 :            :   {
    1003                 :          0 :     i--;
    1004                 :          0 :     p = prime_table[i].p;
    1005                 :            :   }
    1006                 :            :   /* if beyond prime table, then p <= a */
    1007                 :    8667021 :   n = prime_table[i].n;
    1008                 :    8667021 :   d = diffptr + n;
    1009         [ +  + ]:    8667021 :   if (p < a)
    1010                 :            :   {
    1011         [ -  + ]:    4332117 :     if (a > maxp) pari_err_MAXPRIME(a);
    1012         [ +  + ]:   11618154 :     do { n++; NEXT_PRIME_VIADIFF(p,d); } while (p < a);
    1013                 :            :   }
    1014         [ +  + ]:    4334904 :   else if (p != a)
    1015                 :            :   {
    1016         [ +  + ]:    2152451 :     do { n--; PREC_PRIME_VIADIFF(p,d); } while (p > a) ;
    1017         [ +  + ]:       8218 :     if (p < a) { NEXT_PRIME_VIADIFF(p,d); n++; }
    1018                 :            :   }
    1019                 :    8667021 :   *pn = n;
    1020                 :    8667021 :   *pp = p;
    1021                 :    8667021 :   *pd = d;
    1022                 :    8667021 : }
    1023                 :            : 
    1024                 :            : ulong
    1025                 :       9246 : uprimepi(ulong a)
    1026                 :            : {
    1027                 :       9246 :   ulong p, n, maxp = maxprime();
    1028         [ +  + ]:       9246 :   if (a <= maxp)
    1029                 :            :   {
    1030                 :            :     byteptr d;
    1031                 :       9128 :     prime_table_next_p(a, &d, &p, &n);
    1032         [ +  + ]:       9128 :     return p == a? n: n-1;
    1033                 :            :   }
    1034                 :            :   else
    1035                 :            :   {
    1036                 :        118 :     long i = prime_table_closest_p(a);
    1037                 :            :     forprime_t S;
    1038                 :        118 :     p = prime_table[i].p;
    1039         [ +  + ]:        118 :     if (p > a)
    1040                 :            :     {
    1041                 :         28 :       i--;
    1042                 :         28 :       p = prime_table[i].p;
    1043                 :            :     }
    1044                 :            :     /* p = largest prime in table <= a */
    1045                 :        118 :     n = prime_table[i].n;
    1046                 :        118 :     (void)u_forprime_init(&S, p+1, a);
    1047         [ +  + ]:   52168632 :     for (; p; n++) p = u_forprime_next(&S);
    1048                 :       9246 :     return n-1;
    1049                 :            :   }
    1050                 :            : }
    1051                 :            : 
    1052                 :            : GEN
    1053                 :        245 : primepi(GEN x)
    1054                 :            : {
    1055                 :        245 :   pari_sp av = avma;
    1056         [ -  + ]:        245 :   GEN pp, nn, N = typ(x) == t_INT? x: gfloor(x);
    1057                 :            :   forprime_t S;
    1058                 :            :   ulong n, p;
    1059                 :            :   long i, l;
    1060         [ -  + ]:        245 :   if (typ(N) != t_INT) pari_err_TYPE("primepi",N);
    1061         [ -  + ]:        245 :   if (signe(N) <= 0) return gen_0;
    1062                 :        245 :   avma = av; l = lgefint(N);
    1063         [ +  + ]:        245 :   if (l == 3) return utoi(uprimepi(N[2]));
    1064                 :          1 :   i = prime_table_len-1;
    1065                 :          1 :   p = prime_table[i].p;
    1066                 :          1 :   n = prime_table[i].n;
    1067                 :          1 :   (void)forprime_init(&S, utoipos(p+1), N);
    1068                 :          1 :   nn = setloop(utoipos(n));
    1069                 :          1 :   pp = gen_0;
    1070         [ +  + ]:    3280223 :   for (; pp; incloop(nn)) pp = forprime_next(&S);
    1071                 :        245 :   return gerepileuptoint(av, subiu(nn,1));
    1072                 :            : }
    1073                 :            : 
    1074                 :            : /* pi(x) < x/log x * (1 + 1/log x + 2.51/log^2 x)), x>=355991 [ Dusart ]
    1075                 :            :  * pi(x) < x/(log x - 1.1), x >= 60184 [ Dusart ]
    1076                 :            :  * ? \p9
    1077                 :            :  * ? M = 0; for(x = 4, 60184, M = max(M, log(x) - x/primepi(x))); M
    1078                 :            :  * %1 = 1.11196252 */
    1079                 :            : double
    1080                 :      49506 : primepi_upper_bound(double x)
    1081                 :            : {
    1082         [ +  + ]:      49506 :   if (x >= 355991)
    1083                 :            :   {
    1084                 :         14 :     double L = 1/log(x);
    1085                 :         14 :     return x * L * (1 + L + 2.51*L*L);
    1086                 :            :   }
    1087         [ -  + ]:      49492 :   if (x >= 60184) return x / (log(x) - 1.1);
    1088         [ +  + ]:      49492 :   if (x < 5) return 2; /* don't bother */
    1089                 :      49506 :   return x / (log(x) - 1.111963);
    1090                 :            : }
    1091                 :            : /* pi(x) > x/log x (1 + 1/log x), x >= 599 [ Dusart ]
    1092                 :            :  * pi(x) > x / (log x + 2), x >= 55 [ Rosser ] */
    1093                 :            : double
    1094                 :         14 : primepi_lower_bound(double x)
    1095                 :            : {
    1096         [ +  - ]:         14 :   if (x >= 599)
    1097                 :            :   {
    1098                 :         14 :     double L = 1/log(x);
    1099                 :         14 :     return x * L * (1 + L);
    1100                 :            :   }
    1101         [ #  # ]:          0 :   if (x < 55) return 0; /* don't bother */
    1102                 :         14 :   return x / (log(x) + 2.);
    1103                 :            : }
    1104                 :            : GEN
    1105                 :          1 : gprimepi_upper_bound(GEN x)
    1106                 :            : {
    1107                 :          1 :   pari_sp av = avma;
    1108                 :            :   double L;
    1109         [ -  + ]:          1 :   if (typ(x) != t_INT) x = gfloor(x);
    1110         [ +  - ]:          1 :   if (expi(x) <= 1022)
    1111                 :            :   {
    1112                 :          1 :     avma = av;
    1113                 :          1 :     return dbltor(primepi_upper_bound(gtodouble(x)));
    1114                 :            :   }
    1115                 :          0 :   x = itor(x, LOWDEFAULTPREC);
    1116                 :          0 :   L = 1 / rtodbl(logr_abs(x));
    1117                 :          0 :   x = mulrr(x, dbltor(L * (1 + L + 2.51*L*L)));
    1118                 :          1 :   return gerepileuptoleaf(av, x);
    1119                 :            : }
    1120                 :            : GEN
    1121                 :          1 : gprimepi_lower_bound(GEN x)
    1122                 :            : {
    1123                 :          1 :   pari_sp av = avma;
    1124                 :            :   double L;
    1125         [ -  + ]:          1 :   if (typ(x) != t_INT) x = gfloor(x);
    1126         [ -  + ]:          1 :   if (cmpiu(x, 55) <= 0) return gen_0;
    1127         [ +  - ]:          1 :   if (expi(x) <= 1022)
    1128                 :            :   {
    1129                 :          1 :     avma = av;
    1130                 :          1 :     return dbltor(primepi_lower_bound(gtodouble(x)));
    1131                 :            :   }
    1132                 :          0 :   x = itor(x, LOWDEFAULTPREC);
    1133                 :          0 :   L = 1 / rtodbl(logr_abs(x));
    1134                 :          0 :   x = mulrr(x, dbltor(L * (1 + L)));
    1135                 :          1 :   return gerepileuptoleaf(av, x);
    1136                 :            : }
    1137                 :            : 
    1138                 :            : GEN
    1139                 :         63 : primes(long n)
    1140                 :            : {
    1141                 :            :   forprime_t S;
    1142                 :            :   long i;
    1143                 :            :   GEN y;
    1144         [ -  + ]:         63 :   if (n <= 0) return cgetg(1, t_VEC);
    1145                 :         63 :   y = cgetg(n+1, t_VEC);
    1146                 :         63 :   (void)new_chunk(3*n); /*HACK*/
    1147                 :         63 :   u_forprime_init(&S, 2, ULONG_MAX);
    1148                 :         63 :   avma = (pari_sp)y;
    1149         [ +  + ]:       3731 :   for (i = 1; i <= n; i++) gel(y, i) = utoipos( u_forprime_next(&S) );
    1150                 :         63 :   return y;
    1151                 :            : }
    1152                 :            : GEN
    1153                 :          0 : primes_zv(long n)
    1154                 :            : {
    1155                 :            :   forprime_t S;
    1156                 :            :   long i;
    1157                 :            :   GEN y;
    1158         [ #  # ]:          0 :   if (n <= 0) return cgetg(1, t_VECSMALL);
    1159                 :          0 :   y = cgetg(n+1,t_VECSMALL);
    1160                 :          0 :   u_forprime_init(&S, 2, ULONG_MAX);
    1161         [ #  # ]:          0 :   for (i = 1; i <= n; i++) y[i] =  u_forprime_next(&S);
    1162                 :          0 :   avma = (pari_sp)y; return y;
    1163                 :            : }
    1164                 :            : GEN
    1165                 :        119 : primes0(GEN N)
    1166                 :            : {
    1167      [ +  +  - ]:        119 :   switch(typ(N))
    1168                 :            :   {
    1169                 :         63 :     case t_INT: return primes(itos(N));
    1170                 :            :     case t_VEC:
    1171         [ +  - ]:         56 :       if (lg(N) == 3) return primes_interval(gel(N,1),gel(N,2));
    1172                 :            :   }
    1173                 :          0 :   pari_err_TYPE("primes", N);
    1174                 :        112 :   return NULL;
    1175                 :            : }
    1176                 :            : 
    1177                 :            : GEN
    1178                 :         56 : primes_interval(GEN a, GEN b)
    1179                 :            : {
    1180                 :         56 :   pari_sp av = avma;
    1181                 :            :   forprime_t S;
    1182                 :            :   long i, n;
    1183                 :            :   GEN y, d, p;
    1184         [ -  + ]:         56 :   if (typ(a) != t_INT)
    1185                 :            :   {
    1186                 :          0 :     a = gceil(a);
    1187         [ #  # ]:          0 :     if (typ(a) != t_INT) pari_err_TYPE("primes_interval",a);
    1188                 :            :   }
    1189         [ +  + ]:         56 :   if (typ(b) != t_INT)
    1190                 :            :   {
    1191                 :          7 :     b = gfloor(b);
    1192         [ +  - ]:          7 :     if (typ(b) != t_INT) pari_err_TYPE("primes_interval",b);
    1193                 :            :   }
    1194         [ +  + ]:         49 :   if (signe(a) < 0) a = gen_2;
    1195                 :         49 :   d = subii(b, a);
    1196 [ +  - ][ -  + ]:         49 :   if (signe(d) < 0 || signe(b) <= 0) { avma = av; return cgetg(1, t_VEC); }
    1197         [ +  + ]:         49 :   if (lgefint(b) == 3)
    1198                 :            :   {
    1199                 :         33 :     avma = av;
    1200                 :         33 :     y = primes_interval_zv(itou(a), itou(b));
    1201                 :         33 :     n = lg(y); settyp(y, t_VEC);
    1202         [ +  + ]:     468256 :     for (i = 1; i < n; i++) gel(y,i) = utoipos(y[i]);
    1203                 :         33 :     return y;
    1204                 :            :   }
    1205                 :            :   /* at most d+1 primes in [a,b]. If d large, try better bound to lower
    1206                 :            :    * memory use */
    1207         [ +  + ]:         16 :   if (cmpiu(d,100000) > 0)
    1208                 :            :   {
    1209                 :          1 :     GEN D = gsub(gprimepi_upper_bound(b), gprimepi_lower_bound(a));
    1210                 :          1 :     D = ceil_safe(D);
    1211         [ -  + ]:          1 :     if (cmpii(D, d) < 0) d = D;
    1212                 :            :   }
    1213                 :         16 :   n = itos(d)+1;
    1214                 :         16 :   forprime_init(&S, a, b);
    1215                 :         16 :   y = cgetg(n+1, t_VEC); i = 1;
    1216         [ +  + ]:       5854 :   while ((p = forprime_next(&S))) gel(y, i++) = icopy(p);
    1217                 :         49 :   setlg(y, i); return gerepileupto(av, y);
    1218                 :            : }
    1219                 :            : 
    1220                 :            : /* a <= b, at most d primes in [a,b]. Return them */
    1221                 :            : static GEN
    1222                 :        152 : primes_interval_i(ulong a, ulong b, ulong d)
    1223                 :            : {
    1224                 :        152 :   ulong p, i = 1, n = d + 1;
    1225                 :            :   forprime_t S;
    1226                 :        152 :   GEN y = cgetg(n+1, t_VECSMALL);
    1227                 :        152 :   pari_sp av = avma;
    1228                 :        152 :   u_forprime_init(&S, a, b);
    1229         [ +  + ]:     486169 :   while ((p = u_forprime_next(&S))) y[i++] = p;
    1230                 :        152 :   avma = av; setlg(y, i); stackdummy((pari_sp)(y + i), (pari_sp)(y + n+1));
    1231                 :        152 :   return y;
    1232                 :            : }
    1233                 :            : GEN
    1234                 :         33 : primes_interval_zv(ulong a, ulong b)
    1235                 :            : {
    1236                 :            :   ulong d;
    1237         [ -  + ]:         33 :   if (!a) return primes_upto_zv(b);
    1238         [ -  + ]:         33 :   if (b < a) return cgetg(1, t_VECSMALL);
    1239                 :         33 :   d = b - a;
    1240         [ +  + ]:         33 :   if (d > 100000UL)
    1241                 :            :   {
    1242                 :         13 :     ulong D = (ulong)ceil(primepi_upper_bound(b)-primepi_lower_bound(a));
    1243         [ +  + ]:         13 :     if (D < d) d = D;
    1244                 :            :   }
    1245                 :         33 :   return primes_interval_i(a, b, d);
    1246                 :            : }
    1247                 :            : GEN
    1248                 :        119 : primes_upto_zv(ulong b)
    1249                 :            : {
    1250                 :            :   ulong d;
    1251         [ -  + ]:        119 :   if (b < 2) return cgetg(1, t_VECSMALL);
    1252         [ -  + ]:        119 :   d = (b > 100000UL)? (ulong)primepi_upper_bound(b): b;
    1253                 :        119 :   return primes_interval_i(2, b, d);
    1254                 :            : }
    1255                 :            : 
    1256                 :            : /***********************************************************************/
    1257                 :            : /**                                                                   **/
    1258                 :            : /**                       PRIVATE PRIME TABLE                         **/
    1259                 :            : /**                                                                   **/
    1260                 :            : /***********************************************************************/
    1261                 :            : /* delete dummy NULL entries */
    1262                 :            : static void
    1263                 :         21 : cleanprimetab(GEN T)
    1264                 :            : {
    1265                 :         21 :   long i,j, l = lg(T);
    1266         [ +  + ]:         70 :   for (i = j = 1; i < l; i++)
    1267         [ +  + ]:         49 :     if (T[i]) T[j++] = T[i];
    1268                 :         21 :   setlg(T,j);
    1269                 :         21 : }
    1270                 :            : /* remove p from T */
    1271                 :            : static void
    1272                 :         28 : rmprime(GEN T, GEN p)
    1273                 :            : {
    1274                 :            :   long i;
    1275         [ -  + ]:         28 :   if (typ(p) != t_INT) pari_err_TYPE("removeprimes",p);
    1276                 :         28 :   i = ZV_search(T, p);
    1277         [ +  + ]:         28 :   if (!i)
    1278                 :          7 :     pari_err_DOMAIN("removeprime","prime","not in",
    1279                 :            :                     strtoGENstr("primetable"), p);
    1280                 :         21 :   gunclone(gel(T,i)); gel(T,i) = NULL;
    1281                 :         21 :   cleanprimetab(T);
    1282                 :         21 : }
    1283                 :            : 
    1284                 :            : /*stolen from ZV_union_shallow() : clone entries from y */
    1285                 :            : static GEN
    1286                 :         28 : addp_union(GEN x, GEN y)
    1287                 :            : {
    1288                 :         28 :   long i, j, k, lx = lg(x), ly = lg(y);
    1289                 :         28 :   GEN z = cgetg(lx + ly - 1, t_VEC);
    1290                 :         28 :   i = j = k = 1;
    1291 [ +  + ][ +  - ]:         35 :   while (i<lx && j<ly)
    1292                 :            :   {
    1293                 :          7 :     int s = cmpii(gel(x,i), gel(y,j));
    1294         [ -  + ]:          7 :     if (s < 0)
    1295                 :          0 :       gel(z,k++) = gel(x,i++);
    1296         [ -  + ]:          7 :     else if (s > 0)
    1297                 :          0 :       gel(z,k++) = gclone(gel(y,j++));
    1298                 :            :     else {
    1299                 :          7 :       gel(z,k++) = gel(x,i++);
    1300                 :          7 :       j++;
    1301                 :            :     }
    1302                 :            :   }
    1303         [ -  + ]:         28 :   while (i<lx) gel(z,k++) = gel(x,i++);
    1304         [ +  + ]:         77 :   while (j<ly) gel(z,k++) = gclone(gel(y,j++));
    1305                 :         28 :   setlg(z, k); return z;
    1306                 :            : }
    1307                 :            : 
    1308                 :            : /* p is NULL, or a single element or a row vector with "primes" to add to
    1309                 :            :  * prime table. */
    1310                 :            : static GEN
    1311                 :        161 : addp(GEN *T, GEN p)
    1312                 :            : {
    1313                 :        161 :   pari_sp av = avma;
    1314                 :            :   long i, l;
    1315                 :            :   GEN v;
    1316                 :            : 
    1317 [ +  + ][ -  + ]:        161 :   if (!p || lg(p) == 1) return *T;
    1318         [ +  + ]:         42 :   if (!is_vec_t(typ(p))) p = mkvec(p);
    1319                 :            : 
    1320                 :         42 :   RgV_check_ZV(p, "addprimes");
    1321                 :         35 :   v = gen_indexsort_uniq(p, (void*)&cmpii, &cmp_nodata);
    1322                 :         35 :   p = vecpermute(p, v);
    1323         [ +  + ]:         35 :   if (cmpiu(gel(p,1), 2) < 0) pari_err_DOMAIN("addprimes", "p", "<", gen_2,p);
    1324                 :         28 :   p = addp_union(*T, p);
    1325                 :         28 :   l = lg(p);
    1326         [ +  - ]:         28 :   if (l != lg(*T))
    1327                 :            :   {
    1328                 :         28 :     GEN old = *T, t = cgetalloc(t_VEC, l);
    1329         [ +  + ]:         84 :     for (i = 1; i < l; i++) gel(t,i) = gel(p,i);
    1330                 :         28 :     *T = t; free(old);
    1331                 :            :   }
    1332                 :        147 :   avma = av; return *T;
    1333                 :            : }
    1334                 :            : GEN
    1335                 :        161 : addprimes(GEN p) { return addp(&primetab, p); }
    1336                 :            : 
    1337                 :            : static GEN
    1338                 :         28 : rmprimes(GEN T, GEN prime)
    1339                 :            : {
    1340                 :            :   long i,tx;
    1341                 :            : 
    1342         [ -  + ]:         28 :   if (!prime) return T;
    1343                 :         28 :   tx = typ(prime);
    1344         [ +  + ]:         28 :   if (is_vec_t(tx))
    1345                 :            :   {
    1346         [ +  + ]:         14 :     if (prime == T)
    1347                 :            :     {
    1348         [ +  + ]:         14 :       for (i=1; i < lg(prime); i++) gunclone(gel(prime,i));
    1349                 :          7 :       setlg(prime, 1);
    1350                 :            :     }
    1351                 :            :     else
    1352                 :            :     {
    1353         [ +  + ]:         21 :       for (i=1; i < lg(prime); i++) rmprime(T, gel(prime,i));
    1354                 :            :     }
    1355                 :         14 :     return T;
    1356                 :            :   }
    1357                 :         28 :   rmprime(T, prime); return T;
    1358                 :            : }
    1359                 :            : GEN
    1360                 :         28 : removeprimes(GEN prime) { return rmprimes(primetab, prime); }

Generated by: LCOV version 1.9