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

Generated by: LCOV version 1.13