Code coverage tests

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

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

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

LCOV - code coverage report
Current view: top level - basemath - polclass.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.10.0 lcov report (development 22296-eb24562) Lines: 842 867 97.1 %
Date: 2018-04-19 06:16:12 Functions: 59 60 98.3 %
Legend: Lines: hit not hit

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

Generated by: LCOV version 1.11