Code coverage tests

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

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

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

LCOV - code coverage report
Current view: top level - modules - aprcl.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.8.0 lcov report (development 16791-0d1274a) Lines: 510 583 87.5 %
Date: 2014-09-16 Functions: 45 46 97.8 %
Legend: Lines: hit not hit | Branches: + taken - not taken # not executed Branches: 285 472 60.4 %

           Branch data     Line data    Source code
       1                 :            : /* Copyright (C) 2000  The PARI group.
       2                 :            : 
       3                 :            : This file is part of the PARI/GP package.
       4                 :            : 
       5                 :            : PARI/GP is free software; you can redistribute it and/or modify it under the
       6                 :            : terms of the GNU General Public License as published by the Free Software
       7                 :            : Foundation. It is distributed in the hope that it will be useful, but WITHOUT
       8                 :            : ANY WARRANTY WHATSOEVER.
       9                 :            : 
      10                 :            : Check the License for details. You should have received a copy of it, along
      11                 :            : with the package; see the file 'COPYING'. If not, write to the Free Software
      12                 :            : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
      13                 :            : 
      14                 :            : /*******************************************************************/
      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 associated 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                 :     553370 : makepoldeg1(GEN c, GEN d)
      47                 :            : {
      48                 :            :   GEN z;
      49         [ +  - ]:     553370 :   if (signe(c)) {
      50                 :     553370 :     z = cgetg(4,t_POL); z[1] = evalsigne(1);
      51                 :     553370 :     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                 :     553370 :   return z;
      59                 :            : }
      60                 :            : 
      61                 :            : /* T mod polcyclo(p), assume deg(T) < 2p */
      62                 :            : static GEN
      63                 :     847885 : red_cyclop(GEN T, long p)
      64                 :            : {
      65                 :            :   long i, d;
      66                 :            :   GEN y, z;
      67                 :            : 
      68                 :     847885 :   d = degpol(T) - p; /* < p */
      69         [ -  + ]:     847885 :   if (d <= -2) return T;
      70                 :            : 
      71                 :            :   /* reduce mod (x^p - 1) */
      72                 :     847885 :   y = ZX_mod_Xnm1(T, p);
      73                 :     847885 :   z = y+2;
      74                 :            : 
      75                 :            :   /* reduce mod x^(p-1) + ... + 1 */
      76                 :     847885 :   d = p-1;
      77         [ +  - ]:     847885 :   if (degpol(y) == d)
      78         [ +  + ]:    5268485 :     for (i=0; i<d; i++) gel(z,i) = subii(gel(z,i), gel(z,d));
      79                 :     847885 :   return normalizepol_lg(y, d+2);
      80                 :            : }
      81                 :            : 
      82                 :            : /* x t_VECSMALL, as red_cyclo2n_ip */
      83                 :            : static GEN
      84                 :       5105 : u_red_cyclo2n_ip(GEN x, long n)
      85                 :            : {
      86                 :       5105 :   long i, pow2 = 1L<<(n-1);
      87                 :            :   GEN z;
      88                 :            : 
      89         [ +  + ]:      17435 :   for (i = lg(x)-1; i>pow2; i--) x[i-pow2] -= x[i];
      90         [ +  - ]:       5155 :   for (; i>0; i--)
      91         [ +  + ]:       5155 :     if (x[i]) break;
      92                 :       5105 :   i += 2;
      93                 :       5105 :   z = cgetg(i, t_POL); z[1] = evalsigne(1);
      94         [ +  + ]:      17385 :   for (i--; i>=2; i--) gel(z,i) = stoi(x[i-1]);
      95                 :       5105 :   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                 :     117815 : red_cyclo2n_ip(GEN x, long n)
     100                 :            : {
     101                 :     117815 :   long i, pow2 = 1L<<(n-1);
     102         [ +  + ]:    1234960 :   for (i = lg(x)-1; i>pow2+1; i--)
     103         [ +  + ]:    1117145 :     if (signe(gel(x,i))) gel(x,i-pow2) = subii(gel(x,i-pow2), gel(x,i));
     104                 :     117815 :   return normalizepol_lg(x, i+1);
     105                 :            : }
     106                 :            : static GEN
     107                 :     117675 : red_cyclo2n(GEN x, long n) { return red_cyclo2n_ip(leafcopy(x), n); }
     108                 :            : 
     109                 :            : /* x a non-zero VECSMALL */
     110                 :            : static GEN
     111                 :      15855 : smallpolrev(GEN x)
     112                 :            : {
     113                 :      15855 :   long i,j, lx = lg(x);
     114                 :            :   GEN y;
     115                 :            : 
     116 [ +  - ][ +  + ]:      15865 :   while (lx-- && x[lx]==0) /* empty */;
     117                 :      15855 :   i = lx+2; y = cgetg(i,t_POL);
     118                 :      15855 :   y[1] = evalsigne(1);
     119         [ +  + ]:      89980 :   for (j=2; j<i; j++) gel(y,j) = stoi(x[j-1]);
     120                 :      15855 :   return y;
     121                 :            : }
     122                 :            : 
     123                 :            : /* x polynomial in t_VECSMALL form, T t_POL return x mod T */
     124                 :            : static GEN
     125                 :      15855 : u_red(GEN x, GEN T) {
     126                 :      15855 :   return RgX_rem(smallpolrev(x), T);
     127                 :            : }
     128                 :            : 
     129                 :            : /* special case R->C = polcyclo(2^n) */
     130                 :            : static GEN
     131                 :     117675 : _red_cyclo2n(GEN x, Red *R) {
     132                 :     117675 :   return centermod_i(red_cyclo2n(x, R->n), R->N, R->N2);
     133                 :            : }
     134                 :            : /* special case R->C = polcyclo(p) */
     135                 :            : static GEN
     136                 :     847885 : _red_cyclop(GEN x, Red *R) {
     137                 :     847885 :   return centermod_i(red_cyclop(x, R->n), R->N, R->N2);
     138                 :            : }
     139                 :            : static GEN
     140                 :     131835 : _red(GEN x, Red *R) {
     141                 :     131835 :   return centermod_i(grem(x, R->C), R->N, R->N2);
     142                 :            : }
     143                 :            : static GEN
     144                 :    5014740 : _redsimple(GEN x, Red *R) { return centermodii(x, R->N, R->N2); }
     145                 :            : 
     146                 :            : static GEN
     147                 :    4903700 : sqrmod(GEN x, Red *R) {
     148                 :    4903700 :   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                 :     307190 : sqrmod3(GEN pol, Red *R)
     162                 :            : {
     163                 :            :   GEN a,b,bma,A,B;
     164                 :     307190 :   long lv = lg(pol);
     165                 :            : 
     166         [ -  + ]:     307190 :   if (lv==2) return pol;
     167         [ -  + ]:     307190 :   if (lv==3) return sqrconst(pol, R);
     168                 :     307190 :   a = gel(pol,3);
     169                 :     307190 :   b = gel(pol,2); bma = subii(b,a);
     170                 :     307190 :   A = centermodii(mulii(a,addii(b,bma)), R->N, R->N2);
     171                 :     307190 :   B = centermodii(mulii(bma,addii(a,b)), R->N, R->N2);
     172                 :     307190 :   return makepoldeg1(A,B);
     173                 :            : }
     174                 :            : 
     175                 :            : /* pol^2 mod (x^2+1,N) */
     176                 :            : static GEN
     177                 :     246180 : sqrmod4(GEN pol, Red *R)
     178                 :            : {
     179                 :            :   GEN a,b,A,B;
     180                 :     246180 :   long lv = lg(pol);
     181                 :            : 
     182         [ -  + ]:     246180 :   if (lv==2) return pol;
     183         [ -  + ]:     246180 :   if (lv==3) return sqrconst(pol, R);
     184                 :     246180 :   a = gel(pol,3);
     185                 :     246180 :   b = gel(pol,2);
     186                 :     246180 :   A = centermodii(mulii(a, shifti(b,1)), R->N, R->N2);
     187                 :     246180 :   B = centermodii(mulii(subii(b,a),addii(b,a)), R->N, R->N2);
     188                 :     246180 :   return makepoldeg1(A,B);
     189                 :            : }
     190                 :            : 
     191                 :            : /* pol^2 mod (polcyclo(5),N) */
     192                 :            : static GEN
     193                 :     568035 : sqrmod5(GEN pol, Red *R)
     194                 :            : {
     195                 :            :   GEN c2,b,c,d,A,B,C,D;
     196                 :     568035 :   long lv = lg(pol);
     197                 :            : 
     198         [ -  + ]:     568035 :   if (lv==2) return pol;
     199         [ -  + ]:     568035 :   if (lv==3) return sqrconst(pol, R);
     200                 :     568035 :   c = gel(pol,3); c2 = shifti(c,1);
     201                 :     568035 :   d = gel(pol,2);
     202         [ -  + ]:     568035 :   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                 :     568035 :   b = gel(pol,4);
     212         [ +  + ]:     568035 :   if (lv==5)
     213                 :            :   {
     214                 :         50 :     A = mulii(b, subii(c2,b));
     215                 :         50 :     B = addii(sqri(c), mulii(b, subii(shifti(d,1),b)));
     216                 :         50 :     C = subii(mulii(c2,d), sqri(b));
     217                 :         50 :     D = mulii(subii(d,b), addii(d,b));
     218                 :            :   }
     219                 :            :   else
     220                 :            :   { /* lv == 6 */
     221                 :     567985 :     GEN a = gel(pol,5), a2 = shifti(a,1);
     222                 :            :     /* 2a(d - c) + b(2c - b) */
     223                 :     567985 :     A = addii(mulii(a2, subii(d,c)), mulii(b, subii(c2,b)));
     224                 :            :     /* c(c - 2a) + b(2d - b) */
     225                 :     567985 :     B = addii(mulii(c, subii(c,a2)), mulii(b, subii(shifti(d,1),b)));
     226                 :            :     /* (a-b)(a+b) + 2c(d - a) */
     227                 :     567985 :     C = addii(mulii(subii(a,b),addii(a,b)), mulii(c2,subii(d,a)));
     228                 :            :     /* 2a(b - c) + (d-b)(d+b) */
     229                 :     567985 :     D = addii(mulii(a2, subii(b,c)), mulii(subii(d,b), addii(d,b)));
     230                 :            :   }
     231                 :     568035 :   A = centermodii(A, R->N, R->N2);
     232                 :     568035 :   B = centermodii(B, R->N, R->N2);
     233                 :     568035 :   C = centermodii(C, R->N, R->N2);
     234                 :     568035 :   D = centermodii(D, R->N, R->N2); return mkpoln(4,A,B,C,D);
     235                 :            : }
     236                 :            : 
     237                 :            : static GEN
     238                 :    1190185 : _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                 :      34580 : _powpolmod(Cache *C, GEN jac, Red *R, GEN (*_sqr)(GEN, Red *))
     243                 :            : {
     244                 :      34580 :   const GEN taba = C->aall;
     245                 :      34580 :   const GEN tabt = C->tall;
     246                 :      34580 :   const long efin = lg(taba)-1, lv = R->lv;
     247                 :      34580 :   GEN L, res = jac, pol2 = _sqr(res, R);
     248                 :            :   long f;
     249                 :      34580 :   pari_sp av0 = avma, av;
     250                 :            : 
     251                 :      34580 :   L = cgetg(lv+1, t_VEC); gel(L,1) = res;
     252         [ +  + ]:     304780 :   for (f=2; f<=lv; f++) gel(L,f) = _mul(gel(L,f-1), pol2, R);
     253                 :      34580 :   av = avma;
     254         [ +  + ]:     989145 :   for (f = efin; f >= 1; f--)
     255                 :            :   {
     256                 :     954565 :     GEN t = gel(L, taba[f]);
     257                 :     954565 :     long tf = tabt[f];
     258         [ +  + ]:     954565 :     res = (f==efin)? t: _mul(t, res, R);
     259         [ +  + ]:    6940265 :     while (tf--) {
     260                 :    5985700 :       res = _sqr(res, R);
     261         [ +  + ]:    5985700 :       if (gc_needed(av,1)) {
     262                 :        940 :         res = gerepilecopy(av, res);
     263         [ -  + ]:        940 :         if(DEBUGMEM>1) pari_warn(warnmem,"powpolmod: f = %ld",f);
     264                 :            :       }
     265                 :            :     }
     266                 :            :   }
     267                 :      34580 :   return gerepilecopy(av0, res);
     268                 :            : }
     269                 :            : 
     270                 :            : static GEN
     271                 :       8180 : _powpolmodsimple(Cache *C, Red *R, GEN jac)
     272                 :            : {
     273                 :       8180 :   pari_sp av = avma;
     274                 :       8180 :   GEN w = mulmat_pol(C->matvite, jac);
     275                 :       8180 :   long j, ph = lg(w);
     276                 :            : 
     277                 :       8180 :   R->red = &_redsimple;
     278         [ +  + ]:      30120 :   for (j=1; j<ph; j++)
     279                 :      21940 :     gel(w,j) = _powpolmod(C, centermodii(gel(w,j), R->N, R->N2), R, &sqrmod);
     280                 :       8180 :   w = centermod_i( gmul(C->matinvvite, w), R->N, R->N2 );
     281                 :       8180 :   w = gerepileupto(av, w);
     282                 :       8180 :   return RgV_to_RgX(w, 0);
     283                 :            : }
     284                 :            : 
     285                 :            : static GEN
     286                 :      20820 : powpolmod(Cache *C, Red *R, long p, long k, GEN jac)
     287                 :            : {
     288                 :            :   GEN (*_sqr)(GEN, Red *);
     289                 :            : 
     290         [ -  + ]:      20820 :   if (DEBUGLEVEL>2) C->ctsgt++;
     291         [ +  + ]:      20820 :   if (C->matvite) return _powpolmodsimple(C, R, jac);
     292         [ +  + ]:      12640 :   if (p == 2) /* p = 2 */
     293                 :            :   {
     294         [ +  + ]:       2355 :     if (k == 2) _sqr = &sqrmod4;
     295                 :        100 :     else        _sqr = &sqrmod;
     296                 :       2355 :     R->n = k;
     297                 :       2355 :     R->red = &_red_cyclo2n;
     298         [ +  + ]:      10285 :   } else if (k == 1)
     299                 :            :   {
     300         [ +  + ]:      10185 :     if      (p == 3) _sqr = &sqrmod3;
     301         [ +  + ]:       7195 :     else if (p == 5) _sqr = &sqrmod5;
     302                 :       3040 :     else             _sqr = &sqrmod;
     303                 :      10185 :     R->n = p;
     304                 :      10185 :     R->red = &_red_cyclop;
     305                 :            :   } else {
     306                 :        100 :     R->red = &_red; _sqr = &sqrmod;
     307                 :            :   }
     308                 :      20820 :   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                 :       1160 : compute_e(ulong t, GEN *faet)
     315                 :            : {
     316                 :       1160 :   GEN L, P, D = divisorsu(t);
     317                 :       1160 :   long l = lg(D);
     318                 :            :   ulong k;
     319                 :            : 
     320                 :       1160 :   P = vecsmalltrunc_init(l);
     321                 :       1160 :   L = vecsmalltrunc_init(l);
     322         [ +  + ]:      22400 :   for (k=l-1; k>1; k--) /* k != 1: avoid d = 1 */
     323                 :            :   {
     324                 :      21240 :     ulong d = D[k];
     325         [ +  + ]:      21240 :     if (uisprime(++d))
     326                 :            :     { /* we want q = 1 (mod p) prime, not too large */
     327         [ -  + ]:      11620 :       if (d > 50000000) return gen_0;
     328                 :      11620 :       vecsmalltrunc_append(P, d);
     329                 :      11620 :       vecsmalltrunc_append(L, upowuu(d, 1 + u_lval(t,d)));
     330                 :            :     }
     331                 :            :   }
     332         [ +  - ]:       1160 :   if (faet) *faet = P;
     333                 :       1160 :   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                 :       1160 : compute_t_small(double C)
     362                 :            : {
     363         [ -  + ]:       1160 :   if (C < 17.953) return        6;
     364         [ -  + ]:       1160 :   if (C < 31.996) return       12;
     365         [ -  + ]:       1160 :   if (C < 33.996) return       24;
     366         [ -  + ]:       1160 :   if (C < 54.079) return       36;
     367         [ +  + ]:       1160 :   if (C < 65.325) return       60;
     368         [ -  + ]:        660 :   if (C < 68.457) return       72;
     369         [ -  + ]:        660 :   if (C < 70.783) return      108;
     370         [ +  + ]:        660 :   if (C < 78.039) return      120;
     371         [ -  + ]:        655 :   if (C < 102.41) return      180;
     372         [ -  + ]:        655 :   if (C < 127.50) return      360;
     373         [ +  + ]:        655 :   if (C < 136.68) return      420;
     374         [ -  + ]:          5 :   if (C < 153.43) return      540;
     375         [ -  + ]:          5 :   if (C < 165.66) return      840;
     376         [ -  + ]:          5 :   if (C < 169.17) return     1008;
     377         [ -  + ]:          5 :   if (C < 178.52) return     1080;
     378         [ -  + ]:          5 :   if (C < 192.68) return     1200;
     379         [ -  + ]:          5 :   if (C < 206.34) return     1260;
     380         [ -  + ]:          5 :   if (C < 211.94) return     1620;
     381         [ -  + ]:          5 :   if (C < 222.09) return     1680;
     382         [ -  + ]:          5 :   if (C < 225.11) return     2016;
     383         [ -  + ]:          5 :   if (C < 244.19) return     2160;
     384         [ -  + ]:          5 :   if (C < 270.29) return     2520;
     385         [ -  + ]:          5 :   if (C < 279.50) return     3360;
     386         [ -  + ]:          5 :   if (C < 293.62) return     3780;
     387         [ -  + ]:          5 :   if (C < 346.68) return     5040;
     388         [ -  + ]:          5 :   if (C < 348.70) return     6480;
     389         [ -  + ]:          5 :   if (C < 383.34) return     7560;
     390         [ -  + ]:          5 :   if (C < 396.68) return     8400;
     391         [ -  + ]:          5 :   if (C < 426.04) return    10080;
     392         [ -  + ]:          5 :   if (C < 458.34) return    12600;
     393         [ -  + ]:          5 :   if (C < 527.16) return    15120;
     394         [ -  + ]:          5 :   if (C < 595.38) return    25200;
     395         [ -  + ]:          5 :   if (C < 636.29) return    30240;
     396         [ -  + ]:          5 :   if (C < 672.53) return    42840;
     397         [ -  + ]:          5 :   if (C < 684.90) return    45360;
     398         [ -  + ]:          5 :   if (C < 708.78) return    55440;
     399         [ -  + ]:          5 :   if (C < 771.30) return    60480;
     400         [ -  + ]:          5 :   if (C < 775.86) return    75600;
     401         [ -  + ]:          5 :   if (C < 859.62) return    85680;
     402         [ -  + ]:          5 :   if (C < 893.16) return   100800;
     403         [ -  + ]:          5 :   if (C < 912.27) return   110880;
     404         [ -  + ]:          5 :   if (C < 966.13) return   128520;
     405         [ +  - ]:          5 :   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                 :       1160 :   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                 :       1160 : 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                 :       1160 :   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         [ +  - ]:       1160 :   if (C < 3514.6)
     452                 :            :   {
     453                 :       1160 :     t = compute_t_small(C);
     454                 :       1160 :     *e = compute_e(t, faet);
     455                 :       1160 :     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                 :       1160 :   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                 :      16590 : computetabdl(ulong q)
     472                 :            : {
     473                 :      16590 :   ulong g, a, i, qs2 = q>>1; /* (q-1)/2 */
     474                 :      16590 :   GEN T = cgetg(qs2+2,t_VECSMALL);
     475                 :            : 
     476                 :      16590 :   g = pgener_Fl(q); a = 1;
     477         [ +  + ]:    2129200 :   for (i=1; i < qs2; i++) /* g^((q-1)/2) = -1 */
     478                 :            :   {
     479                 :    2112610 :     a = Fl_mul(g,a,q);
     480         [ +  + ]:    2112610 :     if (a > qs2) T[q-a] = i+qs2; else T[a] = i;
     481                 :            :   }
     482                 :      16590 :   T[qs2+1] = T[qs2] + qs2;
     483                 :      16590 :   T[1] = 0; return T;
     484                 :            : }
     485                 :            : 
     486                 :            : /* Return T: T[x] = dl of x(1-x) */
     487                 :            : static GEN
     488                 :      11625 : compute_g(ulong q)
     489                 :            : {
     490                 :      11625 :   const ulong qs2 = q>>1; /* (q-1)/2 */
     491                 :            :   ulong x, a;
     492                 :      11625 :   GEN T = computetabdl(q); /* updated in place to save on memory */
     493                 :      11625 :   a = 0; /* dl[1] */
     494         [ +  + ]:    1161055 :   for (x=2; x<=qs2+1; x++)
     495                 :            :   { /* a = dl(x) */
     496                 :    1149430 :     ulong b = T[x]; /* = dl(x) */
     497                 :    1149430 :     T[x] = b + a + qs2; /* dl(x) + dl(x-1) + dl(-1) */
     498                 :    1149430 :     a = b;
     499                 :            :   }
     500                 :      11625 :   return T;
     501                 :            : }
     502                 :            : 
     503                 :            : /* p odd prime */
     504                 :            : static GEN
     505                 :      15855 : get_jac(Cache *C, ulong q, long pk, GEN tabg)
     506                 :            : {
     507                 :      15855 :   ulong x, qs2 = q>>1; /* (q-1)/2 */
     508                 :      15855 :   GEN vpk = zero_zv(pk);
     509                 :            : 
     510         [ +  + ]:    3879385 :   for (x=2; x<=qs2; x++) vpk[ tabg[x]%pk + 1 ] += 2;
     511                 :      15855 :   vpk[ tabg[x]%pk + 1 ]++; /* x = (q+1)/2 */
     512                 :      15855 :   return u_red(vpk, C->cyc);
     513                 :            : }
     514                 :            : 
     515                 :            : /* p = 2 */
     516                 :            : static GEN
     517                 :       4965 : get_jac2(GEN N, ulong q, long k, GEN *j2q, GEN *j3q)
     518                 :            : {
     519                 :       4965 :   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                 :       4965 :   pk = 1UL << k;;
     525                 :       4965 :   vpk = zero_zv(pk);
     526                 :            : 
     527                 :       4965 :   qs2 = q>>1; /* (q-1)/2 */
     528                 :            : 
     529         [ +  + ]:     979770 :   for (x=2; x<=qs2; x++) vpk[ (T[x]+T[x-1]+qs2)%pk + 1 ] += 2;
     530                 :       4965 :   vpk[ (T[x]+T[x-1]+qs2)%pk + 1 ]++;
     531                 :       4965 :   jpq = u_red_cyclo2n_ip(vpk, k);
     532                 :            : 
     533         [ +  + ]:       4965 :   if (k == 2) return jpq;
     534                 :            : 
     535         [ -  + ]:        140 :   if (mod8(N) >= 5)
     536                 :            :   {
     537                 :          0 :     GEN v8 = cgetg(9,t_VECSMALL);
     538         [ #  # ]:          0 :     for (x=1; x<=8; x++) v8[x] = 0;
     539         [ #  # ]:          0 :     for (x=2; x<=qs2; x++) v8[ ((3*T[x]+T[x-1]+qs2)&7) + 1 ]++;
     540         [ #  # ]:          0 :     for (   ; x<=q-1; x++) v8[ ((3*T[q-x]+T[q-x+1]-3*qs2)&7) + 1 ]++;
     541                 :          0 :     *j2q = RgX_inflate(ZX_sqr(u_red_cyclo2n_ip(v8,3)), pk>>3);
     542                 :          0 :     *j2q = red_cyclo2n_ip(*j2q, k);
     543                 :            :   }
     544         [ +  + ]:       2820 :   for (i=1; i<=pk; i++) vpk[i] = 0;
     545         [ +  + ]:     738380 :   for (x=2; x<=qs2; x++) vpk[ (2*T[x]+T[x-1]+qs2)%pk + 1 ]++;
     546         [ +  + ]:     738520 :   for (   ; x<=q-1; x++) vpk[ (2*T[q-x]+T[q-x+1]-2*qs2)%pk + 1 ]++;
     547                 :        140 :   *j3q = ZX_mul(jpq, u_red_cyclo2n_ip(vpk,k));
     548                 :        140 :   *j3q = red_cyclo2n_ip(*j3q, k);
     549                 :       4965 :   return jpq;
     550                 :            : }
     551                 :            : 
     552                 :            : /* N = 1 mod p^k, return an elt of order p^k in (Z/N)^* */
     553                 :            : static GEN
     554                 :       1615 : finda(Cache *Cp, GEN N, long pk, long p)
     555                 :            : {
     556                 :            :   GEN a, pv;
     557 [ +  - ][ +  + ]:       1615 :   if (Cp && Cp->avite) {
     558                 :          5 :     a  = Cp->avite;
     559                 :          5 :     pv = Cp->pkvite;
     560                 :            :   }
     561                 :            :   else
     562                 :            :   {
     563                 :            :     GEN ph, b, q;
     564                 :       1610 :     ulong u = 2;
     565                 :       1610 :     long v = Z_lvalrem(addis(N,-1), p, &q);
     566                 :       1610 :     ph = powuu(p, v-1); pv = muliu(ph, p); /* N - 1 = p^v q */
     567         [ +  + ]:       1610 :     if (p > 2)
     568                 :            :     {
     569                 :        585 :       for (;;u++)
     570                 :            :       {
     571                 :       1580 :         a = Fp_pow(utoipos(u), q, N);
     572                 :       1580 :         b = Fp_pow(a, ph, N);
     573         [ +  + ]:       1580 :         if (!gequal1(b)) break;
     574                 :       1580 :       }
     575                 :            :     }
     576                 :            :     else
     577                 :            :     {
     578         [ +  + ]:       2220 :       while (krosi(u,N) >= 0) u++;
     579                 :        615 :       a = Fp_pow(utoipos(u), q, N);
     580                 :        615 :       b = Fp_pow(a, ph, N);
     581                 :            :     }
     582                 :            :     /* checking b^p = 1 mod N done economically in caller */
     583                 :       1610 :     b = gcdii(addis(b,-1), N);
     584         [ -  + ]:       1610 :     if (!gequal1(b)) return NULL;
     585                 :            : 
     586         [ +  - ]:       1610 :     if (Cp) {
     587                 :       1610 :       Cp->avite  = a; /* a has order p^v */
     588                 :       1610 :       Cp->pkvite = pv;
     589                 :            :     }
     590                 :            :   }
     591                 :       1615 :   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                 :       5330 : 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                 :       5330 :   C->cyc = polcyclo(pk,0);
     604                 :            : 
     605         [ +  + ]:       5330 :   if (p > 2)
     606                 :            :   {
     607                 :       2990 :     long LE = pk - pk/p + 1;
     608                 :       2990 :     GEN E = cgetg(LE, t_VECSMALL), eta = cgetg(pk+1,t_VEC);
     609         [ +  + ]:      14020 :     for (i=1,j=0; i<pk; i++)
     610         [ +  + ]:      11030 :       if (i%p) E[++j] = i;
     611                 :       2990 :     C->E = E;
     612                 :            : 
     613         [ +  + ]:      17010 :     for (i=1; i<=pk; i++)
     614                 :            :     {
     615                 :      14020 :       GEN z = FpX_rem(monomial(gen_1, i-1, 0), C->cyc, R->N);
     616                 :      14020 :       gel(eta,i) = FpX_center(z, R->N, R->N2);
     617                 :            :     }
     618                 :       2990 :     C->eta = eta;
     619                 :            :   }
     620         [ +  + ]:       2340 :   else if (pk >= 8)
     621                 :            :   {
     622                 :         20 :     long LE = (pk>>2) + 1;
     623                 :         20 :     GEN E = cgetg(LE, t_VECSMALL);
     624         [ +  + ]:        320 :     for (i=1,j=0; i<pk; i++)
     625 [ +  + ][ +  + ]:        300 :       if ((i%8)==1 || (i%8)==3) E[++j] = i;
     626                 :         20 :     C->E = E;
     627                 :            :   }
     628                 :            : 
     629 [ +  + ][ +  + ]:       5330 :   if (pk > 2 && umodiu(R->N,pk) == 1)
     630                 :            :   {
     631                 :       1615 :     GEN vpa, p1, p2, p3, a2 = NULL, a = finda(Cp, R->N, pk, p);
     632                 :            :     long jj, ph;
     633                 :            : 
     634         [ -  + ]:       1615 :     if (!a) return 0;
     635                 :       1615 :     ph = pk - pk/p;
     636                 :       1615 :     vpa = cgetg(ph+1,t_COL); gel(vpa,1) = a;
     637         [ +  + ]:       1615 :     if (pk > p) a2 = centermodii(sqri(a), R->N, R->N2);
     638                 :       1615 :     jj = 1;
     639         [ +  + ]:       4810 :     for (i=2; i<pk; i++) /* vpa = { a^i, (i,p) = 1 } */
     640         [ +  + ]:       3195 :       if (i%p) {
     641         [ +  + ]:       2565 :         GEN z = mulii((i%p==1) ? a2 : a, gel(vpa,jj));
     642                 :       2565 :         gel(vpa,++jj) = centermodii(z , R->N, R->N2);
     643                 :            :       }
     644         [ -  + ]:       1615 :     if (!gequal1( centermodii( mulii(a, gel(vpa,ph)), R->N, R->N2) )) return 0;
     645                 :       1615 :     p1 = cgetg(ph+1,t_MAT);
     646                 :       1615 :     p2 = cgetg(ph+1,t_COL); gel(p1,1) = p2;
     647         [ +  + ]:       5795 :     for (i=1; i<=ph; i++) gel(p2,i) = gen_1;
     648                 :       1615 :     j = 2; gel(p1,j) = vpa; p3 = vpa;
     649         [ +  + ]:       2565 :     for (j++; j <= ph; j++)
     650                 :            :     {
     651                 :        950 :       p2 = cgetg(ph+1,t_COL); gel(p1,j) = p2;
     652         [ +  + ]:       5710 :       for (i=1; i<=ph; i++)
     653                 :       4760 :         gel(p2,i) = centermodii(mulii(gel(vpa,i),gel(p3,i)), R->N, R->N2);
     654                 :        950 :       p3 = p2;
     655                 :            :     }
     656                 :       1615 :     C->matvite = p1;
     657                 :       1615 :     C->matinvvite = FpM_inv(p1, R->N);
     658                 :            :   }
     659                 :            : 
     660                 :       5330 :   tabt = cgetg(ltab+1, t_VECSMALL);
     661                 :       5330 :   taba = cgetg(ltab+1, t_VECSMALL);
     662                 :       5330 :   av = avma; m = divis(R->N, pk);
     663 [ +  - ][ +  + ]:      89495 :   for (e=1; e<=ltab && signe(m); e++)
     664                 :            :   {
     665                 :      84165 :     long s = vali(m); m = shifti(m,-s);
     666         [ +  + ]:      84165 :     tabt[e] = e==1? s: s + R->k;
     667         [ +  - ]:      84165 :     taba[e] = signe(m)? ((mod2BIL(m) & R->mask)+1)>>1: 0;
     668                 :      84165 :     m = shifti(m, -R->k);
     669                 :            :   }
     670                 :       5330 :   setlg(taba, e); C->aall = taba;
     671                 :       5330 :   setlg(tabt, e); C->tall = tabt;
     672                 :       5330 :   avma = av; return 1;
     673                 :            : }
     674                 :            : 
     675                 :            : static Cache *
     676                 :       6485 : alloc_cache(void)
     677                 :            : {
     678                 :       6485 :   Cache *C = (Cache*)stack_malloc(sizeof(Cache));
     679                 :       6485 :   C->matvite = NULL;
     680                 :       6485 :   C->avite   = NULL;
     681                 :       6485 :   C->ctsgt = 0; return C;
     682                 :            : }
     683                 :            : 
     684                 :            : static Cache **
     685                 :       1160 : 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                 :       1160 :   b = expi(R->N)+1;
     692         [ +  + ]:       1825 :   k = 3; while (((k+1)*(k+2) << (k-1)) < b) k++;
     693                 :       1160 :   *pltab = (b/k)+2;
     694                 :       1160 :   R->k  = k;
     695                 :       1160 :   R->lv = 1L << (k-1);
     696                 :       1160 :   R->mask = (1UL << k) - 1;
     697                 :            : 
     698                 :       1160 :   fat = factoru_pow(t);
     699                 :       1160 :   P = gel(fat,1);
     700                 :       1160 :   E = gel(fat,2);
     701                 :       1160 :   PE= gel(fat,3);
     702                 :       1160 :   *plpC = lv = vecsmall_max(PE); /* max(p^e, p^e | t) */
     703                 :       1160 :   pC = (Cache**)stack_malloc((lv+1)*sizeof(Cache*));
     704                 :       1160 :   pC[1] = alloc_cache(); /* to be used as temp in step5() */
     705         [ +  + ]:       7250 :   for (i = 2; i <= lv; i++) pC[i] = NULL;
     706         [ +  + ]:       5300 :   for (i=1; i<lg(P); i++)
     707                 :            :   {
     708                 :       4140 :     long pk, p = P[i], e = E[i];
     709                 :       4140 :     pk = p;
     710         [ +  + ]:       9465 :     for (k=1; k<=e; k++, pk*=p)
     711                 :            :     {
     712                 :       5325 :       pC[pk] = alloc_cache();
     713         [ -  + ]:       5325 :       if (!filltabs(pC[pk], pC[p], R, p,pk, *pltab)) return NULL;
     714                 :            :     }
     715                 :            :   }
     716                 :       1160 :   *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                 :      87525 : aut(long pk, GEN z, long a)
     723                 :            : {
     724                 :            :   GEN v;
     725                 :      87525 :   long b, i, dz = degpol(z);
     726 [ +  + ][ -  + ]:      87525 :   if (a == 1 || dz < 0) return z;
     727                 :      71530 :   v = cgetg(pk+2,t_POL);
     728                 :      71530 :   v[1] = evalvarn(0);
     729                 :      71530 :   b = 0;
     730                 :      71530 :   gel(v,2) = gel(z,2); /* i = 0 */
     731         [ +  + ]:     434820 :   for (i = 1; i < pk; i++)
     732                 :            :   {
     733         [ +  + ]:     363290 :     b += a; if (b > pk) b -= pk; /* b = (a*i) % pk */
     734         [ +  + ]:     363290 :     gel(v,i+2) = b > dz? gen_0: gel(z,b+2);
     735                 :            :   }
     736                 :      87525 :   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                 :      15995 : autvec_TH(long pk, GEN z, GEN v, GEN C)
     742                 :            : {
     743                 :      15995 :   long i, lv = lg(v);
     744                 :      15995 :   GEN s = pol_1(varn(C));
     745         [ +  + ]:      74745 :   for (i=1; i<lv; i++)
     746                 :            :   {
     747                 :      58750 :     long y = v[i];
     748         [ +  - ]:      58750 :     if (y) s = RgXQ_mul(s, RgXQ_powu(aut(pk,z, y), y, C), C);
     749                 :            :   }
     750                 :      15995 :   return s;
     751                 :            : }
     752                 :            : 
     753                 :            : static GEN
     754                 :      15995 : autvec_AL(long pk, GEN z, GEN v, Red *R)
     755                 :            : {
     756                 :      15995 :   const long r = umodiu(R->N, pk);
     757                 :      15995 :   GEN s = pol_1(varn(R->C));
     758                 :      15995 :   long i, lv = lg(v);
     759         [ +  + ]:      74745 :   for (i=1; i<lv; i++)
     760                 :            :   {
     761                 :      58750 :     long y = (r*v[i]) / pk;
     762         [ +  + ]:      58750 :     if (y) s = RgXQ_mul(s, RgXQ_powu(aut(pk,z, v[i]), y, R->C), R->C);
     763                 :            :   }
     764                 :      15995 :   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                 :      15855 : look_eta(GEN eta, long pk, GEN z)
     770                 :            : {
     771                 :            :   long i;
     772         [ +  - ]:      44790 :   for (i=1; i<=pk; i++)
     773         [ +  + ]:      44790 :     if (ZX_equal(z, gel(eta,i))) return i-1;
     774                 :      15855 :   return -1;
     775                 :            : }
     776                 :            : /* same pk = 2^k */
     777                 :            : static long
     778                 :       4965 : look_eta2(long k, GEN z)
     779                 :            : {
     780                 :            :   long d, s;
     781         [ -  + ]:       4965 :   if (typ(z) != t_POL) d = 0; /* t_INT */
     782                 :            :   else
     783                 :            :   {
     784         [ -  + ]:       4965 :     if (!RgX_is_monomial(z)) return -1;
     785                 :       4965 :     d = degpol(z);
     786                 :       4965 :     z = gel(z,d+2); /* leading term */
     787                 :            :   }
     788                 :       4965 :   s = signe(z);
     789 [ +  - ][ -  + ]:       4965 :   if (!s || !is_pm1(z)) return -1;
     790         [ +  + ]:       4965 :   return (s > 0)? d: d + (1L<<(k-1));
     791                 :            : }
     792                 :            : 
     793                 :            : static long
     794                 :      15855 : step4a(Cache *C, Red *R, ulong q, long p, long k, GEN tabg)
     795                 :            : {
     796                 :      15855 :   const long pk = upowuu(p,k);
     797                 :            :   long ind;
     798                 :            :   GEN jpq, s1, s2, s3;
     799                 :            : 
     800         [ +  + ]:      15855 :   if (!tabg) tabg = compute_g(q);
     801                 :      15855 :   jpq = get_jac(C, q, pk, tabg);
     802                 :      15855 :   s1 = autvec_TH(pk, jpq, C->E, C->cyc);
     803                 :      15855 :   s2 = powpolmod(C,R, p,k, s1);
     804                 :      15855 :   s3 = autvec_AL(pk, jpq, C->E, R);
     805                 :      15855 :   s3 = _red(gmul(s3,s2), R);
     806                 :            : 
     807                 :      15855 :   ind = look_eta(C->eta, pk, s3);
     808         [ -  + ]:      15855 :   if (ind < 0) return -1;
     809                 :      15855 :   return (ind%p) != 0;
     810                 :            : }
     811                 :            : 
     812                 :            : /* x == -1 mod N ? */
     813                 :            : static long
     814                 :       5560 : is_m1(GEN x, GEN N)
     815                 :            : {
     816                 :       5560 :   return equalii(addis(x,1), N);
     817                 :            : }
     818                 :            : 
     819                 :            : /* p=2, k>=3 */
     820                 :            : static long
     821                 :        140 : step4b(Cache *C, Red *R, ulong q, long k)
     822                 :            : {
     823                 :        140 :   const long pk = 1L << k;
     824                 :            :   long ind;
     825                 :        140 :   GEN s1, s2, s3, j2q = NULL, j3q = NULL;
     826                 :            : 
     827                 :        140 :   (void)get_jac2(R->N,q,k, &j2q,&j3q);
     828                 :            : 
     829                 :        140 :   s1 = autvec_TH(pk, j3q, C->E, C->cyc);
     830                 :        140 :   s2 = powpolmod(C,R, 2,k, s1);
     831                 :        140 :   s3 = autvec_AL(pk, j3q, C->E, R);
     832                 :        140 :   s3 = _red(gmul(s3,s2), R);
     833         [ -  + ]:        140 :   if (j2q) s3 = _red(gmul(j2q, s3), R);
     834                 :            : 
     835                 :        140 :   ind = look_eta2(k, s3);
     836         [ -  + ]:        140 :   if (ind < 0) return -1;
     837         [ +  + ]:        140 :   if ((ind&1)==0) return 0;
     838         [ -  + ]:         70 :   if (DEBUGLEVEL>2) C->ctsgt++;
     839                 :         70 :   s3 = Fp_pow(utoipos(q), R->N2, R->N);
     840                 :        140 :   return is_m1(s3, R->N);
     841                 :            : }
     842                 :            : 
     843                 :            : /* p=2, k=2 */
     844                 :            : static long
     845                 :       4825 : step4c(Cache *C, Red *R, ulong q)
     846                 :            : {
     847                 :            :   long ind;
     848                 :       4825 :   GEN s0,s1,s3, jpq = get_jac2(R->N,q,2, NULL,NULL);
     849                 :            : 
     850                 :       4825 :   s0 = sqrmod4(jpq, R);
     851                 :       4825 :   s1 = gmulsg(q,s0);
     852                 :       4825 :   s3 = powpolmod(C,R, 2,2, s1);
     853         [ +  + ]:       4825 :   if (mod4(R->N) == 3) s3 = _red(gmul(s0,s3), R);
     854                 :            : 
     855                 :       4825 :   ind = look_eta2(2, s3);
     856         [ -  + ]:       4825 :   if (ind < 0) return -1;
     857         [ +  + ]:       4825 :   if ((ind&1)==0) return 0;
     858         [ -  + ]:       2275 :   if (DEBUGLEVEL>2) C->ctsgt++;
     859                 :       2275 :   s3 = Fp_pow(utoipos(q), R->N2, R->N);
     860                 :       4825 :   return is_m1(s3, R->N);
     861                 :            : }
     862                 :            : 
     863                 :            : /* p=2, k=1 */
     864                 :            : static long
     865                 :       6655 : step4d(Cache *C, Red *R, ulong q)
     866                 :            : {
     867                 :       6655 :   GEN s1 = Fp_pow(utoipos(q), R->N2, R->N);
     868         [ -  + ]:       6655 :   if (DEBUGLEVEL>2) C->ctsgt++;
     869         [ +  + ]:       6655 :   if (is_pm1(s1)) return 0;
     870         [ +  - ]:       3215 :   if (is_m1(s1, R->N)) return (mod4(R->N) == 1);
     871                 :       6655 :   return -1;
     872                 :            : }
     873                 :            : 
     874                 :            : static GEN
     875         [ -  + ]:          5 : _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                 :          5 : step5(Cache **pC, Red *R, long p, GEN et, ulong ltab, long lpC)
     880                 :            : {
     881                 :            :   pari_sp av;
     882                 :            :   ulong q;
     883                 :          5 :   long pk, k, fl = -1;
     884                 :            :   Cache *C, *Cp;
     885                 :            :   forprime_t T;
     886                 :            : 
     887                 :          5 :   (void)u_forprime_arith_init(&T, 3, ULONG_MAX, 1,p);
     888         [ +  - ]:         15 :   while( (q = u_forprime_next(&T)) )
     889                 :            :   { /* q = 1 (mod p) */
     890         [ +  + ]:         15 :     if (umodiu(et,q) == 0) continue;
     891         [ -  + ]:          5 :     if (umodiu(R->N,q) == 0) return _res(1,p);
     892                 :          5 :     k = u_lval(q-1, p);
     893                 :          5 :     pk = upowuu(p,k);
     894 [ -  + ][ #  # ]:          5 :     if (pk <= lpC && pC[pk]) {
     895                 :          0 :       C = pC[pk];
     896                 :          0 :       Cp = pC[p];
     897                 :            :     } else {
     898                 :          5 :       C = pC[1];
     899                 :          5 :       Cp = NULL;
     900                 :          5 :       C->matvite = NULL; /* re-init */
     901                 :            :     }
     902                 :            : 
     903                 :          5 :     av = avma;
     904         [ -  + ]:          5 :     if (!filltabs(C, Cp, R, p, pk, ltab)) return _res(1,0);
     905                 :          5 :     R->C = C->cyc;
     906         [ +  - ]:          5 :     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         [ -  + ]:          5 :     if (fl == -1) return _res(q,p);
     911         [ +  - ]:          5 :     if (fl == 1) return NULL; /*OK*/
     912                 :          0 :     avma = av;
     913                 :            :   }
     914                 :          0 :   pari_err_BUG("aprcl test fails! This is highly improbable");
     915                 :          5 :   return NULL;
     916                 :            : }
     917                 :            : 
     918                 :            : static GEN
     919                 :       1160 : step6(GEN N, ulong t, GEN et)
     920                 :            : {
     921                 :       1160 :   GEN r, N1 = remii(N, et);
     922                 :            :   ulong i;
     923                 :       1160 :   pari_sp av = avma;
     924                 :            : 
     925                 :       1160 :   r = gen_1;
     926         [ +  + ]:     953600 :   for (i=1; i<t; i++)
     927                 :            :   {
     928                 :     952485 :     r = remii(mulii(r,N1), et);
     929         [ +  + ]:     952485 :     if (equali1(r)) break;
     930 [ -  + ][ #  # ]:     952440 :     if (!signe(remii(N,r)) && !equalii(r,N)) return mkvec2(r, gen_0);
     931         [ +  + ]:     952440 :     if ((i & 0x1f) == 0) r = gerepileuptoint(av, r);
     932                 :            :   }
     933                 :       1160 :   return gen_1;
     934                 :            : }
     935                 :            : 
     936                 :            : static GEN
     937                 :       1175 : aprcl(GEN N)
     938                 :            : {
     939                 :       1175 :   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         [ -  + ]:       1175 :   if (typ(N) != t_INT) pari_err_TYPE("aprcl",N);
     946         [ +  + ]:       1175 :   if (cmpis(N,12) <= 0)
     947         [ +  + ]:         15 :     switch(itos(N))
     948                 :            :     {
     949                 :         10 :       case 2: case 3: case 5: case 7: case 11: return gen_1;
     950                 :          5 :       default: return _res(0,0);
     951                 :            :     }
     952         [ -  + ]:       1160 :   if (Z_issquare(N)) return _res(0,0);
     953                 :       1160 :   t = compute_t(N, &et, &faet);
     954         [ -  + ]:       1160 :   if (DEBUGLEVEL) err_printf("Starting APRCL with t = %ld\n",t);
     955         [ -  + ]:       1160 :   if (cmpii(sqri(et),N) < 0) pari_err_BUG("aprcl: e(t) too small");
     956         [ -  + ]:       1160 :   if (!equali1(gcdii(N,mului(t,et)))) return _res(1,0);
     957                 :            : 
     958                 :       1160 :   R.N = N;
     959                 :       1160 :   R.N2= shifti(N, -1);
     960                 :       1160 :   pC = calcglobs(&R, t, &lpC, &ltab, &fat);
     961         [ -  + ]:       1160 :   if (!pC) return _res(1,0);
     962                 :       1160 :   lfat = lg(fat);
     963                 :       1160 :   flaglp = cgetg(lfat, t_VECSMALL);
     964                 :       1160 :   flaglp[1] = 0;
     965         [ +  + ]:       4140 :   for (i=2; i<lfat; i++)
     966                 :            :   {
     967                 :       2980 :     ulong p = fat[i];
     968                 :       2980 :     flaglp[i] = equaliu(Fp_powu(N, p-1, sqru(p)), 1);
     969                 :            :   }
     970                 :       1160 :   vecsmall_sort(faet);
     971                 :            : 
     972                 :       1160 :   l = lg(faet);
     973         [ -  + ]:       1160 :   if (DEBUGLEVEL>2) err_printf("Step4: %ld q-values\n", l-1);
     974         [ +  + ]:      12780 :   for (i=l-1; i>0; i--) /* decreasing order: slowest first */
     975                 :            :   {
     976                 :      11620 :     pari_sp av1 = avma;
     977                 :      11620 :     ulong q = faet[i];
     978                 :      11620 :     GEN faq = factoru_pow(q-1), tabg = compute_g(q);
     979                 :      11620 :     GEN P = gel(faq,1), E = gel(faq,2), PE = gel(faq,3);
     980                 :      11620 :     long lfaq = lg(P);
     981                 :      11620 :     pari_sp av2 = avma;
     982         [ -  + ]:      11620 :     if (DEBUGLEVEL>2) err_printf("testing Jacobi sums for q = %ld...",q);
     983         [ +  + ]:      39090 :     for (j=1; j<lfaq; j++, avma = av2)
     984                 :            :     {
     985                 :      27470 :       long p = P[j], e = E[j], pe = PE[j], fl;
     986                 :      27470 :       Cache *C = pC[pe];
     987                 :      27470 :       R.C = C->cyc;
     988         [ +  + ]:      27470 :       if (p >= 3)      fl = step4a(C,&R, q,p,e, tabg);
     989         [ +  + ]:      11620 :       else if (e >= 3) fl = step4b(C,&R, q,e);
     990         [ +  + ]:      11480 :       else if (e == 2) fl = step4c(C,&R, q);
     991                 :       6655 :       else             fl = step4d(C,&R, q);
     992         [ -  + ]:      27470 :       if (fl == -1) return _res(q,p);
     993         [ +  + ]:      27470 :       if (fl == 1) flaglp[ zv_search(fat, p) ] = 0;
     994                 :            :     }
     995         [ -  + ]:      11620 :     if (DEBUGLEVEL>2) err_printf("OK\n");
     996                 :      11620 :     avma = av1;
     997                 :            :   }
     998         [ -  + ]:       1160 :   if (DEBUGLEVEL>2) err_printf("Step5: testing conditions lp\n");
     999         [ +  + ]:       4140 :   for (i=2; i<lfat; i++) /*skip fat[1] = 2*/
    1000                 :            :   {
    1001                 :       2980 :     pari_sp av = avma;
    1002                 :       2980 :     long p = fat[i];
    1003                 :            :     GEN r;
    1004 [ +  + ][ -  + ]:       2980 :     if (flaglp[i] && (r = step5(pC, &R, p, et, ltab, lpC))) return r;
    1005                 :       2980 :     avma = av;
    1006                 :            :   }
    1007         [ -  + ]:       1160 :   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                 :       1175 :   return step6(N, t, et);
    1020                 :            : }
    1021                 :            : 
    1022                 :            : long
    1023                 :       1175 : isprimeAPRCL(GEN N)
    1024                 :            : {
    1025                 :       1175 :   pari_sp av = avma;
    1026                 :       1175 :   GEN res = aprcl(N);
    1027                 :       1175 :   avma = av; return (typ(res) == t_INT);
    1028                 :            : }

Generated by: LCOV version 1.9