Code coverage tests

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

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

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

LCOV - code coverage report
Current view: top level - basemath - aprcl.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.10.0 lcov report (development 21353-12523aa) Lines: 610 721 84.6 %
Date: 2017-11-24 06:20:58 Functions: 51 52 98.1 %
Legend: Lines: hit not hit

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

Generated by: LCOV version 1.11