Code coverage tests

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

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

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

LCOV - code coverage report
Current view: top level - basemath - ecpp.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.10.0 lcov report (development 21059-cbe0d6a) Lines: 886 958 92.5 %
Date: 2017-09-22 06:24:58 Functions: 98 104 94.2 %
Legend: Lines: hit not hit

          Line data    Source code
       1             : /* Copyright (C) 2014 The PARI group.
       2             : 
       3             : This file is part of the PARI/GP package.
       4             : 
       5             : PARI/GP is free software; you can redistribute it and/or modify it under the
       6             : terms of the GNU General Public License as published by the Free Software
       7             : Foundation. It is distributed in the hope that it will be useful, but WITHOUT
       8             : ANY WARRANTY WHATSOEVER.
       9             : 
      10             : Check the License for details. You should have received a copy of it, along
      11             : with the package; see the file 'COPYING'. If not, write to the Free Software
      12             : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
      13             : 
      14             : #include "pari.h"
      15             : #include "paripriv.h"
      16             : 
      17             : #define dbg_mode() if (DEBUGLEVEL >= 2)
      18             : 
      19             : #define ANSI_COLOR_RED     "\x1b[31m"
      20             : #define ANSI_COLOR_GREEN   "\x1b[32m"
      21             : #define ANSI_COLOR_YELLOW  "\x1b[33m"
      22             : #define ANSI_COLOR_BLUE    "\x1b[34m"
      23             : #define ANSI_COLOR_MAGENTA "\x1b[35m"
      24             : #define ANSI_COLOR_CYAN    "\x1b[36m"
      25             : #define ANSI_COLOR_WHITE   "\x1b[37m"
      26             : 
      27             : #define ANSI_COLOR_BRIGHT_RED     "\x1b[31;1m"
      28             : #define ANSI_COLOR_BRIGHT_GREEN   "\x1b[32;1m"
      29             : #define ANSI_COLOR_BRIGHT_YELLOW  "\x1b[33;1m"
      30             : #define ANSI_COLOR_BRIGHT_BLUE    "\x1b[34;1m"
      31             : #define ANSI_COLOR_BRIGHT_MAGENTA "\x1b[35;1m"
      32             : #define ANSI_COLOR_BRIGHT_CYAN    "\x1b[36;1m"
      33             : #define ANSI_COLOR_BRIGHT_WHITE   "\x1b[37;1m"
      34             : 
      35             : #define ANSI_COLOR_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        2702 : zv_binsearch0(void *E, long (*f)(void* E, GEN x), GEN x)
      46             : {
      47        2702 :   long lo = 1;
      48        2702 :   long hi = lg(x)-1;
      49       15722 :   while (lo < hi) {
      50       10318 :     long mi = lo + (hi - lo)/2 + 1;
      51       10318 :     if (f(E,gel(x,mi))) lo = mi;
      52        6146 :     else hi = mi - 1;
      53             :   }
      54        2702 :   if (f(E,gel(x,lo))) return lo;
      55        1330 :   return 0;
      56             : }
      57             : 
      58             : INLINE long
      59           0 : timer_record(GEN* X0, const char* Xx, pari_timer* ti, long c)
      60             : {
      61           0 :   long t = 0;
      62           0 :   dbg_mode() {
      63             :     long i, j;
      64           0 :     t = timer_delay(ti);
      65           0 :     if (!X0) return t;
      66           0 :     i = Xx[0]-'A'+1;
      67           0 :     j = Xx[1]-'1'+1;
      68           0 :     umael3(*X0, 1, i, j) += t;
      69           0 :     umael3(*X0, 2, i, j) += c;
      70             :   }
      71           0 :   return t;
      72             : }
      73             : 
      74             : INLINE long
      75        3701 : FpJ_is_inf(GEN P)
      76             : {
      77        3701 :   return signe(gel(P, 3)) == 0;
      78             : }
      79             : 
      80             : /* n > 1, D = divisors(n); sets L = 2*lambda(n), S = sigma(n) */
      81             : static void
      82    13999930 : lamsig(GEN D, long *pL, long *pS)
      83             : {
      84    13999930 :   pari_sp av = avma;
      85    13999930 :   long i, l = lg(D), L = 1, S = D[l-1]+1;
      86    54426932 :   for (i = 2; i < l; ++i) /* skip d = 1 */
      87             :   {
      88    54426932 :     long d = D[i], nd = D[l-i]; /* nd = n/d */
      89    54426932 :     if (d < nd) { L += d; S += d + nd; }
      90             :     else
      91             :     {
      92    13999930 :       L <<= 1; if (d == nd) { L += d; S += d; }
      93    13999930 :       break;
      94             :     }
      95             :   }
      96    13999930 :   avma = av; *pL = L; *pS = S;
      97    13999930 : }
      98             : 
      99             : static GEN
     100          35 : allhclassno2(long lim)
     101             : {
     102          35 :   pari_sp av = avma;
     103          35 :   GEN VHDH, CACHE = NULL;
     104          35 :   const long cachestep = 1000; /* don't increase this: RAM cache thrashing */
     105          35 :   long r, N, cachea, cacheb, lim0 = 2, LIM = lim0 << 1;
     106             : 
     107          35 :   if (lim <= 0) lim = 5;
     108          35 :   r = lim&3L; if (r) lim += 4-r;
     109          35 :   VHDH = cgetg_block(lim/2 + 1, t_VECSMALL);
     110          35 :   VHDH[1] = 2;
     111          35 :   VHDH[2] = 3;
     112          35 :   for (N = 3; N <= lim0; N++) VHDH[N] = VHDH[N];
     113          35 :   cachea = cacheb = 0;
     114     7000000 :   for (N = LIM + 3; N <= lim; N += 4)
     115             :   {
     116     6999965 :     long s = 0, limt = usqrt(N>>2), flsq = 0, ind, t, L, S;
     117             :     GEN DN, DN2;
     118             :     {
     119             :       GEN F;
     120     6999965 :       if (N + 2 > cacheb)
     121             :       { /* update local cache (recycle memory) */
     122       27965 :         cachea = N;
     123       27965 :         if (cachea + 2*cachestep > lim)
     124          35 :           cacheb = lim+2; /* fuse last 2 chunks */
     125             :         else
     126       27930 :           cacheb = cachea + cachestep;
     127       27965 :         avma = av; /* FIXME: need only factor odd integers in the range */
     128       27965 :         CACHE = vecfactoru_i(cachea, cacheb);
     129             :       }
     130             :       /* use local cache */
     131     6999965 :       F = gel(CACHE,N - cachea + 1); /* factoru(N) */
     132     6999965 :       DN = divisorsu_fact(gel(F,1), gel(F,2));
     133     6999965 :       F = gel(CACHE,N - cachea + 3); /* factoru(N+2) */
     134     6999965 :       DN2 = divisorsu_fact(gel(F,1), gel(F,2));
     135             :     }
     136     6999965 :     ind = N >> 1;
     137  2774330965 :     for (t = 1; t <= limt; ++t)
     138             :     {
     139  2767331000 :       ind -= (t<<2)-2; /* N/2 - 2t^2 */
     140  2767331000 :       if (ind) s += VHDH[ind]; else flsq = 1;
     141             :     }
     142     6999965 :     lamsig(DN, &L,&S);
     143     6999965 :     VHDH[N >> 1] = 2*S - 3*L - 2*s + flsq;
     144     6999965 :     s = 0; flsq = 0; limt = (usqrt(N+2) - 1) >> 1;
     145     6999965 :     ind = (N+1) >> 1;
     146  2770851930 :     for (t = 1; t <= limt; ++t)
     147             :     {
     148  2763851965 :       ind -= t<<2; /* (N+1)/2 - 2t(t+1) */
     149  2763851965 :       if (ind) s += VHDH[ind]; else flsq = 1;
     150             :     }
     151     6999965 :     lamsig(DN2, &L,&S);
     152     6999965 :     VHDH[(N+1) >> 1] = S - 3*(L >> 1) - s - flsq;
     153             :   }
     154          35 :   return VHDH;
     155             : }
     156             : 
     157             : /*****************************************************************/
     158             : 
     159             : /* Assuming D < 0,
     160             :    these return the number of units in Q(sqrt(D)). */
     161             : INLINE long
     162    73173317 : D_get_wD(long D)
     163             : {
     164    73173317 :   if (D == -4) return 4;
     165    73172463 :   if (D == -3) return 6;
     166    73171441 :   return 2;
     167             : }
     168             : 
     169             : /* Dinfo contains information related to the discirminant
     170             :       D: the discriminant (D < 0)
     171             :       h: the class number associated to D
     172             :      bi: the "best invariant" for computing polclass(D, bi)
     173             :          In particular, we choose this invariant to be the one
     174             :            in which D is compatible
     175             :            and the coefficients of polclass(D, bi) are minimized
     176             :      pd: the degree of polclass
     177             :          This is usually equal to h except when the invariant
     178             :            is a double eta-quotient w_p,q and p|D and q|D.
     179             :            In this special case, the degree is h/2.
     180             :    Dfac: the prime factorization of D
     181             :          Note that D can be factored as D = q0 q1* q2* ... qn*
     182             :            where q0* = 1, 4, -4, 8
     183             :              and qi* = 1 mod 4 and |qi| is a prime
     184             :          The prime factorization is implemented as a vecsmall
     185             :            listing the indices of the qi* as they appear in
     186             :            the primelist. If q0* = 1, the first component of
     187             :            this vecsmall contains the index of q1*. Otherwise,
     188             :            it contains the index of q0*.
     189             : */
     190             : INLINE long
     191    73170020 : Dinfo_get_D(GEN Dinfo)
     192    73170020 : { return gel(Dinfo, 1)[1]; }
     193             : 
     194             : INLINE long
     195    73163076 : Dinfo_get_h(GEN Dinfo)
     196    73163076 : { return gel(Dinfo, 1)[2]; }
     197             : 
     198             : INLINE long
     199    73163741 : Dinfo_get_bi(GEN Dinfo)
     200    73163741 : { return gel(Dinfo, 1)[3]; }
     201             : 
     202             : INLINE ulong
     203    73323152 : Dinfo_get_pd(GEN Dinfo)
     204    73323152 : { return umael(Dinfo, 1, 4); }
     205             : 
     206             : INLINE GEN
     207      194635 : Dinfo_get_Dfac(GEN Dinfo)
     208      194635 : { return gel(Dinfo, 2); }
     209             : 
     210             : /* primelist and indexlist
     211             : 
     212             :    primelist begins with 8, -4, -8, which we consider as "prime"
     213             :    the subsequent elements are the corresponding p^star of the
     214             :    odd primes below the indicated limit (maxsqrt) listed in
     215             :    increasing absolute value
     216             : 
     217             :    indexlist keeps the index of the odd primes. see tables below:
     218             : 
     219             :               i |   1 |   2 |   3 |   4 |   5 |   6 |   7 |   8 |   9 |  10 |  11 |
     220             :    -------------+-----+-----+-----+-----+-----+-----+-----+-----+-----+-----+-----+
     221             :    primelist[i] |   8 |  -4 |  -8 |  -3 |   5 |  -7 | -11 |  13 |  17 | -19 | -23 |
     222             : 
     223             :               i |   1 |   2 |   3 |   4 |   5 |   6 |   7 |   8 |   9 |  10 |  11 |
     224             :               p |   3 |   5 |   7 |   9 |  11 |  13 |  15 |  17 |  19 |  21 |  23 |
     225             :    -------------+-----+-----+-----+-----+-----+-----+-----+-----+-----+-----+-----+
     226             :    indexlist[i] |   4 |   5 |   6 | XXX |   7 |   8 | XXX |   9 |  10 | XXX |  11 |
     227             :    -------------+-----+-----+-----+-----+-----+-----+-----+-----+-----+-----+-----+
     228             : */
     229             : 
     230             : /*  Input: maxsqrt
     231             :    Output: primelist containing primes whose absolute value is at most maxsqrt
     232             : */
     233             : static GEN
     234          35 : ecpp_primelist_init( long maxsqrt)
     235             : {
     236             :   forprime_t T;
     237          35 :   GEN primelist = cgetg(maxsqrt, t_VECSMALL);
     238          35 :   long j = 1;
     239          35 :   u_forprime_init(&T, 3, ULONG_MAX);
     240          35 :   primelist[j] =  8; j++;
     241          35 :   primelist[j] = -4; j++;
     242          35 :   primelist[j] = -8; j++;
     243             :   while (1) {
     244        4844 :     long p = u_forprime_next(&T);
     245        4844 :     if (p > maxsqrt) break;
     246        4809 :     if (p == 0) break;
     247        4809 :     if (p%4 == 3) p = -p;
     248        4809 :     primelist[j] = p;
     249        4809 :     j++;
     250        4809 :   }
     251          35 :   setlg(primelist, j);
     252          35 :   return primelist;
     253             : }
     254             : 
     255             : /*  Input: p, indexlist
     256             :    Output: index of p in the primelist associated to indexlist
     257             :            notice the special case when p is even
     258             : */
     259             : INLINE ulong
     260    22797019 : p_to_index(long p, GEN indexlist)
     261             : {
     262    22797019 :   if (p%2 == 0)
     263             :   {
     264     7000000 :     if (p ==  8) return 1;
     265     6125014 :     if (p == -4) return 2;
     266     2625014 :     if (p == -8) return 3;
     267             :   }
     268    17547019 :   return uel(indexlist, (labs(p))/2);
     269             : }
     270             : 
     271             : /*  Input: i, primelist
     272             :    Output: returns the ith element of primelist
     273             : */
     274             : INLINE long
     275      217511 : index_to_p(long i, GEN primelist) { return primelist[i]; }
     276             : 
     277             : /*  Input: primelist
     278             :    Output: returns a vecsmall indicating at what index
     279             :            of primelist each odd prime occurs
     280             : */
     281             : INLINE GEN
     282          35 : primelist_to_indexlist(GEN primelist)
     283             : {
     284             :   long i;
     285          35 :   long lgprimelist = lg(primelist);
     286          35 :   long maxprime = labs(primelist[lgprimelist-1]);
     287             :   GEN indexlist;
     288             : 
     289          35 :   if (maxprime < 8) maxprime = 8;
     290          35 :   indexlist = zero_zv(maxprime/2);
     291        4844 :   for (i = 4; i < lgprimelist; i++)
     292             :   {
     293        4809 :     long p = labs(primelist[i]);
     294        4809 :     uel(indexlist, p/2) = i;
     295             :   }
     296          35 :   return indexlist;
     297             : }
     298             : 
     299             : INLINE ulong
     300          35 : ecpp_param_get_maxsqrt(GEN param)
     301          35 : { return umael3(param, 1, 1, 1); }
     302             : 
     303             : INLINE ulong
     304          35 : ecpp_param_get_maxdisc(GEN param)
     305          35 : { return umael3(param, 1, 1, 2); }
     306             : 
     307             : INLINE ulong
     308         686 : ecpp_param_get_maxpcdg(GEN param)
     309         686 : { return umael3(param, 1, 1, 3); }
     310             : 
     311             : INLINE GEN
     312      354746 : ecpp_param_get_primelist(GEN param)
     313      354746 : { return gmael3(param, 1, 2, 1); }
     314             : 
     315             : INLINE GEN
     316      158725 : ecpp_param_get_disclist(GEN param)
     317      158725 : { return gmael(param, 1, 3); }
     318             : 
     319             : INLINE GEN
     320        1351 : ecpp_param_get_primorial_vec(GEN param)
     321        1351 : { return gel(param, 2); }
     322             : 
     323             : /*  Input: tree obtained using ZV_producttree(v), and integer a
     324             :    Output: product of v[1] to v[i]
     325             : */
     326             : static GEN
     327          49 : producttree_find_partialprod(GEN tree, GEN v, ulong a)
     328             : {
     329          49 :   GEN b = gen_1;
     330             :   long i;
     331          49 :   long lgtree = lg(tree);
     332        1057 :   for (i = 0; i < lgtree; i++)
     333             :   {
     334        1008 :     if (a%2 != 0)
     335             :     {
     336         350 :       if (i==0) b = muliu(b, uel(v, a));
     337         301 :       else b = mulii(b, gmael(tree, i, a));
     338             :     }
     339        1008 :     a/=2;
     340             :   }
     341          49 :   return b;
     342             : }
     343             : 
     344             : /*  Input: x, 22 <= x <= 26
     345             :    Output: v, a vector whose ith component is the product of all primes below 2^(21+x)
     346             : */
     347             : static GEN
     348          35 : primorial_vec(ulong x)
     349             : {
     350          35 :   pari_sp av = avma;
     351             :   long i;
     352          35 :   GEN v = primes_upto_zv(1UL << x);
     353          35 :   GEN tree = ZV_producttree(v);
     354             :   /* ind[i] is the number such that the ind[i]th prime number is the largest prime number below 2^(21+i) */
     355          35 :   GEN ind = mkvecsmall5(295947, 564163, 1077871, 2063689, 3957809);
     356          35 :   long y = x-21;
     357          35 :   GEN ret = cgetg(y+1, t_VEC);
     358          35 :   for (i = 1; i <= y; i++) gel(ret, i) = producttree_find_partialprod(tree, v, uel(ind, i));
     359          35 :   return gerepilecopy(av, ret);
     360             : }
     361             : 
     362             : INLINE GEN
     363        1351 : ecpp_param_get_primorial(GEN param, long v)
     364             : {
     365        1351 :   GEN primorial_vec = ecpp_param_get_primorial_vec(param);
     366        1351 :   return gel(primorial_vec, v);
     367             : }
     368             : 
     369             : INLINE void
     370          35 : ecpp_param_set_maxdisc(GEN param, ulong x)
     371          35 : { umael3(param, 1, 1, 2) = x; }
     372             : 
     373             : INLINE void
     374          35 : ecpp_param_set_maxpcdg(GEN param, ulong x)
     375          35 : { umael3(param, 1, 1, 3) = x; }
     376             : 
     377             : INLINE void
     378          35 : ecpp_param_set_tdivexp(GEN param, ulong x)
     379          35 : { gel(param, 2) = primorial_vec(x); }
     380             : 
     381             : static GEN
     382             : ecpp_disclist_init( long maxsqrt, ulong maxdisc, GEN primelist);
     383             : 
     384             : INLINE void
     385          35 : ecpp_param_set_disclist(GEN param)
     386             : {
     387          35 :    long maxsqrt = ecpp_param_get_maxsqrt(param);
     388          35 :   ulong maxdisc = ecpp_param_get_maxdisc(param);
     389          35 :   GEN primelist = ecpp_param_get_primelist(param);
     390          35 :   gmael(param, 1, 3) = ecpp_disclist_init(maxsqrt, maxdisc, primelist);
     391          35 : }
     392             : 
     393             : /* NDinfomqg contains
     394             :    N, as in the theorem in ??ecpp
     395             :    Dinfo, as discussed in the comment about Dinfo
     396             :    m, as in the theorem in ??ecpp
     397             :    q, as in the theorem in ??ecpp
     398             :    g, a quadratic non-residue modulo N
     399             : */
     400             : 
     401             : INLINE GEN
     402         665 : NDinfomqg_get_N(GEN x)
     403         665 : { return gel(x,1); }
     404             : 
     405             : INLINE GEN
     406        6944 : NDinfomqg_get_Dinfo(GEN x)
     407        6944 : { return gel(x,2); }
     408             : 
     409             : INLINE GEN
     410        6279 : NDinfomqg_get_D(GEN x)
     411        6279 : { return stoi(Dinfo_get_D(NDinfomqg_get_Dinfo(x))); }
     412             : 
     413             : INLINE long
     414        5614 : NDinfomqg_get_longD(GEN x)
     415        5614 : { return itos(NDinfomqg_get_D(x)); }
     416             : 
     417             : INLINE GEN
     418         665 : NDinfomqg_get_m(GEN x)
     419         665 : { return gel(x,3); }
     420             : 
     421             : INLINE GEN
     422         665 : NDinfomqg_get_q(GEN x)
     423         665 : { return gel(x,4); }
     424             : 
     425             : INLINE GEN
     426         665 : NDinfomqg_get_g(GEN x)
     427         665 : { return gel(x,5); }
     428             : 
     429             : /* COMPARATOR FUNCTIONS */
     430             : 
     431             : static int
     432    36581538 : sort_disclist(void *data, GEN x, GEN y)
     433             : {
     434             :   long d1, h1, bi1, pd1, hf1, wD1, d2, h2, bi2, pd2, hf2, wD2;
     435    36581538 :   if (data != NULL) return itos(closure_callgen2( (GEN)data, x, y) );
     436             : 
     437    36581538 :   d1 = Dinfo_get_D(x); /* discriminant */
     438    36581538 :   h1 = Dinfo_get_h(x); /* class number */
     439    36581538 :   bi1 = Dinfo_get_bi(x); /* best invariant */
     440    36581538 :   pd1 = Dinfo_get_pd(x); /* degree of polclass */
     441    36581538 :   hf1 = modinv_height_factor(bi1); /* height factor */
     442    36581538 :   wD1 = D_get_wD(d1); /* number of units */
     443    36581538 :   d2 = Dinfo_get_D(y);
     444    36581538 :   h2 = Dinfo_get_h(y);
     445    36581538 :   bi2 = Dinfo_get_bi(y);
     446    36581538 :   pd2 = Dinfo_get_pd(y);
     447    36581538 :   hf2 = modinv_height_factor(bi2);
     448    36581538 :   wD2 = D_get_wD(d2);
     449             : 
     450             :   /* higher number of units means more elliptic curves to try */
     451    36581538 :   if (wD1 != wD2) return wD2 > wD1 ? 1 : -1;
     452             :   /* lower polclass degree is better because of faster computation of roots modulo N */
     453    36580523 :   if (pd1 != pd2) return pd1 > pd2 ? 1 : -1;
     454             :   /* lower class number is better because of higher probability of passing cornacchia step */
     455    19031698 :   if (h1 != h2) return h1 > h2 ? 1 : -1;
     456             :   /* higher height factor is better because polclass would have lower coefficients */
     457    17833144 :   if (hf1 != hf2) return hf2 > hf1 ? 1 : -1;
     458             :   /* "higher" discriminant is better since its absolute value is lower */
     459     8233953 :   if (d1 != d2) return d2 > d1 ? 1 : -1;
     460             : 
     461           0 :   return 0;
     462             : }
     463             : 
     464             : static int
     465        2807 : sort_NDmq_by_D(void *data, GEN x, GEN y)
     466             : {
     467        2807 :   long d1 = NDinfomqg_get_longD(x);
     468        2807 :   long d2 = NDinfomqg_get_longD(y);
     469        2807 :   (void)data; return d2 > d1 ? 1 : -1;
     470             : }
     471             : 
     472             : static int
     473       68614 : sort_Dmq_by_q(void *data, GEN x, GEN y)
     474             : {
     475       68614 :   GEN q1 = gel(x, 3);
     476       68614 :   GEN q2 = gel(y, 3);
     477       68614 :   (void)data; return cmpii(q1, q2);
     478             : }
     479             : 
     480             : static int
     481         938 : sort_Dmq_by_cnum(void *data, GEN x, GEN y)
     482             : {
     483         938 :   ulong h1 = umael3(x, 1, 1, 2);
     484         938 :   ulong h2 = umael3(y, 1, 1, 2);
     485         938 :   if (h1 != h2) return h1 > h2 ? 1 : -1;
     486         245 :   else return sort_Dmq_by_q(data, x, y);
     487             : }
     488             : 
     489             : static void
     490          35 : ecpp_param_set_maxsqrt(GEN param, long x)
     491             : {
     492          35 :   umael3(param, 1, 1, 1) = x;
     493          35 :   gmael3(param, 1, 2, 1) = ecpp_primelist_init(x);
     494          35 : }
     495             : 
     496             : static GEN
     497          35 : ecpp_disclist_init( long maxsqrt, ulong maxdisc, GEN primelist)
     498             : {
     499          35 :   pari_sp av = avma;
     500             :   long i, p;
     501             :   forprime_t T;
     502          35 :   GEN Harr = allhclassno2(maxdisc); /* table of class numbers*/
     503             :   GEN ev, od; /* ev: D = 0 mod 4; od: D = 1 mod 4 */
     504          35 :   GEN indexlist = primelist_to_indexlist(primelist); /* for making Dfac */
     505             :   GEN merge; /* return this */
     506          35 :   long lenv = maxdisc/4; /* max length of od/ev FIXME: ev can have length maxdisc/8 */
     507             :   long lgmerge; /* length of merge */
     508             :   long N; /* maximum number of positive prime factors */
     509          35 :   long u = 0; /* keeps track of size of merge */
     510          35 :   long u3 = 1, u4 = 1; /* keeps track of size of od/ev */
     511             : 
     512             :   /* tuning paramaters blatantly copied from vecfactoru */
     513          35 :   if (maxdisc < 510510UL) N = 7;
     514          14 :   else if (maxdisc < 9699690UL) N = 8;
     515             :   #ifdef LONG_IS_64BIT
     516           0 :     else if (maxdisc < 223092870UL) N = 9;
     517           0 :     else if (maxdisc < 6469693230UL) N = 10;
     518           0 :     else if (maxdisc < 200560490130UL) N = 11;
     519           0 :     else if (maxdisc < 7420738134810UL) N = 12;
     520           0 :     else if (maxdisc < 304250263527210UL) N = 13;
     521           0 :     else N = 16; /* don't bother */
     522             :   #else
     523           0 :     else N = 9;
     524             :   #endif
     525             : 
     526             :   /* class number computations */
     527          35 :   Harr[1] = Harr[2] = 6;
     528             : 
     529             :   /* initialization */
     530          35 :   od = cgetg(lenv + 1, t_VEC);
     531          35 :   ev = cgetg(lenv + 1, t_VEC);
     532             : 
     533             :   /* od[i] holds Dinfo of -(4*i-1)
     534             :      ev[i] holds Dinfo of -(4*i)   */
     535     7000035 :   for (i = 1; i <= lenv; i++)
     536             :   {
     537             :     long h;
     538     7000000 :     long x = 1;
     539     7000000 :     h = Harr[2*i-1]/6; /* class number of -(4*i-1) */
     540     7000000 :     gel(od, i) = mkvec2( mkvecsmall4(1, h, 0, h), vecsmalltrunc_init(N) );
     541     7000000 :     switch(i&7)
     542             :     {
     543             :       case 0:
     544     1750000 :       case 4: {x = 0;} break;
     545      875014 :       case 2: {x = -8;} break;
     546      874986 :       case 6: {x = 8;} break;
     547     3500000 :       default: {x = -4;} break;
     548             :     }
     549     7000000 :     h = Harr[2*i]/6; /* class number of -(4*i) */
     550     7000000 :     gel(ev, i) = mkvec2( mkvecsmall4(x, h, 0, h), vecsmalltrunc_init(N) );
     551     7000000 :     vecsmalltrunc_append(gmael(ev, i, 2), p_to_index(x, indexlist));
     552             :   }
     553             : 
     554             :   /* sieve part */
     555          35 :   u_forprime_init(&T, 3, maxsqrt);
     556        4879 :   while ( (p = u_forprime_next(&T)) )
     557             :   {
     558        4809 :     long s = 0; /* sp = number we're looking at, detects nonsquarefree numbers */
     559        4809 :     long t = 0; /* points to which number we're looking at */
     560        4809 :     long q = 1; /* p^star */
     561        4809 :     long tod = 0; /* starting t for od */
     562        4809 :     long tev = 0; /* starting t for ev */
     563        4809 :     switch(p&3)
     564             :     {
     565        2296 :       case 1: {q = p; s = 3; tod = (s*p+1)/4; tev = p;} break;
     566        2513 :       case 3: {q = -p; s = 1; tod = (s*p+1)/4; tev = p;} break;
     567             :     }
     568    11984987 :     for (t = tod; t <= lenv; t += p)
     569             :     {
     570    11980178 :       if (s%p != 0 && umael3(od, t, 1, 1) != 0)
     571             :       {
     572     9025940 :         u++;
     573     9025940 :         umael3(od, t, 1, 1) *= q;
     574     9025940 :         vecsmalltrunc_append( gmael(od, t, 2), p_to_index(q, indexlist) );
     575             :       } else
     576     2954238 :         umael3(od, t, 1, 1) = 0;
     577    11980178 :       s += 4;
     578             :     }
     579        4809 :     s = 4;
     580    11982775 :     for (t = tev; t <= lenv; t += p)
     581             :     {
     582    11977966 :       if (s%p != 0 && umael3(ev, t, 1, 1) != 0)
     583             :       {
     584     6771079 :         u++;
     585     6771079 :         umael3(ev, t, 1, 1) *= q;
     586     6771079 :         vecsmalltrunc_append( gmael(ev, t, 2), p_to_index(q, indexlist) );
     587             :       } else
     588     5206887 :         umael3(ev, t, 1, 1) = 0;
     589    11977966 :       s += 4;
     590             :     }
     591             :   }
     592             : 
     593             :   /* merging the two arrays */
     594          35 :   merge = cgetg(u+1, t_VEC);
     595          35 :   lgmerge = 0;
     596             :   while (1)
     597             :   {
     598             :     long o3, o4;
     599    14000035 :     if (u3 > lenv && u4 > lenv) break;
     600    14000000 :     o3 = -1;
     601    14000000 :     o4 = -1;
     602    14000000 :     if (u3 <= lenv)
     603             :     {
     604    14000000 :       o3 = umael3(od, u3, 1, 1);
     605    14000000 :       if (o3 != -4*u3+1) {u3++; continue;}
     606             :     }
     607     8351378 :     if (u4 <= lenv)
     608             :     {
     609     8351322 :       o4 = umael3(ev, u4, 1, 1);
     610     8351322 :       if (o4 != -4*u4) {u4++; continue;}
     611             :     }
     612     2373518 :     if (o3 == -1)
     613             :     {
     614           0 :       gel(merge, ++lgmerge) = gel(ev, u4++);
     615           0 :       continue;
     616             :     }
     617     2373518 :     if (o4 == -1)
     618             :     {
     619          56 :       gel(merge, ++lgmerge) = gel(od, u3++);
     620          56 :       continue;
     621             :     }
     622     2373462 :     if (o3 > o4)
     623     1351322 :       gel(merge, ++lgmerge) = gel(od, u3++);
     624             :     else
     625     1022140 :       gel(merge, ++lgmerge) = gel(ev, u4++);
     626    14000000 :   }
     627          35 :   setlg(merge, lgmerge);
     628             :   /* filling in bestinv [1][3] and poldegree [1][4] */
     629     2373518 :   for (i = 1; i < lgmerge; i++)
     630             :   {
     631     2373483 :     long D = umael3(merge, i, 1, 1);
     632     2373483 :     long h = umael3(merge, i, 1, 2);
     633     2373483 :     long modinv = umael3(merge, i, 1, 3) = disc_best_modinv(D);
     634     2373483 :     long p1 = 1, p2 = 1;
     635     2373483 :     if (modinv_degree(&p1, &p2, modinv) && (-D)%p1 == 0 && (-D)%p2 == 0)
     636      165424 :       umael3(merge, i, 1, 4) = h/2;
     637             :   }
     638          35 :   merge = gen_sort(merge, NULL, &sort_disclist);
     639          35 :   return gerepilecopy(av, merge);
     640             : }
     641             : 
     642             : /*  Input: a vector tune whose components are vectors of length 3
     643             :    Output: vector param of precomputations
     644             :              let x = [maxsqrt, maxdisc, maxpcdg] be the last component of tune then
     645             :              param[1][1] =   tune[#tune] = [maxsqrt, maxdisc, maxpcdg]
     646             :              param[1][2] =     primelist = [ Vecsmall with the list of primes ]
     647             :              param[1][3] =      disclist = vector whose elements are Dinfos, sorted by disclist
     648             :              param[2]    = primorial_vec
     649             :              param[3]    =          tune
     650             : */
     651             : static GEN
     652          35 : ecpp_param_set(GEN tune)
     653             : {
     654          35 :   pari_sp av = avma;
     655          35 :   long lgtune = lg(tune);
     656          35 :   GEN param1 = mkvec3(zero_zv(3), zerovec(1), zerovec(1));
     657          35 :   GEN param2 = gen_1;
     658          35 :   GEN param3 = tune;
     659          35 :   GEN param = mkvec3(param1, param2, param3);
     660          35 :   GEN x = gel(tune, lgtune-1);
     661          35 :   long  maxsqrt = typ(x) == t_VECSMALL ? uel(x, 1) : itos(gel(x, 1));
     662          35 :   ulong maxpcdg = typ(x) == t_VECSMALL ? uel(x, 2) : itos(gel(x, 2));
     663          35 :   ulong tdivexp = typ(x) == t_VECSMALL ? uel(x, 3) : itos(gel(x, 3));
     664             : 
     665          35 :   if (tune != NULL)
     666             :   {
     667          35 :     ecpp_param_set_maxsqrt(param, maxsqrt);
     668          35 :     ecpp_param_set_maxdisc(param, maxsqrt*maxsqrt);
     669          35 :     ecpp_param_set_maxpcdg(param, maxpcdg);
     670          35 :     ecpp_param_set_disclist(param);
     671          35 :     ecpp_param_set_tdivexp(param, tdivexp);
     672             :   }
     673          35 :   return gerepilecopy(av, param);
     674             : }
     675             : 
     676             : /* cert contains [N, t, s, a4, [x, y]] as documented in ??ecpp
     677             :    the following information can be obtained:
     678             :    D = squarefreepart(t^2-4N)
     679             :    m = (N+1-t)
     680             :    q = m/s
     681             :    a6 = y^3 - x^2 - a4*x
     682             :    J = use formula
     683             : */
     684             : INLINE GEN
     685        2100 : cert_get_N(GEN x)
     686        2100 : { return gel(x,1); }
     687             : 
     688             : INLINE GEN
     689        2023 : cert_get_t(GEN x)
     690        2023 : { return gel(x,2); }
     691             : 
     692             : INLINE GEN
     693          28 : cert_get_D(GEN x)
     694             : {
     695          28 :   GEN N = cert_get_N(x);
     696          28 :   GEN t = cert_get_t(x);
     697          28 :   GEN t2m4N = subii(sqri(t), shifti(N, 2));
     698          28 :   return coredisc(t2m4N);
     699             : }
     700             : 
     701             : INLINE GEN
     702        1036 : cert_get_m(GEN x)
     703             : {
     704        1036 :   GEN N = cert_get_N(x);
     705        1036 :   GEN t = cert_get_t(x);
     706        1036 :   return subii(addiu(N, 1), t);
     707             : }
     708             : 
     709             : INLINE GEN
     710        1036 : cert_get_s(GEN x)
     711             : {
     712        1036 :   return gel(x,3);
     713             : }
     714             : 
     715             : INLINE GEN
     716          77 : cert_get_q(GEN x)
     717             : {
     718          77 :   GEN m = cert_get_m(x);
     719          77 :   GEN s = cert_get_s(x);
     720          77 :   return diviiexact(m, s);
     721             : }
     722             : 
     723             : INLINE GEN
     724        1043 : cert_get_a4(GEN x)
     725             : {
     726        1043 :   return gel(x, 4);
     727             : }
     728             : 
     729             : INLINE GEN
     730        1029 : cert_get_P(GEN x)
     731             : {
     732        1029 :   return gel(x, 5);
     733             : }
     734             : 
     735             : INLINE GEN
     736          14 : cert_get_x(GEN x)
     737             : {
     738          14 :   return gel(cert_get_P(x), 1);
     739             : }
     740             : 
     741             : INLINE GEN
     742          56 : cert_get_a6(GEN z)
     743             : {
     744          56 :   GEN N = cert_get_N(z);
     745          56 :   GEN a = cert_get_a4(z);
     746          56 :   GEN P = cert_get_P(z);
     747             : 
     748          56 :   GEN x = gel(P, 1);
     749          56 :   GEN y = gel(P, 2);
     750          56 :   GEN yy = Fp_sqr(y, N);
     751          56 :   GEN xx = Fp_sqr(x, N);
     752          56 :   GEN xxpa = Fp_add(xx, a, N);
     753          56 :   GEN xxxpax = Fp_mul(x, xxpa, N);
     754          56 :   GEN yymxxxmax = Fp_sub(yy, xxxpax, N);
     755          56 :   return yymxxxmax;
     756             : }
     757             : 
     758             : INLINE GEN
     759          28 : cert_get_E(GEN x)
     760             : {
     761          28 :   GEN a = cert_get_a4(x);
     762          28 :   GEN b = cert_get_a6(x);
     763          28 :   return mkvec2(a, b);
     764             : }
     765             : 
     766             : INLINE GEN
     767           0 : cert_get_J(GEN x)
     768             : {
     769           0 :   GEN N = cert_get_N(x);
     770           0 :   GEN a = cert_get_a4(x);
     771           0 :   GEN b = cert_get_a6(x);
     772           0 :   return Fp_ellj(a, b, N);
     773             : }
     774             : 
     775             : /* "Twist factor"
     776             :    Does not cover J = 0, 1728 */
     777             : static GEN
     778           0 : FpE_get_lambda(GEN a, GEN b, GEN A, GEN B, GEN N)
     779             : {
     780           0 :   GEN aB = Fp_mul(a, B, N);
     781           0 :   GEN bA = Fp_mul(b, A, N);
     782           0 :   return Fp_div(aB, bA, N);
     783             : }
     784             : 
     785             : /*  Input: J, N
     786             :    Output: [a, b]
     787             :            a = 3*J*(1728-J)   mod N
     788             :            b = 2*J*(1728-J)^2 mod N
     789             : */
     790             : static GEN
     791         665 : Fp_ellfromj(GEN j, GEN N)
     792             : {
     793             :   GEN k, kk, jk, a, b;
     794         665 :   j = Fp_red(j, N);
     795         665 :   if (isintzero(j)) return mkvec2(gen_0, gen_1);
     796         518 :   if (equalii(Fp_red(stoi(1728), N), j)) return mkvec2(gen_1, gen_0);
     797         469 :   k = Fp_sub(stoi(1728), j, N);
     798         469 :   kk = Fp_sqr(k, N);
     799         469 :   jk = Fp_mul(j, k, N);
     800         469 :   a = Fp_mulu(jk, 3, N);
     801         469 :   b = Fp_mulu(Fp_mul(j, kk, N), 2, N);
     802         469 :   return mkvec2(a, b);
     803             : }
     804             : 
     805             : static GEN
     806           0 : cert_get_lambda(GEN x)
     807             : {
     808             :   GEN N, J, a, b, A, B, AB;
     809           0 :   J = cert_get_J(x);
     810           0 :   N = cert_get_N(x);
     811           0 :   a = cert_get_a4(x);
     812           0 :   b = cert_get_a6(x);
     813           0 :   AB = Fp_ellfromj(J, N);
     814           0 :   A = gel(AB, 1);
     815           0 :   B = gel(AB, 2);
     816           0 :   return FpE_get_lambda(a, b, A, B, N);
     817             : }
     818             : 
     819             : /* Solves for T such that if
     820             :    [A, B] = [3*J*(1728-J), 2*J*(1728-J)^2]
     821             :    and you let
     822             :    L = T^3 + A*T + B
     823             :    A = A*L^2
     824             :    B = B*L^3
     825             :    then
     826             :    x == TL
     827             :    y == L^2
     828             : */
     829             : static GEN
     830           0 : cert_get_T(GEN z)
     831             : {
     832           0 :   GEN N = cert_get_N(z);
     833           0 :   GEN P = cert_get_P(z);
     834           0 :   GEN x = gel(P, 1);
     835           0 :   GEN l = cert_get_lambda(z); /* l = 1/L */
     836           0 :   GEN T = Fp_mul(x, l, N);
     837           0 :   return T;
     838             : }
     839             : 
     840             : /* database for all invariants
     841             :    stolen from Hamish's code
     842             : */
     843             : static GEN
     844          28 : polmodular_db_init_allinv(void)
     845             : {
     846             :   GEN ret1;
     847          28 :   GEN ret2 = cgetg(39+1, t_VEC);
     848             :   enum { DEFAULT_MODPOL_DB_LEN = 32 };
     849             :   long i;
     850        1120 :   for (i = 1; i < lg(ret2); i++)
     851        1092 :     gel(ret2, i) = zerovec_block(DEFAULT_MODPOL_DB_LEN);
     852          28 :   ret1 = zerovec_block(DEFAULT_MODPOL_DB_LEN);
     853          28 :   return mkvec2(ret1, ret2);
     854             : }
     855             : 
     856             : /*  Input: a discriminant D and a database of previously computed modular polynomials,
     857             :    Output: polclass(D,disc_best_modinv(D))
     858             : */
     859             : static GEN
     860         427 : D_polclass(GEN D, GEN *db)
     861             : {
     862             :   GEN HD;
     863         427 :   long inv = disc_best_modinv(itos(D));
     864         427 :   GEN tmp_db = mkvec2(gel(*db, 1), gmael(*db, 2, inv));
     865         427 :   if (inv == 0) tmp_db = mkvec2(gel(*db, 1), gen_0);
     866         427 :   HD = polclass0( itos(D), inv, 0, &tmp_db);
     867         427 :   gel(*db, 1) = gel(tmp_db, 1);
     868         427 :   if (inv != 0) gmael(*db, 2, inv) = gel(tmp_db, 2);
     869         427 :   return HD;
     870             : }
     871             : 
     872             : /*  Input: N, Dinfo, a root rt mod N of polclass(D,disc_best_modinv(D))
     873             :    Output: a root J mod N of polclass(D)
     874             : */
     875             : INLINE GEN
     876         665 : NDinfor_find_J(GEN N, GEN Dinfo, GEN rt)
     877             : {
     878         665 :   long inv = Dinfo_get_bi(Dinfo);
     879         665 :   return Fp_modinv_to_j(rt, inv, N);
     880             : }
     881             : 
     882             : INLINE long
     883        1385 : NmqEP_check(GEN N, GEN q, GEN E, GEN P, GEN s)
     884             : {
     885        1385 :   GEN a = gel(E, 1);
     886             :   GEN mP, sP;
     887        1385 :   sP = FpJ_mul(P, s, a, N);
     888        1385 :   if (FpJ_is_inf(sP)) return 0;
     889        1385 :   mP = FpJ_mul(sP, q, a, N);
     890        1385 :   if (FpJ_is_inf(mP)) return 1;
     891         720 :   return 0;
     892             : }
     893             : 
     894             : /* This finds an elliptic curve E modulo N and a point P on E
     895             :      which corresponds to the ith element of the certificate.
     896             :    It uses the corresponding N, D, m, q, J obtained in previous steps.
     897             : 
     898             :    All computations are to be done modulo N unless stated otherwise.
     899             : */
     900             : 
     901             : /* g is a quadratic and cubic non-residue modulo N */
     902             : static GEN
     903         268 : j0_find_g(GEN N)
     904             : {
     905             :   while (1)
     906             :   {
     907         268 :     GEN g = randomi(N);
     908         268 :     if (kronecker(g, N) != -1) continue;
     909         171 :     if (isint1(Fp_pow(g, diviuexact(subiu(N, 1), 3), N))) continue;
     910         224 :     return g;
     911         156 :   }
     912             : }
     913             : 
     914             : static GEN
     915        1385 : random_FpJ(GEN A, GEN B, GEN N)
     916             : {
     917        1385 :   GEN P = random_FpE(A, B, N);
     918        1385 :   return FpE_to_FpJ(P);
     919             : }
     920             : 
     921             : static GEN
     922         665 : NDinfomqgJ_find_EP(GEN N, GEN Dinfo, GEN m, GEN q, GEN g, GEN J, GEN s)
     923             : {
     924             :   long i;
     925         665 :   long D = Dinfo_get_D(Dinfo);
     926             :   GEN gg;
     927         665 :   GEN E = Fp_ellfromj(J, N);
     928         665 :   GEN A = gel(E, 1);
     929         665 :   GEN B = gel(E, 2);
     930         665 :   GEN P = random_FpJ(A, B, N);
     931         665 :   if (NmqEP_check(N, q, E, P, s)) return mkvec2(E, P);
     932         357 :   switch (D_get_wD(D))
     933             :   {
     934             :     case 2:
     935         203 :       gg = Fp_sqr(g, N);
     936         203 :       A = Fp_mul(A, gg, N); /* Ag^2 */
     937         203 :       B = Fp_mul(Fp_mul(B, gg, N), g, N); /* Bg^3 */
     938         203 :       E = mkvec2(A, B);
     939         203 :       P = random_FpJ(A, B, N);
     940         203 :       if (NmqEP_check(N, q, E, P, s)) return mkvec2(E, P);
     941           0 :       else return NDinfomqgJ_find_EP(N, Dinfo, m, q, g, J, s);
     942             :     case 4:
     943          56 :       for (i = 1; i < 4; i++)
     944             :       {
     945          56 :         A = Fp_mul(A, g, N); /* Ag */
     946          56 :         E = mkvec2(A, B);
     947          56 :         P = random_FpJ(A, B, N);
     948          56 :         if (NmqEP_check(N, q, E, P, s)) return mkvec2(E, P);
     949             :       }
     950           0 :       return NDinfomqgJ_find_EP(N, Dinfo, m, q, g, J, s);
     951             :     case 6:
     952         126 :       B = Fp_mul(B, g, N); /* Bg */
     953         126 :       E = mkvec2(A, B);
     954         126 :       P = random_FpJ(A, B, N);
     955         126 :       if (NmqEP_check(N, q, E, P, s)) return mkvec2(E, P);
     956         112 :       g = j0_find_g(N);
     957         417 :       for (i = 1; i < 6; i++)
     958             :       {
     959         417 :         B = Fp_mul(B, g, N); /* Bg */
     960         417 :         if (i % 3 == 0) continue;
     961         335 :         E = mkvec2(A, B);
     962         335 :         P = random_FpJ(A, B, N);
     963         335 :         if (NmqEP_check(N, q, E, P, s)) return mkvec2(E, P);
     964             :       }
     965           0 :       return NDinfomqgJ_find_EP(N, Dinfo, m, q, g, J, s);
     966             :   }
     967           0 :   pari_err(e_BUG, "NDinfomqgJ_find_EP");
     968           0 :   return NULL;
     969             : }
     970             : 
     971             : /* This uses the unravelled [N, D, m, q] in step 1
     972             :    to find the appropriate j-invariants
     973             :    to be used in step 3.
     974             :    T2v is a timer used.
     975             :    Step 2 is divided into three substeps 2a, 2b, 2c. */
     976             : static GEN
     977          28 : ecpp_step2(GEN step1, GEN *X0)
     978             : {
     979             :   long j;
     980             :   pari_timer ti;
     981          28 :   GEN perm = gen_indexsort(step1, NULL, &sort_NDmq_by_D);
     982          28 :   GEN step2 = cgetg(lg(step1), t_VEC);
     983          28 :   GEN HD = NULL;
     984          28 :   GEN Dprev = gen_0;
     985          28 :   GEN db = polmodular_db_init_allinv();
     986             : 
     987         693 :   for (j = 1; j < lg(step2); j++)
     988             :   {
     989             : 
     990         665 :     long i = uel(perm, j), tt = 0;
     991         665 :     GEN NDinfomqg_i = gel(step1, i);
     992             : 
     993         665 :     GEN N = NDinfomqg_get_N(NDinfomqg_i);
     994         665 :     GEN Dinfo = NDinfomqg_get_Dinfo(NDinfomqg_i);
     995         665 :     GEN D = NDinfomqg_get_D(NDinfomqg_i);
     996         665 :     GEN m = NDinfomqg_get_m(NDinfomqg_i);
     997         665 :     GEN q = NDinfomqg_get_q(NDinfomqg_i);
     998         665 :     GEN g = NDinfomqg_get_g(NDinfomqg_i);
     999             :     GEN J, t, s, a4, P, EP, rt;
    1000             : 
    1001             : 
    1002             :     /* C1: Find the appropriate class polynomial modulo N. */
    1003         665 :     dbg_mode() timer_start(&ti);
    1004         665 :     if (!equalii(D, Dprev)) HD = D_polclass(D, &db);
    1005         665 :     dbg_mode() tt = timer_record(X0, "C1", &ti, 1);
    1006         665 :     dbg_mode() err_printf(ANSI_COLOR_BRIGHT_GREEN "\n[ %3d | %4ld bits]" ANSI_COLOR_RESET, i, expi(N));
    1007         665 :     dbg_mode() err_printf(ANSI_COLOR_GREEN "      D = %8Ps  poldeg = %4ld" ANSI_COLOR_RESET, D, degpol(HD));
    1008         665 :     dbg_mode() if (equalii(D, Dprev)) err_printf("  %8ld", tt);
    1009         665 :     dbg_mode() if (!equalii(D, Dprev)) err_printf(ANSI_COLOR_BRIGHT_WHITE "  %8ld" ANSI_COLOR_RESET, tt);
    1010             : 
    1011             :     /* C2: Find a root modulo N of the polynomial obtained in the previous step. */
    1012         665 :     dbg_mode() timer_start(&ti);
    1013         665 :     rt = FpX_oneroot_split(HD, N);
    1014         665 :     dbg_mode() tt = timer_record(X0, "C2", &ti, 1);
    1015         665 :     dbg_mode() err_printf("  %8ld", tt);
    1016             : 
    1017             :     /* C3: Convert the root obtained from the previous step into the appropriate j-invariant. */
    1018         665 :     dbg_mode() timer_start(&ti);
    1019         665 :     J = NDinfor_find_J(N, Dinfo, rt);
    1020         665 :     dbg_mode() tt = timer_record(X0, "C3", &ti, 1);
    1021         665 :     dbg_mode() err_printf("  %8ld", tt);
    1022             : 
    1023             :     /* D1: Find an elliptic curve E with a point P satisfying the theorem. */
    1024         665 :     dbg_mode() timer_start(&ti);
    1025         665 :     s = diviiexact(m, q);
    1026         665 :     EP = NDinfomqgJ_find_EP(N, Dinfo, m, q, g, J, s);
    1027         665 :     dbg_mode() tt = timer_record(X0, "D1", &ti, 1);
    1028         665 :     dbg_mode() err_printf("  %8ld", tt);
    1029             : 
    1030             :     /* D2: Compute for t and s */
    1031         665 :     dbg_mode() timer_start(&ti);
    1032         665 :     t = subii(addiu(N, 1), m); /* t = N+1-m */
    1033         665 :     a4 = gmael(EP, 1, 1);
    1034         665 :     P = FpJ_to_FpE(gel(EP, 2), N);
    1035         665 :     dbg_mode() tt = timer_record(X0, "D2", &ti, 1);
    1036         665 :     dbg_mode() err_printf("  %8ld", tt);
    1037             : 
    1038         665 :     gel(step2, i) = mkvec5(N, t, s, a4, P);
    1039             : 
    1040         665 :     Dprev = D;
    1041             :   }
    1042             : 
    1043          28 :   return step2;
    1044             : }
    1045             : /* end of functions for step 2 */
    1046             : 
    1047             : /* start of functions for step 1 */
    1048             : 
    1049             : /* start of macros for step 1 */
    1050             : 
    1051             : /* earlyabort_pcdg is used when the downrun is told to NOT persevere
    1052             :    it returns TRUE (i.e. do an early abort) if the polynomial degree
    1053             :      of the discriminant in Dinfo[i] is larger than the param maxpcdg.
    1054             : */
    1055             : INLINE long
    1056      158039 : earlyabort_pcdg(GEN param, ulong maxpcdg, long i)
    1057             : {
    1058      158039 :   GEN x = ecpp_param_get_disclist(param);
    1059      158039 :   GEN Dinfo = gel(x, i);
    1060      158039 :   ulong pd = Dinfo_get_pd(Dinfo);
    1061      158039 :   return (pd > maxpcdg);
    1062             : }
    1063             : 
    1064             : /*  Input: vector whose components are [D, m]
    1065             :    Output: vector whose components are m
    1066             : */
    1067             : static GEN
    1068        1351 : Dmvec_to_mvec(GEN Dmvec)
    1069             : {
    1070        1351 :   long lgmvec = lg(Dmvec);
    1071        1351 :   GEN mvec = cgetg(lgmvec, t_VEC);
    1072             :   long i;
    1073        1351 :   for (i = 1; i < lgmvec; i++) gel(mvec, i) = gmael(Dmvec, i, 2);
    1074        1351 :   return mvec;
    1075             : }
    1076             : 
    1077             : /*  Input: vector v whose components are [D, m], vector w whose components are q
    1078             :    Output: vector whose components are [D, m, q]
    1079             : */
    1080             : static GEN
    1081        1351 : Dmvec_qvec_to_Dmqvec(GEN Dmvec, GEN qvec)
    1082             : {
    1083        1351 :   long lgDmqvec = lg(Dmvec);
    1084        1351 :   GEN Dmqvec = cgetg(lgDmqvec, t_VEC);
    1085             :   long i;
    1086       23205 :   for (i = 1; i < lgDmqvec; i++)
    1087             :   {
    1088       21854 :     GEN D = gmael(Dmvec, i, 1);
    1089       21854 :     GEN m = gmael(Dmvec, i, 2);
    1090       21854 :     GEN q = gel(qvec, i);
    1091       21854 :     gel(Dmqvec, i) = mkvec3(D, m, q);
    1092             :   }
    1093        1351 :   return Dmqvec;
    1094             : }
    1095             : 
    1096             : /*  Input: vector whose components are [D, m, q]
    1097             :    Output: vector whose components are q
    1098             : */
    1099             : static GEN
    1100        1351 : Dmqvec_to_qvec(GEN Dmqvec)
    1101             : {
    1102        1351 :   long lgqvec = lg(Dmqvec);
    1103        1351 :   GEN qvec = cgetg(lgqvec, t_VEC);
    1104             :   long i;
    1105        1351 :   for (i = 1; i < lgqvec; i++) gel(qvec, i) = gmael(Dmqvec, i, 3);
    1106        1351 :   return qvec;
    1107             : }
    1108             : 
    1109             : /* This initializes sqrtlist as a vector [A, B]
    1110             :      where A is a t_VEC of gen_0's
    1111             :        and B is a t_VECSMALL of 0's,
    1112             :    both of which are as long as primelist. */
    1113             : INLINE void
    1114         686 : primelist_sqrtlist_init(GEN primelist, GEN *sqrtlist)
    1115             : {
    1116         686 :   long l = lg(primelist)-1;
    1117         686 :   *sqrtlist = mkvec2(zerovec(l), zero_zv(l));
    1118         686 : }
    1119             : 
    1120             : /* This returns the square root modulo N
    1121             :      of the ith entry of the primelist.
    1122             :    If this square root is already available on sqrtlist,
    1123             :      then simply return it.
    1124             :    Otherwise, compute it, save it and return it.
    1125             :    y is a quadratic non-residue that is needed
    1126             :      somehow in the algorithm for taking square roots modulo N.
    1127             : */
    1128             : static GEN
    1129       64988 : p_find_primesqrt(GEN N, GEN* X0, GEN primelist, GEN sqrtlist, long i, GEN g)
    1130             : {
    1131             :   pari_timer ti;
    1132       64988 :   long p = uel(primelist, i);
    1133       64988 :   if (isintzero(gmael(sqrtlist,1,i)))
    1134             :   {
    1135       16009 :     dbg_mode() err_printf(ANSI_COLOR_MAGENTA "S" ANSI_COLOR_RESET);
    1136             :     /* A4: Get the square root of a prime factor of D. */
    1137       16009 :     dbg_mode() timer_start(&ti);
    1138       16009 :     gmael(sqrtlist, 1, i) = Fp_sqrt_i(stoi(p), g, N); /* NULL if invalid. */
    1139       16009 :     dbg_mode() timer_record(X0, "A4", &ti, 1);
    1140             :   }
    1141       64988 :   return gmael(sqrtlist, 1, i);
    1142             : }
    1143             : 
    1144             : /* This finds the legit square root of D modulo N where D is the discriminant in Dinfo.
    1145             :    This algorithm finds the square root modulo N of each prime p dividing D.
    1146             :    It then assembles the actual square root of D by multiplying the prime square roots.
    1147             : */
    1148             : static GEN
    1149       34027 : D_find_discsqrt(GEN N, GEN param, GEN *X0, GEN Dinfo, GEN sqrtlist, GEN g)
    1150             : {
    1151       34027 :   GEN discsqrt = gen_1;
    1152       34027 :   GEN Dfac = Dinfo_get_Dfac(Dinfo);
    1153       34027 :   long lgDfac = lg(Dfac);
    1154             :   long i;
    1155       34027 :   GEN primelist = ecpp_param_get_primelist(param);
    1156       99015 :   for (i = 1; i < lgDfac; i++)
    1157             :   {
    1158       64988 :     GEN sqrtfi = p_find_primesqrt(N, X0, primelist, sqrtlist, uel(Dfac, i), g);
    1159       64988 :     if (sqrtfi == NULL) return NULL; /* Only happens when N is composite. */
    1160       64988 :     discsqrt = Fp_mul(discsqrt, sqrtfi, N);
    1161             :   }
    1162       34027 :   return discsqrt;
    1163             : }
    1164             : 
    1165             : /* Given a solution U, V to U^2 + |D|V^2 = 4N, this find all the possible
    1166             :      cardinalities of the elliptic curves over the finite field F_N
    1167             :      whose endomorphism ring is isomorphic to the maximal order
    1168             :      in the imaginary quadratic number field K = Q(sqrt(D)). ???
    1169             : */
    1170             : static GEN
    1171        9884 : NUV_find_mvec(GEN N, GEN U, GEN V, long wD)
    1172             : {
    1173        9884 :   GEN Nplus1 = addiu(N, 1);
    1174             :   GEN m;
    1175        9884 :   GEN u = U;
    1176        9884 :   GEN mvec = cgetg(wD + 1, t_VEC);
    1177             :   long i;
    1178       31738 :   for (i = 0; i < wD; i++)
    1179             :   {
    1180       21854 :     if (i%2 == 0)
    1181             :     {
    1182       10927 :       if (wD == 4 && i==2) u = shifti(V, 1);
    1183       10927 :       if (wD == 6 && i==2) u = shifti(submuliu(U, V, 3), -1);
    1184       10927 :       if (wD == 6 && i==4) u = shifti(addmuliu(U, V, 3), -1);
    1185       10927 :       m = addii(Nplus1, u);
    1186             :     } else
    1187       10927 :       m = subii(Nplus1, u);
    1188       21854 :     gel(mvec, i + 1) = m;
    1189             :   }
    1190        9884 :   return mvec;
    1191             : }
    1192             : 
    1193             : /* This populates Dmbatch with Dmvec's -- vectors whose components are of the form [D, m],
    1194             :      where m is a cardinalities of an elliptic curve with endomorphism ring isomorphic to
    1195             :      the maximal order of the imaginary quadratic number field K = Q(sqrt(D)).
    1196             :    It returns 0 if:
    1197             :      * D (of Dinfo) is not a quadratic residue mod N
    1198             :      * Any of the p* dividing D is not a quadratic residue mod N
    1199             :      * Cornacchia cannot find a solution U^2 + |D|V^2 = 4N.
    1200             :    Otherwise, it returns the number of cardinalities added to Dmbatch.
    1201             :    Moreover, if N is determined to be composite, it sets failflag to 1.
    1202             :    Finally, sqrtlist and y are used to help compute the square root modulo N of D+.
    1203             : */
    1204             : static long
    1205      319998 : D_collectcards(GEN N, GEN param, GEN* X0, GEN Dinfo, GEN sqrtlist, GEN g, GEN Dmbatch, long* failflag)
    1206             : {
    1207             : 
    1208             :   long i, j;
    1209             :   GEN U, V;
    1210             :   long corn_succ;
    1211             :   long wD;
    1212             :   GEN Dfac;
    1213             :   long lgDfac;
    1214             :   GEN sqrtofDmodN;
    1215      319998 :   GEN primelist = ecpp_param_get_primelist(param);
    1216             :   GEN mvec;
    1217             : 
    1218             :   /* Unpacking info on the discriminant. */
    1219      319998 :   long D = gel(Dinfo, 1)[1];
    1220             : 
    1221             :   /* timers / counters */
    1222             :   pari_timer ti;
    1223             :   long kronDN;
    1224             : 
    1225             :   /* A1: Check (D|N) = 1. */
    1226      319998 :   dbg_mode() timer_start(&ti);
    1227      319998 :   kronDN = krosi(D, N);
    1228      319998 :   dbg_mode() timer_record(X0, "A1", &ti, 1);
    1229      319998 :   switch(kronDN) {
    1230           0 :      case 0: *failflag = 1;
    1231      159390 :     case -1: return 0;
    1232             :   }
    1233             : 
    1234             :   /* A2: Check (p*|N) = 1 for all p.
    1235             :      This is equivalent to checking (N|p) = 1. */
    1236      160608 :   dbg_mode() timer_start(&ti);
    1237      160608 :   Dfac = Dinfo_get_Dfac(Dinfo);
    1238      160608 :   lgDfac = lg(Dfac);
    1239      251538 :   for (i = 1; i < lgDfac; i++)
    1240             :   {
    1241      217511 :     long p = index_to_p(uel(Dfac, i), primelist);
    1242      217511 :     if (krosi(p, N) != 1) return 0;
    1243             :   }
    1244       34027 :   dbg_mode() timer_record(X0, "A2", &ti, 1);
    1245             : 
    1246             :   /* A3: Get square root of D mod N. */
    1247             :   /* If sqrtofDmodN is NULL, then N is composite. */
    1248       34027 :   dbg_mode() timer_start(&ti);
    1249       34027 :   sqrtofDmodN = D_find_discsqrt(N, param, X0, Dinfo, sqrtlist, g);
    1250       34027 :   if (sqrtofDmodN == NULL) pari_err_BUG("D_find_discsqrt");
    1251       34027 :   dbg_mode() if (!equalii(Fp_sqr(sqrtofDmodN, N), addis(N, D)) /* D mod N, D < 0*/ ) pari_err_BUG("D_find_discsqrt");
    1252       34027 :   dbg_mode() timer_record(X0, "A3", &ti, 1);
    1253             : 
    1254             :   /* A5: Use the square root to use cornacchia to find the solution to U^2 + |D|V^2 = 4N. */
    1255       34027 :   dbg_mode() timer_start(&ti);
    1256       34027 :   corn_succ = cornacchia2_sqrt( absi(stoi(D)), N, sqrtofDmodN, &U, &V);
    1257       34027 :   dbg_mode() timer_record(X0, "A5", &ti, 1);
    1258       34027 :   if (!corn_succ) {
    1259       24143 :     dbg_mode() err_printf(ANSI_COLOR_YELLOW "c" ANSI_COLOR_RESET);
    1260       24143 :     return 0;
    1261             :   }
    1262             : 
    1263             :   /* We're sticking with this D. */
    1264        9884 :   dbg_mode() err_printf(ANSI_COLOR_BRIGHT_YELLOW "D" ANSI_COLOR_RESET);
    1265             : 
    1266        9884 :   dbg_mode() timer_start(&ti);
    1267             :   /* A6: Collect the w(D) possible cardinalities of elliptic curves over F_N whose discriminant is D. */
    1268        9884 :   wD = D_get_wD(D);
    1269        9884 :   mvec = NUV_find_mvec(N, U, V, wD);
    1270        9884 :   for (j = 1; j < lg(mvec); j++) vectrunc_append(Dmbatch, mkvec2(Dinfo, gel(mvec, j)) );
    1271        9884 :   dbg_mode() timer_record(X0, "A6", &ti, 1);
    1272             : 
    1273        9884 :   return wD;
    1274             : }
    1275             : 
    1276             : /* Compute (S(16N, 2) + S(4096N, 4) + 4)\4 + 1
    1277             :      where S is the nth root rounded down.
    1278             :    This is at most one more than (N^1/4 + 1)^2.
    1279             : */
    1280             : static GEN
    1281        1351 : ecpp_qlo(GEN N)
    1282             : {
    1283        1351 :   GEN a = sqrtnint(shifti(N, 4), 2);
    1284        1351 :   GEN b = sqrtnint(shifti(N, 12), 4);
    1285        1351 :   GEN c = shifti(gen_1, 2);
    1286        1351 :   GEN d = addii(addii(a, b), c);
    1287        1351 :   GEN e = shifti(d, -2);
    1288        1351 :   GEN f = addiu(e, 1);
    1289        1351 :   return f;
    1290             : }
    1291             : 
    1292             : /* E is &qlo, use for zv_binsearch */
    1293             : static long
    1294        6923 : lessthan_qlo(void* E, GEN q)
    1295             : {
    1296        6923 :   GEN qlo = *((GEN*)E);
    1297        6923 :   return (cmpii(q, qlo) < 0);
    1298             : }
    1299             : 
    1300             : /* E is &goal, use for zv_binsearch */
    1301             : static long
    1302        6097 : gained_bits(void* E, GEN q)
    1303             : {
    1304        6097 :   long goal = *((long*)E);
    1305        6097 :   return (expi(q) <= goal);
    1306             : }
    1307             : 
    1308             : /*  Input: Dmqvec
    1309             :    Output: Dmqvec such that q satisfies
    1310             :      (N^1/4 + 1)^2 < q < N/2
    1311             : */
    1312             : static GEN
    1313        1351 : Dmqvec_slice_Dmqvec(GEN N, GEN Dmqvec)
    1314             : {
    1315        1351 :   pari_sp av = avma;
    1316             :   GEN qlo;
    1317             :   GEN qvec;
    1318             :   long lo_ind, bitgain, hi_ind, goal;
    1319             : 
    1320             :   /* Dmqvec is sorted according to q */
    1321        1351 :   Dmqvec = gen_sort(Dmqvec, NULL, &sort_Dmq_by_q);
    1322        1351 :   qvec = Dmqvec_to_qvec(Dmqvec);
    1323             : 
    1324        1351 :   qlo = ecpp_qlo(N);
    1325        1351 :   lo_ind = zv_binsearch0(&qlo, &lessthan_qlo, qvec); lo_ind++;
    1326        1351 :   if (lo_ind >= lg(qvec)) { avma = av; return NULL; }
    1327             : 
    1328        1351 :   bitgain = 1;
    1329        1351 :   goal = expi(N)-bitgain;
    1330        1351 :   hi_ind = zv_binsearch0(&goal, &gained_bits, qvec);
    1331        1351 :   if (hi_ind == 0) { avma = av; return NULL; }
    1332             : 
    1333        1351 :   Dmqvec = vecslice(Dmqvec, lo_ind, hi_ind);
    1334        1351 :   gerepileall(av, 1, &Dmqvec);
    1335        1351 :   return Dmqvec;
    1336             : }
    1337             : 
    1338             : /* Given a vector mvec of mi's,
    1339             :    This simultaneously removes all prime factors of each mi less then BOUND_PRIMORIAL
    1340             :    This implements Franke 2004: Proving the Primality of Very Large Numbers with fastECPP. */
    1341             : static GEN
    1342        1351 : mvec_batchfactor_qvec(GEN mvec, GEN primorial)
    1343             : {
    1344        1351 :   pari_sp av = avma;
    1345             :   long i;
    1346             :   /* Obtain the product tree of mvec. */
    1347        1351 :   GEN leaf = Z_ZV_mod_tree(primorial, mvec, ZV_producttree(mvec));
    1348             : 
    1349             :   /* Go through each leaf and remove small factors. */
    1350       23205 :   for (i = 1; i < lg(mvec); i++)
    1351             :   {
    1352       21854 :     GEN m = gel(mvec, i);
    1353      110117 :     while (!isint1(gel(leaf, i)) )
    1354             :     {
    1355       66409 :       gel(leaf, i) = gcdii( m, gel(leaf,i) );
    1356       66409 :       m = diviiexact( m, gel(leaf,i) );
    1357             :     }
    1358       21854 :     gel(mvec, i) = m;
    1359             :   }
    1360        1351 :   gerepileall(av, 1, &mvec);
    1361        1351 :   return mvec;
    1362             : }
    1363             : 
    1364             : /*  Input: Dmvec
    1365             :    Output: Dmqvec
    1366             :    Uses mvec_batchfactor_qvec to produce the output.
    1367             : */
    1368             : static GEN
    1369        1351 : Dmvec_batchfactor_Dmqvec(GEN Dmvec, GEN primorial)
    1370             : {
    1371        1351 :   pari_sp av = avma;
    1372        1351 :   GEN mvec = Dmvec_to_mvec(Dmvec);
    1373        1351 :   GEN qvec = mvec_batchfactor_qvec(mvec, primorial);
    1374        1351 :   GEN Dmqvec = Dmvec_qvec_to_Dmqvec(Dmvec, qvec);
    1375        1351 :   gerepileall(av, 1, &Dmqvec);
    1376        1351 :   return Dmqvec;
    1377             : }
    1378             : 
    1379             : /* tunevec = [maxsqrt, maxpcdg, tdivexp]
    1380             :      and is a shorthand for specifying the parameters of ecpp
    1381             : */
    1382             : static GEN
    1383        2723 : tunevec(GEN N, GEN param)
    1384             : {
    1385        2723 :   long e = expi(N);
    1386        2723 :   GEN T = gel(param,3);
    1387        2723 :   if      (e <= 1500) return gel(T,1);
    1388          91 :   else if (e <= 2500) return gel(T,2);
    1389           0 :   else if (e <= 3500) return gel(T,3);
    1390           0 :   return gel(T,4);
    1391             : }
    1392             : static long
    1393        2723 : tunevec_tdivbd(GEN N, GEN param)
    1394             : {
    1395        2723 :   return uel(tunevec(N, param), 3);
    1396             : }
    1397             : static long
    1398         686 : tunevec_batchsize(GEN N, GEN param)
    1399             : {
    1400         686 :   return (28-tunevec_tdivbd(N, param) < 0 ? 0 : 28-tunevec_tdivbd(N, param));
    1401             : }
    1402             : 
    1403             : /*  Input: Dmbatch (from the downrun)
    1404             :    Output: Dmqvec
    1405             : */
    1406             : static GEN
    1407        1351 : Dmbatch_factor_Dmqvec(GEN N, GEN* X0, GEN Dmbatch, GEN param)
    1408             : {
    1409             :   pari_timer ti;
    1410        1351 :   GEN curr_primorial = ecpp_param_get_primorial(param, tunevec_tdivbd(N, param) - 21);
    1411             :   GEN Dmqvec;
    1412             : 
    1413             :   /* B1: Factor by batch. */
    1414        1351 :   dbg_mode() timer_start(&ti);
    1415        1351 :   Dmqvec = Dmvec_batchfactor_Dmqvec(Dmbatch, curr_primorial);
    1416        1351 :   dbg_mode() timer_record(X0, "B1", &ti, 1);
    1417             : 
    1418             :   /* B2: For each batch,
    1419             :          remove cardinalities lower than (N^(1/4)+1)^2
    1420             :          and cardinalities in which we didn't win enough bits. */
    1421        1351 :   dbg_mode() timer_start(&ti);
    1422        1351 :   Dmqvec = Dmqvec_slice_Dmqvec(N, Dmqvec);
    1423        1351 :   dbg_mode() timer_record(X0, "B2", &ti, lg(Dmqvec)-1);
    1424             : 
    1425             :   /* If nothing is left after B2, return NULL */
    1426        1351 :   if (Dmqvec == NULL) return NULL;
    1427             : 
    1428        1351 :   return Dmqvec;
    1429             : }
    1430             : 
    1431             : /* Dmq macros */
    1432             : INLINE GEN
    1433       15778 : Dmq_get_Dinfo(GEN Dmq)
    1434       15778 : { return gel(Dmq, 1); }
    1435             : 
    1436             : INLINE GEN
    1437       15778 : Dmq_get_m(GEN Dmq)
    1438       15778 : { return gel(Dmq, 2); }
    1439             : 
    1440             : INLINE GEN
    1441       31556 : Dmq_get_q(GEN Dmq)
    1442       31556 : { return gel(Dmq, 3); }
    1443             : 
    1444             : INLINE long
    1445           0 : Dmq_get_cnum(GEN Dmq)
    1446           0 : { return gmael(Dmq, 1, 1)[2]; }
    1447             : 
    1448             : /*  Input: Dmq (where q does not have small prime factors)
    1449             :    Output: whether q is pseudoprime or not
    1450             : */
    1451             : static long
    1452       15778 : Dmq_isgoodq(GEN Dmq, GEN* X0)
    1453             : {
    1454             :   pari_timer ti;
    1455             :   long is_pseudoprime;
    1456       15778 :   GEN q = Dmq_get_q(Dmq);
    1457             : 
    1458             :   /* B3: Check pseudoprimality of each q on the list. */
    1459       15778 :   dbg_mode() timer_start(&ti);
    1460       15778 :   is_pseudoprime = BPSW_psp_nosmalldiv(q);
    1461       15778 :   dbg_mode() timer_record(X0, "B3", &ti, 1);
    1462       15778 :   return is_pseudoprime; /* did not find for this m */
    1463             : }
    1464             : 
    1465             : /*  Input: N
    1466             :    Output: [ NDmqg, therest ]
    1467             :      NDmqg is a vector of the form [N, D, m, q, g].
    1468             :      therest is a vector of the form the same as the output OR [N].
    1469             :    This is the downrun.
    1470             :    It tries to find [N, D, m, q, g].
    1471             :    If N is small, terminate.
    1472             :    It finds a quadratic non-residue g.
    1473             :    It starts with finding the square root of D modulo N.
    1474             :    It then uses this square root for cornacchia's algorithm.
    1475             :    If there is a solution to U^2 + |D|V^2 = 4N, use it to find candidates for m.
    1476             :    It continues until you get a batch of size at least t = log2(N)/2^log2bsize
    1477             :      or if there's no more discriminants left on the list.
    1478             :    It factors the m's of the batch and produces the q's.
    1479             :    If one of the q's are pseudoprime, then call this function again with N = q.
    1480             : */
    1481             : static GEN
    1482         714 : N_downrun_NDinfomq(GEN N, GEN param, GEN *X0, long *depth, long persevere)
    1483             : {
    1484         714 :   pari_sp ave = avma;
    1485             :   pari_timer T;
    1486             :   long i, j;
    1487         714 :   long expiN = expi(N);
    1488         714 :   long persevere_next = 0;
    1489             :   long lgdisclist, t, numcard;
    1490             :   ulong maxpcdg;
    1491             :   GEN primelist, disclist, sqrtlist, g, Dmbatch;
    1492         714 :   long FAIL = 0;
    1493             : 
    1494         714 :   dbg_mode() timer_start(&T);
    1495             : 
    1496             :   /* Unpack trustbits. */
    1497         714 :   if (expiN < 64) return gerepilecopy(ave, mkvec(N));
    1498             : 
    1499             :   /* This means we're going down the tree. */
    1500         686 :   *depth += 1;
    1501             : 
    1502             :   /* Unpack maxpcdg. */
    1503         686 :   maxpcdg = ecpp_param_get_maxpcdg(param);
    1504             : 
    1505             :   /* Unpack primelist, disclist. */
    1506         686 :   primelist = ecpp_param_get_primelist(param);
    1507         686 :   disclist = ecpp_param_get_disclist(param);
    1508         686 :   lgdisclist = lg(disclist);
    1509             : 
    1510             :   /* Initialize sqrtlist for this N. */
    1511         686 :   primelist_sqrtlist_init(primelist, &sqrtlist);
    1512             : 
    1513             :   /* Precomputation for batch size t. */
    1514             :   /* Tuning! */
    1515         686 :   t = expiN >> tunevec_batchsize(N, param);
    1516         686 :   if (t < 1) t = 1;
    1517             : 
    1518             :   /* Precomputation for taking square roots.
    1519             :        g will be needed for Fp_sqrt_i
    1520             :   */
    1521         686 :   g = Fp_2gener(N);
    1522         686 :   if (g == NULL) return gen_0; /* Composite if this happens. */
    1523             : 
    1524             :   /* Print the start of this iteration. */
    1525         686 :   dbg_mode() if (!persevere) err_printf(ANSI_COLOR_BRIGHT_CYAN "\n[ %3d | %4ld bits] " ANSI_COLOR_RESET, *depth, expiN, persevere);
    1526         686 :   dbg_mode() if (persevere) err_printf(ANSI_COLOR_BRIGHT_CYAN "\n< %3d | %4ld bits> " ANSI_COLOR_RESET, *depth, expiN, persevere);
    1527             : 
    1528             :   /* Initialize Dmbatch
    1529             :        It will be populated with candidate cardinalities on Phase I (until its length reaches at least t).
    1530             :        Its elements will be factored on Phase II.
    1531             :        We allocate a vectrunc_init of t+7.
    1532             :   */
    1533         686 :   Dmbatch = vectrunc_init(t+7);
    1534             : 
    1535             :   /* Number of cardinalities collected so far.
    1536             :      Should always be equal to lg(Dmbatch)-1. */
    1537         686 :   numcard = 0;
    1538             : 
    1539             :   /* i determines which discriminant we are considering. */
    1540         686 :   i = 1;
    1541             : 
    1542        2058 :   while (!FAIL)
    1543             :   {
    1544             :     pari_timer F;
    1545        1358 :     long last_i = i, failflag;
    1546             :     GEN Dmqlist;
    1547             :     long lgDmqlist;
    1548             : 
    1549             :     /* Dmbatch begins "empty", but we keep the allocated memory. */
    1550        1358 :     setlg(Dmbatch, 1);
    1551        1358 :     numcard = 0;
    1552             : 
    1553             :     /* PHASE I:
    1554             :        Go through the D's and search for canidate m's.
    1555             :        We fill up Dmbatch until there are at least t elements.
    1556             :     */
    1557             : 
    1558        1358 :     failflag = 0;
    1559      321377 :     while (i < lgdisclist )
    1560             :     {
    1561             :       GEN Dinfo;
    1562      320012 :       if (!persevere && earlyabort_pcdg(param, maxpcdg, i)) { FAIL = 1; break; }
    1563      319998 :       Dinfo = gel(disclist, i);
    1564      319998 :       numcard += D_collectcards(N, param, X0, Dinfo, sqrtlist, g, Dmbatch, &failflag);
    1565      320663 :       if (failflag) return gen_0;
    1566      319998 :       last_i = i++;
    1567      319998 :       if (numcard >= t) break;
    1568             :     }
    1569             : 
    1570             :     /* If we have exhausted disclist and there are no cardinalities to be factored. */
    1571        1365 :     if (FAIL && numcard <= 0) break;
    1572        1351 :     if (i >= lgdisclist && numcard <= 0) break;
    1573             : 
    1574             :     /* PHASE II:
    1575             :        Find the corresponding q's for the m's found.
    1576             :        We use Dmbatch.
    1577             :     */
    1578             : 
    1579             :     /* Find coresponding q of the candidate m's. */
    1580        1351 :     dbg_mode() timer_start(&F);
    1581        1351 :     Dmqlist = Dmbatch_factor_Dmqvec(N, X0, Dmbatch, param);
    1582             : 
    1583             :     /* If none left, move to the next discriminant. */
    1584        1351 :     if (Dmqlist == NULL) continue;
    1585             : 
    1586        1351 :     lgDmqlist = lg(Dmqlist);
    1587             : 
    1588             :     /* If we are cheating by adding class numbers, sort by class number instead. */
    1589        1351 :     if (Dinfo_get_pd(gel(disclist, last_i)) > maxpcdg)
    1590          28 :       Dmqlist = gen_sort(Dmqlist, NULL, &sort_Dmq_by_cnum);
    1591             : 
    1592             :     /* Check pseudoprimality before going down. */
    1593       16464 :     for (j = 1; j < lgDmqlist; j++)
    1594             :     {
    1595             :       GEN NDinfomq;
    1596       15778 :       GEN Dmq = gel(Dmqlist, j);
    1597       15778 :       GEN Dinfo = Dmq_get_Dinfo(Dmq);
    1598       15778 :       GEN m = Dmq_get_m(Dmq);
    1599       15778 :       GEN q = Dmq_get_q(Dmq);
    1600             :       GEN ret;
    1601             :       long a;
    1602       15778 :       if (expiN - expi(q) < 1)
    1603             :       {
    1604           0 :         dbg_mode() err_printf(ANSI_COLOR_BRIGHT_RED "  x" ANSI_COLOR_RESET);
    1605           0 :         continue;
    1606             :       }
    1607       15778 :       dbg_mode() err_printf(ANSI_COLOR_WHITE "." ANSI_COLOR_RESET);
    1608       15778 :       a = Dmq_isgoodq(Dmq, X0);
    1609       15778 :       if (!a) continue;
    1610             : 
    1611         679 :       dbg_mode() err_printf(ANSI_COLOR_BRIGHT_BLUE "  %ld" ANSI_COLOR_RESET, Dmq_get_cnum(Dmq));
    1612         679 :       dbg_mode() err_printf(ANSI_COLOR_BRIGHT_RED "\n       %5ld bits " ANSI_COLOR_RESET, expi(q)-expiN);
    1613             : 
    1614             :       /* Cardinality is pseudoprime. Call the next downrun! */
    1615         679 :       ret = N_downrun_NDinfomq(q, param, X0, depth, persevere_next);
    1616             : 
    1617             :       /* That downrun failed. */
    1618         679 :       if (ret == NULL) {
    1619          14 :         dbg_mode() if (!persevere) err_printf(ANSI_COLOR_CYAN "\n[ %3d | %4ld bits] " ANSI_COLOR_RESET, *depth, expiN, persevere);
    1620          14 :         dbg_mode() if (persevere) err_printf(ANSI_COLOR_CYAN "\n< %3d | %4ld bits> " ANSI_COLOR_RESET, *depth, expiN, persevere);
    1621          14 :         continue;
    1622             :       }
    1623             : 
    1624             :       /* That downrun succeeded. */
    1625         665 :       NDinfomq = mkvec5(N, Dinfo, m, q, g);
    1626         665 :       return gerepilecopy(ave, mkvec2(NDinfomq, ret));
    1627             :     }
    1628             : 
    1629             :     /* We have exhausted all the discriminants. */
    1630         686 :     if (i >= lgdisclist) FAIL = 1;
    1631             : 
    1632         686 :     if (Dinfo_get_pd(gel(disclist, last_i)) > maxpcdg)
    1633             :     {
    1634          28 :       dbg_mode() err_printf(ANSI_COLOR_BRIGHT_RED "  !" ANSI_COLOR_RESET);
    1635          28 :       persevere_next = 1;
    1636             :     }
    1637             : 
    1638             :   }
    1639             : 
    1640             :   /* FAILED: Out of discriminants. */
    1641          21 :   if (X0) umael(*X0, 3, 1)++; /* FAILS++ */
    1642          21 :   (*depth)--;
    1643          21 :   dbg_mode() err_printf(ANSI_COLOR_BRIGHT_RED "  X" ANSI_COLOR_RESET);
    1644          21 :   return NULL;
    1645             : }
    1646             : 
    1647             : /*  Input: the output of the (recursive) downrun function
    1648             :    Output: a vector whose components are [N, D, m, q, g]
    1649             : */
    1650             : static GEN
    1651          28 : ecpp_flattencert(GEN notflat, long depth)
    1652             : {
    1653          28 :   long i, lgret = depth+1;
    1654          28 :   GEN ret = cgetg(lgret, t_VEC);
    1655          28 :   GEN x = notflat;
    1656          28 :   if (typ(x) != t_VEC) return gen_0;
    1657          28 :   for (i = 1; i < lgret; i++) { gel(ret, i) = gel(x, 1); x = gel(x, 2); }
    1658          28 :   return ret;
    1659             : }
    1660             : 
    1661             : /* This calls the first downrun.
    1662             :    Then unravels the skeleton of the certificate.
    1663             :    This returns one of the following:
    1664             :    * a vector of the form [N, D, m, q, y]
    1665             :    * gen_0 (if N is composite)
    1666             :    * NULL (if FAIL)
    1667             : */
    1668             : static GEN
    1669          35 : ecpp_step1(GEN N, GEN param, GEN* X0)
    1670             : {
    1671          35 :   long depth = 0;
    1672          35 :   GEN downrun = N_downrun_NDinfomq(N, param, X0, &depth, 1);
    1673          35 :   if (downrun == NULL) return NULL;
    1674          28 :   return ecpp_flattencert(downrun, depth);
    1675             : }
    1676             : 
    1677             : /* This is the main function called.
    1678             :    The input is an integer N.
    1679             :    It uses the precomputation step0 done in ecpp_step0. */
    1680             : GEN
    1681          35 : ecpp0(GEN N, GEN param, GEN* X0)
    1682             : {
    1683             : 
    1684          35 :   pari_sp av = avma;
    1685             :   long i, j;
    1686             :   GEN step1;
    1687             :   GEN answer;
    1688             :   GEN Tv, Cv;
    1689             : 
    1690             :   /* Check if N is pseudoprime to begin with. */
    1691          35 :   if (X0 != NULL && !ispseudoprime(N, 0)) return gen_0;
    1692             : 
    1693             :   /* Check if we should even prove it. */
    1694          35 :   if (X0 != NULL && expi(N) < 64) return N;
    1695             : 
    1696             :   /* Timers and Counters */
    1697          35 :   Tv = mkvec4(zero_zv(6), zero_zv(3), zero_zv(3), zero_zv(2));
    1698          35 :   Cv = mkvec4(zero_zv(6), zero_zv(3), zero_zv(3), zero_zv(2));
    1699          35 :   if (X0) *X0 = mkvec3(Tv, Cv, zero_zv(1));
    1700             : 
    1701          35 :   step1 = ecpp_step1(N, param, X0);
    1702          35 :   if (step1 == NULL) return NULL;
    1703          28 :   if (typ(step1) != t_VEC) return gen_0;
    1704             : 
    1705          28 :   answer = ecpp_step2(step1, X0);
    1706          28 :   if (answer == NULL) pari_err(e_BUG, "ecpp");
    1707             : 
    1708         140 :   for (i = 1; i < lg(Tv); i++)
    1709             :   {
    1710         112 :     GEN v = gel(Tv, i);
    1711         112 :     long lgv = lg(v);
    1712         504 :     for (j = 1; j < lgv; j++)
    1713         392 :       dbg_mode() err_printf("\n   %c%ld: %16ld %16ld %16f", 'A'+i-1, j, umael3(*X0, 1, i, j), umael3(*X0, 2, i, j), (double)(umael3(*X0, 1, i, j))/(double)(umael3(*X0, 2, i, j)));
    1714             :   }
    1715          28 :   dbg_mode() err_printf("\n" ANSI_COLOR_BRIGHT_RED "\nFAILS: %16ld" ANSI_COLOR_RESET "\n", umael(*X0, 3, 1));
    1716             : 
    1717          28 :   if (X0) *X0 = gcopy(mkvec3(Tv, Cv, stoi(umael(*X0, 3, 1))));
    1718             : 
    1719          28 :   gerepileall(av, 1, &answer);
    1720          28 :   return answer;
    1721             : }
    1722             : 
    1723             : GEN
    1724          56 : ecpp(GEN N)
    1725             : {
    1726             :   long expiN, tunelen;
    1727             :   GEN param, answer, garbage, tune;
    1728          56 :   if (typ(N) != t_INT) pari_err_TYPE("ecpp", N);
    1729             : 
    1730             :   /* Check if N is pseudoprime to begin with. */
    1731          56 :   if (!ispseudoprime(N, 0)) return gen_0;
    1732             :   /* Check if we should even prove it. */
    1733          49 :   if (expi(N) < 64) return N;
    1734             : 
    1735          28 :   expiN = expi(N);
    1736          28 :   param = NULL;
    1737          28 :   answer = NULL;
    1738          28 :   garbage = NULL;
    1739             : 
    1740             :   /* tuning for 1500, 2500, 3500, 4500 bits */
    1741             :   /* ecpp is supposed to be faster than isprime on average if N is more than 1500 bits */
    1742          28 :   if (expiN <= 1500) tunelen = 1;
    1743           7 :   else if (expiN <= 2500) tunelen = 2;
    1744           0 :   else if (expiN <= 3500) tunelen = 3;
    1745           0 :   else tunelen = 4;
    1746             : 
    1747          28 :   tune = cgetg(tunelen+1, t_VEC);
    1748          28 :   if (tunelen >= 1) gel(tune, 1) = mkvecsmall3( 500, 24, 22);
    1749          28 :   if (tunelen >= 2) gel(tune, 2) = mkvecsmall3(1500, 32, 23);
    1750          28 :   if (tunelen >= 3) gel(tune, 3) = mkvecsmall3(1500, 32, 24);
    1751          28 :   if (tunelen >= 4) gel(tune, 4) = mkvecsmall3(3000, 40, 24);
    1752             : 
    1753          63 :   while (answer == NULL)
    1754             :   {
    1755             :     pari_timer T;
    1756          35 :     dbg_mode() timer_start(&T);
    1757          35 :     param = ecpp_param_set( tune );
    1758          35 :     dbg_mode() err_printf(ANSI_COLOR_BRIGHT_WHITE "\n%Ps" ANSI_COLOR_RESET, gel(tune, tunelen));
    1759          35 :     dbg_mode() err_printf(ANSI_COLOR_WHITE "  %8ld" ANSI_COLOR_RESET, timer_delay(&T));
    1760          35 :     answer = ecpp0(N, param, &garbage);
    1761          35 :     if (answer != NULL) break;
    1762           7 :     umael(tune, tunelen, 1) *= 2;
    1763           7 :     umael(tune, tunelen, 2) *= 2;
    1764           7 :     umael(tune, tunelen, 3)++;
    1765           7 :     if (umael(tune, tunelen, 3) > 26) umael(tune, tunelen, 3) = 26;
    1766             :   }
    1767          28 :   return answer;
    1768             : }
    1769             : 
    1770             : /*  Input: PARI ECPP Certificate
    1771             :    Output: Human-readable format.
    1772             : */
    1773             : static GEN
    1774           7 : cert_out(GEN x)
    1775             : {
    1776             :   long i;
    1777           7 :   long lgx = lg(x);
    1778             :   pari_str str;
    1779           7 :   str_init(&str, 1);
    1780           7 :   if (typ(x) == t_INT)
    1781             :   {
    1782           0 :     str_printf(&str, "%Ps is prime.\nIndeed, ispseudoprime(%Ps) = 1 and %Ps < 2^64.\n", x, x, x);
    1783           0 :     return strtoGENstr(str.string);
    1784             :   }
    1785          21 :   for (i = 1; i < lgx; i++)
    1786             :   {
    1787          14 :     GEN xi = gel(x, i);
    1788          14 :     str_printf(&str, "\n[%ld]\n", i);
    1789          14 :     str_printf(&str, " N = %Ps\n", cert_get_N(xi));
    1790          14 :     str_printf(&str, " t = %Ps\n", cert_get_t(xi));
    1791          14 :     str_printf(&str, " s = %Ps\n", cert_get_s(xi));
    1792          14 :     str_printf(&str, "a4 = %Ps\n", cert_get_a4(xi));
    1793          14 :     str_printf(&str, "D = %Ps\n", cert_get_D(xi));
    1794          14 :     str_printf(&str, "m = %Ps\n", cert_get_m(xi));
    1795          14 :     str_printf(&str, "q = %Ps\n", cert_get_q(xi));
    1796          14 :     str_printf(&str, "E = %Ps\n", cert_get_E(xi));
    1797          14 :     str_printf(&str, "P = %Ps\n", cert_get_P(xi));
    1798             :   }
    1799           7 :   return strtoGENstr(str.string);
    1800             : }
    1801             : 
    1802             : /*  Input: PARI ECPP Certificate
    1803             :    Output: Magma Certificate
    1804             :      Magma certificate looks like this
    1805             :      (newlines and extraneous spaces for clarity)
    1806             :      [*
    1807             :        [*
    1808             :          N,
    1809             :          |D|,
    1810             :          -1,
    1811             :          m,
    1812             :          [* a, b *],
    1813             :          [* x, y, 1 *],
    1814             :          [*
    1815             :            [* q, 1 *]
    1816             :          *]
    1817             :        *], ...
    1818             :      *]
    1819             : */
    1820             : static GEN
    1821           7 : magma_out(GEN x)
    1822             : {
    1823             :   long i;
    1824           7 :   long size = lg(x);
    1825             :   pari_str ret;
    1826           7 :   str_init(&ret, 1);
    1827           7 :   if (typ(x) == t_INT)
    1828             :   {
    1829           0 :     str_printf(&ret, "Operation not supported.");
    1830           0 :     return strtoGENstr(ret.string);
    1831             :   }
    1832           7 :   str_printf(&ret, "[* ");
    1833          21 :   for (i = 1; i<size; i++)
    1834             :   {
    1835          14 :     GEN xi = gel(x, i);
    1836          14 :     GEN E = cert_get_E(xi);
    1837          14 :     GEN P = cert_get_P(xi);
    1838          14 :     str_printf(&ret, "[* ");
    1839          14 :     str_printf(&ret, "%Ps, ", cert_get_N(xi));
    1840          14 :     str_printf(&ret, "%Ps, ", negi(cert_get_D(xi)));
    1841          14 :     str_printf(&ret, "-1, ");
    1842          14 :     str_printf(&ret, "%Ps, ", cert_get_m(xi));
    1843          14 :     str_printf(&ret, "[* %Ps, %Ps *], ", gel(E, 1), gel(E, 2));
    1844          14 :     str_printf(&ret, "[* %Ps, %Ps, 1 *], ", gel(P, 1), gel(P, 2));
    1845          14 :     str_printf(&ret, "[* [* %Ps, 1 *] *]", cert_get_q(xi));
    1846          14 :     str_printf(&ret, " *]");
    1847          14 :     if (i != size-1) str_printf(&ret, ", ");
    1848             :   }
    1849           7 :   str_printf(&ret, " *]");
    1850           7 :   return strtoGENstr(ret.string);
    1851             : }
    1852             : 
    1853             : /*  Input: &ret, label, value
    1854             :    Prints: label=(sign)hexvalue\n
    1855             :      where sign is + or -
    1856             :      hexvalue is value in hex, of the form: 0xABC123
    1857             : */
    1858             : static void
    1859          70 : primo_printval(pari_str *ret, const char* label, GEN value)
    1860             : {
    1861          70 :   str_printf(ret, label);
    1862          70 :   if (signe(value) >= 0) str_printf(ret, "=0x%P0X\n", value);
    1863          14 :   else str_printf(ret, "=-0x%P0X\n", negi(value));
    1864          70 : }
    1865             : 
    1866             : /*  Input: &ret, label, value
    1867             :    Prints: label=(sign)hexvalue\n
    1868             :      where sign is + or -
    1869             :      hexvalue is x in hex, of the form: 0xABC123
    1870             :        where |x| < N/2 and x = value mod N.
    1871             : */
    1872             : static void
    1873          14 : primo_printval_center(pari_str *ret, const char* label, GEN value, GEN N, GEN shif)
    1874             : {
    1875          14 :   primo_printval(ret, label, Fp_center(value, N, shif) );
    1876          14 : }
    1877             : 
    1878             : /*  Input: PARI ECPP Certificate
    1879             :    Output: PRIMO Certificate
    1880             : 
    1881             :    Let Y^2 = X^3 + aX + b be the elliptic curve over N
    1882             :    with the point (x, y) be the ones described in the
    1883             :    PARI certificate.
    1884             : 
    1885             :    For the case where J != 0, 1728, PRIMO asks for [J,T].
    1886             :    We obtain J easily using the well-known formula.
    1887             :    we obtain T by using the formula T = a/A * B/b * x
    1888             :    where A = 3J(1728-J) and B = 2J(1728-J)^2.
    1889             : 
    1890             :    For the case where J=0 or J=1728, PRIMO asks for [A,B,T].
    1891             :    We let A and B be a and b, respectively.
    1892             :    Consequently, T = x.
    1893             : */
    1894             : static GEN
    1895           7 : primo_out(GEN x)
    1896             : {
    1897             :   long i;
    1898           7 :   long size = (typ(x) == t_INT) ? 1 : lg(x);
    1899             :   pari_str ret;
    1900           7 :   str_init(&ret, 1);
    1901           7 :   str_printf(&ret, "[PRIMO - Primality Certificate]\n");
    1902           7 :   str_printf(&ret, "Format=4\n");
    1903           7 :   str_printf(&ret, "TestCount=%ld\n", size-1);
    1904           7 :   str_printf(&ret, "\n");
    1905           7 :   str_printf(&ret, "[Comments]\n");
    1906           7 :   str_printf(&ret, "Generated by PARI/GP\n");
    1907           7 :   str_printf(&ret, "http://pari.math.u-bordeaux.fr/\n");
    1908           7 :   str_printf(&ret, "\n");
    1909           7 :   str_printf(&ret, "[Candidate]\n");
    1910           7 :   if (typ(x) == t_INT)
    1911             :   {
    1912           0 :     str_printf(&ret, "N=0x%P0X\n", x);
    1913           0 :     return strtoGENstr(ret.string);
    1914             :   } else
    1915           7 :     str_printf(&ret, "N=0x%P0X\n", cert_get_N(gel(x, 1)));
    1916           7 :   str_printf(&ret, "\n");
    1917             : 
    1918          21 :   for (i = 1; i<size; i++)
    1919             :   {
    1920             :     GEN xi, N, Nover2;
    1921             :     long Ais0, Bis0;
    1922          14 :     str_printf(&ret, "[%ld]\n", i);
    1923          14 :     xi = gel(x, i);
    1924          14 :     N = cert_get_N(xi);
    1925          14 :     Nover2 = shifti(N, -1);
    1926          14 :     primo_printval(&ret, "S", cert_get_s(xi));
    1927          14 :     primo_printval(&ret, "W", cert_get_t(xi));
    1928          14 :     Ais0 = isintzero(cert_get_a4(xi));
    1929          14 :     Bis0 = isintzero(cert_get_a6(xi));
    1930          14 :     if (!Ais0 && !Bis0) { /* J != 0, 1728 */
    1931           0 :       primo_printval_center(&ret, "J", cert_get_J(xi), N, Nover2);
    1932           0 :       primo_printval(&ret, "T", cert_get_T(xi));
    1933          14 :     } else if (Ais0) { /* J == 0 */
    1934          14 :       primo_printval(&ret, "A", gen_0);
    1935          14 :       primo_printval_center(&ret, "B", cert_get_a6(xi), N, Nover2);
    1936          14 :       primo_printval(&ret, "T", cert_get_x(xi));
    1937             :     }else{ /* J == 1728 */
    1938           0 :       primo_printval_center(&ret, "A", cert_get_a4(xi), N, Nover2);
    1939           0 :       primo_printval(&ret, "B", gen_0);
    1940           0 :       primo_printval(&ret, "T", cert_get_x(xi));
    1941             :     }
    1942          14 :     str_printf(&ret, "\n");
    1943             :   }
    1944           7 :   return strtoGENstr(ret.string);
    1945             : }
    1946             : 
    1947             : /*  Input: N, q
    1948             :    Output: 1 if q > (N^1/4 + 1)^2, 0 otherwise.
    1949             : */
    1950             : static long
    1951         931 : Nq_isvalid(GEN N, GEN q)
    1952             : {
    1953         931 :   GEN qm1sqr = sqri(subiu(q, 1));
    1954         931 :   if (cmpii(qm1sqr,N) > 0) /* (q-1)^2 > N */
    1955             :   { /* (q-1)^4 + N^2 > 16Nq + 2N(q-1)^2 */
    1956         931 :     GEN a = addii(sqri(qm1sqr), sqri(N));
    1957         931 :     GEN b = addii(shifti(mulii(N, q), 4 ), mulii(shifti(N, 1), qm1sqr));
    1958         931 :     return (cmpii(a, b) > 0);
    1959             :   }
    1960           0 :   return 0;
    1961             : }
    1962             : 
    1963             : /*  Input: cert, trustbits=64
    1964             :    Output: 1 if cert is a valid PARI ECPP certificate
    1965             :      assuming that all q below 2^trustbits is prime
    1966             : */
    1967             : static long
    1968          56 : ecppisvalid0(GEN cert, long trustbits)
    1969             : {
    1970          56 :   pari_sp av = avma;
    1971             :   long i;
    1972          56 :   long lgcert = lg(cert);
    1973             :   GEN ql;
    1974          56 :   GEN N, W, m, s, q = gen_0, P, a;
    1975             :   GEN mP, sP;
    1976             :   GEN r;
    1977             : 
    1978          56 :   if (typ(cert) == t_INT)
    1979             :   {
    1980           7 :     if (expi(cert) > trustbits) return 0;
    1981           7 :     if (cmpii(cert, shifti(gen_1, 64)) > 0) return 0;
    1982           7 :     return ispseudoprime(cert, 0);
    1983             :   }
    1984             : 
    1985          49 :   if (typ(cert) != t_VEC) pari_err_TYPE("ecppisvalid", cert);
    1986             : 
    1987          49 :   if (lgcert <= 1) return 0;
    1988          49 :   if (lg(gel(cert, lgcert-1))-1 != 5) return 0;
    1989          49 :   ql = cert_get_q(gel(cert, lgcert-1));
    1990          49 :   if (expi(ql) > trustbits) return 0;
    1991          49 :   if (!ispseudoprime(ql, 0)) return 0;
    1992             : 
    1993         980 :   for (i = 1; i < lgcert; i++)
    1994             :   {
    1995         931 :     GEN certi = gel(cert, i);
    1996         931 :     if (lg(certi)-1 != 5) return 0;
    1997             : 
    1998         931 :     N = cert_get_N(certi);
    1999         931 :     if (typ(N) != t_INT) return 0;
    2000             :     /* Check if N > 1 */
    2001         931 :     if (signe(N) != 1) return 0;
    2002             :     /* Check if N is not divisible by 2 or 3 */
    2003         931 :     if (!isint1(gcdii(N, stoi(6)))) return 0;
    2004             :     /* Check if N matches the q from the previous entry. */
    2005         931 :     if (i > 1 && !equalii(N, q)) return 0;
    2006             : 
    2007         931 :     W = cert_get_t(certi);
    2008         931 :     if (typ(W) != t_INT) return 0;
    2009             :     /* Check if W^2 < 4N (Hasse interval) */
    2010         931 :     if (!(cmpii(sqri(W), shifti(N,2)) < 0)) return 0;
    2011             : 
    2012         931 :     s = cert_get_s(certi);
    2013         931 :     if (typ(s) != t_INT) return 0;
    2014         931 :     if (signe(s) < 1) return 0;
    2015             : 
    2016         931 :     m = cert_get_m(certi);
    2017         931 :     q = dvmdii(m, s, &r);
    2018             : 
    2019             :     /* Check m%s == 0 */
    2020         931 :     if (!isintzero(r)) return 0;
    2021             : 
    2022             :     /* Check q > (N^(1/4) + 1)^2 */
    2023         931 :     if (!Nq_isvalid(N, q)) return 0;
    2024             : 
    2025         931 :     a = cert_get_a4(certi);
    2026         931 :     if (typ(a) != t_INT) return 0;
    2027             : 
    2028         931 :     P = cert_get_P(certi);
    2029         931 :     if (typ(P) != t_VEC) return 0;
    2030         931 :     if (lg(P)-1 != 2) return 0;
    2031         931 :     P = FpE_to_FpJ(P);
    2032             : 
    2033             :     /* Check mP == 0 */
    2034         931 :     mP = FpJ_mul(P, m, a, N);
    2035         931 :     if (!FpJ_is_inf(mP)) return 0;
    2036             : 
    2037             :     /* Check sP != 0 and third component is coprime to N */
    2038         931 :     sP = FpJ_mul(P, s, a, N);
    2039         931 :     if (!isint1(gcdii(gel(sP, 3), N))) return 0;
    2040             :   }
    2041          49 :   avma = av; return 1;
    2042             : }
    2043             : 
    2044             : long
    2045          56 : ecppisvalid(GEN cert)
    2046          56 : { return ecppisvalid0(cert, 64); }
    2047             : 
    2048             : GEN
    2049          21 : ecppexport(GEN cert, long flag)
    2050             : {
    2051          21 :   if (!ecppisvalid(cert))
    2052           0 :     pari_err_TYPE("ecppexport - certificate not valid", cert);
    2053          21 :   switch(flag)
    2054             :   {
    2055           7 :     case 1: return magma_out(cert);
    2056           7 :     case 2: return primo_out(cert);
    2057           7 :     default: return cert_out(cert);
    2058             :   }
    2059             : }

Generated by: LCOV version 1.11