Code coverage tests

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

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

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

LCOV - code coverage report
Current view: top level - basemath - ecpp.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.12.0 lcov report (development 23010-740a36cf0) Lines: 615 680 90.4 %
Date: 2018-09-21 05:37:29 Functions: 81 86 94.2 %
Legend: Lines: hit not hit

          Line data    Source code
       1             : /* Copyright (C) 2014 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             : #include "pari.h"
      15             : #include "paripriv.h"
      16             : 
      17             : #define dbg_mode() if (DEBUGLEVEL >= 2)
      18             : 
      19             : #define ANSI_RED     "\x1b[31m"
      20             : #define ANSI_GREEN   "\x1b[32m"
      21             : #define ANSI_YELLOW  "\x1b[33m"
      22             : #define ANSI_BLUE    "\x1b[34m"
      23             : #define ANSI_MAGENTA "\x1b[35m"
      24             : #define ANSI_CYAN    "\x1b[36m"
      25             : #define ANSI_WHITE   "\x1b[37m"
      26             : 
      27             : #define ANSI_BRIGHT_RED     "\x1b[31;1m"
      28             : #define ANSI_BRIGHT_GREEN   "\x1b[32;1m"
      29             : #define ANSI_BRIGHT_YELLOW  "\x1b[33;1m"
      30             : #define ANSI_BRIGHT_BLUE    "\x1b[34;1m"
      31             : #define ANSI_BRIGHT_MAGENTA "\x1b[35;1m"
      32             : #define ANSI_BRIGHT_CYAN    "\x1b[36;1m"
      33             : #define ANSI_BRIGHT_WHITE   "\x1b[37;1m"
      34             : 
      35             : #define ANSI_RESET   "\x1b[0m"
      36             : 
      37             : /* THINGS THAT DON'T BELONG */
      38             : 
      39             : /* Assume that x is a vector such that
      40             :    f(x) = 1 for x <= k
      41             :    f(x) = 0 otherwise.
      42             :    Return k.
      43             : */
      44             : static long
      45        6818 : zv_binsearch0(void *E, long (*f)(void* E, GEN x), GEN x)
      46             : {
      47        6818 :   long lo = 1;
      48        6818 :   long hi = lg(x)-1;
      49       29946 :   while (lo < hi) {
      50       16310 :     long mi = lo + (hi - lo)/2 + 1;
      51       16310 :     if (f(E,gel(x,mi))) lo = mi;
      52        9450 :     else hi = mi - 1;
      53             :   }
      54        6818 :   if (f(E,gel(x,lo))) return lo;
      55        3388 :   return 0;
      56             : }
      57             : 
      58             : INLINE long
      59           0 : timer_record(GEN* X0, const char* Xx, pari_timer* ti)
      60             : {
      61           0 :   long t = timer_delay(ti), i = Xx[0]-'A'+1, j = Xx[1]-'1'+1;
      62           0 :   umael3(*X0, 1, i, j) += t;
      63           0 :   umael3(*X0, 2, i, j) ++; return t;
      64             : }
      65             : 
      66             : INLINE long
      67        4624 : FpJ_is_inf(GEN P) { return signe(gel(P, 3)) == 0; }
      68             : 
      69             : /*****************************************************************/
      70             : 
      71             : /* D < 0 fundamental: return the number of units in Q(sqrt(D)) */
      72             : INLINE long
      73    57465009 : D_get_wD(long D)
      74             : {
      75    57465009 :   if (D == -4) return 4;
      76    57463693 :   if (D == -3) return 6;
      77    57462132 :   return 2;
      78             : }
      79             : 
      80             : /* Dinfo contains information related to the discirminant
      81             :  *    D: the discriminant (D < 0)
      82             :  *    h: the class number associated to D
      83             :  *   bi: the "best invariant" for computing polclass(D, bi)
      84             :  *   pd: the degree of polclass; equal to h except when bi is a double
      85             :  *       eta-quotient w_p,q with p|D and q|D, where pd = h/2.
      86             :  * Dfac: the prime factorization of D; we have D = q0 q1* q2* ... qn*
      87             :  *       where q0 = 1, 4, -4, 8, qi* = 1 mod 4 and |qi| is a prime.
      88             :  *       The factorization is a vecsmall listing the indices of the qi* as
      89             :  *       they appear in the primelist (q0 = 1 is omitted) */
      90             : INLINE long
      91    57997688 : Dinfo_get_D(GEN Dinfo) { return gel(Dinfo, 1)[1]; }
      92             : INLINE long
      93    28061362 : Dinfo_get_h(GEN Dinfo) { return gel(Dinfo, 1)[2]; }
      94             : INLINE long
      95    21023184 : Dinfo_get_bi(GEN Dinfo) { return gel(Dinfo, 1)[3]; }
      96             : INLINE ulong
      97    57609895 : Dinfo_get_pd(GEN Dinfo) { return umael(Dinfo, 1, 4); }
      98             : INLINE GEN
      99    28594398 : Dinfo_get_Dfac(GEN Dinfo) { return gel(Dinfo, 2); }
     100             : 
     101             : /* primelist and indexlist
     102             :  *
     103             :  * primelist begins with 8, -4, -8; subsequent elements are the p^* for
     104             :  * successive odd primes <= maxsqrt, by increasing absolute value
     105             :  *
     106             :  *        i |  1 |  2 |  3 |  4 |  5 |  6 |  7 |  8 |  9 |  10 |  11 |
     107             :  * ---------+----+----+----+----+----+----+----+----+----+-----+-----+
     108             :  *        i |  1 |  2 |  3 |  4 |  5 |  6 |  7 |  8 |  9 |  10 |  11 |
     109             :  * Plist[i] |  8 | -4 | -8 | -3 |  5 | -7 |-11 | 13 | 17 | -19 | -23 |
     110             :  *        p |  3 |  5 |  7 |  9 | 11 | 13 | 15 | 17 | 19 |  21 |  23 |
     111             :  * ---------+----+----+----+----+----+----+----+----+----+-----+-----+*/
     112             : 
     113             : /* primelist containing primes whose absolute value is at most maxsqrt */
     114             : static GEN
     115          63 : ecpp_primelist_init(long maxsqrt)
     116             : {
     117          63 :   GEN P = cgetg(maxsqrt, t_VECSMALL);
     118             :   forprime_t T;
     119          63 :   long p, j = 4;
     120          63 :   u_forprime_init(&T, 3, maxsqrt);
     121          63 :   P[1] =  8; P[2] = -4; P[3] = -8;
     122          63 :   while((p = u_forprime_next(&T))) P[j++] = ((p & 3UL) == 1)? p: -p;
     123          63 :   setlg(P,j); return P;
     124             : }
     125             : 
     126             : static GEN
     127         924 : Dfac_to_p(GEN x, GEN P) { pari_APPLY_long(uel(P,x[i])); }
     128             : static GEN
     129         924 : Dfac_to_roots(GEN x, GEN P) { pari_APPLY_type(t_VEC, gel(P,x[i])); }
     130             : 
     131             : #if 0
     132             : INLINE ulong
     133             : ecpp_param_get_maxsqrt(GEN param) { return umael3(param, 1, 1, 1); }
     134             : INLINE ulong
     135             : ecpp_param_get_maxdisc(GEN param) { return umael3(param, 1, 1, 2); }
     136             : #endif
     137             : INLINE ulong
     138         980 : ecpp_param_get_maxpcdg(GEN param) { return umael3(param, 1, 1, 3); }
     139             : INLINE GEN
     140      536606 : ecpp_param_get_primelist(GEN param) { return gmael(param, 1, 2); }
     141             : INLINE GEN
     142         980 : ecpp_param_get_disclist(GEN param) { return gmael(param, 1, 3); }
     143             : INLINE GEN
     144        3409 : ecpp_param_get_primorial_vec(GEN param) { return gel(param, 2); }
     145             : INLINE GEN
     146        4389 : ecpp_param_get_tune(GEN param) { return gel(param, 3); }
     147             : 
     148             : /*  Input: x, 20 <= x <= 35
     149             :  * Output: a vector whose ith entry is the product of all primes below 2^x */
     150             : static GEN
     151          63 : primorial_vec(ulong x)
     152             : {
     153          63 :   pari_sp av = avma;
     154          63 :   long i, y = x-19;
     155          63 :   GEN v = primes_upto_zv(1UL << x), w = cgetg(y+1, t_VEC);
     156             :   /* ind[i]th prime number is the largest prime <= 2^(20+i) */
     157          63 :   long ind[] = { 0, 82025L, 155611L, 295947L, 564163L, 1077871L, 2063689L,
     158             :                  3957809L, 7603553L, 14630843L, 28192750L, 54400028L,
     159             :                  105097565L, 203280221L, 393615806L, 762939111L, 1480206279L};
     160          63 :   gel(w,1) = zv_prod_Z(vecslice(v,1,ind[1]));
     161         126 :   for (i = 2; i <= y; i++)
     162             :   {
     163          63 :     pari_sp av2 = avma;
     164          63 :     GEN q = mulii(gel(w,i-1), zv_prod_Z(vecslice(v,ind[i-1]+1,ind[i])));
     165          63 :     gel(w,i) = gerepileuptoint(av2, q);
     166             :   }
     167          63 :   return gerepileupto(av, w);
     168             : }
     169             : 
     170             : /* NDmqg contains
     171             :    N, as in the theorem in ??ecpp
     172             :    Dinfo, as discussed in the comment about Dinfo
     173             :    m, as in the theorem in ??ecpp
     174             :    q, as in the theorem in ??ecpp
     175             :    g, a quadratic non-residue modulo N
     176             :    sqrt, a list of squareroots
     177             : */
     178             : INLINE GEN
     179         924 : NDmqg_get_N(GEN x) { return gel(x,1); }
     180             : INLINE GEN
     181        8106 : NDmqg_get_Dinfo(GEN x) { return gel(x,2); }
     182             : INLINE GEN
     183         924 : NDmqg_get_m(GEN x) { return gel(x,3); }
     184             : INLINE GEN
     185         924 : NDmqg_get_q(GEN x) { return gel(x,4); }
     186             : INLINE GEN
     187         924 : NDmqg_get_g(GEN x) { return gel(x,5); }
     188             : INLINE GEN
     189         924 : NDmqg_get_sqrt(GEN x) { return gel(x,6); }
     190             : 
     191             : /* COMPARATOR FUNCTIONS */
     192             : 
     193             : static int
     194    28726999 : sort_disclist(void *data, GEN x, GEN y)
     195             : {
     196             :   long d1, h1, g1, o1, pd1, hf1, wD1, d2, h2, g2, o2, pd2, hf2, wD2;
     197             :   (void)data;
     198    28726999 :   d1 = Dinfo_get_D(x); wD1 = D_get_wD(d1);
     199    28726999 :   d2 = Dinfo_get_D(y); wD2 = D_get_wD(d2);
     200             :   /* higher number of units means more elliptic curves to try */
     201    28726999 :   if (wD1 != wD2) return wD2 > wD1 ? 1 : -1;
     202             :   /* lower polclass degree is better: faster computation of roots modulo N */
     203    28725529 :   pd1 = Dinfo_get_pd(x); /* degree of polclass */
     204    28725529 :   pd2 = Dinfo_get_pd(y);
     205    28725529 :   if (pd1 != pd2) return pd1 > pd2 ? 1 : -1;
     206    14028483 :   g1 = lg(Dinfo_get_Dfac(x))-1; /* genus number */
     207    14028483 :   h1 = Dinfo_get_h(x); /* class number */
     208    14028483 :   o1 = h1 >> (g1-1); /* odd class number */
     209    14028483 :   g2 = lg(Dinfo_get_Dfac(y))-1;
     210    14028483 :   h2 = Dinfo_get_h(y);
     211    14028483 :   o2 = h2 >> (g2-1);
     212    14028483 :   if (o1 != o2) return g1 > g2 ? 1 : -1;
     213             :   /* lower class number is better: higher probability of succesful cornacchia */
     214    11024825 :   if (h1 != h2) return h1 > h2 ? 1 : -1;
     215             :   /* higher height factor is better: polclass would have lower coefficients */
     216    10511130 :   hf1 = modinv_height_factor(Dinfo_get_bi(x)); /* height factor for best inv. */
     217    10511130 :   hf2 = modinv_height_factor(Dinfo_get_bi(y));
     218    10511130 :   if (hf1 != hf2) return hf2 > hf1 ? 1 : -1;
     219             :   /* "higher" discriminant is better since its absolute value is lower */
     220     4854640 :   if (d1 != d2) return d2 > d1 ? 1 : -1;
     221           0 :   return 0;
     222             : }
     223             : 
     224             : INLINE long
     225        7182 : NDmqg_get_D(GEN x) { return Dinfo_get_D(NDmqg_get_Dinfo(x)); }
     226             : static int
     227        3591 : sort_NDmq_by_D(void *data, GEN x, GEN y)
     228             : {
     229        3591 :   long d1 = NDmqg_get_D(x);
     230        3591 :   long d2 = NDmqg_get_D(y);
     231        3591 :   (void)data; return d2 > d1 ? 1 : -1;
     232             : }
     233             : 
     234             : static int
     235       55678 : sort_Dmq_by_q(void *data, GEN x, GEN y)
     236       55678 : { (void)data; return cmpii(gel(x,3), gel(y,3)); }
     237             : 
     238             : INLINE long
     239           0 : Dmq_get_D(GEN Dmq) { return Dinfo_get_D(gel(Dmq,1)); }
     240             : INLINE long
     241        4396 : Dmq_get_h(GEN Dmq) { return Dinfo_get_h(gel(Dmq,1)); }
     242             : static int
     243        2198 : sort_Dmq_by_cnum(void *data, GEN x, GEN y)
     244             : {
     245        2198 :   ulong h1 = Dmq_get_h(x);
     246        2198 :   ulong h2 = Dmq_get_h(y);
     247        2198 :   if (h1 != h2) return h1 > h2 ? 1 : -1;
     248         805 :   return sort_Dmq_by_q(data, x, y);
     249             : }
     250             : 
     251             : /* return H s.t if -maxD <= D < 0 is fundamental then H[(-D)>>1] is the
     252             :  * ordinary class number of Q(sqrt(D)); junk at other entries. */
     253             : static GEN
     254          63 : allh(ulong maxD)
     255             : {
     256          63 :   ulong a, A = usqrt(maxD/3), maxD2 = maxD >> 1;
     257          63 :   GEN H = zero_zv(maxD2);
     258       15393 :   for (a = 1; a <= A; a++)
     259             :   {
     260       15330 :     ulong a2 = a << 1, aa = a*a, aa4 = aa << 2, b, c;
     261             :     { /* b = 0 */
     262       15330 :       ulong D = aa << 1;
     263       15330 :       for (c = a; D <= maxD2; c++) { H[D]++; D += a2; }
     264             :     }
     265     2821805 :     for (b = 1; b < a; b++)
     266             :     {
     267     2808505 :       ulong B = b*b, D = (aa4 - B) >> 1;
     268     2808505 :       if (D > maxD2) break;
     269     2806475 :       H[D]++; D += a2; /* c = a */
     270     2806475 :       for (c = a+1; D <= maxD2; c++) { H[D] += 2; D += a2; }
     271             :     }
     272             :     { /* b = a */
     273       15330 :       ulong D = (aa4 - aa) >> 1;
     274       15330 :       for (c = a; D <= maxD2; c++) { H[D]++; D += a2; }
     275             :     }
     276             :   }
     277          63 :   return H;
     278             : }
     279             : 
     280             : static GEN
     281     1895915 : mkDinfo(GEN c, long D, long h)
     282             : {
     283             :   long bi, pd, p1, p2;
     284     1895915 :   bi = disc_best_modinv(D);
     285     1895915 :   pd = (modinv_degree(&p1,&p2,bi) && (-D)%p1==0 && (-D)%p2==0)? h/2: h;
     286     1895915 :   gel(c,1) = mkvecsmall4(D, h, bi, pd); return c;
     287             : }
     288             : 
     289             : static GEN
     290          63 : ecpp_disclist_init(ulong maxdisc, GEN primelist)
     291             : {
     292          63 :   pari_sp av = avma;
     293             :   long t, ip, u, lp, lmerge;
     294          63 :   GEN merge, ev, od, Harr = allh(maxdisc); /* table of class numbers*/
     295          63 :   long lenv = maxdisc/4; /* max length of od/ev */
     296             :   long N; /* maximum number of positive prime factors */
     297             : 
     298             :   /* tuning paramaters blatantly copied from vecfactoru */
     299          63 :   if (maxdisc < 510510UL) N = 7;
     300          14 :   else if (maxdisc < 9699690UL) N = 8;
     301             :   #ifdef LONG_IS_64BIT
     302           0 :     else if (maxdisc < 223092870UL) N = 9;
     303           0 :     else if (maxdisc < 6469693230UL) N = 10;
     304           0 :     else if (maxdisc < 200560490130UL) N = 11;
     305           0 :     else if (maxdisc < 7420738134810UL) N = 12;
     306           0 :     else if (maxdisc < 304250263527210UL) N = 13;
     307           0 :     else N = 16; /* don't bother */
     308             :   #else
     309           0 :     else N = 9;
     310             :   #endif
     311             : 
     312             :   /* od[t] attached to discriminant 1-4*t, ev[t] attached to -4*t */
     313          63 :   od = cgetg(lenv + 1, t_VEC); /* contains 'long' at first: save memory */
     314          63 :   ev = cgetg(lenv + 1, t_VEC);
     315             :   /* first pass: squarefree sieve and restrict to maxsqrt-smooth disc */
     316     5626313 :   for (t = 1; t <= lenv; t++)
     317             :   {
     318     5626250 :     od[t] = 1;
     319     5626250 :     switch(t&7)
     320             :     {
     321     1406559 :       case 0: case 4: ev[t] = 0; break;
     322      703297 :       case 2: ev[t] =-8; break;
     323      703262 :       case 6: ev[t] = 8; break;
     324     2813132 :       default:ev[t] =-4; break;
     325             :     }
     326             :   }
     327          63 :   lp = lg(primelist);
     328        4865 :   for (ip = 4; ip < lp; ip++) /* skip 8,-4,-8 */
     329             :   { /* sieve by squares of primes */
     330        4802 :     long s, q = primelist[ip], p = labs(q);
     331        4802 :     s = (q == p)? 3: 1;
     332     9551948 :     for (t = (s*p+1)>>2; t <= lenv; t += p, s += 4)
     333             :     {
     334     9547146 :       long c = od[t];
     335     9547146 :       if (c) { if (s%p == 0) od[t] = 0; else  od[t] = c*q; }
     336             :     }
     337        4802 :     s = 1;
     338     9549512 :     for (t = p; t <= lenv; t += p, s++)
     339             :     {
     340     9544710 :       long c = ev[t];
     341     9544710 :       if (c) { if (s%p == 0) ev[t] = 0; else  ev[t] = c*q; }
     342             :     }
     343             :   }
     344     5626313 :   for (u = 0, t = 1; t <= lenv; t++)
     345             :   { /* restrict to maxsqrt-smooth disc */
     346     5626250 :     if (od[t] != -4*t+1) od[t] = 0; else u++;
     347     5626250 :     if (ev[t] != -4*t)   ev[t] = 0; else u++;
     348             :   }
     349             : 
     350             :   /* second pass: sieve by primes and record factorization */
     351     5626313 :   for (t = 1; t <= lenv; t++)
     352             :   {
     353     5626250 :     if (od[t]) gel(od,t) = mkvec2(NULL, vecsmalltrunc_init(N));
     354     5626250 :     if (ev[t])
     355             :     {
     356      820435 :       GEN F = vecsmalltrunc_init(N);
     357             :       long id;
     358      820435 :       switch(t&7)
     359             :       {
     360      220339 :         case 2: id = 3; break;
     361      220311 :         case 6: id = 1; break;
     362      379785 :         default:id = 2; break;
     363             :       }
     364      820435 :       vecsmalltrunc_append(F, id);
     365      820435 :       gel(ev,t) = mkvec2(NULL, F);
     366             :     }
     367             :   }
     368          63 :   lp = lg(primelist);
     369        4865 :   for (ip = 4; ip < lp; ip++) /* skip 8,-4,-8 */
     370             :   {
     371        4802 :     long s, q = primelist[ip], p = labs(q);
     372        4802 :     s = (q == p)? 3: 1;
     373     9551948 :     for (t = (s*p+1)>>2; t <= lenv; t += p, s += 4)
     374             :     {
     375     9547146 :       GEN c = gel(od,t);
     376     9547146 :       if (c) vecsmalltrunc_append(gel(c,2), ip);
     377             :     }
     378        4802 :     s = 1;
     379     9549512 :     for (t = p; t <= lenv; t += p, s++)
     380             :     {
     381     9544710 :       GEN c = gel(ev,t);
     382     9544710 :       if (c) vecsmalltrunc_append(gel(c,2), ip);
     383             :     }
     384             :   }
     385             :   /* merging the two arrays */
     386          63 :   merge = cgetg(u+1, t_VEC); lmerge = 0;
     387     5626313 :   for (t = 1; t <= lenv; t++)
     388             :   {
     389             :     GEN c;
     390     5626250 :     c = gel(od,t); if (c) gel(merge, ++lmerge) = mkDinfo(c,1-4*t, Harr[2*t-1]);
     391     5626250 :     c = gel(ev,t); if (c) gel(merge, ++lmerge) = mkDinfo(c, -4*t, Harr[2*t]);
     392             :   }
     393          63 :   setlg(merge, lmerge);
     394          63 :   gen_sort_inplace(merge, NULL, &sort_disclist, NULL);
     395          63 :   return gerepilecopy(av, merge);
     396             : }
     397             : 
     398             : /*  Input: a vector tune whose components are [maxsqrt,maxpcdg,tdivexp,expiN]
     399             :  * Output: vector param of precomputations
     400             :  *   let x =  be a component of tune then
     401             :  *   param[1][1] = [maxsqrt, maxdisc, maxpcdg]
     402             :  *   param[1][2] = primelist = Vecsmall of primes
     403             :  *   param[1][3] = disclist  = vector of Dinfos, sorted by quality
     404             :  *   param[2]    = primorial_vec
     405             :  *   param[3]    = tune */
     406             : static GEN
     407          63 : ecpp_param_set(GEN tune, GEN x)
     408             : {
     409          63 :   pari_sp av = avma;
     410          63 :   ulong maxsqrt = uel(x,1), maxpcdg = uel(x,2), tdivexp = uel(x,3);
     411          63 :   ulong maxdisc = maxsqrt * maxsqrt;
     412          63 :   GEN T = mkvecsmall3(maxsqrt, maxdisc, maxpcdg);
     413          63 :   GEN Plist = ecpp_primelist_init(maxsqrt);
     414          63 :   GEN Dlist = ecpp_disclist_init(maxdisc, Plist);
     415          63 :   GEN primorial = primorial_vec(tdivexp);
     416          63 :   return gerepilecopy(av, mkvec3(mkvec3(T,Plist,Dlist), primorial, tune));
     417             : }
     418             : 
     419             : /* cert contains [N, t, s, a4, [x, y]] as documented in ??ecpp; the following
     420             :  * information can be obtained:
     421             :  * D = squarefreepart(t^2-4N)
     422             :  * m = (N+1-t), q = m/s, a6 = y^3 - x^2 - a4*x, J */
     423             : INLINE GEN
     424         854 : cert_get_N(GEN x) { return gel(x,1); }
     425             : INLINE GEN
     426         784 : cert_get_t(GEN x) { return gel(x,2); }
     427             : INLINE GEN
     428          28 : cert_get_D(GEN x)
     429             : {
     430          28 :   GEN N = cert_get_N(x), t = cert_get_t(x);
     431          28 :   return coredisc(subii(sqri(t), shifti(N, 2)));
     432             : }
     433             : INLINE GEN
     434         406 : cert_get_m(GEN x)
     435             : {
     436         406 :   GEN N = cert_get_N(x), t = cert_get_t(x);
     437         406 :   return subii(addiu(N, 1), t);
     438             : }
     439             : INLINE GEN
     440         406 : cert_get_s(GEN x) { return gel(x,3); }
     441             : INLINE GEN
     442          56 : cert_get_q(GEN x) { return diviiexact(cert_get_m(x), cert_get_s(x)); }
     443             : INLINE GEN
     444         420 : cert_get_a4(GEN x) { return gel(x,4); }
     445             : INLINE GEN
     446         406 : cert_get_P(GEN x) { return gel(x,5); }
     447             : INLINE GEN
     448          14 : cert_get_x(GEN x) { return gel(cert_get_P(x), 1); }
     449             : INLINE GEN
     450          42 : cert_get_a6(GEN z)
     451             : {
     452          42 :   GEN N = cert_get_N(z), a = cert_get_a4(z), P = cert_get_P(z);
     453          42 :   GEN x = gel(P,1), xx = Fp_sqr(x, N);
     454          42 :   GEN y = gel(P,2), yy = Fp_sqr(y, N);
     455          42 :   return Fp_sub(yy, Fp_mul(x, Fp_add(xx,a,N), N), N);
     456             : }
     457             : INLINE GEN
     458          28 : cert_get_E(GEN x) { return mkvec2(cert_get_a4(x), cert_get_a6(x)); }
     459             : INLINE GEN
     460           0 : cert_get_J(GEN x)
     461             : {
     462           0 :   GEN a = cert_get_a4(x), b = cert_get_a6(x), N = cert_get_N(x);
     463           0 :   return Fp_ellj(a, b, N);
     464             : }
     465             : 
     466             : /* Given J, N, set A = 3*J*(1728-J) mod N, B = 2*J*(1728-J)^2 mod N */
     467             : static void
     468         924 : Fp_ellfromj(GEN j, GEN N, GEN *A, GEN *B)
     469             : {
     470             :   GEN k, jk;
     471         924 :   j = Fp_red(j, N);
     472         924 :   if (isintzero(j)) { *A = gen_0; *B = gen_1; return; }
     473         665 :   if (absequalui(umodui(1728,N), j)) { *A = gen_1; *B = gen_0; return; }
     474         539 :   k = Fp_sub(utoi(1728), j, N);
     475         539 :   jk = Fp_mul(j, k, N);
     476         539 :   *A = Fp_mulu(jk, 3, N);
     477         539 :   *B = Fp_mulu(Fp_mul(j, Fp_sqr(k,N), N), 2, N);
     478             : }
     479             : 
     480             : /* "Twist factor". Does not cover J = 0, 1728 */
     481             : static GEN
     482           0 : cert_get_lambda(GEN z)
     483             : {
     484             :   GEN N, J, a, b, A, B;
     485           0 :   J = cert_get_J(z);
     486           0 :   N = cert_get_N(z);
     487           0 :   a = cert_get_a4(z);
     488           0 :   b = cert_get_a6(z);
     489           0 :   Fp_ellfromj(J, N, &A, &B);
     490           0 :   return Fp_div(Fp_mul(a,B,N), Fp_mul(b,A,N), N);
     491             : }
     492             : 
     493             : /* Solves for T such that if
     494             :    [A, B] = [3*J*(1728-J), 2*J*(1728-J)^2]
     495             :    and you let
     496             :    L = T^3 + A*T + B, A = A*L^2, B = B*L^3
     497             :    then
     498             :    x == TL and y == L^2
     499             : */
     500             : static GEN
     501           0 : cert_get_T(GEN z)
     502             : {
     503           0 :   GEN N = cert_get_N(z), x = cert_get_x(z);
     504           0 :   GEN l = cert_get_lambda(z); /* l = 1/L */
     505           0 :   return Fp_mul(x, l, N);
     506             : }
     507             : 
     508             : /* database for all invariants, stolen from Hamish's code */
     509             : static GEN
     510          42 : polmodular_db_init_allinv(void)
     511             : {
     512          42 :   const long DEFAULT_MODPOL_DB_LEN = 32;
     513          42 :   GEN a, b = cgetg(39+1, t_VEC);
     514             :   long i;
     515          42 :   for (i = 1; i < 40; i++) gel(b,i) = zerovec_block(DEFAULT_MODPOL_DB_LEN);
     516          42 :   a = zerovec_block(DEFAULT_MODPOL_DB_LEN);
     517          42 :   return mkvec2(a, b);
     518             : }
     519             : 
     520             : /* Given D and a database of modular polynomials, return polclass(D, inv) */
     521             : static GEN
     522         525 : D_polclass(long D, long inv, GEN *db)
     523             : {
     524         525 :   GEN HD, t = mkvec2(gel(*db, 1), inv == 0? gen_0: gmael(*db, 2, inv));
     525         525 :   HD = polclass0(D, inv, 0, &t);
     526         525 :   gel(*db, 1) = gel(t,1);
     527         525 :   if (inv != 0) gmael(*db, 2, inv) = gel(t,2);
     528         525 :   return HD;
     529             : }
     530             : 
     531             : static GEN
     532        2151 : NqE_check(GEN N, GEN q, GEN a, GEN b, GEN s)
     533             : {
     534        2151 :   GEN mP, sP, P = random_FpE(a, b, N);
     535        2151 :   GEN PJ = mkvec3(gel(P,1), gel(P,2), gen_1);
     536        2151 :   sP = FpJ_mul(PJ, s, a, N); if (FpJ_is_inf(sP)) return NULL;
     537        2151 :   mP = FpJ_mul(sP, q, a, N); return FpJ_is_inf(mP)? mkvec2(a, P): NULL;
     538             : }
     539             : 
     540             : /* Find an elliptic curve E modulo N and a point P on E which corresponds to
     541             :  * the ith element of the certificate; uses N, D, m, q, J from previous steps.
     542             :  * All computations are modulo N unless stated otherwise */
     543             : 
     544             : /* g is a quadratic and cubic non-residue modulo N */
     545             : static GEN
     546         182 : j0_find_g(GEN N)
     547             : {
     548         182 :   GEN n = diviuexact(subiu(N, 1), 3);
     549             :   for(;;)
     550         421 :   {
     551         603 :     GEN g = randomi(N);
     552         785 :     if (kronecker(g, N) == -1 && !equali1(Fp_pow(g, n, N))) return g;
     553             :   }
     554             : }
     555             : 
     556             : static GEN
     557         924 : find_EP(GEN N, long D, GEN q, GEN g, GEN J, GEN s)
     558             : {
     559         924 :   GEN A0, B0; Fp_ellfromj(J, N, &A0, &B0);
     560             :   for(;;)
     561           0 :   { /* expect one iteration: not worth saving the A's and B's */
     562         924 :     GEN gg, v, A = A0, B = B0;
     563             :     long i;
     564         924 :     if ((v = NqE_check(N, q, A, B, s))) return v;
     565         595 :     switch (D_get_wD(D))
     566             :     {
     567             :       case 2:
     568         273 :         gg = Fp_sqr(g, N);
     569         273 :         A = Fp_mul(A, gg, N);
     570         273 :         B = Fp_mul(Fp_mul(B, gg, N), g, N);
     571         273 :         if ((v = NqE_check(N, q, A, B, s))) return v;
     572           0 :         break;
     573             :       case 4:
     574         224 :         for (i = 1; i < 4; i++)
     575             :         {
     576         224 :           A = Fp_mul(A, g, N);
     577         224 :           if ((v = NqE_check(N, q, A, B, s))) return v;
     578             :         }
     579           0 :         break;
     580             :       default: /* 6 */
     581         217 :         B = Fp_mul(B, g, N);
     582         217 :         if ((v = NqE_check(N, q, A, B, s))) return v;
     583         182 :         g = j0_find_g(N);
     584         646 :         for (i = 1; i < 6; i++)
     585             :         {
     586         646 :           B = Fp_mul(B, g, N);
     587         646 :           if (i == 3) continue;
     588         513 :           if ((v = NqE_check(N, q, A, B, s))) return v;
     589             :         }
     590           0 :         break;
     591             :     }
     592             :   }
     593             : }
     594             : 
     595             : /* Convert the disc. factorisation of a genus field to the
     596             :  * disc. factorisation of its real part. */
     597             : static GEN
     598         399 : realgenusfield(GEN Dfac, GEN sq, GEN p)
     599             : {
     600         399 :   long i, j, l = lg(Dfac), dn, n = 0;
     601         399 :   GEN sn, s = gen_0, R = cgetg(l-1, t_VECSMALL);
     602         448 :   for (i = 1; i < l; i++)
     603         448 :     if (Dfac[i] < 0) { n = i; break; }
     604         399 :   if (n == 0) pari_err_BUG("realgenusfield");
     605         399 :   dn = Dfac[n]; sn = gel(sq, n);
     606        1211 :   for (j = i = 1; i < l; i++)
     607         812 :     if (i != n)
     608             :     {
     609         413 :       long d = Dfac[i];
     610         413 :       GEN si = gel(sq,i);
     611         413 :       R[j++] = d > 0 ? d : d * dn;
     612         413 :       s = Fp_add(s, d > 0? si: Fp_mul(si, sn, p), p);
     613             :     }
     614         399 :   return mkvec2(R, s);
     615             : }
     616             : 
     617             : static GEN
     618         924 : FpX_classtower_oneroot(GEN P, GEN Dfac, GEN sq, GEN p)
     619             : {
     620             :   pari_timer ti;
     621         924 :   pari_sp av = avma;
     622             :   GEN C;
     623         924 :   dbg_mode() timer_start(&ti);
     624         924 :   if (degpol(P) > 1)
     625             :   {
     626         399 :     GEN N = NULL, V = realgenusfield(Dfac, sq, p), v = gel(V,1), R = gel(V,2);
     627         399 :     long i, l = lg(v);
     628         812 :     for (i = 1; i < l; i++)
     629             :     {
     630         413 :       GEN Q = deg2pol_shallow(gen_1, gen_0, stoi(-v[i]), 1);
     631         413 :       if (!N) N = Q;
     632             :       else
     633             :       {
     634          77 :         GEN cm = polcompositum0(N,Q,3);
     635          77 :         N = gel(cm,1); P = gsubst(P, 1, gel(cm,2));
     636             :       }
     637         413 :       P = liftpol_shallow(gmael(nffactor(N,P),1,1));
     638             :     }
     639         399 :     if (N)
     640             :     {
     641         336 :       P = FpXY_evalx(Q_primpart(P), R, p);
     642         336 :       dbg_mode() err_printf(ANSI_BRIGHT_GREEN " %6ld[%ld]" ANSI_RESET,
     643             :                             timer_delay(&ti), l-1);
     644             :     } else
     645          63 :       dbg_mode() err_printf("          ");
     646             :   }
     647         924 :   C = FpX_oneroot_split(P, p);
     648         924 :   dbg_mode() err_printf(" %6ld", timer_delay(&ti));
     649         924 :   return gerepileupto(av, C);
     650             : }
     651             : 
     652             : /* This uses [N, D, m, q] from step 1 to find the appropriate j-invariants
     653             :  * to be used in step 3. Step 2 is divided into substeps 2a, 2b, 2c */
     654             : static GEN
     655          42 : ecpp_step2(GEN step1, GEN *X0, GEN primelist)
     656             : {
     657          42 :   long j, Dprev = 0;
     658             :   pari_timer ti;
     659          42 :   GEN perm = gen_indexsort(step1, NULL, &sort_NDmq_by_D);
     660          42 :   GEN step2 = cgetg(lg(step1), t_VEC);
     661          42 :   GEN HD = NULL, db = polmodular_db_init_allinv();
     662             : 
     663         966 :   for (j = 1; j < lg(step2); j++)
     664             :   {
     665         924 :     long i = uel(perm, j);
     666         924 :     GEN J, t, s, EP, rt, S = gel(step1, i);
     667         924 :     GEN N = NDmqg_get_N(S), Dinfo = NDmqg_get_Dinfo(S);
     668         924 :     GEN m = NDmqg_get_m(S), q = NDmqg_get_q(S);
     669         924 :     GEN g = NDmqg_get_g(S), sq = NDmqg_get_sqrt(S);
     670         924 :     long D = Dinfo_get_D(Dinfo), inv = Dinfo_get_bi(Dinfo);
     671         924 :     GEN Dfacp = Dfac_to_p(Dinfo_get_Dfac(Dinfo), primelist);
     672             : 
     673             :     /* C1: Find the appropriate class polynomial modulo N */
     674         924 :     dbg_mode() timer_start(&ti);
     675         924 :     if (D != Dprev) HD = D_polclass(D, inv, &db);
     676         924 :     dbg_mode() {
     677           0 :       long tt = timer_record(X0, "C1", &ti);
     678           0 :       err_printf(ANSI_BRIGHT_GREEN "\n[ %3d | %4ld bits]" ANSI_RESET, i, expi(N));
     679           0 :       err_printf(ANSI_GREEN " D = %8ld poldeg = %4ld" ANSI_RESET, D, degpol(HD));
     680           0 :       if (D == Dprev) err_printf(" %6ld", tt);
     681           0 :       else err_printf(ANSI_BRIGHT_WHITE " %6ld" ANSI_RESET, tt);
     682             :     }
     683             :     /* C2: Find a root modulo N of polclass(D,inv) */
     684         924 :     dbg_mode() timer_start(&ti);
     685         924 :     rt = FpX_classtower_oneroot(HD, Dfacp, sq, N);
     686         924 :     dbg_mode() err_printf(" %6ld", timer_record(X0, "C2", &ti));
     687             : 
     688             :     /* C3: Convert root from previous step into the appropriate j-invariant */
     689         924 :     dbg_mode() timer_start(&ti);
     690         924 :     J = Fp_modinv_to_j(rt, inv, N); /* root of polclass(D) */
     691         924 :     dbg_mode() err_printf(" %6ld", timer_record(X0, "C3", &ti));
     692             : 
     693             :     /* D1: Find an elliptic curve E with a point P satisfying the theorem */
     694         924 :     dbg_mode() timer_start(&ti);
     695         924 :     s = diviiexact(m, q);
     696         924 :     EP = find_EP(N, D, q, g, J, s);
     697         924 :     dbg_mode() err_printf(" %6ld", timer_record(X0, "D1", &ti));
     698             : 
     699             :     /* D2: Compute for t and s */
     700         924 :     t = subii(addiu(N, 1), m); /* t = N+1-m */
     701         924 :     gel(step2, i) = mkvec5(N, t, s, gel(EP,1), gel(EP,2));
     702         924 :     Dprev = D;
     703             :   }
     704          42 :   gunclone_deep(db); return step2;
     705             : }
     706             : /* end of functions for step 2 */
     707             : 
     708             : /* start of functions for step 1 */
     709             : 
     710             : /* This finds the square root of D modulo N [given by Dfac]: find the square
     711             :  * root modulo N of each prime p dividing D and multiply them out */
     712             : static GEN
     713       40040 : D_find_discsqrt(GEN N, GEN primelist, GEN Dfac, GEN sqrtlist, GEN g)
     714             : {
     715       40040 :   GEN s = NULL;
     716       40040 :   long i, l = lg(Dfac);
     717      122178 :   for (i = 1; i < l; i++)
     718             :   {
     719       82138 :     long j = Dfac[i];
     720       82138 :     GEN sj = gel(sqrtlist,j);
     721       82138 :     if (!signe(sj))
     722             :     {
     723       14427 :       GEN p = stoi(primelist[j]);
     724       14427 :       dbg_mode() err_printf(ANSI_MAGENTA "S" ANSI_RESET);
     725             :       /* A4: Get the square root of a prime factor of D. */
     726       14427 :       sj = gel(sqrtlist, j) = Fp_sqrt_i(p, g, N);
     727       14427 :       if (!sj) pari_err_BUG("D_find_discsqrt"); ; /* possible if N composite */
     728             :     }
     729       82138 :     s = s? Fp_mul(s, sj, N): sj;
     730             :   }
     731       40040 :   return s;/* != NULL */
     732             : }
     733             : 
     734             : /* Given a solution U, V to U^2 + |D|V^2 = 4N, this find all the possible
     735             :      cardinalities of the elliptic curves over the finite field F_N
     736             :      whose endomorphism ring is isomorphic to the maximal order
     737             :      in the imaginary quadratic number field K = Q(sqrt(D)). ???
     738             : */
     739             : static GEN
     740       10416 : NUV_find_m(GEN Dinfo, GEN N, GEN U, GEN V, long wD)
     741             : {
     742       10416 :   GEN m, Nplus1 = addiu(N, 1), u = U, mvec = cgetg(wD+1, t_VEC);
     743             :   long i;
     744       34482 :   for (i = 0; i < wD; i++)
     745             :   {
     746       24066 :     if (odd(i)) m = subii(Nplus1, u);
     747             :     else
     748             :     {
     749       12033 :       if (wD == 4 && i==2) u = shifti(V, 1);
     750       11606 :       else if (wD == 6 && i==2) u = shifti(submuliu(U, V, 3), -1);
     751       11011 :       else if (wD == 6 && i==4) u = shifti(addmuliu(U, V, 3), -1);
     752       12033 :       m = addii(Nplus1, u);
     753             :     }
     754       24066 :     gel(mvec, i+1) = mkvec2(Dinfo, m);
     755             :   }
     756       10416 :   return mvec;
     757             : }
     758             : 
     759             : /* Populates Dmbatch with Dmvec's -- whose components are of the form [D,m],
     760             :    where m is a cardinalities of an elliptic curve with endomorphism ring
     761             :    isomorphic to the maximal order of imaginary quadratic K = Q(sqrt(D)).
     762             :    It returns 0 if:
     763             :      * Any of the p* dividing D is not a quadratic residue mod N
     764             :      * Cornacchia cannot find a solution U^2 + |D|V^2 = 4N.
     765             :    Otherwise, it returns the number of cardinalities added to Dmbatch.
     766             :    Finally, sqrtlist and g help compute the square root modulo N of D.
     767             : */
     768             : static long
     769      535584 : D_collectcards(GEN N, GEN param, GEN* X0, GEN Dinfo, GEN sqrtlist, GEN g, GEN Dmbatch, GEN kroP)
     770             : {
     771      535584 :   long i, l, corn_succ, wD, D = Dinfo_get_D(Dinfo);
     772      535584 :   GEN U, V, sqrtofDmodN, Dfac = Dinfo_get_Dfac(Dinfo);
     773      535584 :   GEN primelist = ecpp_param_get_primelist(param);
     774             :   pari_timer ti;
     775             : 
     776             :   /* A1: Check (p*|N) = 1 for all primes dividing D */
     777      535584 :   l = lg(Dfac);
     778      764946 :   for (i = 1; i < l; i++)
     779             :   {
     780      724906 :     long j = Dfac[i], s = kroP[j];
     781      724906 :     if (s > 1) kroP[j] = s = krosi(primelist[j], N); /* update cache */
     782      724906 :     if (s == 0) return -1; /* N is composite */
     783      724906 :     if (s < 0) return 0;
     784             :   }
     785             :   /* A3: Get square root of D mod N */
     786       40040 :   dbg_mode() timer_start(&ti);
     787       40040 :   sqrtofDmodN = D_find_discsqrt(N, primelist, Dfac, sqrtlist, g);
     788       40040 :   dbg_mode() timer_record(X0, "A3", &ti);
     789             :   /* A5: Use square root with Cornacchia to solve U^2 + |D|V^2 = 4N */
     790       40040 :   dbg_mode() timer_start(&ti);
     791       40040 :   corn_succ = cornacchia2_sqrt( utoi(labs(D)), N, sqrtofDmodN, &U, &V);
     792       40040 :   dbg_mode() timer_record(X0, "A5", &ti);
     793       40040 :   if (!corn_succ) {
     794       29624 :     dbg_mode() err_printf(ANSI_YELLOW "c" ANSI_RESET);
     795       29624 :     return 0;
     796             :   }
     797             :   /* A6: Collect the w(D) cardinalities of E/(F_N) with CM by D */
     798       10416 :   dbg_mode() err_printf(ANSI_BRIGHT_YELLOW "D" ANSI_RESET);
     799       10416 :   wD = D_get_wD(D);
     800       10416 :   vectrunc_append_batch(Dmbatch, NUV_find_m(Dinfo,N,U,V,wD));
     801       10416 :   return wD;
     802             : }
     803             : 
     804             : /* Compute (S(16N, 2) + S(4096N, 4) + 4)\4 + 1,  where S is the nth root
     805             :  * rounded down. This is at most one more than (N^1/4 + 1)^2. */
     806             : static GEN
     807        3409 : ecpp_qlo(GEN N)
     808             : {
     809        3409 :   GEN a = sqrtnint(shifti(N, 4), 2);
     810        3409 :   GEN b = sqrtnint(shifti(N, 12), 4);
     811        3409 :   return addiu(shifti(addii(a, b), -2), 2);
     812             : }
     813             : 
     814             : static long
     815       12194 : lessthan_qlo(void* E, GEN Dmq) { return (cmpii(gel(Dmq,3), (GEN)E) < 0); }
     816             : static long
     817       10934 : gained_bits(void* E, GEN Dmq) { return (expi(gel(Dmq,3)) <= (long)E); }
     818             : 
     819             : /*  Input: Dmqvec
     820             :  * Output: Dmqvec such that q satisfies (N^1/4 + 1)^2 < q < N/2 */
     821             : static GEN
     822        3409 : Dmqvec_slice(GEN N, GEN Dmqvec)
     823             : {
     824             :   long lo, hi;
     825             : 
     826        3409 :   gen_sort_inplace(Dmqvec, NULL, &sort_Dmq_by_q, NULL); /* sort wrt q */
     827        3409 :   lo = zv_binsearch0((void*)ecpp_qlo(N), &lessthan_qlo, Dmqvec); lo++;
     828        3409 :   if (lo >= lg(Dmqvec)) return NULL;
     829             : 
     830        3409 :   hi = zv_binsearch0((void*)(expi(N)-1), &gained_bits, Dmqvec);
     831        3409 :   if (hi == 0) return NULL;
     832             : 
     833        3409 :   return vecslice(Dmqvec, lo, hi);
     834             : }
     835             : 
     836             : /* Given a Dmvec of [D,m]'s, simultaneously remove all prime factors of each
     837             :  * m less then BOUND_PRIMORIAL. This implements Franke 2004: Proving the
     838             :  * Primality of Very Large Numbers with fastECPP */
     839             : static GEN
     840        3409 : Dmvec_batchfactor(GEN Dmvec, GEN primorial)
     841             : {
     842        3409 :   long i, l = lg(Dmvec);
     843        3409 :   GEN leaf, v = cgetg(l, t_VEC);
     844       27475 :   for (i = 1; i < l; i++)
     845             :   { /* cheaply remove powers of 2 */
     846       24066 :     GEN m = gmael(Dmvec, i, 2);
     847       24066 :     long v2 = vali(m);
     848       24066 :     gel(v,i) = v2? shifti(m,-v2): m;
     849             :   }
     850        3409 :   leaf = Z_ZV_mod_tree(primorial, v, ZV_producttree(v));
     851             :   /* Go through each leaf and remove small factors. */
     852       27475 :   for (i = 1; i < l; i++)
     853             :   {
     854       24066 :     GEN q = gel(v,i), Dm = gel(Dmvec,i), L = gel(leaf,i);
     855       24066 :     while (!equali1(L)) { L = gcdii(q, L); q = diviiexact(q, L); }
     856       24066 :     gel(v,i) = mkvec3(gel(Dm,1), gel(Dm,2), q);
     857             :   }
     858        3409 :   return v;
     859             : }
     860             : 
     861             : /* return good parameters [maxsqrt, maxpcdg, tdivexp] for ecpp(N) */
     862             : static GEN
     863        4389 : tunevec(long expiN, GEN param)
     864             : {
     865        4389 :   GEN T = ecpp_param_get_tune(param);
     866        4389 :   long i, n = lg(T)-1;
     867        4389 :   for (i = 1; i < n; i++) { GEN x = gel(T,i); if (expiN <= x[4]) return x; }
     868        1631 :   return gel(T,n);
     869             : }
     870             : static long
     871        4389 : tunevec_tdivbd(long expiN, GEN param) { return uel(tunevec(expiN, param), 3); }
     872             : static long
     873         980 : tunevec_batchsize(long expiN, GEN param)
     874             : {
     875         980 :   long t, b = 28 - tunevec_tdivbd(expiN, param);
     876         980 :   if (b < 0) return expiN;
     877         980 :   t = expiN >> b; return t < 1? 1: t;
     878             : }
     879             : 
     880             : static GEN
     881        3409 : Dmbatch_factor_Dmqvec(GEN N, long expiN, GEN* X0, GEN Dmbatch, GEN param)
     882             : {
     883        3409 :   pari_sp av = avma;
     884             :   pari_timer ti;
     885        3409 :   GEN Dmqvec, primorial_vec = ecpp_param_get_primorial_vec(param);
     886        3409 :   GEN primorial = gel(primorial_vec, tunevec_tdivbd(expiN, param)-19);
     887             : 
     888             :   /* B1: Factor by batch */
     889        3409 :   dbg_mode() timer_start(&ti);
     890        3409 :   Dmqvec = Dmvec_batchfactor(Dmbatch, primorial);
     891        3409 :   dbg_mode() timer_record(X0, "B1", &ti);
     892             : 
     893             :   /* B2: For each batch, remove cardinalities lower than (N^(1/4)+1)^2
     894             :    *     and cardinalities in which we didn't win enough bits */
     895        3409 :   Dmqvec = Dmqvec_slice(N, Dmqvec);
     896        3409 :   if (!Dmqvec) return gc_NULL(av); /* nothing is left */
     897        3409 :   return gerepilecopy(av, Dmqvec);
     898             : }
     899             : 
     900             : /* Dmq (where q has no small prime factors): determine if q is pseudoprime */
     901             : static long
     902       19250 : Dmq_isgoodq(GEN Dmq, GEN* X0)
     903             : {
     904             :   pari_timer ti;
     905             :   long s;
     906             :   /* B3: Check pseudoprimality of each q on the list. */
     907       19250 :   dbg_mode() timer_start(&ti);
     908       19250 :   s = BPSW_psp_nosmalldiv(gel(Dmq,3));
     909       19250 :   dbg_mode() timer_record(X0, "B3", &ti);
     910       19250 :   return s; /* did not find for this m */
     911             : }
     912             : static GEN
     913         924 : mkNDmqg(GEN z, GEN N, GEN Dmq, GEN g, GEN sqrtlist)
     914             : {
     915         924 :   GEN Dinfo = gel(Dmq,1), sq =  Dfac_to_roots(Dinfo_get_Dfac(Dinfo),sqrtlist);
     916         924 :   GEN NDmqg = mkcol6(N, Dinfo, gel(Dmq,2), gel(Dmq,3), g, sq);
     917         924 :   return mkvec2(NDmqg, z);
     918             : }
     919             : /* expi(N) > 64. Return [ NDmqg, N_downrun(q) ]; NDmqg is a vector of the form
     920             :  * [N,D,m,q,g,sqrt]. For successive D, find a square root of D mod N (g is a
     921             :  * quadratic non-residue), solve U^2 + |D|V^2 = 4N, then find candidates for m.
     922             :  * When enough m's, batch-factor them to produce the q's. If one of the q's is
     923             :  * pseudoprime, recursive call with N = q. May return gen_0 at toplevel
     924             :  * => N has a small prime divisor */
     925             : static GEN
     926         980 : N_downrun(GEN N, GEN param, GEN *X0, long *depth, long persevere)
     927             : {
     928             :   pari_timer T, ti;
     929         980 :   long lgdisclist, lprimelist, t, i, j, expiN = expi(N);
     930         980 :   long persevere_next = 0, FAIL = 0;
     931             :   ulong maxpcdg;
     932             :   GEN primelist, disclist, sqrtlist, g, Dmbatch, kroP;
     933             : 
     934         980 :   dbg_mode() timer_start(&T);
     935         980 :   (*depth)++; /* we're going down the tree. */
     936         980 :   maxpcdg = ecpp_param_get_maxpcdg(param);
     937         980 :   primelist = ecpp_param_get_primelist(param);
     938         980 :   disclist = ecpp_param_get_disclist(param);
     939         980 :   lgdisclist = lg(disclist);
     940         980 :   t = tunevec_batchsize(expiN, param);
     941             : 
     942             :   /* Precomputation for taking square roots, g needed for Fp_sqrt_i */
     943         980 :   g = Fp_2gener(N); if (!g) return gen_0; /* Composite if this happens. */
     944             : 
     945             :   /* Initialize sqrtlist for this N. */
     946         980 :   lprimelist = lg(primelist);
     947         980 :   sqrtlist = zerovec(lprimelist-1);
     948             : 
     949             :   /* A2: Check (p*|N) = 1 for all p */
     950         980 :   dbg_mode() timer_start(&ti);
     951             :   /* kroP[i] will contain (primelist[i] | N) */
     952         980 :   kroP = const_vecsmall(lprimelist-1, 2/*sentinel*/);
     953         980 :   switch(mod8(N))
     954             :   { /* primelist[1,2,3] = [8,-4,-8]; N is odd */
     955         280 :     case 1: kroP[1] = kroP[2] = kroP[3] = 1; break;
     956         245 :     case 3: kroP[1] = -1; kroP[2] =-1; kroP[3] = 1; break;
     957         294 :     case 5: kroP[1] = -1; kroP[2] = 1; kroP[3] =-1; break;
     958         161 :     case 7: kroP[1] =  1; kroP[2] =-1; kroP[3] =-1; break;
     959             :   }
     960         980 :   dbg_mode() timer_record(X0, "A2", &ti);
     961             : 
     962             :   /* Print the start of this iteration. */
     963         980 :   dbg_mode() {
     964           0 :     char o = persevere? '<': '[';
     965           0 :     char c = persevere? '>': ']';
     966           0 :     err_printf(ANSI_BRIGHT_CYAN "\n%c %3d | %4ld bits%c "
     967             :                ANSI_RESET, o, *depth, expiN, c);
     968             :   }
     969             :   /* Initialize Dmbatch, populated with candidate cardinalities in Phase I
     970             :    * (until #Dmbatch >= t); its elements will be factored on Phase II */
     971         980 :   Dmbatch = vectrunc_init(t+7);
     972             : 
     973             :   /* Number of cardinalities so far; should always be equal to lg(Dmbatch)-1. */
     974             :   /* i determines which discriminant we are considering. */
     975         980 :   i = 1;
     976        4438 :   while (!FAIL)
     977             :   {
     978             :     pari_timer F;
     979        3451 :     long lgDmqlist, last_i = i, numcard = 0; /* #Dmbatch */
     980             :     GEN Dmqlist;
     981             : 
     982             :     /* Dmbatch begins "empty", but keep the allocated memory. */
     983        3451 :     setlg(Dmbatch, 1);
     984             : 
     985             :     /* PHASE I: Go through the D's and search for candidate m's.
     986             :      * Fill up Dmbatch until there are at least t elements. */
     987      539098 :     while (i < lgdisclist )
     988             :     {
     989      535619 :       GEN Dinfo = gel(disclist, i);
     990             :       long n;
     991      535619 :       if (!persevere && Dinfo_get_pd(Dinfo) > maxpcdg) { FAIL = 1; break; }
     992      535584 :       n = D_collectcards(N,param, X0, Dinfo, sqrtlist, g, Dmbatch, kroP);
     993      536508 :       if (n < 0) return gen_0;
     994      535584 :       last_i = i++;
     995      535584 :       numcard += n; if (numcard >= t) break;
     996             :     }
     997             : 
     998             :     /* We have exhausted disclist and there are no card. to be factored */
     999        3451 :     if (numcard <= 0 && (FAIL || i >= lgdisclist)) break;
    1000             : 
    1001             :     /* PHASE II: Find the corresponding q's for the m's found. Use Dmbatch */
    1002             :     /* Find coresponding q of the candidate m's. */
    1003        3409 :     dbg_mode() timer_start(&F);
    1004        3409 :     Dmqlist = Dmbatch_factor_Dmqvec(N, expiN, X0, Dmbatch, param);
    1005        3409 :     if (Dmqlist == NULL) continue; /* none left => next discriminant. */
    1006             : 
    1007             :     /* If we are cheating by adding class numbers, sort by class number */
    1008        3409 :     if (Dinfo_get_pd(gel(disclist, last_i)) > maxpcdg)
    1009         266 :       gen_sort_inplace(Dmqlist, NULL, &sort_Dmq_by_cnum, NULL);
    1010             : 
    1011             :     /* Check pseudoprimality before going down. */
    1012        3409 :     lgDmqlist = lg(Dmqlist);
    1013       21735 :     for (j = 1; j < lgDmqlist; j++)
    1014             :     {
    1015       19250 :       GEN z, Dmq = gel(Dmqlist,j), q = gel(Dmq,3);
    1016       19250 :       dbg_mode() err_printf(ANSI_WHITE "." ANSI_RESET);
    1017       19250 :       if (!Dmq_isgoodq(Dmq, X0)) continue;
    1018             : 
    1019         959 :       dbg_mode() {
    1020           0 :         err_printf(ANSI_BRIGHT_RED "\n       %5ld bits " ANSI_RESET,
    1021           0 :                    expi(q)-expiN);
    1022           0 :         err_printf("  D = %ld", Dmq_get_D(Dmq));
    1023           0 :         err_printf(ANSI_BRIGHT_BLUE "  h = %ld" ANSI_RESET, Dmq_get_h(Dmq));
    1024             :       }
    1025             :       /* q is pseudoprime */
    1026         959 :       if (expi(q) < 64) z = gen_1; /* q is prime; sentinel */
    1027             :       else
    1028             :       {
    1029         917 :         z = N_downrun(q, param, X0, depth, persevere_next);
    1030         917 :         if (!z) {
    1031          35 :           dbg_mode() { char o = persevere? '<': '[', c = persevere? '>': ']';
    1032           0 :                        err_printf(ANSI_CYAN "\n%c %3d | %4ld bits%c "
    1033             :                                   ANSI_RESET, o, *depth, expiN, c); }
    1034          35 :           continue; /* downrun failed */
    1035             :         }
    1036             :       }
    1037         924 :       return mkNDmqg(z, N, Dmq, g, sqrtlist); /* SUCCESS */
    1038             :     }
    1039        2485 :     if (i >= lgdisclist) break; /* discriminants exhausted: FAIL */
    1040        2478 :     if (Dinfo_get_pd(gel(disclist, last_i)) > maxpcdg)
    1041             :     {
    1042         252 :       dbg_mode() err_printf(ANSI_BRIGHT_RED "  !" ANSI_RESET);
    1043         252 :       persevere_next = 1;
    1044             :     }
    1045             :   }
    1046             :   /* FAILED: Out of discriminants. */
    1047          56 :   umael(*X0, 3, 1)++; /* FAILS++ */
    1048          56 :   dbg_mode() err_printf(ANSI_BRIGHT_RED "  X" ANSI_RESET);
    1049          56 :   (*depth)--; return NULL;
    1050             : }
    1051             : 
    1052             : /* x: the output of the (recursive) downrun function; return a vector
    1053             :  * whose components are [N, D, m, q, g] */
    1054             : static GEN
    1055          42 : ecpp_flattencert(GEN x, long depth)
    1056             : {
    1057          42 :   long i, l = depth+1;
    1058          42 :   GEN ret = cgetg(l, t_VEC);
    1059          42 :   if (typ(x) != t_VEC) return gen_0;
    1060          42 :   for (i = 1; i < l; i++) { gel(ret, i) = gel(x,1); x = gel(x,2); }
    1061          42 :   return ret;
    1062             : }
    1063             : 
    1064             : /* Call the first downrun then unravel the skeleton of the certificate.
    1065             :  * Assume expi(N) < 64. This returns one of the following:
    1066             :  * - a vector of the form [N, D, m, q, y]
    1067             :  * - gen_0 (if N is composite)
    1068             :  * - NULL (if FAIL) */
    1069             : static GEN
    1070          63 : ecpp_step1(GEN N, GEN param, GEN* X0)
    1071             : {
    1072          63 :   pari_sp av = avma;
    1073          63 :   long depth = 0;
    1074          63 :   GEN downrun = N_downrun(N, param, X0, &depth, 1);
    1075          63 :   if (downrun == NULL) return gc_NULL(av);
    1076          42 :   return gerepilecopy(av, ecpp_flattencert(downrun, depth));
    1077             : }
    1078             : 
    1079             : /* The input is an integer N.
    1080             :    It uses the precomputation step0 done in ecpp_step0. */
    1081             : static GEN
    1082          63 : ecpp0(GEN N, GEN param)
    1083             : {
    1084             : 
    1085             :   GEN step1, answer, Tv, Cv, X0;
    1086          63 :   pari_sp av = avma;
    1087             :   long i, j;
    1088             : 
    1089             :   /* Check if N is pseudoprime to begin with. */
    1090          63 :   if (!BPSW_psp(N)) return gen_0;
    1091             : 
    1092             :   /* Check if we should even prove it. */
    1093          63 :   if (expi(N) < 64) return N;
    1094             : 
    1095             :   /* Timers and Counters */
    1096          63 :   Tv = mkvec4(zero_zv(5), zero_zv(3), zero_zv(3), zero_zv(1));
    1097          63 :   Cv = mkvec4(zero_zv(5), zero_zv(3), zero_zv(3), zero_zv(1));
    1098          63 :   X0 = mkvec3(Tv, Cv, zero_zv(1));
    1099             : 
    1100          63 :   step1 = ecpp_step1(N, param, &X0);
    1101          63 :   if (step1 == NULL) return gc_NULL(av);
    1102          42 :   if (typ(step1) != t_VEC) { set_avma(av); return gen_0; }
    1103             : 
    1104          42 :   answer = ecpp_step2(step1, &X0, ecpp_param_get_primelist(param));
    1105             : 
    1106          42 :   dbg_mode() {
    1107           0 :   for (i = 1; i < lg(Tv); i++)
    1108             :   {
    1109           0 :     GEN v = gel(Tv, i);
    1110           0 :     long l = lg(v);
    1111           0 :     for (j = 1; j < l; j++)
    1112             :     {
    1113           0 :       long t = umael3(X0,1, i,j), n = umael3(X0,2, i,j);
    1114           0 :       if (!t) continue;
    1115           0 :       err_printf("\n  %c%ld: %16ld %16ld %16.3f", 'A'+i-1,j, t,n, (double)t/n);
    1116             :     }
    1117             :   }
    1118           0 :     err_printf("\n" ANSI_BRIGHT_RED "\nFAILS: %16ld" ANSI_RESET "\n", umael(X0, 3, 1));
    1119             :   }
    1120          42 :   return gerepilecopy(av, answer);
    1121             : }
    1122             : 
    1123             : static const long ecpp_tune[][4]=
    1124             : {
    1125             :   {  100,  2,  20,   500 },
    1126             :   {  350, 14,  21,  1000 },
    1127             :   {  450, 18,  21,  1500 },
    1128             :   {  750, 22,  22,  2000 },
    1129             :   {  900, 26,  23,  2500 },
    1130             :   { 1000, 32,  23,  3000 },
    1131             :   { 1200, 38,  24,  3500 },
    1132             :   { 1400, 44,  24,  4000 },
    1133             :   { 1600, 50,  24,  4500 },
    1134             :   { 1800, 56,  25,  5000 },
    1135             :   { 2000, 62,  25,  5500 },
    1136             :   { 2200, 68,  26,  6000 },
    1137             :   { 2400, 74,  26,  6500 },
    1138             :   { 2600, 80,  26,  7000 },
    1139             :   { 2800, 86,  26,  7500 },
    1140             :   { 3000, 92,  27,  8000 },
    1141             :   { 3200, 98,  27,  8500 },
    1142             :   { 3400, 104, 28,  9000 },
    1143             :   { 3600, 110, 28,  9500 },
    1144             :   { 3800, 116, 29, 10000 },
    1145             :   { 4000, 122, 29,     0 }
    1146             : };
    1147             : 
    1148             : /* assume N BPSW-pseudoprime */
    1149             : GEN
    1150          84 : ecpp(GEN N)
    1151             : {
    1152             :   long expiN, i, tunelen;
    1153             :   GEN tune;
    1154             : 
    1155             :   /* Check if we should even prove it. */
    1156          84 :   expiN = expi(N);
    1157          84 :   if (expiN < 64) return N;
    1158             : 
    1159          42 :   tunelen = (expiN+499)/500;
    1160          42 :   tune = cgetg(tunelen+1, t_VEC);
    1161         112 :   for (i = 1; i <= tunelen && ecpp_tune[i-1][3]; i++)
    1162         210 :     gel(tune,i) = mkvecsmall4(ecpp_tune[i-1][0], ecpp_tune[i-1][1],
    1163         140 :                               ecpp_tune[i-1][2], ecpp_tune[i-1][3]);
    1164          42 :   for (; i <= tunelen; i++) gel(tune,i) = mkvecsmall4(200*(i-1),6*i-4,30,500*i);
    1165             :   for(;;)
    1166          21 :   {
    1167          63 :     GEN C, param, x = gel(tune, tunelen);
    1168             :     pari_timer T;
    1169          63 :     dbg_mode() timer_start(&T);
    1170          63 :     param = ecpp_param_set(tune, x);
    1171          63 :     dbg_mode() {
    1172           0 :       err_printf(ANSI_BRIGHT_WHITE "\n%Ps" ANSI_RESET, x);
    1173           0 :       err_printf(ANSI_WHITE "  %8ld" ANSI_RESET, timer_delay(&T));
    1174             :     }
    1175          63 :     if ((C = ecpp0(N, param))) return C;
    1176          21 :     x[1] *= 2;
    1177          21 :     x[2] *= 2;
    1178          21 :     x[3] = minss(x[3]+1, 30);
    1179             :   }
    1180             : }
    1181             : long
    1182          42 : isprimeECPP(GEN N)
    1183             : {
    1184          42 :   pari_sp av = avma;
    1185          42 :   if (!BPSW_psp(N)) return 0;
    1186          28 :   return gc_long(av, !isintzero(ecpp(N)));
    1187             : }
    1188             : 
    1189             : /* PARI ECPP Certificate -> Human-readable format */
    1190             : static GEN
    1191           7 : cert_out(GEN x)
    1192             : {
    1193           7 :   long i, l = lg(x);
    1194             :   pari_str str;
    1195           7 :   str_init(&str, 1);
    1196           7 :   if (typ(x) == t_INT)
    1197             :   {
    1198           0 :     str_printf(&str, "%Ps is prime.\nIndeed, ispseudoprime(%Ps) = 1 and %Ps < 2^64.\n", x, x, x);
    1199           0 :     return strtoGENstr(str.string);
    1200             :   }
    1201          21 :   for (i = 1; i < l; i++)
    1202             :   {
    1203          14 :     GEN xi = gel(x, i);
    1204          14 :     str_printf(&str, "[%ld]\n", i);
    1205          14 :     str_printf(&str, " N = %Ps\n", cert_get_N(xi));
    1206          14 :     str_printf(&str, " t = %Ps\n", cert_get_t(xi));
    1207          14 :     str_printf(&str, " s = %Ps\n", cert_get_s(xi));
    1208          14 :     str_printf(&str, "a4 = %Ps\n", cert_get_a4(xi));
    1209          14 :     str_printf(&str, "D = %Ps\n", cert_get_D(xi));
    1210          14 :     str_printf(&str, "m = %Ps\n", cert_get_m(xi));
    1211          14 :     str_printf(&str, "q = %Ps\n", cert_get_q(xi));
    1212          14 :     str_printf(&str, "E = %Ps\n", cert_get_E(xi));
    1213          14 :     str_printf(&str, "P = %Ps\n", cert_get_P(xi));
    1214             :   }
    1215           7 :   return strtoGENstr(str.string);
    1216             : }
    1217             : 
    1218             : /* PARI ECPP Certificate -> Magma Certificate
    1219             :  * [* [*
    1220             :  *     N, |D|, -1, m,
    1221             :  *     [* a, b *],
    1222             :  *     [* x, y, 1 *],
    1223             :  *     [* [* q, 1 *] *]
    1224             :  *   *], ... *] */
    1225             : static GEN
    1226           7 : magma_out(GEN x)
    1227             : {
    1228           7 :   long i, n = lg(x)-1;
    1229             :   pari_str ret;
    1230           7 :   str_init(&ret, 1);
    1231           7 :   if (typ(x) == t_INT)
    1232             :   {
    1233           0 :     str_printf(&ret, "Operation not supported.");
    1234           0 :     return strtoGENstr(ret.string);
    1235             :   }
    1236           7 :   str_printf(&ret, "[* ");
    1237          21 :   for (i = 1; i<=n; i++)
    1238             :   {
    1239          14 :     GEN xi = gel(x,i), E = cert_get_E(xi), P = cert_get_P(xi);
    1240          14 :     str_printf(&ret, "[* %Ps, %Ps, -1, %Ps, ",
    1241             :               cert_get_N(xi), negi(cert_get_D(xi)), cert_get_m(xi));
    1242          14 :     str_printf(&ret, "[* %Ps, %Ps *], ", gel(E,1), gel(E,2));
    1243          14 :     str_printf(&ret, "[* %Ps, %Ps, 1 *], ", gel(P,1), gel(P,2));
    1244          14 :     str_printf(&ret, "[* [* %Ps, 1 *] *] *]", cert_get_q(xi));
    1245          14 :     if (i != n) str_printf(&ret, ", ");
    1246             :   }
    1247           7 :   str_printf(&ret, " *]");
    1248           7 :   return strtoGENstr(ret.string);
    1249             : }
    1250             : 
    1251             : /* Prints: label=(sign)hexvalue\n
    1252             :  *   where sign is + or -
    1253             :  *   hexvalue is value in hex, of the form: 0xABC123 */
    1254             : static void
    1255          70 : primo_printval(pari_str *ret, const char* label, GEN value)
    1256             : {
    1257          70 :   str_printf(ret, label);
    1258          70 :   if (signe(value) >= 0) str_printf(ret, "=0x%P0X\n", value);
    1259          14 :   else str_printf(ret, "=-0x%P0X\n", negi(value));
    1260          70 : }
    1261             : 
    1262             : /* Input: PARI ECPP Certificate / Output: PRIMO Certificate
    1263             :  *
    1264             :  * Let Y^2 = X^3 + aX + b be the elliptic curve over N with the point (x,y)
    1265             :  * as described in the PARI certificate.
    1266             :  *
    1267             :  * If J != 0, 1728, PRIMO asks for [J,T], where T = a/A * B/b * x,
    1268             :  * A = 3J(1728-J) and B = 2J(1728-J)^2.
    1269             :  *
    1270             :  * If J=0 or J=1728, PRIMO asks for [A,B,T]; we let A=a and B=b => T = x */
    1271             : static GEN
    1272           7 : primo_out(GEN x)
    1273             : {
    1274           7 :   long i, l = (typ(x) == t_INT) ? 1 : lg(x);
    1275             :   pari_str ret;
    1276           7 :   str_init(&ret, 1);
    1277           7 :   str_printf(&ret, "[PRIMO - Primality Certificate]\nFormat=4\n");
    1278           7 :   str_printf(&ret, "TestCount=%ld\n\n[Comments]\n", l-1);
    1279           7 :   str_printf(&ret, "Generated by %s\n",paricfg_version);
    1280           7 :   str_printf(&ret, "https://pari.math.u-bordeaux.fr/\n\n[Candidate]\n");
    1281           7 :   if (typ(x) == t_INT)
    1282             :   {
    1283           0 :     str_printf(&ret, "N=0x%P0X\n", x);
    1284           0 :     return strtoGENstr(ret.string);
    1285             :   }
    1286           7 :   str_printf(&ret, "N=0x%P0X\n\n", cert_get_N(gel(x,1)));
    1287          21 :   for (i = 1; i < l; i++)
    1288             :   {
    1289          14 :     GEN a4, a6, N, Nover2, xi = gel(x,i);
    1290             :     long Ais0, Bis0;
    1291          14 :     str_printf(&ret, "[%ld]\n", i);
    1292          14 :     N = cert_get_N(xi);
    1293          14 :     Nover2 = shifti(N, -1);
    1294          14 :     primo_printval(&ret, "S", cert_get_s(xi));
    1295          14 :     primo_printval(&ret, "W", cert_get_t(xi));
    1296          14 :     a4 = cert_get_a4(xi); Ais0 = isintzero(a4);
    1297          14 :     a6 = cert_get_a6(xi); Bis0 = isintzero(a6);
    1298          14 :     if (!Ais0 && !Bis0) { /* J != 0, 1728 */
    1299           0 :       primo_printval(&ret, "J", Fp_center_i(cert_get_J(xi), N, Nover2));
    1300           0 :       primo_printval(&ret, "T", cert_get_T(xi));
    1301             :     } else {
    1302          14 :       if (!Ais0) a4 = Fp_center_i(a4, N, Nover2);
    1303          14 :       if (!Bis0) a6 = Fp_center_i(a6, N, Nover2);
    1304          14 :       primo_printval(&ret, "A", a4);
    1305          14 :       primo_printval(&ret, "B", a6);
    1306          14 :       primo_printval(&ret, "T", cert_get_x(xi));
    1307             :     }
    1308          14 :     str_printf(&ret, "\n");
    1309             :   }
    1310           7 :   return strtoGENstr(ret.string);
    1311             : }
    1312             : 
    1313             : /* return 1 if q > (N^1/4 + 1)^2, 0 otherwise */
    1314             : static long
    1315         322 : Nq_isvalid(GEN N, GEN q)
    1316             : {
    1317         322 :   GEN c = subii(sqri(subiu(q,1)), N);
    1318         322 :   if (signe(c) <= 0) return 0;
    1319             :   /* (q-1)^2 > N; check that (N - (q-1)^2)^2 > 16Nq */
    1320         322 :   return (cmpii(sqri(c), shifti(mulii(N,q), 4)) > 0);
    1321             : }
    1322             : 
    1323             : /* return 1 if 'cert' is a valid PARI ECPP certificate */
    1324             : static long
    1325          28 : ecppisvalid_i(GEN cert)
    1326             : {
    1327          28 :   const long trustbits = 64;/* a pseudoprime < 2^trustbits is prime */
    1328          28 :   long i, lgcert = lg(cert);
    1329          28 :   GEN ql, q = gen_0;
    1330             : 
    1331          28 :   if (typ(cert) == t_INT)
    1332             :   {
    1333           0 :     if (expi(cert) >= trustbits) return 0; /* >= 2^trustbits */
    1334           0 :     return BPSW_psp(cert);
    1335             :   }
    1336          28 :   if (typ(cert) != t_VEC || lgcert <= 1) return 0;
    1337          28 :   ql = gel(cert, lgcert-1); if (lg(ql)-1 != 5) return 0;
    1338          28 :   ql = cert_get_q(ql);
    1339          28 :   if (expi(ql) >= trustbits || !BPSW_psp(ql)) return 0;
    1340             : 
    1341         700 :   for (i = 1; i < lgcert; i++)
    1342             :   {
    1343         329 :     GEN N, W, mP, sP, r, m, s, a, P, certi = gel(cert, i);
    1344         336 :     if (lg(certi)-1 != 5) return 0;
    1345             : 
    1346         329 :     N = cert_get_N(certi);
    1347         329 :     if (typ(N) != t_INT || signe(N) <= 0) return 0;
    1348         329 :     switch(umodiu(N,6))
    1349             :     {
    1350         322 :       case 1: case 5: break; /* Check if N is not divisible by 2 or 3 */
    1351           7 :       default: return 0;
    1352             :     }
    1353             :     /* Check if N matches the q from the previous entry. */
    1354         322 :     if (i > 1 && !equalii(N, q)) return 0;
    1355             : 
    1356         322 :     W = cert_get_t(certi);
    1357         322 :     if (typ(W) != t_INT) return 0;
    1358             :     /* Check if W^2 < 4N (Hasse interval) */
    1359         322 :     if (!(cmpii(sqri(W), shifti(N,2)) < 0)) return 0;
    1360             : 
    1361         322 :     s = cert_get_s(certi);
    1362         322 :     if (typ(s) != t_INT || signe(s) < 0) return 0;
    1363             : 
    1364         322 :     m = cert_get_m(certi);
    1365         322 :     q = dvmdii(m, s, &r);
    1366             :     /* Check m%s == 0 */
    1367         322 :     if (!isintzero(r)) return 0;
    1368             :     /* Check q > (N^(1/4) + 1)^2 */
    1369         322 :     if (!Nq_isvalid(N, q)) return 0;
    1370             : 
    1371         322 :     a = cert_get_a4(certi);
    1372         322 :     if (typ(a) != t_INT) return 0;
    1373             : 
    1374         322 :     P = cert_get_P(certi);
    1375         322 :     if (typ(P) != t_VEC || lg(P) != 3) return 0;
    1376         322 :     P = FpE_to_FpJ(P);
    1377             : 
    1378             :     /* Check mP == 0 */
    1379         322 :     mP = FpJ_mul(P, m, a, N);
    1380         322 :     if (!FpJ_is_inf(mP)) return 0;
    1381             : 
    1382             :     /* Check sP != 0 and third component is coprime to N */
    1383         322 :     sP = FpJ_mul(P, s, a, N);
    1384         322 :     if (!isint1(gcdii(gel(sP, 3), N))) return 0;
    1385             :   }
    1386          21 :   return 1;
    1387             : }
    1388             : long
    1389          28 : ecppisvalid(GEN cert)
    1390          28 : { pari_sp av = avma; return gc_long(av, ecppisvalid_i(cert)); }
    1391             : 
    1392             : GEN
    1393          21 : ecppexport(GEN cert, long flag)
    1394             : {
    1395          21 :   switch(flag)
    1396             :   {
    1397           7 :     case 0: return cert_out(cert);
    1398           7 :     case 1: return primo_out(cert);
    1399           7 :     case 2: return magma_out(cert);
    1400             :   }
    1401           0 :   pari_err_FLAG("primecertexport");
    1402             :   return NULL;/*LCOV_EXCL_LINE*/
    1403             : }

Generated by: LCOV version 1.13