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 21343-6216058) Lines: 845 917 92.1 %
Date: 2017-11-19 06:21:17 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        3701 : FpJ_is_inf(GEN P)
      76             : {
      77        3701 :   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        1385 : NmqEP_check(GEN N, GEN q, GEN E, GEN P, GEN s)
     830             : {
     831        1385 :   GEN a = gel(E, 1);
     832             :   GEN mP, sP;
     833        1385 :   sP = FpJ_mul(P, s, a, N);
     834        1385 :   if (FpJ_is_inf(sP)) return 0;
     835        1385 :   mP = FpJ_mul(sP, q, a, N);
     836        1385 :   if (FpJ_is_inf(mP)) return 1;
     837         720 :   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         268 : j0_find_g(GEN N)
     850             : {
     851             :   while (1)
     852             :   {
     853         268 :     GEN g = randomi(N);
     854         268 :     if (kronecker(g, N) != -1) continue;
     855         171 :     if (isint1(Fp_pow(g, diviuexact(subiu(N, 1), 3), N))) continue;
     856         224 :     return g;
     857         156 :   }
     858             : }
     859             : 
     860             : static GEN
     861        1385 : random_FpJ(GEN A, GEN B, GEN N)
     862             : {
     863        1385 :   GEN P = random_FpE(A, B, N);
     864        1385 :   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         417 :       for (i = 1; i < 6; i++)
     904             :       {
     905         417 :         B = Fp_mul(B, g, N); /* Bg */
     906         417 :         if (i % 3 == 0) continue;
     907         335 :         E = mkvec2(A, B);
     908         335 :         P = random_FpJ(A, B, N);
     909         335 :         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             : 
     948             :     /* C1: Find the appropriate class polynomial modulo N. */
     949         665 :     dbg_mode() timer_start(&ti);
     950         665 :     if (!equalii(D, Dprev)) HD = D_polclass(D, &db);
     951         665 :     dbg_mode() tt = timer_record(X0, "C1", &ti, 1);
     952         665 :     dbg_mode() err_printf(ANSI_COLOR_BRIGHT_GREEN "\n[ %3d | %4ld bits]" ANSI_COLOR_RESET, i, expi(N));
     953         665 :     dbg_mode() err_printf(ANSI_COLOR_GREEN "      D = %8Ps  poldeg = %4ld" ANSI_COLOR_RESET, D, degpol(HD));
     954         665 :     dbg_mode() if (equalii(D, Dprev)) err_printf("  %8ld", tt);
     955         665 :     dbg_mode() 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 the previous step. */
     958         665 :     dbg_mode() timer_start(&ti);
     959         665 :     rt = FpX_oneroot_split(HD, N);
     960         665 :     dbg_mode() tt = timer_record(X0, "C2", &ti, 1);
     961         665 :     dbg_mode() err_printf("  %8ld", tt);
     962             : 
     963             :     /* C3: Convert the root obtained from the previous step into the appropriate j-invariant. */
     964         665 :     dbg_mode() timer_start(&ti);
     965         665 :     J = NDinfor_find_J(N, Dinfo, rt);
     966         665 :     dbg_mode() tt = timer_record(X0, "C3", &ti, 1);
     967         665 :     dbg_mode() err_printf("  %8ld", tt);
     968             : 
     969             :     /* D1: Find an elliptic curve E with a point P satisfying the theorem. */
     970         665 :     dbg_mode() timer_start(&ti);
     971         665 :     s = diviiexact(m, q);
     972         665 :     EP = NDinfomqgJ_find_EP(N, Dinfo, m, q, g, J, s);
     973         665 :     dbg_mode() tt = timer_record(X0, "D1", &ti, 1);
     974         665 :     dbg_mode() err_printf("  %8ld", tt);
     975             : 
     976             :     /* D2: Compute for t and s */
     977         665 :     dbg_mode() timer_start(&ti);
     978         665 :     t = subii(addiu(N, 1), m); /* t = N+1-m */
     979         665 :     a4 = gmael(EP, 1, 1);
     980         665 :     P = FpJ_to_FpE(gel(EP, 2), N);
     981         665 :     dbg_mode() tt = timer_record(X0, "D2", &ti, 1);
     982         665 :     dbg_mode() err_printf("  %8ld", tt);
     983             : 
     984         665 :     gel(step2, i) = mkvec5(N, t, s, a4, P);
     985             : 
     986         665 :     Dprev = D;
     987             :   }
     988             : 
     989          28 :   return step2;
     990             : }
     991             : /* end of functions for step 2 */
     992             : 
     993             : /* start of functions for step 1 */
     994             : 
     995             : /* start of macros for step 1 */
     996             : 
     997             : /* earlyabort_pcdg is used when the downrun is told to NOT persevere
     998             :    it returns TRUE (i.e. do an early abort) if the polynomial degree
     999             :      of the discriminant in Dinfo[i] is larger than the param maxpcdg.
    1000             : */
    1001             : INLINE long
    1002      158039 : earlyabort_pcdg(GEN param, ulong maxpcdg, long i)
    1003             : {
    1004      158039 :   GEN x = ecpp_param_get_disclist(param);
    1005      158039 :   GEN Dinfo = gel(x, i);
    1006      158039 :   ulong pd = Dinfo_get_pd(Dinfo);
    1007      158039 :   return (pd > maxpcdg);
    1008             : }
    1009             : 
    1010             : /*  Input: vector whose components are [D, m]
    1011             :    Output: vector whose components are m
    1012             : */
    1013             : static GEN
    1014        1351 : Dmvec_to_mvec(GEN Dmvec)
    1015             : {
    1016        1351 :   long lgmvec = lg(Dmvec);
    1017        1351 :   GEN mvec = cgetg(lgmvec, t_VEC);
    1018             :   long i;
    1019        1351 :   for (i = 1; i < lgmvec; i++) gel(mvec, i) = gmael(Dmvec, i, 2);
    1020        1351 :   return mvec;
    1021             : }
    1022             : 
    1023             : /*  Input: vector v whose components are [D, m], vector w whose components are q
    1024             :    Output: vector whose components are [D, m, q]
    1025             : */
    1026             : static GEN
    1027        1351 : Dmvec_qvec_to_Dmqvec(GEN Dmvec, GEN qvec)
    1028             : {
    1029        1351 :   long lgDmqvec = lg(Dmvec);
    1030        1351 :   GEN Dmqvec = cgetg(lgDmqvec, t_VEC);
    1031             :   long i;
    1032       23205 :   for (i = 1; i < lgDmqvec; i++)
    1033             :   {
    1034       21854 :     GEN D = gmael(Dmvec, i, 1);
    1035       21854 :     GEN m = gmael(Dmvec, i, 2);
    1036       21854 :     GEN q = gel(qvec, i);
    1037       21854 :     gel(Dmqvec, i) = mkvec3(D, m, q);
    1038             :   }
    1039        1351 :   return Dmqvec;
    1040             : }
    1041             : 
    1042             : /*  Input: vector whose components are [D, m, q]
    1043             :    Output: vector whose components are q
    1044             : */
    1045             : static GEN
    1046        1351 : Dmqvec_to_qvec(GEN Dmqvec)
    1047             : {
    1048        1351 :   long lgqvec = lg(Dmqvec);
    1049        1351 :   GEN qvec = cgetg(lgqvec, t_VEC);
    1050             :   long i;
    1051        1351 :   for (i = 1; i < lgqvec; i++) gel(qvec, i) = gmael(Dmqvec, i, 3);
    1052        1351 :   return qvec;
    1053             : }
    1054             : 
    1055             : /* This initializes sqrtlist as a vector [A, B]
    1056             :      where A is a t_VEC of gen_0's
    1057             :        and B is a t_VECSMALL of 0's,
    1058             :    both of which are as long as primelist. */
    1059             : INLINE void
    1060         686 : primelist_sqrtlist_init(GEN primelist, GEN *sqrtlist)
    1061             : {
    1062         686 :   long l = lg(primelist)-1;
    1063         686 :   *sqrtlist = mkvec2(zerovec(l), zero_zv(l));
    1064         686 : }
    1065             : 
    1066             : /* This returns the square root modulo N
    1067             :      of the ith entry of the primelist.
    1068             :    If this square root is already available on sqrtlist,
    1069             :      then simply return it.
    1070             :    Otherwise, compute it, save it and return it.
    1071             :    y is a quadratic non-residue that is needed
    1072             :      somehow in the algorithm for taking square roots modulo N.
    1073             : */
    1074             : static GEN
    1075       64988 : p_find_primesqrt(GEN N, GEN* X0, GEN primelist, GEN sqrtlist, long i, GEN g)
    1076             : {
    1077             :   pari_timer ti;
    1078       64988 :   long p = uel(primelist, i);
    1079       64988 :   if (isintzero(gmael(sqrtlist,1,i)))
    1080             :   {
    1081       16009 :     dbg_mode() err_printf(ANSI_COLOR_MAGENTA "S" ANSI_COLOR_RESET);
    1082             :     /* A4: Get the square root of a prime factor of D. */
    1083       16009 :     dbg_mode() timer_start(&ti);
    1084       16009 :     gmael(sqrtlist, 1, i) = Fp_sqrt_i(stoi(p), g, N); /* NULL if invalid. */
    1085       16009 :     dbg_mode() timer_record(X0, "A4", &ti, 1);
    1086             :   }
    1087       64988 :   return gmael(sqrtlist, 1, i);
    1088             : }
    1089             : 
    1090             : /* This finds the legit square root of D modulo N where D is the discriminant in Dinfo.
    1091             :    This algorithm finds the square root modulo N of each prime p dividing D.
    1092             :    It then assembles the actual square root of D by multiplying the prime square roots.
    1093             : */
    1094             : static GEN
    1095       34027 : D_find_discsqrt(GEN N, GEN param, GEN *X0, GEN Dinfo, GEN sqrtlist, GEN g)
    1096             : {
    1097       34027 :   GEN discsqrt = gen_1;
    1098       34027 :   GEN Dfac = Dinfo_get_Dfac(Dinfo);
    1099       34027 :   long lgDfac = lg(Dfac);
    1100             :   long i;
    1101       34027 :   GEN primelist = ecpp_param_get_primelist(param);
    1102       99015 :   for (i = 1; i < lgDfac; i++)
    1103             :   {
    1104       64988 :     GEN sqrtfi = p_find_primesqrt(N, X0, primelist, sqrtlist, uel(Dfac, i), g);
    1105       64988 :     if (sqrtfi == NULL) return NULL; /* Only happens when N is composite. */
    1106       64988 :     discsqrt = Fp_mul(discsqrt, sqrtfi, N);
    1107             :   }
    1108       34027 :   return discsqrt;
    1109             : }
    1110             : 
    1111             : /* Given a solution U, V to U^2 + |D|V^2 = 4N, this find all the possible
    1112             :      cardinalities of the elliptic curves over the finite field F_N
    1113             :      whose endomorphism ring is isomorphic to the maximal order
    1114             :      in the imaginary quadratic number field K = Q(sqrt(D)). ???
    1115             : */
    1116             : static GEN
    1117        9884 : NUV_find_mvec(GEN N, GEN U, GEN V, long wD)
    1118             : {
    1119        9884 :   GEN Nplus1 = addiu(N, 1);
    1120             :   GEN m;
    1121        9884 :   GEN u = U;
    1122        9884 :   GEN mvec = cgetg(wD + 1, t_VEC);
    1123             :   long i;
    1124       31738 :   for (i = 0; i < wD; i++)
    1125             :   {
    1126       21854 :     if (i%2 == 0)
    1127             :     {
    1128       10927 :       if (wD == 4 && i==2) u = shifti(V, 1);
    1129       10927 :       if (wD == 6 && i==2) u = shifti(submuliu(U, V, 3), -1);
    1130       10927 :       if (wD == 6 && i==4) u = shifti(addmuliu(U, V, 3), -1);
    1131       10927 :       m = addii(Nplus1, u);
    1132             :     } else
    1133       10927 :       m = subii(Nplus1, u);
    1134       21854 :     gel(mvec, i + 1) = m;
    1135             :   }
    1136        9884 :   return mvec;
    1137             : }
    1138             : 
    1139             : /* This populates Dmbatch with Dmvec's -- vectors whose components are of the form [D, m],
    1140             :      where m is a cardinalities of an elliptic curve with endomorphism ring isomorphic to
    1141             :      the maximal order of the imaginary quadratic number field K = Q(sqrt(D)).
    1142             :    It returns 0 if:
    1143             :      * D (of Dinfo) is not a quadratic residue mod N
    1144             :      * Any of the p* dividing D is not a quadratic residue mod N
    1145             :      * Cornacchia cannot find a solution U^2 + |D|V^2 = 4N.
    1146             :    Otherwise, it returns the number of cardinalities added to Dmbatch.
    1147             :    Moreover, if N is determined to be composite, it sets failflag to 1.
    1148             :    Finally, sqrtlist and y are used to help compute the square root modulo N of D+.
    1149             : */
    1150             : static long
    1151      319998 : D_collectcards(GEN N, GEN param, GEN* X0, GEN Dinfo, GEN sqrtlist, GEN g, GEN Dmbatch, long* failflag)
    1152             : {
    1153             : 
    1154             :   long i, j;
    1155             :   GEN U, V;
    1156             :   long corn_succ;
    1157             :   long wD;
    1158             :   GEN Dfac;
    1159             :   long lgDfac;
    1160             :   GEN sqrtofDmodN;
    1161      319998 :   GEN primelist = ecpp_param_get_primelist(param);
    1162             :   GEN mvec;
    1163             : 
    1164             :   /* Unpacking info on the discriminant. */
    1165      319998 :   long D = gel(Dinfo, 1)[1];
    1166             : 
    1167             :   /* timers / counters */
    1168             :   pari_timer ti;
    1169             :   long kronDN;
    1170             : 
    1171             :   /* A1: Check (D|N) = 1. */
    1172      319998 :   dbg_mode() timer_start(&ti);
    1173      319998 :   kronDN = krosi(D, N);
    1174      319998 :   dbg_mode() timer_record(X0, "A1", &ti, 1);
    1175      319998 :   switch(kronDN) {
    1176           0 :      case 0: *failflag = 1;
    1177      159390 :     case -1: return 0;
    1178             :   }
    1179             : 
    1180             :   /* A2: Check (p*|N) = 1 for all p.
    1181             :      This is equivalent to checking (N|p) = 1. */
    1182      160608 :   dbg_mode() timer_start(&ti);
    1183      160608 :   Dfac = Dinfo_get_Dfac(Dinfo);
    1184      160608 :   lgDfac = lg(Dfac);
    1185      251538 :   for (i = 1; i < lgDfac; i++)
    1186             :   {
    1187      217511 :     long p = index_to_p(uel(Dfac, i), primelist);
    1188      217511 :     if (krosi(p, N) != 1) return 0;
    1189             :   }
    1190       34027 :   dbg_mode() timer_record(X0, "A2", &ti, 1);
    1191             : 
    1192             :   /* A3: Get square root of D mod N. */
    1193             :   /* If sqrtofDmodN is NULL, then N is composite. */
    1194       34027 :   dbg_mode() timer_start(&ti);
    1195       34027 :   sqrtofDmodN = D_find_discsqrt(N, param, X0, Dinfo, sqrtlist, g);
    1196       34027 :   if (sqrtofDmodN == NULL) pari_err_BUG("D_find_discsqrt");
    1197       34027 :   dbg_mode() if (!equalii(Fp_sqr(sqrtofDmodN, N), addis(N, D)) /* D mod N, D < 0*/ ) pari_err_BUG("D_find_discsqrt");
    1198       34027 :   dbg_mode() timer_record(X0, "A3", &ti, 1);
    1199             : 
    1200             :   /* A5: Use the square root to use cornacchia to find the solution to U^2 + |D|V^2 = 4N. */
    1201       34027 :   dbg_mode() timer_start(&ti);
    1202       34027 :   corn_succ = cornacchia2_sqrt( absi(stoi(D)), N, sqrtofDmodN, &U, &V);
    1203       34027 :   dbg_mode() timer_record(X0, "A5", &ti, 1);
    1204       34027 :   if (!corn_succ) {
    1205       24143 :     dbg_mode() err_printf(ANSI_COLOR_YELLOW "c" ANSI_COLOR_RESET);
    1206       24143 :     return 0;
    1207             :   }
    1208             : 
    1209             :   /* We're sticking with this D. */
    1210        9884 :   dbg_mode() err_printf(ANSI_COLOR_BRIGHT_YELLOW "D" ANSI_COLOR_RESET);
    1211             : 
    1212        9884 :   dbg_mode() timer_start(&ti);
    1213             :   /* A6: Collect the w(D) possible cardinalities of elliptic curves over F_N whose discriminant is D. */
    1214        9884 :   wD = D_get_wD(D);
    1215        9884 :   mvec = NUV_find_mvec(N, U, V, wD);
    1216        9884 :   for (j = 1; j < lg(mvec); j++) vectrunc_append(Dmbatch, mkvec2(Dinfo, gel(mvec, j)) );
    1217        9884 :   dbg_mode() timer_record(X0, "A6", &ti, 1);
    1218             : 
    1219        9884 :   return wD;
    1220             : }
    1221             : 
    1222             : /* Compute (S(16N, 2) + S(4096N, 4) + 4)\4 + 1
    1223             :      where S is the nth root rounded down.
    1224             :    This is at most one more than (N^1/4 + 1)^2.
    1225             : */
    1226             : static GEN
    1227        1351 : ecpp_qlo(GEN N)
    1228             : {
    1229        1351 :   GEN a = sqrtnint(shifti(N, 4), 2);
    1230        1351 :   GEN b = sqrtnint(shifti(N, 12), 4);
    1231        1351 :   GEN c = shifti(gen_1, 2);
    1232        1351 :   GEN d = addii(addii(a, b), c);
    1233        1351 :   GEN e = shifti(d, -2);
    1234        1351 :   GEN f = addiu(e, 1);
    1235        1351 :   return f;
    1236             : }
    1237             : 
    1238             : /* E is &qlo, use for zv_binsearch */
    1239             : static long
    1240        6923 : lessthan_qlo(void* E, GEN q)
    1241             : {
    1242        6923 :   GEN qlo = *((GEN*)E);
    1243        6923 :   return (cmpii(q, qlo) < 0);
    1244             : }
    1245             : 
    1246             : /* E is &goal, use for zv_binsearch */
    1247             : static long
    1248        6097 : gained_bits(void* E, GEN q)
    1249             : {
    1250        6097 :   long goal = *((long*)E);
    1251        6097 :   return (expi(q) <= goal);
    1252             : }
    1253             : 
    1254             : /*  Input: Dmqvec
    1255             :    Output: Dmqvec such that q satisfies
    1256             :      (N^1/4 + 1)^2 < q < N/2
    1257             : */
    1258             : static GEN
    1259        1351 : Dmqvec_slice_Dmqvec(GEN N, GEN Dmqvec)
    1260             : {
    1261        1351 :   pari_sp av = avma;
    1262             :   GEN qlo;
    1263             :   GEN qvec;
    1264             :   long lo_ind, bitgain, hi_ind, goal;
    1265             : 
    1266             :   /* Dmqvec is sorted according to q */
    1267        1351 :   Dmqvec = gen_sort(Dmqvec, NULL, &sort_Dmq_by_q);
    1268        1351 :   qvec = Dmqvec_to_qvec(Dmqvec);
    1269             : 
    1270        1351 :   qlo = ecpp_qlo(N);
    1271        1351 :   lo_ind = zv_binsearch0(&qlo, &lessthan_qlo, qvec); lo_ind++;
    1272        1351 :   if (lo_ind >= lg(qvec)) { avma = av; return NULL; }
    1273             : 
    1274        1351 :   bitgain = 1;
    1275        1351 :   goal = expi(N)-bitgain;
    1276        1351 :   hi_ind = zv_binsearch0(&goal, &gained_bits, qvec);
    1277        1351 :   if (hi_ind == 0) { avma = av; return NULL; }
    1278             : 
    1279        1351 :   Dmqvec = vecslice(Dmqvec, lo_ind, hi_ind);
    1280        1351 :   gerepileall(av, 1, &Dmqvec);
    1281        1351 :   return Dmqvec;
    1282             : }
    1283             : 
    1284             : /* Given a vector mvec of mi's,
    1285             :    This simultaneously removes all prime factors of each mi less then BOUND_PRIMORIAL
    1286             :    This implements Franke 2004: Proving the Primality of Very Large Numbers with fastECPP. */
    1287             : static GEN
    1288        1351 : mvec_batchfactor_qvec(GEN mvec, GEN primorial)
    1289             : {
    1290        1351 :   pari_sp av = avma;
    1291             :   long i;
    1292             :   /* Obtain the product tree of mvec. */
    1293        1351 :   GEN leaf = Z_ZV_mod_tree(primorial, mvec, ZV_producttree(mvec));
    1294             : 
    1295             :   /* Go through each leaf and remove small factors. */
    1296       23205 :   for (i = 1; i < lg(mvec); i++)
    1297             :   {
    1298       21854 :     GEN m = gel(mvec, i);
    1299      110117 :     while (!isint1(gel(leaf, i)) )
    1300             :     {
    1301       66409 :       gel(leaf, i) = gcdii( m, gel(leaf,i) );
    1302       66409 :       m = diviiexact( m, gel(leaf,i) );
    1303             :     }
    1304       21854 :     gel(mvec, i) = m;
    1305             :   }
    1306        1351 :   gerepileall(av, 1, &mvec);
    1307        1351 :   return mvec;
    1308             : }
    1309             : 
    1310             : /*  Input: Dmvec
    1311             :    Output: Dmqvec
    1312             :    Uses mvec_batchfactor_qvec to produce the output.
    1313             : */
    1314             : static GEN
    1315        1351 : Dmvec_batchfactor_Dmqvec(GEN Dmvec, GEN primorial)
    1316             : {
    1317        1351 :   pari_sp av = avma;
    1318        1351 :   GEN mvec = Dmvec_to_mvec(Dmvec);
    1319        1351 :   GEN qvec = mvec_batchfactor_qvec(mvec, primorial);
    1320        1351 :   GEN Dmqvec = Dmvec_qvec_to_Dmqvec(Dmvec, qvec);
    1321        1351 :   gerepileall(av, 1, &Dmqvec);
    1322        1351 :   return Dmqvec;
    1323             : }
    1324             : 
    1325             : /* tunevec = [maxsqrt, maxpcdg, tdivexp]
    1326             :      and is a shorthand for specifying the parameters of ecpp
    1327             : */
    1328             : static GEN
    1329        2723 : tunevec(GEN N, GEN param)
    1330             : {
    1331        2723 :   long e = expi(N);
    1332        2723 :   GEN T = gel(param,3);
    1333        2723 :   if      (e <= 1500) return gel(T,1);
    1334          91 :   else if (e <= 2500) return gel(T,2);
    1335           0 :   else if (e <= 3500) return gel(T,3);
    1336           0 :   return gel(T,4);
    1337             : }
    1338             : static long
    1339        2723 : tunevec_tdivbd(GEN N, GEN param)
    1340             : {
    1341        2723 :   return uel(tunevec(N, param), 3);
    1342             : }
    1343             : static long
    1344         686 : tunevec_batchsize(GEN N, GEN param)
    1345             : {
    1346         686 :   return (28-tunevec_tdivbd(N, param) < 0 ? 0 : 28-tunevec_tdivbd(N, param));
    1347             : }
    1348             : 
    1349             : /*  Input: Dmbatch (from the downrun)
    1350             :    Output: Dmqvec
    1351             : */
    1352             : static GEN
    1353        1351 : Dmbatch_factor_Dmqvec(GEN N, GEN* X0, GEN Dmbatch, GEN param)
    1354             : {
    1355             :   pari_timer ti;
    1356        1351 :   GEN curr_primorial = ecpp_param_get_primorial(param, tunevec_tdivbd(N, param) - 21);
    1357             :   GEN Dmqvec;
    1358             : 
    1359             :   /* B1: Factor by batch. */
    1360        1351 :   dbg_mode() timer_start(&ti);
    1361        1351 :   Dmqvec = Dmvec_batchfactor_Dmqvec(Dmbatch, curr_primorial);
    1362        1351 :   dbg_mode() timer_record(X0, "B1", &ti, 1);
    1363             : 
    1364             :   /* B2: For each batch,
    1365             :          remove cardinalities lower than (N^(1/4)+1)^2
    1366             :          and cardinalities in which we didn't win enough bits. */
    1367        1351 :   dbg_mode() timer_start(&ti);
    1368        1351 :   Dmqvec = Dmqvec_slice_Dmqvec(N, Dmqvec);
    1369        1351 :   dbg_mode() timer_record(X0, "B2", &ti, lg(Dmqvec)-1);
    1370             : 
    1371             :   /* If nothing is left after B2, return NULL */
    1372        1351 :   if (Dmqvec == NULL) return NULL;
    1373             : 
    1374        1351 :   return Dmqvec;
    1375             : }
    1376             : 
    1377             : /* Dmq macros */
    1378             : INLINE GEN
    1379       15778 : Dmq_get_Dinfo(GEN Dmq)
    1380       15778 : { return gel(Dmq, 1); }
    1381             : 
    1382             : INLINE GEN
    1383       15778 : Dmq_get_m(GEN Dmq)
    1384       15778 : { return gel(Dmq, 2); }
    1385             : 
    1386             : INLINE GEN
    1387       31556 : Dmq_get_q(GEN Dmq)
    1388       31556 : { return gel(Dmq, 3); }
    1389             : 
    1390             : INLINE long
    1391           0 : Dmq_get_cnum(GEN Dmq)
    1392           0 : { return gmael(Dmq, 1, 1)[2]; }
    1393             : 
    1394             : /*  Input: Dmq (where q does not have small prime factors)
    1395             :    Output: whether q is pseudoprime or not
    1396             : */
    1397             : static long
    1398       15778 : Dmq_isgoodq(GEN Dmq, GEN* X0)
    1399             : {
    1400             :   pari_timer ti;
    1401             :   long is_pseudoprime;
    1402       15778 :   GEN q = Dmq_get_q(Dmq);
    1403             : 
    1404             :   /* B3: Check pseudoprimality of each q on the list. */
    1405       15778 :   dbg_mode() timer_start(&ti);
    1406       15778 :   is_pseudoprime = BPSW_psp_nosmalldiv(q);
    1407       15778 :   dbg_mode() timer_record(X0, "B3", &ti, 1);
    1408       15778 :   return is_pseudoprime; /* did not find for this m */
    1409             : }
    1410             : 
    1411             : /*  Input: N
    1412             :    Output: [ NDmqg, therest ]
    1413             :      NDmqg is a vector of the form [N, D, m, q, g].
    1414             :      therest is a vector of the form the same as the output OR [N].
    1415             :    This is the downrun.
    1416             :    It tries to find [N, D, m, q, g].
    1417             :    If N is small, terminate.
    1418             :    It finds a quadratic non-residue g.
    1419             :    It starts with finding the square root of D modulo N.
    1420             :    It then uses this square root for cornacchia's algorithm.
    1421             :    If there is a solution to U^2 + |D|V^2 = 4N, use it to find candidates for m.
    1422             :    It continues until you get a batch of size at least t = log2(N)/2^log2bsize
    1423             :      or if there's no more discriminants left on the list.
    1424             :    It factors the m's of the batch and produces the q's.
    1425             :    If one of the q's are pseudoprime, then call this function again with N = q.
    1426             : */
    1427             : static GEN
    1428         714 : N_downrun_NDinfomq(GEN N, GEN param, GEN *X0, long *depth, long persevere)
    1429             : {
    1430         714 :   pari_sp ave = avma;
    1431             :   pari_timer T;
    1432             :   long i, j;
    1433         714 :   long expiN = expi(N);
    1434         714 :   long persevere_next = 0;
    1435             :   long lgdisclist, t, numcard;
    1436             :   ulong maxpcdg;
    1437             :   GEN primelist, disclist, sqrtlist, g, Dmbatch;
    1438         714 :   long FAIL = 0;
    1439             : 
    1440         714 :   dbg_mode() timer_start(&T);
    1441             : 
    1442             :   /* Unpack trustbits. */
    1443         714 :   if (expiN < 64) return gerepilecopy(ave, mkvec(N));
    1444             : 
    1445             :   /* This means we're going down the tree. */
    1446         686 :   *depth += 1;
    1447             : 
    1448             :   /* Unpack maxpcdg. */
    1449         686 :   maxpcdg = ecpp_param_get_maxpcdg(param);
    1450             : 
    1451             :   /* Unpack primelist, disclist. */
    1452         686 :   primelist = ecpp_param_get_primelist(param);
    1453         686 :   disclist = ecpp_param_get_disclist(param);
    1454         686 :   lgdisclist = lg(disclist);
    1455             : 
    1456             :   /* Initialize sqrtlist for this N. */
    1457         686 :   primelist_sqrtlist_init(primelist, &sqrtlist);
    1458             : 
    1459             :   /* Precomputation for batch size t. */
    1460             :   /* Tuning! */
    1461         686 :   t = expiN >> tunevec_batchsize(N, param);
    1462         686 :   if (t < 1) t = 1;
    1463             : 
    1464             :   /* Precomputation for taking square roots.
    1465             :        g will be needed for Fp_sqrt_i
    1466             :   */
    1467         686 :   g = Fp_2gener(N);
    1468         686 :   if (g == NULL) return gen_0; /* Composite if this happens. */
    1469             : 
    1470             :   /* Print the start of this iteration. */
    1471         686 :   dbg_mode() if (!persevere) err_printf(ANSI_COLOR_BRIGHT_CYAN "\n[ %3d | %4ld bits] " ANSI_COLOR_RESET, *depth, expiN, persevere);
    1472         686 :   dbg_mode() if (persevere) err_printf(ANSI_COLOR_BRIGHT_CYAN "\n< %3d | %4ld bits> " ANSI_COLOR_RESET, *depth, expiN, persevere);
    1473             : 
    1474             :   /* Initialize Dmbatch
    1475             :        It will be populated with candidate cardinalities on Phase I (until its length reaches at least t).
    1476             :        Its elements will be factored on Phase II.
    1477             :        We allocate a vectrunc_init of t+7.
    1478             :   */
    1479         686 :   Dmbatch = vectrunc_init(t+7);
    1480             : 
    1481             :   /* Number of cardinalities collected so far.
    1482             :      Should always be equal to lg(Dmbatch)-1. */
    1483         686 :   numcard = 0;
    1484             : 
    1485             :   /* i determines which discriminant we are considering. */
    1486         686 :   i = 1;
    1487             : 
    1488        2058 :   while (!FAIL)
    1489             :   {
    1490             :     pari_timer F;
    1491        1358 :     long last_i = i, failflag;
    1492             :     GEN Dmqlist;
    1493             :     long lgDmqlist;
    1494             : 
    1495             :     /* Dmbatch begins "empty", but we keep the allocated memory. */
    1496        1358 :     setlg(Dmbatch, 1);
    1497        1358 :     numcard = 0;
    1498             : 
    1499             :     /* PHASE I:
    1500             :        Go through the D's and search for canidate m's.
    1501             :        We fill up Dmbatch until there are at least t elements.
    1502             :     */
    1503             : 
    1504        1358 :     failflag = 0;
    1505      321377 :     while (i < lgdisclist )
    1506             :     {
    1507             :       GEN Dinfo;
    1508      320012 :       if (!persevere && earlyabort_pcdg(param, maxpcdg, i)) { FAIL = 1; break; }
    1509      319998 :       Dinfo = gel(disclist, i);
    1510      319998 :       numcard += D_collectcards(N, param, X0, Dinfo, sqrtlist, g, Dmbatch, &failflag);
    1511      320663 :       if (failflag) return gen_0;
    1512      319998 :       last_i = i++;
    1513      319998 :       if (numcard >= t) break;
    1514             :     }
    1515             : 
    1516             :     /* If we have exhausted disclist and there are no cardinalities to be factored. */
    1517        1365 :     if (FAIL && numcard <= 0) break;
    1518        1351 :     if (i >= lgdisclist && numcard <= 0) break;
    1519             : 
    1520             :     /* PHASE II:
    1521             :        Find the corresponding q's for the m's found.
    1522             :        We use Dmbatch.
    1523             :     */
    1524             : 
    1525             :     /* Find coresponding q of the candidate m's. */
    1526        1351 :     dbg_mode() timer_start(&F);
    1527        1351 :     Dmqlist = Dmbatch_factor_Dmqvec(N, X0, Dmbatch, param);
    1528             : 
    1529             :     /* If none left, move to the next discriminant. */
    1530        1351 :     if (Dmqlist == NULL) continue;
    1531             : 
    1532        1351 :     lgDmqlist = lg(Dmqlist);
    1533             : 
    1534             :     /* If we are cheating by adding class numbers, sort by class number instead. */
    1535        1351 :     if (Dinfo_get_pd(gel(disclist, last_i)) > maxpcdg)
    1536          28 :       Dmqlist = gen_sort(Dmqlist, NULL, &sort_Dmq_by_cnum);
    1537             : 
    1538             :     /* Check pseudoprimality before going down. */
    1539       16464 :     for (j = 1; j < lgDmqlist; j++)
    1540             :     {
    1541             :       GEN NDinfomq;
    1542       15778 :       GEN Dmq = gel(Dmqlist, j);
    1543       15778 :       GEN Dinfo = Dmq_get_Dinfo(Dmq);
    1544       15778 :       GEN m = Dmq_get_m(Dmq);
    1545       15778 :       GEN q = Dmq_get_q(Dmq);
    1546             :       GEN ret;
    1547             :       long a;
    1548       15778 :       if (expiN - expi(q) < 1)
    1549             :       {
    1550           0 :         dbg_mode() err_printf(ANSI_COLOR_BRIGHT_RED "  x" ANSI_COLOR_RESET);
    1551           0 :         continue;
    1552             :       }
    1553       15778 :       dbg_mode() err_printf(ANSI_COLOR_WHITE "." ANSI_COLOR_RESET);
    1554       15778 :       a = Dmq_isgoodq(Dmq, X0);
    1555       15778 :       if (!a) continue;
    1556             : 
    1557         679 :       dbg_mode() err_printf(ANSI_COLOR_BRIGHT_BLUE "  %ld" ANSI_COLOR_RESET, Dmq_get_cnum(Dmq));
    1558         679 :       dbg_mode() err_printf(ANSI_COLOR_BRIGHT_RED "\n       %5ld bits " ANSI_COLOR_RESET, expi(q)-expiN);
    1559             : 
    1560             :       /* Cardinality is pseudoprime. Call the next downrun! */
    1561         679 :       ret = N_downrun_NDinfomq(q, param, X0, depth, persevere_next);
    1562             : 
    1563             :       /* That downrun failed. */
    1564         679 :       if (ret == NULL) {
    1565          14 :         dbg_mode() if (!persevere) err_printf(ANSI_COLOR_CYAN "\n[ %3d | %4ld bits] " ANSI_COLOR_RESET, *depth, expiN, persevere);
    1566          14 :         dbg_mode() if (persevere) err_printf(ANSI_COLOR_CYAN "\n< %3d | %4ld bits> " ANSI_COLOR_RESET, *depth, expiN, persevere);
    1567          14 :         continue;
    1568             :       }
    1569             : 
    1570             :       /* That downrun succeeded. */
    1571         665 :       NDinfomq = mkvec5(N, Dinfo, m, q, g);
    1572         665 :       return gerepilecopy(ave, mkvec2(NDinfomq, ret));
    1573             :     }
    1574             : 
    1575             :     /* We have exhausted all the discriminants. */
    1576         686 :     if (i >= lgdisclist) FAIL = 1;
    1577             : 
    1578         686 :     if (Dinfo_get_pd(gel(disclist, last_i)) > maxpcdg)
    1579             :     {
    1580          28 :       dbg_mode() err_printf(ANSI_COLOR_BRIGHT_RED "  !" ANSI_COLOR_RESET);
    1581          28 :       persevere_next = 1;
    1582             :     }
    1583             : 
    1584             :   }
    1585             : 
    1586             :   /* FAILED: Out of discriminants. */
    1587          21 :   if (X0) umael(*X0, 3, 1)++; /* FAILS++ */
    1588          21 :   (*depth)--;
    1589          21 :   dbg_mode() err_printf(ANSI_COLOR_BRIGHT_RED "  X" ANSI_COLOR_RESET);
    1590          21 :   return NULL;
    1591             : }
    1592             : 
    1593             : /*  Input: the output of the (recursive) downrun function
    1594             :    Output: a vector whose components are [N, D, m, q, g]
    1595             : */
    1596             : static GEN
    1597          28 : ecpp_flattencert(GEN notflat, long depth)
    1598             : {
    1599          28 :   long i, lgret = depth+1;
    1600          28 :   GEN ret = cgetg(lgret, t_VEC);
    1601          28 :   GEN x = notflat;
    1602          28 :   if (typ(x) != t_VEC) return gen_0;
    1603          28 :   for (i = 1; i < lgret; i++) { gel(ret, i) = gel(x, 1); x = gel(x, 2); }
    1604          28 :   return ret;
    1605             : }
    1606             : 
    1607             : /* This calls the first downrun.
    1608             :    Then unravels the skeleton of the certificate.
    1609             :    This returns one of the following:
    1610             :    * a vector of the form [N, D, m, q, y]
    1611             :    * gen_0 (if N is composite)
    1612             :    * NULL (if FAIL)
    1613             : */
    1614             : static GEN
    1615          35 : ecpp_step1(GEN N, GEN param, GEN* X0)
    1616             : {
    1617          35 :   long depth = 0;
    1618          35 :   GEN downrun = N_downrun_NDinfomq(N, param, X0, &depth, 1);
    1619          35 :   if (downrun == NULL) return NULL;
    1620          28 :   return ecpp_flattencert(downrun, depth);
    1621             : }
    1622             : 
    1623             : /* This is the main function called.
    1624             :    The input is an integer N.
    1625             :    It uses the precomputation step0 done in ecpp_step0. */
    1626             : GEN
    1627          35 : ecpp0(GEN N, GEN param, GEN* X0)
    1628             : {
    1629             : 
    1630          35 :   pari_sp av = avma;
    1631             :   long i, j;
    1632             :   GEN step1;
    1633             :   GEN answer;
    1634             :   GEN Tv, Cv;
    1635             : 
    1636             :   /* Check if N is pseudoprime to begin with. */
    1637          35 :   if (X0 != NULL && !ispseudoprime(N, 0)) return gen_0;
    1638             : 
    1639             :   /* Check if we should even prove it. */
    1640          35 :   if (X0 != NULL && expi(N) < 64) return N;
    1641             : 
    1642             :   /* Timers and Counters */
    1643          35 :   Tv = mkvec4(zero_zv(6), zero_zv(3), zero_zv(3), zero_zv(2));
    1644          35 :   Cv = mkvec4(zero_zv(6), zero_zv(3), zero_zv(3), zero_zv(2));
    1645          35 :   if (X0) *X0 = mkvec3(Tv, Cv, zero_zv(1));
    1646             : 
    1647          35 :   step1 = ecpp_step1(N, param, X0);
    1648          35 :   if (step1 == NULL) return NULL;
    1649          28 :   if (typ(step1) != t_VEC) return gen_0;
    1650             : 
    1651          28 :   answer = ecpp_step2(step1, X0);
    1652          28 :   if (answer == NULL) pari_err(e_BUG, "ecpp");
    1653             : 
    1654         140 :   for (i = 1; i < lg(Tv); i++)
    1655             :   {
    1656         112 :     GEN v = gel(Tv, i);
    1657         112 :     long lgv = lg(v);
    1658         504 :     for (j = 1; j < lgv; j++)
    1659         392 :       dbg_mode() err_printf("\n   %c%ld: %16ld %16ld %16f", 'A'+i-1, j, umael3(*X0, 1, i, j), umael3(*X0, 2, i, j), (double)(umael3(*X0, 1, i, j))/(double)(umael3(*X0, 2, i, j)));
    1660             :   }
    1661          28 :   dbg_mode() err_printf("\n" ANSI_COLOR_BRIGHT_RED "\nFAILS: %16ld" ANSI_COLOR_RESET "\n", umael(*X0, 3, 1));
    1662             : 
    1663          28 :   if (X0) *X0 = gcopy(mkvec3(Tv, Cv, stoi(umael(*X0, 3, 1))));
    1664             : 
    1665          28 :   gerepileall(av, 1, &answer);
    1666          28 :   return answer;
    1667             : }
    1668             : 
    1669             : GEN
    1670          56 : ecpp(GEN N)
    1671             : {
    1672             :   long expiN, tunelen;
    1673             :   GEN param, answer, garbage, tune;
    1674          56 :   if (typ(N) != t_INT) pari_err_TYPE("ecpp", N);
    1675             : 
    1676             :   /* Check if N is pseudoprime to begin with. */
    1677          56 :   if (!ispseudoprime(N, 0)) return gen_0;
    1678             :   /* Check if we should even prove it. */
    1679          49 :   if (expi(N) < 64) return N;
    1680             : 
    1681          28 :   expiN = expi(N);
    1682          28 :   param = NULL;
    1683          28 :   answer = NULL;
    1684          28 :   garbage = NULL;
    1685             : 
    1686             :   /* tuning for 1500, 2500, 3500, 4500 bits */
    1687             :   /* ecpp is supposed to be faster than isprime on average if N is more than 1500 bits */
    1688          28 :   if (expiN <= 1500) tunelen = 1;
    1689           7 :   else if (expiN <= 2500) tunelen = 2;
    1690           0 :   else if (expiN <= 3500) tunelen = 3;
    1691           0 :   else tunelen = 4;
    1692             : 
    1693          28 :   tune = cgetg(tunelen+1, t_VEC);
    1694          28 :   if (tunelen >= 1) gel(tune, 1) = mkvecsmall3( 500, 24, 22);
    1695          28 :   if (tunelen >= 2) gel(tune, 2) = mkvecsmall3(1500, 32, 23);
    1696          28 :   if (tunelen >= 3) gel(tune, 3) = mkvecsmall3(1500, 32, 24);
    1697          28 :   if (tunelen >= 4) gel(tune, 4) = mkvecsmall3(3000, 40, 24);
    1698             : 
    1699          63 :   while (answer == NULL)
    1700             :   {
    1701             :     pari_timer T;
    1702          35 :     dbg_mode() timer_start(&T);
    1703          35 :     param = ecpp_param_set( tune );
    1704          35 :     dbg_mode() err_printf(ANSI_COLOR_BRIGHT_WHITE "\n%Ps" ANSI_COLOR_RESET, gel(tune, tunelen));
    1705          35 :     dbg_mode() err_printf(ANSI_COLOR_WHITE "  %8ld" ANSI_COLOR_RESET, timer_delay(&T));
    1706          35 :     answer = ecpp0(N, param, &garbage);
    1707          35 :     if (answer != NULL) break;
    1708           7 :     umael(tune, tunelen, 1) *= 2;
    1709           7 :     umael(tune, tunelen, 2) *= 2;
    1710           7 :     umael(tune, tunelen, 3)++;
    1711           7 :     if (umael(tune, tunelen, 3) > 26) umael(tune, tunelen, 3) = 26;
    1712             :   }
    1713          28 :   return answer;
    1714             : }
    1715             : 
    1716             : /*  Input: PARI ECPP Certificate
    1717             :    Output: Human-readable format.
    1718             : */
    1719             : static GEN
    1720           7 : cert_out(GEN x)
    1721             : {
    1722             :   long i;
    1723           7 :   long lgx = lg(x);
    1724             :   pari_str str;
    1725           7 :   str_init(&str, 1);
    1726           7 :   if (typ(x) == t_INT)
    1727             :   {
    1728           0 :     str_printf(&str, "%Ps is prime.\nIndeed, ispseudoprime(%Ps) = 1 and %Ps < 2^64.\n", x, x, x);
    1729           0 :     return strtoGENstr(str.string);
    1730             :   }
    1731          21 :   for (i = 1; i < lgx; i++)
    1732             :   {
    1733          14 :     GEN xi = gel(x, i);
    1734          14 :     str_printf(&str, "\n[%ld]\n", i);
    1735          14 :     str_printf(&str, " N = %Ps\n", cert_get_N(xi));
    1736          14 :     str_printf(&str, " t = %Ps\n", cert_get_t(xi));
    1737          14 :     str_printf(&str, " s = %Ps\n", cert_get_s(xi));
    1738          14 :     str_printf(&str, "a4 = %Ps\n", cert_get_a4(xi));
    1739          14 :     str_printf(&str, "D = %Ps\n", cert_get_D(xi));
    1740          14 :     str_printf(&str, "m = %Ps\n", cert_get_m(xi));
    1741          14 :     str_printf(&str, "q = %Ps\n", cert_get_q(xi));
    1742          14 :     str_printf(&str, "E = %Ps\n", cert_get_E(xi));
    1743          14 :     str_printf(&str, "P = %Ps\n", cert_get_P(xi));
    1744             :   }
    1745           7 :   return strtoGENstr(str.string);
    1746             : }
    1747             : 
    1748             : /*  Input: PARI ECPP Certificate
    1749             :    Output: Magma Certificate
    1750             :      Magma certificate looks like this
    1751             :      (newlines and extraneous spaces for clarity)
    1752             :      [*
    1753             :        [*
    1754             :          N,
    1755             :          |D|,
    1756             :          -1,
    1757             :          m,
    1758             :          [* a, b *],
    1759             :          [* x, y, 1 *],
    1760             :          [*
    1761             :            [* q, 1 *]
    1762             :          *]
    1763             :        *], ...
    1764             :      *]
    1765             : */
    1766             : static GEN
    1767           7 : magma_out(GEN x)
    1768             : {
    1769             :   long i;
    1770           7 :   long size = lg(x);
    1771             :   pari_str ret;
    1772           7 :   str_init(&ret, 1);
    1773           7 :   if (typ(x) == t_INT)
    1774             :   {
    1775           0 :     str_printf(&ret, "Operation not supported.");
    1776           0 :     return strtoGENstr(ret.string);
    1777             :   }
    1778           7 :   str_printf(&ret, "[* ");
    1779          21 :   for (i = 1; i<size; i++)
    1780             :   {
    1781          14 :     GEN xi = gel(x, i);
    1782          14 :     GEN E = cert_get_E(xi);
    1783          14 :     GEN P = cert_get_P(xi);
    1784          14 :     str_printf(&ret, "[* ");
    1785          14 :     str_printf(&ret, "%Ps, ", cert_get_N(xi));
    1786          14 :     str_printf(&ret, "%Ps, ", negi(cert_get_D(xi)));
    1787          14 :     str_printf(&ret, "-1, ");
    1788          14 :     str_printf(&ret, "%Ps, ", cert_get_m(xi));
    1789          14 :     str_printf(&ret, "[* %Ps, %Ps *], ", gel(E, 1), gel(E, 2));
    1790          14 :     str_printf(&ret, "[* %Ps, %Ps, 1 *], ", gel(P, 1), gel(P, 2));
    1791          14 :     str_printf(&ret, "[* [* %Ps, 1 *] *]", cert_get_q(xi));
    1792          14 :     str_printf(&ret, " *]");
    1793          14 :     if (i != size-1) str_printf(&ret, ", ");
    1794             :   }
    1795           7 :   str_printf(&ret, " *]");
    1796           7 :   return strtoGENstr(ret.string);
    1797             : }
    1798             : 
    1799             : /*  Input: &ret, label, value
    1800             :    Prints: label=(sign)hexvalue\n
    1801             :      where sign is + or -
    1802             :      hexvalue is value in hex, of the form: 0xABC123
    1803             : */
    1804             : static void
    1805          70 : primo_printval(pari_str *ret, const char* label, GEN value)
    1806             : {
    1807          70 :   str_printf(ret, label);
    1808          70 :   if (signe(value) >= 0) str_printf(ret, "=0x%P0X\n", value);
    1809          14 :   else str_printf(ret, "=-0x%P0X\n", negi(value));
    1810          70 : }
    1811             : 
    1812             : /*  Input: &ret, label, value
    1813             :    Prints: label=(sign)hexvalue\n
    1814             :      where sign is + or -
    1815             :      hexvalue is x in hex, of the form: 0xABC123
    1816             :        where |x| < N/2 and x = value mod N.
    1817             : */
    1818             : static void
    1819          14 : primo_printval_center(pari_str *ret, const char* label, GEN value, GEN N, GEN shif)
    1820             : {
    1821          14 :   primo_printval(ret, label, Fp_center(value, N, shif) );
    1822          14 : }
    1823             : 
    1824             : /*  Input: PARI ECPP Certificate
    1825             :    Output: PRIMO Certificate
    1826             : 
    1827             :    Let Y^2 = X^3 + aX + b be the elliptic curve over N
    1828             :    with the point (x, y) be the ones described in the
    1829             :    PARI certificate.
    1830             : 
    1831             :    For the case where J != 0, 1728, PRIMO asks for [J,T].
    1832             :    We obtain J easily using the well-known formula.
    1833             :    we obtain T by using the formula T = a/A * B/b * x
    1834             :    where A = 3J(1728-J) and B = 2J(1728-J)^2.
    1835             : 
    1836             :    For the case where J=0 or J=1728, PRIMO asks for [A,B,T].
    1837             :    We let A and B be a and b, respectively.
    1838             :    Consequently, T = x.
    1839             : */
    1840             : static GEN
    1841           7 : primo_out(GEN x)
    1842             : {
    1843             :   long i;
    1844           7 :   long size = (typ(x) == t_INT) ? 1 : lg(x);
    1845             :   pari_str ret;
    1846           7 :   str_init(&ret, 1);
    1847           7 :   str_printf(&ret, "[PRIMO - Primality Certificate]\n");
    1848           7 :   str_printf(&ret, "Format=4\n");
    1849           7 :   str_printf(&ret, "TestCount=%ld\n", size-1);
    1850           7 :   str_printf(&ret, "\n");
    1851           7 :   str_printf(&ret, "[Comments]\n");
    1852           7 :   str_printf(&ret, "Generated by PARI/GP\n");
    1853           7 :   str_printf(&ret, "http://pari.math.u-bordeaux.fr/\n");
    1854           7 :   str_printf(&ret, "\n");
    1855           7 :   str_printf(&ret, "[Candidate]\n");
    1856           7 :   if (typ(x) == t_INT)
    1857             :   {
    1858           0 :     str_printf(&ret, "N=0x%P0X\n", x);
    1859           0 :     return strtoGENstr(ret.string);
    1860             :   } else
    1861           7 :     str_printf(&ret, "N=0x%P0X\n", cert_get_N(gel(x, 1)));
    1862           7 :   str_printf(&ret, "\n");
    1863             : 
    1864          21 :   for (i = 1; i<size; i++)
    1865             :   {
    1866             :     GEN xi, N, Nover2;
    1867             :     long Ais0, Bis0;
    1868          14 :     str_printf(&ret, "[%ld]\n", i);
    1869          14 :     xi = gel(x, i);
    1870          14 :     N = cert_get_N(xi);
    1871          14 :     Nover2 = shifti(N, -1);
    1872          14 :     primo_printval(&ret, "S", cert_get_s(xi));
    1873          14 :     primo_printval(&ret, "W", cert_get_t(xi));
    1874          14 :     Ais0 = isintzero(cert_get_a4(xi));
    1875          14 :     Bis0 = isintzero(cert_get_a6(xi));
    1876          14 :     if (!Ais0 && !Bis0) { /* J != 0, 1728 */
    1877           0 :       primo_printval_center(&ret, "J", cert_get_J(xi), N, Nover2);
    1878           0 :       primo_printval(&ret, "T", cert_get_T(xi));
    1879          14 :     } else if (Ais0) { /* J == 0 */
    1880          14 :       primo_printval(&ret, "A", gen_0);
    1881          14 :       primo_printval_center(&ret, "B", cert_get_a6(xi), N, Nover2);
    1882          14 :       primo_printval(&ret, "T", cert_get_x(xi));
    1883             :     }else{ /* J == 1728 */
    1884           0 :       primo_printval_center(&ret, "A", cert_get_a4(xi), N, Nover2);
    1885           0 :       primo_printval(&ret, "B", gen_0);
    1886           0 :       primo_printval(&ret, "T", cert_get_x(xi));
    1887             :     }
    1888          14 :     str_printf(&ret, "\n");
    1889             :   }
    1890           7 :   return strtoGENstr(ret.string);
    1891             : }
    1892             : 
    1893             : /*  Input: N, q
    1894             :    Output: 1 if q > (N^1/4 + 1)^2, 0 otherwise.
    1895             : */
    1896             : static long
    1897         931 : Nq_isvalid(GEN N, GEN q)
    1898             : {
    1899         931 :   GEN qm1sqr = sqri(subiu(q, 1));
    1900         931 :   if (cmpii(qm1sqr,N) > 0) /* (q-1)^2 > N */
    1901             :   { /* (q-1)^4 + N^2 > 16Nq + 2N(q-1)^2 */
    1902         931 :     GEN a = addii(sqri(qm1sqr), sqri(N));
    1903         931 :     GEN b = addii(shifti(mulii(N, q), 4 ), mulii(shifti(N, 1), qm1sqr));
    1904         931 :     return (cmpii(a, b) > 0);
    1905             :   }
    1906           0 :   return 0;
    1907             : }
    1908             : 
    1909             : /*  Input: cert, trustbits=64
    1910             :    Output: 1 if cert is a valid PARI ECPP certificate
    1911             :      assuming that all q below 2^trustbits is prime
    1912             : */
    1913             : static long
    1914          56 : ecppisvalid0(GEN cert, long trustbits)
    1915             : {
    1916          56 :   pari_sp av = avma;
    1917             :   long i;
    1918          56 :   long lgcert = lg(cert);
    1919             :   GEN ql;
    1920          56 :   GEN N, W, m, s, q = gen_0, P, a;
    1921             :   GEN mP, sP;
    1922             :   GEN r;
    1923             : 
    1924          56 :   if (typ(cert) == t_INT)
    1925             :   {
    1926           7 :     if (expi(cert) > trustbits) return 0;
    1927           7 :     if (cmpii(cert, shifti(gen_1, 64)) > 0) return 0;
    1928           7 :     return ispseudoprime(cert, 0);
    1929             :   }
    1930             : 
    1931          49 :   if (typ(cert) != t_VEC) pari_err_TYPE("ecppisvalid", cert);
    1932             : 
    1933          49 :   if (lgcert <= 1) return 0;
    1934          49 :   if (lg(gel(cert, lgcert-1))-1 != 5) return 0;
    1935          49 :   ql = cert_get_q(gel(cert, lgcert-1));
    1936          49 :   if (expi(ql) > trustbits) return 0;
    1937          49 :   if (!ispseudoprime(ql, 0)) return 0;
    1938             : 
    1939         980 :   for (i = 1; i < lgcert; i++)
    1940             :   {
    1941         931 :     GEN certi = gel(cert, i);
    1942         931 :     if (lg(certi)-1 != 5) return 0;
    1943             : 
    1944         931 :     N = cert_get_N(certi);
    1945         931 :     if (typ(N) != t_INT) return 0;
    1946             :     /* Check if N > 1 */
    1947         931 :     if (signe(N) != 1) return 0;
    1948             :     /* Check if N is not divisible by 2 or 3 */
    1949         931 :     if (!isint1(gcdii(N, stoi(6)))) return 0;
    1950             :     /* Check if N matches the q from the previous entry. */
    1951         931 :     if (i > 1 && !equalii(N, q)) return 0;
    1952             : 
    1953         931 :     W = cert_get_t(certi);
    1954         931 :     if (typ(W) != t_INT) return 0;
    1955             :     /* Check if W^2 < 4N (Hasse interval) */
    1956         931 :     if (!(cmpii(sqri(W), shifti(N,2)) < 0)) return 0;
    1957             : 
    1958         931 :     s = cert_get_s(certi);
    1959         931 :     if (typ(s) != t_INT) return 0;
    1960         931 :     if (signe(s) < 1) return 0;
    1961             : 
    1962         931 :     m = cert_get_m(certi);
    1963         931 :     q = dvmdii(m, s, &r);
    1964             : 
    1965             :     /* Check m%s == 0 */
    1966         931 :     if (!isintzero(r)) return 0;
    1967             : 
    1968             :     /* Check q > (N^(1/4) + 1)^2 */
    1969         931 :     if (!Nq_isvalid(N, q)) return 0;
    1970             : 
    1971         931 :     a = cert_get_a4(certi);
    1972         931 :     if (typ(a) != t_INT) return 0;
    1973             : 
    1974         931 :     P = cert_get_P(certi);
    1975         931 :     if (typ(P) != t_VEC) return 0;
    1976         931 :     if (lg(P)-1 != 2) return 0;
    1977         931 :     P = FpE_to_FpJ(P);
    1978             : 
    1979             :     /* Check mP == 0 */
    1980         931 :     mP = FpJ_mul(P, m, a, N);
    1981         931 :     if (!FpJ_is_inf(mP)) return 0;
    1982             : 
    1983             :     /* Check sP != 0 and third component is coprime to N */
    1984         931 :     sP = FpJ_mul(P, s, a, N);
    1985         931 :     if (!isint1(gcdii(gel(sP, 3), N))) return 0;
    1986             :   }
    1987          49 :   avma = av; return 1;
    1988             : }
    1989             : 
    1990             : long
    1991          56 : ecppisvalid(GEN cert)
    1992          56 : { return ecppisvalid0(cert, 64); }
    1993             : 
    1994             : GEN
    1995          21 : ecppexport(GEN cert, long flag)
    1996             : {
    1997          21 :   if (!ecppisvalid(cert))
    1998           0 :     pari_err_TYPE("ecppexport - certificate not valid", cert);
    1999          21 :   switch(flag)
    2000             :   {
    2001           7 :     case 1: return magma_out(cert);
    2002           7 :     case 2: return primo_out(cert);
    2003           7 :     default: return cert_out(cert);
    2004             :   }
    2005             : }

Generated by: LCOV version 1.11