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

Generated by: LCOV version 1.13