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 - polclass.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.16.2 lcov report (development 29115-f22e516b23) Lines: 825 846 97.5 %
Date: 2024-03-29 08:06:26 Functions: 72 73 98.6 %
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_polclass
      19             : 
      20             : #define dbg_printf(lvl) if (DEBUGLEVEL >= (lvl) + 3) err_printf
      21             : 
      22             : static GEN
      23        9299 : pcp_get_L(GEN G) { return gmael(G,1,1)+1; }
      24             : static GEN
      25        6439 : pcp_get_n(GEN G) { return gmael(G,1,2)+1; }
      26             : static GEN
      27        6712 : pcp_get_m(GEN G) { return gmael(G,1,4)+1; }
      28             : static GEN
      29        6439 : pcp_get_o(GEN G) { return gmael(G,1,3)+1; }
      30             : static GEN
      31        6439 : pcp_get_r(GEN G) { return gmael(G,1,5)+1; }
      32             : static GEN
      33        6712 : pcp_get_orient_p(GEN G) { return gmael(G,3,1)+1; }
      34             : static GEN
      35        6712 : pcp_get_orient_q(GEN G) { return gmael(G,3,2)+1; }
      36             : static GEN
      37        6439 : pcp_get_orient_reps(GEN G) { return gmael(G,3,3)+1; }
      38             : static long
      39        9299 : pcp_get_L0(GEN G) { return mael(G,2,1); }
      40             : static long
      41       11300 : pcp_get_k(GEN G) { return mael(G,2,2); }
      42             : static long
      43      187788 : pcp_get_h(GEN G) { return mael(G,4,1); }
      44             : static long
      45      147085 : pcp_get_inv(GEN G) { return mael(G,4,2); }
      46             : static long
      47      146760 : pcp_get_D(GEN G) { return mael(G,4,3); }
      48             : static long
      49      187788 : pcp_get_D0(GEN G) { return mael(G,4,4); }
      50             : static long
      51      137070 : pcp_get_u(GEN G) { return mael(G,4,5); }
      52             : static GEN
      53        4861 : pcp_get_fau(GEN G) { return gel(G,5); }
      54             : 
      55             : /* SECTION: Functions dedicated to finding a j-invariant with a given
      56             :  * trace. */
      57             : 
      58             : static void
      59      138690 : hasse_bounds(long *low, long *high, long p)
      60      138690 : { long u = usqrt(4*p); *low = p + 1 - u; *high = p + 1 + u; }
      61             : 
      62             : /* a / b : a and b are from factoru and b must divide a exactly */
      63             : INLINE GEN
      64        2542 : famatsmall_divexact(GEN a, GEN b)
      65             : {
      66        2542 :   GEN a1 = gel(a,1), a2 = gel(a,2), c1, c2;
      67        2542 :   GEN b1 = gel(b,1), b2 = gel(b,2);
      68        2542 :   long i, j, k, la = lg(a1);
      69        2542 :   c1 = cgetg(la, t_VECSMALL);
      70        2542 :   c2 = cgetg(la, t_VECSMALL);
      71        8984 :   for (i = j = k = 1; j < la; j++)
      72             :   {
      73        6442 :     c1[k] = a1[j];
      74        6442 :     c2[k] = a2[j];
      75        6442 :     if (a1[j] == b1[i]) { c2[k] -= b2[i++]; if (!c2[k]) continue; }
      76        4290 :     k++;
      77             :   }
      78        2542 :   setlg(c1, k);
      79        2542 :   setlg(c2, k); return mkvec2(c1,c2);
      80             : }
      81             : 
      82             : /* This is Sutherland, 2009, TestCurveOrder.
      83             :  *
      84             :  * [a4, a6] and p specify an elliptic curve over FF_p.  N0,N1 are the two
      85             :  * possible curve orders (N0+N1 = 2p+2), and n0,n1 their factoru */
      86             : static long
      87      164729 : test_curve_order(norm_eqn_t ne, ulong a4, ulong a6,
      88             :   long N0, long N1, GEN n0, GEN n1, long hasse_low, long hasse_high)
      89             : {
      90      164729 :   pari_sp ltop = avma, av;
      91      164729 :   ulong a4t, a6t, p = ne->p, pi = ne->pi, T = ne->T;
      92             :   long m0, m1;
      93      164729 :   int swapped = 0, first = 1;
      94             : 
      95      164729 :   if (p <= 11) {
      96         414 :     long card = (long)p + 1 - Fl_elltrace(a4, a6, p);
      97         414 :     return card == N0 || card == N1;
      98             :   }
      99      164315 :   m0 = m1 = 1;
     100      164315 :   for (av = avma;;)
     101      127114 :   {
     102             :     long a1, x, n_s;
     103      291429 :     GEN Q = random_Flj_pre(a4, a6, p, pi);
     104      291421 :     if (m0 == 1)
     105      288760 :       n_s = Flj_order_ufact(Q, N0, n0, a4, p, pi);
     106        2661 :     else if (N0 % m0) n_s = 0;
     107             :     else
     108             :     { /* m0 | N0 */
     109        2542 :       GEN fa0 = famatsmall_divexact(n0, factoru(m0));
     110        2542 :       Q = Flj_mulu_pre(Q, m0, a4, p, pi);
     111        2542 :       n_s = Flj_order_ufact(Q, N0 / m0, fa0, a4, p, pi);
     112             :     }
     113      291430 :     if (n_s == 0)
     114             :     { /* If m0 divides N1 and m1 divides N0 and N0 < N1, then swap */
     115      126450 :       if (swapped || N1 % m0 || N0 % m1) return gc_long(ltop,0);
     116      100595 :       swapspec(n0, n1, N0, N1); swapped = 1; continue;
     117             :     }
     118             : 
     119      164980 :     m0 *= n_s; a1 = (2 * p + 2) % m1;
     120      323278 :     for (x = ceildivuu(hasse_low, m0) * m0; x <= hasse_high; x += m0)
     121      184819 :       if ((x % m1) == a1 && x != N0 && x != N1) break;
     122             :     /* every x was either N0 or N1, so we return true */
     123      164982 :     if (x > hasse_high) return gc_long(ltop,1);
     124             : 
     125             :     /* [a4, a6] is the given curve and [a4t, a6t] is its quadratic twist */
     126       26518 :     if (first) { Fl_elltwist_disc(a4, a6, T, p, &a4t, &a6t); first = 0; }
     127       26518 :     lswap(a4, a4t);
     128       26518 :     lswap(a6, a6t);
     129       26518 :     lswap(m0, m1); set_avma(av);
     130             :   }
     131             : }
     132             : 
     133             : static GEN
     134      854175 : random_FleV(GEN x, GEN a6, ulong p, ulong pi)
     135     3971779 : { pari_APPLY_type(t_VEC, random_Fle_pre(uel(x,i), uel(a6,i), p, pi)) }
     136             : 
     137             : /* START Code from AVSs "torcosts.h" */
     138             : struct torctab_rec {
     139             :   int m;
     140             :   int fix2, fix3;
     141             :   int N;
     142             :   int s2_flag;
     143             :   int t3_flag;
     144             :   double rating;
     145             : };
     146             : 
     147             : /* These costs assume p=2 mod 3, 3 mod 4 and not 1 mod N */
     148             : static struct torctab_rec torctab1[] = {
     149             : { 11, 1, 1, 11, 1, 1, 0.047250 },
     150             : { 33, 1, 0, 11, 1, 2, 0.047250 },
     151             : { 22, 1, 1, 11, 3, 1, 0.055125 },
     152             : { 66, 1, 0, 11, 3, 2, 0.055125 },
     153             : { 11, 1, 0, 11, 1, 0, 0.058000 },
     154             : { 13, 1, 1, 13, 1, 1, 0.058542 },
     155             : { 39, 1, 0, 13, 1, 2, 0.058542 },
     156             : { 22, 0, 1, 11, 2, 1, 0.061333 },
     157             : { 66, 0, 0, 11, 2, 2, 0.061333 },
     158             : { 22, 1, 0, 11, 3, 0, 0.061750 },
     159             : { 14, 1, 1, 14, 3, 1, 0.062500 },
     160             : { 42, 1, 0, 14, 3, 2, 0.062500 },
     161             : { 26, 1, 1, 13, 3, 1, 0.064583 },
     162             : { 78, 1, 0, 13, 3, 2, 0.064583 },
     163             : { 28, 0, 1, 14, 4, 1, 0.065625 },
     164             : { 84, 0, 0, 14, 4, 2, 0.065625 },
     165             : { 7, 1, 1, 7, 1, 1, 0.068750 },
     166             : { 13, 1, 0, 13, 1, 0, 0.068750 },
     167             : { 21, 1, 0, 7, 1, 2, 0.068750 },
     168             : { 26, 1, 0, 13, 3, 0, 0.069583 },
     169             : { 17, 1, 1, 17, 1, 1, 0.069687 },
     170             : { 51, 1, 0, 17, 1, 2, 0.069687 },
     171             : { 11, 0, 1, 11, 0, 1, 0.072500 },
     172             : { 33, 0, 0, 11, 0, 2, 0.072500 },
     173             : { 44, 1, 0, 11, 130, 0, 0.072667 },
     174             : { 52, 0, 1, 13, 4, 1, 0.073958 },
     175             : { 156, 0, 0, 13, 4, 2, 0.073958 },
     176             : { 34, 1, 1, 17, 3, 1, 0.075313 },
     177             : { 102, 1, 0, 17, 3, 2, 0.075313 },
     178             : { 15, 1, 0, 15, 1, 0, 0.075625 },
     179             : { 13, 0, 1, 13, 0, 1, 0.076667 },
     180             : { 39, 0, 0, 13, 0, 2, 0.076667 },
     181             : { 44, 0, 0, 11, 4, 0, 0.076667 },
     182             : { 30, 1, 0, 15, 3, 0, 0.077188 },
     183             : { 22, 0, 0, 11, 2, 0, 0.077333 },
     184             : { 34, 1, 0, 17, 3, 0, 0.077969 },
     185             : { 17, 1, 0, 17, 1, 0, 0.078750 },
     186             : { 14, 0, 1, 14, 0, 1, 0.080556 },
     187             : { 28, 0, 0, 14, 4, 0, 0.080556 },
     188             : { 42, 0, 0, 14, 0, 2, 0.080556 },
     189             : { 7, 1, 0, 7, 1, 0, 0.080833 },
     190             : { 9, 1, 0, 9, 1, 0, 0.080833 },
     191             : { 68, 0, 1, 17, 4, 1, 0.081380 },
     192             : { 204, 0, 0, 17, 4, 2, 0.081380 },
     193             : { 52, 0, 0, 13, 4, 0, 0.082292 },
     194             : { 10, 1, 1, 10, 3, 1, 0.084687 },
     195             : { 17, 0, 1, 17, 0, 1, 0.084687 },
     196             : { 51, 0, 0, 17, 0, 2, 0.084687 },
     197             : { 20, 0, 1, 10, 4, 1, 0.085938 },
     198             : { 60, 0, 0, 10, 4, 2, 0.085938 },
     199             : { 19, 1, 1, 19, 1, 1, 0.086111 },
     200             : { 57, 1, 0, 19, 1, 2, 0.086111 },
     201             : { 68, 0, 0, 17, 4, 0, 0.088281 },
     202             : { 38, 1, 1, 19, 3, 1, 0.089514 },
     203             : { 114, 1, 0, 19, 3, 2, 0.089514 },
     204             : { 20, 0, 0, 10, 4, 0, 0.090625 },
     205             : { 36, 0, 0, 18, 4, 0, 0.090972 },
     206             : { 26, 0, 0, 13, 2, 0, 0.091667 },
     207             : { 11, 0, 0, 11, 0, 0, 0.092000 },
     208             : { 19, 1, 0, 19, 1, 0, 0.092778 },
     209             : { 38, 1, 0, 19, 3, 0, 0.092778 },
     210             : { 14, 1, 0, 7, 3, 0, 0.092917 },
     211             : { 18, 1, 0, 9, 3, 0, 0.092917 },
     212             : { 76, 0, 1, 19, 4, 1, 0.095255 },
     213             : { 228, 0, 0, 19, 4, 2, 0.095255 },
     214             : { 10, 0, 1, 10, 0, 1, 0.096667 },
     215             : { 13, 0, 0, 13, 0, 0, 0.096667 },
     216             : { 30, 0, 0, 10, 0, 2, 0.096667 },
     217             : { 19, 0, 1, 19, 0, 1, 0.098333 },
     218             : { 57, 0, 0, 19, 0, 2, 0.098333 },
     219             : { 17, 0, 0, 17, 0, 0, 0.100000 },
     220             : { 23, 1, 1, 23, 1, 1, 0.100227 },
     221             : { 69, 1, 0, 23, 1, 2, 0.100227 },
     222             : { 7, 0, 1, 7, 0, 1, 0.100833 },
     223             : { 21, 0, 0, 7, 0, 2, 0.100833 },
     224             : { 76, 0, 0, 19, 4, 0, 0.102083 },
     225             : { 14, 0, 0, 14, 0, 0, 0.102222 },
     226             : { 18, 0, 0, 9, 2, 0, 0.102222 },
     227             : { 5, 1, 1, 5, 1, 1, 0.103125 },
     228             : { 46, 1, 1, 23, 3, 1, 0.104318 },
     229             : { 138, 1, 0, 23, 3, 2, 0.104318 },
     230             : { 23, 1, 0, 23, 1, 0, 0.105682 },
     231             : { 46, 1, 0, 23, 3, 0, 0.106705 },
     232             : { 92, 0, 1, 23, 4, 1, 0.109091 },
     233             : { 276, 0, 0, 23, 4, 2, 0.109091 },
     234             : { 19, 0, 0, 19, 0, 0, 0.110000 },
     235             : { 23, 0, 1, 23, 0, 1, 0.112273 },
     236             : { 69, 0, 0, 23, 0, 2, 0.112273 },
     237             : { 7, 0, 0, 7, 0, 0, 0.113333 },
     238             : { 9, 0, 0, 9, 0, 0, 0.113333 },
     239             : { 92, 0, 0, 23, 4, 0, 0.113826 },
     240             : { 16, 0, 1, 16, 0, 1, 0.118125 },
     241             : { 48, 0, 0, 16, 0, 2, 0.118125 },
     242             : { 5, 1, 0, 5, 1, 0, 0.121250 },
     243             : { 15, 0, 0, 15, 0, 0, 0.121250 },
     244             : { 10, 0, 0, 10, 0, 0, 0.121667 },
     245             : { 23, 0, 0, 23, 0, 0, 0.123182 },
     246             : { 12, 0, 0, 12, 0, 0, 0.141667 },
     247             : { 5, 0, 1, 5, 0, 1, 0.145000 },
     248             : { 16, 0, 0, 16, 0, 0, 0.145000 },
     249             : { 8, 0, 1, 8, 0, 1, 0.151250 },
     250             : { 29, 1, 1, 29, 1, 1, 0.153036 },
     251             : { 87, 1, 0, 29, 1, 2, 0.153036 },
     252             : { 25, 0, 0, 25, 0, 0, 0.155000 },
     253             : { 58, 1, 1, 29, 3, 1, 0.156116 },
     254             : { 174, 1, 0, 29, 3, 2, 0.156116 },
     255             : { 29, 1, 0, 29, 1, 0, 0.157500 },
     256             : { 58, 1, 0, 29, 3, 0, 0.157500 },
     257             : { 116, 0, 1, 29, 4, 1, 0.161086 },
     258             : { 29, 0, 1, 29, 0, 1, 0.163393 },
     259             : { 87, 0, 0, 29, 0, 2, 0.163393 },
     260             : { 116, 0, 0, 29, 4, 0, 0.163690 },
     261             : { 5, 0, 0, 5, 0, 0, 0.170000 },
     262             : { 8, 0, 0, 8, 0, 0, 0.170000 },
     263             : { 29, 0, 0, 29, 0, 0, 0.171071 },
     264             : { 31, 1, 1, 31, 1, 1, 0.186583 },
     265             : { 93, 1, 0, 31, 1, 2, 0.186583 },
     266             : { 62, 1, 1, 31, 3, 1, 0.189750 },
     267             : { 186, 1, 0, 31, 3, 2, 0.189750 },
     268             : { 31, 1, 0, 31, 1, 0, 0.191333 },
     269             : { 62, 1, 0, 31, 3, 0, 0.192167 },
     270             : { 124, 0, 1, 31, 4, 1, 0.193056 },
     271             : { 31, 0, 1, 31, 0, 1, 0.195333 },
     272             : { 93, 0, 0, 31, 0, 2, 0.195333 },
     273             : { 124, 0, 0, 31, 4, 0, 0.197917 },
     274             : { 2, 1, 1, 2, 3, 1, 0.200000 },
     275             : { 6, 1, 0, 2, 3, 2, 0.200000 },
     276             : { 31, 0, 0, 31, 0, 0, 0.206667 },
     277             : { 4, 1, 1, 4, 130, 1, 0.214167 },
     278             : { 6, 0, 0, 6, 0, 0, 0.226667 },
     279             : { 3, 1, 0, 3, 1, 0, 0.230000 },
     280             : { 4, 0, 1, 4, 0, 1, 0.241667 },
     281             : { 4, 1, 0, 2, 130, 0, 0.266667 },
     282             : { 4, 0, 0, 4, 0, 0, 0.283333 },
     283             : { 3, 0, 0, 3, 0, 0, 0.340000 },
     284             : { 1, 1, 1, 1, 1, 1, 0.362500 },
     285             : { 2, 0, 1, 2, 0, 1, 0.386667 },
     286             : { 1, 1, 0, 1, 1, 0, 0.410000 },
     287             : { 2, 0, 0, 2, 0, 0, 0.453333 },
     288             : };
     289             : 
     290             : static struct torctab_rec torctab2[] = {
     291             : { 11, 1, 1, 11, 1, 1, 0.047250 },
     292             : { 33, 1, 0, 11, 1, 2, 0.047250 },
     293             : { 22, 1, 1, 11, 3, 1, 0.055125 },
     294             : { 66, 1, 0, 11, 3, 2, 0.055125 },
     295             : { 13, 1, 1, 13, 1, 1, 0.057500 },
     296             : { 39, 1, 0, 13, 1, 2, 0.057500 },
     297             : { 11, 1, 0, 11, 1, 0, 0.058000 },
     298             : { 22, 0, 1, 11, 2, 1, 0.061333 },
     299             : { 66, 0, 0, 11, 2, 2, 0.061333 },
     300             : { 14, 1, 1, 14, 3, 1, 0.061458 },
     301             : { 42, 1, 0, 14, 3, 2, 0.061458 },
     302             : { 22, 1, 0, 11, 3, 0, 0.061750 },
     303             : { 26, 1, 1, 13, 3, 1, 0.064062 },
     304             : { 78, 1, 0, 13, 3, 2, 0.064062 },
     305             : { 28, 0, 1, 14, 4, 1, 0.065625 },
     306             : { 84, 0, 0, 14, 4, 2, 0.065625 },
     307             : { 13, 1, 0, 13, 1, 0, 0.066667 },
     308             : { 26, 1, 0, 13, 3, 0, 0.069583 },
     309             : { 17, 1, 1, 17, 1, 1, 0.069687 },
     310             : { 51, 1, 0, 17, 1, 2, 0.069687 },
     311             : { 11, 0, 1, 11, 0, 1, 0.070000 },
     312             : { 33, 0, 0, 11, 0, 2, 0.070000 },
     313             : { 7, 1, 1, 7, 1, 1, 0.070417 },
     314             : { 21, 1, 0, 7, 1, 2, 0.070417 },
     315             : { 15, 1, 0, 15, 1, 0, 0.072500 },
     316             : { 52, 0, 1, 13, 4, 1, 0.073090 },
     317             : { 156, 0, 0, 13, 4, 2, 0.073090 },
     318             : { 34, 1, 1, 17, 3, 1, 0.074219 },
     319             : { 102, 1, 0, 17, 3, 2, 0.074219 },
     320             : { 7, 1, 0, 7, 1, 0, 0.076667 },
     321             : { 13, 0, 1, 13, 0, 1, 0.076667 },
     322             : { 39, 0, 0, 13, 0, 2, 0.076667 },
     323             : { 44, 0, 0, 11, 4, 0, 0.076667 },
     324             : { 17, 1, 0, 17, 1, 0, 0.077188 },
     325             : { 22, 0, 0, 11, 2, 0, 0.077333 },
     326             : { 34, 1, 0, 17, 3, 0, 0.077969 },
     327             : { 30, 1, 0, 15, 3, 0, 0.080312 },
     328             : { 14, 0, 1, 14, 0, 1, 0.080556 },
     329             : { 28, 0, 0, 14, 4, 0, 0.080556 },
     330             : { 42, 0, 0, 14, 0, 2, 0.080556 },
     331             : { 9, 1, 0, 9, 1, 0, 0.080833 },
     332             : { 68, 0, 1, 17, 4, 1, 0.081380 },
     333             : { 204, 0, 0, 17, 4, 2, 0.081380 },
     334             : { 52, 0, 0, 13, 4, 0, 0.082292 },
     335             : { 10, 1, 1, 10, 3, 1, 0.083125 },
     336             : { 20, 0, 1, 10, 4, 1, 0.083333 },
     337             : { 60, 0, 0, 10, 4, 2, 0.083333 },
     338             : { 17, 0, 1, 17, 0, 1, 0.084687 },
     339             : { 51, 0, 0, 17, 0, 2, 0.084687 },
     340             : { 19, 1, 1, 19, 1, 1, 0.084722 },
     341             : { 57, 1, 0, 19, 1, 2, 0.084722 },
     342             : { 11, 0, 0, 11, 0, 0, 0.087000 },
     343             : { 68, 0, 0, 17, 4, 0, 0.088281 },
     344             : { 38, 1, 1, 19, 3, 1, 0.090139 },
     345             : { 114, 1, 0, 19, 3, 2, 0.090139 },
     346             : { 36, 0, 0, 18, 4, 0, 0.090972 },
     347             : { 19, 1, 0, 19, 1, 0, 0.091389 },
     348             : { 26, 0, 0, 13, 2, 0, 0.091667 },
     349             : { 13, 0, 0, 13, 0, 0, 0.092500 },
     350             : { 38, 1, 0, 19, 3, 0, 0.092778 },
     351             : { 14, 1, 0, 7, 3, 0, 0.092917 },
     352             : { 18, 1, 0, 9, 3, 0, 0.092917 },
     353             : { 20, 0, 0, 10, 4, 0, 0.095833 },
     354             : { 76, 0, 1, 19, 4, 1, 0.096412 },
     355             : { 228, 0, 0, 19, 4, 2, 0.096412 },
     356             : { 17, 0, 0, 17, 0, 0, 0.096875 },
     357             : { 19, 0, 1, 19, 0, 1, 0.098056 },
     358             : { 57, 0, 0, 19, 0, 2, 0.098056 },
     359             : { 23, 1, 1, 23, 1, 1, 0.100682 },
     360             : { 69, 1, 0, 23, 1, 2, 0.100682 },
     361             : { 7, 0, 1, 7, 0, 1, 0.100833 },
     362             : { 21, 0, 0, 7, 0, 2, 0.100833 },
     363             : { 30, 0, 0, 15, 2, 0, 0.100833 },
     364             : { 76, 0, 0, 19, 4, 0, 0.102083 },
     365             : { 14, 0, 0, 14, 0, 0, 0.102222 },
     366             : { 5, 1, 1, 5, 1, 1, 0.103125 },
     367             : { 46, 1, 1, 23, 3, 1, 0.104034 },
     368             : { 138, 1, 0, 23, 3, 2, 0.104034 },
     369             : { 23, 1, 0, 23, 1, 0, 0.104545 },
     370             : { 7, 0, 0, 7, 0, 0, 0.105000 },
     371             : { 10, 0, 1, 10, 0, 1, 0.105000 },
     372             : { 16, 0, 1, 16, 0, 1, 0.105417 },
     373             : { 48, 0, 0, 16, 0, 2, 0.105417 },
     374             : { 46, 1, 0, 23, 3, 0, 0.106705 },
     375             : { 18, 0, 0, 9, 2, 0, 0.107778 },
     376             : { 92, 0, 1, 23, 4, 1, 0.108239 },
     377             : { 276, 0, 0, 23, 4, 2, 0.108239 },
     378             : { 19, 0, 0, 19, 0, 0, 0.110000 },
     379             : { 23, 0, 1, 23, 0, 1, 0.111136 },
     380             : { 69, 0, 0, 23, 0, 2, 0.111136 },
     381             : { 9, 0, 0, 9, 0, 0, 0.113333 },
     382             : { 10, 0, 0, 10, 0, 0, 0.113333 },
     383             : { 92, 0, 0, 23, 4, 0, 0.113826 },
     384             : { 5, 1, 0, 5, 1, 0, 0.115000 },
     385             : { 15, 0, 0, 15, 0, 0, 0.115000 },
     386             : { 23, 0, 0, 23, 0, 0, 0.120909 },
     387             : { 8, 0, 1, 8, 0, 1, 0.126042 },
     388             : { 24, 0, 0, 8, 0, 2, 0.126042 },
     389             : { 16, 0, 0, 16, 0, 0, 0.127188 },
     390             : { 8, 0, 0, 8, 0, 0, 0.141667 },
     391             : { 25, 0, 1, 25, 0, 1, 0.144000 },
     392             : { 5, 0, 1, 5, 0, 1, 0.151250 },
     393             : { 12, 0, 0, 12, 0, 0, 0.152083 },
     394             : { 29, 1, 1, 29, 1, 1, 0.153929 },
     395             : { 87, 1, 0, 29, 1, 2, 0.153929 },
     396             : { 25, 0, 0, 25, 0, 0, 0.155000 },
     397             : { 58, 1, 1, 29, 3, 1, 0.155045 },
     398             : { 174, 1, 0, 29, 3, 2, 0.155045 },
     399             : { 29, 1, 0, 29, 1, 0, 0.156429 },
     400             : { 58, 1, 0, 29, 3, 0, 0.157857 },
     401             : { 116, 0, 1, 29, 4, 1, 0.158631 },
     402             : { 116, 0, 0, 29, 4, 0, 0.163542 },
     403             : { 29, 0, 1, 29, 0, 1, 0.164286 },
     404             : { 87, 0, 0, 29, 0, 2, 0.164286 },
     405             : { 29, 0, 0, 29, 0, 0, 0.169286 },
     406             : { 5, 0, 0, 5, 0, 0, 0.170000 },
     407             : { 31, 1, 1, 31, 1, 1, 0.187000 },
     408             : { 93, 1, 0, 31, 1, 2, 0.187000 },
     409             : { 62, 1, 1, 31, 3, 1, 0.188500 },
     410             : { 186, 1, 0, 31, 3, 2, 0.188500 },
     411             : { 31, 1, 0, 31, 1, 0, 0.191333 },
     412             : { 62, 1, 0, 31, 3, 0, 0.192083 },
     413             : { 124, 0, 1, 31, 4, 1, 0.193472 },
     414             : { 31, 0, 1, 31, 0, 1, 0.196167 },
     415             : { 93, 0, 0, 31, 0, 2, 0.196167 },
     416             : { 124, 0, 0, 31, 4, 0, 0.197083 },
     417             : { 2, 1, 1, 2, 3, 1, 0.200000 },
     418             : { 6, 1, 0, 2, 3, 2, 0.200000 },
     419             : { 31, 0, 0, 31, 0, 0, 0.205000 },
     420             : { 6, 0, 0, 6, 0, 0, 0.226667 },
     421             : { 3, 1, 0, 3, 1, 0, 0.230000 },
     422             : { 4, 0, 1, 4, 0, 1, 0.241667 },
     423             : { 4, 0, 0, 4, 0, 0, 0.283333 },
     424             : { 3, 0, 0, 3, 0, 0, 0.340000 },
     425             : { 1, 1, 1, 1, 1, 1, 0.362500 },
     426             : { 2, 0, 1, 2, 0, 1, 0.370000 },
     427             : { 1, 1, 0, 1, 1, 0, 0.385000 },
     428             : { 2, 0, 0, 2, 0, 0, 0.453333 },
     429             : };
     430             : 
     431             : static struct torctab_rec torctab3[] = {
     432             : { 66, 1, 0, 11, 3, 2, 0.040406 },
     433             : { 33, 1, 0, 11, 1, 2, 0.043688 },
     434             : { 78, 1, 0, 13, 3, 2, 0.045391 },
     435             : { 132, 1, 0, 11, 130, 2, 0.046938 },
     436             : { 39, 1, 0, 13, 1, 2, 0.047656 },
     437             : { 102, 1, 0, 17, 3, 2, 0.049922 },
     438             : { 42, 1, 0, 14, 3, 2, 0.050000 },
     439             : { 51, 1, 0, 17, 1, 2, 0.051680 },
     440             : { 132, 0, 0, 11, 4, 2, 0.052188 },
     441             : { 156, 1, 0, 13, 130, 2, 0.053958 },
     442             : { 156, 0, 0, 13, 4, 2, 0.054818 },
     443             : { 84, 1, 0, 14, 130, 2, 0.055000 },
     444             : { 15, 1, 0, 15, 1, 0, 0.056719 },
     445             : { 204, 0, 0, 17, 4, 2, 0.057227 },
     446             : { 114, 1, 0, 19, 3, 2, 0.057500 },
     447             : { 11, 1, 0, 11, 1, 0, 0.058000 },
     448             : { 66, 0, 0, 11, 2, 2, 0.058000 },
     449             : { 57, 1, 0, 19, 1, 2, 0.059062 },
     450             : { 30, 1, 0, 15, 3, 0, 0.059063 },
     451             : { 84, 0, 0, 14, 4, 2, 0.060677 },
     452             : { 22, 1, 0, 11, 3, 0, 0.061750 },
     453             : { 78, 0, 0, 13, 2, 2, 0.063542 },
     454             : { 228, 0, 0, 19, 4, 2, 0.063889 },
     455             : { 21, 1, 0, 7, 1, 2, 0.065000 },
     456             : { 138, 1, 0, 23, 3, 2, 0.065028 },
     457             : { 69, 1, 0, 23, 1, 2, 0.066903 },
     458             : { 13, 1, 0, 13, 1, 0, 0.068750 },
     459             : { 102, 0, 0, 17, 2, 2, 0.068906 },
     460             : { 26, 1, 0, 13, 3, 0, 0.069583 },
     461             : { 51, 0, 0, 17, 0, 2, 0.070312 },
     462             : { 60, 1, 0, 15, 130, 0, 0.071094 },
     463             : { 276, 0, 0, 23, 4, 2, 0.071236 },
     464             : { 39, 0, 0, 13, 0, 2, 0.071250 },
     465             : { 33, 0, 0, 11, 0, 2, 0.072750 },
     466             : { 44, 1, 0, 11, 130, 0, 0.073500 },
     467             : { 60, 0, 0, 15, 4, 0, 0.073828 },
     468             : { 9, 1, 0, 9, 1, 0, 0.074097 },
     469             : { 30, 0, 0, 15, 2, 0, 0.075625 },
     470             : { 57, 0, 0, 19, 0, 2, 0.075625 },
     471             : { 7, 1, 0, 7, 1, 0, 0.076667 },
     472             : { 44, 0, 0, 11, 4, 0, 0.076667 },
     473             : { 22, 0, 0, 11, 2, 0, 0.077333 },
     474             : { 17, 1, 0, 17, 1, 0, 0.078750 },
     475             : { 34, 1, 0, 17, 3, 0, 0.078750 },
     476             : { 69, 0, 0, 23, 0, 2, 0.079943 },
     477             : { 28, 0, 0, 14, 4, 0, 0.080556 },
     478             : { 42, 0, 0, 14, 0, 2, 0.080833 },
     479             : { 52, 0, 0, 13, 4, 0, 0.082292 },
     480             : { 14, 1, 1, 14, 3, 1, 0.083333 },
     481             : { 36, 0, 0, 18, 4, 0, 0.083391 },
     482             : { 18, 1, 0, 9, 3, 0, 0.085174 },
     483             : { 68, 0, 0, 17, 4, 0, 0.089583 },
     484             : { 15, 0, 0, 15, 0, 0, 0.090938 },
     485             : { 19, 1, 0, 19, 1, 0, 0.091389 },
     486             : { 26, 0, 0, 13, 2, 0, 0.091667 },
     487             : { 11, 0, 0, 11, 0, 0, 0.092000 },
     488             : { 13, 0, 0, 13, 0, 0, 0.092500 },
     489             : { 38, 1, 0, 19, 3, 0, 0.092778 },
     490             : { 14, 1, 0, 7, 3, 0, 0.092917 },
     491             : { 18, 0, 0, 9, 2, 0, 0.093704 },
     492             : { 174, 1, 0, 29, 3, 2, 0.095826 },
     493             : { 20, 0, 0, 10, 4, 0, 0.095833 },
     494             : { 96, 1, 0, 16, 133, 2, 0.096562 },
     495             : { 21, 0, 0, 21, 0, 0, 0.096875 },
     496             : { 87, 1, 0, 29, 1, 2, 0.096964 },
     497             : { 17, 0, 0, 17, 0, 0, 0.100000 },
     498             : { 348, 0, 0, 29, 4, 2, 0.100558 },
     499             : { 76, 0, 0, 19, 4, 0, 0.100926 },
     500             : { 14, 0, 0, 14, 0, 0, 0.102222 },
     501             : { 9, 0, 0, 9, 0, 0, 0.103889 },
     502             : { 46, 1, 0, 23, 3, 0, 0.105114 },
     503             : { 23, 1, 0, 23, 1, 0, 0.105682 },
     504             : { 48, 0, 0, 16, 0, 2, 0.106406 },
     505             : { 87, 0, 0, 29, 0, 2, 0.107545 },
     506             : { 19, 0, 0, 19, 0, 0, 0.107778 },
     507             : { 7, 0, 0, 7, 0, 0, 0.113333 },
     508             : { 10, 0, 0, 10, 0, 0, 0.113333 },
     509             : { 92, 0, 0, 23, 4, 0, 0.113636 },
     510             : { 12, 0, 0, 12, 0, 0, 0.114062 },
     511             : { 5, 1, 0, 5, 1, 0, 0.115000 },
     512             : { 186, 1, 0, 31, 3, 2, 0.115344 },
     513             : { 93, 1, 0, 31, 1, 2, 0.118125 },
     514             : { 23, 0, 0, 23, 0, 0, 0.120909 },
     515             : { 93, 0, 0, 31, 0, 2, 0.128250 },
     516             : { 16, 0, 0, 16, 0, 0, 0.138750 },
     517             : { 25, 0, 0, 25, 0, 0, 0.155000 },
     518             : { 58, 1, 0, 29, 3, 0, 0.155714 },
     519             : { 29, 1, 0, 29, 1, 0, 0.158214 },
     520             : { 3, 1, 0, 3, 1, 0, 0.163125 },
     521             : { 116, 0, 0, 29, 4, 0, 0.163690 },
     522             : { 5, 0, 0, 5, 0, 0, 0.170000 },
     523             : { 6, 0, 0, 6, 0, 0, 0.170000 },
     524             : { 8, 0, 0, 8, 0, 0, 0.170000 },
     525             : { 29, 0, 0, 29, 0, 0, 0.172857 },
     526             : { 31, 1, 0, 31, 1, 0, 0.191333 },
     527             : { 62, 1, 0, 31, 3, 0, 0.191750 },
     528             : { 124, 0, 0, 31, 4, 0, 0.197917 },
     529             : { 31, 0, 0, 31, 0, 0, 0.201667 },
     530             : { 3, 0, 0, 3, 0, 0, 0.236250 },
     531             : { 4, 0, 0, 4, 0, 0, 0.262500 },
     532             : { 2, 1, 1, 2, 3, 1, 0.317187 },
     533             : { 1, 1, 0, 1, 1, 0, 0.410000 },
     534             : { 2, 0, 0, 2, 0, 0, 0.453333 },
     535             : };
     536             : 
     537             : static struct torctab_rec torctab4[] = {
     538             : { 66, 1, 0, 11, 3, 2, 0.041344 },
     539             : { 33, 1, 0, 11, 1, 2, 0.042750 },
     540             : { 78, 1, 0, 13, 3, 2, 0.045781 },
     541             : { 39, 1, 0, 13, 1, 2, 0.046875 },
     542             : { 264, 1, 0, 11, 131, 2, 0.049043 },
     543             : { 42, 1, 0, 14, 3, 2, 0.050000 },
     544             : { 102, 1, 0, 17, 3, 2, 0.050508 },
     545             : { 51, 1, 0, 17, 1, 2, 0.051094 },
     546             : { 528, 1, 0, 11, 132, 2, 0.052891 },
     547             : { 132, 0, 0, 11, 4, 2, 0.052969 },
     548             : { 168, 1, 0, 14, 131, 2, 0.053965 },
     549             : { 156, 0, 0, 13, 4, 2, 0.054948 },
     550             : { 336, 1, 0, 14, 132, 2, 0.056120 },
     551             : { 15, 1, 0, 15, 1, 0, 0.056719 },
     552             : { 66, 0, 0, 11, 2, 2, 0.057000 },
     553             : { 114, 1, 0, 19, 3, 2, 0.057812 },
     554             : { 11, 1, 0, 11, 1, 0, 0.058000 },
     555             : { 204, 0, 0, 17, 4, 2, 0.058203 },
     556             : { 57, 1, 0, 19, 1, 2, 0.058542 },
     557             : { 84, 0, 0, 14, 4, 2, 0.059375 },
     558             : { 30, 1, 0, 15, 3, 0, 0.061406 },
     559             : { 22, 1, 0, 11, 3, 0, 0.063000 },
     560             : { 78, 0, 0, 13, 2, 2, 0.063542 },
     561             : { 138, 1, 0, 23, 3, 2, 0.064815 },
     562             : { 21, 1, 0, 7, 1, 2, 0.065000 },
     563             : { 228, 0, 0, 19, 4, 2, 0.065104 },
     564             : { 69, 1, 0, 23, 1, 2, 0.066477 },
     565             : { 13, 1, 0, 13, 1, 0, 0.068750 },
     566             : { 102, 0, 0, 17, 2, 2, 0.068906 },
     567             : { 51, 0, 0, 17, 0, 2, 0.069141 },
     568             : { 26, 1, 0, 13, 3, 0, 0.070625 },
     569             : { 276, 0, 0, 23, 4, 2, 0.071236 },
     570             : { 39, 0, 0, 13, 0, 2, 0.071250 },
     571             : { 33, 0, 0, 11, 0, 2, 0.072750 },
     572             : { 60, 0, 0, 15, 4, 0, 0.073828 },
     573             : { 9, 1, 0, 9, 1, 0, 0.074097 },
     574             : { 57, 0, 0, 19, 0, 2, 0.074583 },
     575             : { 30, 0, 0, 15, 2, 0, 0.075625 },
     576             : { 44, 0, 0, 11, 4, 0, 0.076667 },
     577             : { 17, 1, 0, 17, 1, 0, 0.077188 },
     578             : { 22, 0, 0, 11, 2, 0, 0.077333 },
     579             : { 69, 0, 0, 23, 0, 2, 0.080114 },
     580             : { 36, 0, 0, 18, 4, 0, 0.080208 },
     581             : { 34, 1, 0, 17, 3, 0, 0.080312 },
     582             : { 28, 0, 0, 14, 4, 0, 0.080556 },
     583             : { 7, 1, 0, 7, 1, 0, 0.080833 },
     584             : { 52, 0, 0, 13, 4, 0, 0.082292 },
     585             : { 42, 0, 0, 14, 0, 2, 0.082500 },
     586             : { 14, 1, 1, 14, 3, 1, 0.083333 },
     587             : { 15, 0, 0, 15, 0, 0, 0.086250 },
     588             : { 18, 1, 0, 9, 3, 0, 0.087083 },
     589             : { 26, 0, 0, 13, 2, 0, 0.088889 },
     590             : { 68, 0, 0, 17, 4, 0, 0.089583 },
     591             : { 48, 1, 0, 16, 132, 2, 0.089844 },
     592             : { 19, 1, 0, 19, 1, 0, 0.091389 },
     593             : { 11, 0, 0, 11, 0, 0, 0.092000 },
     594             : { 38, 1, 0, 19, 3, 0, 0.092917 },
     595             : { 18, 0, 0, 9, 2, 0, 0.093704 },
     596             : { 14, 1, 0, 7, 3, 0, 0.095000 },
     597             : { 96, 1, 0, 16, 133, 2, 0.095391 },
     598             : { 20, 0, 0, 10, 4, 0, 0.095833 },
     599             : { 174, 1, 0, 29, 3, 2, 0.095893 },
     600             : { 13, 0, 0, 13, 0, 0, 0.096667 },
     601             : { 17, 0, 0, 17, 0, 0, 0.096875 },
     602             : { 21, 0, 0, 21, 0, 0, 0.096875 },
     603             : { 87, 1, 0, 29, 1, 2, 0.097366 },
     604             : { 48, 0, 0, 16, 0, 2, 0.097969 },
     605             : { 24, 1, 0, 12, 131, 0, 0.098789 },
     606             : { 76, 0, 0, 19, 4, 0, 0.100926 },
     607             : { 348, 0, 0, 29, 4, 2, 0.101116 },
     608             : { 14, 0, 0, 14, 0, 0, 0.102222 },
     609             : { 9, 0, 0, 9, 0, 0, 0.103889 },
     610             : { 23, 1, 0, 23, 1, 0, 0.104545 },
     611             : { 46, 1, 0, 23, 3, 0, 0.105682 },
     612             : { 12, 0, 0, 12, 0, 0, 0.106250 },
     613             : { 87, 0, 0, 29, 0, 2, 0.108348 },
     614             : { 19, 0, 0, 19, 0, 0, 0.110000 },
     615             : { 7, 0, 0, 7, 0, 0, 0.113333 },
     616             : { 10, 0, 0, 10, 0, 0, 0.113333 },
     617             : { 92, 0, 0, 23, 4, 0, 0.113826 },
     618             : { 186, 1, 0, 31, 3, 2, 0.116094 },
     619             : { 93, 1, 0, 31, 1, 2, 0.116813 },
     620             : { 23, 0, 0, 23, 0, 0, 0.120909 },
     621             : { 5, 1, 0, 5, 1, 0, 0.121250 },
     622             : { 93, 0, 0, 31, 0, 2, 0.127625 },
     623             : { 16, 0, 0, 16, 0, 0, 0.132917 },
     624             : { 8, 0, 0, 8, 0, 0, 0.141667 },
     625             : { 25, 0, 0, 25, 0, 0, 0.152500 },
     626             : { 58, 1, 0, 29, 3, 0, 0.157946 },
     627             : { 29, 1, 0, 29, 1, 0, 0.158393 },
     628             : { 116, 0, 0, 29, 4, 0, 0.162946 },
     629             : { 3, 1, 0, 3, 1, 0, 0.163125 },
     630             : { 29, 0, 0, 29, 0, 0, 0.169286 },
     631             : { 5, 0, 0, 5, 0, 0, 0.170000 },
     632             : { 6, 0, 0, 6, 0, 0, 0.170000 },
     633             : { 31, 1, 0, 31, 1, 0, 0.191333 },
     634             : { 62, 1, 0, 31, 3, 0, 0.192083 },
     635             : { 124, 0, 0, 31, 4, 0, 0.196389 },
     636             : { 31, 0, 0, 31, 0, 0, 0.205000 },
     637             : { 3, 0, 0, 3, 0, 0, 0.255000 },
     638             : { 4, 0, 0, 4, 0, 0, 0.262500 },
     639             : { 2, 1, 1, 2, 3, 1, 0.325000 },
     640             : { 1, 1, 0, 1, 1, 0, 0.385000 },
     641             : { 2, 0, 0, 2, 0, 0, 0.420000 },
     642             : };
     643             : 
     644             : #define TWIST_DOUBLE_RATIO              (9.0/16.0)
     645             : 
     646             : static long
     647      415121 : torsion_constraint(struct torctab_rec *torctab, long ltorc, double tormod[], long n, long m)
     648             : {
     649      415121 :   long i, b = -1;
     650      415121 :   double rb = -1.;
     651    49971890 :   for (i = 0 ; i < ltorc ; i++)
     652             :   {
     653    49556769 :     struct torctab_rec *ti = torctab + i;
     654    49556769 :     if ( ! (n%ti->m) && ( !ti->fix2 || (n%(2*ti->m)) ) && ( ! ti->fix3 || (n%(3*ti->m)) ) )
     655     4496526 :       if ( n == m || ( ! (m%ti->m) && ( !ti->fix2 || (m%(2*ti->m)) ) && ( ! ti->fix3 || (m%(3*ti->m)) ) ) )
     656             :       {
     657     3618064 :         double ri = ti->rating*tormod[ti->N];
     658     3618064 :         if ( b < 0 || ri < rb ) {  b = i; rb = ri; }
     659             :       }
     660             :   }
     661      415121 :   if (b < 0) pari_err_BUG("find_rating");
     662      414501 :   return b;
     663             : }
     664             : 
     665             : /* p > 3 prime */
     666             : static void
     667      138523 : best_torsion_constraint(ulong p, long t, int *ptwist, ulong *ptor)
     668             : {
     669             :   struct torctab_rec *torctab;
     670             :   double tormod[32];
     671             :   long ltorc, n1, n2, b, b1, b2, b12, i;
     672             : 
     673      138523 :   switch(p % 12)
     674             :   {
     675       33226 :     case 11:torctab = torctab1; ltorc = numberof(torctab1); break;
     676       33225 :     case 5: torctab = torctab2; ltorc = numberof(torctab2); break;
     677       39644 :     case 7: torctab = torctab3; ltorc = numberof(torctab3); break;
     678       32428 :     default/*1*/: torctab = torctab4; ltorc = numberof(torctab4); break;
     679             :   }
     680     4574821 :   for ( i = 0 ; i < 32 ; i++ ) tormod[i] = 1.0;
     681      138523 :   if (p%5  == 1) tormod[5] = tormod[10] = tormod[15] = 6.0/5.0;
     682      138523 :   if (p%7  == 1) tormod[7] = tormod[14] = 8.0/7.0;
     683      138523 :   if (p%11 == 1) tormod[11] = 12.0/11.0;
     684      138523 :   if (p%13 == 1) tormod[13] = 14.0/13.0;
     685      138523 :   if (p%17 == 1) tormod[17] = 18.0/17.0;
     686      138523 :   if (p%19 == 1) tormod[19] = 20.0/19.0;
     687      138523 :   if (p%23 == 1) tormod[23] = 24.0/23.0;
     688      138523 :   if (p%29 == 1) tormod[29] = 30.0/29.0;
     689      138523 :   if (p%31 == 1) tormod[31] = 32.0/31.0;
     690             : 
     691      138523 :   n1 = p+1-t;
     692      138523 :   n2 = p+1+t;
     693      138523 :   b1  = torsion_constraint(torctab, ltorc, tormod, n1, n1);
     694      138453 :   b2  = torsion_constraint(torctab, ltorc, tormod, n2, n2);
     695      138523 :   b12 = torsion_constraint(torctab, ltorc, tormod, n1, n2);
     696      138775 :   if (b1 > b2) {
     697       64262 :     if (torctab[b2].rating / TWIST_DOUBLE_RATIO > torctab[b12].rating)
     698       17280 :       *ptwist = 3;
     699             :     else
     700       46982 :       *ptwist = 2;
     701             :   } else
     702       74513 :     if (torctab[b1].rating / TWIST_DOUBLE_RATIO > torctab[b12].rating)
     703       24166 :       *ptwist = 3;
     704             :     else
     705       50347 :       *ptwist = 1;
     706      138775 :   b = *ptwist ==1 ? b1: *ptwist ==2 ? b2: b12;
     707      138775 :   *ptor = torctab[b].N;
     708      138775 : }
     709             : 
     710             : /* This is Sutherland 2009 Algorithm 1.1 */
     711             : static long
     712      138712 : find_j_inv_with_given_trace(
     713             :   ulong *j_t, norm_eqn_t ne, long rho_inv, long max_curves)
     714             : {
     715      138712 :   pari_sp ltop = avma, av;
     716      138712 :   long tested = 0, t = ne->t, batch_size, N0, N1, hasse_low, hasse_high, i;
     717             :   GEN n0, n1, A4, A6, tx, ty;
     718      138712 :   ulong p = ne->p, pi = ne->pi, p1 = p + 1, a4, a6, m, N;
     719             :   int twist;
     720             : 
     721      138712 :   if (p == 2 || p == 3) {
     722          84 :     if (t == 0) pari_err_BUG("find_j_inv_with_given_trace");
     723          98 :     *j_t = t; return 1;
     724             :   }
     725             : 
     726      138656 :   N0 = (long)p1 - t; n0 = factoru(N0);
     727      138503 :   N1 = (long)p1 + t; n1 = factoru(N1);
     728      138517 :   best_torsion_constraint(p, t, &twist, &m);
     729      138756 :   switch(twist)
     730             :   {
     731       50347 :     case 1: N = N0; break;
     732       46979 :     case 2: N = N1; break;
     733       41430 :     default: N = p1;
     734             :   }
     735             : 
     736             :   /* Select batch size so that we have roughly a 50% chance of finding
     737             :    * a good curve in a batch. */
     738      138756 :   batch_size = 1.0 + rho_inv / (2.0 * m);
     739      138756 :   A4 = cgetg(batch_size + 1, t_VECSMALL);
     740      138608 :   A6 = cgetg(batch_size + 1, t_VECSMALL);
     741      138660 :   tx = cgetg(batch_size + 1, t_VECSMALL);
     742      138678 :   ty = cgetg(batch_size + 1, t_VECSMALL);
     743             : 
     744      138686 :   dbg_printf(2)("  Selected torsion constraint m = %lu and batch "
     745             :                 "size = %ld\n", m, batch_size);
     746      138686 :   hasse_bounds(&hasse_low, &hasse_high, p);
     747      138699 :   av = avma;
     748      853717 :   while (max_curves <= 0 || tested < max_curves)
     749             :   {
     750             :     GEN Pp1, Pt;
     751      853717 :     random_curves_with_m_torsion((ulong *)(A4 + 1), (ulong *)(A6 + 1),
     752      853717 :                                  (ulong *)(tx + 1), (ulong *)(ty + 1),
     753             :                                  batch_size, m, p, pi);
     754      854155 :     Pp1 = random_FleV(A4, A6, p, pi);
     755      854342 :     Pt = gcopy(Pp1);
     756      854530 :     FleV_mulu_pre_inplace(Pp1, N, A4, p, pi);
     757      853679 :     if (twist >= 3) FleV_mulu_pre_inplace(Pt, t, A4,  p, pi);
     758     3601284 :     for (i = 1; i <= batch_size; ++i) {
     759     2886291 :       ++tested;
     760     2886291 :       a4 = A4[i];
     761     2886291 :       a6 = A6[i]; if (a4 == 0 || a6 == 0) continue;
     762             : 
     763     2884298 :       if (( (twist >= 3 && mael(Pp1,i,1) == mael(Pt,i,1))
     764     2838435 :          || (twist < 3 && umael(Pp1,i,1) == p))
     765      164575 :           && test_curve_order(ne, a4, a6, N0,N1, n0,n1, hasse_low,hasse_high)) {
     766      138844 :         *j_t = Fl_ellj_pre(a4, a6, p, pi);
     767      138851 :         return gc_long(ltop, tested);
     768             :       }
     769             :     }
     770      714993 :     set_avma(av);
     771             :   }
     772           0 :   return gc_long(ltop, tested);
     773             : }
     774             : 
     775             : /* SECTION: Functions for dealing with polycyclic presentations. */
     776             : 
     777             : static GEN
     778        5655 : next_generator(GEN DD, long D, ulong u, long filter, GEN *genred, long *P)
     779             : {
     780        5655 :   pari_sp av = avma;
     781        5655 :   ulong p = (ulong)*P;
     782             :   while (1)
     783             :   {
     784        9171 :     p = unextprime(p + 1);
     785        9171 :     if (p > LONG_MAX) pari_err_BUG("next_generator");
     786        9171 :     if (kross(D, (long)p) != -1 && u % p != 0 && filter % p != 0)
     787             :     {
     788        5655 :       GEN gen = primeform_u(DD, p);
     789             :       /* If gen is in the principal class, skip it */
     790        5655 :       *genred = qfbred_i(gen);
     791        5655 :       if (!equali1(gel(*genred,1))) { *P = (long)p; return gen; }
     792           0 :       set_avma(av);
     793             :     }
     794             :   }
     795             : }
     796             : 
     797             : INLINE long *
     798        5796 : evec_ri_mutate(long r[], long i)
     799        5796 : { return r + (i * (i - 1) >> 1); }
     800             : 
     801             : INLINE const long *
     802       13523 : evec_ri(const long r[], long i)
     803       13523 : { return r + (i * (i - 1) >> 1); }
     804             : 
     805             : /* Reduces evec e so that e[i] < n[i] (assume e[i] >= 0) using pcp(n,r,k).
     806             :  * No check for overflow, this could be an issue for large groups */
     807             : INLINE void
     808       21554 : evec_reduce(long e[], const long n[], const long r[], long k)
     809             : {
     810             :   long i, j, q;
     811             :   const long *ri;
     812       21554 :   if (!k) return;
     813       51288 :   for (i = k - 1; i > 0; i--) {
     814       29734 :     if (e[i] >= n[i]) {
     815        8645 :       q = e[i] / n[i];
     816        8645 :       ri = evec_ri(r, i);
     817       21642 :       for (j = 0; j < i; j++) e[j] += q * ri[j];
     818        8645 :       e[i] -= q * n[i];
     819             :     }
     820             :   }
     821       21554 :   e[0] %= n[0];
     822             : }
     823             : 
     824             : /* Computes e3 = log(a^e1*a^e2) in terms of the given polycyclic
     825             :  * presentation (here a denotes the implicit vector of generators) */
     826             : INLINE void
     827         518 : evec_compose(long e3[],
     828             :   const long e1[], const long e2[], const long n[],const long r[], long k)
     829             : {
     830             :     long i;
     831        2128 :     for (i = 0; i < k; i++) e3[i] = e1[i] + e2[i];
     832         518 :     evec_reduce(e3, n, r, k);
     833         518 : }
     834             : 
     835             : /* Converts an evec to an integer index corresponding to the
     836             :  * multi-radix representation of the evec with moduli corresponding to
     837             :  * the subgroup orders m[i] */
     838             : INLINE long
     839       10797 : evec_to_index(const long e[], const long m[], long k)
     840             : {
     841       10797 :   long i, index = e[0];
     842       25802 :   for (i = 1; i < k; i++) index += e[i] * m[i - 1];
     843       10797 :   return index;
     844             : }
     845             : 
     846             : INLINE void
     847       12213 : evec_copy(long f[], const long e[], long k)
     848             : {
     849             :   long i;
     850       28886 :   for (i = 0; i < k; ++i) f[i] = e[i];
     851       12213 : }
     852             : 
     853             : INLINE void
     854        5042 : evec_clear(long e[], long k)
     855             : {
     856             :   long i;
     857       13000 :   for (i = 0; i < k; ++i) e[i] = 0;
     858        5042 : }
     859             : 
     860             : /* e1 and e2 may overlap */
     861             : /* Note that this function is not very efficient because it does not know the
     862             :  * orders of the elements in the presentation, only the relative orders */
     863             : INLINE void
     864         203 : evec_inverse(long e2[], const long e1[], const long n[], const long r[], long k)
     865             : {
     866         203 :   pari_sp av = avma;
     867             :   long i, *e3, *e4;
     868             : 
     869         203 :   e3 = new_chunk(k);
     870         203 :   e4 = new_chunk(k);
     871         203 :   evec_clear(e4, k);
     872         203 :   evec_copy(e3, e1, k);
     873             :   /* We have e1 + e4 = e3 which we maintain throughout while making e1
     874             :    * the zero vector */
     875         917 :   for (i = k - 1; i >= 0; i--) if (e3[i])
     876             :   {
     877          35 :     e4[i] += n[i] - e3[i];
     878          35 :     evec_reduce(e4, n, r, k);
     879          35 :     e3[i] = n[i];
     880          35 :     evec_reduce(e3, n, r, k);
     881             :   }
     882         203 :   evec_copy(e2, e4, k);
     883         203 :   set_avma(av);
     884         203 : }
     885             : 
     886             : /* e1 and e2 may overlap */
     887             : /* This is a faster way to compute inverses, if the presentation
     888             :  * element orders are known (these are specified in the array o, the
     889             :  * array n holds the relative orders) */
     890             : INLINE void
     891         735 : evec_inverse_o(
     892             :   long e2[],
     893             :   const long e1[], const long n[], const long o[], const long r[], long k)
     894             : {
     895             :   long j;
     896        3311 :   for (j = 0; j < k; j++) e2[j] = (e1[j] ? o[j] - e1[j] : 0);
     897         735 :   evec_reduce(e2, n, r, k);
     898         735 : }
     899             : 
     900             : /* Computes the order of the group element a^e using the pcp (n,r,k) */
     901             : INLINE long
     902        5529 : evec_order(const long e[], const long n[], const long r[], long k)
     903             : {
     904        5529 :   pari_sp av = avma;
     905        5529 :   long *f = new_chunk(k);
     906             :   long i, j, o, m;
     907             : 
     908        5529 :   evec_copy(f, e, k);
     909       14230 :   for (o = 1, i = k - 1; i >= 0; i--) if (f[i])
     910             :   {
     911        6485 :     m = n[i] / ugcd(f[i], n[i]);
     912       17784 :     for (j = 0; j < k; j++) f[j] *= m;
     913        6485 :     evec_reduce(f, n, r, k);
     914        6485 :     o *= m;
     915             :   }
     916        5529 :   return gc_long(av,o);
     917             : }
     918             : 
     919             : /* Computes orders o[] for each generator using relative orders n[]
     920             :  * and power relations r[] */
     921             : INLINE void
     922        4328 : evec_orders(long o[], const long n[], const long r[], long k)
     923             : {
     924        4328 :   pari_sp av = avma;
     925        4328 :   long i, *e = new_chunk(k);
     926             : 
     927        4328 :   evec_clear(e, k);
     928        9857 :   for (i = 0; i < k; i++) {
     929        5529 :     e[i] = 1;
     930        5529 :     if (i) e[i - 1] = 0;
     931        5529 :     o[i] = evec_order(e, n, r, k);
     932             :   }
     933        4328 :   set_avma(av);
     934        4328 : }
     935             : 
     936             : INLINE int
     937         259 : evec_equal(const long e1[], const long e2[], long k)
     938             : {
     939             :   long j;
     940         371 :   for (j = 0; j < k; ++j)
     941         343 :     if (e1[j] != e2[j]) break;
     942         259 :   return j == k;
     943             : }
     944             : 
     945             : INLINE void
     946        1386 : index_to_evec(long e[], long index, const long m[], long k)
     947             : {
     948             :   long i;
     949        4158 :   for (i = k - 1; i > 0; --i) {
     950        2772 :     e[i] = index / m[i - 1];
     951        2772 :     index -= e[i] * m[i - 1];
     952             :   }
     953        1386 :   e[0] = index;
     954        1386 : }
     955             : 
     956             : INLINE void
     957        4328 : evec_n_to_m(long m[], const long n[], long k)
     958             : {
     959             :   long i;
     960        4328 :   m[0] = n[0];
     961        5529 :   for (i = 1; i < k; ++i) m[i] = m[i - 1] * n[i];
     962        4328 : }
     963             : 
     964             : /* Based on logfac() in Sutherland's classpoly package.
     965             :  * Ramanujan approximation to log(n!), accurate to O(1/n^3) */
     966             : INLINE double
     967           0 : logfac(long n)
     968             : {
     969           0 :   const double HALFLOGPI = 0.57236494292470008707171367567653;
     970           0 :   return n * log((double) n) - (double) n +
     971           0 :     log((double) n * (1.0 + 4.0 * n * (1.0 + 2.0 * n))) / 6.0 + HALFLOGPI;
     972             : }
     973             : 
     974             : /* This is based on Sutherland 2009, Lemma 8 (p31). */
     975             : static double
     976        4861 : upper_bound_on_classpoly_coeffs(long D, long h, GEN qfinorms)
     977             : {
     978        4861 :   double B, logbinom, lnMk, lnMh = 0, C = 2114.567, t = M_PI * sqrt((double)-D);
     979        4861 :   ulong maxak = 0;
     980             :   long k, m;
     981             : 
     982       80185 :   for (k = 1, B = 0.0; k <= h; ++k)
     983             :   {
     984       75324 :     ulong ak = uel(qfinorms, k);
     985       75324 :     double tk = t / ak;
     986       75324 :     lnMk = tk + log(1.0 + C * exp(-tk));
     987       75324 :     B += lnMk;
     988       75324 :     if (ak > maxak) { maxak = ak; lnMh = lnMk; }
     989             :   }
     990        4861 :   m = floor((h + 1)/(exp(lnMh) + 1.0));
     991             :   /* log(binom(h, m)); 0 unless D <= -1579751 */
     992        4861 :   logbinom = (m > 0 && m < h)? logfac(h) - logfac(m) - logfac(h - m): 0;
     993        4861 :   return (B + logbinom - m * lnMh) * (1 / M_LN2) + 2.0;
     994             : }
     995             : 
     996             : INLINE long
     997        1001 : distinct_inverses(const long f[], const long ef[], const long ei[],
     998             :   const long n[], const long o[], const long r[], long k, long L0, long i)
     999             : {
    1000        1001 :   pari_sp av = avma;
    1001             :   long j, *e2, *e3;
    1002             : 
    1003        1001 :   if ( ! ef[i] || (L0 && ef[0])) return 0;
    1004         574 :   for (j = i + 1; j < k; ++j)
    1005         455 :     if (ef[j]) break;
    1006         490 :   if (j < k) return 0;
    1007             : 
    1008         119 :   e2 = new_chunk(k);
    1009         119 :   evec_copy(e2, ef, i);
    1010         119 :   e2[i] = o[i] - ef[i];
    1011         175 :   for (j = i + 1; j < k; ++j) e2[j] = 0;
    1012         119 :   evec_reduce(e2, n, r, k);
    1013             : 
    1014         119 :   if (evec_equal(ef, e2, k)) return gc_long(av,0);
    1015             : 
    1016         112 :   e3 = new_chunk(k);
    1017         112 :   evec_inverse_o(e3, ef, n, o, r, k);
    1018         112 :   if (evec_equal(e2, e3, k)) return gc_long(av,0);
    1019             : 
    1020          91 :   if (f) {
    1021          14 :     evec_compose(e3, f, ei, n, r, k);
    1022          14 :     if (evec_equal(e2, e3, k)) return gc_long(av,0);
    1023             : 
    1024          14 :     evec_inverse_o(e3, e3, n, o, r, k);
    1025          14 :     if (evec_equal(e2, e3, k)) return gc_long(av,0);
    1026             :   }
    1027          91 :   return gc_long(av,1);
    1028             : }
    1029             : 
    1030             : INLINE long
    1031        1239 : next_prime_evec(long *qq, long f[], const long m[], long k,
    1032             :   hashtable *tbl, long D, GEN DD, long u, long lvl, long ubound)
    1033             : {
    1034        1239 :   pari_sp av = avma;
    1035             :   hashentry *he;
    1036             :   GEN P;
    1037        1239 :   long idx, q = *qq;
    1038             : 
    1039        3164 :   do q = unextprime(q + 1);
    1040        3164 :   while (!(u % q) || kross(D, q) == -1 || !(lvl % q) || !(D % (q * q)));
    1041        1239 :   if (q > ubound) return 0;
    1042        1113 :   *qq = q;
    1043             : 
    1044             :   /* Get evec f corresponding to q */
    1045        1113 :   P = qfbred_i(primeform_u(DD, q));
    1046        1113 :   he = hash_search(tbl, P);
    1047        1113 :   if (!he) pari_err_BUG("next_prime_evec");
    1048        1113 :   idx = (long)he->val;
    1049        1113 :   index_to_evec(f, idx, m, k);
    1050        1113 :   return gc_long(av,1);
    1051             : }
    1052             : 
    1053             : /* Return 1 on success, 0 on failure. */
    1054             : static int
    1055         280 : orient_pcp(GEN G, long *ni, long D, long u, hashtable *tbl)
    1056             : {
    1057         280 :   pari_sp av = avma;
    1058             :   /* 199 seems to suffice, but can be increased if necessary */
    1059             :   enum { MAX_ORIENT_P = 199 };
    1060         280 :   const long *L = pcp_get_L(G), *n = pcp_get_n(G), *r = pcp_get_r(G), *m = pcp_get_m(G), *o = pcp_get_o(G);
    1061         280 :   long i, *ps = pcp_get_orient_p(G), *qs = pcp_get_orient_q(G), *reps = pcp_get_orient_reps(G);
    1062         280 :   long *ef, *e, *ei, *f, k = pcp_get_k(G), lvl = modinv_level(pcp_get_inv(G));
    1063         280 :   GEN DD = stoi(D);
    1064         280 :   long L0 = pcp_get_L0(G);
    1065             : 
    1066         280 :   memset(ps, 0, k * sizeof(long));
    1067         280 :   memset(qs, 0, k * sizeof(long));
    1068         280 :   memset(reps, 0, k * k * sizeof(long));
    1069             : 
    1070         357 :   for (i = 0; i < k; ++i) { ps[i] = -1; if (o[i] > 2) break; }
    1071         385 :   for (++i; i < k; ++i) ps[i] = (o[i] > 2) ? 0 : -1; /* ps[i] = -!(o[i] > 2); */
    1072             : 
    1073         280 :   e = new_chunk(k);
    1074         280 :   ei = new_chunk(k);
    1075         280 :   f = new_chunk(k);
    1076             : 
    1077         728 :   for (i = 0; i < k; ++i) {
    1078             :     long p;
    1079         532 :     if (ps[i]) continue;
    1080          98 :     p = L[i];
    1081          98 :     ef = &reps[i * k];
    1082         595 :     while (!ps[i]) {
    1083         518 :       if (!next_prime_evec(&p, ef, m, k, tbl, D, DD, u, lvl, MAX_ORIENT_P))
    1084          21 :         break;
    1085         497 :       evec_inverse_o(ei, ef, n, o, r, k);
    1086         497 :       if (!distinct_inverses(NULL, ef, ei, n, o, r, k, L0, i)) continue;
    1087          77 :       ps[i] = p;
    1088          77 :       qs[i] = 1;
    1089             :     }
    1090          98 :     if (ps[i]) continue;
    1091             : 
    1092          21 :     p = unextprime(L[i] + 1);
    1093         133 :     while (!ps[i]) {
    1094             :       long q;
    1095             : 
    1096         119 :       if (!next_prime_evec(&p, e, m, k, tbl, D, DD, u, lvl, MAX_ORIENT_P))
    1097           7 :         break;
    1098         112 :       evec_inverse_o(ei, e, n, o, r, k);
    1099             : 
    1100         112 :       q = L[i];
    1101         616 :       while (!qs[i]) {
    1102         602 :         if (!next_prime_evec(&q, f, m, k, tbl, D, DD, u, lvl, p - 1)) break;
    1103         504 :         evec_compose(ef, e, f, n, r, k);
    1104         504 :         if (!distinct_inverses(f, ef, ei, n, o, r, k, L0, i)) continue;
    1105          14 :         ps[i] = p;
    1106          14 :         qs[i] = q;
    1107             :       }
    1108             :     }
    1109          21 :     if (!ps[i]) return 0;
    1110             :   }
    1111         273 :   if (ni) {
    1112         273 :     GEN N = qfb_nform(D, *ni);
    1113         273 :     hashentry *he = hash_search(tbl, N);
    1114         273 :     if (!he) pari_err_BUG("orient_pcp");
    1115         273 :     *ni = (long)he->val;
    1116             :   }
    1117         273 :   return gc_bool(av,1);
    1118             : }
    1119             : 
    1120             : /* We must avoid situations where L_i^{+/-2} = L_j^2 (or = L_0*L_j^2
    1121             :  * if ell0 flag is set), with |L_i| = |L_j| = 4 (or have 4th powers in
    1122             :  * <L0> but not 2nd powers in <L0>) and j < i */
    1123             : /* These cases cause problems when enumerating roots via gcds */
    1124             : /* returns the index of the first bad generator, or -1 if no bad
    1125             :  * generators are found */
    1126             : static long
    1127        4349 : classgp_pcp_check_generators(const long *n, long *r, long k, long L0)
    1128             : {
    1129        4349 :   pari_sp av = avma;
    1130             :   long *e1, i, i0, j, s;
    1131             :   const long *ei;
    1132             : 
    1133        4349 :   s = !!L0;
    1134        4349 :   e1 = new_chunk(k);
    1135             : 
    1136        5445 :   for (i = s + 1; i < k; i++) {
    1137        1117 :     if (n[i] != 2) continue;
    1138         689 :     ei = evec_ri(r, i);
    1139         983 :     for (j = s; j < i; j++)
    1140         752 :       if (ei[j]) break;
    1141         689 :     if (j == i) continue;
    1142        1077 :     for (i0 = s; i0 < i; i0++) {
    1143         640 :       if ((4 % n[i0])) continue;
    1144         287 :       evec_clear(e1, k);
    1145         287 :       e1[i0] = 4;
    1146         287 :       evec_reduce(e1, n, r, k);
    1147         735 :       for (j = s; j < i; j++)
    1148         511 :         if (e1[j]) break;
    1149         287 :       if (j < i) continue; /* L_i0^4 is not trivial or in <L_0> */
    1150         224 :       evec_clear(e1, k);
    1151         224 :       e1[i0] = 2;
    1152         224 :       evec_reduce(e1, n, r, k); /* compute L_i0^2 */
    1153         357 :       for (j = s; j < i; j++)
    1154         336 :         if (e1[j] != ei[j]) break;
    1155         224 :       if (j == i) return i;
    1156         203 :       evec_inverse(e1, e1, n, r, k); /* compute L_i0^{-2} */
    1157         315 :       for (j = s; j < i; j++)
    1158         315 :         if (e1[j] != ei[j]) break;
    1159         203 :       if (j == i) return i;
    1160             :     }
    1161             :   }
    1162        4328 :   return gc_long(av,-1);
    1163             : }
    1164             : 
    1165             : /* This is Sutherland 2009, Algorithm 2.2 (p16). */
    1166             : static GEN
    1167        4861 : classgp_make_pcp(double *height, long *ni, long h, long D, long D0, ulong u,
    1168             :   GEN Pu, GEN Eu, long inv, long Lfilter, long orient)
    1169             : {
    1170        4861 :   const long MAX_GENS = 16, lvl = modinv_level(inv);
    1171             :   pari_sp av, av2;
    1172             :   long curr_p, nelts, L0, enum_cnt, GLfilter, i, k, L1, L2;
    1173             :   GEN G, DD, ident, T, v, L, m, n, o, r, orient_p, orient_q, orient_reps;
    1174             :   hashtable *tbl;
    1175             : 
    1176        4861 :   L = zero_zv(MAX_GENS);
    1177        4861 :   m = zero_zv(MAX_GENS);
    1178        4861 :   n = zero_zv(MAX_GENS);
    1179        4861 :   o = zero_zv(MAX_GENS);
    1180        4861 :   r = zero_zv(MAX_GENS * (MAX_GENS-1) / 2);
    1181        4861 :   orient_p = zero_zv(MAX_GENS);
    1182        4861 :   orient_q = zero_zv(MAX_GENS);
    1183        4861 :   orient_reps = zero_zv(MAX_GENS*MAX_GENS);
    1184        4861 :   G = mkvec5(mkvec5(L, n, o, m, r), mkvecsmall3(0, 0, 0),
    1185             :              mkvec3(orient_p, orient_q, orient_reps),
    1186             :              mkvecsmall5(h, inv, D, D0, u), mkmat2(Pu, Eu));
    1187        4861 :   av = avma;
    1188        4861 :   if (h == 1)
    1189             :   {
    1190         540 :     *height = upper_bound_on_classpoly_coeffs(D, h, mkvecsmall(1));
    1191         540 :     return gc_const(av, G); /* no need to set *ni when h = 1 */
    1192             :   }
    1193        4321 :   if (!modinv_is_double_eta(inv) || !modinv_ramified(D, inv, &L0)) L0 = 0;
    1194        4321 :   enum_cnt = h / (1 + !!L0);
    1195        4321 :   GLfilter = ulcm(Lfilter, lvl);
    1196        4321 :   DD = stoi(D); av2 = avma;
    1197             :   while (1) {
    1198        4349 :     k = 0;
    1199             :     /* Hash table has an imaginary QFB as a key and the index of that
    1200             :      * form in T as its value */
    1201        4349 :     tbl = hash_create(h, (ulong(*)(void*)) hash_GEN,
    1202             :                          (int(*)(void*,void*))&gidentical, 1);
    1203        4349 :     ident = primeform(DD, gen_1);
    1204        4349 :     hash_insert(tbl, ident, (void*)0);
    1205             : 
    1206        4349 :     T = vectrunc_init(h + 1);
    1207        4349 :     vectrunc_append(T, ident);
    1208        4349 :     nelts = 1;
    1209        4349 :     curr_p = 1;
    1210             : 
    1211       10109 :     while (nelts < h) {
    1212        5760 :       GEN gamma_i = NULL, beta;
    1213             :       hashentry *e;
    1214        5760 :       long N = lg(T)-1, ri = 1;
    1215             : 
    1216        5760 :       if (k == MAX_GENS) pari_err_IMPL("classgp_pcp");
    1217             : 
    1218        5760 :       if (nelts == 1 && L0) {
    1219         105 :         curr_p = L0;
    1220         105 :         gamma_i = qfb_nform(D, curr_p);
    1221         105 :         beta = qfbred_i(gamma_i);
    1222         105 :         if (equali1(gel(beta, 1))) { curr_p = 1; gamma_i  = NULL; }
    1223             :       }
    1224        5760 :       if (!gamma_i)
    1225        5655 :         gamma_i = next_generator(DD, D, u, GLfilter, &beta, &curr_p);
    1226       62459 :       while ((e = hash_search(tbl, beta)) == NULL) {
    1227             :         long j;
    1228      130564 :         for (j = 1; j <= N; ++j) {
    1229       73865 :           GEN t = qfbcomp_i(beta, gel(T, j));
    1230       73865 :           vectrunc_append(T, t);
    1231       73865 :           hash_insert(tbl, t, (void*)tbl->nb);
    1232             :         }
    1233       56699 :         beta = qfbcomp_i(beta, gamma_i);
    1234       56699 :         ++ri;
    1235             :       }
    1236        5760 :       if (ri > 1) {
    1237             :         long j, si;
    1238        5571 :         L[k+1] = curr_p;
    1239        5571 :         n[k+1] = ri;
    1240        5571 :         nelts *= ri;
    1241             : 
    1242             :         /* This is to reset the curr_p counter when we have L0 != 0
    1243             :          * in the first position of L. */
    1244        5571 :         if (curr_p == L0) curr_p = 1;
    1245             : 
    1246        5571 :         N = 1;
    1247        5571 :         si = (long)e->val;
    1248        7178 :         for (j = 1; j <= k; ++j) {
    1249        1607 :           evec_ri_mutate(r, k)[j] = (si / N) % n[j];
    1250        1607 :           N *= n[j];
    1251             :         }
    1252        5571 :         ++k;
    1253             :       }
    1254             :     }
    1255             : 
    1256        4349 :     if ((i = classgp_pcp_check_generators(n+1, r+1, k, L0)) < 0) {
    1257        4328 :       mael(G,2,1) = L0;
    1258        4328 :       mael(G,2,2) = k;
    1259        4328 :       mael(G,2,3) = enum_cnt;
    1260        4328 :       evec_orders(o+1, n+1, r+1, k);
    1261        4328 :       evec_n_to_m(m+1, n+1, k);
    1262        4328 :       if (!orient || orient_pcp(G, ni, D, u, tbl)) break;
    1263           7 :       GLfilter *= L[1];
    1264             :     } else {
    1265          21 :       GLfilter = umuluu_or_0(GLfilter, L[i+1]);
    1266          21 :       if (!GLfilter) pari_err_IMPL("classgp_pcp");
    1267             :     }
    1268          28 :     set_avma(av2);
    1269             :   }
    1270        4321 :   v = cgetg(h + 1, t_VECSMALL);
    1271        4321 :   v[1] = 1;
    1272       76366 :   for (i = 2; i <= h; ++i) uel(v,i) = itou(gmael(T,i,1));
    1273        4321 :   *height = upper_bound_on_classpoly_coeffs(D, enum_cnt, v);
    1274             : 
    1275             :   /* The norms of the last one or two generators. */
    1276        4321 :   L1 = L[k];
    1277        4321 :   L2 = k > 1 ? L[k - 1] : 1;
    1278             :   /* 4 * L1^2 * L2^2 must fit in a ulong */
    1279        4321 :   if (2 * (1 + log2(L1) + log2(L2)) >= BITS_IN_LONG)
    1280           0 :     pari_err_IMPL("classgp_pcp");
    1281        4321 :   if (L0 && (L[1] != L0 || o[1] != 2)) pari_err_BUG("classgp_pcp");
    1282        4321 :   return gc_const(av, G);
    1283             : }
    1284             : 
    1285             : /* SECTION: Functions for calculating class polynomials. */
    1286             : static const long SMALL_PRIMES[11] = { 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31 };
    1287             : 
    1288             : /* Encodes prime divisors of smooth integers <= 1200
    1289             :   P = primes(11); V = vector(31, i, vecsearch(P, i));
    1290             :   vfactor(v) =
    1291             :   { if (v == 1, return (0));
    1292             :     my(q = factor(v)[,1]);
    1293             :     if (vecmax(q) > 31, return (-1)); \\ not smooth
    1294             :     sum(i = 1, #P, 1 << (V[q[i]]-1)); }
    1295             :   vector(1200, v, vfactor(v)) */
    1296             : static const long
    1297             : SMOOTH_INTS[] = { -1, /* 0 */
    1298             : 0,1,2,1,4,3,8,1,2,5,16,3,32,9,6,1,64,3,128,5,10,17,256,3,4,33,2,9,512,7,1024,1,
    1299             : 18,65,12,3,-1,129,34,5,-1,11,-1,17,6,257,-1,3,8,5,66,33,-1,3,20,9,130,513,-1,7,
    1300             : -1,1025,10,1,36,19,-1,65,258,13,-1,3,-1,-1,6,129,24,35,-1,5,2,-1,-1,11,68,-1,
    1301             : 514,17,-1,7,40,257,1026,-1,132,3,-1,9,18,5,-1,67,-1,33,14,-1,-1,3,-1,21,-1,9,-1,
    1302             : 131,260,513,34,-1,72,7,16,-1,-1,1025,4,11,-1,1,-1,37,-1,19,136,-1,6,65,-1,259,
    1303             : -1,13,-1,-1,48,3,516,-1,10,-1,-1,7,-1,129,66,25,1028,35,-1,-1,-1,5,264,3,-1,-1,
    1304             : 22,-1,-1,11,32,69,130,-1,-1,515,12,17,-1,-1,-1,7,-1,41,-1,257,-1,1027,80,-1,10,
    1305             : 133,-1,3,-1,-1,38,9,-1,19,-1,5,-1,-1,520,67,-1,-1,258,33,144,15,-1,-1,-1,-1,-1,
    1306             : 3,1032,-1,-1,21,96,-1,-1,9,6,-1,-1,131,-1,261,26,513,-1,35,-1,-1,-1,73,-1,7,-1,
    1307             : 17,2,-1,12,-1,160,1025,-1,5,-1,11,272,-1,70,1,-1,-1,-1,37,514,-1,-1,19,-1,137,
    1308             : -1,-1,-1,7,-1,65,42,-1,20,259,-1,-1,1026,13,-1,-1,-1,-1,134,49,-1,3,64,517,-1,
    1309             : -1,-1,11,-1,-1,18,-1,288,7,-1,-1,-1,129,-1,67,-1,25,-1,1029,-1,35,-1,-1,14,-1,
    1310             : -1,-1,528,5,-1,265,192,3,36,-1,-1,-1,-1,23,-1,-1,-1,-1,-1,11,-1,33,-1,69,1040,
    1311             : 131,8,-1,262,-1,-1,515,-1,13,34,17,-1,-1,-1,-1,74,-1,-1,7,128,-1,18,41,-1,-1,-1,
    1312             : 257,-1,-1,-1,1027,-1,81,6,-1,544,11,-1,133,-1,-1,-1,3,28,-1,-1,-1,-1,39,320,9,
    1313             : -1,-1,-1,19,-1,-1,138,5,-1,-1,1056,-1,6,521,-1,67,-1,-1,-1,-1,-1,259,-1,33,-1,
    1314             : 145,-1,15,-1,-1,-1,-1,68,-1,-1,-1,50,-1,-1,3,-1,1033,518,-1,384,-1,-1,21,10,97,
    1315             : -1,-1,-1,-1,-1,9,-1,7,-1,-1,-1,-1,44,131,-1,-1,66,261,-1,27,-1,513,1030,-1,-1,
    1316             : 35,-1,-1,-1,-1,-1,-1,132,73,-1,-1,-1,7,-1,-1,266,17,-1,3,-1,-1,-1,13,-1,-1,576,
    1317             : 161,22,1025,-1,-1,-1,5,-1,-1,-1,11,-1,273,34,-1,-1,71,-1,1,130,-1,-1,-1,-1,-1,
    1318             : -1,37,-1,515,-1,-1,14,-1,1088,19,256,-1,-1,137,-1,-1,-1,-1,-1,-1,24,7,-1,-1,-1,
    1319             : 65,-1,43,-1,-1,-1,21,640,259,-1,-1,-1,-1,-1,1027,-1,13,82,-1,-1,-1,-1,-1,10,-1,
    1320             : -1,135,-1,49,-1,-1,260,3,-1,65,-1,517,-1,-1,-1,-1,38,-1,-1,11,1152,-1,-1,-1,-1,
    1321             : 19,76,-1,-1,289,-1,7,-1,-1,-1,-1,20,-1,-1,129,522,-1,-1,67,-1,-1,-1,25,-1,-1,-1,
    1322             : 1029,258,-1,-1,35,4,-1,146,-1,-1,15,-1,-1,-1,-1,-1,-1,40,529,-1,5,-1,-1,-1,265,
    1323             : -1,193,-1,3,-1,37,1034,-1,-1,-1,-1,-1,-1,-1,-1,23,-1,-1,98,-1,140,-1,768,-1,-1,
    1324             : -1,-1,11,-1,-1,6,33,-1,-1,-1,69,-1,1041,-1,131,-1,9,-1,-1,-1,263,-1,-1,26,-1,-1,
    1325             : 515,-1,-1,-1,13,-1,35,-1,17,-1,-1,-1,-1,-1,-1,-1,-1,1280,75,52,-1,-1,-1,-1,7,-1,
    1326             : 129,-1,-1,516,19,-1,41,2,-1,-1,-1,-1,-1,14,257,-1,-1,-1,-1,162,-1,-1,1027,-1,-1,
    1327             : -1,81,-1,7,-1,-1,-1,545,-1,11,-1,-1,274,133,-1,-1,-1,-1,70,-1,-1,3,-1,29,-1,-1,
    1328             : -1,-1,1028,-1,-1,-1,-1,39,-1,321,514,9,-1,-1,-1,-1,-1,-1,-1,19,-1,-1,-1,-1,-1,
    1329             : 139,-1,5,-1,-1,-1,-1,268,1057,-1,-1,-1,7,-1,521,-1,-1,-1,67,-1,-1,42,-1,-1,-1,
    1330             : -1,-1,22,-1,-1,259,-1,-1,-1,33,72,-1,-1,145,1026,-1,-1,15,512,-1,-1,-1,36,-1,24,
    1331             : -1,-1,69,-1,-1,-1,-1,134,-1,-1,51,-1,-1,-1,-1,-1,3,-1,-1,66,1033,-1,519,-1,-1,
    1332             : -1,385,12,-1,-1,-1,-1,21,-1,11,-1,97,-1,-1,-1,-1,-1,-1,18,-1,-1,-1,-1,9,290,-1,
    1333             : 1536,7,-1,-1,-1,-1,-1,-1,-1,-1,-1,45,-1,131,-1,-1,-1,-1,-1,67,-1,261,-1,-1,-1,
    1334             : 27,-1,-1,-1,513,-1,1031,136,-1,-1,-1,84,35,-1,-1,-1,-1,-1,-1,-1,-1,14,-1,-1,-1,
    1335             : -1,133,-1,73,-1,-1,-1,-1,530,-1,-1,7,1024,-1,-1,-1,-1,267,-1,17,194,-1,-1,3,-1,
    1336             : -1,38,-1,-1,-1,-1,13,-1,-1,-1,-1,-1,577,-1,161,-1,23,-1,1025,-1,-1,-1,-1,-1,-1,
    1337             : -1,5,56,-1,-1,-1,-1,-1,-1,11,-1,-1,-1,273,-1,35,524,-1,-1,-1,-1,71,-1,-1,1042,1,
    1338             : -1,131,-1,-1,10,-1,-1,-1,-1,-1,262,-1,-1,-1,-1,37,-1,-1,-1,515,148,-1,-1,-1,-1,
    1339             : 15,-1,-1,34,1089,-1,19,-1,257,-1,-1,-1,-1,-1,137,-1,-1,-1,-1,-1,-1,74,-1,-1,-1,
    1340             : -1,-1,-1,25,-1,7,-1,-1,130,-1,1036,-1,-1,65,18,-1,-1,43,-1,-1,-1,-1,-1,-1,-1,21,
    1341             : -1,641,-1,259,100,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,1027,-1,-1,-1,13,-1,83,-1,-1,6,
    1342             : -1,264,-1,-1,-1,546,-1,-1,11,-1,-1,-1,-1,-1,135,-1,-1,-1,49,-1,-1,-1,-1,-1,261,
    1343             : -1,3,-1,-1,30,65,-1,-1,-1,517,-1,-1,-1,-1,-1,-1,-1,-1,-1,39,-1,-1,322,-1,-1,11,
    1344             : -1,1153,-1,-1,-1,-1,40,-1,-1,-1,-1,19,-1,77,-1,-1,-1,-1,-1,289,138,-1,-1,7
    1345             : };
    1346             : 
    1347             : /* Upper bound for H(v^2 d) / H(d) = \prod_{p | v} (p + 1) / (p - 1)
    1348             :  * We actually store ceil(128 * bound)
    1349             : P = primes(11); V = vector(31, i, vecsearch(P, i));
    1350             : hbound(v) =
    1351             : { my(q = factor(v)[,1]); if (q && vecmax(q) > 31, return (0));
    1352             :    ceil(prod(i = 1, #q, (q[i] + 1)/(q[i] - 1), 1.) * 128); }
    1353             : vector(1200, v, hbound(v))
    1354             : */
    1355             : static const long HURWITZ_RATIO[] = {
    1356             : 128,384,256,384,192,768,171,384,256,576,154,768,150,512,384,
    1357             : 384,144,768,143,576,342,461,140,768,192,448,256,512,138,1152,
    1358             : 137,384,308,432,256,768,0,427,299,576,0,1024,0,461,384,419,0,
    1359             : 768,171,576,288,448,0,768,231,512,285,412,0,1152,0,410,342,
    1360             : 384,224,922,0,432,280,768,0,768,0,0,384,427,205,896,0,576,
    1361             : 256,0,0,1024,216,0,275,461,0,1152,200,419,274,0,214,768,0,
    1362             : 512,308,576,0,864,0,448,512,0,0,768,0,692,0,512,0,854,210,
    1363             : 412,299,0,192,1152,154,0,0,410,192,1024,0,384,0,672,0,922,
    1364             : 190,0,384,432,0,838,0,768,0,0,180,768,206,0,342,0,0,1152,0,
    1365             : 427,288,615,205,896,0,0,0,576,187,768,0,0,461,0,0,1024,150,
    1366             : 648,285,0,0,823,256,461,0,0,0,1152,0,598,0,419,0,820,173,0,
    1367             : 342,640,0,768,0,0,448,512,0,922,0,576,0,0,183,864,0,0,280,
    1368             : 448,171,1536,0,0,0,0,0,768,183,0,0,692,168,0,0,512,384,0,
    1369             : 0,854,0,629,410,412,0,896,0,0,0,576,0,1152,0,461,256,0,256,
    1370             : 0,166,410,0,576,0,1024,168,0,432,384,0,0,0,672,275,0,0,922,
    1371             : 0,569,0,0,0,1152,0,432,399,0,231,838,0,0,274,768,0,0,0,0,
    1372             : 427,538,0,768,144,618,0,0,0,1024,0,0,308,0,163,1152,0,0,0,
    1373             : 427,0,864,0,615,0,615,0,896,0,0,512,0,0,0,165,576,0,559,
    1374             : 160,768,224,0,0,0,0,1383,0,0,0,0,0,1024,0,448,0,648,164,
    1375             : 854,171,0,419,0,0,823,0,768,299,461,0,0,0,0,384,0,0,1152,
    1376             : 143,0,308,598,0,0,0,419,0,0,0,820,0,519,384,0,160,1024,0,
    1377             : 640,0,0,0,768,308,0,0,0,0,1344,158,512,0,0,0,922,0,0,380,
    1378             : 576,0,0,160,0,384,549,0,864,0,0,0,0,0,838,0,448,0,512,0,
    1379             : 1536,0,0,0,0,216,0,0,0,359,0,0,768,0,547,412,0,156,0,0,
    1380             : 692,342,504,0,0,0,0,0,512,0,1152,0,0,0,0,299,854,0,0,288,
    1381             : 629,0,1229,0,412,410,0,0,896,0,0,0,0,0,0,214,576,0,0,0,
    1382             : 1152,0,0,373,461,0,768,0,0,0,768,0,0,155,498,461,410,0,0,
    1383             : 0,576,0,0,0,1024,0,503,299,0,0,1296,0,384,285,0,0,0,0,0,
    1384             : 0,672,0,823,0,0,512,0,154,922,140,0,0,569,0,0,0,0,0,0,
    1385             : 205,1152,0,0,0,432,0,1195,0,0,0,692,153,838,0,0,0,0,0,820,
    1386             : 0,768,346,0,0,0,0,0,342,0,0,1280,0,538,0,0,210,768,0,432,
    1387             : 0,618,0,0,0,0,448,0,0,1024,152,0,0,0,0,922,288,0,0,489,0,
    1388             : 1152,0,0,0,0,231,0,0,427,366,0,0,864,0,0,0,615,0,0,0,615,
    1389             : 280,0,0,896,192,0,342,0,0,1536,0,0,0,0,0,0,200,494,0,576,
    1390             : 0,0,0,559,0,480,0,768,0,672,365,0,0,0,0,0,0,0,0,1383,0,
    1391             : 0,336,0,285,0,150,0,0,0,0,1024,0,0,384,448,0,0,0,648,0,
    1392             : 492,0,854,0,512,0,0,0,1257,0,0,410,0,0,823,0,0,0,768,0,
    1393             : 896,0,461,0,0,0,0,0,0,0,0,149,1152,269,0,0,0,0,1152,0,
    1394             : 427,0,0,206,922,0,598,256,0,0,0,0,0,512,419,0,0,0,0,332,
    1395             : 0,0,820,0,0,0,519,0,1152,0,0,0,480,0,1024,0,0,336,640,0,
    1396             : 0,0,0,432,0,0,768,0,922,0,0,0,0,205,0,0,0,0,1344,0,472,
    1397             : 275,512,0,0,0,0,0,0,0,922,0,0,0,0,0,1138,0,576,0,0,0,0,
    1398             : 280,478,0,0,0,1152,0,549,0,0,0,864,0,0,399,0,0,0,0,0,461,
    1399             : 0,0,838,0,0,0,448,192,0,0,512,274,0,0,1536,138,0,0,0,224,
    1400             : 0,205,0,0,648,0,0,0,0,427,0,0,1076,0,0,0,0,0,768,0,0,
    1401             : 288,547,0,1235,0,0,0,466,256,0,0,0,0,692,0,1024,0,504,0,0,
    1402             : 0,0,0,0,308,0,0,0,0,512,326,0,147,1152,0,0,0,0,0,0,0,0,
    1403             : 0,896,0,854,0,0,0,0,0,864,0,629,0,0,0,1229,0,0,0,412,0,
    1404             : 1229,190,0,0,0,260,896,0,0,0,0,0,0,0,0,512,0,0,0,0,640,
    1405             : 0,576,0,0,0,0,330,0,0,1152,137,0,0,0,0,1118,0,461,320,0,
    1406             : 0,768,0,0,448,0,0,0,0,768,0,0,0,0,0,463,0,498,0,1383,0,
    1407             : 410,0,0,0,0,0,0,0,576,239,0,0,0,0,0,0,1024,0,0,0,503,0,
    1408             : 896,275,0,0,0,0,1296,0,0,328,384,0,854,0,0,342,0,0,0,0,0,
    1409             : 419,0,0,0,0,672,0,0,0,823,256,0,0,0,0,1536,0,0,299,461,0,
    1410             : 922,0,419,0,0,0,0,0,569,0,0,0,0,0,0,384,0,0,0,0,0,0,
    1411             : 615,0,1152,0,0,285,0,274,0,0,432,308,0,0,1195,0,0,0,0,0,
    1412             : 0,0,692,0,458,0,838,252,0,0,0,0,0,0,0,0,0,0,820,0,0,0,
    1413             : 768,0,1037,0,0,384,0,187,0,0,0,320,0,0,1024,0,0,0,0,0,
    1414             : 1280,0,0,0,538,0,0,0,0,0,629,0,768,0,0,615,432,0,0,0,618,
    1415             : 0,0,0,0,0,0,0,0,0,1344,0,0,315,0,0,1024,0,456,0,0,0,0,
    1416             : 200,0,0,0,0,922,0,864,0,0,0,0,0,489,380,0,0,1152
    1417             : };
    1418             : 
    1419             : /* Hurwitz class number of Df^2, D < 0; h = classno(D), Faf = factoru(f) */
    1420             : static double
    1421      187788 : hclassno_wrapper(long h, long D, GEN Faf)
    1422             : {
    1423             :   pari_sp av;
    1424      187788 :   if (lg(gel(Faf,1)) == 1) switch(D)
    1425             :   {
    1426           0 :     case -3: return 1. / 3;
    1427           0 :     case -4: return 1. / 2;
    1428       27365 :     default: return (double)h;
    1429             :   }
    1430      160423 :   av = avma;
    1431      160423 :   return (double)gc_long(av,  h * uhclassnoF_fact(Faf, D));
    1432             : }
    1433             : 
    1434             : /* return factor(u*v) */
    1435             : static GEN
    1436      182927 : factor_uv(GEN fau, ulong v, ulong vfactors)
    1437             : {
    1438             :   GEN P, E;
    1439             :   long i;
    1440      182927 :   if (!vfactors) return fau;
    1441      159106 :   P = gel(fau,1);
    1442      159106 :   E = gel(fau,2);
    1443      353888 :   for (i = 0; vfactors; i++, vfactors >>= 1)
    1444      353888 :     if (vfactors & 1UL)
    1445             :     {
    1446      185731 :       long p = SMALL_PRIMES[i];
    1447      185731 :       P = vecsmall_append(P, p);
    1448      185731 :       E = vecsmall_append(E, u_lvalrem(v, p, &v));
    1449      185731 :       if (v == 1) break;
    1450             :     }
    1451      159106 :   return famatsmall_reduce(mkmat2(P, E));
    1452             : }
    1453             : 
    1454             : /* This is Sutherland 2009, Algorithm 2.1 (p8); delta > 0 */
    1455             : static GEN
    1456        4861 : select_classpoly_prime_pool(double min_bits, double delta, GEN G)
    1457             : { /* Sutherland defines V_MAX to be 1200 without saying why */
    1458        4861 :   const long V_MAX = 1200;
    1459        4861 :   double bits = 0.0, hurwitz, z;
    1460        4861 :   ulong t_size_lim, d = (ulong)-pcp_get_D(G);
    1461        4861 :   long ires, inv = pcp_get_inv(G);
    1462        4861 :   GEN fau = pcp_get_fau(G);
    1463             :   GEN res, t_min; /* t_min[v] = lower bound for the t we look at for that v */
    1464             : #ifdef LONG_IS_64BIT
    1465        4158 :   long L = pcp_get_L(G)[!!pcp_get_L0(G)];
    1466             : #endif
    1467             : 
    1468        4861 :   hurwitz = hclassno_wrapper(pcp_get_h(G), pcp_get_D0(G), fau);
    1469             : 
    1470        4861 :   res = cgetg(128+1, t_VEC);
    1471        4861 :   ires = 1;
    1472             :   /* Initialise t_min to be all 2's.  This avoids trace 0 and trace 1 curves */
    1473        4861 :   t_min = const_vecsmall(V_MAX-1, 2);
    1474             : 
    1475             :   /* maximum possible trace = sqrt(2^BIL + D) */
    1476        4861 :   t_size_lim = 2.0 * sqrt((double)((1UL << (BITS_IN_LONG - 2)) - (d >> 2)));
    1477             : 
    1478        4861 :   for (z = d / (2.0 * hurwitz); ; z *= delta + 1.0)
    1479       18960 :   {
    1480       23821 :     double v_bound_aux = 4.0 * z * hurwitz; /* = 4 z H(d) */
    1481             :     ulong v;
    1482       23821 :     dbg_printf(1)("z = %.2f\n", z);
    1483      206777 :     for (v = 1; v < V_MAX; v++)
    1484             : #ifdef LONG_IS_64BIT
    1485      177540 :       if (L<=2 || v%L)
    1486             : #endif
    1487             :       {
    1488             :         ulong p, t, t_max, vfactors, v2d, vd;
    1489             :         double hurwitz_ratio_bound, max_p, H;
    1490             :         long ires0;
    1491             :         GEN faw;
    1492             : 
    1493      201887 :         if ((long)(vfactors = SMOOTH_INTS[v]) < 0) continue;
    1494      201887 :         hurwitz_ratio_bound = HURWITZ_RATIO[v] / 128.0;
    1495      201887 :         vd = v * d;
    1496      201887 :         if (vd >= v_bound_aux * hurwitz_ratio_bound) break;
    1497      182927 :         v2d = v * vd;
    1498      182927 :         faw = factor_uv(fau, v, vfactors);
    1499      182927 :         H = hclassno_wrapper(pcp_get_h(G), pcp_get_D0(G), faw);
    1500             :         /* t <= 2 sqrt(p) and p <= z H(v^2 d) and
    1501             :          *   H(v^2 d) < v H(d) \prod_{p | v} (p+1)/(p-1)
    1502             :          * This last term is v * hurwitz * hurwitz_ratio_bound. */
    1503      182927 :         max_p = z * v * hurwitz * hurwitz_ratio_bound;
    1504      182927 :         t_max = 2.0 * sqrt(mindd((1UL<<(BITS_IN_LONG-2)) - (v2d>>2), max_p));
    1505      182927 :         t = t_min[v]; if ((t & 1) != (v2d & 1)) t++;
    1506      182927 :         p = (t * t + v2d) >> 2;
    1507      182927 :         ires0 = ires;
    1508    10304893 :         for (; t <= t_max; p += t+1, t += 2) /* 4p = t^2 + v^2*d */
    1509    10121966 :           if (modinv_good_prime(inv,p) && uisprime(p))
    1510             :           {
    1511      284626 :             if (ires == lg(res)) res = vec_lengthen(res, lg(res) << 1);
    1512      284626 :             gel(res, ires++) = mkvec2(mkvecsmall5(p,t,v,(long)(p/H),vfactors),
    1513             :                 faw);
    1514      284626 :             bits += log2(p);
    1515             :           }
    1516      182927 :         t_min[v] = t;
    1517             : 
    1518      182927 :         if (ires - ires0) {
    1519       68937 :           dbg_printf(2)("  Found %lu primes for v = %lu.\n", ires - ires0, v);
    1520             :         }
    1521      182927 :         if (bits > min_bits) {
    1522        4861 :           dbg_printf(1)("Found %ld primes; total size %.2f bits.\n", ires-1,bits);
    1523        4861 :           setlg(res, ires); return res;
    1524             :         }
    1525             :       }
    1526       18960 :     if (uel(t_min,1) >= t_size_lim) {
    1527             :       /* exhausted all solutions that fit in ulong */
    1528           0 :       char *err = stack_sprintf("class polynomial of discriminant %ld", pcp_get_D(G));
    1529           0 :       pari_err(e_ARCH, err);
    1530             :     }
    1531             :   }
    1532             : }
    1533             : 
    1534             : static int
    1535     1216473 : primecmp(void *data, GEN v1, GEN v2)
    1536     1216473 : { (void)data; return cmpss(gel(v1,1)[4], gel(v2,1)[4]); }
    1537             : 
    1538             : static long
    1539        4861 : height_margin(long inv, long D)
    1540             : {
    1541             :   (void)D;
    1542             :   /* NB: avs just uses a height margin of 256 for everyone and everything. */
    1543        4861 :   if (inv == INV_F) return 64;  /* checked for discriminants up to -350000 */
    1544        4742 :   if (inv == INV_G2) return 5;
    1545        4441 :   if (inv != INV_J) return 256; /* TODO: This should be made more accurate */
    1546        3916 :   return 0;
    1547             : }
    1548             : 
    1549             : static GEN
    1550        4861 : select_classpoly_primes(ulong *vfactors, ulong *biggest_v, double delta,
    1551             :   GEN G, double height)
    1552             : {
    1553        4861 :   const long k = 2;
    1554        4861 :   pari_sp av = avma;
    1555        4861 :   long i, s, D = pcp_get_D(G), inv = pcp_get_inv(G);
    1556             :   ulong biggest_p;
    1557             :   double prime_bits, min_prime_bits, b;
    1558             :   GEN prime_pool;
    1559             : 
    1560        4861 :   s = modinv_height_factor(inv);
    1561        4861 :   b = height / s + height_margin(inv, D);
    1562        4861 :   dbg_printf(1)("adjusted height = %.2f\n", b);
    1563        4861 :   min_prime_bits = k * b;
    1564             : 
    1565        4861 :   prime_pool = select_classpoly_prime_pool(min_prime_bits, delta, G);
    1566             : 
    1567             :   /* FIXME: Apply torsion constraints */
    1568             :   /* FIXME: Rank elts of res according to cost/benefit ratio */
    1569        4861 :   gen_sort_inplace(prime_pool, NULL, primecmp, NULL);
    1570        4861 :   prime_bits = 0.0; biggest_p = *biggest_v = *vfactors = 0;
    1571      137319 :   for (i = 1; i < lg(prime_pool); i++)
    1572             :   {
    1573      137319 :     GEN q = gmael(prime_pool, i, 1);
    1574      137319 :     ulong p = q[1], v = q[3];
    1575      137319 :     *vfactors |= q[5];
    1576      137319 :     prime_bits += log2(p);
    1577      137319 :     if (p > biggest_p) biggest_p = p;
    1578      137319 :     if (v > *biggest_v) *biggest_v = v;
    1579      137319 :     if (prime_bits > b) break;
    1580             :   }
    1581        4861 :   dbg_printf(1)("Selected %ld primes; largest is %lu ~ 2^%.2f\n",
    1582             :              i, biggest_p, log2(biggest_p));
    1583        4861 :   return gerepilecopy(av, vecslice(prime_pool, 1, i));
    1584             : }
    1585             : 
    1586             : /* This is Sutherland 2009 Algorithm 1.2. */
    1587             : static long
    1588      138949 : oneroot_of_classpoly(
    1589             :   ulong *j_endo, int *endo_cert, ulong j, norm_eqn_t ne, GEN jdb)
    1590             : {
    1591      138949 :   pari_sp av = avma;
    1592             :   long nfactors, L_bound, i;
    1593      138949 :   ulong p = ne->p, pi = ne->pi;
    1594             :   GEN factors, vdepths;
    1595             : 
    1596      138949 :   if (j == 0 || j == 1728 % p) pari_err_BUG("oneroot_of_classpoly");
    1597             : 
    1598      138950 :   *endo_cert = 1;
    1599      138950 :   factors = gel(ne->faw, 1); nfactors = lg(factors) - 1;
    1600      138950 :   if (!nfactors) { *j_endo = j; return 1; }
    1601      135093 :   vdepths = gel(ne->faw, 2);
    1602             : 
    1603             :   /* FIXME: This should be bigger */
    1604      135093 :   L_bound = maxdd(log((double) -ne->D), (double)ne->v);
    1605             : 
    1606             :   /* Iterate over the primes L dividing w */
    1607      337447 :   for (i = 1; i <= nfactors; ++i) {
    1608      211140 :     pari_sp av2 = avma;
    1609             :     GEN phi;
    1610      211140 :     long jlvl, lvl_diff, depth = vdepths[i], L = factors[i];
    1611      211140 :     if (L > L_bound) { *endo_cert = 0; break; }
    1612             : 
    1613      202352 :     phi = polmodular_db_getp(jdb, L, p);
    1614             : 
    1615             :     /* TODO: See if I can reuse paths created in j_level_in_volcano()
    1616             :      * later in {ascend,descend}_volcano(), perhaps by combining the
    1617             :      * functions into one "adjust_level" function. */
    1618      202354 :     jlvl = j_level_in_volcano(phi, j, p, pi, L, depth);
    1619      202350 :     lvl_diff = z_lval(ne->u, L) - jlvl;
    1620      202359 :     if (lvl_diff < 0)
    1621             :       /* j's level is less than v(u) so we must ascend */
    1622      125229 :       j = ascend_volcano(phi, j, p, pi, jlvl, L, depth, -lvl_diff);
    1623       77130 :     else if (lvl_diff > 0)
    1624             :       /* otherwise j's level is greater than v(u) so we descend */
    1625        1265 :       j = descend_volcano(phi, j, p, pi, jlvl, L, depth, lvl_diff);
    1626      202360 :     set_avma(av2);
    1627             :   }
    1628             :   /* Prob(j has the wrong endomorphism ring) ~ \sum_{p|u_compl} 1/p
    1629             :    * (and u_compl > L_bound), so just return it and rely on detection code in
    1630             :    * enum_j_with_endo_ring().  Detection is that we hit a previously found
    1631             :    * j-invariant earlier than expected.  OR, we evaluate class polynomials of
    1632             :    * the suborders at j and if any are zero then j must be chosen again. */
    1633      135095 :   set_avma(av); *j_endo = j; return j != 0 && j != 1728 % p;
    1634             : }
    1635             : 
    1636             : INLINE long
    1637        7621 : vecsmall_isin_skip(GEN v, long x, long k)
    1638             : {
    1639        7621 :   long i, l = lg(v);
    1640      458409 :   for (i = k; i < l; ++i)
    1641      451258 :     if (v[i] == x) return i;
    1642        7151 :   return 0;
    1643             : }
    1644             : 
    1645             : void
    1646      174702 : norm_eqn_set(norm_eqn_t ne, long D, long t, long u, long v, GEN faw, ulong p)
    1647             : {
    1648      174702 :   ne->D = D;
    1649      174702 :   ne->u = u;
    1650      174702 :   ne->t = t;
    1651      174702 :   ne->v = v;
    1652      174702 :   ne->faw = faw;
    1653      174702 :   ne->p = p;
    1654      174702 :   ne->pi = get_Fl_red(ne->p);
    1655      174704 :   ne->s2 = Fl_2gener_pre(ne->p, ne->pi);
    1656             :   /* select twisting parameter */
    1657      349831 :   do ne->T = random_Fl(p); while (krouu(ne->T, p) != -1);
    1658      174726 : }
    1659             : 
    1660             : INLINE ulong
    1661       15978 : Flv_powsum_pre(GEN v, ulong n, ulong p, ulong pi)
    1662             : {
    1663       15978 :   long i, l = lg(v);
    1664       15978 :   ulong psum = 0;
    1665      375896 :   for (i = 1; i < l; ++i)
    1666      359918 :     psum = Fl_add(psum, Fl_powu_pre(uel(v,i), n, p, pi), p);
    1667       15978 :   return psum;
    1668             : }
    1669             : 
    1670             : INLINE int
    1671        4861 : modinv_has_sign_ambiguity(long inv)
    1672             : {
    1673        4861 :   switch (inv) {
    1674         553 :   case INV_F:
    1675             :   case INV_F3:
    1676             :   case INV_W2W3E2:
    1677             :   case INV_W2W7E2:
    1678             :   case INV_W2W3:
    1679             :   case INV_W2W5:
    1680             :   case INV_W2W7:
    1681             :   case INV_W3W3:
    1682             :   case INV_W2W13:
    1683         553 :   case INV_W3W7: return 1;
    1684             :   }
    1685        4308 :   return 0;
    1686             : }
    1687             : 
    1688             : INLINE int
    1689        3472 : modinv_units(int inv)
    1690        3472 : { return modinv_is_double_eta(inv) || modinv_is_Weber(inv); }
    1691             : 
    1692             : INLINE int
    1693       10843 : adjust_signs(GEN js, ulong p, ulong pi, long inv, GEN T, long e)
    1694             : {
    1695       10843 :   long negate = 0;
    1696       10843 :   long h = lg(js) - 1;
    1697       14098 :   if ((h & 1) && modinv_units(inv)) {
    1698        3255 :     ulong prod = Flv_prod_pre(js, p, pi);
    1699        3255 :     if (prod != p - 1) {
    1700        1743 :       if (prod != 1) pari_err_BUG("adjust_signs: constant term is not +/-1");
    1701        1743 :       negate = 1;
    1702             :     }
    1703             :   } else {
    1704             :     ulong tp, t;
    1705        7588 :     tp = umodiu(T, p);
    1706        7588 :     t = Flv_powsum_pre(js, e, p, pi);
    1707        7588 :     if (t == 0) return 0;
    1708        7574 :     if (t != tp) {
    1709        3767 :       if (Fl_neg(t, p) != tp) pari_err_BUG("adjust_signs: incorrect trace");
    1710        3767 :       negate = 1;
    1711             :     }
    1712             :   }
    1713       10829 :   if (negate) Flv_neg_inplace(js, p);
    1714       10829 :   return 1;
    1715             : }
    1716             : 
    1717             : static ulong
    1718      138375 : find_jinv(
    1719             :   long *trace_tries, long *endo_tries, int *cert,
    1720             :   norm_eqn_t ne, long inv, long rho_inv, GEN jdb)
    1721             : {
    1722      138375 :   long found, ok = 1;
    1723             :   ulong j, r;
    1724             :   do {
    1725             :     do {
    1726             :       long tries;
    1727      138668 :       ulong j_t = 0;
    1728             :       /* TODO: Set batch size according to expected number of tries and
    1729             :        * experimental cost/benefit analysis. */
    1730      138668 :       tries = find_j_inv_with_given_trace(&j_t, ne, rho_inv, 0);
    1731      138949 :       if (j_t == 0)
    1732           0 :         pari_err_BUG("polclass0: Couldn't find j-invariant with given trace.");
    1733      138949 :       dbg_printf(2)("  j-invariant %ld has trace +/-%ld (%ld tries, 1/rho = %ld)\n",
    1734             :           j_t, ne->t, tries, rho_inv);
    1735      138949 :       *trace_tries += tries;
    1736             : 
    1737      138949 :       found = oneroot_of_classpoly(&j, cert, j_t, ne, jdb);
    1738      138952 :       ++*endo_tries;
    1739      138952 :     } while (!found);
    1740             : 
    1741      138904 :     if (modinv_is_double_eta(inv))
    1742       12294 :       ok = modfn_unambiguous_root(&r, inv, j, ne, jdb);
    1743             :     else
    1744      126610 :       r = modfn_root(j, ne, inv);
    1745      138903 :   } while (!ok);
    1746      138658 :   return r;
    1747             : }
    1748             : 
    1749             : static GEN
    1750      137069 : polclass_roots_modp(long *n_trace_curves,
    1751             :   norm_eqn_t ne, long rho_inv, GEN G, GEN db)
    1752             : {
    1753             :   pari_sp av;
    1754      137069 :   ulong j = 0;
    1755      137069 :   long inv = pcp_get_inv(G), endo_tries = 0;
    1756             :   int endo_cert;
    1757      137080 :   GEN res, jdb, fdb, vshape = factoru(ne->v);
    1758             : 
    1759      136934 :   jdb = polmodular_db_for_inv(db, INV_J);
    1760      136995 :   fdb = polmodular_db_for_inv(db, inv);
    1761      137018 :   dbg_printf(2)("p = %ld, t = %ld, v = %ld\n", ne->p, ne->t, ne->v);
    1762      137019 :   av = avma;
    1763             :   do {
    1764      138363 :     j = find_jinv(n_trace_curves,&endo_tries,&endo_cert, ne, inv, rho_inv, jdb);
    1765      138658 :     res = enum_roots(j, ne, fdb, G, vshape);
    1766      138653 :     if (!res && endo_cert) pari_err_BUG("polclass_roots_modp");
    1767      138653 :     if (res && !endo_cert && vecsmall_isin_skip(res, res[1], 2))
    1768             :     {
    1769         470 :       set_avma(av);
    1770         470 :       res = NULL;
    1771             :     }
    1772      138653 :   } while (!res);
    1773             : 
    1774      137309 :   dbg_printf(2)("  j-invariant %ld has correct endomorphism ring "
    1775             :              "(%ld tries)\n", j, endo_tries);
    1776      137309 :   dbg_printf(4)("  all such j-invariants: %Ps\n", res);
    1777      137309 :   return res;
    1778             : }
    1779             : 
    1780             : INLINE int
    1781        2225 : modinv_inverted_involution(long inv)
    1782        2225 : { return modinv_is_double_eta(inv); }
    1783             : 
    1784             : INLINE int
    1785        2225 : modinv_negated_involution(long inv)
    1786             : { /* determined by trial and error */
    1787        2225 :   return inv == INV_F || inv == INV_W3W5 || inv == INV_W3W7
    1788        4450 :     || inv == INV_W3W3 || inv == INV_W5W7;
    1789             : }
    1790             : 
    1791             : /* Return true iff Phi_L(j0, j1) = 0. */
    1792             : INLINE long
    1793        3966 : verify_edge(ulong j0, ulong j1, ulong p, ulong pi, long L, GEN fdb)
    1794             : {
    1795        3966 :   pari_sp av = avma;
    1796        3966 :   GEN phi = polmodular_db_getp(fdb, L, p);
    1797        3966 :   GEN f = Flm_Fl_polmodular_evalx(phi, L, j1, p, pi);
    1798        3966 :   return gc_long(av, Flx_eval_pre(f, j0, p, pi) == 0);
    1799             : }
    1800             : 
    1801             : INLINE long
    1802         672 : verify_2path(
    1803             :   ulong j1, ulong j2, ulong p, ulong pi, long L1, long L2, GEN fdb)
    1804             : {
    1805         672 :   pari_sp av = avma;
    1806         672 :   GEN phi1 = polmodular_db_getp(fdb, L1, p);
    1807         672 :   GEN phi2 = polmodular_db_getp(fdb, L2, p);
    1808         672 :   GEN f = Flm_Fl_polmodular_evalx(phi1, L1, j1, p, pi);
    1809         672 :   GEN g = Flm_Fl_polmodular_evalx(phi2, L2, j2, p, pi);
    1810         672 :   GEN d = Flx_gcd(f, g, p);
    1811         672 :   long n = degpol(d);
    1812         672 :   if (n >= 2) n = Flx_nbroots(d, p);
    1813         672 :   return gc_long(av, n);
    1814             : }
    1815             : 
    1816             : static long
    1817        6159 : oriented_n_action(
    1818             :   const long *ni, GEN G, GEN v, ulong p, ulong pi, GEN fdb)
    1819             : {
    1820        6159 :   pari_sp av = avma;
    1821        6159 :   long i, j, k = pcp_get_k(G);
    1822        6159 :   long nr = k * (k - 1) / 2;
    1823        6159 :   const long *n = pcp_get_n(G), *m = pcp_get_m(G), *o = pcp_get_o(G), *r = pcp_get_r(G),
    1824        6159 :     *ps = pcp_get_orient_p(G), *qs = pcp_get_orient_q(G), *reps = pcp_get_orient_reps(G);
    1825        6159 :   long *signs = new_chunk(k);
    1826        6159 :   long *e = new_chunk(k);
    1827        6159 :   long *rels = new_chunk(nr);
    1828             : 
    1829        6159 :   evec_copy(rels, r, nr);
    1830             : 
    1831       16635 :   for (i = 0; i < k; ++i) {
    1832             :     /* If generator doesn't require orientation, continue; power rels already
    1833             :      * copied to *rels in initialisation */
    1834       10476 :     if (ps[i] <= 0) { signs[i] = 1; continue; }
    1835             :     /* Get rep of orientation element and express it in terms of the
    1836             :      * (partially) oriented presentation */
    1837        6508 :     for (j = 0; j < i; ++j) {
    1838        4189 :       long t = reps[i * k + j];
    1839        4189 :       e[j] = (signs[j] < 0 ? o[j] - t : t);
    1840             :     }
    1841        2319 :     e[j] = reps[i * k + j];
    1842        3474 :     for (++j; j < k; ++j) e[j] = 0;
    1843        2319 :     evec_reduce(e, n, rels, k);
    1844        2319 :     j = evec_to_index(e, m, k);
    1845             : 
    1846             :     /* FIXME: These calls to verify_edge recalculate powers of v[0]
    1847             :      * and v[j] over and over again, they also reduce Phi_{ps[i]} modulo p over
    1848             :      * and over again.  Need to cache these things! */
    1849        2319 :     if (qs[i] > 1)
    1850         336 :       signs[i] =
    1851         336 :         (verify_2path(uel(v,1), uel(v,j+1), p, pi, ps[i], qs[i], fdb) ? 1 : -1);
    1852             :     else
    1853             :       /* Verify ps[i]-edge to orient ith generator */
    1854        1983 :       signs[i] =
    1855        1983 :         (verify_edge(uel(v,1), uel(v,j+1), p, pi, ps[i], fdb) ? 1 : -1);
    1856             :     /* Update power relation */
    1857        6508 :     for (j = 0; j < i; ++j) {
    1858        4189 :       long t = evec_ri(r, i)[j];
    1859        4189 :       e[j] = (signs[i] * signs[j] < 0 ? o[j] - t : t);
    1860             :     }
    1861        5793 :     while (j < k) e[j++] = 0;
    1862        2319 :     evec_reduce(e, n, rels, k);
    1863        6508 :     for (j = 0; j < i; ++j) evec_ri_mutate(rels, i)[j] = e[j];
    1864             :     /* TODO: This is a sanity check, can be removed if everything is working */
    1865        8827 :     for (j = 0; j <= i; ++j) {
    1866        6508 :       long t = reps[i * k + j];
    1867        6508 :       e[j] = (signs[j] < 0 ? o[j] - t : t);
    1868             :     }
    1869        3474 :     while (j < k) e[j++] = 0;
    1870        2319 :     evec_reduce(e, n, rels, k);
    1871        2319 :     j = evec_to_index(e, m, k);
    1872        2319 :     if (qs[i] > 1) {
    1873         336 :       if (!verify_2path(uel(v,1), uel(v, j+1), p, pi, ps[i], qs[i], fdb))
    1874           0 :         pari_err_BUG("oriented_n_action");
    1875             :     } else {
    1876        1983 :       if (!verify_edge(uel(v,1), uel(v, j+1), p, pi, ps[i], fdb))
    1877           0 :         pari_err_BUG("oriented_n_action");
    1878             :     }
    1879             :   }
    1880             : 
    1881             :   /* Orient representation of [N] relative to the torsor <signs, rels> */
    1882       16635 :   for (i = 0; i < k; ++i) e[i] = (signs[i] < 0 ? o[i] - ni[i] : ni[i]);
    1883        6159 :   evec_reduce(e, n, rels, k);
    1884        6159 :   return gc_long(av, evec_to_index(e,m,k));
    1885             : }
    1886             : 
    1887             : /* F = double_eta_raw(inv) */
    1888             : INLINE void
    1889        6159 : adjust_orientation(GEN F, long inv, GEN v, long e, ulong p, ulong pi)
    1890             : {
    1891        6159 :   ulong j0 = uel(v, 1), je = uel(v, e);
    1892             : 
    1893        6159 :   if (!modinv_j_from_2double_eta(F, inv, j0, je, p, pi)) {
    1894        2225 :     if (modinv_inverted_involution(inv)) Flv_inv_pre_inplace(v, p, pi);
    1895        2225 :     if (modinv_negated_involution(inv)) Flv_neg_inplace(v, p);
    1896             :   }
    1897        6159 : }
    1898             : 
    1899             : static void
    1900         553 : polclass_psum(
    1901             :   GEN *psum, long *d, GEN roots, GEN primes, GEN pilist, ulong h, long inv)
    1902             : {
    1903             :   /* Number of consecutive CRT stabilisations before we assume we have
    1904             :    * the correct answer. */
    1905             :   enum { MIN_STAB_CNT = 3 };
    1906         553 :   pari_sp av = avma, btop;
    1907             :   GEN ps, psum_sqr, P;
    1908         553 :   long i, e, stabcnt, nprimes = lg(primes) - 1;
    1909             : 
    1910         553 :   if ((h & 1) && modinv_units(inv)) { *psum = gen_1; *d = 0; return; }
    1911         336 :   e = -1;
    1912         336 :   ps = cgetg(nprimes+1, t_VECSMALL);
    1913             :   do {
    1914         378 :     e += 2;
    1915        8768 :     for (i = 1; i <= nprimes; ++i)
    1916             :     {
    1917        8390 :       GEN roots_modp = gel(roots, i);
    1918        8390 :       ulong p = uel(primes, i), pi = uel(pilist, i);
    1919        8390 :       uel(ps, i) = Flv_powsum_pre(roots_modp, e, p, pi);
    1920             :     }
    1921         378 :     btop = avma;
    1922         378 :     psum_sqr = Z_init_CRT(0, 1);
    1923         378 :     P = gen_1;
    1924        2149 :     for (i = 1, stabcnt = 0; stabcnt < MIN_STAB_CNT && i <= nprimes; ++i)
    1925             :     {
    1926        1771 :       ulong p = uel(primes, i), pi = uel(pilist, i);
    1927        1771 :       ulong ps2 = Fl_sqr_pre(uel(ps, i), p, pi);
    1928        1771 :       ulong stab = Z_incremental_CRT(&psum_sqr, ps2, &P, p);
    1929             :       /* stabcnt = stab * (stabcnt + 1) */
    1930        1771 :       if (stab) ++stabcnt; else stabcnt = 0;
    1931        1771 :       if (gc_needed(av, 2)) gerepileall(btop, 2, &psum_sqr, &P);
    1932             :     }
    1933         378 :     if (stabcnt == 0 && nprimes >= MIN_STAB_CNT)
    1934           0 :       pari_err_BUG("polclass_psum");
    1935         378 :   } while (!signe(psum_sqr));
    1936             : 
    1937         336 :   if ( ! Z_issquareall(psum_sqr, psum)) pari_err_BUG("polclass_psum");
    1938             : 
    1939         336 :   dbg_printf(1)("Classpoly power sum (e = %ld) is %Ps; found with %.2f%% of the primes\n",
    1940           0 :       e, *psum, 100 * (i - 1) / (double) nprimes);
    1941         336 :   *psum = gerepileupto(av, *psum);
    1942         336 :   *d = e;
    1943             : }
    1944             : 
    1945             : static GEN
    1946        1491 : polclass_small_disc(long D, long inv, long vx)
    1947             : {
    1948        1491 :   if (D == -3) return pol_x(vx);
    1949         546 :   if (D == -4) {
    1950         546 :     switch (inv) {
    1951         364 :     case INV_J: return deg1pol(gen_1, stoi(-1728), vx);
    1952         182 :     case INV_G2:return deg1pol(gen_1, stoi(-12), vx);
    1953           0 :     default: /* no other invariants for which we can calculate H_{-4}(X) */
    1954           0 :       pari_err_BUG("polclass_small_disc");
    1955             :     }
    1956             :   }
    1957           0 :   return NULL;
    1958             : }
    1959             : 
    1960             : static ulong
    1961        4861 : quadnegclassnou(long D, long *pD0, GEN *pP, GEN *pE)
    1962             : {
    1963        4861 :   ulong d = (ulong)-D, d0 = coredisc2u_fact(factoru(d), -1, pP, pE);
    1964        4861 :   ulong h = uquadclassnoF_fact(d0, -1, *pP, *pE) * quadclassnos(-(long)d0);
    1965        4861 :   *pD0 = -(long)d0; return h;
    1966             : }
    1967             : 
    1968             : GEN
    1969      137167 : polclass_worker(GEN q, GEN G, GEN db)
    1970             : {
    1971      137167 :   GEN T = cgetg(3, t_VEC), z, P = gel(q,1), faw = gel(q,2);
    1972      137051 :   long n_curves_tested = 0, t = P[2], v = P[3], rho_inv = P[4];
    1973      137051 :   ulong p = (ulong)P[1];
    1974      137051 :   pari_sp av = avma;
    1975      137051 :   norm_eqn_t ne; norm_eqn_set(ne, pcp_get_D(G), t, pcp_get_u(G), v, faw, p);
    1976      137081 :   z = polclass_roots_modp(&n_curves_tested, ne, rho_inv, G, db);
    1977      137310 :   gel(T,1) = gerepileuptoleaf(av, z);
    1978      137315 :   gel(T,2) = mkvecsmall3(ne->p, ne->pi, n_curves_tested); return T;
    1979             : }
    1980             : 
    1981             : GEN
    1982        6352 : polclass0(long D, long inv, long vx, GEN *db)
    1983             : {
    1984        6352 :   pari_sp av = avma;
    1985             :   GEN primes, H, plist, pilist, Pu, Eu;
    1986        6352 :   long n_curves_tested = 0, filter = 1, pending = 0, cnt = 0;
    1987             :   long D0, nprimes, s, i, j, del, ni, orient, h, p1, p2, k;
    1988             :   ulong u, vfactors, biggest_v;
    1989             :   GEN G, worker, vec;
    1990             :   double height;
    1991             :   static const double delta = 0.5;
    1992             :   struct pari_mt pt;
    1993             : 
    1994        6352 :   if (D >= -4) return polclass_small_disc(D, inv, vx);
    1995             : 
    1996        4861 :   h = quadnegclassnou(D, &D0, &Pu, &Eu);
    1997        4861 :   u = D == D0? 1: (ulong)sqrt(D / D0); /* u = \prod Pu[i]^Eu[i] */
    1998        4861 :   dbg_printf(1)("D = %ld, conductor = %ld, inv = %ld\n", D, u, inv);
    1999             : 
    2000        4861 :   ni = modinv_degree(&p1, &p2, inv);
    2001        4861 :   orient = modinv_is_double_eta(inv) && kross(D, p1) && kross(D, p2);
    2002             : 
    2003        4861 :   G = classgp_make_pcp(&height, &ni, h, D, D0, u, Pu, Eu, inv, filter, orient);
    2004        4861 :   primes = select_classpoly_primes(&vfactors, &biggest_v, delta, G, height);
    2005             : 
    2006             :   /* Prepopulate *db with all the modpolys we might need */
    2007        4861 :   if (u > 1)
    2008             :   { /* L_bound in oneroot_of_classpoly() */
    2009         322 :     long l = lg(Pu), maxL = maxdd(log((double) -D), (double)biggest_v);
    2010         550 :     for (i = 1; i < l; i++)
    2011             :     {
    2012         364 :       long L = Pu[i];
    2013         364 :       if (L > maxL) break;
    2014         228 :       polmodular_db_add_level(db, L, INV_J);
    2015             :     }
    2016             :   }
    2017       19255 :   for (i = 0; vfactors; ++i) {
    2018       14394 :     if (vfactors & 1UL)
    2019       13664 :       polmodular_db_add_level(db, SMALL_PRIMES[i], INV_J);
    2020       14394 :     vfactors >>= 1;
    2021             :   }
    2022        4861 :   if (p1 > 1) polmodular_db_add_level(db, p1, INV_J);
    2023        4861 :   if (p2 > 1) polmodular_db_add_level(db, p2, INV_J);
    2024        4861 :   s = !!pcp_get_L0(G); k = pcp_get_k(G);
    2025        4861 :   polmodular_db_add_levels(db, pcp_get_L(G) + s, k - s, inv);
    2026        4861 :   if (orient)
    2027             :   {
    2028         273 :     GEN orient_p = pcp_get_orient_p(G);
    2029         273 :     GEN orient_q = pcp_get_orient_q(G);
    2030         714 :     for (i = 0; i < k; ++i)
    2031             :     {
    2032         441 :       if (orient_p[i] > 1) polmodular_db_add_level(db, orient_p[i], inv);
    2033         441 :       if (orient_q[i] > 1) polmodular_db_add_level(db, orient_q[i], inv);
    2034             :     }
    2035             :   }
    2036        4861 :   nprimes = lg(primes) - 1;
    2037        4861 :   worker = snm_closure(is_entry("_polclass_worker"), mkvec2(G, *db));
    2038        4861 :   H = cgetg(nprimes + 1, t_VEC);
    2039        4861 :   plist = cgetg(nprimes + 1, t_VECSMALL);
    2040        4861 :   pilist = cgetg(nprimes + 1, t_VECSMALL);
    2041        4861 :   vec = cgetg(2, t_VEC);
    2042        4861 :   dbg_printf(0)("Calculating class polynomial of disc %ld and inv %ld:", D, inv);
    2043        4861 :   mt_queue_start_lim(&pt, worker, nprimes);
    2044      150757 :   for (i = 1; i <= nprimes || pending; i++)
    2045             :   {
    2046             :     long workid;
    2047             :     GEN done;
    2048      145896 :     if (i <= nprimes) gel(vec,1) = gel(primes,i);
    2049      145896 :     mt_queue_submit(&pt, i, i <= nprimes? vec: NULL);
    2050      145896 :     done = mt_queue_get(&pt, &workid, &pending);
    2051      145896 :     if (done)
    2052             :     {
    2053      137319 :       gel(H, workid) = gel(done,1);
    2054      137319 :       uel(plist, workid) = mael(done,2,1);
    2055      137319 :       uel(pilist, workid) = mael(done,2,2);
    2056      137319 :       n_curves_tested += mael(done,2,3);
    2057      137319 :       cnt++;
    2058      137319 :       if (DEBUGLEVEL>2 && (cnt & 3L)==0)
    2059           0 :         err_printf(" %ld%%", cnt*100/nprimes);
    2060             :     }
    2061             :   }
    2062        4861 :   mt_queue_end(&pt);
    2063        4861 :   dbg_printf(0)(" done\n");
    2064             : 
    2065        4861 :   if (orient) {
    2066         273 :     GEN nvec = new_chunk(k);
    2067         273 :     GEN fdb = polmodular_db_for_inv(*db, inv);
    2068         273 :     GEN F = double_eta_raw(inv);
    2069         273 :     index_to_evec((long *)nvec, ni, pcp_get_m(G), k);
    2070        6432 :     for (i = 1; i <= nprimes; ++i) {
    2071        6159 :       GEN v = gel(H, i);
    2072        6159 :       ulong p = uel(plist, i), pi = uel(pilist, i);
    2073        6159 :       long oni = oriented_n_action(nvec, G, v, p, pi, fdb);
    2074        6159 :       adjust_orientation(F, inv, v, oni + 1, p, pi);
    2075             :     }
    2076             :   }
    2077             : 
    2078        4861 :   if (modinv_has_sign_ambiguity(inv)) {
    2079             :     GEN psum;
    2080             :     long e;
    2081         553 :     polclass_psum(&psum, &e, H, plist, pilist, h, inv);
    2082       11396 :     for (i = 1; i <= nprimes; ++i) {
    2083       10843 :       GEN v = gel(H, i);
    2084       10843 :       ulong p = uel(plist, i), pi = uel(pilist, i);
    2085       10843 :       if (!adjust_signs(v, p, pi, inv, psum, e))
    2086          14 :         uel(plist, i) = 0;
    2087             :     }
    2088             :   }
    2089             : 
    2090      142180 :   for (i = 1, j = 1, del = 0; i <= nprimes; ++i)
    2091             :   {
    2092      137319 :     pari_sp av2 = avma;
    2093      137319 :     ulong p = uel(plist, i);
    2094             :     long l;
    2095             :     GEN v;
    2096      137319 :     if (!p) { del++; continue; }
    2097      137305 :     v = Flv_roots_to_pol(gel(H,i), p, vx); l = lg(v);
    2098      137305 :     *++v = evaltyp(t_VECSMALL) | _evallg(l-1); /* Flx_to_Flv inplace */
    2099      137305 :     uel(plist, j) = p;
    2100      137305 :     gel(H, j++) = gerepileuptoleaf(av2, v);
    2101             :   }
    2102        4861 :   setlg(H,nprimes+1-del);
    2103        4861 :   setlg(plist,nprimes+1-del);
    2104             : 
    2105        4861 :   dbg_printf(1)("Total number of curves tested: %ld\n", n_curves_tested);
    2106        4861 :   H = ncV_chinese_center(H, plist, NULL);
    2107        4861 :   dbg_printf(1)("Result height: %.2f\n", dbllog2(gsupnorm(H, DEFAULTPREC)));
    2108        4861 :   return gerepilecopy(av, RgV_to_RgX(H, vx));
    2109             : }
    2110             : 
    2111             : void
    2112        2786 : check_modinv(long inv)
    2113             : {
    2114        2786 :   switch (inv) {
    2115        2772 :   case INV_J:
    2116             :   case INV_F:
    2117             :   case INV_F2:
    2118             :   case INV_F3:
    2119             :   case INV_F4:
    2120             :   case INV_G2:
    2121             :   case INV_W2W3:
    2122             :   case INV_F8:
    2123             :   case INV_W3W3:
    2124             :   case INV_W2W5:
    2125             :   case INV_W2W7:
    2126             :   case INV_W3W5:
    2127             :   case INV_W3W7:
    2128             :   case INV_W2W3E2:
    2129             :   case INV_W2W5E2:
    2130             :   case INV_W2W13:
    2131             :   case INV_W2W7E2:
    2132             :   case INV_W3W3E2:
    2133             :   case INV_W5W7:
    2134             :   case INV_W3W13:
    2135        2772 :     break;
    2136          14 :   default:
    2137          14 :     pari_err_DOMAIN("polmodular", "inv", "invalid invariant", stoi(inv), gen_0);
    2138             :   }
    2139        2772 : }
    2140             : 
    2141             : GEN
    2142        2170 : polclass(GEN DD, long inv, long vx)
    2143             : {
    2144             :   GEN db, H;
    2145             :   long D;
    2146             : 
    2147        2170 :   if (vx < 0) vx = 0;
    2148        2170 :   check_quaddisc_imag(DD, NULL, "polclass");
    2149        2163 :   check_modinv(inv);
    2150             : 
    2151        2156 :   D = itos(DD);
    2152        2156 :   if (!modinv_good_disc(inv, D))
    2153           0 :     pari_err_DOMAIN("polclass", "D", "incompatible with given invariant", stoi(inv), DD);
    2154             : 
    2155        2156 :   db = polmodular_db_init(inv);
    2156        2156 :   H = polclass0(D, inv, vx, &db);
    2157        2156 :   gunclone_deep(db); return H;
    2158             : }

Generated by: LCOV version 1.14