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

Generated by: LCOV version 1.13