Code coverage tests

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

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

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

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

Generated by: LCOV version 1.13