Code coverage tests

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

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

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

LCOV - code coverage report
Current view: top level - basemath - aprcl.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.14.0 lcov report (development 27775-aca467eab2) Lines: 600 706 85.0 %
Date: 2022-07-03 07:33:15 Functions: 50 51 98.0 %
Legend: Lines: hit not hit

          Line data    Source code
       1             : /* Copyright (C) 2000  The PARI group.
       2             : 
       3             : This file is part of the PARI/GP package.
       4             : 
       5             : PARI/GP is free software; you can redistribute it and/or modify it under the
       6             : terms of the GNU General Public License as published by the Free Software
       7             : Foundation; either version 2 of the License, or (at your option) any later
       8             : version. It is distributed in the hope that it will be useful, but WITHOUT
       9             : ANY WARRANTY WHATSOEVER.
      10             : 
      11             : Check the License for details. You should have received a copy of it, along
      12             : with the package; see the file 'COPYING'. If not, write to the Free Software
      13             : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
      14             : 
      15             : /*******************************************************************/
      16             : /*                    THE APRCL PRIMALITY TEST                     */
      17             : /*******************************************************************/
      18             : #include "pari.h"
      19             : #include "paripriv.h"
      20             : 
      21             : #define DEBUGLEVEL DEBUGLEVEL_isprime
      22             : 
      23             : typedef struct Red {
      24             : /* global data */
      25             :   GEN N; /* prime we are certifying */
      26             :   GEN N2; /* floor(N/2) */
      27             : /* global data for flexible window */
      28             :   long k, lv;
      29             :   ulong mask;
      30             : /* reduction data */
      31             :   long n;
      32             :   GEN C; /* polcyclo(n) */
      33             :   GEN (*red)(GEN x, struct Red*);
      34             : } Red;
      35             : 
      36             : #define cache_aall(C)  (gel((C),1))
      37             : #define cache_tall(C)  (gel((C),2))
      38             : #define cache_cyc(C)  (gel((C),3))
      39             : #define cache_E(C)  (gel((C),4))
      40             : #define cache_eta(C)  (gel((C),5))
      41             : /* the last 4 only when N = 1 mod p^k */
      42             : #define cache_mat(C)  (gel((C),6))
      43             : #define cache_matinv(C)  (gel((C),7))
      44             : #define cache_a(C)  (gel((C),8)) /* a has order p^k in (Z/NZ)^* */
      45             : #define cache_pk(C)  (gel((C),9)) /* p^k */
      46             : 
      47             : static GEN
      48     1289444 : makepoldeg1(GEN c, GEN d)
      49             : {
      50             :   GEN z;
      51     1289444 :   if (signe(c)) {
      52     1289466 :     z = cgetg(4,t_POL); z[1] = evalsigne(1);
      53     1291031 :     gel(z,2) = d; gel(z,3) = c;
      54           0 :   } else if (signe(d)) {
      55           0 :     z = cgetg(3,t_POL); z[1] = evalsigne(1);
      56           0 :     gel(z,2) = d;
      57             :   } else {
      58           0 :     z = cgetg(2,t_POL); z[1] = evalsigne(0);
      59             :   }
      60     1291031 :   return z;
      61             : }
      62             : 
      63             : /* T mod polcyclo(p), assume deg(T) < 2p */
      64             : static GEN
      65     1800741 : red_cyclop(GEN T, long p)
      66             : {
      67             :   long i, d;
      68             :   GEN y, z;
      69             : 
      70     1800741 :   d = degpol(T) - p; /* < p */
      71     1800469 :   if (d <= -2) return T;
      72             : 
      73             :   /* reduce mod (x^p - 1) */
      74     1800469 :   y = ZX_mod_Xnm1(T, p);
      75     1798729 :   z = y+2;
      76             : 
      77             :   /* reduce mod x^(p-1) + ... + 1 */
      78     1798729 :   d = p-1;
      79     1798729 :   if (degpol(y) == d)
      80    11132899 :     for (i=0; i<d; i++) gel(z,i) = subii(gel(z,i), gel(z,d));
      81     1793490 :   return normalizepol_lg(y, d+2);
      82             : }
      83             : 
      84             : /* x t_VECSMALL, as red_cyclo2n_ip */
      85             : static GEN
      86        9537 : u_red_cyclo2n_ip(GEN x, long n)
      87             : {
      88        9537 :   long i, pow2 = 1L<<(n-1);
      89             :   GEN z;
      90             : 
      91       43883 :   for (i = lg(x)-1; i>pow2; i--) x[i-pow2] -= x[i];
      92        9841 :   for (; i>0; i--)
      93        9841 :     if (x[i]) break;
      94        9537 :   i += 2;
      95        9537 :   z = cgetg(i, t_POL); z[1] = evalsigne(1);
      96       43579 :   for (i--; i>=2; i--) gel(z,i) = stoi(x[i-1]);
      97        9537 :   return z;
      98             : }
      99             : /* x t_POL, n > 0. Return x mod polcyclo(2^n) = (x^(2^(n-1)) + 1). IN PLACE */
     100             : static GEN
     101      567983 : red_cyclo2n_ip(GEN x, long n)
     102             : {
     103      567983 :   long i, pow2 = 1L<<(n-1);
     104     5639898 :   for (i = lg(x)-1; i>pow2+1; i--)
     105     5073300 :     if (signe(gel(x,i))) gel(x,i-pow2) = subii(gel(x,i-pow2), gel(x,i));
     106      566598 :   return normalizepol_lg(x, i+1);
     107             : }
     108             : static GEN
     109      566563 : red_cyclo2n(GEN x, long n) { return red_cyclo2n_ip(leafcopy(x), n); }
     110             : 
     111             : /* special case R->C = polcyclo(2^n) */
     112             : static GEN
     113      566561 : _red_cyclo2n(GEN x, Red *R)
     114      566561 : { return centermod_i(red_cyclo2n(x, R->n), R->N, R->N2); }
     115             : /* special case R->C = polcyclo(p) */
     116             : static GEN
     117     1800812 : _red_cyclop(GEN x, Red *R)
     118     1800812 : { return centermod_i(red_cyclop(x, R->n), R->N, R->N2); }
     119             : static GEN
     120      808777 : _red(GEN x, Red *R)
     121      808777 : { return centermod_i(ZX_rem(x, R->C), R->N, R->N2); }
     122             : static GEN
     123    16700398 : modZ(GEN x, Red *R) { return centermodii(x, R->N, R->N2); }
     124             : static GEN
     125     7683866 : sqrmodZ(GEN x, Red *R) { return modZ(sqri(x), R); }
     126             : static GEN
     127     2246601 : sqrmod(GEN x, Red *R) { return R->red(ZX_sqr(x), R); }
     128             : static GEN
     129     2228629 : _mul(GEN x, GEN y, Red *R)
     130     1331405 : { return typ(x) == t_INT? modZ(mulii(x,y), R)
     131     3560259 :                         : R->red(ZX_mul(x,y), R); }
     132             : 
     133             : static GEN
     134           0 : sqrconst(GEN pol, Red *R)
     135             : {
     136           0 :   GEN z = cgetg(3,t_POL);
     137           0 :   gel(z,2) = modZ(sqri(gel(pol,2)), R);
     138           0 :   z[1] = pol[1]; return z;
     139             : }
     140             : 
     141             : /* pol^2 mod (x^2+x+1, N) */
     142             : static GEN
     143      756796 : sqrmod3(GEN pol, Red *R)
     144             : {
     145             :   GEN a,b,bma,A,B;
     146      756796 :   long lv = lg(pol);
     147             : 
     148      756796 :   if (lv==2) return pol;
     149      756796 :   if (lv==3) return sqrconst(pol, R);
     150      756796 :   a = gel(pol,3);
     151      756796 :   b = gel(pol,2); bma = subii(b,a);
     152      754995 :   A = modZ(mulii(a,addii(b,bma)), R);
     153      755857 :   B = modZ(mulii(bma,addii(a,b)), R); return makepoldeg1(A,B);
     154             : }
     155             : 
     156             : /* pol^2 mod (x^2+1,N) */
     157             : static GEN
     158      536436 : sqrmod4(GEN pol, Red *R)
     159             : {
     160             :   GEN a,b,A,B;
     161      536436 :   long lv = lg(pol);
     162             : 
     163      536436 :   if (lv==2) return pol;
     164      536436 :   if (lv==3) return sqrconst(pol, R);
     165      536436 :   a = gel(pol,3);
     166      536436 :   b = gel(pol,2);
     167      536436 :   A = modZ(mulii(a, shifti(b,1)), R);
     168      536091 :   B = modZ(mulii(subii(b,a),addii(b,a)), R); return makepoldeg1(A,B);
     169             : }
     170             : 
     171             : /* pol^2 mod (polcyclo(5),N) */
     172             : static GEN
     173     1313723 : sqrmod5(GEN pol, Red *R)
     174             : {
     175             :   GEN c2,b,c,d,A,B,C,D;
     176     1313723 :   long lv = lg(pol);
     177             : 
     178     1313723 :   if (lv==2) return pol;
     179     1313723 :   if (lv==3) return sqrconst(pol, R);
     180     1313723 :   c = gel(pol,3); c2 = shifti(c,1);
     181     1311894 :   d = gel(pol,2);
     182     1311894 :   if (lv==4)
     183             :   {
     184          42 :     A = modZ(sqri(d), R);
     185           0 :     B = modZ(mulii(c2, d), R);
     186           0 :     C = modZ(sqri(c), R); return mkpoln(3,A,B,C);
     187             :   }
     188     1311852 :   b = gel(pol,4);
     189     1311852 :   if (lv==5)
     190             :   {
     191          70 :     A = mulii(b, subii(c2,b));
     192          70 :     B = addii(sqri(c), mulii(b, subii(shifti(d,1),b)));
     193          70 :     C = subii(mulii(c2,d), sqri(b));
     194          70 :     D = mulii(subii(d,b), addii(d,b));
     195             :   }
     196             :   else
     197             :   { /* lv == 6 */
     198     1311782 :     GEN a = gel(pol,5), a2 = shifti(a,1);
     199             :     /* 2a(d - c) + b(2c - b) */
     200     1310913 :     A = addii(mulii(a2, subii(d,c)), mulii(b, subii(c2,b)));
     201             :     /* c(c - 2a) + b(2d - b) */
     202     1311243 :     B = addii(mulii(c, subii(c,a2)), mulii(b, subii(shifti(d,1),b)));
     203             :     /* (a-b)(a+b) + 2c(d - a) */
     204     1311138 :     C = addii(mulii(subii(a,b),addii(a,b)), mulii(c2,subii(d,a)));
     205             :     /* 2a(b - c) + (d-b)(d+b) */
     206     1311114 :     D = addii(mulii(a2, subii(b,c)), mulii(subii(d,b), addii(d,b)));
     207             :   }
     208     1311296 :   A = modZ(A, R);
     209     1312584 :   B = modZ(B, R);
     210     1312673 :   C = modZ(C, R);
     211     1312716 :   D = modZ(D, R); return mkpoln(4,A,B,C,D);
     212             : }
     213             : 
     214             : /* jac^floor(N/pk) mod (N, polcyclo(pk)), flexible window */
     215             : static GEN
     216       55523 : _powpolmod(GEN C, GEN jac, Red *R, GEN (*_sqr)(GEN, Red *))
     217             : {
     218       55523 :   const GEN taba = cache_aall(C);
     219       55523 :   const GEN tabt = cache_tall(C);
     220       55523 :   const long efin = lg(taba)-1, lv = R->lv;
     221       55523 :   GEN L, res = jac, pol2 = _sqr(res, R);
     222             :   long f;
     223       55517 :   pari_sp av0 = avma, av;
     224             : 
     225       55517 :   L = cgetg(lv+1, t_VEC); gel(L,1) = res;
     226      539565 :   for (f=2; f<=lv; f++) gel(L,f) = _mul(gel(L,f-1), pol2, R);
     227       55523 :   av = avma;
     228     1857828 :   for (f = efin; f >= 1; f--)
     229             :   {
     230     1802527 :     GEN t = gel(L, taba[f]);
     231     1802527 :     long tf = tabt[f];
     232     1802527 :     res = (f==efin)? t: _mul(t, res, R);
     233    14248228 :     while (tf--) {
     234    12445923 :       res = _sqr(res, R);
     235    12448256 :       if (gc_needed(av,1)) {
     236        1150 :         res = gerepilecopy(av, res);
     237        1150 :         if(DEBUGMEM>1) pari_warn(warnmem,"powpolmod: f = %ld",f);
     238             :       }
     239             :     }
     240             :   }
     241       55301 :   return gerepilecopy(av0, res);
     242             : }
     243             : 
     244             : static GEN
     245       12072 : _powpolmodsimple(GEN C, Red *R, GEN jac)
     246             : {
     247       12072 :   pari_sp av = avma;
     248       12072 :   GEN w = ZM_ZX_mul(cache_mat(C), jac);
     249       12070 :   long j, ph = lg(w);
     250             : 
     251       12070 :   R->red = &modZ;
     252       45973 :   for (j=1; j<ph; j++) gel(w,j) = _powpolmod(C, modZ(gel(w,j), R), R, &sqrmodZ);
     253       12072 :   w = centermod_i( ZM_ZC_mul(cache_matinv(C), w), R->N, R->N2);
     254       12072 :   w = gerepileupto(av, w);
     255       12072 :   return RgV_to_RgX(w, 0);
     256             : }
     257             : 
     258             : static GEN
     259       33696 : powpolmod(GEN C, Red *R, long p, long k, GEN jac)
     260             : {
     261             :   GEN (*_sqr)(GEN, Red *);
     262             : 
     263       33696 :   if (!isintzero(cache_mat(C))) return _powpolmodsimple(C, R, jac);
     264       21624 :   if (p == 2) /* p = 2 */
     265             :   {
     266        4277 :     if (k == 2) _sqr = &sqrmod4;
     267         790 :     else        _sqr = &sqrmod;
     268        4277 :     R->n = k;
     269        4277 :     R->red = &_red_cyclo2n;
     270       17347 :   } else if (k == 1)
     271             :   {
     272       16325 :     if      (p == 3) _sqr = &sqrmod3;
     273       11587 :     else if (p == 5) _sqr = &sqrmod5;
     274        4896 :     else             _sqr = &sqrmod;
     275       16325 :     R->n = p;
     276       16325 :     R->red = &_red_cyclop;
     277             :   } else {
     278        1022 :     R->red = &_red; _sqr = &sqrmod;
     279             :   }
     280       21624 :   return _powpolmod(C, jac, R, _sqr);
     281             : }
     282             : 
     283             : /* Return e(t) = \prod_{p-1 | t} p^{1+v_p(t)}}
     284             :  * faet contains the odd prime divisors of e(t) */
     285             : static GEN
     286        1672 : compute_e(ulong t, GEN *faet)
     287             : {
     288        1672 :   GEN L, P, D = divisorsu(t);
     289        1672 :   long l = lg(D);
     290             :   ulong k;
     291             : 
     292        1672 :   P = vecsmalltrunc_init(l);
     293        1672 :   L = vecsmalltrunc_init(l);
     294       35824 :   for (k=l-1; k>1; k--) /* k != 1: avoid d = 1 */
     295             :   {
     296       34152 :     ulong d = D[k];
     297       34152 :     if (uisprime(++d))
     298             :     { /* we want q = 1 (mod p) prime, not too large */
     299             : #ifdef LONG_IS_64BIT
     300       15808 :       if (d > 5000000000UL) return gen_0;
     301             : #endif
     302       18132 :       vecsmalltrunc_append(P, d);
     303       18132 :       vecsmalltrunc_append(L, upowuu(d, 1 + u_lval(t,d)));
     304             :     }
     305             :   }
     306        1672 :   if (faet) *faet = P;
     307        1672 :   return shifti(zv_prod_Z(L), 2 + u_lval(t,2));
     308             : }
     309             : 
     310             : /* Table obtained by the following script:
     311             : 
     312             : install(compute_e, "LD&"); \\ remove 'static' first
     313             : table(first = 6, step = 6, MAXT = 6983776800)=
     314             : {
     315             :   emax = 0;
     316             :   forstep(t = first, MAXT, step,
     317             :     e = compute_e(t);
     318             :     if (e > 1.9*emax, emax = e;
     319             :       printf("  if (C < %7.2f) return %8d;\n", 2*log(e)/log(2) - 1e-2, t)
     320             :     );
     321             :   );
     322             : }
     323             : table(,, 147026880);
     324             : table(147026880,5040, 6983776800);
     325             : */
     326             : 
     327             : /* assume C < 20003.8 */
     328             : static ulong
     329        1672 : compute_t_small(double C)
     330             : {
     331        1672 :   if (C <   17.94) return        6;
     332        1672 :   if (C <   31.99) return       12;
     333        1672 :   if (C <   33.99) return       24;
     334        1672 :   if (C <   54.07) return       36;
     335        1672 :   if (C <   65.32) return       60;
     336         972 :   if (C <   68.45) return       72;
     337         972 :   if (C <   70.78) return      108;
     338         972 :   if (C <   78.04) return      120;
     339         965 :   if (C <  102.41) return      180;
     340         965 :   if (C <  127.50) return      360;
     341         963 :   if (C <  136.68) return      420;
     342          53 :   if (C <  153.44) return      540;
     343          53 :   if (C <  165.67) return      840;
     344          53 :   if (C <  169.18) return     1008;
     345          53 :   if (C <  178.53) return     1080;
     346          53 :   if (C <  192.69) return     1200;
     347          53 :   if (C <  206.35) return     1260;
     348          53 :   if (C <  211.96) return     1620;
     349          53 :   if (C <  222.10) return     1680;
     350          53 :   if (C <  225.12) return     2016;
     351          53 :   if (C <  244.20) return     2160;
     352          53 :   if (C <  270.31) return     2520;
     353          53 :   if (C <  279.52) return     3360;
     354          53 :   if (C <  293.64) return     3780;
     355          53 :   if (C <  346.70) return     5040;
     356          53 :   if (C <  348.73) return     6480;
     357          53 :   if (C <  383.37) return     7560;
     358          53 :   if (C <  396.71) return     8400;
     359          53 :   if (C <  426.08) return    10080;
     360          53 :   if (C <  458.38) return    12600;
     361          53 :   if (C <  527.20) return    15120;
     362          53 :   if (C <  595.43) return    25200;
     363          53 :   if (C <  636.34) return    30240;
     364           7 :   if (C <  672.58) return    42840;
     365           7 :   if (C <  684.96) return    45360;
     366           7 :   if (C <  708.84) return    55440;
     367           7 :   if (C <  771.37) return    60480;
     368           7 :   if (C <  775.93) return    75600;
     369           7 :   if (C <  859.69) return    85680;
     370           7 :   if (C <  893.24) return   100800;
     371           7 :   if (C <  912.35) return   110880;
     372           7 :   if (C <  966.22) return   128520;
     373           7 :   if (C < 1009.18) return   131040;
     374           0 :   if (C < 1042.04) return   166320;
     375           0 :   if (C < 1124.98) return   196560;
     376           0 :   if (C < 1251.09) return   257040;
     377           0 :   if (C < 1375.13) return   332640;
     378           0 :   if (C < 1431.11) return   393120;
     379           0 :   if (C < 1483.46) return   514080;
     380           0 :   if (C < 1546.46) return   655200;
     381           0 :   if (C < 1585.94) return   665280;
     382           0 :   if (C < 1661.44) return   786240;
     383           0 :   if (C < 1667.67) return   831600;
     384           0 :   if (C < 1677.07) return   917280;
     385           0 :   if (C < 1728.17) return   982800;
     386           0 :   if (C < 1747.57) return  1081080;
     387           0 :   if (C < 1773.76) return  1179360;
     388           0 :   if (C < 1810.81) return  1285200;
     389           0 :   if (C < 1924.66) return  1310400;
     390           0 :   if (C < 2001.27) return  1441440;
     391           0 :   if (C < 2096.51) return  1663200;
     392           0 :   if (C < 2166.02) return  1965600;
     393           0 :   if (C < 2321.86) return  2162160;
     394           0 :   if (C < 2368.45) return  2751840;
     395           0 :   if (C < 2377.39) return  2827440;
     396           0 :   if (C < 2514.97) return  3326400;
     397           0 :   if (C < 2588.72) return  3341520;
     398           0 :   if (C < 2636.84) return  3603600;
     399           0 :   if (C < 2667.46) return  3931200;
     400           0 :   if (C < 3028.92) return  4324320;
     401           0 :   if (C < 3045.76) return  5654880;
     402           0 :   if (C < 3080.78) return  6652800;
     403           0 :   if (C < 3121.88) return  6683040;
     404           0 :   if (C < 3283.38) return  7207200;
     405           0 :   if (C < 3514.94) return  8648640;
     406           0 :   if (C < 3725.71) return 10810800;
     407           0 :   if (C < 3817.49) return 12972960;
     408           0 :   if (C < 3976.57) return 14414400;
     409           0 :   if (C < 3980.72) return 18378360;
     410           0 :   if (C < 4761.70) return 21621600;
     411           0 :   if (C < 5067.62) return 36756720;
     412           0 :   if (C < 5657.30) return 43243200;
     413           0 :   if (C < 5959.24) return 64864800;
     414           0 :   if (C < 6423.60) return 73513440;
     415           0 :   if (C < 6497.01) return 86486400;
     416           0 :   if (C < 6529.89) return 113097600;
     417           0 :   if (C < 6899.19) return 122522400;
     418           0 :   if (C < 7094.26) return 129729600;
     419           0 :   if (C < 7494.60) return 147026880;
     420           0 :   if (C < 7606.21) return 172972800;
     421           0 :   if (C < 7785.10) return 183783600;
     422           0 :   if (C < 7803.68) return 216216000;
     423           0 :   if (C < 8024.18) return 220540320;
     424           0 :   if (C < 8278.12) return 245044800;
     425           0 :   if (C < 8316.48) return 273873600;
     426           0 :   if (C < 8544.02) return 294053760;
     427           0 :   if (C < 8634.14) return 302702400;
     428           0 :   if (C < 9977.69) return 367567200;
     429           0 :   if (C < 10053.06) return 514594080;
     430           0 :   if (C < 10184.29) return 551350800;
     431           0 :   if (C < 11798.33) return 735134400;
     432           0 :   if (C < 11812.60) return 821620800;
     433           0 :   if (C < 11935.31) return 1029188160;
     434           0 :   if (C < 12017.99) return 1074427200;
     435           0 :   if (C < 12723.99) return 1102701600;
     436           0 :   if (C < 13702.71) return 1470268800;
     437           0 :   if (C < 13748.76) return 1643241600;
     438           0 :   if (C < 13977.37) return 2058376320;
     439           0 :   if (C < 14096.03) return 2148854400UL;
     440           0 :   if (C < 15082.25) return 2205403200UL;
     441           0 :   if (C < 15344.18) return 2572970400UL;
     442           0 :   if (C < 15718.37) return 2940537600UL;
     443           0 :   if (C < 15868.65) return 3491888400UL;
     444           0 :   if (C < 15919.88) return 3675672000UL;
     445           0 :   if (C < 16217.23) return 4108104000UL;
     446             : #ifdef LONG_IS_64BIT
     447           0 :   if (C < 17510.32) return 4410806400UL;
     448           0 :   if (C < 18312.87) return 5145940800UL;
     449           0 :   return 6983776800UL;
     450             : #else
     451           0 :   pari_err_IMPL("APRCL for large numbers on 32bit arch");
     452           0 :   return 0;
     453             : #endif
     454             : }
     455             : 
     456             : /* return t such that e(t) > sqrt(N), set *faet = odd prime divisors of e(t) */
     457             : static ulong
     458        1672 : compute_t(GEN N, GEN *e, GEN *faet)
     459             : { /* 2^e b <= N < 2^e (b+1), where b >= 2^52. Approximating log_2 N by
     460             :    * log2(gtodouble(N)) ~ e+log2(b), the error is less than log(1+1/b) < 1e-15*/
     461        1672 :   double C = dbllog2(N) + 1e-10; /* > log_2 N at least for N < 2^(2^21) */
     462             :   ulong t;
     463             :   /* Return "smallest" t such that f(t) >= C, which implies e(t) > sqrt(N) */
     464             :   /* For N < 2^20003.8 ~ 5.5 10^6021 */
     465        1672 :   if (C < 20003.8)
     466             :   {
     467        1672 :     t = compute_t_small(C);
     468        1672 :     *e = compute_e(t, faet);
     469             :   }
     470             :   else
     471             :   {
     472             : #ifdef LONG_IS_64BIT
     473           0 :     GEN B = sqrti(N);
     474           0 :     for (t = 6983776800UL+5040UL;; t+=5040)
     475           0 :     {
     476           0 :       pari_sp av = avma;
     477           0 :       *e = compute_e(t, faet);
     478           0 :       if (cmpii(*e, B) > 0) break;
     479           0 :       set_avma(av);
     480             :     }
     481             : #else
     482             :     *e = NULL; /* LCOV_EXCL_LINE */
     483             :     t = 0; /* LCOV_EXCL_LINE */
     484             : #endif
     485             :   }
     486        1672 :   return t;
     487             : }
     488             : 
     489             : /* T[i] = discrete log of i in (Z/q)^*, q odd prime
     490             :  * To save on memory, compute half the table: T[q-x] = T[x] + (q-1)/2 */
     491             : static GEN
     492       26392 : computetabdl(ulong q)
     493             : {
     494       26392 :   ulong g, a, i, qs2 = q>>1; /* (q-1)/2 */
     495       26392 :   GEN T = cgetg(qs2+2,t_VECSMALL);
     496             : 
     497       26392 :   g = pgener_Fl(q); a = 1;
     498     5965419 :   for (i=1; i < qs2; i++) /* g^((q-1)/2) = -1 */
     499             :   {
     500     5939027 :     a = Fl_mul(g,a,q);
     501     5939236 :     if (a > qs2) T[q-a] = i+qs2; else T[a] = i;
     502             :   }
     503       26392 :   T[qs2+1] = T[qs2] + qs2;
     504       26392 :   T[1] = 0; return T;
     505             : }
     506             : 
     507             : /* Return T: T[x] = dl of x(1-x) */
     508             : static GEN
     509       18139 : compute_g(ulong q)
     510             : {
     511       18139 :   const ulong qs2 = q>>1; /* (q-1)/2 */
     512             :   ulong x, a;
     513       18139 :   GEN T = computetabdl(q); /* updated in place to save on memory */
     514       18139 :   a = 0; /* dl[1] */
     515     3242779 :   for (x=2; x<=qs2+1; x++)
     516             :   { /* a = dl(x) */
     517     3224640 :     ulong b = T[x]; /* = dl(x) */
     518     3224640 :     T[x] = b + a + qs2; /* dl(x) + dl(x-1) + dl(-1) */
     519     3224640 :     a = b;
     520             :   }
     521       18139 :   return T;
     522             : }
     523             : 
     524             : /* x a nonzero VECSMALL with >= 0 entries */
     525             : static GEN
     526       25443 : zv_to_ZX(GEN x)
     527             : {
     528       25443 :   long i,j, lx = lg(x);
     529             :   GEN y;
     530             : 
     531       25505 :   while (lx-- && x[lx]==0) /* empty */;
     532       25443 :   i = lx+2; y = cgetg(i,t_POL); y[1] = evalsigne(1);
     533      157915 :   for (j=2; j<i; j++) gel(y,j) = utoi(x[j-1]);
     534       25440 :   return y;
     535             : }
     536             : /* p odd prime */
     537             : static GEN
     538       25442 : get_jac(GEN C, ulong q, long pk, GEN tabg)
     539             : {
     540       25442 :   ulong x, qs2 = q>>1; /* (q-1)/2 */
     541       25442 :   GEN vpk = zero_zv(pk);
     542             : 
     543    10019267 :   for (x=2; x<=qs2; x++) vpk[ tabg[x]%pk + 1 ] += 2;
     544       25442 :   vpk[ tabg[x]%pk + 1 ]++; /* x = (q+1)/2 */
     545       25442 :   return ZX_rem(zv_to_ZX(vpk), cache_cyc(C));
     546             : }
     547             : 
     548             : /* p = 2 */
     549             : static GEN
     550        8253 : get_jac2(GEN N, ulong q, long k, GEN *j2q, GEN *j3q)
     551             : {
     552        8253 :   GEN jpq, vpk, T = computetabdl(q);
     553             :   ulong x, pk, i, qs2;
     554             : 
     555             :   /* could store T[x+1] + T[x] + qs2 (cf compute_g).
     556             :    * Recompute instead, saving half the memory. */
     557        8253 :   pk = 1UL << k;;
     558        8253 :   vpk = zero_zv(pk);
     559             : 
     560        8253 :   qs2 = q>>1; /* (q-1)/2 */
     561             : 
     562     2935939 :   for (x=2; x<=qs2; x++) vpk[ (T[x]+T[x-1]+qs2)%pk + 1 ] += 2;
     563        8253 :   vpk[ (T[x]+T[x-1]+qs2)%pk + 1 ]++;
     564        8253 :   jpq = u_red_cyclo2n_ip(vpk, k);
     565             : 
     566        8253 :   if (k == 2) return jpq;
     567             : 
     568        1028 :   if (mod8(N) >= 5)
     569             :   {
     570         256 :     GEN v8 = cgetg(9,t_VECSMALL);
     571        2304 :     for (x=1; x<=8; x++) v8[x] = 0;
     572      464072 :     for (x=2; x<=qs2; x++) v8[ ((3*T[x]+T[x-1]+qs2)&7) + 1 ]++;
     573      464328 :     for (   ; x<=q-1; x++) v8[ ((3*T[q-x]+T[q-x+1]-3*qs2)&7) + 1 ]++;
     574         256 :     *j2q = RgX_inflate(ZX_sqr(u_red_cyclo2n_ip(v8,3)), pk>>3);
     575         256 :     *j2q = red_cyclo2n_ip(*j2q, k);
     576             :   }
     577       19900 :   for (i=1; i<=pk; i++) vpk[i] = 0;
     578     2545505 :   for (x=2; x<=qs2; x++) vpk[ (2*T[x]+T[x-1]+qs2)%pk + 1 ]++;
     579     2544598 :   for (   ; x<=q-1; x++) vpk[ (2*T[q-x]+T[q-x+1]-2*qs2)%pk + 1 ]++;
     580        1028 :   *j3q = ZX_mul(jpq, u_red_cyclo2n_ip(vpk,k));
     581        1028 :   *j3q = red_cyclo2n_ip(*j3q, k);
     582        1028 :   return jpq;
     583             : }
     584             : 
     585             : /* N = 1 mod p^k, return an elt of order p^k in (Z/N)^* */
     586             : static GEN
     587        2321 : finda(GEN Cp, GEN N, long pk, long p)
     588             : {
     589             :   GEN a, pv;
     590        2321 :   if (Cp && !isintzero(cache_a(Cp))) {
     591          35 :     a  = cache_a(Cp);
     592          35 :     pv = cache_pk(Cp);
     593             :   }
     594             :   else
     595             :   {
     596             :     GEN ph, b, q;
     597        2286 :     ulong u = 2;
     598        2286 :     long v = Z_lvalrem(subiu(N,1), p, &q);
     599        2286 :     ph = powuu(p, v-1); pv = muliu(ph, p); /* N - 1 = p^v q */
     600        2286 :     if (p > 2)
     601             :     {
     602         819 :       for (;;u++)
     603             :       {
     604        2230 :         a = Fp_pow(utoipos(u), q, N);
     605        2230 :         b = Fp_pow(a, ph, N);
     606        2230 :         if (!gequal1(b)) break;
     607             :       }
     608             :     }
     609             :     else
     610             :     {
     611        3136 :       while (krosi(u,N) >= 0) u++;
     612         875 :       a = Fp_pow(utoipos(u), q, N);
     613         875 :       b = Fp_pow(a, ph, N);
     614             :     }
     615             :     /* checking b^p = 1 mod N done economically in caller */
     616        2286 :     b = gcdii(subiu(b,1), N);
     617        2286 :     if (!gequal1(b)) return NULL;
     618             : 
     619        2286 :     if (Cp) {
     620        2286 :       cache_a(Cp)  = a; /* a has order p^v */
     621        2286 :       cache_pk(Cp) = pv;
     622             :     }
     623             :   }
     624        2321 :   return Fp_pow(a, divis(pv, pk), N);
     625             : }
     626             : 
     627             : /* return 0: N not a prime, 1: no problem so far */
     628             : static long
     629        7934 : filltabs(GEN C, GEN Cp, Red *R, long p, long pk, long ltab)
     630             : {
     631             :   pari_sp av;
     632             :   long i, j;
     633             :   long e;
     634             :   GEN tabt, taba, m;
     635             : 
     636        7934 :   cache_cyc(C) = polcyclo(pk,0);
     637             : 
     638        7934 :   if (p > 2)
     639             :   {
     640        4422 :     long LE = pk - pk/p + 1;
     641        4422 :     GEN E = cgetg(LE, t_VECSMALL), eta = cgetg(pk+1,t_VEC);
     642       22008 :     for (i=1,j=0; i<pk; i++)
     643       17586 :       if (i%p) E[++j] = i;
     644        4422 :     cache_E(C) = E;
     645             : 
     646       26430 :     for (i=1; i<=pk; i++)
     647             :     {
     648       22008 :       GEN z = FpX_rem(pol_xn(i-1, 0), cache_cyc(C), R->N);
     649       22008 :       gel(eta,i) = FpX_center_i(z, R->N, R->N2);
     650             :     }
     651        4422 :     cache_eta(C) = eta;
     652             :   }
     653        3512 :   else if (pk >= 8)
     654             :   {
     655         168 :     long LE = (pk>>2) + 1;
     656         168 :     GEN E = cgetg(LE, t_VECSMALL);
     657        3040 :     for (i=1,j=0; i<pk; i++)
     658        2872 :       if ((i%8)==1 || (i%8)==3) E[++j] = i;
     659         168 :     cache_E(C) = E;
     660             :   }
     661             : 
     662        7934 :   if (pk > 2 && umodiu(R->N,pk) == 1)
     663             :   {
     664        2321 :     GEN vpa, p1, p2, p3, a2 = NULL, a = finda(Cp, R->N, pk, p);
     665             :     long jj, ph;
     666             : 
     667        2321 :     if (!a) return 0;
     668        2321 :     ph = pk - pk/p;
     669        2321 :     vpa = cgetg(ph+1,t_COL); gel(vpa,1) = a;
     670        2321 :     if (pk > p) a2 = modZ(sqri(a), R);
     671        2321 :     jj = 1;
     672        7180 :     for (i=2; i<pk; i++) /* vpa = { a^i, (i,p) = 1 } */
     673        4859 :       if (i%p) {
     674        3823 :         GEN z = mulii((i%p==1) ? a2 : a, gel(vpa,jj));
     675        3823 :         gel(vpa,++jj) = modZ(z , R);
     676             :       }
     677        2321 :     if (!gequal1( modZ( mulii(a, gel(vpa,ph)), R) )) return 0;
     678        2321 :     p1 = cgetg(ph+1,t_MAT);
     679        2321 :     p2 = cgetg(ph+1,t_COL); gel(p1,1) = p2;
     680        8465 :     for (i=1; i<=ph; i++) gel(p2,i) = gen_1;
     681        2321 :     j = 2; gel(p1,j) = vpa; p3 = vpa;
     682        3823 :     for (j++; j <= ph; j++)
     683             :     {
     684        1502 :       p2 = cgetg(ph+1,t_COL); gel(p1,j) = p2;
     685        9302 :       for (i=1; i<=ph; i++)
     686        7800 :         gel(p2,i) = modZ(mulii(gel(vpa,i),gel(p3,i)), R);
     687        1502 :       p3 = p2;
     688             :     }
     689        2321 :     cache_mat(C) = p1;
     690        2321 :     cache_matinv(C) = FpM_inv(p1, R->N);
     691             :   }
     692             : 
     693        7934 :   tabt = cgetg(ltab+1, t_VECSMALL);
     694        7934 :   taba = cgetg(ltab+1, t_VECSMALL);
     695        7934 :   av = avma; m = divis(R->N, pk);
     696      149576 :   for (e=1; e<=ltab && signe(m); e++)
     697             :   {
     698      141642 :     long s = vali(m); m = shifti(m,-s);
     699      141642 :     tabt[e] = e==1? s: s + R->k;
     700      141642 :     taba[e] = signe(m)? ((mod2BIL(m) & R->mask)+1)>>1: 0;
     701      141642 :     m = shifti(m, -R->k);
     702             :   }
     703        7934 :   setlg(taba, e); cache_aall(C) = taba;
     704        7934 :   setlg(tabt, e); cache_tall(C) = tabt;
     705        7934 :   return gc_long(av,1);
     706             : }
     707             : 
     708             : static GEN
     709        1672 : calcglobs(Red *R, ulong t, long *plpC, long *pltab, GEN *pP)
     710             : {
     711             :   GEN fat, P, E, PE;
     712             :   long lv, i, k, b;
     713             :   GEN pC;
     714             : 
     715        1672 :   b = expi(R->N)+1;
     716        2697 :   k = 3; while (((k+1)*(k+2) << (k-1)) < b) k++;
     717        1672 :   *pltab = (b/k)+2;
     718        1672 :   R->k  = k;
     719        1672 :   R->lv = 1L << (k-1);
     720        1672 :   R->mask = (1UL << k) - 1;
     721             : 
     722        1672 :   fat = factoru_pow(t);
     723        1672 :   P = gel(fat,1);
     724        1672 :   E = gel(fat,2);
     725        1672 :   PE= gel(fat,3);
     726        1672 :   *plpC = lv = vecsmall_max(PE); /* max(p^e, p^e | t) */
     727        1672 :   pC = zerovec(lv);
     728        1672 :   gel(pC,1) = zerovec(9); /* to be used as temp in step5() */
     729       11640 :   for (i = 2; i <= lv; i++) gel(pC,i) = gen_0;
     730        7658 :   for (i=1; i<lg(P); i++)
     731             :   {
     732        5986 :     long pk, p = P[i], e = E[i];
     733        5986 :     pk = p;
     734       13913 :     for (k=1; k<=e; k++, pk*=p)
     735             :     {
     736        7927 :       gel(pC,pk) = zerovec(9);
     737        7927 :       if (!filltabs(gel(pC,pk), gel(pC,p), R, p,pk, *pltab)) return NULL;
     738             :     }
     739             :   }
     740        1672 :   *pP = P; return pC;
     741             : }
     742             : 
     743             : /* sig_a^{-1}(z) for z in Q(zeta_pk) and sig_a: zeta -> zeta^a. Assume
     744             :  * a reduced mod pk := p^k*/
     745             : static GEN
     746      163351 : aut(long pk, GEN z, long a)
     747             : {
     748             :   GEN v;
     749      163351 :   long b, i, dz = degpol(z);
     750      163352 :   if (a == 1 || dz < 0) return z;
     751      136889 :   v = cgetg(pk+2,t_POL); v[1] = evalvarn(0); b = 0;
     752      136886 :   gel(v,2) = gel(z,2); /* i = 0 */
     753     1236111 :   for (i = 1; i < pk; i++)
     754             :   {
     755     1099225 :     b += a; if (b > pk) b -= pk; /* b = (a*i) % pk */
     756     1099225 :     gel(v,i+2) = b > dz? gen_0: gel(z,b+2);
     757             :   }
     758      136886 :   return normalizepol_lg(v, pk+2);
     759             : }
     760             : 
     761             : /* z^v for v in Z[G], represented by couples [sig_x^{-1},x] */
     762             : static GEN
     763       26470 : autvec_TH(long pk, GEN z, GEN v, GEN C)
     764             : {
     765       26470 :   pari_sp av = avma;
     766       26470 :   long i, lv = lg(v);
     767       26470 :   GEN s = pol_1(varn(C));
     768      133477 :   for (i=1; i<lv; i++)
     769             :   {
     770      107009 :     long y = v[i];
     771      107009 :     if (y) s = ZXQ_mul(s, ZXQ_powu(aut(pk,z, y), y, C), C);
     772             :   }
     773       26468 :   return gerepileupto(av,s);
     774             : }
     775             : 
     776             : static GEN
     777       26471 : autvec_AL(long pk, GEN z, GEN v, Red *R)
     778             : {
     779       26471 :   pari_sp av = avma;
     780       26471 :   const long r = umodiu(R->N, pk);
     781       26471 :   GEN s = pol_1(varn(R->C));
     782       26473 :   long i, lv = lg(v);
     783      133487 :   for (i=1; i<lv; i++)
     784             :   {
     785      107016 :     long y = (r*v[i]) / pk;
     786      107016 :     if (y) s = RgXQ_mul(s, ZXQ_powu(aut(pk,z, v[i]), y, R->C), R->C);
     787             :   }
     788       26471 :   return gerepileupto(av,s);
     789             : }
     790             : 
     791             : /* 0 <= i < pk, such that x^i = z mod polcyclo(pk),  -1 if no such i exist */
     792             : static long
     793       25442 : look_eta(GEN eta, long pk, GEN z)
     794             : {
     795             :   long i;
     796       78713 :   for (i=1; i<=pk; i++)
     797       78715 :     if (ZX_equal(z, gel(eta,i))) return i-1;
     798           0 :   return -1;
     799             : }
     800             : /* same pk = 2^k */
     801             : static long
     802        8252 : look_eta2(long k, GEN z)
     803             : {
     804             :   long d, s;
     805        8252 :   if (typ(z) != t_POL) d = 0; /* t_INT */
     806             :   else
     807             :   {
     808        8252 :     if (!RgX_is_monomial(z)) return -1;
     809        8252 :     d = degpol(z);
     810        8252 :     z = gel(z,d+2); /* leading term */
     811             :   }
     812        8252 :   s = signe(z);
     813        8252 :   if (!s || !is_pm1(z)) return -1;
     814        8252 :   return (s > 0)? d: d + (1L<<(k-1));
     815             : }
     816             : 
     817             : static long
     818       25441 : step4a(GEN C, Red *R, ulong q, long p, long k, GEN tabg)
     819             : {
     820       25441 :   const long pk = upowuu(p,k);
     821             :   long ind;
     822             :   GEN jpq, s1, s2, s3;
     823             : 
     824       25441 :   if (!tabg) tabg = compute_g(q);
     825       25441 :   jpq = get_jac(C, q, pk, tabg);
     826       25442 :   s1 = autvec_TH(pk, jpq, cache_E(C), cache_cyc(C));
     827       25443 :   s2 = powpolmod(C,R, p,k, s1);
     828       25443 :   s3 = autvec_AL(pk, jpq, cache_E(C), R);
     829       25443 :   s3 = _red(gmul(s3,s2), R);
     830             : 
     831       25442 :   ind = look_eta(cache_eta(C), pk, s3);
     832       25442 :   if (ind < 0) return -1;
     833       25442 :   return (ind%p) != 0;
     834             : }
     835             : 
     836             : /* x == -1 mod N ? */
     837             : static long
     838        8700 : is_m1(GEN x, GEN N)
     839             : {
     840        8700 :   return equalii(addiu(x,1), N);
     841             : }
     842             : 
     843             : /* p=2, k>=3 */
     844             : static long
     845        1028 : step4b(GEN C, Red *R, ulong q, long k)
     846             : {
     847        1028 :   const long pk = 1L << k;
     848             :   long ind;
     849        1028 :   GEN s1, s2, s3, j2q = NULL, j3q = NULL;
     850             : 
     851        1028 :   (void)get_jac2(R->N,q,k, &j2q,&j3q);
     852             : 
     853        1028 :   s1 = autvec_TH(pk, j3q, cache_E(C), cache_cyc(C));
     854        1028 :   s2 = powpolmod(C,R, 2,k, s1);
     855        1028 :   s3 = autvec_AL(pk, j3q, cache_E(C), R);
     856        1028 :   s3 = _red(gmul(s3,s2), R);
     857        1028 :   if (j2q) s3 = _red(gmul(j2q, s3), R);
     858             : 
     859        1028 :   ind = look_eta2(k, s3);
     860        1028 :   if (ind < 0) return -1;
     861        1028 :   if ((ind&1)==0) return 0;
     862         448 :   s3 = Fp_pow(utoipos(q), R->N2, R->N);
     863         448 :   return is_m1(s3, R->N);
     864             : }
     865             : 
     866             : /* p=2, k=2 */
     867             : static long
     868        7225 : step4c(GEN C, Red *R, ulong q)
     869             : {
     870             :   long ind;
     871        7225 :   GEN s0,s1,s3, jpq = get_jac2(R->N,q,2, NULL,NULL);
     872             : 
     873        7225 :   s0 = sqrmod4(jpq, R);
     874        7225 :   s1 = gmulsg(q,s0);
     875        7225 :   s3 = powpolmod(C,R, 2,2, s1);
     876        7225 :   if (mod4(R->N) == 3) s3 = _red(gmul(s0,s3), R);
     877             : 
     878        7224 :   ind = look_eta2(2, s3);
     879        7224 :   if (ind < 0) return -1;
     880        7224 :   if ((ind&1)==0) return 0;
     881        3451 :   s3 = Fp_pow(utoipos(q), R->N2, R->N);
     882        3451 :   return is_m1(s3, R->N);
     883             : }
     884             : 
     885             : /* p=2, k=1 */
     886             : static long
     887        9879 : step4d(Red *R, ulong q)
     888             : {
     889        9879 :   GEN s1 = Fp_pow(utoipos(q), R->N2, R->N);
     890        9879 :   if (is_pm1(s1)) return 0;
     891        4801 :   if (is_m1(s1, R->N)) return (mod4(R->N) == 1);
     892           0 :   return -1;
     893             : }
     894             : 
     895             : /* return 1 [OK so far] or 0 [not a prime] */
     896             : static long
     897           7 : step5(GEN pC, Red *R, long p, GEN et, ulong ltab, long lpC)
     898             : {
     899             :   pari_sp av;
     900             :   ulong q;
     901           7 :   long pk, k, fl = -1;
     902             :   GEN C, Cp;
     903             :   forprime_t T;
     904             : 
     905           7 :   (void)u_forprime_arith_init(&T, 3, ULONG_MAX, 1,p);
     906          21 :   while( (q = u_forprime_next(&T)) )
     907             :   { /* q = 1 (mod p) */
     908          21 :     if (umodiu(et,q) == 0) continue;
     909           7 :     if (umodiu(R->N,q) == 0) return 0;
     910           7 :     k = u_lval(q-1, p);
     911           7 :     pk = upowuu(p,k);
     912           7 :     if (pk <= lpC && !isintzero(gel(pC,pk))) {
     913           0 :       C = gel(pC,pk);
     914           0 :       Cp = gel(pC,p);
     915             :     } else {
     916           7 :       C = gel(pC,1);
     917           7 :       Cp = NULL;
     918           7 :       cache_mat(C) = gen_0; /* re-init */
     919             :     }
     920           7 :     av = avma;
     921           7 :     if (!filltabs(C, Cp, R, p, pk, ltab)) return 0;
     922           7 :     R->C = cache_cyc(C);
     923           7 :     if (p >= 3)      fl = step4a(C,R, q,p,k, NULL);
     924           0 :     else if (k >= 3) fl = step4b(C,R, q,k);
     925           0 :     else if (k == 2) fl = step4c(C,R, q);
     926           0 :     else             fl = step4d(R, q);
     927           7 :     if (fl == -1) return 0;
     928           7 :     if (fl == 1) return 1; /*OK*/
     929           0 :     set_avma(av);
     930             :   }
     931           0 :   pari_err_BUG("aprcl test fails! This is highly improbable");
     932             :   return 0;/*LCOV_EXCL_LINE*/
     933             : }
     934             : 
     935             : GEN
     936        1901 : aprcl_step6_worker(GEN r, long t, GEN N, GEN N1, GEN et)
     937             : {
     938             :   long i;
     939        1901 :   pari_sp av = avma;
     940     2695831 :   for (i=1; i<=t; i++)
     941             :   {
     942     2694014 :     r = remii(mulii(r,N1), et);
     943     2700247 :     if (equali1(r)) break;
     944     2695405 :     if (dvdii(N,r) && !equalii(r,N)) return gen_0; /* not prime */
     945     2695087 :     if ((i & 0x1f) == 0) r = gerepileuptoint(av, r);
     946             :   }
     947        1882 :   return cgetg(1,t_VECSMALL);
     948             : }
     949             : 
     950             : /* return 1 if N prime, else 0 */
     951             : static long
     952        1672 : step6(GEN N, ulong t, GEN et)
     953             : {
     954        1672 :   GEN worker, r, rk, N1 = remii(N, et);
     955        1672 :   ulong k = 10000;
     956             :   ulong i;
     957        1672 :   long pending = 0, res = 1;
     958             :   struct pari_mt pt;
     959             :   pari_sp btop;
     960        1672 :   worker = snm_closure(is_entry("_aprcl_step6_worker"), mkvec3(N, N1, et));
     961        1672 :   r = gen_1;
     962        1672 :   rk = Fp_powu(N1, k, et);
     963        1672 :   mt_queue_start_lim(&pt, worker, (t-1+k-1)/k);
     964        1672 :   btop = avma;
     965        3586 :   for (i=1; (i<t && res) || pending; i+=k)
     966             :   {
     967             :     GEN done;
     968        1914 :     mt_queue_submit(&pt, i, (i<t && res)? mkvec2(r, utoi(minss(k,t-i))): NULL);
     969        1914 :     done = mt_queue_get(&pt, NULL, &pending);
     970        1914 :     if (res && done && typ(done) != t_VECSMALL) res = 0; /* not a prime */
     971        1914 :     r = Fp_mul(r, rk, et);
     972        1914 :     if (gc_needed(btop, 2))
     973             :     {
     974           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"APRCL: i = %ld",i);
     975           0 :       r = gerepileupto(btop, r);
     976             :     }
     977             :   }
     978        1672 :   mt_queue_end(&pt);
     979        1672 :   return res;
     980             : }
     981             : 
     982             : GEN
     983       18132 : aprcl_step4_worker(ulong q, GEN pC, GEN N, GEN v)
     984             : {
     985       18132 :   pari_sp av1 = avma, av2 = avma;
     986             :   long j, k;
     987             :   Red R;
     988       18132 :   GEN faq = factoru_pow(q-1), tabg = compute_g(q);
     989       18132 :   GEN P = gel(faq,1), E = gel(faq,2), PE = gel(faq,3);
     990       18132 :   long lfaq = lg(P);
     991       18132 :   GEN flags = cgetg(lfaq, t_VECSMALL);
     992       18132 :   R.N = N;
     993       18132 :   R.N2= shifti(N, -1);
     994       18132 :   R.k = v[1]; R.lv = v[2]; R.mask = uel(v,3); R.n = v[4];
     995       18132 :   av2 = avma;
     996       61697 :   for (j=1, k=1; j<lfaq; j++, set_avma(av2))
     997             :   {
     998       43566 :     long p = P[j], e = E[j], pe = PE[j], fl;
     999       43566 :     GEN C = gel(pC,pe);
    1000       43566 :     R.C = cache_cyc(C);
    1001       43566 :     if (p >= 3)      fl = step4a(C,&R, q,p,e, tabg);
    1002       18132 :     else if (e >= 3) fl = step4b(C,&R, q,e);
    1003       17104 :     else if (e == 2) fl = step4c(C,&R, q);
    1004        9879 :     else             fl = step4d(&R, q);
    1005       43565 :     if (fl == -1) return gen_0; /* not prime */
    1006       43565 :     if (fl == 1) flags[k++] = p;
    1007             :   }
    1008       18131 :   setlg(flags, k);
    1009       18131 :   return gerepileuptoleaf(av1, flags); /* OK so far */
    1010             : }
    1011             : 
    1012             : /* return 1 if prime, else 0 */
    1013             : static long
    1014        1693 : aprcl(GEN N)
    1015             : {
    1016        1693 :   GEN pC, worker, et, fat, flaglp, faet = NULL; /*-Wall*/
    1017        1693 :   long i, j, l, ltab, lfat, lpC, workid, pending = 0, res = 1;
    1018             :   ulong t;
    1019             :   Red R;
    1020             :   struct pari_mt pt;
    1021             : 
    1022        1693 :   if (typ(N) != t_INT) pari_err_TYPE("aprcl",N);
    1023        1693 :   if (cmpis(N,12) <= 0)
    1024          21 :     switch(itos(N))
    1025             :     {
    1026          14 :       case 2: case 3: case 5: case 7: case 11: return 1;
    1027           7 :       default: return 0;
    1028             :     }
    1029        1672 :   if (Z_issquare(N)) return 0;
    1030        1672 :   t = compute_t(N, &et, &faet);
    1031        1672 :   if (DEBUGLEVEL) err_printf("Starting APRCL with t = %ld\n",t);
    1032        1672 :   if (cmpii(sqri(et),N) < 0) pari_err_BUG("aprcl: e(t) too small");
    1033        1672 :   if (!equali1(gcdii(N,mului(t,et)))) return 0;
    1034             : 
    1035        1672 :   R.N = N;
    1036        1672 :   R.N2= shifti(N, -1);
    1037        1672 :   pC = calcglobs(&R, t, &lpC, &ltab, &fat);
    1038        1672 :   if (!pC) return 0;
    1039        1672 :   lfat = lg(fat);
    1040        1672 :   flaglp = cgetg(lfat, t_VECSMALL);
    1041        1672 :   flaglp[1] = 0;
    1042        5986 :   for (i=2; i<lfat; i++)
    1043             :   {
    1044        4314 :     ulong p = fat[i];
    1045        4314 :     flaglp[i] = absequaliu(Fp_powu(N, p-1, sqru(p)), 1);
    1046             :   }
    1047        1672 :   vecsmall_sort(faet);
    1048        1672 :   l = lg(faet);
    1049        1672 :   worker = snm_closure(is_entry("_aprcl_step4_worker"),
    1050        1672 :                        mkvec3(pC, N, mkvecsmall4(R.k, R.lv, R.mask, R.n)));
    1051        1672 :   if (DEBUGLEVEL>2) err_printf("Step4: %ld q-values\n", l-1);
    1052        1672 :   mt_queue_start_lim(&pt, worker, l-1);
    1053       21872 :   for (i=l-1; (i>0 && res) || pending; i--)
    1054             :   {
    1055             :     GEN done;
    1056       20200 :     ulong q = i>0 ? faet[i]: 0;
    1057       20200 :     mt_queue_submit(&pt, q, q>0? mkvec(utoi(q)): NULL);
    1058       20200 :     done = mt_queue_get(&pt, &workid, &pending);
    1059       20200 :     if (res && done)
    1060             :     {
    1061       18132 :       long lf = lg(done);
    1062       18132 :       if (typ(done) != t_VECSMALL) { res = 0; continue; } /* not prime */
    1063       43823 :       for (j=1; j<lf; j++) flaglp[ zv_search(fat, done[j]) ] = 0;
    1064       18132 :       if (DEBUGLEVEL>2) err_printf("testing Jacobi sums for q = %ld...OK\n", workid);
    1065             :     }
    1066             :   }
    1067        1672 :   mt_queue_end(&pt);
    1068        1672 :   if (!res) return 0;
    1069        1672 :   if (DEBUGLEVEL>2) err_printf("Step5: testing conditions lp\n");
    1070        5986 :   for (i=2; i<lfat; i++) /*skip fat[1] = 2*/
    1071             :   {
    1072        4314 :     pari_sp av = avma;
    1073        4314 :     long p = fat[i];
    1074        4314 :     if (flaglp[i] && !step5(pC, &R, p, et, ltab, lpC)) return 0;
    1075        4314 :     set_avma(av);
    1076             :   }
    1077        1672 :   if (DEBUGLEVEL>2)
    1078           0 :     err_printf("Step6: testing potential divisors\n");
    1079        1672 :   return step6(N, t, et);
    1080             : }
    1081             : 
    1082             : long
    1083        1693 : isprimeAPRCL(GEN N)
    1084        1693 : { pari_sp av = avma; return gc_long(av, aprcl(N)); }
    1085             : 
    1086             : /*******************************************************************/
    1087             : /*              DIVISORS IN RESIDUE CLASSES (LENSTRA)              */
    1088             : /*******************************************************************/
    1089             : /* This would allow to replace e(t) > N^(1/2) by e(t) > N^(1/3), but step6
    1090             :  * becomes so expensive that, at least up to 6000 bits, this is useless
    1091             :  * in this application. */
    1092             : static void
    1093         385 : set_add(hashtable *H, void *d)
    1094             : {
    1095         385 :   ulong h = H->hash(d);
    1096         385 :   if (!hash_search2(H, d, h)) hash_insert2(H, d, NULL, h);
    1097         385 : }
    1098             : static GEN
    1099          63 : GEN_hash_keys(hashtable *H)
    1100          63 : { GEN v = hash_keys(H); settyp(v, t_VEC); return ZV_sort(v); }
    1101             : static void
    1102         336 : add(hashtable *H, GEN t1, GEN t2, GEN a, GEN b, GEN r, GEN s)
    1103             : {
    1104         336 :   GEN ra, qa = dvmdii(t1, a, &ra);
    1105         336 :   if (!signe(ra) && dvdii(t2, b) && equalii(modii(qa, s), r))
    1106         266 :     set_add(H, (void*)qa);
    1107         336 : }
    1108             : /* T^2 - B*T + C has integer roots ? */
    1109             : static void
    1110         301 : check_t(hashtable *H, GEN B, GEN C4, GEN a, GEN b, GEN r, GEN s)
    1111             : {
    1112         301 :   GEN d, t1, t2, D = subii(sqri(B), C4);
    1113         301 :   if (!Z_issquareall(D, &d)) return;
    1114         245 :   t1 = shifti(addii(B, d), -1); /* >= 0 */
    1115         245 :   t2 = subii(B, t1);
    1116         245 :   add(H, t1,t2, a,b,r,s);
    1117         245 :   if (signe(t2) >= 0) add(H, t2,t1, a,b,r,s);
    1118             : }
    1119             : /* N > s > r >= 0, (r,s) = 1 */
    1120             : GEN
    1121          63 : divisorslenstra(GEN N, GEN r, GEN s)
    1122             : {
    1123          63 :   pari_sp av = avma;
    1124             :   GEN u, Ns2, rp, a0, a1, b0, b1, c0, c1, s2;
    1125          63 :   hashtable *H = hash_create(11, (ulong(*)(void*))&hash_GEN,
    1126             :                                  (int(*)(void*,void*))&equalii, 1);
    1127             :   long j;
    1128          63 :   if (typ(N) != t_INT) pari_err_TYPE("divisorslenstra", N);
    1129          63 :   if (typ(r) != t_INT) pari_err_TYPE("divisorslenstra", r);
    1130          63 :   if (typ(s) != t_INT) pari_err_TYPE("divisorslenstra", s);
    1131          63 :   u = Fp_inv(r, s);
    1132          63 :   rp = Fp_mul(u, N, s); /* r' */
    1133          63 :   s2 = sqri(s);
    1134          63 :   a0 = s;
    1135          63 :   b0 = gen_0;
    1136          63 :   c0 = gen_0;
    1137          63 :   if (dvdii(N, r)) set_add(H, (void*)r); /* case i = 0 */
    1138          63 :   a1 = Fp_mul(u, rp, s); if (!signe(a1)) a1 = s; /* 0 < a1 <= s */
    1139          63 :   b1 = gen_1;
    1140          63 :   c1 = Fp_mul(u, diviiexact(subii(N,mulii(r,rp)), s), s);
    1141          63 :   Ns2 = divii(N, s2);
    1142          63 :   j = 1;
    1143             :   for (;;)
    1144         273 :   {
    1145         336 :     GEN Cs, q, c, ab = mulii(a1,b1);
    1146             :     long i, lC;
    1147         336 :     if (j == 0) /* i even, a1 >= 0 */
    1148             :     {
    1149         168 :       if (!signe(c1)) Cs = mkvec(gen_0);
    1150             :       else
    1151             :       {
    1152         112 :         GEN cs = mulii(c1, s);
    1153         112 :         Cs = mkvec2(subii(cs,s2), cs);
    1154             :       }
    1155             :     }
    1156             :     else
    1157             :     { /* i odd, a1 > 0 */
    1158         168 :       GEN X = shifti(ab,1);
    1159         168 :       c = c1;
    1160             :       /* smallest c >= 2ab, c = c1 (mod s) */
    1161         168 :       if (cmpii(c, X) < 0)
    1162             :       {
    1163         147 :         GEN rX, qX = dvmdii(subii(X,c), s, &rX);
    1164         147 :         if (signe(rX)) qX = addiu(qX,1); /* ceil((X-c)/s) */
    1165         147 :         c = addii(c, mulii(s, qX));
    1166             :       }
    1167         168 :       Cs = (cmpii(c, addii(Ns2,ab)) <= 0)? mkvec(mulii(c,s)): cgetg(1,t_VEC);
    1168             :     }
    1169         336 :     lC = lg(Cs);
    1170         336 :     if (signe(a1))
    1171             :     {
    1172         273 :       GEN abN4 = shifti(mulii(ab, N), 2);
    1173         273 :       GEN B = addii(mulii(a1,r), mulii(b1,rp));
    1174         574 :       for (i = 1; i < lC; i++)
    1175         301 :         check_t(H, addii(B, gel(Cs,i)), abN4, a1, b1, r, s);
    1176             :     }
    1177             :     else
    1178             :     { /* a1 = 0, last batch */
    1179         133 :       for (i = 1; i < lC; i++)
    1180             :       {
    1181          70 :         GEN ry, ys = dvmdii(gel(Cs,i), b1, &ry);
    1182          70 :         if (!signe(ry))
    1183             :         {
    1184          70 :           GEN d = addii(ys, rp);
    1185          70 :           if (signe(d) > 0)
    1186             :           {
    1187          56 :             d = dvmdii(N, d, &ry);
    1188          56 :             if (!signe(ry)) set_add(H, (void*)d);
    1189             :           }
    1190             :         }
    1191             :       }
    1192          63 :       break; /* DONE */
    1193             :     }
    1194         273 :     j = 1-j;
    1195         273 :     q = dvmdii(a0, a1, &c);
    1196         273 :     if (j == 1 && !signe(c)) { q = subiu(q,1); c = a1; }
    1197         273 :     a0 = a1; a1 = c;
    1198         273 :     if (equali1(q)) /* frequent */
    1199             :     {
    1200         119 :       c = subii(b0, b1); b0 = b1; b1 = c;
    1201         119 :       c = Fp_sub(c0, c1, s); c0 = c1; c1 = c;
    1202             :     }
    1203             :     else
    1204             :     {
    1205         154 :       c = subii(b0, mulii(q,b1)); b0 = b1; b1 = c;
    1206         154 :       c = modii(subii(c0, mulii(q,c1)), s); c0 = c1; c1 = c;
    1207             :     }
    1208             :   }
    1209          63 :   return gerepileupto(av, GEN_hash_keys(H));
    1210             : }

Generated by: LCOV version 1.13