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

Generated by: LCOV version 1.11