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.16.2 lcov report (development 29115-f22e516b23) Lines: 748 805 92.9 %
Date: 2024-03-28 08:06:56 Functions: 95 98 96.9 %
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        5350 : 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    57654507 : D_get_wD(long D)
      82             : {
      83    57654507 :   if (D == -4) return 4;
      84    57652946 :   if (D == -3) return 6;
      85    57651133 :   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        2553 : 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       15521 : 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       14345 : 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        1057 : NDmqg_get_N(GEN x) { return gel(x,1); }
     192             : INLINE GEN
     193        9786 : 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        1057 : NDmqg_get_g(GEN x) { return gel(x,5); }
     200             : INLINE GEN
     201        1057 : 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(veclast(x)); }
     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             : /* "Twist factor". Does not cover J = 0, 1728 */
     520             : static GEN
     521          49 : cert_get_lambda(GEN z)
     522             : {
     523             :   GEN N, J, a, b, A, B;
     524          49 :   J = cert_get_J(z);
     525          49 :   N = cert_get_N(z);
     526          49 :   a = cert_get_a4(z);
     527          49 :   b = cert_get_a6(z);
     528          49 :   Fp_ellj_to_a4a6(J, N, &A, &B);
     529          49 :   return Fp_div(Fp_mul(a,B,N), Fp_mul(b,A,N), N);
     530             : }
     531             : 
     532             : /* Solves for T such that if
     533             :    [A, B] = [3*J*(1728-J), 2*J*(1728-J)^2]
     534             :    and you let
     535             :    L = T^3 + A*T + B, A = A*L^2, B = B*L^3
     536             :    then
     537             :    x == TL and y == L^2
     538             : */
     539             : static GEN
     540          49 : cert_get_T(GEN z)
     541             : {
     542          49 :   GEN N = cert_get_N(z), x = cert_get_x(z);
     543          49 :   GEN l = cert_get_lambda(z); /* l = 1/L */
     544          49 :   return Fp_mul(x, l, N);
     545             : }
     546             : 
     547             : /* database for all invariants, stolen from Hamish's code */
     548             : static GEN
     549          56 : polmodular_db_init_allinv(void)
     550             : {
     551          56 :   const long DEFAULT_MODPOL_DB_LEN = 32;
     552          56 :   GEN a, b = cgetg(39+1, t_VEC);
     553             :   long i;
     554        2240 :   for (i = 1; i < 40; i++) gel(b,i) = zerovec_block(DEFAULT_MODPOL_DB_LEN);
     555          56 :   a = zerovec_block(DEFAULT_MODPOL_DB_LEN);
     556          56 :   return mkvec2(a, b);
     557             : }
     558             : 
     559             : /* Given D and a database of modular polynomials, return polclass(D, inv) */
     560             : static GEN
     561        1057 : D_polclass(long D, long inv, GEN *db)
     562             : {
     563        1057 :   GEN HD, t = mkvec2(gel(*db, 1), inv == 0? gen_0: gmael(*db, 2, inv));
     564        1057 :   HD = polclass0(D, inv, 0, &t);
     565        1057 :   gel(*db, 1) = gel(t,1);
     566        1057 :   if (inv != 0) gmael(*db, 2, inv) = gel(t,2);
     567        1057 :   return HD;
     568             : }
     569             : 
     570             : static GEN
     571        2423 : NqE_check(GEN N, GEN q, GEN a, GEN b, GEN s)
     572             : {
     573        2423 :   GEN mP, sP, P = random_FpE(a, b, N);
     574        2423 :   GEN PJ = mkvec3(gel(P,1), gel(P,2), gen_1);
     575        2423 :   sP = FpJ_mul(PJ, s, a, N); if (FpJ_is_inf(sP)) return NULL;
     576        2423 :   mP = FpJ_mul(sP, q, a, N); return FpJ_is_inf(mP)? mkvec2(a, P): NULL;
     577             : }
     578             : 
     579             : /* Find an elliptic curve E modulo N and a point P on E which corresponds to
     580             :  * the ith element of the certificate; uses N, D, m, q, J from previous steps.
     581             :  * All computations are modulo N unless stated otherwise */
     582             : 
     583             : /* g is a quadratic and cubic nonresidue modulo N */
     584             : static GEN
     585         203 : j0_find_g(GEN N)
     586             : {
     587         203 :   GEN n = diviuexact(subiu(N, 1), 3);
     588             :   for(;;)
     589         490 :   {
     590         693 :     GEN g = randomi(N);
     591         693 :     if (kronecker(g, N) == -1 && !equali1(Fp_pow(g, n, N))) return g;
     592             :   }
     593             : }
     594             : 
     595             : static GEN
     596        1057 : find_EP(GEN N, long D, GEN q, GEN g, GEN J, GEN s)
     597             : {
     598             :   GEN A0, B0;
     599        1057 :   Fp_ellj_to_a4a6(J, N, &A0, &B0);
     600             :   for(;;)
     601           0 :   { /* expect one iteration: not worth saving the A's and B's */
     602        1057 :     GEN gg, v, A = A0, B = B0;
     603             :     long i;
     604        1057 :     if ((v = NqE_check(N, q, A, B, s))) return v;
     605         687 :     switch (D_get_wD(D))
     606             :     {
     607         302 :       case 2:
     608         302 :         gg = Fp_sqr(g, N);
     609         302 :         A = Fp_mul(A, gg, N);
     610         302 :         B = Fp_mul(Fp_mul(B, gg, N), g, N);
     611         302 :         if ((v = NqE_check(N, q, A, B, s))) return v;
     612           0 :         break;
     613         140 :       case 4:
     614         294 :         for (i = 1; i < 4; i++)
     615             :         {
     616         294 :           A = Fp_mul(A, g, N);
     617         294 :           if ((v = NqE_check(N, q, A, B, s))) return v;
     618             :         }
     619           0 :         break;
     620         245 :       default: /* 6 */
     621         245 :         B = Fp_mul(B, g, N);
     622         245 :         if ((v = NqE_check(N, q, A, B, s))) return v;
     623         203 :         g = j0_find_g(N);
     624         630 :         for (i = 1; i < 6; i++)
     625             :         {
     626         630 :           B = Fp_mul(B, g, N);
     627         630 :           if (i == 3) continue;
     628         525 :           if ((v = NqE_check(N, q, A, B, s))) return v;
     629             :         }
     630           0 :         break;
     631             :     }
     632             :   }
     633             : }
     634             : 
     635             : /* Return the genus field. If real is set only the real part */
     636             : static GEN
     637         413 : genusfield(GEN Dfac, GEN sq, GEN p, long real)
     638             : {
     639         413 :   long i, j, l = lg(Dfac), dn, n = 0;
     640         413 :   GEN sn, s = gen_0, R = cgetg(l-1+!real, t_VECSMALL);
     641         476 :   for (i = 1; i < l; i++)
     642         476 :     if (Dfac[i] < 0) { n = i; break; }
     643         413 :   if (n == 0) pari_err_BUG("realgenusfield");
     644         413 :   dn = Dfac[n]; sn = gel(sq, n);
     645        1267 :   for (j = i = 1; i < l; i++)
     646         854 :     if (i != n)
     647             :     {
     648         441 :       long d = Dfac[i];
     649         441 :       GEN si = gel(sq,i);
     650         441 :       R[j++] = d > 0 ? d : d * dn;
     651         441 :       s = Fp_add(s, d > 0? si: Fp_mul(si, sn, p), p);
     652             :     }
     653         413 :   if (!real)
     654             :   {
     655         161 :     R[j++] = dn;
     656         161 :     s = Fp_add(s,sn, p);
     657             :   }
     658         413 :   return mkvec2(R, s);
     659             : }
     660             : 
     661             : static GEN
     662          35 : realpart(GEN P, GEN R)
     663             : {
     664          35 :   GEN G = galoisinit(rnfequation(P,R), NULL);
     665          35 :   GEN T = galoisfixedfield(G, gel(gal_get_gen(G),2), 1, -1);
     666          35 :   setvarn(T,1); return polredbest(T,0);
     667             : }
     668             : 
     669             : static GEN
     670          35 : dihedralpol(long D, long l)
     671             : {
     672          35 :   GEN P = quadpoly0(utoineg(-D), 1);
     673          35 :   GEN bnf = Buchall(P, 1, BIGDEFAULTPREC);
     674          35 :   GEN x = bnrclassfield(bnf, utoipos(l), 0, BIGDEFAULTPREC);
     675          70 :   pari_APPLY_same(realpart(P, gel(x,i)))
     676             : }
     677             : 
     678             : static GEN
     679        1057 : FpX_classtower_oneroot(GEN P, GEN Dinfo, GEN Dfac, GEN sq, GEN p)
     680             : {
     681        1057 :   long D = Dinfo_get_D(Dinfo), h = Dinfo_get_h(Dinfo);
     682        1057 :   if (degpol(P) > 1)
     683             :   {
     684         413 :     long real = !modinv_is_double_eta(Dinfo_get_bi(Dinfo));
     685         413 :     GEN V = genusfield(Dfac, sq, p, real), v = gel(V,1), R = gel(V,2);
     686         413 :     GEN N = NULL;
     687         413 :     long i, l = lg(v);
     688        1015 :     for (i = 1; i < l; i++)
     689             :     {
     690         602 :       GEN Q = deg2pol_shallow(gen_1, gen_0, stoi(-v[i]), 1);
     691         602 :       if (!N) N = Q;
     692             :       else
     693             :       {
     694         231 :         GEN cm = polcompositum0(N,Q,3);
     695         231 :         N = gel(cm,1); P = QXY_QXQ_evalx(P,gmael(cm,2,2),N);
     696             :       }
     697         602 :       P = liftpol_shallow(gmael(nffactor(N,P),1,1));
     698             :     }
     699         413 :     if (h%3 == 0 && (!N || degpol(N)*3<degpol(P)))
     700             :     {
     701          35 :       GEN v = dihedralpol(D,3);
     702          35 :       long i, l = lg(v);
     703          35 :       for (i = 1; i < l; i++)
     704             :       {
     705          35 :         GEN Q = gel(v,i);
     706          35 :         if (!N) N = Q;
     707             :         else
     708             :         {
     709           0 :           GEN cm = polcompositum0(N, Q, 3);
     710           0 :           N = gel(cm,1); P = QXY_QXQ_evalx(P,gmael(cm,2,2),N);
     711           0 :           R = Fp_mul(R, gel(cm,4),p);
     712             :         }
     713          35 :         R = Fp_add(R, FpX_oneroot_split(Q, p), p);
     714          35 :         P = liftpol_shallow(gmael(nffactor(N,P),1,1));
     715          35 :         if (degpol(N)*3>=degpol(P)) break;
     716             :       }
     717             :     }
     718         413 :     if (N)
     719         406 :       P = FpXY_evalx(Q_primpart(P), R, p);
     720             :   }
     721        1057 :   return P;
     722             : }
     723             : 
     724             : GEN
     725        1057 : ecpp_step2_worker(GEN S, GEN HD, GEN primelist, long dbg)
     726             : {
     727        1057 :   pari_sp av = avma;
     728             :   pari_timer ti;
     729             :   GEN J, t, s, EP, rt, res;
     730        1057 :   GEN N = NDmqg_get_N(S), Dinfo = NDmqg_get_Dinfo(S);
     731        1057 :   GEN m = NDmqg_get_m(S), q = NDmqg_get_q(S);
     732        1057 :   GEN g = NDmqg_get_g(S), sq = NDmqg_get_sqrt(S);
     733        1056 :   long D = Dinfo_get_D(Dinfo), inv = Dinfo_get_bi(Dinfo);
     734        1056 :   GEN Dfacp = Dfac_to_p(Dinfo_get_Dfac(Dinfo), primelist);
     735        1056 :   long C2 = 0, C3 = 0, C4 = 0, D1 = 0;
     736        1056 :   GEN P, c = getrand();
     737        1057 :   setrand(gen_1); /* for reproducibility */
     738             :   /* C2: Find a root modulo N of polclass(D,inv) */
     739        1057 :   if (dbg >= 2) timer_start(&ti);
     740        1057 :   P = FpX_classtower_oneroot(HD, Dinfo, Dfacp, sq, N);
     741        1057 :   if (dbg >= 2) C2 = timer_delay(&ti);
     742        1057 :   rt = FpX_oneroot_split(P, N);
     743        1057 :   if (dbg >= 2) C3 = timer_delay(&ti);
     744             :   /* C3: Convert root from previous step into the appropriate j-invariant */
     745        1057 :   J = Fp_modinv_to_j(rt, inv, N); /* root of polclass(D) */
     746        1057 :   if (dbg >= 2) C4 = timer_delay(&ti);
     747             :   /* D1: Find an elliptic curve E with a point P satisfying the theorem */
     748        1057 :   s = diviiexact(m, q);
     749        1057 :   EP = find_EP(N, D, q, g, J, s);
     750        1057 :   if (dbg >= 2) D1 = timer_delay(&ti);
     751             : 
     752             :   /* D2: Compute for t and s */
     753        1057 :   t = subii(addiu(N, 1), m); /* t = N+1-m */
     754        1057 :   res = mkvec2(mkvec5(N, t, s, gel(EP,1), gel(EP,2)),mkvecsmall4(C2,C3,C4,D1));
     755        1057 :   setrand(c);
     756        1057 :   return gerepilecopy(av, res);
     757             : }
     758             : 
     759             : /* This uses [N, D, m, q] from step 1 to find the appropriate j-invariants
     760             :  * to be used in step 3. Step 2 is divided into substeps 2a, 2b, 2c */
     761             : static GEN
     762          56 : ecpp_step2(GEN step1, GEN *X0, GEN primelist)
     763             : {
     764          56 :   long j, Dprev = 0;
     765             :   pari_timer ti;
     766          56 :   GEN perm = gen_indexsort(step1, NULL, &sort_NDmq_by_D);
     767          56 :   GEN step2 = cgetg(lg(step1), t_VEC);
     768          56 :   GEN HD = NULL, db = polmodular_db_init_allinv();
     769          56 :   long ls = lg(step2), pending = 0;
     770             :   GEN worker;
     771             :   struct pari_mt pt;
     772          56 :   GEN Hi = cgetg(ls, t_VEC);
     773        1113 :   for (j = 1; j < ls; j++)
     774             :   {
     775        1057 :     long i = uel(perm, j);
     776        1057 :     GEN S = gel(step1, i);
     777        1057 :     GEN Dinfo = NDmqg_get_Dinfo(S);
     778        1057 :     long D = Dinfo_get_D(Dinfo), inv = Dinfo_get_bi(Dinfo);
     779             :     /* C1: Find the appropriate class polynomial modulo N */
     780        1057 :     dbg_mode() timer_start(&ti);
     781        1057 :     if (D != Dprev) HD = D_polclass(D, inv, &db);
     782        1057 :     dbg_mode() {
     783           0 :       GEN N = NDmqg_get_N(S);
     784           0 :       long tt = timer_record(X0, "C1", &ti);
     785           0 :       err_printf(ANSI_BRIGHT_GREEN "\n[ %3d | %4ld bits]" ANSI_RESET, i, expi(N));
     786           0 :       err_printf(ANSI_GREEN " D = %8ld poldeg = %4ld" ANSI_RESET, D, degpol(HD));
     787           0 :       if (D == Dprev) err_printf(" %6ld", tt);
     788           0 :       else err_printf(ANSI_BRIGHT_WHITE " %6ld" ANSI_RESET, tt);
     789             :     }
     790        1057 :     gel(Hi, i) = HD;
     791             :   }
     792          56 :   gunclone_deep(db);
     793          56 :   worker = snm_closure(is_entry("_ecpp_step2_worker"),
     794             :                        mkvec2(primelist,stoi(DEBUGLEVEL)));
     795          56 :   mt_queue_start_lim(&pt, worker, ls-1);
     796        1182 :   for (j = 1; j < ls || pending; j++)
     797             :   {
     798             :     long jj;
     799        1126 :     GEN S = gel(step1, j), HD = gel(Hi, j), done;
     800        1126 :     mt_queue_submit(&pt, j, j < ls ? mkvec2(S, HD) : NULL);
     801        1126 :     done = mt_queue_get(&pt, &jj, &pending);
     802        1126 :     if (!done) continue;
     803        1057 :     dbg_mode() {
     804           0 :       GEN T = gel(done,2);
     805           0 :       GEN N = NDmqg_get_N(gel(step1, jj));
     806           0 :       err_printf(ANSI_BRIGHT_GREEN "\n[ %3d | %4ld bits]" ANSI_RESET, jj, expi(N));
     807           0 :       err_printf(" %6ld", time_record(X0, "C2", T[1]));
     808           0 :       err_printf(" %6ld", time_record(X0, "C3", T[2]));
     809           0 :       err_printf(" %6ld", time_record(X0, "C4", T[3]));
     810           0 :       err_printf(" %6ld", time_record(X0, "D1", T[4]));
     811             :     }
     812        1057 :     gel(step2, jj) = gel(done,1);
     813             :   }
     814          56 :   mt_queue_end(&pt);
     815          56 :   return step2;
     816             : }
     817             : /* end of functions for step 2 */
     818             : 
     819             : GEN
     820       17997 : ecpp_sqrt_worker(GEN p, GEN g, GEN N)
     821             : {
     822       17997 :   GEN r = Fp_sqrt_i(p, g, N);
     823       18002 :   return r ? r: gen_0;
     824             : }
     825             : 
     826             : static void
     827       14345 : ecpp_parsqrt(GEN N, GEN param, GEN kroP, GEN sqrtlist, GEN g, long nb, long *nbsqrt)
     828             : {
     829       14345 :   GEN primelist = ecpp_param_get_primelist(param);
     830       14345 :   GEN ordinv = ecpp_param_get_primeord(param);
     831       14345 :   long i, k, l = lg(ordinv), n = *nbsqrt+1;
     832       14345 :   GEN worker = snm_closure(is_entry("_ecpp_sqrt_worker"), mkvec2(g, N));
     833       14345 :   GEN W = cgetg(nb+1,t_VEC), V;
     834       50301 :   for (i=n, k=1; k<=nb && i<l; i++)
     835             :   {
     836       35956 :     long pi = ordinv[i];
     837       35956 :     if (kroP[pi] > 0)
     838       18003 :       gel(W,k++) = stoi(primelist[pi]);
     839             :   }
     840       14345 :   *nbsqrt=i-1;
     841       14345 :   setlg(W,k);
     842       14345 :   V = gen_parapply(worker, W);
     843       50301 :   for (i=n, k=1; k<=nb && i<l; i++)
     844             :   {
     845       35956 :     long pi = ordinv[i];
     846       35956 :     if (kroP[pi] > 0)
     847             :     {
     848       18003 :       dbg_mode() err_printf(ANSI_MAGENTA "S" ANSI_RESET);
     849       18003 :       if (!signe(gel(V,k)))
     850           0 :         pari_err_BUG("D_find_discsqrt"); /* possible if N composite */
     851       18003 :       gel(sqrtlist,pi) = gel(V,k++);
     852             :     }
     853             :   }
     854       14345 : }
     855             : 
     856             : /* start of functions for step 1 */
     857             : 
     858             : /* This finds the square root of D modulo N [given by Dfac]: find the square
     859             :  * root modulo N of each prime p dividing D and multiply them out */
     860             : static GEN
     861       40677 : D_find_discsqrt(GEN N, GEN param, GEN Dfac, GEN kroP, GEN sqrtlist, GEN g, long *nbsqrt)
     862             : {
     863       40677 :   GEN s = NULL, sj;
     864       40677 :   long i, l = lg(Dfac);
     865      123914 :   for (i = 1; i < l; i++)
     866             :   {
     867       83237 :     long j = Dfac[i];
     868       97582 :     while(!signe(gel(sqrtlist,j)))
     869       14345 :       ecpp_parsqrt(N, param, kroP, sqrtlist, g, mt_nbthreads(), nbsqrt);
     870       83237 :     sj = gel(sqrtlist,j);
     871       83237 :     s = s? Fp_mul(s, sj, N): sj;
     872             :   }
     873       40677 :   return s;/* != NULL */
     874             : }
     875             : 
     876             : /* Given a solution U, V to U^2 + |D|V^2 = 4N, this find all the possible
     877             :      cardinalities of the elliptic curves over the finite field F_N
     878             :      whose endomorphism ring is isomorphic to the maximal order
     879             :      in the imaginary quadratic number field K = Q(sqrt(D)). ???
     880             : */
     881             : static GEN
     882       10864 : NUV_find_m(GEN Dinfo, GEN N, GEN U, GEN V, long wD)
     883             : {
     884       10864 :   GEN m, Nplus1 = addiu(N, 1), u = U, mvec = cgetg(wD+1, t_VEC);
     885             :   long i;
     886       36358 :   for (i = 0; i < wD; i++)
     887             :   {
     888       25494 :     if (odd(i)) m = subii(Nplus1, u);
     889             :     else
     890             :     {
     891       12747 :       if (wD == 4 && i==2) u = shifti(V, 1);
     892       12250 :       else if (wD == 6 && i==2) u = shifti(submuliu(U, V, 3), -1);
     893       11557 :       else if (wD == 6 && i==4) u = shifti(addmuliu(U, V, 3), -1);
     894       12747 :       m = addii(Nplus1, u);
     895             :     }
     896       25494 :     gel(mvec, i+1) = mkvec2(Dinfo, m);
     897             :   }
     898       10864 :   return mvec;
     899             : }
     900             : 
     901             : /* Populates Dmbatch with Dmvec's -- whose components are of the form [D,m],
     902             :    where m is a cardinalities of an elliptic curve with endomorphism ring
     903             :    isomorphic to the maximal order of imaginary quadratic K = Q(sqrt(D)).
     904             :    It returns 0 if:
     905             :      * Any of the p* dividing D is not a quadratic residue mod N
     906             :      * Cornacchia cannot find a solution U^2 + |D|V^2 = 4N.
     907             :    Otherwise, it returns the number of cardinalities added to Dmbatch.
     908             :    Finally, sqrtlist and g help compute the square root modulo N of D.
     909             : */
     910             : static long
     911      538034 : D_collectcards(GEN N, GEN param, GEN* X0, GEN Dinfo, GEN sqrtlist, GEN g, GEN Dmbatch, GEN kroP, long *nbsqrt)
     912             : {
     913      538034 :   long i, l, corn_succ, wD, D = Dinfo_get_D(Dinfo);
     914      538034 :   GEN U, V, sqrtofDmodN, Dfac = Dinfo_get_Dfac(Dinfo);
     915             :   pari_timer ti;
     916             : 
     917             :   /* A1: Check (p*|N) = 1 for all primes dividing D */
     918      538034 :   l = lg(Dfac);
     919      768992 :   for (i = 1; i < l; i++)
     920             :   {
     921      728315 :     long j = Dfac[i];
     922      728315 :     if (kroP[j] < 0) return 0;
     923             :   }
     924             :   /* A3: Get square root of D mod N */
     925       40677 :   dbg_mode() timer_start(&ti);
     926       40677 :   sqrtofDmodN = D_find_discsqrt(N, param, Dfac, kroP, sqrtlist, g, nbsqrt);
     927       40677 :   dbg_mode() timer_record(X0, "A3", &ti);
     928             :   /* A5: Use square root with Cornacchia to solve U^2 + |D|V^2 = 4N */
     929       40677 :   dbg_mode() timer_start(&ti);
     930       40677 :   corn_succ = cornacchia2_sqrt( utoi(labs(D)), N, sqrtofDmodN, &U, &V);
     931       40677 :   dbg_mode() timer_record(X0, "A5", &ti);
     932       40677 :   if (!corn_succ) {
     933       29813 :     dbg_mode1() err_printf(ANSI_YELLOW "c" ANSI_RESET);
     934       29813 :     return 0;
     935             :   }
     936             :   /* A6: Collect the w(D) cardinalities of E/(F_N) with CM by D */
     937       10864 :   dbg_mode() err_printf(ANSI_BRIGHT_YELLOW "D" ANSI_RESET);
     938       10864 :   wD = D_get_wD(D);
     939       10864 :   vectrunc_append_batch(Dmbatch, NUV_find_m(Dinfo,N,U,V,wD));
     940       10864 :   return wD;
     941             : }
     942             : 
     943             : /* Compute (S(16N, 2) + S(4096N, 4) + 4)\4 + 1,  where S is the nth root
     944             :  * rounded down. This is at most one more than (N^1/4 + 1)^2. */
     945             : static GEN
     946        3857 : ecpp_qlo(GEN N)
     947             : {
     948        3857 :   GEN a = sqrtnint(shifti(N, 4), 2);
     949        3857 :   GEN b = sqrtnint(shifti(N, 12), 4);
     950        3857 :   return addiu(shifti(addii(a, b), -2), 2);
     951             : }
     952             : 
     953             : static long
     954       13356 : lessthan_qlo(void* E, GEN Dmq) { return (cmpii(gel(Dmq,3), (GEN)E) < 0); }
     955             : static long
     956       12047 : gained_bits(void* E, GEN Dmq) { return (expi(gel(Dmq,3)) <= (long)E); }
     957             : 
     958             : /*  Input: Dmqvec
     959             :  * Output: Dmqvec such that q satisfies (N^1/4 + 1)^2 < q < N/2 */
     960             : static GEN
     961        3857 : Dmqvec_slice(GEN N, GEN Dmqvec)
     962             : {
     963             :   long lo, hi;
     964             : 
     965        3857 :   gen_sort_inplace(Dmqvec, NULL, &sort_Dmq_by_q, NULL); /* sort wrt q */
     966        3857 :   lo = zv_binsearch0((void*)ecpp_qlo(N), &lessthan_qlo, Dmqvec); lo++;
     967        3857 :   if (lo >= lg(Dmqvec)) return NULL;
     968             : 
     969        3857 :   hi = zv_binsearch0((void*)(expi(N)-1), &gained_bits, Dmqvec);
     970        3857 :   if (hi == 0) return NULL;
     971             : 
     972        3857 :   return vecslice(Dmqvec, lo, hi);
     973             : }
     974             : 
     975             : /* Given a Dmvec of [D,m]'s, simultaneously remove all prime factors of each
     976             :  * m less then BOUND_PRIMORIAL. This implements Franke 2004: Proving the
     977             :  * Primality of Very Large Numbers with fastECPP */
     978             : static GEN
     979        3857 : Dmvec_batchfactor(GEN Dmvec, GEN primorial)
     980             : {
     981        3857 :   long i, l = lg(Dmvec);
     982        3857 :   GEN leaf, v = cgetg(l, t_VEC);
     983       29351 :   for (i = 1; i < l; i++)
     984             :   { /* cheaply remove powers of 2 */
     985       25494 :     GEN m = gmael(Dmvec, i, 2);
     986       25494 :     long v2 = vali(m);
     987       25494 :     gel(v,i) = v2? shifti(m,-v2): m;
     988             :   }
     989        3857 :   leaf = Z_ZV_mod_tree(primorial, v, ZV_producttree(v));
     990             :   /* Go through each leaf and remove small factors. */
     991       29351 :   for (i = 1; i < l; i++)
     992             :   {
     993       25494 :     GEN q = gel(v,i), Dm = gel(Dmvec,i), L = gel(leaf,i);
     994       84315 :     while (!equali1(L)) { L = gcdii(q, L); q = diviiexact(q, L); }
     995       25494 :     gel(v,i) = mkvec3(gel(Dm,1), gel(Dm,2), q);
     996             :   }
     997        3857 :   return v;
     998             : }
     999             : 
    1000             : /* return good parameters [maxsqrt, maxpcdg, tdivexp] for ecpp(N) */
    1001             : static GEN
    1002        4977 : tunevec(long expiN, GEN param)
    1003             : {
    1004        4977 :   GEN T = ecpp_param_get_tune(param);
    1005        4977 :   long i, n = lg(T)-1;
    1006        8106 :   for (i = 1; i < n; i++) { GEN x = gel(T,i); if (expiN <= x[4]) return x; }
    1007        2219 :   return gel(T,n);
    1008             : }
    1009             : static long
    1010        4977 : tunevec_tdivbd(long expiN, GEN param) { return uel(tunevec(expiN, param), 3); }
    1011             : static long
    1012        1120 : tunevec_batchsize(long expiN, GEN param)
    1013             : {
    1014        1120 :   long t, b = 28 - tunevec_tdivbd(expiN, param);
    1015        1120 :   if (b < 0) return expiN;
    1016        1120 :   t = expiN >> b; return t < 1? 1: t;
    1017             : }
    1018             : 
    1019             : static GEN
    1020        3857 : Dmbatch_factor_Dmqvec(GEN N, long expiN, GEN* X0, GEN Dmbatch, GEN param)
    1021             : {
    1022        3857 :   pari_sp av = avma;
    1023             :   pari_timer ti;
    1024        3857 :   GEN Dmqvec, primorial_vec = ecpp_param_get_primorial_vec(param);
    1025        3857 :   GEN primorial = gel(primorial_vec, tunevec_tdivbd(expiN, param)-19);
    1026             : 
    1027             :   /* B1: Factor by batch */
    1028        3857 :   dbg_mode() timer_start(&ti);
    1029        3857 :   Dmqvec = Dmvec_batchfactor(Dmbatch, primorial);
    1030        3857 :   dbg_mode() timer_record(X0, "B1", &ti);
    1031             : 
    1032             :   /* B2: For each batch, remove cardinalities lower than (N^(1/4)+1)^2
    1033             :    *     and cardinalities in which we didn't win enough bits */
    1034        3857 :   Dmqvec = Dmqvec_slice(N, Dmqvec);
    1035        3857 :   if (!Dmqvec) return gc_NULL(av); /* nothing is left */
    1036        3857 :   return gerepilecopy(av, Dmqvec);
    1037             : }
    1038             : 
    1039             : GEN
    1040       20938 : ecpp_ispsp_worker(GEN N)
    1041       20938 : { return mkvecsmall(BPSW_psp_nosmalldiv(N)); }
    1042             : 
    1043             : static GEN
    1044        1057 : mkNDmqg(GEN z, GEN N, GEN Dmq, GEN g, GEN sqrtlist)
    1045             : {
    1046        1057 :   GEN Dinfo = gel(Dmq,1), sq =  Dfac_to_roots(Dinfo_get_Dfac(Dinfo),sqrtlist);
    1047        1057 :   GEN NDmqg = mkcol6(N, Dinfo, gel(Dmq,2), gel(Dmq,3), g, sq);
    1048        1057 :   return mkvec2(NDmqg, z);
    1049             : }
    1050             : /* expi(N) > 64. Return [ NDmqg, N_downrun(q) ]; NDmqg is a vector of the form
    1051             :  * [N,D,m,q,g,sqrt]. For successive D, find a square root of D mod N (g is a
    1052             :  * quadratic nonresidue), solve U^2 + |D|V^2 = 4N, then find candidates for m.
    1053             :  * When enough m's, batch-factor them to produce the q's. If one of the q's is
    1054             :  * pseudoprime, recursive call with N = q. May return gen_0 at toplevel
    1055             :  * => N has a small prime divisor */
    1056             : static GEN
    1057        1120 : N_downrun(GEN N, GEN param, GEN worker, GEN *X0, long *depth, long persevere, long stopat)
    1058             : {
    1059             :   pari_timer T, ti;
    1060        1120 :   long lgdisclist, lprimelist, nbsqrt = 0, t, i, j, expiN = expi(N);
    1061        1120 :   long persevere_next = 0, FAIL = 0;
    1062             :   ulong maxpcdg;
    1063             :   GEN primelist, uprimelist, disclist, sqrtlist, g, Dmbatch, kroP, Np;
    1064             : 
    1065        1120 :   dbg_mode() timer_start(&T);
    1066        1120 :   (*depth)++; /* we're going down the tree. */
    1067        1120 :   maxpcdg = ecpp_param_get_maxpcdg(param);
    1068        1120 :   primelist = ecpp_param_get_primelist(param);
    1069        1120 :   uprimelist = ecpp_param_get_uprimelist(param);
    1070        1120 :   disclist = ecpp_param_get_disclist(param);
    1071        1120 :   lgdisclist = lg(disclist);
    1072        1120 :   t = tunevec_batchsize(expiN, param);
    1073             : 
    1074             :   /* Precomputation for taking square roots, g needed for Fp_sqrt_i */
    1075        1120 :   g = Fp_2gener(N); if (!g) return gen_0; /* Composite if this happens. */
    1076             : 
    1077             :   /* Initialize sqrtlist for this N. */
    1078        1120 :   lprimelist = lg(primelist);
    1079        1120 :   sqrtlist = zerovec(lprimelist-1);
    1080             : 
    1081             :   /* A2: Check (p*|N) = 1 for all p */
    1082        1120 :   dbg_mode() timer_start(&ti);
    1083             :   /* kroP[i] will contain (primelist[i] | N) */
    1084        1120 :   kroP = cgetg(lprimelist,t_VECSMALL);
    1085        1120 :   switch(mod8(N))
    1086             :   { /* primelist[1,2,3] = [8,-4,-8]; N is odd */
    1087         350 :     case 1: kroP[1] = kroP[2] = kroP[3] = 1; break;
    1088         259 :     case 3: kroP[1] = -1; kroP[2] =-1; kroP[3] = 1; break;
    1089         336 :     case 5: kroP[1] = -1; kroP[2] = 1; kroP[3] =-1; break;
    1090         175 :     case 7: kroP[1] =  1; kroP[2] =-1; kroP[3] =-1; break;
    1091             :   }
    1092        1120 :   Np = Z_ZV_mod_tree(N, gel(uprimelist,1), gel(uprimelist,2));
    1093      133161 :   for(i=4; i<lprimelist; i++)
    1094             :   {
    1095      132041 :     if (Np[i]==0) return gen_0;
    1096      132041 :     kroP[i] = kross(Np[i],primelist[i]);
    1097             :   }
    1098        1120 :   dbg_mode() timer_record(X0, "A2", &ti);
    1099             : 
    1100             :   /* Print the start of this iteration. */
    1101        1120 :   dbg_mode() {
    1102           0 :     char o = persevere? '<': '[';
    1103           0 :     char c = persevere? '>': ']';
    1104           0 :     err_printf(ANSI_BRIGHT_CYAN "\n%c %3d | %4ld bits%c "
    1105             :                ANSI_RESET, o, *depth, expiN, c);
    1106             :   }
    1107             :   /* Initialize Dmbatch, populated with candidate cardinalities in Phase I
    1108             :    * (until #Dmbatch >= t); its elements will be factored on Phase II */
    1109        1120 :   Dmbatch = vectrunc_init(t+7);
    1110             : 
    1111             :   /* Number of cardinalities so far; should always be equal to lg(Dmbatch)-1. */
    1112             :   /* i determines which discriminant we are considering. */
    1113        1120 :   i = 1;
    1114        3913 :   while (!FAIL)
    1115             :   {
    1116             :     pari_timer F;
    1117        3906 :     long lgDmqlist, last_i = i, numcard = 0; /* #Dmbatch */
    1118             :     GEN Dmqlist;
    1119             : 
    1120             :     /* Dmbatch begins "empty", but keep the allocated memory. */
    1121        3906 :     setlg(Dmbatch, 1);
    1122             : 
    1123             :     /* PHASE I: Go through the D's and search for candidate m's.
    1124             :      * Fill up Dmbatch until there are at least t elements. */
    1125      538104 :     while (i < lgdisclist )
    1126             :     {
    1127      538076 :       GEN Dinfo = gel(disclist, i);
    1128             :       long n;
    1129      538076 :       if (!persevere && Dinfo_get_pd(Dinfo) > maxpcdg) { FAIL = 1; break; }
    1130      538034 :       n = D_collectcards(N,param, X0, Dinfo, sqrtlist, g, Dmbatch, kroP, &nbsqrt);
    1131      539091 :       if (n < 0) return gen_0;
    1132      538034 :       last_i = i++;
    1133      538034 :       numcard += n; if (numcard >= t) break;
    1134             :     }
    1135             : 
    1136             :     /* We have exhausted disclist and there are no card. to be factored */
    1137        3906 :     if (numcard <= 0 && (FAIL || i >= lgdisclist)) break;
    1138             : 
    1139             :     /* PHASE II: Find the corresponding q's for the m's found. Use Dmbatch */
    1140             :     /* Find coresponding q of the candidate m's. */
    1141        3857 :     dbg_mode() timer_start(&F);
    1142        3857 :     Dmqlist = Dmbatch_factor_Dmqvec(N, expiN, X0, Dmbatch, param);
    1143        3857 :     if (Dmqlist == NULL) continue; /* none left => next discriminant. */
    1144             : 
    1145             :     /* If we are cheating by adding class numbers, sort by class number */
    1146        3857 :     if (Dinfo_get_pd(gel(disclist, last_i)) > maxpcdg)
    1147         385 :       gen_sort_inplace(Dmqlist, NULL, &sort_Dmq_by_cnum, NULL);
    1148             : 
    1149             :     /* Check pseudoprimality before going down. */
    1150        3857 :     lgDmqlist = lg(Dmqlist);
    1151        3899 :     for (j = 1; j < lgDmqlist; j++)
    1152             :     {
    1153             :       GEN z, Dmq, q;
    1154             :       struct pari_mt pt;
    1155             :       pari_timer ti2;
    1156        3899 :       long running, stop = 0, pending = 0;
    1157        3899 :       mt_queue_start_lim(&pt, worker, lgDmqlist-j);
    1158        3899 :       dbg_mode() timer_start(&ti2);
    1159       27515 :       while ((running = (!stop && j < lgDmqlist)) || pending)
    1160             :       {
    1161             :         long jj;
    1162             :         GEN done;
    1163       23616 :         mt_queue_submit(&pt, j, running ? mkvec(gmael(Dmqlist,j,3)) : NULL);
    1164       23616 :         done = mt_queue_get(&pt, &jj, &pending);
    1165       23616 :         if (done)
    1166             :         {
    1167       20942 :           if (done[1] == 0)
    1168       19784 :           { dbg_mode() err_printf(ANSI_WHITE "." ANSI_RESET); }
    1169        1158 :           else if (!stop || jj < stop)
    1170        1103 :             stop = jj;
    1171             :         }
    1172       23616 :         j++;
    1173             :       }
    1174        3899 :       mt_queue_end(&pt);
    1175        3899 :       dbg_mode() timer_record(X0, "B3", &ti2);
    1176        3899 :       if (!stop && j >= lgDmqlist) break;
    1177        1099 :       j = stop; Dmq = gel(Dmqlist,j); q = gel(Dmq,3);
    1178             : 
    1179        1099 :       dbg_mode() {
    1180           0 :         err_printf(ANSI_BRIGHT_RED "\n       %5ld bits " ANSI_RESET,
    1181           0 :                    expi(q)-expiN);
    1182           0 :         err_printf("  D = %ld", Dmq_get_D(Dmq));
    1183           0 :         err_printf(ANSI_BRIGHT_BLUE "  h = %ld" ANSI_RESET, Dmq_get_h(Dmq));
    1184             :       }
    1185             :       /* q is pseudoprime */
    1186        1099 :       if (expi(q) < stopat) z = gen_1; /* q is prime; sentinel */
    1187             :       else
    1188             :       {
    1189        1043 :         z = N_downrun(q, param, worker, X0, depth, persevere_next, stopat);
    1190        1043 :         if (!z) {
    1191          42 :           dbg_mode() { char o = persevere? '<': '[', c = persevere? '>': ']';
    1192           0 :                        err_printf(ANSI_CYAN "\n%c %3d | %4ld bits%c "
    1193             :                                   ANSI_RESET, o, *depth, expiN, c); }
    1194          42 :           continue; /* downrun failed */
    1195             :         }
    1196        2058 :         if (isintzero(z)) return z;
    1197             :       }
    1198        1057 :       return mkNDmqg(z, N, Dmq, g, sqrtlist); /* SUCCESS */
    1199             :     }
    1200        2800 :     if (i >= lgdisclist) break; /* discriminants exhausted: FAIL */
    1201        2793 :     if (Dinfo_get_pd(gel(disclist, last_i)) > maxpcdg)
    1202             :     {
    1203         364 :       dbg_mode() err_printf(ANSI_BRIGHT_RED "  !" ANSI_RESET);
    1204         364 :       persevere_next = 1;
    1205             :     }
    1206             :   }
    1207             :   /* FAILED: Out of discriminants. */
    1208          63 :   umael(*X0, 3, 1)++; /* FAILS++ */
    1209          63 :   dbg_mode() err_printf(ANSI_BRIGHT_RED "  X" ANSI_RESET);
    1210          63 :   (*depth)--; return NULL;
    1211             : }
    1212             : 
    1213             : /* x: the output of the (recursive) downrun function; return a vector
    1214             :  * whose components are [N, D, m, q, g] */
    1215             : static GEN
    1216          56 : ecpp_flattencert(GEN x, long depth)
    1217             : {
    1218          56 :   long i, l = depth+1;
    1219          56 :   GEN ret = cgetg(l, t_VEC);
    1220        1113 :   for (i = 1; i < l; i++) { gel(ret, i) = gel(x,1); x = gel(x,2); }
    1221          56 :   return ret;
    1222             : }
    1223             : 
    1224             : /* Call the first downrun then unravel the skeleton of the certificate.
    1225             :  * Assume expi(N) < 64. This returns one of the following:
    1226             :  * - a vector of the form [N, D, m, q, y]
    1227             :  * - gen_0 (if N is composite)
    1228             :  * - NULL (if FAIL) */
    1229             : static GEN
    1230          77 : ecpp_step1(GEN N, GEN param, GEN* X0, long stopat)
    1231             : {
    1232          77 :   pari_sp av = avma;
    1233          77 :   long depth = 0;
    1234          77 :   GEN worker = strtofunction("_ecpp_ispsp_worker");
    1235          77 :   GEN downrun = N_downrun(N, param, worker, X0, &depth, 1, stopat);
    1236          77 :   if (downrun == NULL) return gc_NULL(av);
    1237          56 :   if (typ(downrun) != t_VEC) return gen_0;
    1238          56 :   return gerepilecopy(av, ecpp_flattencert(downrun, depth));
    1239             : }
    1240             : 
    1241             : /* The input is an integer N.
    1242             :    It uses the precomputation step0 done in ecpp_step0. */
    1243             : static GEN
    1244          77 : ecpp_param(GEN N, GEN param, long stopat)
    1245             : {
    1246             : 
    1247             :   GEN step1, answer, Tv, Cv, X0;
    1248          77 :   pari_sp av = avma;
    1249             :   long i, j;
    1250             : 
    1251             :   /* Check if N is pseudoprime to begin with. */
    1252          77 :   if (!BPSW_psp(N)) return gen_0;
    1253             : 
    1254             :   /* Check if we should even prove it. */
    1255          77 :   if (expi(N) < stopat) return N;
    1256             : 
    1257             :   /* Timers and Counters */
    1258          77 :   Tv = mkvec4(zero_zv(5), zero_zv(3), zero_zv(4), zero_zv(1));
    1259          77 :   Cv = mkvec4(zero_zv(5), zero_zv(3), zero_zv(4), zero_zv(1));
    1260          77 :   X0 = mkvec3(Tv, Cv, zero_zv(1));
    1261             : 
    1262          77 :   step1 = ecpp_step1(N, param, &X0, stopat);
    1263          77 :   if (step1 == NULL) return gc_NULL(av);
    1264          56 :   if (typ(step1) != t_VEC) { set_avma(av); return gen_0; }
    1265             : 
    1266          56 :   answer = ecpp_step2(step1, &X0, ecpp_param_get_primelist(param));
    1267             : 
    1268          56 :   dbg_mode() {
    1269           0 :   for (i = 1; i < lg(Tv); i++)
    1270             :   {
    1271           0 :     GEN v = gel(Tv, i);
    1272           0 :     long l = lg(v);
    1273           0 :     for (j = 1; j < l; j++)
    1274             :     {
    1275           0 :       long t = umael3(X0,1, i,j), n = umael3(X0,2, i,j);
    1276           0 :       if (!t) continue;
    1277           0 :       err_printf("\n  %c%ld: %16ld %16ld %16.3f", 'A'+i-1,j, t,n, (double)t/n);
    1278             :     }
    1279             :   }
    1280           0 :     err_printf("\n" ANSI_BRIGHT_RED "\nFAILS: %16ld" ANSI_RESET "\n", umael(X0, 3, 1));
    1281             :   }
    1282          56 :   return gerepilecopy(av, answer);
    1283             : }
    1284             : 
    1285             : static const long ecpp_tune[][4]=
    1286             : {
    1287             :   {  100,  2,  20,   500 },
    1288             :   {  350, 14,  21,  1000 },
    1289             :   {  450, 18,  21,  1500 },
    1290             :   {  750, 22,  22,  2000 },
    1291             :   {  900, 26,  23,  2500 },
    1292             :   { 1000, 32,  23,  3000 },
    1293             :   { 1200, 38,  24,  3500 },
    1294             :   { 1400, 44,  24,  4000 },
    1295             :   { 1600, 50,  24,  4500 },
    1296             :   { 1800, 56,  25,  5000 },
    1297             :   { 2000, 62,  25,  5500 },
    1298             :   { 2200, 68,  26,  6000 },
    1299             :   { 2400, 74,  26,  6500 },
    1300             :   { 2600, 80,  26,  7000 },
    1301             :   { 2800, 86,  26,  7500 },
    1302             :   { 3000, 92,  27,  8000 },
    1303             :   { 3200, 98,  27,  8500 },
    1304             :   { 3400, 104, 28,  9000 },
    1305             :   { 3600, 110, 28,  9500 },
    1306             :   { 3800, 116, 29, 10000 },
    1307             :   { 4000, 122, 29,     0 }
    1308             : };
    1309             : 
    1310             : /* assume N BPSW-pseudoprime */
    1311             : GEN
    1312          98 : ecpp0(GEN C, long stopat)
    1313             : {
    1314          98 :   pari_sp av = avma;
    1315             :   long expiN, i, tunelen;
    1316          98 :   GEN N = C, tune;
    1317             : 
    1318             :   /* Is it a partial certificate? */
    1319          98 :   if (typ(C) == t_VEC && check_ecppcert(C))
    1320           7 :     N = cert_get_qlast(C);
    1321          98 :   if (typ(N) != t_INT)
    1322           0 :     pari_err_TYPE("ecpp",C);
    1323             : 
    1324             :   /* Check if we should even prove it. */
    1325          98 :   expiN = expi(N);
    1326          98 :   if (stopat < 64) stopat = 64;
    1327          98 :   if (expiN < stopat) return C;
    1328             : 
    1329          56 :   tunelen = (expiN+499)/500;
    1330          56 :   tune = cgetg(tunelen+1, t_VEC);
    1331         140 :   for (i = 1; i <= tunelen && ecpp_tune[i-1][3]; i++)
    1332          84 :     gel(tune,i) = mkvecsmall4(ecpp_tune[i-1][0], ecpp_tune[i-1][1],
    1333          84 :                               ecpp_tune[i-1][2], ecpp_tune[i-1][3]);
    1334          56 :   for (; i <= tunelen; i++) gel(tune,i) = mkvecsmall4(100*(i+20),4*i+42,30,500*i);
    1335             :   for(;;)
    1336          21 :   {
    1337          77 :     GEN C2, param, x = gel(tune, tunelen);
    1338          77 :     param = ecpp_param_set(tune, x);
    1339          77 :     C2 = ecpp_param(N, param, stopat);
    1340          77 :     if (C2)
    1341          56 :       return typ(C)==t_INT? C2: gerepilecopy(av, shallowconcat(C,C2));
    1342          21 :     x[1] *= 2;
    1343          21 :     x[2] *= 2;
    1344          21 :     x[3] = minss(x[3]+1, 30);
    1345             :   }
    1346             : }
    1347             : 
    1348             : GEN
    1349          28 : ecpp(GEN N)
    1350          28 : { return ecpp0(N, 0); }
    1351             : 
    1352             : long
    1353          42 : isprimeECPP(GEN N)
    1354             : {
    1355          42 :   pari_sp av = avma;
    1356          42 :   if (!BPSW_psp(N)) return 0;
    1357          28 :   return gc_long(av, !isintzero(ecpp(N)));
    1358             : }
    1359             : 
    1360             : /* PARI ECPP Certificate -> Human-readable format */
    1361             : static GEN
    1362           7 : cert_out(GEN x)
    1363             : {
    1364           7 :   long i, l = lg(x);
    1365             :   pari_str str;
    1366           7 :   str_init(&str, 1);
    1367           7 :   if (typ(x) == t_INT)
    1368             :   {
    1369           0 :     str_printf(&str, "%Ps is prime.\nIndeed, ispseudoprime(%Ps) = 1 and %Ps < 2^64.\n", x, x, x);
    1370           0 :     return strtoGENstr(str.string);
    1371             :   }
    1372          21 :   for (i = 1; i < l; i++)
    1373             :   {
    1374          14 :     GEN xi = gel(x, i);
    1375          14 :     str_printf(&str, "[%ld]\n", i);
    1376          14 :     str_printf(&str, " N = %Ps\n", cert_get_N(xi));
    1377          14 :     str_printf(&str, " t = %Ps\n", cert_get_t(xi));
    1378          14 :     str_printf(&str, " s = %Ps\n", cert_get_s(xi));
    1379          14 :     str_printf(&str, "a = %Ps\n", cert_get_a4(xi));
    1380          14 :     str_printf(&str, "D = %Ps\n", cert_get_D(xi));
    1381          14 :     str_printf(&str, "m = %Ps\n", cert_get_m(xi));
    1382          14 :     str_printf(&str, "q = %Ps\n", cert_get_q(xi));
    1383          14 :     str_printf(&str, "E = %Ps\n", cert_get_E(xi));
    1384          14 :     str_printf(&str, "P = %Ps\n", cert_get_P(xi));
    1385             :   }
    1386           7 :   return strtoGENstr(str.string);
    1387             : }
    1388             : 
    1389             : /* PARI ECPP Certificate -> Magma Certificate
    1390             :  * [* [*
    1391             :  *     N, |D|, -1, m,
    1392             :  *     [* a, b *],
    1393             :  *     [* x, y, 1 *],
    1394             :  *     [* [* q, 1 *] *]
    1395             :  *   *], ... *] */
    1396             : static GEN
    1397          14 : magma_out(GEN x)
    1398             : {
    1399          14 :   long i, n = lg(x)-1;
    1400             :   pari_str ret;
    1401          14 :   str_init(&ret, 1);
    1402          14 :   if (typ(x) == t_INT)
    1403             :   {
    1404           0 :     str_printf(&ret, "Operation not supported.");
    1405           0 :     return strtoGENstr(ret.string);
    1406             :   }
    1407          14 :   str_printf(&ret, "[* ");
    1408         168 :   for (i = 1; i<=n; i++)
    1409             :   {
    1410         154 :     GEN xi = gel(x,i), E = cert_get_E(xi), P = cert_get_P(xi);
    1411         154 :     str_printf(&ret, "[* %Ps, %Ps, -1, %Ps, ",
    1412             :               cert_get_N(xi), negi(cert_get_D(xi)), cert_get_m(xi));
    1413         154 :     str_printf(&ret, "[* %Ps, %Ps *], ", gel(E,1), gel(E,2));
    1414         154 :     str_printf(&ret, "[* %Ps, %Ps, 1 *], ", gel(P,1), gel(P,2));
    1415         154 :     str_printf(&ret, "[* [* %Ps, 1 *] *] *]", cert_get_q(xi));
    1416         154 :     if (i != n) str_printf(&ret, ", ");
    1417             :   }
    1418          14 :   str_printf(&ret, " *]");
    1419          14 :   return strtoGENstr(ret.string);
    1420             : }
    1421             : 
    1422             : /* Prints: label=(sign)hexvalue\n
    1423             :  *   where sign is + or -
    1424             :  *   hexvalue is value in hex, of the form: 0xABC123 */
    1425             : static void
    1426         721 : primo_printval(pari_str *ret, const char* label, GEN value)
    1427             : {
    1428         721 :   str_printf(ret, label);
    1429         721 :   if (signe(value) >= 0) str_printf(ret, "=0x%P0X\n", value);
    1430         175 :   else str_printf(ret, "=-0x%P0X\n", negi(value));
    1431         721 : }
    1432             : 
    1433             : /* Input: PARI ECPP Certificate / Output: PRIMO Certificate
    1434             :  *
    1435             :  * Let Y^2 = X^3 + aX + b be the elliptic curve over N with the point (x,y)
    1436             :  * as described in the PARI certificate.
    1437             :  *
    1438             :  * If J != 0, 1728, PRIMO asks for [J,T], where T = a/A * B/b * x,
    1439             :  * A = 3J(1728-J) and B = 2J(1728-J)^2.
    1440             :  *
    1441             :  * If J=0 or J=1728, PRIMO asks for [A,B,T]; we let A=a and B=b => T = x */
    1442             : static GEN
    1443          14 : primo_out(GEN x)
    1444             : {
    1445          14 :   long i, l = (typ(x) == t_INT) ? 1 : lg(x);
    1446             :   pari_str ret;
    1447          14 :   str_init(&ret, 1);
    1448          14 :   str_printf(&ret, "[PRIMO - Primality Certificate]\nFormat=4\n");
    1449          14 :   str_printf(&ret, "TestCount=%ld\n\n[Comments]\n", l-1);
    1450          14 :   str_printf(&ret, "Generated by %s\n",paricfg_version);
    1451          14 :   str_printf(&ret, "https://pari.math.u-bordeaux.fr/\n\n[Candidate]\n");
    1452          14 :   if (typ(x) == t_INT)
    1453             :   {
    1454           0 :     str_printf(&ret, "N=0x%P0X\n", x);
    1455           0 :     return strtoGENstr(ret.string);
    1456             :   }
    1457          14 :   str_printf(&ret, "N=0x%P0X\n\n", cert_get_N(gel(x,1)));
    1458         168 :   for (i = 1; i < l; i++)
    1459             :   {
    1460         154 :     GEN a4, a6, N, Nover2, xi = gel(x,i);
    1461             :     long Ais0, Bis0;
    1462         154 :     str_printf(&ret, "[%ld]\n", i);
    1463         154 :     N = cert_get_N(xi);
    1464         154 :     Nover2 = shifti(N, -1);
    1465         154 :     primo_printval(&ret, "S", cert_get_s(xi));
    1466         154 :     primo_printval(&ret, "W", cert_get_t(xi));
    1467         154 :     a4 = cert_get_a4(xi); Ais0 = isintzero(a4);
    1468         154 :     a6 = cert_get_a6(xi); Bis0 = isintzero(a6);
    1469         154 :     if (!Ais0 && !Bis0) { /* J != 0, 1728 */
    1470          49 :       primo_printval(&ret, "J", Fp_center_i(cert_get_J(xi), N, Nover2));
    1471          49 :       primo_printval(&ret, "T", cert_get_T(xi));
    1472             :     } else {
    1473         105 :       if (!Ais0) a4 = Fp_center_i(a4, N, Nover2);
    1474         105 :       if (!Bis0) a6 = Fp_center_i(a6, N, Nover2);
    1475         105 :       primo_printval(&ret, "A", a4);
    1476         105 :       primo_printval(&ret, "B", a6);
    1477         105 :       primo_printval(&ret, "T", cert_get_x(xi));
    1478             :     }
    1479         154 :     str_printf(&ret, "\n");
    1480             :   }
    1481          14 :   return strtoGENstr(ret.string);
    1482             : }
    1483             : 
    1484             : /* return 1 if q > (N^1/4 + 1)^2, 0 otherwise */
    1485             : static long
    1486         503 : Nq_isvalid(GEN N, GEN q)
    1487             : {
    1488         503 :   GEN c = subii(sqri(subiu(q,1)), N);
    1489         503 :   if (signe(c) <= 0) return 0;
    1490             :   /* (q-1)^2 > N; check that (N - (q-1)^2)^2 > 16Nq */
    1491         503 :   return (cmpii(sqri(c), shifti(mulii(N,q), 4)) > 0);
    1492             : }
    1493             : 
    1494             : GEN
    1495         511 : primecertisvalid_ecpp_worker(GEN certi)
    1496             : {
    1497             :   GEN N, W, mP, sP, r, m, s, a, P, q;
    1498         511 :   if (lg(certi)-1 != 5) return gen_0;
    1499             : 
    1500         511 :   N = cert_get_N(certi);
    1501         511 :   if (typ(N) != t_INT || signe(N) <= 0) return gen_0;
    1502         511 :   switch(umodiu(N,6))
    1503             :   {
    1504         504 :     case 1: case 5: break; /* Check if N is not divisible by 2 or 3 */
    1505           7 :     default: return gen_0;
    1506             :   }
    1507             : 
    1508         504 :   W = cert_get_t(certi);
    1509         504 :   if (typ(W) != t_INT) return gen_0;
    1510             :   /* Check if W^2 < 4N (Hasse interval) */
    1511         504 :   if (!(cmpii(sqri(W), shifti(N,2)) < 0)) return gen_0;
    1512             : 
    1513         504 :   s = cert_get_s(certi);
    1514         504 :   if (typ(s) != t_INT || signe(s) < 0) return gen_0;
    1515             : 
    1516         504 :   m = cert_get_m(certi);
    1517         503 :   q = dvmdii(m, s, &r);
    1518             :   /* Check m%s == 0 */
    1519         503 :   if (!isintzero(r)) return gen_0;
    1520             :   /* Check q > (N^(1/4) + 1)^2 */
    1521         503 :   if (!Nq_isvalid(N, q)) return gen_0;
    1522             : 
    1523         504 :   a = cert_get_a4(certi);
    1524         504 :   if (typ(a) != t_INT) return gen_0;
    1525             : 
    1526         504 :   P = cert_get_P(certi);
    1527         504 :   if (typ(P) != t_VEC || lg(P) != 3) return gen_0;
    1528         504 :   P = FpE_to_FpJ(P);
    1529             : 
    1530             :   /* Check mP == 0 */
    1531         504 :   mP = FpJ_mul(P, m, a, N);
    1532         504 :   if (!FpJ_is_inf(mP)) return gen_0;
    1533             : 
    1534             :   /* Check sP != 0 and third component is coprime to N */
    1535         504 :   sP = FpJ_mul(P, s, a, N);
    1536         504 :   if (!isint1(gcdii(gel(sP, 3), N))) return gen_0;
    1537         504 :   return q;
    1538             : }
    1539             : 
    1540             : /* return 1 if 'cert' is a valid PARI ECPP certificate */
    1541             : static long
    1542          35 : ecppisvalid_i(GEN cert)
    1543             : {
    1544          35 :   const long trustbits = 64;/* a pseudoprime < 2^trustbits is prime */
    1545          35 :   long i, lgcert = lg(cert);
    1546          35 :   GEN ql, q = gen_0, worker, check;
    1547             : 
    1548          35 :   if (typ(cert) == t_INT)
    1549             :   {
    1550           0 :     if (expi(cert) >= trustbits) return 0; /* >= 2^trustbits */
    1551           0 :     return BPSW_psp(cert);
    1552             :   }
    1553          35 :   if (typ(cert) != t_VEC || lgcert <= 1) return 0;
    1554          35 :   ql = veclast(cert); if (lg(ql)-1 != 5) return 0;
    1555          35 :   ql = cert_get_q(ql);
    1556          35 :   if (expi(ql) >= trustbits || !BPSW_psp(ql)) return 0;
    1557          35 :   worker = strtofunction("_primecertisvalid_ecpp_worker");
    1558          35 :   check = gen_parapply(worker, cert);
    1559         490 :   for (i = 1; i < lgcert; i++)
    1560             :   {
    1561         462 :     GEN certi = gel(cert, i);
    1562         462 :     GEN qq = gel(check,i), N = cert_get_N(certi);
    1563         462 :     if (isintzero(qq)) return 0;
    1564         455 :     if (i > 1 && !equalii(N, q)) return 0; /* N of this entry doesn't match q of previous */
    1565         455 :     q = qq;
    1566             :   }
    1567          28 :   return 1;
    1568             : }
    1569             : long
    1570          35 : ecppisvalid(GEN cert)
    1571          35 : { pari_sp av = avma; return gc_long(av, ecppisvalid_i(cert)); }
    1572             : 
    1573             : GEN
    1574          35 : ecppexport(GEN cert, long flag)
    1575             : {
    1576          35 :   switch(flag)
    1577             :   {
    1578           7 :     case 0: return cert_out(cert);
    1579          14 :     case 1: return primo_out(cert);
    1580          14 :     case 2: return magma_out(cert);
    1581             :   }
    1582           0 :   pari_err_FLAG("primecertexport");
    1583             :   return NULL;/*LCOV_EXCL_LINE*/
    1584             : }

Generated by: LCOV version 1.14