Code coverage tests

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

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

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

LCOV - code coverage report
Current view: top level - basemath - polclass.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.10.0 lcov report (development 21059-cbe0d6a) Lines: 925 961 96.3 %
Date: 2017-09-22 06:24:58 Functions: 57 58 98.3 %
Legend: Lines: hit not hit

          Line data    Source code
       1             : /* Copyright (C) 2014  The PARI group.
       2             : 
       3             : This file is part of the PARI/GP package.
       4             : 
       5             : PARI/GP is free software; you can redistribute it and/or modify it under the
       6             : terms of the GNU General Public License as published by the Free Software
       7             : Foundation. It is distributed in the hope that it will be useful, but WITHOUT
       8             : ANY WARRANTY WHATSOEVER.
       9             : 
      10             : Check the License for details. You should have received a copy of it, along
      11             : with the package; see the file 'COPYING'. If not, write to the Free Software
      12             : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
      13             : 
      14             : #include "pari.h"
      15             : #include "paripriv.h"
      16             : 
      17             : #define dbg_printf(lvl) if (DEBUGLEVEL >= (lvl) + 3) err_printf
      18             : 
      19             : /**
      20             :  * SECTION: Functions dedicated to finding a j-invariant with a given
      21             :  * trace.
      22             :  */
      23             : 
      24             : /* TODO: This code is shared with
      25             :  * torsion_compatible_with_characteristic() in 'torsion.c'. */
      26             : static void
      27      115046 : hasse_bounds(long *low, long *high, long p)
      28             : {
      29      115046 :   long two_sqrt_p = usqrt(4*p);
      30      115046 :   *low = p + 1 - two_sqrt_p;
      31      115046 :   *high = p + 1 + two_sqrt_p;
      32      115046 : }
      33             : 
      34             : 
      35             : /*
      36             :  * a and b must be the result of factoru_pow(), and b must divide a
      37             :  * exactly.
      38             :  */
      39             : INLINE void
      40      233761 : famatsmall_divexact(GEN a, GEN b)
      41             : {
      42             :   long i, j;
      43      235085 :   for (i = j = 1; j < lg(gel(a, 1)) && i < lg(gel(b, 1)); ++j)
      44        1324 :     if (gel(a, 1)[j] == gel(b, 1)[i])
      45        1171 :       gel(a, 2)[j] -= gel(b, 2)[i++];
      46             : 
      47      939397 :   for (i = j = 1; j < lg(gel(a, 1)); ++j) {
      48      705636 :     if (gel(a, 2)[j]) {
      49      705180 :       gel(a, 1)[i] = gel(a, 1)[j];
      50      705180 :       gel(a, 2)[i] = gel(a, 2)[j];
      51      705180 :       ++i;
      52             :     }
      53             :   }
      54      233761 :   if (i == 1) {
      55             :     /* b == a, so a must now be 1. */
      56           6 :     gel(a, 1)[1] = 1;
      57           6 :     gel(a, 2)[1] = 0;
      58           6 :     setlg(gel(a, 1), 2);
      59           6 :     setlg(gel(a, 2), 2);
      60             :   } else {
      61      233755 :     setlg(gel(a, 1), i);
      62      233755 :     setlg(gel(a, 2), i);
      63             :   }
      64      233761 : }
      65             : 
      66             : 
      67             : /*
      68             :  * This is Sutherland, 2009, TestCurveOrder.
      69             :  *
      70             :  * [a4, a6] and p specify an elliptic curve over FF_p.  N is a
      71             :  * two-element array containing the two possible curve orders, and n
      72             :  * is a two-element array containg the corresponding factorisations as
      73             :  * famats.
      74             :  */
      75             : static long
      76      136901 : test_curve_order(
      77             :   norm_eqn_t ne, ulong a4, ulong a6,
      78             :   long N0, long N1, GEN n0, GEN n1,
      79             :   const long hasse[2])
      80             : {
      81      136901 :   pari_sp ltop = avma, av;
      82             :   ulong a4t, a6t;
      83             :   long m0, m1;
      84             :   long hasse_low, hasse_high;
      85      136901 :   ulong p = ne->p, pi = ne->pi, T = ne->T;
      86      136901 :   ulong swapped = 0;
      87             : 
      88      136901 :   if (p <= 11) {
      89          93 :     long card = (long)p + 1 - Fl_elltrace(a4, a6, p);
      90          93 :     return card == N0 || card == N1;
      91             :   }
      92             : 
      93             :   /* [a4, a6] is the given curve and [a4t, a6t] is its quadratic
      94             :    * twist. */
      95      136808 :   Fl_elltwist_disc(a4, a6, T, p, &a4t, &a6t);
      96             : 
      97      136808 :   m0 = m1 = 1;
      98             : 
      99      136808 :   if (N0 + N1 != 2 * (long)p + 2) pari_err_BUG("test_curve_order");
     100             : 
     101      136808 :   hasse_low = hasse[0];
     102      136808 :   hasse_high = hasse[1];
     103             : 
     104      136808 :   av = avma;
     105             :   for ( ; ; ) {
     106             :     GEN pt, Q, tmp;
     107             :     long a1, x, n_s;
     108             : 
     109      233761 :     pt = random_Flj_pre(a4, a6, p, pi);
     110      233761 :     Q = Flj_mulu_pre(pt, m0, a4, p, pi);
     111             :     /* TODO: Work out how to avoid this copying. */
     112      233761 :     tmp = gcopy(n0);
     113      233761 :     famatsmall_divexact(tmp, factoru(m0));
     114      233761 :     n_s = Flj_order_ufact(Q, N0 / m0, tmp, a4, p, pi);
     115             : 
     116      233761 :     if (n_s == 0) {
     117             :       /* If m0 divides N1 and m1 divides N0 and N0 < N1,
     118             :        * then swap. */
     119      103187 :       if ( ! swapped && N1 % m0 == 0 && N0 % m1 == 0) {
     120       81347 :         swapspec(n0, n1, N0, N1);
     121       81347 :         swapped = 1;
     122       81347 :         continue;
     123             :       } else {
     124       21840 :         avma = ltop;
     125       21840 :         return 0;
     126             :       }
     127             :     }
     128             : 
     129      130574 :     m0 *= n_s;
     130      130574 :     a1 = (2 * p + 2) % m1;
     131             :     /* Using ceil(n/d) = (n + d - 1)/d */
     132      130574 :     x = (hasse_low + m0 - 1) / m0;
     133      130574 :     x *= m0;
     134      255404 :     for ( ; x <= hasse_high; x += m0) {
     135      140436 :       if ((x % m1) == a1 && x != N0 && x != N1)
     136       15606 :         break;
     137             :     }
     138             :     /* We exited the loop because we finished iterating, not because
     139             :      * of the break.  That means every x in N was either N0 or N1, so
     140             :      * we return true. */
     141      130574 :     if (x > hasse_high) {
     142      114968 :       avma = ltop;
     143      114968 :       return 1;
     144             :     }
     145             : 
     146       15606 :     lswap(a4, a4t);
     147       15606 :     lswap(a6, a6t);
     148       15606 :     lswap(m0, m1);
     149       15606 :     avma = av;
     150       96953 :   }
     151             : }
     152             : 
     153             : INLINE int
     154     3324418 : jac_eq_or_opp(GEN P, GEN Q, ulong p, ulong pi)
     155             : {
     156             :   /* (X1:Y1:Z1) and (X2:Y2:Z2) in Jacobian coordinates are equal
     157             :    * or opposite iff X1 Z2^2 = X2 Z1^2. */
     158     6648836 :   return ! Fl_sub(Fl_mul_pre(P[1], Fl_sqr_pre(Q[3], p, pi), p, pi),
     159     6648836 :                   Fl_mul_pre(Q[1], Fl_sqr_pre(P[3], p, pi), p, pi), p);
     160             : }
     161             : 
     162             : 
     163             : /* This is Sutherland 2009 Algorithm 1.1 */
     164             : static long
     165      115053 : find_j_inv_with_given_trace(
     166             :   ulong *j_t, norm_eqn_t ne, long rho_inv, long max_curves)
     167             : {
     168      115053 :   pari_sp ltop = avma, av;
     169      115053 :   long curves_tested = 0, batch_size;
     170             :   long N0, N1, hasse[2];
     171             :   GEN n0, n1;
     172      115053 :   long i, found = 0;
     173      115053 :   ulong p = ne->p, pi = ne->pi;
     174      115053 :   long t = ne->t;
     175      115053 :   ulong p1 = p + 1, j = 0, m, c_1728 = 1728 % p;
     176             :   GEN A4, A6, tx, ty;
     177             :   /* This number must be the same as LAST_X1_LEVEL in 'torsion.c', */
     178             :   enum { MAX_X1_CURVE_LVL = 39 };
     179             : 
     180             :   /* ellap(ellinit(ellfromj(Mod(1,2)))) == -1
     181             :    * ellap(ellinit(ellfromj(Mod(1,3)))) ==  1
     182             :    * ellap(ellinit(ellfromj(Mod(2,3)))) ==  2 */
     183      115053 :   if (p == 2 || p == 3) {
     184           7 :     if (t == 0) pari_err_BUG("find_j_inv_with_given_trace");
     185           7 :     *j_t = t;
     186           7 :     return 1;
     187             :   }
     188             : 
     189      115046 :   N0 = (long)p1 - t;
     190      115046 :   N1 = (long)p1 + t;
     191      115046 :   n0 = factoru(N0);
     192      115046 :   n1 = factoru(N1);
     193             : 
     194             :   /* FIXME: Select m more intelligently.  Currently just the biggest
     195             :    * common divisor of N0 and N1 less than 39. */
     196      115046 :   m = cgcd(N0, N1);
     197      115046 :   av = avma;
     198      115046 :   if (m > MAX_X1_CURVE_LVL) {
     199       16651 :     GEN factm = factoru(m);
     200       16651 :     long nfactors = lg(gel(factm, 1)) - 1;
     201       45559 :     for (i = 1; i <= nfactors; ) {
     202       28908 :       m /= gel(factm, 1)[i];
     203       28908 :       if (m <= MAX_X1_CURVE_LVL)
     204       16651 :         break;
     205       12257 :       gel(factm, 2)[i] -= 1;
     206       12257 :       if (gel(factm, 2)[i] == 0)
     207        1058 :         ++i;
     208             :     }
     209       16651 :     avma = av;
     210             :   }
     211             : 
     212             :   /* Select batch size so that we have roughly a 50% chance of finding
     213             :    * a good curve in a batch. */
     214      115046 :   batch_size = 1.0 + rho_inv / (2.0 * m);
     215      115046 :   A4 = cgetg(batch_size + 1, t_VECSMALL);
     216      115046 :   A6 = cgetg(batch_size + 1, t_VECSMALL);
     217      115046 :   tx = cgetg(batch_size + 1, t_VECSMALL);
     218      115046 :   ty = cgetg(batch_size + 1, t_VECSMALL);
     219             : 
     220      115046 :   dbg_printf(2)("  Selected torsion constraint m = %lu and batch "
     221             :              "size = %ld\n", m, batch_size);
     222             : 
     223      115046 :   hasse_bounds(&hasse[0], &hasse[1], p);
     224             : 
     225      115046 :   av = avma;
     226      687814 :   while ( ! found && (max_curves <= 0 || curves_tested < max_curves)) {
     227      457722 :     random_curves_with_m_torsion((ulong *)(A4 + 1), (ulong *)(A6 + 1),
     228             :                                  (ulong *)(tx + 1), (ulong *)(ty + 1),
     229             :                                  batch_size, m, p);
     230     3669089 :     for (i = 1; i <= batch_size; ++i) {
     231             :       ulong a4, a6;
     232             :       GEN P, p1P, tP;
     233     3326413 :       ++curves_tested;
     234     3326413 :       a4 = A4[i];
     235     3326413 :       a6 = A6[i];
     236     3326413 :       j = Fl_ellj_pre(a4, a6, p, pi);
     237     3326413 :       if (j == 0 || j == c_1728)
     238        1995 :         continue;
     239             : 
     240     3324418 :       P = random_Flj_pre(a4, a6, p, pi);
     241     3324418 :       p1P = Flj_mulu_pre(P, p1, a4, p, pi);
     242     3324418 :       tP = Flj_mulu_pre(P, t, a4, p, pi);
     243             : 
     244     3324418 :       if (jac_eq_or_opp(p1P, tP, p, pi)
     245      136901 :           && test_curve_order(ne, a4, a6, N0, N1, n0, n1, hasse)) {
     246      115046 :         found = 1;
     247      115046 :         break;
     248             :       }
     249     3209372 :       avma = av;
     250             :     }
     251             :   }
     252      115046 :   *j_t = j;
     253      115046 :   avma = ltop;
     254      115046 :   return curves_tested;
     255             : }
     256             : 
     257             : 
     258             : /**
     259             :  * SECTION: Functions for dealing with polycyclic presentations.
     260             :  */
     261             : 
     262             : static GEN
     263        4533 : next_generator(GEN DD, long D, ulong u, long filter, long *P)
     264             : {
     265        4533 :   pari_sp av = avma, av1;
     266             :   GEN gen, genred;
     267             :   long norm;
     268        4533 :   ulong p = (ulong)*P;
     269             :   while (1) {
     270        7622 :     p = unextprime(p + 1);
     271        7622 :     if (p > LONG_MAX) pari_err_BUG("next_generator");
     272        7622 :     if (kross(D, (long)p) != -1 && u % p != 0 && filter % p != 0) {
     273        4533 :       gen = primeform_u(DD, p);
     274        4533 :       av1 = avma;
     275             : 
     276             :       /* If the reduction of gen has norm = 1, then it is the identity
     277             :        * form and is therefore skipped. */
     278        4533 :       genred = redimag(gen);
     279        4533 :       norm = itos(gel(genred, 1));
     280        4533 :       avma = av1;
     281        4533 :       if (norm != 1)
     282        4533 :         break;
     283           0 :       avma = av;
     284             :     }
     285        3089 :   }
     286        4533 :   *P = (long)p;
     287        4533 :   return gen;
     288             : }
     289             : 
     290             : INLINE long *
     291        5359 : evec_ri_mutate(long r[], long i)
     292             : {
     293        5359 :   return r + (i * (i - 1) >> 1);
     294             : }
     295             : 
     296             : INLINE const long *
     297       12647 : evec_ri(const long r[], long i)
     298             : {
     299       12647 :   return r + (i * (i - 1) >> 1);
     300             : }
     301             : 
     302             : /* Reduces evec e so that e[i] < n[i] (e[i] are assumed to be
     303             :  * nonnegative) using pcp (n,r,k). Does not check for overflow -- this
     304             :  * could be an issue for large groups. */
     305             : INLINE void
     306       18579 : evec_reduce(long e[], const long n[], const long r[], long k)
     307             : {
     308             :   long i, j, q;
     309             :   const long *ri;
     310       18579 :   if (!k)
     311       18579 :     return;
     312       46017 :   for (i = k - 1; i > 0; i--) {
     313       27438 :     if (e[i] >= n[i]) {
     314        8094 :       q = e[i] / n[i];
     315        8094 :       ri = evec_ri(r, i);
     316       20578 :       for (j = 0; j < i; j++)
     317       12484 :         e[j] += q * ri[j];
     318        8094 :       e[i] -= q * n[i];
     319             :     }
     320             :   }
     321       18579 :   e[0] %= n[0];
     322             : }
     323             : 
     324             : /* Computes e3 = log(a^e1*a^e2) in terms of the given polycyclic
     325             :  * presentation (here a denotes the implicit vector of generators) */
     326             : INLINE void
     327         518 : evec_compose(
     328             :   long e3[],
     329             :   const long e1[], const long e2[], const long n[],const long r[], long k)
     330             : {
     331             :     long i;
     332        2128 :     for (i = 0; i < k; i++)
     333        1610 :       e3[i] = e1[i] + e2[i];
     334         518 :     evec_reduce(e3, n, r, k);
     335         518 : }
     336             : 
     337             : /* Converts an evec to an integer index corresponding to the
     338             :  * multi-radix representation of the evec with moduli corresponding to
     339             :  * the subgroup orders m[i] */
     340             : INLINE long
     341        9408 : evec_to_index(const long e[], const long m[], long k)
     342             : {
     343        9408 :   long i, index = e[0];
     344       23282 :   for (i = 1; i < k; i++)
     345       13874 :     index += e[i] * m[i - 1];
     346        9408 :   return index;
     347             : }
     348             : 
     349             : INLINE void
     350       20297 : evec_copy(long f[], const long e[], long k)
     351             : {
     352             :   long i;
     353       44324 :   for (i = 0; i < k; ++i)
     354       24027 :     f[i] = e[i];
     355       20297 : }
     356             : 
     357             : INLINE void
     358        4025 : evec_clear(long e[], long k)
     359             : {
     360             :   long i;
     361       10553 :   for (i = 0; i < k; ++i)
     362        6528 :     e[i] = 0;
     363        4025 : }
     364             : 
     365             : /* e1 and e2 may overlap */
     366             : /* Note that this function is not very efficient because it does not
     367             :  * know the orders of the elements in the presentation, only the
     368             :  * relative orders */
     369             : INLINE void
     370         168 : evec_inverse(
     371             :   long e2[], const long e1[], const long n[], const long r[], long k)
     372             : {
     373         168 :   pari_sp av = avma;
     374             :   long i, *e3, *e4;
     375             : 
     376         168 :   e3 = new_chunk(k);
     377         168 :   e4 = new_chunk(k);
     378         168 :   evec_clear(e4, k);
     379         168 :   evec_copy(e3, e1, k);
     380             :   /* We have e1 + e4 = e3 which we maintain throughout while making e1
     381             :    * the zero vector */
     382         805 :   for (i = k - 1; i >= 0; i--) {
     383         637 :     if (e3[i]) {
     384          28 :       e4[i] += n[i] - e3[i];
     385          28 :       evec_reduce(e4, n, r, k);
     386          28 :       e3[i] = n[i];
     387          28 :       evec_reduce(e3, n, r, k);
     388             :     }
     389             :   }
     390         168 :   evec_copy(e2, e4, k);
     391         168 :   avma = av;
     392         168 : }
     393             : 
     394             : /* e1 and e2 may overlap */
     395             : /* This is a faster way to compute inverses, if the presentation
     396             :  * element orders are known (these are specified in the array o, the
     397             :  * array n holds the relative orders) */
     398             : INLINE void
     399         707 : evec_inverse_o(
     400             :   long e2[],
     401             :   const long e1[], const long n[], const long o[], const long r[], long k)
     402             : {
     403             :     long j;
     404        3227 :     for (j = 0; j < k; j++)
     405        2520 :       e2[j] = (e1[j] ? o[j] - e1[j] : 0);
     406         707 :     evec_reduce(e2, n, r, k);
     407         707 : }
     408             : 
     409             : /* Computes the order of the group element a^e using the pcp (n,r,k) */
     410             : INLINE long
     411        4421 : evec_order(const long e[], const long n[], const long r[], long k)
     412             : {
     413        4421 :   pari_sp av = avma;
     414        4421 :   long *f = new_chunk(k);
     415             :   long i, j, o, m;
     416             : 
     417        4421 :   evec_copy(f, e, k);
     418       11468 :   for (o = 1, i = k - 1; i >= 0; i--) {
     419        7047 :     if (f[i]) {
     420        5216 :       m = n[i] / cgcd(f[i], n[i]);
     421       14434 :       for (j = 0; j < k; j++)
     422        9218 :         f[j] *= m;
     423        5216 :       evec_reduce(f, n, r, k);
     424        5216 :       o *= m;
     425             :     }
     426             :   }
     427        4421 :   avma = av;
     428        4421 :   return o;
     429             : }
     430             : 
     431             : /* Computes orders o[] for each generator using relative orders n[]
     432             :  * and power relations r[] */
     433             : INLINE void
     434        3451 : evec_orders(long o[], const long n[], const long r[], long k)
     435             : {
     436        3451 :   pari_sp av = avma;
     437        3451 :   long i, *e = new_chunk(k);
     438             : 
     439        3451 :   evec_clear(e, k);
     440        7872 :   for (i = 0; i < k; i++) {
     441        4421 :     e[i] = 1;
     442        4421 :     if (i)
     443         970 :       e[i - 1] = 0;
     444        4421 :     o[i] = evec_order(e, n, r, k);
     445             :   }
     446        3451 :   avma = av;
     447        3451 : }
     448             : 
     449             : INLINE int
     450         231 : evec_equal(const long e1[], const long e2[], long k)
     451             : {
     452             :   long j;
     453         329 :   for (j = 0; j < k; ++j)
     454         308 :     if (e1[j] != e2[j])
     455         210 :       break;
     456         231 :   return j == k;
     457             : }
     458             : 
     459             : INLINE void
     460        1323 : index_to_evec(long e[], long index, const long m[], long k)
     461             : {
     462             :   long i;
     463        4046 :   for (i = k - 1; i > 0; --i) {
     464        2723 :     e[i] = index / m[i - 1];
     465        2723 :     index -= e[i] * m[i - 1];
     466             :   }
     467        1323 :   e[0] = index;
     468        1323 : }
     469             : 
     470             : INLINE void
     471        3451 : evec_n_to_m(long m[], const long n[], long k)
     472             : {
     473             :   long i;
     474        3451 :   m[0] = n[0];
     475        4421 :   for (i = 1; i < k; ++i)
     476         970 :     m[i] = m[i - 1] * n[i];
     477        3451 : }
     478             : 
     479             : #define HALFLOGPI 0.57236494292470008707171367567653
     480             : 
     481             : /*
     482             :  * Based on logfac() in Sutherland's classpoly package.
     483             :  *
     484             :  * Ramanujan approximation to log(n!), accurate to O(1/n^3)
     485             :  */
     486             : INLINE double
     487           0 : logfac(long n)
     488             : {
     489           0 :   return n * log((double) n) - (double) n +
     490           0 :     log((double) n * (1.0 + 4.0 * n * (1.0 + 2.0 * n))) / 6.0 +
     491             :     HALFLOGPI;
     492             : }
     493             : 
     494             : 
     495             : #define LOG2E 1.44269504088896340735992468100189
     496             : 
     497             : 
     498             : /* This is based on Sutherland 2009, Lemma 8 (p31). */
     499             : static double
     500        3550 : upper_bound_on_classpoly_coeffs(long D, long h, GEN qfinorms)
     501             : {
     502        3550 :   pari_sp ltop = avma;
     503        3550 :   GEN C = dbltor(2114.567);
     504             :   double Mk, m, logbinom;
     505        3550 :   GEN tmp = mulrr(mppi(LOWDEFAULTPREC), sqrtr(stor(-D, LOWDEFAULTPREC)));
     506             :   /* We treat this case separately since the table is not initialised
     507             :    * when h = 1. This is the same as in the for loop below but with ak
     508             :    * = 1. */
     509        3550 :   double log2Mk = dbllog2r(mpadd(mpexp(tmp), C));
     510        3550 :   double res = log2Mk;
     511        3550 :   ulong maxak = 1;
     512        3550 :   double log2Mh = log2Mk;
     513             : 
     514        3550 :   pari_sp btop = avma;
     515             :   long k;
     516       65568 :   for (k = 2; k <= h; ++k) {
     517       62018 :     ulong ak = uel(qfinorms, k);
     518             :     /* Unfortunately exp(tmp/a[k]) can overflow for even moderate
     519             :      * discriminants, so we need to do this calculation with t_REALs
     520             :      * instead of just doubles.  Sutherland has a (much more
     521             :      * complicated) implementation in the classpoly package which
     522             :      * should be consulted if this ever turns out to be a bottleneck.
     523             :      *
     524             :      * [Note that one idea to avoid t_REALs is the following: we have
     525             :      * log(e^x + C) - x <= log(2) ~ 0.69 for x >= log(C) ~ 0.44 and
     526             :      * the difference is basically zero for x slightly bigger than
     527             :      * log(C).  Hence for large discriminants, we will always have x =
     528             :      * \pi\sqrt{-D}/ak >> log(C) and so we could approximate log(e^x +
     529             :      * C) by x.] */
     530       62018 :     log2Mk = dbllog2r(mpadd(mpexp(divru(tmp, ak)), C));
     531       62018 :     res += log2Mk;
     532       62018 :     if (ak > maxak) {
     533       12599 :       maxak = ak;
     534       12599 :       log2Mh = log2Mk;
     535             :     }
     536       62018 :     avma = btop;
     537             :   }
     538             : 
     539        3550 :   Mk = pow(2.0, log2Mh);
     540        3550 :   m = floor((h + 1)/(Mk + 1.0));
     541             :   /* This line computes "log2(itos(binomialuu(h, m)))".  The smallest
     542             :    * fundamental discriminant for which logbinom is not zero is
     543             :    * -1579751. */
     544        3550 :   logbinom = (m > 0 && m < h)
     545           0 :     ? LOG2E * (logfac(h) - logfac(m) - logfac(h - m))
     546        3550 :     : 0;
     547        3550 :   avma = ltop;
     548        3550 :   return res + logbinom - m * log2Mh + 2.0;
     549             : }
     550             : 
     551             : INLINE long
     552         987 : distinct_inverses(
     553             :   const long f[], const long ef[], const long ei[],
     554             :   const long n[], const long o[], const long r[], long k,
     555             :   long L0, long i)
     556             : {
     557         987 :   pari_sp av = avma;
     558             :   long j, *e2, *e3;
     559             : 
     560         987 :   if ( ! ef[i] || (L0 && ef[0]))
     561         511 :     return 0;
     562         560 :   for (j = i + 1; j < k; ++j)
     563         455 :     if (ef[j])
     564         371 :       break;
     565         476 :   if (j < k)
     566         371 :     return 0;
     567             : 
     568         105 :   e2 = new_chunk(k);
     569         105 :   evec_copy(e2, ef, i);
     570         105 :   e2[i] = o[i] - ef[i];
     571         161 :   for (j = i + 1; j < k; ++j)
     572          56 :     e2[j] = 0;
     573         105 :   evec_reduce(e2, n, r, k);
     574             : 
     575         105 :   if (evec_equal(ef, e2, k)) {
     576           7 :     avma = av;
     577           7 :     return 0;
     578             :   }
     579             : 
     580          98 :   e3 = new_chunk(k);
     581          98 :   evec_inverse_o(e3, ef, n, o, r, k);
     582          98 :   if (evec_equal(e2, e3, k)) {
     583          14 :     avma = av;
     584          14 :     return 0;
     585             :   }
     586             : 
     587          84 :   if (f) {
     588          14 :     evec_compose(e3, f, ei, n, r, k);
     589          14 :     if (evec_equal(e2, e3, k)) {
     590           0 :       avma = av;
     591           0 :       return 0;
     592             :     }
     593             : 
     594          14 :     evec_inverse_o(e3, e3, n, o, r, k);
     595          14 :     if (evec_equal(e2, e3, k)) {
     596           0 :       avma = av;
     597           0 :       return 0;
     598             :     }
     599             :   }
     600          84 :   avma = av;
     601          84 :   return 1;
     602             : }
     603             : 
     604             : INLINE long
     605        1225 : next_prime_evec(
     606             :   long *qq, long f[], const long m[], long k,
     607             :   hashtable *tbl, long D, GEN DD, long u, long lvl, long ubound)
     608             : {
     609        1225 :   pari_sp av = avma;
     610             :   hashentry *he;
     611             :   GEN P;
     612        1225 :   long idx, q = *qq;
     613             : 
     614        3150 :   do q = unextprime(q + 1);
     615        6300 :   while (!(u % q) || kross(D, q) == -1
     616        4375 :       || !(lvl % q) || !(D % (q * q)));
     617        1225 :   if (q > ubound)
     618         126 :     return 0;
     619        1099 :   *qq = q;
     620             : 
     621             :   /* Get evec f corresponding to q */
     622        1099 :   P = redimag(primeform_u(DD, q));
     623        1099 :   he = hash_search(tbl, P);
     624        1099 :   if ( ! he) pari_err_BUG("next_prime_evec");
     625        1099 :   idx = itos((GEN) he->val);
     626        1099 :   index_to_evec(f, idx, m, k);
     627             : 
     628        1099 :   avma = av;
     629        1099 :   return 1;
     630             : }
     631             : 
     632             : /* Return 1 on success, 0 on failure. */
     633             : static int
     634         231 : orient_pcp(classgp_pcp_t G, long *ni, long D, long u, hashtable *tbl)
     635             : {
     636         231 :   pari_sp av = avma;
     637             :   /* NB: MAX_ORIENT_P = 199 seems to suffice, but can be increased if
     638             :    * necessary. */
     639             :   enum { MAX_ORIENT_P = 199 };
     640         231 :   const long *L = G->L, *n = G->n, *r = G->r, *m = G->m, *o = G->o;
     641         231 :   long i, *ps = G->orient_p, *qs = G->orient_q, *reps = G->orient_reps;
     642         231 :   long *ef, *e, *ei, *f, k = G->k, lvl = modinv_level(G->inv);
     643         231 :   GEN DD = stoi(D);
     644             : 
     645         231 :   memset(ps, 0, k * sizeof(long));
     646         231 :   memset(qs, 0, k * sizeof(long));
     647         231 :   memset(reps, 0, k * k * sizeof(long));
     648             : 
     649         280 :   for (i = 0; i < k; ++i) {
     650         280 :     ps[i] = -1;
     651         280 :     if (o[i] > 2)
     652         231 :       break;
     653             :   }
     654         329 :   for (++i; i < k; ++i)
     655          98 :     ps[i] = (o[i] > 2) ? 0 : -1; /* ps[i] = -!(o[i] > 2); */
     656             : 
     657         231 :   e = new_chunk(k);
     658         231 :   ei = new_chunk(k);
     659         231 :   f = new_chunk(k);
     660             : 
     661         595 :   for (i = 0; i < k; ++i) {
     662             :     long p;
     663         371 :     if (ps[i])
     664         630 :       continue;
     665          91 :     p = L[i];
     666          91 :     ef = &reps[i * k];
     667         665 :     while ( ! ps[i]) {
     668         504 :       if ( ! next_prime_evec(&p, ef, m, k, tbl, D, DD, u, lvl, MAX_ORIENT_P))
     669          21 :         break;
     670         483 :       evec_inverse_o(ei, ef, n, o, r, k);
     671         483 :       if ( ! distinct_inverses(NULL, ef, ei, n, o, r, k, G->L0, i))
     672         413 :         continue;
     673          70 :       ps[i] = p;
     674          70 :       qs[i] = 1;
     675             :     }
     676          91 :     if (ps[i])
     677          70 :       continue;
     678             : 
     679          21 :     p = unextprime(L[i] + 1);
     680         154 :     while ( ! ps[i]) {
     681             :       long q;
     682             : 
     683         119 :       if ( ! next_prime_evec(&p, e, m, k, tbl, D, DD, u, lvl, MAX_ORIENT_P))
     684           7 :         break;
     685         112 :       evec_inverse_o(ei, e, n, o, r, k);
     686             : 
     687         112 :       q = L[i];
     688         728 :       while ( ! qs[i]) {
     689         602 :         if ( ! next_prime_evec(&q, f, m, k, tbl, D, DD, u, lvl, p - 1))
     690          98 :           break;
     691         504 :         evec_compose(ef, e, f, n, r, k);
     692         504 :         if ( ! distinct_inverses(f, ef, ei, n, o, r, k, G->L0, i))
     693         490 :           continue;
     694             : 
     695          14 :         ps[i] = p;
     696          14 :         qs[i] = q;
     697             :       }
     698             :     }
     699          21 :     if ( ! ps[i])
     700           7 :       return 0;
     701             :   }
     702             : 
     703         224 :   if (ni) {
     704         224 :     GEN N = qfb_nform(D, *ni);
     705         224 :     hashentry *he = hash_search(tbl, N);
     706         224 :     if ( ! he) pari_err_BUG("orient_pcp");
     707         224 :     *ni = itos((GEN) he->val);
     708             :   }
     709         224 :   avma = av;
     710         224 :   return 1;
     711             : }
     712             : 
     713             : /* We must avoid situations where L_i^{+/-2} = L_j^2 (or = L_0*L_j^2
     714             :  * if ell0 flag is set), with |L_i| = |L_j| = 4 (or have 4th powers in
     715             :  * <L0> but not 2nd powers in <L0>) and j < i */
     716             : /* These cases cause problems when enumerating roots via gcds */
     717             : /* returns the index of the first bad generator, or -1 if no bad
     718             :  * generators are found */
     719             : static long
     720        3465 : classgp_pcp_check_generators(const long *n, long *r, long k, long L0)
     721             : {
     722        3465 :   pari_sp av = avma;
     723             :   long *e1, i, i0, j, s;
     724             :   const long *ei;
     725             : 
     726        3465 :   s = !!L0;
     727        3465 :   e1 = new_chunk(k);
     728             : 
     729        4337 :   for (i = s + 1; i < k; i++) {
     730         886 :     if (n[i] != 2)
     731         365 :       continue;
     732         521 :     ei = evec_ri(r, i);
     733         717 :     for (j = s; j < i; j++)
     734         577 :       if (ei[j])
     735         381 :         break;
     736         521 :     if (j == i)
     737         140 :       continue;
     738         909 :     for (i0 = s; i0 < i; i0++) {
     739         542 :       if ((4 % n[i0]))
     740         318 :         continue;
     741         224 :       evec_clear(e1, k);
     742         224 :       e1[i0] = 4;
     743         224 :       evec_reduce(e1, n, r, k);
     744         616 :       for (j = s; j < i; j++)
     745         434 :         if (e1[j])
     746          42 :           break;
     747         224 :       if (j < i)
     748          42 :         continue;       /* L_i0^4 is not trivial or in <L_0> */
     749         182 :       evec_clear(e1, k);
     750         182 :       e1[i0] = 2;
     751         182 :       evec_reduce(e1, n, r, k);   /* compute L_i0^2 */
     752         287 :       for (j = s; j < i; j++)
     753         273 :         if (e1[j] != ei[j])
     754         168 :           break;
     755         182 :       if (j == i)
     756          14 :         return i;
     757         168 :       evec_inverse(e1, e1, n, r, k);   /* compute L_i0^{-2} */
     758         259 :       for (j = s; j < i; j++)
     759         259 :         if (e1[j] != ei[j])
     760         168 :           break;
     761         168 :       if (j == i)
     762           0 :         return i;
     763             :     }
     764             :   }
     765        3451 :   avma = av;
     766        3451 :   return -1;
     767             : }
     768             : 
     769             : static void
     770        3451 : pcp_alloc_and_set(
     771             :   classgp_pcp_t G, const long *L, const long *n, const long *r, long k)
     772             : {
     773             :   /* classgp_pcp contains 6 arrays of length k (L, m, n, o, orient_p,
     774             :    * orient_q), one of length binom(k, 2) (r) and one of length k^2
     775             :    * (orient_reps). */
     776        3451 :   long rlen = k * (k - 1) / 2;
     777        3451 :   long datalen = 6 * k + rlen + k * k;
     778        3451 :   G->_data = newblock(datalen);
     779        3451 :   G->L = G->_data;
     780        3451 :   G->m = G->L + k;
     781        3451 :   G->n = G->m + k;
     782        3451 :   G->o = G->n + k;
     783        3451 :   G->r = G->o + k;
     784        3451 :   G->orient_p = G->r + rlen;
     785        3451 :   G->orient_q = G->orient_p + k;
     786        3451 :   G->orient_reps = G->orient_q + k;
     787        3451 :   G->k = k;
     788             : 
     789        3451 :   evec_copy(G->L, L, k);
     790        3451 :   evec_copy(G->n, n, k);
     791        3451 :   evec_copy(G->r, r, rlen);
     792        3451 :   evec_orders(G->o, n, r, k);
     793        3451 :   evec_n_to_m(G->m, n, k);
     794        3451 : }
     795             : 
     796             : static void
     797        3557 : classgp_pcp_clear(classgp_pcp_t G)
     798             : {
     799        3557 :   if (G->_data)
     800        3451 :     killblock(G->_data);
     801        3557 : }
     802             : 
     803             : /*
     804             :  * This is Sutherland 2009, Algorithm 2.2 (p16).
     805             :  */
     806             : static void
     807        3550 : classgp_make_pcp(
     808             :   classgp_pcp_t G, double *height, long *ni,
     809             :   long h, long D, ulong u, long inv, long Lfilter, long orient)
     810             : {
     811             :   enum { MAX_GENS = 16, MAX_RLEN = MAX_GENS * (MAX_GENS - 1) / 2 };
     812        3550 :   pari_sp av = avma, bv;
     813             :   long curr_p;
     814        3550 :   long h2, nelts, lvl = modinv_level(inv);
     815             :   GEN DD, ident, T, v;
     816             :   hashtable *tbl;
     817             :   long i, L1, L2;
     818             :   long k, L[MAX_GENS], n[MAX_GENS], r[MAX_RLEN];
     819             : 
     820        3550 :   memset(G, 0, sizeof *G);
     821             : 
     822        3550 :   G->D = D;
     823        3550 :   G->h = h;
     824        3550 :   G->inv = inv;
     825        7541 :   G->L0 = (modinv_is_double_eta(inv) && modinv_ramified(D, inv))
     826        3648 :     ? modinv_degree(NULL, NULL, inv) : 0;
     827        3550 :   G->enum_cnt = h / (1 + !!G->L0);
     828        3550 :   G->Lfilter = clcm(Lfilter, lvl);
     829             : 
     830        3550 :   if (h == 1) {
     831         106 :     if (G->L0) pari_err_BUG("classgp_pcp");
     832         106 :     G->k = 0;
     833         106 :     G->_data = NULL;
     834         106 :     v = const_vecsmall(1, 1);
     835         106 :     *height = upper_bound_on_classpoly_coeffs(D, h, v);
     836             :     /* NB: No need to set *ni when h = 1 */
     837         106 :     avma = av;
     838         106 :     return;
     839             :   }
     840             : 
     841        3444 :   DD = stoi(D);
     842        3444 :   bv = avma;
     843             :   while (1) {
     844        3465 :     k = 0;
     845             :     /* Hash table has a QFI as a key and the (boxed) index of that QFI
     846             :      * in T as its value */
     847        3465 :     tbl = hash_create(h, (ulong(*)(void*)) hash_GEN,
     848             :                          (int(*)(void*,void*))&gequal, 1);
     849        3465 :     ident = redimag(primeform_u(DD, 1));
     850        3465 :     hash_insert(tbl, ident, gen_0);
     851             : 
     852        3465 :     T = vectrunc_init(h + 1);
     853        3465 :     vectrunc_append(T, ident);
     854        3465 :     nelts = 1;
     855        3465 :     curr_p = 1;
     856             : 
     857       11561 :     while (nelts < h) {
     858             :       GEN gamma_i, beta;
     859             :       hashentry *e;
     860        4631 :       long N = glength(T), Tlen = N, ri = 1;
     861             : 
     862        4631 :       if (k == MAX_GENS) pari_err_IMPL("classgp_pcp");
     863             : 
     864        4631 :       if (nelts == 1 && G->L0) {
     865          98 :         curr_p = G->L0;
     866          98 :         gamma_i = qfb_nform(D, curr_p);
     867          98 :         beta = redimag(gamma_i);
     868         196 :         if (gequal1(gel(beta, 1))) {
     869           0 :           curr_p = 1;
     870           0 :           gamma_i = next_generator(DD, D, u, G->Lfilter, &curr_p);
     871           0 :           beta = redimag(gamma_i);
     872             :         }
     873             :       } else {
     874        4533 :         gamma_i = next_generator(DD, D, u, G->Lfilter, &curr_p);
     875        4533 :         beta = redimag(gamma_i);
     876             :       }
     877       59525 :       while ((e = hash_search(tbl, beta)) == NULL) {
     878             :         long j;
     879      115130 :         for (j = 1; j <= N; ++j) {
     880       64867 :           GEN t = qficomp(beta, gel(T, j));
     881       64867 :           vectrunc_append(T, t);
     882       64867 :           hash_insert(tbl, t, stoi(Tlen++));
     883             :         }
     884       50263 :         beta = qficomp(beta, gamma_i);
     885       50263 :         ++ri;
     886             :       }
     887        4631 :       if (ri > 1) {
     888             :         long j, si;
     889        4449 :         L[k] = curr_p;
     890        4449 :         n[k] = ri;
     891        4449 :         nelts *= ri;
     892             : 
     893             :         /* This is to reset the curr_p counter when we have G->L0 != 0
     894             :          * in the first position of L. */
     895        4449 :         if (curr_p == G->L0)
     896          98 :           curr_p = 1;
     897             : 
     898        4449 :         N = 1;
     899        4449 :         si = itos((GEN) e->val);
     900        5776 :         for (j = 0; j < k; ++j) {
     901        1327 :           evec_ri_mutate(r, k)[j] = (si / N) % n[j];
     902        1327 :           N *= n[j];
     903             :         }
     904        4449 :         ++k;
     905             :       }
     906             :     }
     907             : 
     908        3465 :     if ((i = classgp_pcp_check_generators(n, r, k, G->L0)) < 0) {
     909        3451 :       pcp_alloc_and_set(G, L, n, r, k);
     910        3451 :       if ( ! orient || orient_pcp(G, ni, D, u, tbl))
     911             :         break;
     912           7 :       G->Lfilter *= G->L[0];
     913           7 :       classgp_pcp_clear(G);
     914          14 :     } else if (log2(G->Lfilter) + log2(L[i]) >= BITS_IN_LONG)
     915           0 :       pari_err_IMPL("classgp_pcp");
     916             :     else
     917          14 :       G->Lfilter *= L[i];
     918          21 :     avma = bv;
     919          21 :   }
     920             : 
     921        3444 :   v = cgetg(h + 1, t_VECSMALL);
     922        3444 :   v[1] = 1;
     923       66540 :   for (i = 2; i <= h; ++i)
     924       63096 :     uel(v, i) = itou(gmael(T, i, 1));
     925             : 
     926        3444 :   h2 = G->L0 ? h / 2 : h;
     927        3444 :   *height = upper_bound_on_classpoly_coeffs(D, h2, v);
     928             : 
     929             :   /* The norms of the last one or two generators. */
     930        3444 :   L1 = L[k - 1];
     931        3444 :   L2 = k > 1 ? L[k - 2] : 1;
     932             :   /* 4 * L1^2 * L2^2 must fit in a ulong */
     933        3444 :   if (2 * (1 + log2(L1) + log2(L2)) >= BITS_IN_LONG)
     934           0 :     pari_err_IMPL("classgp_pcp");
     935             : 
     936        3444 :   if (G->L0 && (G->L[0] != G->L0 || G->o[0] != 2))
     937           0 :     pari_err_BUG("classgp_pcp");
     938             : 
     939        3444 :   avma = av;
     940        3444 :   return;
     941             : }
     942             : 
     943             : INLINE ulong
     944        3550 : classno_wrapper(long D)
     945             : {
     946        3550 :   pari_sp av = avma;
     947             :   GEN clsgp;
     948             :   ulong h;
     949        3550 :   clsgp = quadclassunit0(stoi(D), 0, NULL, DEFAULTPREC);
     950        3550 :   h = itou(gel(clsgp, 1));
     951        3550 :   avma = av;
     952        3550 :   return h;
     953             : }
     954             : 
     955             : 
     956             : /**
     957             :  * SECTION: Functions for calculating class polynomials.
     958             :  */
     959             : 
     960             : /* NB: Sutherland defines V_MAX to be 1200 with saying why. */
     961             : #define V_MAX 1200
     962             : 
     963             : #define NSMALL_PRIMES 11
     964             : static const long SMALL_PRIMES[11] = {
     965             :   2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31
     966             : };
     967             : 
     968             : static long
     969       44401 : is_smooth_enough(ulong *factors, long v)
     970             : {
     971             :   long i;
     972       44401 :   *factors = 0;
     973      112685 :   for (i = 0; i < NSMALL_PRIMES; ++i) {
     974      112678 :     long p = SMALL_PRIMES[i];
     975      112678 :     if (v % p == 0)
     976       53233 :       *factors |= 1UL << i;
     977      297933 :     while (v % p == 0)
     978       72577 :       v /= p;
     979      112678 :     if (v == 1)
     980       44394 :       break;
     981             :   }
     982       44401 :   return v == 1;
     983             : }
     984             : 
     985             : 
     986             : /* Hurwitz class number of |D| assuming hclassno() and attached
     987             :  * conversion to double costs much more than unegisfundamental(). */
     988             : INLINE double
     989       47944 : hclassno_wrapper(long D, long h)
     990             : {
     991             :   /* TODO: Can probably calculate hurwitz faster using -D, factor(u)
     992             :    * and classno(D). */
     993       47944 :   pari_sp av = avma;
     994       47944 :   ulong abs_D = D < 0 ? -D : D;
     995             :   double hurwitz;
     996             : 
     997       47944 :   if (h && unegisfundamental(abs_D))
     998        3375 :     hurwitz = (double) h;
     999             :   else
    1000       44569 :     hurwitz = rtodbl(gtofp(hclassno(utoi(abs_D)), DEFAULTPREC));
    1001       47944 :   avma = av;
    1002       47944 :   return hurwitz;
    1003             : }
    1004             : 
    1005             : 
    1006             : /*
    1007             :  * This is Sutherland 2009, Algorithm 2.1 (p8).
    1008             :  *
    1009             :  * NB: This function is not gerepileupto-safe.
    1010             :  */
    1011             : static GEN
    1012        3550 : select_classpoly_prime_pool(
    1013             :   double min_prime_bits, double delta, classgp_pcp_t G)
    1014             : {
    1015             :   pari_sp av;
    1016        3550 :   double prime_bits = 0.0, hurwitz, z;
    1017             :   ulong i;
    1018             :   /* t_min[v] will hold the lower bound of the t we need to look at
    1019             :    * for a given v. */
    1020             :   ulong t_min[V_MAX], t_size_lim;
    1021             :   GEN res;
    1022        3550 :   long D = G->D, inv = G->inv;
    1023             : 
    1024        3550 :   if (delta <= 0) pari_err_BUG("select_suitable_primes");
    1025        3550 :   hurwitz = hclassno_wrapper(D, G->h);
    1026             : 
    1027        3550 :   res = cgetg(1, t_VEC);
    1028             :   /* Initialise t_min to be all 2's.  This avoids trace 0 and trace
    1029             :    * 1 curves. */
    1030     4263550 :   for (i = 0; i < V_MAX; ++i)
    1031     4260000 :     t_min[i] = 2;
    1032             : 
    1033             :   /* maximum possible trace = sqrt(2^BIL - D) */
    1034        3550 :   t_size_lim = 2.0 * sqrt((double)((1UL << (BITS_IN_LONG - 2)) - (((ulong)-D) >> 2)));
    1035             : 
    1036        3550 :   av = avma;
    1037        3557 :   for (z = -D / (2.0 * hurwitz); ; z *= delta + 1.0) {
    1038             :     /* v_bound_aux = -4 z H(-D). */
    1039        3557 :     double v_bound_aux = -4.0 * z * hurwitz;
    1040             :     ulong v;
    1041        3557 :     dbg_printf(1)("z = %.2f\n", z);
    1042       44408 :     for (v = 1; ; ++v) {
    1043       44408 :       ulong pcount = 0, t, t_max, vfactors;
    1044       44408 :       ulong m_vsqr_D = v * v * (ulong)(-D);
    1045             :       /* hurwitz_ratio_bound = 11 * log(log(v + 4))^2 */
    1046       44408 :       double hurwitz_ratio_bound = log(log(v + 4.0)), max_p, H;
    1047       44408 :       hurwitz_ratio_bound *= 11.0 * hurwitz_ratio_bound;
    1048             : 
    1049       44408 :       if (v >= v_bound_aux * hurwitz_ratio_bound / D || v >= V_MAX)
    1050             :         break;
    1051             : 
    1052       44401 :       if ( ! is_smooth_enough(&vfactors, v))
    1053           7 :         continue;
    1054       44394 :       H = hclassno_wrapper(m_vsqr_D, 0);
    1055             : 
    1056             :       /* t <= 2 sqrt(p) and p <= z H(-v^2 D) and
    1057             :        *
    1058             :        *   H(-v^2 D) < vH(-D) (11 log(log(v + 4))^2)
    1059             :        *
    1060             :        * This last term is v * hurwitz * hurwitz_ratio_bound. */
    1061             : 
    1062       44394 :       max_p = z * v * hurwitz * hurwitz_ratio_bound;
    1063       44394 :       t_max = 2.0 * mindd(sqrt((double)((1UL << (BITS_IN_LONG - 2)) - (m_vsqr_D >> 2))),
    1064             :                           sqrt(max_p));
    1065    13287753 :       for (t = t_min[v]; t <= t_max; ++t) {
    1066    13243359 :         ulong possible_4p = t * t + m_vsqr_D;
    1067    13243359 :         if (possible_4p % 4 == 0) {
    1068     6622920 :           ulong possible_p = possible_4p / 4;
    1069     6622920 :           if (uisprime(possible_p) && modinv_good_prime(inv, possible_p)) {
    1070      235220 :             long p = possible_p;
    1071      235220 :             double rho_inv = p / H;
    1072             :             GEN hit;
    1073             : 
    1074      235220 :             hit = mkvecsmall5(p, t, v, (long)rho_inv, vfactors);
    1075             :             /* FIXME: Avoid doing GC for every prime as here. */
    1076      235220 :             res = gerepileupto(av, gconcat(res, hit));
    1077      235220 :             prime_bits += log2(p);
    1078      235220 :             ++pcount;
    1079             :           }
    1080             :         }
    1081             :       }
    1082       44394 :       t_min[v] = t_max + 1;
    1083             : 
    1084       44394 :       if (pcount) {
    1085       22316 :         dbg_printf(2)("  Found %lu primes for v = %lu.\n", pcount, v);
    1086       22316 :         if (gc_needed(av, 2))
    1087           0 :           res = gerepilecopy(av, res);
    1088             :       }
    1089       44394 :       if (prime_bits > min_prime_bits) {
    1090        3550 :         dbg_printf(1)("Found %ld primes; total size %.2f bits.\n",
    1091             :                     glength(res), prime_bits);
    1092        7100 :         return gerepilecopy(av, res);
    1093             :       }
    1094       40851 :     }
    1095             : 
    1096             :     /* Have we exhausted all possible solutions that fit in machine words? */
    1097           7 :     if (t_min[1] >= t_size_lim) {
    1098           0 :       char *err = stack_sprintf("class polynomial of discriminant %ld", D);
    1099           0 :       pari_err(e_ARCH, err);
    1100             :     }
    1101           7 :   }
    1102             : }
    1103             : 
    1104             : INLINE int
    1105     1031698 : cmp_small(long a, long b)
    1106             : {
    1107     1031698 :   return a>b? 1: (a<b? -1: 0);
    1108             : }
    1109             : 
    1110             : static int
    1111     1031698 : primecmp(void *data, GEN v1, GEN v2)
    1112             : {
    1113             :   (void)data;
    1114     1031698 :   return cmp_small(v1[4], v2[4]);
    1115             : }
    1116             : 
    1117             : 
    1118             : static long
    1119        3550 : height_margin(long inv, long D)
    1120             : {
    1121             :   (void)D;
    1122             :   /* NB: avs just uses a height margin of 256 for everyone and everything. */
    1123        3550 :   if (inv == INV_F)
    1124         133 :     return 64;  /* Verified for "many" discriminants up to about -350000 */
    1125        3417 :   if (inv == INV_G2)
    1126         182 :     return 5;
    1127             : 
    1128             :   /* TODO: This should be made more accurate */
    1129        3235 :   if (inv != INV_J)
    1130         462 :     return 256;
    1131             : 
    1132        2773 :   return 0;
    1133             : }
    1134             : 
    1135             : 
    1136             : static GEN
    1137        3550 : select_classpoly_primes(
    1138             :   ulong *vfactors, ulong *biggest_v,
    1139             :   long k, double delta, classgp_pcp_t G, double height)
    1140             : {
    1141        3550 :   pari_sp av = avma;
    1142        3550 :   long i, s, D = G->D, inv = G->inv;
    1143             :   ulong biggest_p;
    1144             :   double prime_bits, min_prime_bits, b;
    1145             :   GEN prime_pool;
    1146             : 
    1147             : 
    1148        3550 :   if (k < 2) pari_err_BUG("select_suitable_primes");
    1149             : 
    1150        3550 :   s = modinv_height_factor(inv);
    1151        3550 :   b = height / s + height_margin(inv, D);
    1152        3550 :   dbg_printf(1)("adjusted height = %.2f\n", b);
    1153        3550 :   min_prime_bits = k * b;
    1154             : 
    1155        3550 :   prime_pool = select_classpoly_prime_pool(min_prime_bits, delta, G);
    1156             : 
    1157             :   /* FIXME: Apply torsion constraints */
    1158             :   /* FIXME: Rank elts of res according to cost/benefit ratio */
    1159        3550 :   gen_sort_inplace(prime_pool, NULL, primecmp, NULL);
    1160             : 
    1161        3550 :   prime_bits = 0.0;
    1162        3550 :   biggest_p = gel(prime_pool, 1)[1];
    1163        3550 :   *biggest_v = gel(prime_pool, 1)[3];
    1164        3550 :   *vfactors = 0;
    1165      114784 :   for (i = 1; i < lg(prime_pool); ++i) {
    1166      114784 :     ulong p = gel(prime_pool, i)[1];
    1167      114784 :     ulong v = gel(prime_pool, i)[3];
    1168      114784 :     prime_bits += log2(p);
    1169      114784 :     *vfactors |= gel(prime_pool, i)[5];
    1170      114784 :     if (p > biggest_p)
    1171       56124 :       biggest_p = p;
    1172      114784 :     if (v > *biggest_v)
    1173       13423 :       *biggest_v = v;
    1174      114784 :     if (prime_bits > b)
    1175        3550 :       break;
    1176             :   }
    1177        3550 :   dbg_printf(1)("Selected %ld primes; largest is %lu ~ 2^%.2f\n",
    1178             :              i, biggest_p, log2(biggest_p));
    1179        3550 :   return gerepilecopy(av, vecslice0(prime_pool, 1, i));
    1180             : }
    1181             : 
    1182             : /*
    1183             :  * This is Sutherland 2009 Algorithm 1.2.
    1184             :  */
    1185             : static long
    1186      115053 : oneroot_of_classpoly(
    1187             :   ulong *j_endo, int *endo_cert, ulong j, norm_eqn_t ne, GEN jdb)
    1188             : {
    1189      115053 :   pari_sp av = avma;
    1190             :   long nfactors, L_bound, i;
    1191      115053 :   ulong p = ne->p, pi = ne->pi;
    1192             :   GEN factw, factors, u_levels, vdepths;
    1193             : 
    1194      115053 :   if (j == 0 || j == 1728 % p) pari_err_BUG("oneroot_of_classpoly");
    1195             : 
    1196      115053 :   *endo_cert = 1;
    1197      115053 :   if (ne->u * ne->v == 1) {
    1198        2506 :     *j_endo = j;
    1199        2506 :     return 1;
    1200             :   }
    1201             : 
    1202             :   /* TODO: Precalculate all this data further up */
    1203      112547 :   factw = factoru(ne->u * ne->v);
    1204      112547 :   factors = gel(factw, 1);
    1205      112547 :   nfactors = lg(factors) - 1;
    1206      112547 :   u_levels = cgetg(nfactors + 1, t_VECSMALL);
    1207      290149 :   for (i = 1; i <= nfactors; ++i)
    1208      177602 :     u_levels[i] = z_lval(ne->u, gel(factw, 1)[i]);
    1209      112547 :   vdepths = gel(factw, 2);
    1210             : 
    1211             :   /* FIXME: This should be bigger */
    1212      112547 :   L_bound = maxdd(log((double) -ne->D), (double)ne->v);
    1213             : 
    1214             :   /* Iterate over the primes L dividing w */
    1215      286576 :   for (i = 1; i <= nfactors; ++i) {
    1216      177602 :     pari_sp bv = avma;
    1217             :     GEN phi;
    1218      177602 :     long jlvl, lvl_diff, depth = vdepths[i];
    1219      177602 :     long L = factors[i];
    1220      177602 :     if (L > L_bound) {
    1221        3573 :       *endo_cert = 0;
    1222        3573 :       break;
    1223             :     }
    1224             : 
    1225      174029 :     phi = polmodular_db_getp(jdb, L, p);
    1226             : 
    1227             :     /* TODO: See if I can reuse paths created in j_level_in_volcano()
    1228             :      * later in {ascend,descend}_volcano(), perhaps by combining the
    1229             :      * functions into one "adjust_level" function. */
    1230      174029 :     jlvl = j_level_in_volcano(phi, j, p, pi, L, depth);
    1231      174029 :     lvl_diff = u_levels[i] - jlvl;
    1232             : 
    1233      174029 :     if (lvl_diff < 0) {
    1234             :       /* j's level is less than v(u) so we must ascend */
    1235      116820 :       j = ascend_volcano(phi, j, p, pi, jlvl, L, depth, -lvl_diff);
    1236       57209 :     } else if (lvl_diff > 0) {
    1237             :       /* Otherwise j's level is greater than v(u) so we descend */
    1238         678 :       j = descend_volcano(phi, j, p, pi, jlvl, L, depth, lvl_diff);
    1239             :     }
    1240      174029 :     avma = bv;
    1241             :   }
    1242      112547 :   avma = av;
    1243             :   /* At this point the probability that j has the wrong endomorphism
    1244             :    * ring is about \sum_{p|u_compl} 1/p (and u_compl must be bigger
    1245             :    * than L_bound, so pretty big), so just return it and rely on
    1246             :    * detection code in enum_j_with_endo_ring().  Detection is that we
    1247             :    * hit a previously found j-invariant earlier than expected.  OR, we
    1248             :    * evaluate class polynomials of the suborders at j and if any are
    1249             :    * zero then j must be chosen again.  */
    1250      112547 :   *j_endo = j;
    1251      112547 :   return j != 0 && j != 1728 % p;
    1252             : }
    1253             : 
    1254             : INLINE long
    1255        3304 : carray_isin(ulong *v, long n, ulong x)
    1256             : {
    1257             :   long i;
    1258      121163 :   for (i = 0; i < n; ++i)
    1259      117859 :     if (v[i] == x)
    1260           0 :       break;
    1261        3304 :   return i;
    1262             : }
    1263             : 
    1264             : INLINE ulong
    1265      228641 : select_twisting_param(ulong p)
    1266             : {
    1267             :   ulong T;
    1268             :   do
    1269      228641 :     T = random_Fl(p);
    1270      228641 :   while (krouu(T, p) != -1);
    1271      114784 :   return T;
    1272             : }
    1273             : 
    1274             : 
    1275             : INLINE void
    1276      114784 : setup_norm_eqn(norm_eqn_t ne, long D, long u, GEN norm_eqn)
    1277             : {
    1278      114784 :   ne->D = D;
    1279      114784 :   ne->u = u;
    1280      114784 :   ne->t = norm_eqn[2];
    1281      114784 :   ne->v = norm_eqn[3];
    1282      114784 :   ne->p = (ulong) norm_eqn[1];
    1283      114784 :   ne->pi = get_Fl_red(ne->p);
    1284      114784 :   ne->s2 = Fl_2gener_pre(ne->p, ne->pi);
    1285      114784 :   ne->T = select_twisting_param(ne->p);
    1286      114784 : }
    1287             : 
    1288             : INLINE ulong
    1289       13629 : Flv_powsum_pre(GEN v, ulong n, ulong p, ulong pi)
    1290             : {
    1291       13629 :   long i, l = lg(v);
    1292       13629 :   ulong psum = 0;
    1293      322217 :   for (i = 1; i < l; ++i)
    1294      308588 :     psum = Fl_add(psum, Fl_powu_pre(uel(v,i), n, p, pi), p);
    1295       13629 :   return psum;
    1296             : }
    1297             : 
    1298             : INLINE int
    1299        3550 : modinv_has_sign_ambiguity(long inv)
    1300             : {
    1301        3550 :   switch (inv) {
    1302             :   case INV_F:
    1303             :   case INV_F3:
    1304             :   case INV_W2W3E2:
    1305             :   case INV_W2W7E2:
    1306             :   case INV_W2W3:
    1307             :   case INV_W2W5:
    1308             :   case INV_W2W7:
    1309             :   case INV_W3W3:
    1310             :   case INV_W2W13:
    1311             :   case INV_W3W7:
    1312         518 :     return 1;
    1313             :   }
    1314        3032 :   return 0;
    1315             : }
    1316             : 
    1317             : INLINE int
    1318        3178 : modinv_units(int inv)
    1319             : {
    1320        3178 :   return modinv_is_double_eta(inv) || modinv_is_Weber(inv);
    1321             : }
    1322             : 
    1323             : INLINE void
    1324        9730 : adjust_signs(GEN js, ulong p, ulong pi, long inv, GEN T, long e)
    1325             : {
    1326        9730 :   long negate = 0;
    1327        9730 :   long h = lg(js) - 1;
    1328       12705 :   if ((h & 1) && modinv_units(inv)) {
    1329        2975 :     ulong prod = Flv_prod_pre(js, p, pi);
    1330        2975 :     if (prod != p - 1) {
    1331        1467 :       if (prod != 1)
    1332           0 :         pari_err_BUG("adjust_signs: constant term is not +/-1");
    1333        1467 :       negate = 1;
    1334             :     }
    1335             :   } else {
    1336             :     ulong tp, t;
    1337        6755 :     tp = umodiu(T, p);
    1338        6755 :     t = Flv_powsum_pre(js, e, p, pi);
    1339        6755 :     if (t != tp) {
    1340        3302 :       if (Fl_neg(t, p) != tp)
    1341           0 :         pari_err_BUG("adjust_signs: incorrect trace");
    1342        3302 :       negate = 1;
    1343             :     }
    1344             :   }
    1345        9730 :   if (negate)
    1346        4769 :     Flv_neg_inplace(js, p);
    1347        9730 : }
    1348             : 
    1349             : static ulong
    1350      115014 : find_jinv(
    1351             :   long *trace_tries, long *endo_tries, int *cert,
    1352             :   norm_eqn_t ne, long inv, long rho_inv, GEN jdb)
    1353             : {
    1354      115014 :   long found, ok = 1;
    1355             :   ulong j, r;
    1356             :   do {
    1357             :     do {
    1358             :       long tries;
    1359      115053 :       ulong j_t = 0;
    1360             :       /* TODO: Set batch size according to expected number of tries and
    1361             :        * experimental cost/benefit analysis. */
    1362      115053 :       tries = find_j_inv_with_given_trace(&j_t, ne, rho_inv, 0);
    1363      115053 :       if (j_t == 0)
    1364           0 :         pari_err_BUG("polclass0: Couldn't find j-invariant with given trace.");
    1365      115053 :       dbg_printf(2)("  j-invariant %ld has trace +/-%ld (%ld tries, 1/rho = %ld)\n",
    1366             :           j_t, ne->t, tries, rho_inv);
    1367      115053 :       *trace_tries += tries;
    1368             : 
    1369      115053 :       found = oneroot_of_classpoly(&j, cert, j_t, ne, jdb);
    1370      115053 :       ++*endo_tries;
    1371      115053 :     } while ( ! found);
    1372             : 
    1373      115046 :     if (modinv_is_double_eta(inv))
    1374       10452 :       ok = modfn_unambiguous_root(&r, inv, j, ne, jdb);
    1375             :     else
    1376      104594 :       r = modfn_root(j, ne, inv);
    1377      115046 :   } while ( ! ok);
    1378      115014 :   return r;
    1379             : }
    1380             : 
    1381             : 
    1382             : static GEN
    1383      114784 : polclass_roots_modp(
    1384             :   long *n_trace_curves,
    1385             :   norm_eqn_t ne, long rho_inv, classgp_pcp_t G, GEN db)
    1386             : {
    1387      114784 :   pari_sp av = avma;
    1388      114784 :   ulong j = 0;
    1389      114784 :   long inv = G->inv, h = G->h, endo_tries = 0;
    1390             :   int endo_cert;
    1391             :   GEN res, jdb, fdb;
    1392             : 
    1393      114784 :   jdb = polmodular_db_for_inv(db, INV_J);
    1394      114784 :   fdb = polmodular_db_for_inv(db, inv);
    1395             : 
    1396      114784 :   dbg_printf(2)("p = %ld, t = %ld, v = %ld\n", ne->p, ne->t, ne->v);
    1397             : 
    1398             :   do {
    1399      115014 :     j = find_jinv(n_trace_curves, &endo_tries, &endo_cert, ne, inv, rho_inv, jdb);
    1400             : 
    1401      115014 :     res = enum_roots(j, ne, fdb, G);
    1402      115014 :     if ( ! res && endo_cert) pari_err_BUG("polclass_roots_modp");
    1403      115014 :     if (res && ! endo_cert
    1404        3304 :         && carray_isin((ulong *)&res[2], h - 1, res[1]) < h - 1) {
    1405           0 :       avma = av;
    1406           0 :       res = NULL;
    1407             :     }
    1408      115014 :   } while ( ! res);
    1409             : 
    1410      114784 :   dbg_printf(2)("  j-invariant %ld has correct endomorphism ring "
    1411             :              "(%ld tries)\n", j, endo_tries);
    1412      114784 :   dbg_printf(4)("  all such j-invariants: %Ps\n", res);
    1413      114784 :   return gerepileupto(av, res);
    1414             : }
    1415             : 
    1416             : INLINE int
    1417        1869 : modinv_inverted_involution(long inv)
    1418             : {
    1419        1869 :   return modinv_is_double_eta(inv);
    1420             : }
    1421             : 
    1422             : INLINE int
    1423        1869 : modinv_negated_involution(long inv)
    1424             : {
    1425             :   /* determined by trial and error */
    1426        3738 :   return inv == INV_F || inv == INV_W3W5 || inv == INV_W3W7
    1427        3494 :     || inv == INV_W3W3 || inv == INV_W5W7;
    1428             : }
    1429             : 
    1430             : /* Return true iff Phi_L(j0, j1) = 0. */
    1431             : INLINE long
    1432        3654 : verify_edge(ulong j0, ulong j1, ulong p, ulong pi, long L, GEN fdb)
    1433             : {
    1434        3654 :   pari_sp av = avma;
    1435        3654 :   GEN phi = polmodular_db_getp(fdb, L, p);
    1436        3654 :   GEN f = Flm_Fl_polmodular_evalx(phi, L, j1, p, pi);
    1437        3654 :   ulong r = Flx_eval_pre(f, j0, p, pi);
    1438        3654 :   avma = av;
    1439        3654 :   return !r;
    1440             : }
    1441             : 
    1442             : INLINE long
    1443         672 : verify_2path(
    1444             :   ulong j1, ulong j2, ulong p, ulong pi, long L1, long L2, GEN fdb)
    1445             : {
    1446         672 :   pari_sp av = avma;
    1447         672 :   GEN phi1 = polmodular_db_getp(fdb, L1, p);
    1448         672 :   GEN phi2 = polmodular_db_getp(fdb, L2, p);
    1449         672 :   GEN f = Flm_Fl_polmodular_evalx(phi1, L1, j1, p, pi);
    1450         672 :   GEN g = Flm_Fl_polmodular_evalx(phi2, L2, j2, p, pi);
    1451         672 :   GEN d = Flx_gcd(f, g, p);
    1452         672 :   long n = degpol(d);
    1453         672 :   if (n < 2) {
    1454         672 :     avma = av;
    1455         672 :     return n;
    1456             :   }
    1457           0 :   n = Flx_nbroots(d, p);
    1458           0 :   avma = av;
    1459           0 :   return n;
    1460             : }
    1461             : 
    1462             : static long
    1463        5082 : oriented_n_action(
    1464             :   const long *ni, classgp_pcp_t G, GEN v, ulong p, ulong pi, GEN fdb)
    1465             : {
    1466        5082 :   pari_sp av = avma;
    1467        5082 :   long i, j, k = G->k;
    1468        5082 :   long nr = k * (k - 1) / 2;
    1469        5082 :   const long *n = G->n, *m = G->m, *o = G->o, *r = G->r,
    1470        5082 :     *ps = G->orient_p, *qs = G->orient_q, *reps = G->orient_reps;
    1471        5082 :   long *signs = new_chunk(k);
    1472        5082 :   long *e = new_chunk(k);
    1473        5082 :   long *rels = new_chunk(nr);
    1474             : 
    1475        5082 :   evec_copy(rels, r, nr);
    1476             : 
    1477       13664 :   for (i = 0; i < k; ++i) {
    1478             :     /* If generator doesn't require orientation, just copy its power
    1479             :      * relations and continue. */
    1480        8582 :     if (ps[i] <= 0) {
    1481        6419 :       signs[i] = 1;
    1482             :       /* power rels already copied to *rels in initialisation */
    1483        6419 :       continue;
    1484             :     }
    1485             :     /* Get rep of orientation element and express it in terms of the
    1486             :      * (partially) oriented presentation */
    1487        6195 :     for (j = 0; j < i; ++j) {
    1488        4032 :       long t = reps[i * k + j];
    1489        4032 :       e[j] = (signs[j] < 0 ? o[j] - t : t);
    1490             :     }
    1491        2163 :     e[j] = reps[i * k + j];
    1492        3318 :     for (++j; j < k; ++j)
    1493        1155 :       e[j] = 0;
    1494        2163 :     evec_reduce(e, n, rels, k);
    1495        2163 :     j = evec_to_index(e, m, k);
    1496             : 
    1497             :     /* FIXME: These calls to verify_edge recalculate powers of v[0]
    1498             :      * and v[j] over and over again, they also reduce
    1499             :      * Phi_{ps[i]} modulo p over and over again.  Need to cache
    1500             :      * these things! */
    1501        2163 :     if (qs[i] > 1) {
    1502         672 :       signs[i] =
    1503         336 :         (verify_2path(uel(v, 1), uel(v, j + 1), p, pi, ps[i], qs[i], fdb)
    1504             :             ? 1 : -1);
    1505             :     } else {
    1506             :       /* Verify ps[i]-edge to orient ith generator */
    1507        3654 :       signs[i] =
    1508        1827 :         (verify_edge(uel(v, 1), uel(v, j + 1), p, pi, ps[i], fdb)
    1509             :             ? 1 : -1);
    1510             :     }
    1511             :     /* Update power relation */
    1512        6195 :     for (j = 0; j < i; ++j) {
    1513        4032 :       long t = evec_ri(r, i)[j];
    1514        4032 :       e[j] = (signs[i] * signs[j] < 0 ? o[j] - t : t);
    1515             :     }
    1516        7644 :     while (j < k)
    1517        3318 :       e[j++] = 0;
    1518        2163 :     evec_reduce(e, n, rels, k);
    1519        6195 :     for (j = 0; j < i; ++j)
    1520        4032 :       evec_ri_mutate(rels, i)[j] = e[j];
    1521             :     /* TODO: This is a sanity check and can be removed if everything
    1522             :      * is working */
    1523        8358 :     for (j = 0; j <= i; ++j) {
    1524        6195 :       long t = reps[i * k + j];
    1525        6195 :       e[j] = (signs[j] < 0 ? o[j] - t : t);
    1526             :     }
    1527        5481 :     while (j < k)
    1528        1155 :       e[j++] = 0;
    1529        2163 :     evec_reduce(e, n, rels, k);
    1530        2163 :     j = evec_to_index(e, m, k);
    1531        2163 :     if (qs[i] > 1) {
    1532         336 :       if ( ! verify_2path(uel(v, 1), uel(v, j + 1), p, pi, ps[i], qs[i], fdb))
    1533           0 :         pari_err_BUG("oriented_n_action");
    1534             :     } else {
    1535        1827 :       if ( ! verify_edge(uel(v, 1), uel(v, j + 1), p, pi, ps[i], fdb))
    1536           0 :         pari_err_BUG("oriented_n_action");
    1537             :     }
    1538             :   }
    1539             : 
    1540             :   /* Orient representation of [N] relative to the torsor <signs, rels> */
    1541       13664 :   for (i = 0; i < k; ++i)
    1542        8582 :     e[i] = (signs[i] < 0 ? o[i] - ni[i] : ni[i]);
    1543        5082 :   evec_reduce(e, n, rels, k);
    1544        5082 :   avma = av;
    1545        5082 :   return evec_to_index(e, m, k);
    1546             : }
    1547             : 
    1548             : /* F = double_eta_raw(inv) */
    1549             : INLINE void
    1550        5082 : adjust_orientation(GEN F, long inv, GEN v, long e, ulong p, ulong pi)
    1551             : {
    1552        5082 :   ulong j0 = uel(v, 1), je = uel(v, e);
    1553             : 
    1554        5082 :   if ( ! modinv_j_from_2double_eta(F, inv, NULL, j0, je, p, pi)) {
    1555        1869 :     if (modinv_inverted_involution(inv)) {
    1556        1869 :       Flv_inv_pre_inplace(v, p, pi);
    1557             :     }
    1558        1869 :     if (modinv_negated_involution(inv))
    1559         520 :       Flv_neg_inplace(v, p);
    1560             :   }
    1561        5082 : }
    1562             : 
    1563             : 
    1564             : static void
    1565         518 : polclass_psum(
    1566             :   GEN *psum, long *d, GEN roots, GEN primes, GEN pilist, ulong h, long inv)
    1567             : {
    1568             :   /* Number of consecutive CRT stabilisations before we assume we have
    1569             :    * the correct answer. */
    1570             :   enum { MIN_STAB_CNT = 3 };
    1571         518 :   pari_sp av = avma;
    1572             :   GEN ps, psum_sqr, P;
    1573         518 :   long i, e, stabcnt, nprimes = lg(primes) - 1;
    1574             : 
    1575         518 :   if ((h & 1) && modinv_units(inv)) {
    1576         203 :     *psum = gen_1;
    1577         203 :     *d = 0;
    1578         721 :     return;
    1579             :   }
    1580             : 
    1581         315 :   e = -1;
    1582         315 :   ps = cgetg(nprimes+1, t_VECSMALL);
    1583             :   do {
    1584         364 :     e += 2;
    1585        7189 :     for (i = 1; i <= nprimes; ++i)
    1586             :     {
    1587        6874 :       GEN roots_modp = gel(roots, i);
    1588        6874 :       ulong p = uel(primes, i), pi = uel(pilist, i);
    1589        6874 :       uel(ps, i) = Flv_powsum_pre(roots_modp, e, p, pi);
    1590        6874 :       if (uel(ps, i) == 0) break;
    1591             :     }
    1592         364 :     if (i > nprimes) break;
    1593          49 :   } while (1);
    1594         315 :   psum_sqr = Z_init_CRT(0, 1);
    1595         315 :   P = gen_1;
    1596        1827 :   for (i = 1, stabcnt = 0; stabcnt < MIN_STAB_CNT && i <= nprimes; ++i)
    1597             :   {
    1598        1512 :     ulong p = uel(primes, i), pi = uel(pilist, i);
    1599        1512 :     ulong ps2 = Fl_sqr_pre(uel(ps, i), p, pi);
    1600        1512 :     ulong stab = Z_incremental_CRT(&psum_sqr, ps2, &P, p);
    1601             : 
    1602             :     /* stabcnt = stab * (stabcnt + 1) */
    1603        1512 :     if (stab)
    1604         945 :       ++stabcnt;
    1605             :     else
    1606         567 :       stabcnt = 0;
    1607        1512 :     if (gc_needed(av, 2))
    1608           0 :       gerepileall(av, 2, &psum_sqr, &P);
    1609             :   }
    1610         315 :   if (stabcnt < MIN_STAB_CNT && nprimes >= MIN_STAB_CNT)
    1611           0 :     pari_err_BUG("polclass_psum");
    1612             : 
    1613         315 :   if ( ! Z_issquareall(psum_sqr, psum)) pari_err_BUG("polclass_psum");
    1614             : 
    1615         315 :   dbg_printf(1)("Classpoly power sum (e = %ld) is %Ps; found with %.2f%% of the primes\n",
    1616           0 :       e, *psum, 100 * (i - 1) / (double) nprimes);
    1617         315 :   *psum = gerepileupto(av, *psum);
    1618         315 :   *d = e;
    1619             : }
    1620             : 
    1621             : static GEN
    1622          63 : polclass_small_disc(long D, long inv, long xvar)
    1623             : {
    1624          63 :   if (D == -3) /* x */
    1625          21 :     return pol_x(xvar);
    1626          42 :   if (D == -4) {
    1627          42 :     switch (inv) {
    1628             :     case INV_J:
    1629           7 :       return deg1pol(gen_1, stoi(-1728), xvar);
    1630             :     case INV_G2:
    1631          35 :       return deg1pol(gen_1, stoi(-12), xvar);
    1632             :     default:
    1633             :       /* There are no other invariants for which we can calculate
    1634             :        * H_{-4}(X). */
    1635           0 :       pari_err_BUG("polclass_small_disc");
    1636             :     }
    1637             :   }
    1638           0 :   return NULL;
    1639             : }
    1640             : 
    1641             : GEN
    1642        3613 : polclass0(long D, long inv, long xvar, GEN *db)
    1643             : {
    1644        3613 :   pari_sp av = avma;
    1645             :   GEN primes;
    1646        3613 :   long n_curves_tested = 0;
    1647             :   long nprimes, s, i, ni, orient;
    1648             :   GEN P, H, plist, pilist;
    1649             :   ulong u, L, maxL, vfactors, biggest_v;
    1650        3613 :   long h, p1, p2, filter = 1;
    1651             :   classgp_pcp_t G;
    1652             :   double height;
    1653             :   static const long k = 2;
    1654             :   static const double delta = 0.5;
    1655             : 
    1656        3613 :   if (D >= -4)
    1657          63 :     return polclass_small_disc(D, inv, xvar);
    1658             : 
    1659        3550 :   (void) corediscs(D, &u);
    1660        3550 :   h = classno_wrapper(D);
    1661             : 
    1662        3550 :   dbg_printf(1)("D = %ld, conductor = %ld, inv = %ld\n", D, u, inv);
    1663             : 
    1664        3550 :   ni = modinv_degree(&p1, &p2, inv);
    1665        3550 :   orient = modinv_is_double_eta(inv) && kross(D, p1) && kross(D, p2);
    1666             : 
    1667        3550 :   classgp_make_pcp(G, &height, &ni, h, D, u, inv, filter, orient);
    1668        3550 :   primes = select_classpoly_primes(&vfactors, &biggest_v, k, delta, G, height);
    1669             : 
    1670             :   /* Prepopulate *db with all the modpolys we might need */
    1671             :   /* TODO: Clean this up; in particular, note that u is factored later on. */
    1672             :   /* This comes from L_bound in oneroot_of_classpoly() above */
    1673        3550 :   maxL = maxdd(log((double) -D), (double)biggest_v);
    1674        3550 :   if (u > 1) {
    1675         847 :     for (L = 2; L <= maxL; L = unextprime(L + 1))
    1676         672 :       if ( ! (u % L))
    1677          91 :         polmodular_db_add_level(db, L, INV_J);
    1678             :   }
    1679       14613 :   for (i = 0; vfactors; ++i) {
    1680       11063 :     if (vfactors & 1UL)
    1681       10727 :       polmodular_db_add_level(db, SMALL_PRIMES[i], INV_J);
    1682       11063 :     vfactors >>= 1;
    1683             :   }
    1684        3550 :   if (p1 > 1)
    1685         441 :     polmodular_db_add_level(db, p1, INV_J);
    1686        3550 :   if (p2 > 1)
    1687         441 :     polmodular_db_add_level(db, p2, INV_J);
    1688        3550 :   s = !!G->L0;
    1689        3550 :   polmodular_db_add_levels(db, G->L + s, G->k - s, inv);
    1690        3550 :   if (orient) {
    1691         581 :     for (i = 0; i < G->k; ++i)
    1692             :     {
    1693         357 :       if (G->orient_p[i] > 1)
    1694          84 :         polmodular_db_add_level(db, G->orient_p[i], inv);
    1695         357 :       if (G->orient_q[i] > 1)
    1696          14 :         polmodular_db_add_level(db, G->orient_q[i], inv);
    1697             :     }
    1698             :   }
    1699             : 
    1700        3550 :   nprimes = lg(primes) - 1;
    1701        3550 :   H = cgetg(nprimes + 1, t_VEC);
    1702        3550 :   plist = cgetg(nprimes + 1, t_VECSMALL);
    1703        3550 :   pilist = cgetg(nprimes + 1, t_VECSMALL);
    1704      118334 :   for (i = 1; i <= nprimes; ++i) {
    1705      114784 :     long rho_inv = gel(primes, i)[4];
    1706             :     norm_eqn_t ne;
    1707      114784 :     setup_norm_eqn(ne, D, u, gel(primes, i));
    1708             : 
    1709      229568 :     gel(H, i) =
    1710      114784 :       polclass_roots_modp(&n_curves_tested, ne, rho_inv, G, *db);
    1711      114784 :     uel(plist, i) = ne->p;
    1712      114784 :     uel(pilist, i) = ne->pi;
    1713      114784 :     if (DEBUGLEVEL>2 && (i & 3L)==0) err_printf(" %ld%%", i*100/nprimes);
    1714             :   }
    1715        3550 :   dbg_printf(0)("\n");
    1716             : 
    1717        3550 :   if (orient) {
    1718         224 :     GEN nvec = new_chunk(G->k);
    1719         224 :     GEN fdb = polmodular_db_for_inv(*db, inv);
    1720         224 :     GEN F = double_eta_raw(inv);
    1721         224 :     index_to_evec((long *)nvec, ni, G->m, G->k);
    1722        5306 :     for (i = 1; i <= nprimes; ++i) {
    1723        5082 :       GEN v = gel(H, i);
    1724        5082 :       ulong p = uel(plist, i), pi = uel(pilist, i);
    1725        5082 :       long oni = oriented_n_action(nvec, G, v, p, pi, fdb);
    1726        5082 :       adjust_orientation(F, inv, v, oni + 1, p, pi);
    1727             :     }
    1728             :   }
    1729             : 
    1730        3550 :   if (modinv_has_sign_ambiguity(inv)) {
    1731             :     GEN psum;
    1732             :     long e;
    1733         518 :     polclass_psum(&psum, &e, H, plist, pilist, h, inv);
    1734       10248 :     for (i = 1; i <= nprimes; ++i) {
    1735        9730 :       GEN v = gel(H, i);
    1736        9730 :       ulong p = uel(plist, i), pi = uel(pilist, i);
    1737        9730 :       adjust_signs(v, p, pi, inv, psum, e);
    1738             :     }
    1739             :   }
    1740             : 
    1741      118334 :   for (i = 1; i <= nprimes; ++i) {
    1742      114784 :     GEN v = gel(H, i), pol;
    1743      114784 :     ulong p = uel(plist, i);
    1744             : 
    1745      114784 :     pol = Flv_roots_to_pol(v, p, xvar);
    1746      114784 :     gel(H, i) = Flx_to_Flv(pol, lg(pol) - 2);
    1747             :   }
    1748             : 
    1749        3550 :   classgp_pcp_clear(G);
    1750             : 
    1751        3550 :   dbg_printf(1)("Total number of curves tested: %ld\n", n_curves_tested);
    1752        3550 :   H = ncV_chinese_center(H, plist, &P);
    1753        3550 :   dbg_printf(1)("Result height: %.2f\n",
    1754             :              dbllog2r(itor(gsupnorm(H, DEFAULTPREC), DEFAULTPREC)));
    1755        3550 :   return gerepilecopy(av, RgV_to_RgX(H, xvar));
    1756             : }
    1757             : 
    1758             : void
    1759        1218 : check_modinv(long inv)
    1760             : {
    1761        1218 :   switch (inv) {
    1762             :   case INV_J:
    1763             :   case INV_F:
    1764             :   case INV_F2:
    1765             :   case INV_F3:
    1766             :   case INV_F4:
    1767             :   case INV_G2:
    1768             :   case INV_W2W3:
    1769             :   case INV_F8:
    1770             :   case INV_W3W3:
    1771             :   case INV_W2W5:
    1772             :   case INV_W2W7:
    1773             :   case INV_W3W5:
    1774             :   case INV_W3W7:
    1775             :   case INV_W2W3E2:
    1776             :   case INV_W2W5E2:
    1777             :   case INV_W2W13:
    1778             :   case INV_W2W7E2:
    1779             :   case INV_W3W3E2:
    1780             :   case INV_W5W7:
    1781             :   case INV_W3W13:
    1782        1204 :     break;
    1783             :   default:
    1784          14 :     pari_err_DOMAIN("polmodular", "inv", "invalid invariant", stoi(inv), gen_0);
    1785             :   }
    1786        1204 : }
    1787             : 
    1788             : GEN
    1789         616 : polclass(GEN DD, long inv, long xvar)
    1790             : {
    1791             :   GEN db, H;
    1792             :   long dummy, D;
    1793             : 
    1794         616 :   if (xvar < 0)
    1795         609 :     xvar = 0;
    1796         616 :   check_quaddisc_imag(DD, &dummy, "polclass");
    1797         609 :   check_modinv(inv);
    1798             : 
    1799         602 :   D = itos(DD);
    1800         602 :   if ( ! modinv_good_disc(inv, D))
    1801           0 :     pari_err_DOMAIN("polclass", "D", "incompatible with given invariant", stoi(inv), DD);
    1802             : 
    1803         602 :   db = polmodular_db_init(inv);
    1804         602 :   H = polclass0(D, inv, xvar, &db);
    1805         602 :   gunclone_deep(db);
    1806         602 :   return H;
    1807             : }

Generated by: LCOV version 1.11