Code coverage tests

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

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

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

LCOV - code coverage report
Current view: top level - modules - aprcl.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.8.0 lcov report (development 19361-ac4f238) Lines: 511 583 87.7 %
Date: 2016-08-28 06:11:39 Functions: 45 46 97.8 %
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             : /*                                                                 */
      16             : /*                    THE APRCL PRIMALITY TEST                     */
      17             : /*                                                                 */
      18             : /*******************************************************************/
      19             : #include "pari.h"
      20             : #include "paripriv.h"
      21             : 
      22             : typedef struct Red {
      23             : /* global data */
      24             :   GEN N; /* prime we are certifying */
      25             :   GEN N2; /* floor(N/2) */
      26             : /* globa data for flexible window */
      27             :   long k, lv;
      28             :   ulong mask;
      29             : /* reduction data */
      30             :   long n;
      31             :   GEN C; /* polcyclo(n) */
      32             :   GEN (*red)(GEN x, struct Red*);
      33             : } Red;
      34             : 
      35             : typedef struct Cache { /* data attached to p^k */
      36             :   GEN aall, tall;
      37             :   GEN cyc;
      38             :   GEN E;
      39             :   GEN eta;
      40             :   GEN matvite, matinvvite;
      41             :   GEN avite, pkvite;
      42             :   ulong ctsgt; /* DEBUG */
      43             : } Cache;
      44             : 
      45             : static GEN
      46     1109986 : makepoldeg1(GEN c, GEN d)
      47             : {
      48             :   GEN z;
      49     1109986 :   if (signe(c)) {
      50     1109986 :     z = cgetg(4,t_POL); z[1] = evalsigne(1);
      51     1109968 :     gel(z,2) = d; gel(z,3) = c;
      52           0 :   } else if (signe(d)) {
      53           0 :     z = cgetg(3,t_POL); z[1] = evalsigne(1);
      54           0 :     gel(z,2) = d;
      55             :   } else {
      56           0 :     z = cgetg(2,t_POL); z[1] = evalsigne(0);
      57             :   }
      58     1109968 :   return z;
      59             : }
      60             : 
      61             : /* T mod polcyclo(p), assume deg(T) < 2p */
      62             : static GEN
      63     1483088 : red_cyclop(GEN T, long p)
      64             : {
      65             :   long i, d;
      66             :   GEN y, z;
      67             : 
      68     1483088 :   d = degpol(T) - p; /* < p */
      69     1483083 :   if (d <= -2) return T;
      70             : 
      71             :   /* reduce mod (x^p - 1) */
      72     1483083 :   y = ZX_mod_Xnm1(T, p);
      73     1483099 :   z = y+2;
      74             : 
      75             :   /* reduce mod x^(p-1) + ... + 1 */
      76     1483099 :   d = p-1;
      77     1483099 :   if (degpol(y) == d)
      78     1483095 :     for (i=0; i<d; i++) gel(z,i) = subii(gel(z,i), gel(z,d));
      79     1483061 :   return normalizepol_lg(y, d+2);
      80             : }
      81             : 
      82             : /* x t_VECSMALL, as red_cyclo2n_ip */
      83             : static GEN
      84       10783 : u_red_cyclo2n_ip(GEN x, long n)
      85             : {
      86       10783 :   long i, pow2 = 1L<<(n-1);
      87             :   GEN z;
      88             : 
      89       10783 :   for (i = lg(x)-1; i>pow2; i--) x[i-pow2] -= x[i];
      90       11385 :   for (; i>0; i--)
      91       11385 :     if (x[i]) break;
      92       10783 :   i += 2;
      93       10783 :   z = cgetg(i, t_POL); z[1] = evalsigne(1);
      94       10783 :   for (i--; i>=2; i--) gel(z,i) = stoi(x[i-1]);
      95       10783 :   return z;
      96             : }
      97             : /* x t_POL, n > 0. Return x mod polcyclo(2^n) = (x^(2^(n-1)) + 1). IN PLACE */
      98             : static GEN
      99      395229 : red_cyclo2n_ip(GEN x, long n)
     100             : {
     101      395229 :   long i, pow2 = 1L<<(n-1);
     102     3498166 :   for (i = lg(x)-1; i>pow2+1; i--)
     103     3102947 :     if (signe(gel(x,i))) gel(x,i-pow2) = subii(gel(x,i-pow2), gel(x,i));
     104      395219 :   return normalizepol_lg(x, i+1);
     105             : }
     106             : static GEN
     107      393807 : red_cyclo2n(GEN x, long n) { return red_cyclo2n_ip(leafcopy(x), n); }
     108             : 
     109             : /* x a non-zero VECSMALL */
     110             : static GEN
     111       27310 : smallpolrev(GEN x)
     112             : {
     113       27310 :   long i,j, lx = lg(x);
     114             :   GEN y;
     115             : 
     116       27310 :   while (lx-- && x[lx]==0) /* empty */;
     117       27310 :   i = lx+2; y = cgetg(i,t_POL);
     118       27310 :   y[1] = evalsigne(1);
     119       27310 :   for (j=2; j<i; j++) gel(y,j) = stoi(x[j-1]);
     120       27310 :   return y;
     121             : }
     122             : 
     123             : /* x polynomial in t_VECSMALL form, T t_POL return x mod T */
     124             : static GEN
     125       27310 : u_red(GEN x, GEN T) {
     126       27310 :   return RgX_rem(smallpolrev(x), T);
     127             : }
     128             : 
     129             : /* special case R->C = polcyclo(2^n) */
     130             : static GEN
     131      393807 : _red_cyclo2n(GEN x, Red *R) {
     132      393807 :   return centermod_i(red_cyclo2n(x, R->n), R->N, R->N2);
     133             : }
     134             : /* special case R->C = polcyclo(p) */
     135             : static GEN
     136     1483089 : _red_cyclop(GEN x, Red *R) {
     137     1483089 :   return centermod_i(red_cyclop(x, R->n), R->N, R->N2);
     138             : }
     139             : static GEN
     140      584593 : _red(GEN x, Red *R) {
     141      584593 :   return centermod_i(grem(x, R->C), R->N, R->N2);
     142             : }
     143             : static GEN
     144     8648050 : _redsimple(GEN x, Red *R) { return centermodii(x, R->N, R->N2); }
     145             : 
     146             : static GEN
     147     8921888 : sqrmod(GEN x, Red *R) {
     148     8921888 :   return R->red(gsqr(x), R);
     149             : }
     150             : 
     151             : static GEN
     152           0 : sqrconst(GEN pol, Red *R)
     153             : {
     154           0 :   GEN z = cgetg(3,t_POL);
     155           0 :   gel(z,2) = centermodii(sqri(gel(pol,2)), R->N, R->N2);
     156           0 :   z[1] = pol[1]; return z;
     157             : }
     158             : 
     159             : /* pol^2 mod (x^2+x+1, N) */
     160             : static GEN
     161      609268 : sqrmod3(GEN pol, Red *R)
     162             : {
     163             :   GEN a,b,bma,A,B;
     164      609268 :   long lv = lg(pol);
     165             : 
     166      609268 :   if (lv==2) return pol;
     167      609268 :   if (lv==3) return sqrconst(pol, R);
     168      609268 :   a = gel(pol,3);
     169      609268 :   b = gel(pol,2); bma = subii(b,a);
     170      609261 :   A = centermodii(mulii(a,addii(b,bma)), R->N, R->N2);
     171      609267 :   B = centermodii(mulii(bma,addii(a,b)), R->N, R->N2);
     172      609272 :   return makepoldeg1(A,B);
     173             : }
     174             : 
     175             : /* pol^2 mod (x^2+1,N) */
     176             : static GEN
     177      500733 : sqrmod4(GEN pol, Red *R)
     178             : {
     179             :   GEN a,b,A,B;
     180      500733 :   long lv = lg(pol);
     181             : 
     182      500733 :   if (lv==2) return pol;
     183      500733 :   if (lv==3) return sqrconst(pol, R);
     184      500733 :   a = gel(pol,3);
     185      500733 :   b = gel(pol,2);
     186      500733 :   A = centermodii(mulii(a, shifti(b,1)), R->N, R->N2);
     187      500728 :   B = centermodii(mulii(subii(b,a),addii(b,a)), R->N, R->N2);
     188      500731 :   return makepoldeg1(A,B);
     189             : }
     190             : 
     191             : /* pol^2 mod (polcyclo(5),N) */
     192             : static GEN
     193     1117583 : sqrmod5(GEN pol, Red *R)
     194             : {
     195             :   GEN c2,b,c,d,A,B,C,D;
     196     1117583 :   long lv = lg(pol);
     197             : 
     198     1117583 :   if (lv==2) return pol;
     199     1117583 :   if (lv==3) return sqrconst(pol, R);
     200     1117583 :   c = gel(pol,3); c2 = shifti(c,1);
     201     1117557 :   d = gel(pol,2);
     202     1117557 :   if (lv==4)
     203             :   {
     204           0 :     A = sqri(d);
     205           0 :     B = mulii(c2, d);
     206           0 :     C = sqri(c);
     207           0 :     A = centermodii(A, R->N, R->N2);
     208           0 :     B = centermodii(B, R->N, R->N2);
     209           0 :     C = centermodii(C, R->N, R->N2); return mkpoln(3,A,B,C);
     210             :   }
     211     1117557 :   b = gel(pol,4);
     212     1117557 :   if (lv==5)
     213             :   {
     214          98 :     A = mulii(b, subii(c2,b));
     215          98 :     B = addii(sqri(c), mulii(b, subii(shifti(d,1),b)));
     216          98 :     C = subii(mulii(c2,d), sqri(b));
     217          98 :     D = mulii(subii(d,b), addii(d,b));
     218             :   }
     219             :   else
     220             :   { /* lv == 6 */
     221     1117459 :     GEN a = gel(pol,5), a2 = shifti(a,1);
     222             :     /* 2a(d - c) + b(2c - b) */
     223     1117432 :     A = addii(mulii(a2, subii(d,c)), mulii(b, subii(c2,b)));
     224             :     /* c(c - 2a) + b(2d - b) */
     225     1117491 :     B = addii(mulii(c, subii(c,a2)), mulii(b, subii(shifti(d,1),b)));
     226             :     /* (a-b)(a+b) + 2c(d - a) */
     227     1117498 :     C = addii(mulii(subii(a,b),addii(a,b)), mulii(c2,subii(d,a)));
     228             :     /* 2a(b - c) + (d-b)(d+b) */
     229     1117515 :     D = addii(mulii(a2, subii(b,c)), mulii(subii(d,b), addii(d,b)));
     230             :   }
     231     1117607 :   A = centermodii(A, R->N, R->N2);
     232     1117577 :   B = centermodii(B, R->N, R->N2);
     233     1117556 :   C = centermodii(C, R->N, R->N2);
     234     1117582 :   D = centermodii(D, R->N, R->N2); return mkpoln(4,A,B,C,D);
     235             : }
     236             : 
     237             : static GEN
     238     2154063 : _mul(GEN x, GEN y, Red *R) { return R->red(gmul(x,y), R); }
     239             : 
     240             : /* jac^floor(N/pk) mod (N, polcyclo(pk)), flexible window */
     241             : static GEN
     242       60272 : _powpolmod(Cache *C, GEN jac, Red *R, GEN (*_sqr)(GEN, Red *))
     243             : {
     244       60272 :   const GEN taba = C->aall;
     245       60272 :   const GEN tabt = C->tall;
     246       60272 :   const long efin = lg(taba)-1, lv = R->lv;
     247       60272 :   GEN L, res = jac, pol2 = _sqr(res, R);
     248             :   long f;
     249       60272 :   pari_sp av0 = avma, av;
     250             : 
     251       60272 :   L = cgetg(lv+1, t_VEC); gel(L,1) = res;
     252       60270 :   for (f=2; f<=lv; f++) gel(L,f) = _mul(gel(L,f-1), pol2, R);
     253       60270 :   av = avma;
     254     1787213 :   for (f = efin; f >= 1; f--)
     255             :   {
     256     1726941 :     GEN t = gel(L, taba[f]);
     257     1726941 :     long tf = tabt[f];
     258     1726941 :     res = (f==efin)? t: _mul(t, res, R);
     259    14534310 :     while (tf--) {
     260    11080425 :       res = _sqr(res, R);
     261    11080429 :       if (gc_needed(av,1)) {
     262        1694 :         res = gerepilecopy(av, res);
     263        1694 :         if(DEBUGMEM>1) pari_warn(warnmem,"powpolmod: f = %ld",f);
     264             :       }
     265             :     }
     266             :   }
     267       60272 :   return gerepilecopy(av0, res);
     268             : }
     269             : 
     270             : static GEN
     271       13717 : _powpolmodsimple(Cache *C, Red *R, GEN jac)
     272             : {
     273       13717 :   pari_sp av = avma;
     274       13717 :   GEN w = mulmat_pol(C->matvite, jac);
     275       13717 :   long j, ph = lg(w);
     276             : 
     277       13717 :   R->red = &_redsimple;
     278       50995 :   for (j=1; j<ph; j++)
     279       37278 :     gel(w,j) = _powpolmod(C, centermodii(gel(w,j), R->N, R->N2), R, &sqrmod);
     280       13717 :   w = centermod_i( gmul(C->matinvvite, w), R->N, R->N2 );
     281       13717 :   w = gerepileupto(av, w);
     282       13717 :   return RgV_to_RgX(w, 0);
     283             : }
     284             : 
     285             : static GEN
     286       36711 : powpolmod(Cache *C, Red *R, long p, long k, GEN jac)
     287             : {
     288             :   GEN (*_sqr)(GEN, Red *);
     289             : 
     290       36711 :   if (DEBUGLEVEL>2) C->ctsgt++;
     291       36711 :   if (C->matvite) return _powpolmodsimple(C, R, jac);
     292       22994 :   if (p == 2) /* p = 2 */
     293             :   {
     294        4811 :     if (k == 2) _sqr = &sqrmod4;
     295         754 :     else        _sqr = &sqrmod;
     296        4811 :     R->n = k;
     297        4811 :     R->red = &_red_cyclo2n;
     298       18183 :   } else if (k == 1)
     299             :   {
     300       16628 :     if      (p == 3) _sqr = &sqrmod3;
     301       11750 :     else if (p == 5) _sqr = &sqrmod5;
     302        4559 :     else             _sqr = &sqrmod;
     303       16628 :     R->n = p;
     304       16628 :     R->red = &_red_cyclop;
     305             :   } else {
     306        1555 :     R->red = &_red; _sqr = &sqrmod;
     307             :   }
     308       22994 :   return _powpolmod(C, jac, R, _sqr);
     309             : }
     310             : 
     311             : /* Return e(t) = \prod_{p-1 | t} p^{1+v_p(t)}}
     312             :  * faet contains the odd prime divisors of e(t) */
     313             : static GEN
     314        1940 : compute_e(ulong t, GEN *faet)
     315             : {
     316        1940 :   GEN L, P, D = divisorsu(t);
     317        1940 :   long l = lg(D);
     318             :   ulong k;
     319             : 
     320        1940 :   P = vecsmalltrunc_init(l);
     321        1940 :   L = vecsmalltrunc_init(l);
     322       39904 :   for (k=l-1; k>1; k--) /* k != 1: avoid d = 1 */
     323             :   {
     324       37964 :     ulong d = D[k];
     325       37964 :     if (uisprime(++d))
     326             :     { /* we want q = 1 (mod p) prime, not too large */
     327       20459 :       if (d > 50000000) return gen_0;
     328       20459 :       vecsmalltrunc_append(P, d);
     329       20459 :       vecsmalltrunc_append(L, upowuu(d, 1 + u_lval(t,d)));
     330             :     }
     331             :   }
     332        1940 :   if (faet) *faet = P;
     333        1940 :   return shifti(zv_prod_Z(L), 2 + u_lval(t,2));
     334             : }
     335             : 
     336             : /* Table obtained by the following script:
     337             : 
     338             : install(compute_e, "LD&"); \\ remove 'static' first
     339             : 
     340             : table(first = 6, step = 6, MAXT = 8648640)=
     341             : {
     342             :   emax = 0;
     343             :   forstep(t = first, MAXT, step,
     344             :     e = compute_e(t);
     345             :     if (e > 1.9*emax, emax = e;
     346             :       printf("  if (C < %5.5g) return %8d;\n", 2*log(e)/log(2)*0.9999, t)
     347             :     );
     348             :   );
     349             : }
     350             : 
     351             : Previous values can be recovered using the following table:
     352             : 
     353             : T=[6,12,24,48,36,60,120,180,240,504,360,420,720,840,1440,1260,1680,2520,3360,5040,13860,10080,16380,21840,18480,27720,32760,36960,55440,65520,98280,110880,131040,166320,196560,262080,277200,360360,480480,332640,554400,720720,665280,831600,1113840,1441440,1663200,2227680,2162160,2827440,3326400,3603600,6126120,4324320,6683040,7207200,11138400,8648640];
     354             : f(t) = 2*log(compute_e(t))/log(2);
     355             : for (i=1,#T, t=T[i]; printf("  if (C < %5.5g) return %8d;\n", f(t),t));
     356             : 
     357             : */
     358             : 
     359             : /* assume C < 3514.6 */
     360             : static ulong
     361        1940 : compute_t_small(double C)
     362             : {
     363        1940 :   if (C < 17.953) return        6;
     364        1940 :   if (C < 31.996) return       12;
     365        1940 :   if (C < 33.996) return       24;
     366        1940 :   if (C < 54.079) return       36;
     367        1940 :   if (C < 65.325) return       60;
     368        1240 :   if (C < 68.457) return       72;
     369        1240 :   if (C < 70.783) return      108;
     370        1240 :   if (C < 78.039) return      120;
     371        1233 :   if (C < 102.41) return      180;
     372        1149 :   if (C < 127.50) return      360;
     373         942 :   if (C < 136.68) return      420;
     374          32 :   if (C < 153.43) return      540;
     375          32 :   if (C < 165.66) return      840;
     376          32 :   if (C < 169.17) return     1008;
     377          32 :   if (C < 178.52) return     1080;
     378          32 :   if (C < 192.68) return     1200;
     379          32 :   if (C < 206.34) return     1260;
     380          32 :   if (C < 211.94) return     1620;
     381          32 :   if (C < 222.09) return     1680;
     382          32 :   if (C < 225.11) return     2016;
     383          32 :   if (C < 244.19) return     2160;
     384          32 :   if (C < 270.29) return     2520;
     385          25 :   if (C < 279.50) return     3360;
     386          25 :   if (C < 293.62) return     3780;
     387          25 :   if (C < 346.68) return     5040;
     388          25 :   if (C < 348.70) return     6480;
     389          25 :   if (C < 383.34) return     7560;
     390          25 :   if (C < 396.68) return     8400;
     391          25 :   if (C < 426.04) return    10080;
     392          25 :   if (C < 458.34) return    12600;
     393          25 :   if (C < 527.16) return    15120;
     394          25 :   if (C < 595.38) return    25200;
     395          25 :   if (C < 636.29) return    30240;
     396           7 :   if (C < 672.53) return    42840;
     397           7 :   if (C < 684.90) return    45360;
     398           7 :   if (C < 708.78) return    55440;
     399           7 :   if (C < 771.30) return    60480;
     400           7 :   if (C < 775.86) return    75600;
     401           7 :   if (C < 859.62) return    85680;
     402           7 :   if (C < 893.16) return   100800;
     403           7 :   if (C < 912.27) return   110880;
     404           7 :   if (C < 966.13) return   128520;
     405           7 :   if (C < 1009.1) return   131040;
     406           0 :   if (C < 1041.9) return   166320;
     407           0 :   if (C < 1124.9) return   196560;
     408           0 :   if (C < 1251.0) return   257040;
     409           0 :   if (C < 1375.0) return   332640;
     410           0 :   if (C < 1431.0) return   393120;
     411           0 :   if (C < 1483.3) return   514080;
     412           0 :   if (C < 1546.3) return   655200;
     413           0 :   if (C < 1585.8) return   665280;
     414           0 :   if (C < 1661.3) return   786240;
     415           0 :   if (C < 1667.5) return   831600;
     416           0 :   if (C < 1676.9) return   917280;
     417           0 :   if (C < 1728.0) return   982800;
     418           0 :   if (C < 1747.4) return  1081080;
     419           0 :   if (C < 1773.6) return  1179360;
     420           0 :   if (C < 1810.6) return  1285200;
     421           0 :   if (C < 1924.5) return  1310400;
     422           0 :   if (C < 2001.1) return  1441440;
     423           0 :   if (C < 2096.3) return  1663200;
     424           0 :   if (C < 2165.8) return  1965600;
     425           0 :   if (C < 2321.6) return  2162160;
     426           0 :   if (C < 2368.2) return  2751840;
     427           0 :   if (C < 2377.2) return  2827440;
     428           0 :   if (C < 2514.7) return  3326400;
     429           0 :   if (C < 2588.5) return  3341520;
     430           0 :   if (C < 2636.6) return  3603600;
     431           0 :   if (C < 2667.2) return  3931200;
     432           0 :   if (C < 3028.6) return  4324320;
     433           0 :   if (C < 3045.5) return  5654880;
     434           0 :   if (C < 3080.5) return  6652800;
     435           0 :   if (C < 3121.6) return  6683040;
     436           0 :   if (C < 3283.1) return  7207200;
     437           0 :   return  8648640;
     438             : }
     439             : 
     440             : /* return t such that e(t) > sqrt(N), set *faet = odd prime divisors of e(t) */
     441             : static ulong
     442        1940 : compute_t(GEN N, GEN *e, GEN *faet)
     443             : {
     444             :   /* 2^e b <= N < 2^e (b+1), where b >= 2^52. Approximating log_2 N by
     445             :    * log2(gtodouble(N)) ~ e+log2(b), the error is less than log(1+1/b) < 1e-15*/
     446        1940 :   double C = dbllog2(N) + 1e-6; /* > log_2 N */
     447             :   ulong t;
     448             :   GEN B;
     449             :   /* Return "smallest" t such that f(t) >= C, which implies e(t) > sqrt(N) */
     450             :   /* For N < 2^3515 ~ 10^1058 */
     451        1940 :   if (C < 3514.6)
     452             :   {
     453        1940 :     t = compute_t_small(C);
     454        1940 :     *e = compute_e(t, faet);
     455        1940 :     return t;
     456             :   }
     457           0 :   B = sqrti(N);
     458           0 :   for (t = 8648640+840;; t+=840)
     459             :   {
     460           0 :     pari_sp av = avma;
     461           0 :     *e = compute_e(t, faet);
     462           0 :     if (cmpii(*e, B) > 0) break;
     463           0 :     avma = av;
     464           0 :   }
     465           0 :   return t;
     466             : }
     467             : 
     468             : /* T[i] = discrete log of i in (Z/q)^*, q odd prime
     469             :  * To save on memory, compute half the table: T[q-x] = T[x] + (q-1)/2 */
     470             : static GEN
     471       29867 : computetabdl(ulong q)
     472             : {
     473       29867 :   ulong g, a, i, qs2 = q>>1; /* (q-1)/2 */
     474       29867 :   GEN T = cgetg(qs2+2,t_VECSMALL);
     475             : 
     476       29867 :   g = pgener_Fl(q); a = 1;
     477     4375967 :   for (i=1; i < qs2; i++) /* g^((q-1)/2) = -1 */
     478             :   {
     479     4346100 :     a = Fl_mul(g,a,q);
     480     4346100 :     if (a > qs2) T[q-a] = i+qs2; else T[a] = i;
     481             :   }
     482       29867 :   T[qs2+1] = T[qs2] + qs2;
     483       29867 :   T[1] = 0; return T;
     484             : }
     485             : 
     486             : /* Return T: T[x] = dl of x(1-x) */
     487             : static GEN
     488       20466 : compute_g(ulong q)
     489             : {
     490       20466 :   const ulong qs2 = q>>1; /* (q-1)/2 */
     491             :   ulong x, a;
     492       20466 :   GEN T = computetabdl(q); /* updated in place to save on memory */
     493       20466 :   a = 0; /* dl[1] */
     494     2342079 :   for (x=2; x<=qs2+1; x++)
     495             :   { /* a = dl(x) */
     496     2321613 :     ulong b = T[x]; /* = dl(x) */
     497     2321613 :     T[x] = b + a + qs2; /* dl(x) + dl(x-1) + dl(-1) */
     498     2321613 :     a = b;
     499             :   }
     500       20466 :   return T;
     501             : }
     502             : 
     503             : /* p odd prime */
     504             : static GEN
     505       27310 : get_jac(Cache *C, ulong q, long pk, GEN tabg)
     506             : {
     507       27310 :   ulong x, qs2 = q>>1; /* (q-1)/2 */
     508       27310 :   GEN vpk = zero_zv(pk);
     509             : 
     510       27310 :   for (x=2; x<=qs2; x++) vpk[ tabg[x]%pk + 1 ] += 2;
     511       27310 :   vpk[ tabg[x]%pk + 1 ]++; /* x = (q+1)/2 */
     512       27310 :   return u_red(vpk, C->cyc);
     513             : }
     514             : 
     515             : /* p = 2 */
     516             : static GEN
     517        9401 : get_jac2(GEN N, ulong q, long k, GEN *j2q, GEN *j3q)
     518             : {
     519        9401 :   GEN jpq, vpk, T = computetabdl(q);
     520             :   ulong x, pk, i, qs2;
     521             : 
     522             :   /* could store T[x+1] + T[x] + qs2 (cf compute_g).
     523             :    * Recompute instead, saving half the memory. */
     524        9401 :   pk = 1UL << k;;
     525        9401 :   vpk = zero_zv(pk);
     526             : 
     527        9401 :   qs2 = q>>1; /* (q-1)/2 */
     528             : 
     529        9401 :   for (x=2; x<=qs2; x++) vpk[ (T[x]+T[x-1]+qs2)%pk + 1 ] += 2;
     530        9401 :   vpk[ (T[x]+T[x-1]+qs2)%pk + 1 ]++;
     531        9401 :   jpq = u_red_cyclo2n_ip(vpk, k);
     532             : 
     533        9401 :   if (k == 2) return jpq;
     534             : 
     535         962 :   if (mod8(N) >= 5)
     536             :   {
     537         420 :     GEN v8 = cgetg(9,t_VECSMALL);
     538         420 :     for (x=1; x<=8; x++) v8[x] = 0;
     539         420 :     for (x=2; x<=qs2; x++) v8[ ((3*T[x]+T[x-1]+qs2)&7) + 1 ]++;
     540         420 :     for (   ; x<=q-1; x++) v8[ ((3*T[q-x]+T[q-x+1]-3*qs2)&7) + 1 ]++;
     541         420 :     *j2q = RgX_inflate(ZX_sqr(u_red_cyclo2n_ip(v8,3)), pk>>3);
     542         420 :     *j2q = red_cyclo2n_ip(*j2q, k);
     543             :   }
     544         962 :   for (i=1; i<=pk; i++) vpk[i] = 0;
     545         962 :   for (x=2; x<=qs2; x++) vpk[ (2*T[x]+T[x-1]+qs2)%pk + 1 ]++;
     546         962 :   for (   ; x<=q-1; x++) vpk[ (2*T[q-x]+T[q-x+1]-2*qs2)%pk + 1 ]++;
     547         962 :   *j3q = ZX_mul(jpq, u_red_cyclo2n_ip(vpk,k));
     548         962 :   *j3q = red_cyclo2n_ip(*j3q, k);
     549         962 :   return jpq;
     550             : }
     551             : 
     552             : /* N = 1 mod p^k, return an elt of order p^k in (Z/N)^* */
     553             : static GEN
     554        2743 : finda(Cache *Cp, GEN N, long pk, long p)
     555             : {
     556             :   GEN a, pv;
     557        2743 :   if (Cp && Cp->avite) {
     558          70 :     a  = Cp->avite;
     559          70 :     pv = Cp->pkvite;
     560             :   }
     561             :   else
     562             :   {
     563             :     GEN ph, b, q;
     564        2673 :     ulong u = 2;
     565        2673 :     long v = Z_lvalrem(addis(N,-1), p, &q);
     566        2673 :     ph = powuu(p, v-1); pv = muliu(ph, p); /* N - 1 = p^v q */
     567        2673 :     if (p > 2)
     568             :     {
     569         931 :       for (;;u++)
     570             :       {
     571        2595 :         a = Fp_pow(utoipos(u), q, N);
     572        2595 :         b = Fp_pow(a, ph, N);
     573        2595 :         if (!gequal1(b)) break;
     574         931 :       }
     575             :     }
     576             :     else
     577             :     {
     578        1009 :       while (krosi(u,N) >= 0) u++;
     579        1009 :       a = Fp_pow(utoipos(u), q, N);
     580        1009 :       b = Fp_pow(a, ph, N);
     581             :     }
     582             :     /* checking b^p = 1 mod N done economically in caller */
     583        2673 :     b = gcdii(addis(b,-1), N);
     584        2673 :     if (!gequal1(b)) return NULL;
     585             : 
     586        2673 :     if (Cp) {
     587        2673 :       Cp->avite  = a; /* a has order p^v */
     588        2673 :       Cp->pkvite = pv;
     589             :     }
     590             :   }
     591        2743 :   return Fp_pow(a, divis(pv, pk), N);
     592             : }
     593             : 
     594             : /* return 0: N not a prime, 1: no problem so far */
     595             : static long
     596        9353 : filltabs(Cache *C, Cache *Cp, Red *R, long p, long pk, long ltab)
     597             : {
     598             :   pari_sp av;
     599             :   long i, j;
     600             :   long e;
     601             :   GEN tabt, taba, m;
     602             : 
     603        9353 :   C->cyc = polcyclo(pk,0);
     604             : 
     605        9353 :   if (p > 2)
     606             :   {
     607        5177 :     long LE = pk - pk/p + 1;
     608        5177 :     GEN E = cgetg(LE, t_VECSMALL), eta = cgetg(pk+1,t_VEC);
     609       25661 :     for (i=1,j=0; i<pk; i++)
     610       20484 :       if (i%p) E[++j] = i;
     611        5177 :     C->E = E;
     612             : 
     613       30838 :     for (i=1; i<=pk; i++)
     614             :     {
     615       25661 :       GEN z = FpX_rem(monomial(gen_1, i-1, 0), C->cyc, R->N);
     616       25661 :       gel(eta,i) = FpX_center(z, R->N, R->N2);
     617             :     }
     618        5177 :     C->eta = eta;
     619             :   }
     620        4176 :   else if (pk >= 8)
     621             :   {
     622         296 :     long LE = (pk>>2) + 1;
     623         296 :     GEN E = cgetg(LE, t_VECSMALL);
     624        3168 :     for (i=1,j=0; i<pk; i++)
     625        2872 :       if ((i%8)==1 || (i%8)==3) E[++j] = i;
     626         296 :     C->E = E;
     627             :   }
     628             : 
     629        9353 :   if (pk > 2 && umodiu(R->N,pk) == 1)
     630             :   {
     631        2743 :     GEN vpa, p1, p2, p3, a2 = NULL, a = finda(Cp, R->N, pk, p);
     632             :     long jj, ph;
     633             : 
     634        2743 :     if (!a) return 0;
     635        2743 :     ph = pk - pk/p;
     636        2743 :     vpa = cgetg(ph+1,t_COL); gel(vpa,1) = a;
     637        2743 :     if (pk > p) a2 = centermodii(sqri(a), R->N, R->N2);
     638        2743 :     jj = 1;
     639        8419 :     for (i=2; i<pk; i++) /* vpa = { a^i, (i,p) = 1 } */
     640        5676 :       if (i%p) {
     641        4447 :         GEN z = mulii((i%p==1) ? a2 : a, gel(vpa,jj));
     642        4447 :         gel(vpa,++jj) = centermodii(z , R->N, R->N2);
     643             :       }
     644        2743 :     if (!gequal1( centermodii( mulii(a, gel(vpa,ph)), R->N, R->N2) )) return 0;
     645        2743 :     p1 = cgetg(ph+1,t_MAT);
     646        2743 :     p2 = cgetg(ph+1,t_COL); gel(p1,1) = p2;
     647        2743 :     for (i=1; i<=ph; i++) gel(p2,i) = gen_1;
     648        2743 :     j = 2; gel(p1,j) = vpa; p3 = vpa;
     649        4447 :     for (j++; j <= ph; j++)
     650             :     {
     651        1704 :       p2 = cgetg(ph+1,t_COL); gel(p1,j) = p2;
     652       10168 :       for (i=1; i<=ph; i++)
     653        8464 :         gel(p2,i) = centermodii(mulii(gel(vpa,i),gel(p3,i)), R->N, R->N2);
     654        1704 :       p3 = p2;
     655             :     }
     656        2743 :     C->matvite = p1;
     657        2743 :     C->matinvvite = FpM_inv(p1, R->N);
     658             :   }
     659             : 
     660        9353 :   tabt = cgetg(ltab+1, t_VECSMALL);
     661        9353 :   taba = cgetg(ltab+1, t_VECSMALL);
     662        9353 :   av = avma; m = divis(R->N, pk);
     663      174161 :   for (e=1; e<=ltab && signe(m); e++)
     664             :   {
     665      164809 :     long s = vali(m); m = shifti(m,-s);
     666      164803 :     tabt[e] = e==1? s: s + R->k;
     667      164803 :     taba[e] = signe(m)? ((mod2BIL(m) & R->mask)+1)>>1: 0;
     668      164805 :     m = shifti(m, -R->k);
     669             :   }
     670        9352 :   setlg(taba, e); C->aall = taba;
     671        9353 :   setlg(tabt, e); C->tall = tabt;
     672        9353 :   avma = av; return 1;
     673             : }
     674             : 
     675             : static Cache *
     676       11286 : alloc_cache(void)
     677             : {
     678       11286 :   Cache *C = (Cache*)stack_malloc(sizeof(Cache));
     679       11286 :   C->matvite = NULL;
     680       11286 :   C->avite   = NULL;
     681       11286 :   C->ctsgt = 0; return C;
     682             : }
     683             : 
     684             : static Cache **
     685        1940 : calcglobs(Red *R, ulong t, long *plpC, long *pltab, GEN *pP)
     686             : {
     687             :   GEN fat, P, E, PE;
     688             :   long lv, i, k, b;
     689             :   Cache **pC;
     690             : 
     691        1940 :   b = expi(R->N)+1;
     692        1940 :   k = 3; while (((k+1)*(k+2) << (k-1)) < b) k++;
     693        1940 :   *pltab = (b/k)+2;
     694        1940 :   R->k  = k;
     695        1940 :   R->lv = 1L << (k-1);
     696        1940 :   R->mask = (1UL << k) - 1;
     697             : 
     698        1940 :   fat = factoru_pow(t);
     699        1940 :   P = gel(fat,1);
     700        1940 :   E = gel(fat,2);
     701        1940 :   PE= gel(fat,3);
     702        1940 :   *plpC = lv = vecsmall_max(PE); /* max(p^e, p^e | t) */
     703        1940 :   pC = (Cache**)stack_malloc((lv+1)*sizeof(Cache*));
     704        1940 :   pC[1] = alloc_cache(); /* to be used as temp in step5() */
     705        1940 :   for (i = 2; i <= lv; i++) pC[i] = NULL;
     706        8709 :   for (i=1; i<lg(P); i++)
     707             :   {
     708        6769 :     long pk, p = P[i], e = E[i];
     709        6769 :     pk = p;
     710       16115 :     for (k=1; k<=e; k++, pk*=p)
     711             :     {
     712        9346 :       pC[pk] = alloc_cache();
     713        9346 :       if (!filltabs(pC[pk], pC[p], R, p,pk, *pltab)) return NULL;
     714             :     }
     715             :   }
     716        1940 :   *pP = P; return pC;
     717             : }
     718             : 
     719             : /* sig_a^{-1}(z) for z in Q(zeta_pk) and sig_a: zeta -> zeta^a. Assume
     720             :  * a reduced mod pk := p^k*/
     721             : static GEN
     722      163505 : aut(long pk, GEN z, long a)
     723             : {
     724             :   GEN v;
     725      163505 :   long b, i, dz = degpol(z);
     726      163504 :   if (a == 1 || dz < 0) return z;
     727      135232 :   v = cgetg(pk+2,t_POL);
     728      135233 :   v[1] = evalvarn(0);
     729      135233 :   b = 0;
     730      135233 :   gel(v,2) = gel(z,2); /* i = 0 */
     731      998567 :   for (i = 1; i < pk; i++)
     732             :   {
     733      863334 :     b += a; if (b > pk) b -= pk; /* b = (a*i) % pk */
     734      863334 :     gel(v,i+2) = b > dz? gen_0: gel(z,b+2);
     735             :   }
     736      135233 :   return normalizepol_lg(v, pk+2);
     737             : }
     738             : 
     739             : /* z^v for v in Z[G], represented by couples [sig_x^{-1},x] */
     740             : static GEN
     741       28272 : autvec_TH(long pk, GEN z, GEN v, GEN C)
     742             : {
     743       28272 :   long i, lv = lg(v);
     744       28272 :   GEN s = pol_1(varn(C));
     745      136424 :   for (i=1; i<lv; i++)
     746             :   {
     747      108152 :     long y = v[i];
     748      108152 :     if (y) s = RgXQ_mul(s, RgXQ_powu(aut(pk,z, y), y, C), C);
     749             :   }
     750       28272 :   return s;
     751             : }
     752             : 
     753             : static GEN
     754       28272 : autvec_AL(long pk, GEN z, GEN v, Red *R)
     755             : {
     756       28272 :   const long r = umodiu(R->N, pk);
     757       28272 :   GEN s = pol_1(varn(R->C));
     758       28272 :   long i, lv = lg(v);
     759      136424 :   for (i=1; i<lv; i++)
     760             :   {
     761      108152 :     long y = (r*v[i]) / pk;
     762      108152 :     if (y) s = RgXQ_mul(s, RgXQ_powu(aut(pk,z, v[i]), y, R->C), R->C);
     763             :   }
     764       28272 :   return s;
     765             : }
     766             : 
     767             : /* 0 <= i < pk, such that x^i = z mod polcyclo(pk),  -1 if no such i exist */
     768             : static long
     769       27310 : look_eta(GEN eta, long pk, GEN z)
     770             : {
     771             :   long i;
     772       82341 :   for (i=1; i<=pk; i++)
     773       82341 :     if (ZX_equal(z, gel(eta,i))) return i-1;
     774           0 :   return -1;
     775             : }
     776             : /* same pk = 2^k */
     777             : static long
     778        9401 : look_eta2(long k, GEN z)
     779             : {
     780             :   long d, s;
     781        9401 :   if (typ(z) != t_POL) d = 0; /* t_INT */
     782             :   else
     783             :   {
     784        9401 :     if (!RgX_is_monomial(z)) return -1;
     785        9401 :     d = degpol(z);
     786        9401 :     z = gel(z,d+2); /* leading term */
     787             :   }
     788        9401 :   s = signe(z);
     789        9401 :   if (!s || !is_pm1(z)) return -1;
     790        9401 :   return (s > 0)? d: d + (1L<<(k-1));
     791             : }
     792             : 
     793             : static long
     794       27310 : step4a(Cache *C, Red *R, ulong q, long p, long k, GEN tabg)
     795             : {
     796       27310 :   const long pk = upowuu(p,k);
     797             :   long ind;
     798             :   GEN jpq, s1, s2, s3;
     799             : 
     800       27310 :   if (!tabg) tabg = compute_g(q);
     801       27310 :   jpq = get_jac(C, q, pk, tabg);
     802       27310 :   s1 = autvec_TH(pk, jpq, C->E, C->cyc);
     803       27310 :   s2 = powpolmod(C,R, p,k, s1);
     804       27310 :   s3 = autvec_AL(pk, jpq, C->E, R);
     805       27310 :   s3 = _red(gmul(s3,s2), R);
     806             : 
     807       27310 :   ind = look_eta(C->eta, pk, s3);
     808       27310 :   if (ind < 0) return -1;
     809       27310 :   return (ind%p) != 0;
     810             : }
     811             : 
     812             : /* x == -1 mod N ? */
     813             : static long
     814        9651 : is_m1(GEN x, GEN N)
     815             : {
     816        9651 :   return equalii(addis(x,1), N);
     817             : }
     818             : 
     819             : /* p=2, k>=3 */
     820             : static long
     821         962 : step4b(Cache *C, Red *R, ulong q, long k)
     822             : {
     823         962 :   const long pk = 1L << k;
     824             :   long ind;
     825         962 :   GEN s1, s2, s3, j2q = NULL, j3q = NULL;
     826             : 
     827         962 :   (void)get_jac2(R->N,q,k, &j2q,&j3q);
     828             : 
     829         962 :   s1 = autvec_TH(pk, j3q, C->E, C->cyc);
     830         962 :   s2 = powpolmod(C,R, 2,k, s1);
     831         962 :   s3 = autvec_AL(pk, j3q, C->E, R);
     832         962 :   s3 = _red(gmul(s3,s2), R);
     833         962 :   if (j2q) s3 = _red(gmul(j2q, s3), R);
     834             : 
     835         962 :   ind = look_eta2(k, s3);
     836         962 :   if (ind < 0) return -1;
     837         962 :   if ((ind&1)==0) return 0;
     838         429 :   if (DEBUGLEVEL>2) C->ctsgt++;
     839         429 :   s3 = Fp_pow(utoipos(q), R->N2, R->N);
     840         429 :   return is_m1(s3, R->N);
     841             : }
     842             : 
     843             : /* p=2, k=2 */
     844             : static long
     845        8439 : step4c(Cache *C, Red *R, ulong q)
     846             : {
     847             :   long ind;
     848        8439 :   GEN s0,s1,s3, jpq = get_jac2(R->N,q,2, NULL,NULL);
     849             : 
     850        8439 :   s0 = sqrmod4(jpq, R);
     851        8439 :   s1 = gmulsg(q,s0);
     852        8439 :   s3 = powpolmod(C,R, 2,2, s1);
     853        8439 :   if (mod4(R->N) == 3) s3 = _red(gmul(s0,s3), R);
     854             : 
     855        8439 :   ind = look_eta2(2, s3);
     856        8439 :   if (ind < 0) return -1;
     857        8439 :   if ((ind&1)==0) return 0;
     858        4018 :   if (DEBUGLEVEL>2) C->ctsgt++;
     859        4018 :   s3 = Fp_pow(utoipos(q), R->N2, R->N);
     860        4018 :   return is_m1(s3, R->N);
     861             : }
     862             : 
     863             : /* p=2, k=1 */
     864             : static long
     865       11058 : step4d(Cache *C, Red *R, ulong q)
     866             : {
     867       11058 :   GEN s1 = Fp_pow(utoipos(q), R->N2, R->N);
     868       11058 :   if (DEBUGLEVEL>2) C->ctsgt++;
     869       11058 :   if (is_pm1(s1)) return 0;
     870        5204 :   if (is_m1(s1, R->N)) return (mod4(R->N) == 1);
     871           0 :   return -1;
     872             : }
     873             : 
     874             : static GEN
     875           7 : _res(long a, long b) { return b? mkvec2s(a, b): mkvecs(a); }
     876             : 
     877             : /* return 1 [OK so far] or <= 0 [not a prime] */
     878             : static GEN
     879           7 : step5(Cache **pC, Red *R, long p, GEN et, ulong ltab, long lpC)
     880             : {
     881             :   pari_sp av;
     882             :   ulong q;
     883           7 :   long pk, k, fl = -1;
     884             :   Cache *C, *Cp;
     885             :   forprime_t T;
     886             : 
     887           7 :   (void)u_forprime_arith_init(&T, 3, ULONG_MAX, 1,p);
     888          28 :   while( (q = u_forprime_next(&T)) )
     889             :   { /* q = 1 (mod p) */
     890          21 :     if (umodiu(et,q) == 0) continue;
     891           7 :     if (umodiu(R->N,q) == 0) return _res(1,p);
     892           7 :     k = u_lval(q-1, p);
     893           7 :     pk = upowuu(p,k);
     894           7 :     if (pk <= lpC && pC[pk]) {
     895           0 :       C = pC[pk];
     896           0 :       Cp = pC[p];
     897             :     } else {
     898           7 :       C = pC[1];
     899           7 :       Cp = NULL;
     900           7 :       C->matvite = NULL; /* re-init */
     901             :     }
     902             : 
     903           7 :     av = avma;
     904           7 :     if (!filltabs(C, Cp, R, p, pk, ltab)) return _res(1,0);
     905           7 :     R->C = C->cyc;
     906           7 :     if (p >= 3)      fl = step4a(C,R, q,p,k, NULL);
     907           0 :     else if (k >= 3) fl = step4b(C,R, q,k);
     908           0 :     else if (k == 2) fl = step4c(C,R, q);
     909           0 :     else             fl = step4d(C,R, q);
     910           7 :     if (fl == -1) return _res(q,p);
     911           7 :     if (fl == 1) return NULL; /*OK*/
     912           0 :     avma = av;
     913             :   }
     914           0 :   pari_err_BUG("aprcl test fails! This is highly improbable");
     915           0 :   return NULL;
     916             : }
     917             : 
     918             : static GEN
     919        1940 : step6(GEN N, ulong t, GEN et)
     920             : {
     921        1940 :   GEN r, N1 = remii(N, et);
     922             :   ulong i;
     923        1940 :   pari_sp av = avma;
     924             : 
     925        1940 :   r = gen_1;
     926     1982080 :   for (i=1; i<t; i++)
     927             :   {
     928     1980233 :     r = remii(mulii(r,N1), et);
     929     1980233 :     if (equali1(r)) break;
     930     1980140 :     if (!signe(remii(N,r)) && !equalii(r,N)) return mkvec2(r, gen_0);
     931     1980140 :     if ((i & 0x1f) == 0) r = gerepileuptoint(av, r);
     932             :   }
     933        1940 :   return gen_1;
     934             : }
     935             : 
     936             : static GEN
     937        1961 : aprcl(GEN N)
     938             : {
     939        1961 :   GEN et, fat, flaglp, faet = NULL; /*-Wall*/
     940             :   long i, j, l, ltab, lfat, lpC;
     941             :   ulong t;
     942             :   Red R;
     943             :   Cache **pC;
     944             : 
     945        1961 :   if (typ(N) != t_INT) pari_err_TYPE("aprcl",N);
     946        1961 :   if (cmpis(N,12) <= 0)
     947          21 :     switch(itos(N))
     948             :     {
     949          14 :       case 2: case 3: case 5: case 7: case 11: return gen_1;
     950           7 :       default: return _res(0,0);
     951             :     }
     952        1940 :   if (Z_issquare(N)) return _res(0,0);
     953        1940 :   t = compute_t(N, &et, &faet);
     954        1940 :   if (DEBUGLEVEL) err_printf("Starting APRCL with t = %ld\n",t);
     955        1940 :   if (cmpii(sqri(et),N) < 0) pari_err_BUG("aprcl: e(t) too small");
     956        1940 :   if (!equali1(gcdii(N,mului(t,et)))) return _res(1,0);
     957             : 
     958        1940 :   R.N = N;
     959        1940 :   R.N2= shifti(N, -1);
     960        1940 :   pC = calcglobs(&R, t, &lpC, &ltab, &fat);
     961        1940 :   if (!pC) return _res(1,0);
     962        1940 :   lfat = lg(fat);
     963        1940 :   flaglp = cgetg(lfat, t_VECSMALL);
     964        1940 :   flaglp[1] = 0;
     965        6769 :   for (i=2; i<lfat; i++)
     966             :   {
     967        4829 :     ulong p = fat[i];
     968        4829 :     flaglp[i] = absequaliu(Fp_powu(N, p-1, sqru(p)), 1);
     969             :   }
     970        1940 :   vecsmall_sort(faet);
     971             : 
     972        1940 :   l = lg(faet);
     973        1940 :   if (DEBUGLEVEL>2) err_printf("Step4: %ld q-values\n", l-1);
     974       22399 :   for (i=l-1; i>0; i--) /* decreasing order: slowest first */
     975             :   {
     976       20459 :     pari_sp av1 = avma;
     977       20459 :     ulong q = faet[i];
     978       20459 :     GEN faq = factoru_pow(q-1), tabg = compute_g(q);
     979       20459 :     GEN P = gel(faq,1), E = gel(faq,2), PE = gel(faq,3);
     980       20459 :     long lfaq = lg(P);
     981       20459 :     pari_sp av2 = avma;
     982       20459 :     if (DEBUGLEVEL>2) err_printf("testing Jacobi sums for q = %ld...",q);
     983       68221 :     for (j=1; j<lfaq; j++, avma = av2)
     984             :     {
     985       47762 :       long p = P[j], e = E[j], pe = PE[j], fl;
     986       47762 :       Cache *C = pC[pe];
     987       47762 :       R.C = C->cyc;
     988       47762 :       if (p >= 3)      fl = step4a(C,&R, q,p,e, tabg);
     989       20459 :       else if (e >= 3) fl = step4b(C,&R, q,e);
     990       19497 :       else if (e == 2) fl = step4c(C,&R, q);
     991       11058 :       else             fl = step4d(C,&R, q);
     992       47762 :       if (fl == -1) return _res(q,p);
     993       47762 :       if (fl == 1) flaglp[ zv_search(fat, p) ] = 0;
     994             :     }
     995       20459 :     if (DEBUGLEVEL>2) err_printf("OK\n");
     996       20459 :     avma = av1;
     997             :   }
     998        1940 :   if (DEBUGLEVEL>2) err_printf("Step5: testing conditions lp\n");
     999        6769 :   for (i=2; i<lfat; i++) /*skip fat[1] = 2*/
    1000             :   {
    1001        4829 :     pari_sp av = avma;
    1002        4829 :     long p = fat[i];
    1003             :     GEN r;
    1004        4829 :     if (flaglp[i] && (r = step5(pC, &R, p, et, ltab, lpC))) return r;
    1005        4829 :     avma = av;
    1006             :   }
    1007        1940 :   if (DEBUGLEVEL>2)
    1008             :   {
    1009           0 :     ulong sc = pC[1]->ctsgt;
    1010           0 :     err_printf("Individual Fermat powerings:\n");
    1011           0 :     for (i=2; i<lpC; i++)
    1012           0 :       if (pC[i]) {
    1013           0 :         err_printf("  %-3ld: %3ld\n", i, pC[i]->ctsgt);
    1014           0 :         sc += pC[i]->ctsgt;
    1015             :       }
    1016           0 :     err_printf("Number of Fermat powerings = %lu\n",sc);
    1017           0 :     err_printf("Step6: testing potential divisors\n");
    1018             :   }
    1019        1940 :   return step6(N, t, et);
    1020             : }
    1021             : 
    1022             : long
    1023        1961 : isprimeAPRCL(GEN N)
    1024             : {
    1025        1961 :   pari_sp av = avma;
    1026        1961 :   GEN res = aprcl(N);
    1027        1961 :   avma = av; return (typ(res) == t_INT);
    1028             : }

Generated by: LCOV version 1.11