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

Generated by: LCOV version 1.11