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 - polmodular.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.11.0 lcov report (development 22860-5579deb0b) Lines: 2276 2348 96.9 %
Date: 2018-07-19 05:36:41 Functions: 142 142 100.0 %
Legend: Lines: hit not hit

          Line data    Source code
       1             : /* Copyright (C) 2014  The PARI group.
       2             : 
       3             : This file is part of the PARI/GP package.
       4             : 
       5             : PARI/GP is free software; you can redistribute it and/or modify it under the
       6             : terms of the GNU General Public License as published by the Free Software
       7             : Foundation. 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             :  * START Code from AVSs "class_inv.h"
      21             :  */
      22             : 
      23             : /* actually just returns the square-free part of the level, which is
      24             :  * all we care about */
      25             : long
      26       31717 : modinv_level(long inv)
      27             : {
      28       31717 :   switch (inv) {
      29       23391 :     case INV_J:     return 1;
      30             :     case INV_G2:
      31        1239 :     case INV_W3W3E2:return 3;
      32             :     case INV_F:
      33             :     case INV_F2:
      34             :     case INV_F4:
      35        1091 :     case INV_F8:    return 6;
      36          70 :     case INV_F3:    return 2;
      37         448 :     case INV_W3W3:  return 6;
      38             :     case INV_W2W7E2:
      39        1498 :     case INV_W2W7:  return 14;
      40         269 :     case INV_W3W5:  return 15;
      41             :     case INV_W2W3E2:
      42         301 :     case INV_W2W3:  return 6;
      43             :     case INV_W2W5E2:
      44         644 :     case INV_W2W5:  return 30;
      45         441 :     case INV_W2W13: return 26;
      46        1725 :     case INV_W3W7:  return 42;
      47         544 :     case INV_W5W7:  return 35;
      48          56 :     case INV_W3W13: return 39;
      49             :   }
      50             :   pari_err_BUG("modinv_level"); return 0;/*LCOV_EXCL_LINE*/
      51             : }
      52             : 
      53             : /* Where applicable, returns N=p1*p2 (possibly p2=1) s.t. two j's
      54             :  * related to the same f are N-isogenous, and 0 otherwise.  This is
      55             :  * often (but not necessarily) equal to the level. */
      56             : long
      57     7183622 : modinv_degree(long *P1, long *P2, long inv)
      58             : {
      59             :   long *p1, *p2, ignored;
      60     7183622 :   p1 = P1? P1: &ignored;
      61     7183622 :   p2 = P2? P2: &ignored;
      62     7183622 :   switch (inv) {
      63      295707 :     case INV_W3W5:  return (*p1 = 3) * (*p2 = 5);
      64             :     case INV_W2W3E2:
      65      424830 :     case INV_W2W3:  return (*p1 = 2) * (*p2 = 3);
      66             :     case INV_W2W5E2:
      67     1536269 :     case INV_W2W5:  return (*p1 = 2) * (*p2 = 5);
      68             :     case INV_W2W7E2:
      69      931882 :     case INV_W2W7:  return (*p1 = 2) * (*p2 = 7);
      70     1461348 :     case INV_W2W13: return (*p1 = 2) * (*p2 = 13);
      71      508139 :     case INV_W3W7:  return (*p1 = 3) * (*p2 = 7);
      72             :     case INV_W3W3E2:
      73      770735 :     case INV_W3W3:  return (*p1 = 3) * (*p2 = 3);
      74      347823 :     case INV_W5W7:  return (*p1 = 5) * (*p2 = 7);
      75      194054 :     case INV_W3W13: return (*p1 = 3) * (*p2 = 13);
      76             :   }
      77      712835 :   *p1 = *p2 = 1; return 0;
      78             : }
      79             : 
      80             : /* Certain invariants require that D not have 2 in it's conductor, but
      81             :  * this doesn't apply to every invariant with even level so we handle
      82             :  * it separately */
      83             : INLINE int
      84      457168 : modinv_odd_conductor(long inv)
      85             : {
      86      457168 :   switch (inv) {
      87             :     case INV_F:
      88             :     case INV_W3W3:
      89       65433 :     case INV_W3W7: return 1;
      90             :   }
      91      391735 :   return 0;
      92             : }
      93             : 
      94             : long
      95    21028719 : modinv_height_factor(long inv)
      96             : {
      97    21028719 :   switch (inv) {
      98     1746628 :     case INV_J:     return 1;
      99     1617497 :     case INV_G2:    return 3;
     100     2905184 :     case INV_F:     return 72;
     101          28 :     case INV_F2:    return 36;
     102      473893 :     case INV_F3:    return 24;
     103          49 :     case INV_F4:    return 18;
     104          49 :     case INV_F8:    return 9;
     105          63 :     case INV_W2W3:  return 72;
     106     2217376 :     case INV_W3W3:  return 36;
     107     3265507 :     case INV_W2W5:  return 54;
     108     1184226 :     case INV_W2W7:  return 48;
     109        1190 :     case INV_W3W5:  return 36;
     110     3587794 :     case INV_W2W13: return 42;
     111     1048005 :     case INV_W3W7:  return 32;
     112     1059933 :     case INV_W2W3E2:return 36;
     113      132874 :     case INV_W2W5E2:return 27;
     114      995421 :     case INV_W2W7E2:return 24;
     115          49 :     case INV_W3W3E2:return 18;
     116      792939 :     case INV_W5W7:  return 24;
     117          14 :     case INV_W3W13: return 28;
     118             :     default: pari_err_BUG("modinv_height_factor"); return 0;/*LCOV_EXCL_LINE*/
     119             :   }
     120             : }
     121             : 
     122             : long
     123     1895915 : disc_best_modinv(long D)
     124             : {
     125             :   long ret;
     126     1895915 :   ret = INV_F;     if (modinv_good_disc(ret, D)) return ret;
     127     1524481 :   ret = INV_W2W3;  if (modinv_good_disc(ret, D)) return ret;
     128     1524481 :   ret = INV_W2W5;  if (modinv_good_disc(ret, D)) return ret;
     129     1231083 :   ret = INV_W2W7;  if (modinv_good_disc(ret, D)) return ret;
     130     1132859 :   ret = INV_W2W13; if (modinv_good_disc(ret, D)) return ret;
     131      832846 :   ret = INV_W3W3;  if (modinv_good_disc(ret, D)) return ret;
     132      647521 :   ret = INV_W2W3E2;if (modinv_good_disc(ret, D)) return ret;
     133      575603 :   ret = INV_W3W5;  if (modinv_good_disc(ret, D)) return ret;
     134      575477 :   ret = INV_W3W7;  if (modinv_good_disc(ret, D)) return ret;
     135      507647 :   ret = INV_W3W13; if (modinv_good_disc(ret, D)) return ret;
     136      507647 :   ret = INV_W2W5E2;if (modinv_good_disc(ret, D)) return ret;
     137      491463 :   ret = INV_F3;    if (modinv_good_disc(ret, D)) return ret;
     138      461419 :   ret = INV_W2W7E2;if (modinv_good_disc(ret, D)) return ret;
     139      374248 :   ret = INV_W5W7;  if (modinv_good_disc(ret, D)) return ret;
     140      306495 :   ret = INV_W3W3E2;if (modinv_good_disc(ret, D)) return ret;
     141      306495 :   ret = INV_G2;    if (modinv_good_disc(ret, D)) return ret;
     142      159369 :   return INV_J;
     143             : }
     144             : 
     145             : INLINE long
     146       37541 : modinv_sparse_factor(long inv)
     147             : {
     148       37541 :   switch (inv) {
     149             :   case INV_G2:
     150             :   case INV_F8:
     151             :   case INV_W3W5:
     152             :   case INV_W2W5E2:
     153             :   case INV_W3W3E2:
     154        4855 :     return 3;
     155             :   case INV_F:
     156         625 :     return 24;
     157             :   case INV_F2:
     158             :   case INV_W2W3:
     159         357 :     return 12;
     160             :   case INV_F3:
     161         112 :     return 8;
     162             :   case INV_F4:
     163             :   case INV_W2W3E2:
     164             :   case INV_W2W5:
     165             :   case INV_W3W3:
     166        1553 :     return 6;
     167             :   case INV_W2W7:
     168        1046 :     return 4;
     169             :   case INV_W2W7E2:
     170             :   case INV_W2W13:
     171             :   case INV_W3W7:
     172        2873 :     return 2;
     173             :   }
     174       26120 :   return 1;
     175             : }
     176             : 
     177             : #define IQ_FILTER_1MOD3 1
     178             : #define IQ_FILTER_2MOD3 2
     179             : #define IQ_FILTER_1MOD4 4
     180             : #define IQ_FILTER_3MOD4 8
     181             : 
     182             : INLINE long
     183       12530 : modinv_pfilter(long inv)
     184             : {
     185       12530 :   switch (inv) {
     186             :   case INV_G2:
     187             :   case INV_W3W3:
     188             :   case INV_W3W3E2:
     189             :   case INV_W3W5:
     190             :   case INV_W2W5:
     191             :   case INV_W2W3E2:
     192             :   case INV_W2W5E2:
     193             :   case INV_W5W7:
     194             :   case INV_W3W13:
     195        2796 :     return IQ_FILTER_1MOD3; /* ensure unique cube roots */
     196             :   case INV_W2W7:
     197             :   case INV_F3:
     198         529 :     return IQ_FILTER_1MOD4; /* ensure at most two 4th/8th roots */
     199             :   case INV_F:
     200             :   case INV_F2:
     201             :   case INV_F4:
     202             :   case INV_F8:
     203             :   case INV_W2W3:
     204             :     /* Ensure unique cube roots and at most two 4th/8th roots */
     205         972 :     return IQ_FILTER_1MOD3 | IQ_FILTER_1MOD4;
     206             :   }
     207        8233 :   return 0;
     208             : }
     209             : 
     210             : int
     211      311858 : modinv_good_prime(long inv, long p)
     212             : {
     213      311858 :   switch (inv) {
     214             :   case INV_G2:
     215             :   case INV_W2W3E2:
     216             :   case INV_W3W3:
     217             :   case INV_W3W3E2:
     218             :   case INV_W3W5:
     219             :   case INV_W2W5E2:
     220             :   case INV_W2W5:
     221       32345 :     return (p % 3) == 2;
     222             :   case INV_W2W7:
     223             :   case INV_F3:
     224       10065 :     return (p & 3) != 1;
     225             :   case INV_F2:
     226             :   case INV_F4:
     227             :   case INV_F8:
     228             :   case INV_F:
     229             :   case INV_W2W3:
     230       17061 :     return ((p % 3) == 2) && (p & 3) != 1;
     231             :   }
     232      252387 :   return 1;
     233             : }
     234             : 
     235             : /* Returns true if the prime p does not divide the conductor of D */
     236             : INLINE int
     237     3227224 : prime_to_conductor(long D, long p)
     238             : {
     239             :   long b;
     240     3227224 :   if (p > 2) return (D % (p * p));
     241     1239090 :   b = D & 0xF;
     242     1239090 :   return (b && b != 4); /* 2 divides the conductor of D <=> D=0,4 mod 16 */
     243             : }
     244             : 
     245             : INLINE GEN
     246     3227224 : red_primeform(long D, long p)
     247             : {
     248     3227224 :   pari_sp av = avma;
     249             :   GEN P;
     250     3227224 :   if (!prime_to_conductor(D, p)) return NULL;
     251     3227224 :   P = primeform_u(stoi(D), p); /* primitive since p \nmid conductor */
     252     3227224 :   return gerepileupto(av, redimag(P));
     253             : }
     254             : 
     255             : /* Computes product of primeforms over primes appearing in the prime
     256             :  * factorization of n (including multiplicity) */
     257             : GEN
     258      134554 : qfb_nform(long D, long n)
     259             : {
     260      134554 :   pari_sp av = avma;
     261      134554 :   GEN N = NULL, fa = factoru(n), P = gel(fa,1), E = gel(fa,2);
     262      134554 :   long i, l = lg(P);
     263             : 
     264      403585 :   for (i = 1; i < l; ++i)
     265             :   {
     266             :     long j, e;
     267      269031 :     GEN Q = red_primeform(D, P[i]);
     268      269031 :     if (!Q) { avma = av; return NULL; }
     269      269031 :     e = E[i];
     270      269031 :     if (i == 1) { N = Q; j = 1; } else j = 0;
     271      269031 :     for (; j < e; ++j) N = qficomp(Q, N);
     272             :   }
     273      134554 :   return gerepileupto(av, N);
     274             : }
     275             : 
     276             : INLINE int
     277     1668590 : qfb_is_two_torsion(GEN x)
     278             : {
     279     3337180 :   return equali1(gel(x,1)) || !signe(gel(x,2))
     280     3336802 :     || equalii(gel(x,1), gel(x,2)) || equalii(gel(x,1), gel(x,3));
     281             : }
     282             : 
     283             : /* Returns true iff the products p1*p2, p1*p2^-1, p1^-1*p2, and
     284             :  * p1^-1*p2^-1 are all distinct in cl(D) */
     285             : INLINE int
     286      230217 : qfb_distinct_prods(long D, long p1, long p2)
     287             : {
     288             :   GEN P1, P2;
     289             : 
     290      230217 :   P1 = red_primeform(D, p1);
     291      230217 :   if (!P1) return 0;
     292      230217 :   P1 = qfisqr(P1);
     293             : 
     294      230217 :   P2 = red_primeform(D, p2);
     295      230217 :   if (!P2) return 0;
     296      230217 :   P2 = qfisqr(P2);
     297             : 
     298      230217 :   return !(equalii(gel(P1,1), gel(P2,1)) && absequalii(gel(P1,2), gel(P2,2)));
     299             : }
     300             : 
     301             : /* By Corollary 3.1 of Enge-Schertz Constructing elliptic curves over finite
     302             :  * fields using double eta-quotients, we need p1 != p2 to both be non-inert
     303             :  * and prime to the conductor, and if p1=p2=p we want p split and prime to the
     304             :  * conductor. We exclude the case that p1=p2 divides the conductor, even
     305             :  * though this does yield class invariants */
     306             : INLINE int
     307     5268947 : modinv_double_eta_good_disc(long D, long inv)
     308             : {
     309     5268947 :   pari_sp av = avma;
     310             :   GEN P;
     311             :   long i1, i2, p1, p2, N;
     312             : 
     313     5268947 :   N = modinv_degree(&p1, &p2, inv);
     314     5268947 :   if (! N) return 0;
     315     5268947 :   i1 = kross(D, p1);
     316     5268947 :   if (i1 < 0) return 0;
     317             :   /* Exclude ramified case for w_{p,p} */
     318     2385561 :   if (p1 == p2 && !i1) return 0;
     319     2385561 :   i2 = kross(D, p2);
     320     2385561 :   if (i2 < 0) return 0;
     321             :   /* this also verifies that p1 is prime to the conductor */
     322     1357523 :   P = red_primeform(D, p1);
     323     1357523 :   if (!P || gequal1(gel(P,1)) /* don't allow p1 to be principal */
     324             :       /* if p1 is unramified, require it to have order > 2 */
     325     1356753 :       || (i1 && qfb_is_two_torsion(P))) { avma = av; return 0; }
     326     1356179 :   if (p1 == p2)
     327             :   { /* if p1=p2 we need p1*p1 to be distinct from its inverse */
     328      215943 :     int ok = ! qfb_is_two_torsion(qfisqr(P));
     329      215943 :     avma = av; return ok;
     330             :   }
     331             :   /* this also verifies that p2 is prime to the conductor */
     332     1140236 :   P = red_primeform(D, p2);
     333     1140236 :   if (!P || gequal1(gel(P,1)) /* don't allow p2 to be principal */
     334             :       /* if p2 is unramified, require it to have order > 2 */
     335     1140110 :       || (i2 && qfb_is_two_torsion(P))) { avma = av; return 0; }
     336     1139039 :   avma = av;
     337             : 
     338             :   /* if p1 and p2 are split, we also require p1*p2, p1*p2^-1, p1^-1*p2,
     339             :    * and p1^-1*p2^-1 to be distinct */
     340     1139039 :   if (i1>0 && i2>0 && !qfb_distinct_prods(D, p1, p2)) { avma = av; return 0; }
     341     1136154 :   if (!i1 && !i2) {
     342             :     /* if both p1 and p2 are ramified, make sure their product is not
     343             :      * principal */
     344      134197 :     P = qfb_nform(D, N);
     345      134197 :     if (equali1(gel(P,1))) { avma = av; return 0; }
     346      133987 :     avma = av;
     347             :   }
     348     1135944 :   return 1;
     349             : }
     350             : 
     351             : /* Assumes D is a good discriminant for inv, which implies that the
     352             :  * level is prime to the conductor */
     353             : long
     354         483 : modinv_ramified(long D, long inv)
     355             : {
     356         483 :   long p1, p2, N = modinv_degree(&p1, &p2, inv);
     357         483 :   if (N <= 1) return 0;
     358         483 :   return !(D % p1) && !(D % p2);
     359             : }
     360             : 
     361             : int
     362    14550590 : modinv_good_disc(long inv, long D)
     363             : {
     364    14550590 :   switch (inv) {
     365             :   case INV_J:
     366      688754 :     return 1;
     367             :   case INV_G2:
     368      461531 :     return !!(D % 3);
     369             :   case INV_F3:
     370      499555 :     return (-D & 7) == 7;
     371             :   case INV_F:
     372             :   case INV_F2:
     373             :   case INV_F4:
     374             :   case INV_F8:
     375     2051159 :     return ((-D & 7) == 7) && (D % 3);
     376             :   case INV_W3W5:
     377      618219 :     return (D % 3) && modinv_double_eta_good_disc(D, inv);
     378             :   case INV_W3W3E2:
     379      333578 :     return (D % 3) && modinv_double_eta_good_disc(D, inv);
     380             :   case INV_W3W3:
     381      883239 :     return (D & 1) && (D % 3) && modinv_double_eta_good_disc(D, inv);
     382             :   case INV_W2W3E2:
     383      663404 :     return (D % 3) && modinv_double_eta_good_disc(D, inv);
     384             :   case INV_W2W3:
     385     1545145 :     return ((-D & 7) == 7) && (D % 3) && modinv_double_eta_good_disc(D, inv);
     386             :   case INV_W2W5:
     387     1572109 :     return ((-D % 80) != 20) && (D % 3) && modinv_double_eta_good_disc(D, inv);
     388             :   case INV_W2W5E2:
     389      545727 :     return (D % 3) && modinv_double_eta_good_disc(D, inv);
     390             :   case INV_W2W7E2:
     391      552699 :     return ((-D % 112) != 84) && modinv_double_eta_good_disc(D, inv);
     392             :   case INV_W2W7:
     393     1316935 :     return ((-D & 7) == 7) && modinv_double_eta_good_disc(D, inv);
     394             :   case INV_W2W13:
     395     1189741 :     return ((-D % 208) != 52) && modinv_double_eta_good_disc(D, inv);
     396             :   case INV_W3W7:
     397      662984 :     return (D & 1) && (-D % 21) && modinv_double_eta_good_disc(D, inv);
     398             :   case INV_W5W7: /* NB: This is a guess; avs doesn't have an entry */
     399      448567 :     return (D % 3) && modinv_double_eta_good_disc(D, inv);
     400             :   case INV_W3W13: /* NB: This is a guess; avs doesn't have an entry */
     401      517244 :     return (D & 1) && (D % 3) && modinv_double_eta_good_disc(D, inv);
     402             :   }
     403           0 :   pari_err_BUG("modinv_good_discriminant");
     404             :   return 0;/*LCOV_EXCL_LINE*/
     405             : }
     406             : 
     407             : int
     408         672 : modinv_is_Weber(long inv)
     409             : {
     410           0 :   return inv == INV_F || inv == INV_F2 || inv == INV_F3 || inv == INV_F4
     411         672 :     || inv == INV_F8;
     412             : }
     413             : 
     414             : int
     415      198701 : modinv_is_double_eta(long inv)
     416             : {
     417      198701 :   switch (inv) {
     418             :   case INV_W2W3:
     419             :   case INV_W2W3E2:
     420             :   case INV_W2W5:
     421             :   case INV_W2W5E2:
     422             :   case INV_W2W7:
     423             :   case INV_W2W7E2:
     424             :   case INV_W2W13:
     425             :   case INV_W3W3:
     426             :   case INV_W3W3E2:
     427             :   case INV_W3W5:
     428             :   case INV_W3W7:
     429             :   case INV_W5W7:
     430             :   case INV_W3W13:
     431       30138 :     return 1;
     432             :   }
     433      168563 :   return 0;
     434             : }
     435             : 
     436             : /* END Code from "class_inv.h" */
     437             : 
     438             : INLINE int
     439        9637 : safe_abs_sqrt(ulong *r, ulong x, ulong p, ulong pi, ulong s2)
     440             : {
     441        9637 :   if (krouu(x, p) == -1)
     442             :   {
     443        4627 :     if (p%4 == 1) return 0;
     444        4627 :     x = Fl_neg(x, p);
     445             :   }
     446        9637 :   *r = Fl_sqrt_pre_i(x, s2, p, pi);
     447        9637 :   return 1;
     448             : }
     449             : 
     450             : INLINE int
     451        4988 : eighth_root(ulong *r, ulong x, ulong p, ulong pi, ulong s2)
     452             : {
     453             :   ulong s;
     454        4988 :   if (krouu(x, p) == -1) return 0;
     455        2631 :   s = Fl_sqrt_pre_i(x, s2, p, pi);
     456        2631 :   return safe_abs_sqrt(&s, s, p, pi, s2) && safe_abs_sqrt(r, s, p, pi, s2);
     457             : }
     458             : 
     459             : INLINE ulong
     460        2890 : modinv_f_from_j(ulong j, ulong p, ulong pi, ulong s2, long only_residue)
     461             : {
     462        2890 :   pari_sp av = avma;
     463             :   GEN pol, r;
     464             :   long i;
     465        2890 :   ulong g2, f = ULONG_MAX;
     466             : 
     467             :   /* f^8 must be a root of X^3 - \gamma_2 X - 16 */
     468        2890 :   g2 = Fl_sqrtl_pre(j, 3, p, pi);
     469             : 
     470        2890 :   pol = mkvecsmall5(0UL, Fl_neg(16 % p, p), Fl_neg(g2, p), 0UL, 1UL);
     471        2890 :   r = Flx_roots(pol, p);
     472        5491 :   for (i = 1; i < lg(r); ++i)
     473        5491 :     if (only_residue) { if (krouu(r[i], p) != -1) { avma = av; return r[i]; } }
     474        4284 :     else if (eighth_root(&f, r[i], p, pi, s2)) { avma = av; return f; }
     475           0 :   pari_err_BUG("modinv_f_from_j");
     476             :   return 0;/*LCOV_EXCL_LINE*/
     477             : }
     478             : 
     479             : INLINE ulong
     480         357 : modinv_f3_from_j(ulong j, ulong p, ulong pi, ulong s2)
     481             : {
     482         357 :   pari_sp av = avma;
     483             :   GEN pol, r;
     484             :   long i;
     485         357 :   ulong f = ULONG_MAX;
     486             : 
     487        1071 :   pol = mkvecsmall5(0UL,
     488        1071 :       Fl_neg(4096 % p, p), Fl_sub(768 % p, j, p), Fl_neg(48 % p, p), 1UL);
     489         357 :   r = Flx_roots(pol, p);
     490         704 :   for (i = 1; i < lg(r); ++i)
     491         704 :     if (eighth_root(&f, r[i], p, pi, s2)) { avma = av; return f; }
     492           0 :   pari_err_BUG("modinv_f3_from_j");
     493             :   return 0;/*LCOV_EXCL_LINE*/
     494             : }
     495             : 
     496             : /* Return the exponent e for the double-eta "invariant" w such that
     497             :  * w^e is a class invariant.  For example w2w3^12 is a class
     498             :  * invariant, so double_eta_exponent(INV_W2W3) is 12 and
     499             :  * double_eta_exponent(INV_W2W3E2) is 6. */
     500             : INLINE ulong
     501       52527 : double_eta_exponent(long inv)
     502             : {
     503       52527 :   switch (inv) {
     504             :   case INV_W2W3:
     505        2517 :     return 12;
     506             :   case INV_W2W3E2:
     507             :   case INV_W2W5:
     508             :   case INV_W3W3:
     509       12445 :     return 6;
     510             :   case INV_W2W7:
     511        9451 :     return 4;
     512             :   case INV_W3W5:
     513             :   case INV_W2W5E2:
     514             :   case INV_W3W3E2:
     515        5805 :     return 3;
     516             :   case INV_W2W7E2:
     517             :   case INV_W2W13:
     518             :   case INV_W3W7:
     519       15089 :     return 2;
     520             :   default:
     521        7220 :     return 1;
     522             :   }
     523             : }
     524             : 
     525             : INLINE ulong
     526          42 : weber_exponent(long inv)
     527             : {
     528          42 :   switch (inv)
     529             :   {
     530          35 :   case INV_F:  return 24;
     531           0 :   case INV_F2: return 12;
     532           7 :   case INV_F3: return 8;
     533           0 :   case INV_F4: return 6;
     534           0 :   case INV_F8: return 3;
     535           0 :   default:     return 1;
     536             :   }
     537             : }
     538             : 
     539             : INLINE ulong
     540       28591 : double_eta_power(long inv, ulong w, ulong p, ulong pi)
     541             : {
     542       28591 :   return Fl_powu_pre(w, double_eta_exponent(inv), p, pi);
     543             : }
     544             : 
     545             : static GEN
     546         147 : double_eta_raw_to_Fp(GEN f, GEN p)
     547             : {
     548         147 :   GEN u = FpX_red(RgV_to_RgX(gel(f,1), 0), p);
     549         147 :   GEN v = FpX_red(RgV_to_RgX(gel(f,2), 0), p);
     550         147 :   return mkvec3(u, v, gel(f,3));
     551             : }
     552             : 
     553             : /* Given a root x of polclass(D, inv) modulo N, returns a root of polclass(D,0)
     554             :  * modulo N by plugging x to a modular polynomial. For double-eta quotients,
     555             :  * this is done by plugging x into the modular polynomial Phi(INV_WpWq, j)
     556             :  * Enge, Morain 2013: Generalised Weber Functions. */
     557             : GEN
     558         924 : Fp_modinv_to_j(GEN x, long inv, GEN p)
     559             : {
     560         924 :   switch(inv)
     561             :   {
     562         336 :     case INV_J: return Fp_red(x, p);
     563         399 :     case INV_G2: return Fp_powu(x, 3, p);
     564             :     case INV_F: case INV_F2: case INV_F3: case INV_F4: case INV_F8:
     565             :     {
     566          42 :       GEN xe = Fp_powu(x, weber_exponent(inv), p);
     567          42 :       return Fp_div(Fp_powu(subiu(xe, 16), 3, p), xe, p);
     568             :     }
     569             :     default:
     570         147 :     if (modinv_is_double_eta(inv))
     571             :     {
     572         147 :       GEN xe = Fp_powu(x, double_eta_exponent(inv), p);
     573         147 :       GEN uvk = double_eta_raw_to_Fp(double_eta_raw(inv), p);
     574         147 :       GEN J0 = FpX_eval(gel(uvk,1), xe, p);
     575         147 :       GEN J1 = FpX_eval(gel(uvk,2), xe, p);
     576         147 :       GEN J2 = Fp_pow(xe, gel(uvk,3), p);
     577         147 :       GEN phi = mkvec3(J0, J1, J2);
     578         147 :       return FpX_oneroot(RgX_to_FpX(RgV_to_RgX(phi,1), p),p);
     579             :     }
     580             :     pari_err_BUG("Fp_modinv_to_j"); return NULL;/* LCOV_EXCL_LINE */
     581             :   }
     582             : }
     583             : 
     584             : /* Assuming p = 2 (mod 3) and p = 3 (mod 4): if the two 12th roots of
     585             :  * x (mod p) exist, set *r to one of them and return 1, otherwise
     586             :  * return 0 (without touching *r). */
     587             : INLINE int
     588         928 : twelth_root(ulong *r, ulong x, ulong p, ulong pi, ulong s2)
     589             : {
     590             :   register ulong t;
     591             : 
     592         928 :   t = Fl_sqrtl_pre(x, 3, p, pi);
     593         928 :   if (krouu(t, p) == -1)
     594          60 :     return 0;
     595         868 :   t = Fl_sqrt_pre_i(t, s2, p, pi);
     596         868 :   return safe_abs_sqrt(r, t, p, pi, s2);
     597             : }
     598             : 
     599             : INLINE int
     600        5221 : sixth_root(ulong *r, ulong x, ulong p, ulong pi, ulong s2)
     601             : {
     602             :   register ulong t;
     603             : 
     604        5221 :   t = Fl_sqrtl_pre(x, 3, p, pi);
     605        5221 :   if (krouu(t, p) == -1)
     606         209 :     return 0;
     607        5012 :   *r = Fl_sqrt_pre_i(t, s2, p, pi);
     608        5012 :   return 1;
     609             : }
     610             : 
     611             : INLINE int
     612        3813 : fourth_root(ulong *r, ulong x, ulong p, ulong pi, ulong s2)
     613             : {
     614             :   register ulong s;
     615        3813 :   if (krouu(x, p) == -1)
     616         306 :     return 0;
     617        3507 :   s = Fl_sqrt_pre_i(x, s2, p, pi);
     618        3507 :   return safe_abs_sqrt(r, s, p, pi, s2);
     619             : }
     620             : 
     621             : INLINE int
     622       23789 : double_eta_root(long inv, ulong *r, ulong w, ulong p, ulong pi, ulong s2)
     623             : {
     624       23789 :   switch (double_eta_exponent(inv)) {
     625         928 :   case 12: return twelth_root(r, w, p, pi, s2);
     626        5221 :   case 6: return sixth_root(r, w, p, pi, s2);
     627        3813 :   case 4: return fourth_root(r, w, p, pi, s2);
     628        2625 :   case 3: *r = Fl_sqrtl_pre(w, 3, p, pi); return 1;
     629        8117 :   case 2: return krouu(w, p) != -1 && !!(*r = Fl_sqrt_pre_i(w, s2, p, pi));
     630        3084 :   case 1: *r = w; return 1;
     631             :   }
     632             :   pari_err_BUG("double_eta_root"); return 0;/*LCOV_EXCL_LINE*/
     633             : }
     634             : 
     635             : /* F = double_eta_Fl(inv, p) */
     636             : static GEN
     637       40099 : Flx_double_eta_xpoly(GEN F, ulong j, ulong p, ulong pi)
     638             : {
     639       40099 :   GEN u = gel(F,1), v = gel(F,2), w;
     640       40099 :   long i, k = itos(gel(F,3)), lu = lg(u), lv = lg(v), lw = lu + 1;
     641             : 
     642       40099 :   w = cgetg(lw, t_VECSMALL); /* lu >= max(lv,k) */
     643       40097 :   w[1] = 0; /* variable number */
     644     1133688 :   for (i = 1; i < lv; i++)
     645     1093588 :     uel(w, i + 1) = Fl_add(uel(u, i), Fl_mul_pre(j, uel(v, i), p, pi), p);
     646       80200 :   for (     ; i < lu; i++)
     647       40100 :     uel(w, i + 1) = uel(u, i);
     648       40100 :   uel(w, k + 2) = Fl_add(uel(w, k + 2), Fl_sqr_pre(j, p, pi), p);
     649       40100 :   return Flx_renormalize(w, lw);
     650             : }
     651             : 
     652             : /* F = double_eta_Fl(inv, p) */
     653             : static GEN
     654       28591 : Flx_double_eta_jpoly(GEN F, ulong x, ulong p, ulong pi)
     655             : {
     656       28591 :   pari_sp av = avma;
     657       28591 :   GEN u = gel(F,1), v = gel(F,2), xs;
     658       28591 :   long k = itos(gel(F,3));
     659             :   ulong a, b, c;
     660             : 
     661             :   /* u is always longest and the length is bigger than k */
     662       28591 :   xs = Fl_powers_pre(x, lg(u) - 1, p, pi);
     663       28591 :   c = Flv_dotproduct_pre(u, xs, p, pi);
     664       28592 :   b = Flv_dotproduct_pre(v, xs, p, pi);
     665       28592 :   a = uel(xs, k + 1);
     666       28592 :   avma = av;
     667       28592 :   return mkvecsmall4(0, c, b, a);
     668             : }
     669             : 
     670             : /* reduce F = double_eta_raw(inv) mod p */
     671             : static GEN
     672       28753 : double_eta_raw_to_Fl(GEN f, ulong p)
     673             : {
     674       28753 :   GEN u = ZV_to_Flv(gel(f,1), p);
     675       28754 :   GEN v = ZV_to_Flv(gel(f,2), p);
     676       28754 :   return mkvec3(u, v, gel(f,3));
     677             : }
     678             : /* double_eta_raw(inv) mod p */
     679             : static GEN
     680       22846 : double_eta_Fl(long inv, ulong p)
     681       22846 : { return double_eta_raw_to_Fl(double_eta_raw(inv), p); }
     682             : 
     683             : /* Go through roots of Psi(X,j) until one has an double_eta_exponent(inv)-th
     684             :  * root, and return that root. F = double_eta_Fl(inv,p) */
     685             : INLINE ulong
     686        5591 : modinv_double_eta_from_j(GEN F, long inv, ulong j, ulong p, ulong pi, ulong s2)
     687             : {
     688        5591 :   pari_sp av = avma;
     689             :   long i;
     690        5591 :   ulong f = ULONG_MAX;
     691        5591 :   GEN a = Flx_double_eta_xpoly(F, j, p, pi);
     692        5592 :   a = Flx_roots(a, p);
     693        6536 :   for (i = 1; i < lg(a); ++i)
     694        6536 :     if (double_eta_root(inv, &f, uel(a, i), p, pi, s2)) break;
     695        5592 :   if (i == lg(a)) pari_err_BUG("modinv_double_eta_from_j");
     696        5592 :   avma = av; return f;
     697             : }
     698             : 
     699             : /* assume j1 != j2 */
     700             : static long
     701       11662 : modinv_double_eta_from_2j(
     702             :   ulong *r, long inv, ulong j1, ulong j2, ulong p, ulong pi, ulong s2)
     703             : {
     704       11662 :   pari_sp av = avma;
     705       11662 :   GEN f, g, d, F = double_eta_Fl(inv, p);
     706             : 
     707       11662 :   f = Flx_double_eta_xpoly(F, j1, p, pi);
     708       11662 :   g = Flx_double_eta_xpoly(F, j2, p, pi);
     709       11662 :   d = Flx_gcd(f, g, p);
     710             :   /* NB: Morally the next conditional should be written as follows, but,
     711             :    * because of the case when j1 or j2 may not have the correct endomorphism
     712             :    * ring, we need to use the less strict conditional underneath */
     713             : #if 0
     714             :   if (degpol(d) != 1
     715             :       || (*r = Flx_oneroot(d, p)) == p
     716             :       || ! double_eta_root(inv, r, *r, p, pi, s2))
     717             :     pari_err_BUG("modinv_double_eta_from_2j");
     718             : #endif
     719       11662 :   if (degpol(d) > 2
     720       11662 :       || (*r = Flx_oneroot(d, p)) == p
     721       11662 :       || ! double_eta_root(inv, r, *r, p, pi, s2)) return 0;
     722       11647 :   avma = av; return 1;
     723             : }
     724             : 
     725             : long
     726       11718 : modfn_unambiguous_root(ulong *r, long inv, ulong j0, norm_eqn_t ne, GEN jdb)
     727             : {
     728       11718 :   pari_sp av = avma;
     729       11718 :   long p1, p2, v = ne->v, p1_depth;
     730       11718 :   ulong j1, p = ne->p, pi = ne->pi, s2 = ne->s2;
     731             :   GEN phi;
     732             : 
     733       11718 :   (void) modinv_degree(&p1, &p2, inv);
     734       11718 :   p1_depth = u_lval(v, p1);
     735             : 
     736       11718 :   phi = polmodular_db_getp(jdb, p1, p);
     737       11718 :   if (!next_surface_nbr(&j1, phi, p1, p1_depth, j0, NULL, p, pi))
     738           0 :     pari_err_BUG("modfn_unambiguous_root");
     739       11718 :   if (p2 == p1) {
     740        1813 :     if (!next_surface_nbr(&j1, phi, p1, p1_depth, j1, &j0, p, pi))
     741           0 :       pari_err_BUG("modfn_unambiguous_root");
     742             :   } else {
     743        9905 :     long p2_depth = u_lval(v, p2);
     744        9905 :     phi = polmodular_db_getp(jdb, p2, p);
     745        9905 :     if (!next_surface_nbr(&j1, phi, p2, p2_depth, j1, NULL, p, pi))
     746           0 :       pari_err_BUG("modfn_unambiguous_root");
     747             :   }
     748       11718 :   avma = av;
     749       11718 :   return j1 != j0 && modinv_double_eta_from_2j(r, inv, j0, j1, p, pi, s2);
     750             : }
     751             : 
     752             : ulong
     753      160675 : modfn_root(ulong j, norm_eqn_t ne, long inv)
     754             : {
     755      160675 :   ulong f, p = ne->p, pi = ne->pi, s2 = ne->s2;
     756      160675 :   switch (inv) {
     757      149672 :     case INV_J:  return j;
     758        7756 :     case INV_G2: return Fl_sqrtl_pre(j, 3, p, pi);
     759        1525 :     case INV_F:  return modinv_f_from_j(j, p, pi, s2, 0);
     760             :     case INV_F2:
     761         196 :       f = modinv_f_from_j(j, p, pi, s2, 0);
     762         196 :       return Fl_sqr_pre(f, p, pi);
     763         357 :     case INV_F3: return modinv_f3_from_j(j, p, pi, s2);
     764             :     case INV_F4:
     765         553 :       f = modinv_f_from_j(j, p, pi, s2, 0);
     766         553 :       return Fl_sqr_pre(Fl_sqr_pre(f, p, pi), p, pi);
     767         616 :     case INV_F8: return modinv_f_from_j(j, p, pi, s2, 1);
     768             :   }
     769           0 :   if (modinv_is_double_eta(inv))
     770             :   {
     771           0 :     pari_sp av = avma;
     772           0 :     ulong f = modinv_double_eta_from_j(double_eta_Fl(inv,p), inv, j, p, pi, s2);
     773           0 :     avma = av; return f;
     774             :   }
     775             :   pari_err_BUG("modfn_root"); return ULONG_MAX;/*LCOV_EXCL_LINE*/
     776             : }
     777             : 
     778             : INLINE ulong
     779        1858 : modinv_j_from_f(ulong x, ulong n, ulong p, ulong pi)
     780             : { /* If x satisfies (X^24 - 16)^3 - X^24 * j = 0
     781             :    * then j = (x^24 - 16)^3 / x^24 */
     782        1858 :   ulong x24 = Fl_powu_pre(x, 24 / n, p, pi);
     783        1858 :   return Fl_div(Fl_powu_pre(Fl_sub(x24, 16 % p, p), 3, p, pi), x24, p);
     784             : }
     785             : 
     786             : /* TODO: Check whether I can use this to refactor something
     787             :  * F = double_eta_raw(inv) */
     788             : long
     789        5908 : modinv_j_from_2double_eta(
     790             :   GEN F, long inv, ulong *j, ulong x0, ulong x1, ulong p, ulong pi)
     791             : {
     792             :   GEN f, g, d;
     793             : 
     794        5908 :   x0 = double_eta_power(inv, x0, p, pi);
     795        5908 :   x1 = double_eta_power(inv, x1, p, pi);
     796        5908 :   F = double_eta_raw_to_Fl(F, p);
     797        5908 :   f = Flx_double_eta_jpoly(F, x0, p, pi);
     798        5908 :   g = Flx_double_eta_jpoly(F, x1, p, pi);
     799        5908 :   d = Flx_gcd(f, g, p);
     800        5908 :   if (degpol(d) > 1)
     801           0 :     pari_err_BUG("modinv_j_from_2double_eta");
     802        5908 :   else if (degpol(d) < 1)
     803        2134 :     return 0;
     804        3774 :   if (j) *j = Flx_deg1_root(d, p);
     805        3774 :   return 1;
     806             : }
     807             : 
     808             : INLINE ulong
     809       53393 : modfn_preimage(ulong x, norm_eqn_t ne, long inv)
     810             : {
     811       53393 :   ulong p = ne->p, pi = ne->pi;
     812       53393 :   switch (inv) {
     813       45678 :     case INV_J:  return x;
     814        5857 :     case INV_G2: return Fl_powu_pre(x, 3, p, pi);
     815             :     /* NB: could replace these with a single call modinv_j_from_f(x,inv,p,pi)
     816             :      * but avoid the dependence on the actual value of inv */
     817         654 :     case INV_F:  return modinv_j_from_f(x, 1, p, pi);
     818         196 :     case INV_F2: return modinv_j_from_f(x, 2, p, pi);
     819         168 :     case INV_F3: return modinv_j_from_f(x, 3, p, pi);
     820         392 :     case INV_F4: return modinv_j_from_f(x, 4, p, pi);
     821         448 :     case INV_F8: return modinv_j_from_f(x, 8, p, pi);
     822             :   }
     823             :   /* NB: This function should never be called if modinv_double_eta(inv) is
     824             :    * true */
     825             :   pari_err_BUG("modfn_preimage"); return ULONG_MAX;/*LCOV_EXCL_LINE*/
     826             : }
     827             : 
     828             : /**
     829             :  * SECTION: class group bb_group.
     830             :  */
     831             : 
     832             : INLINE GEN
     833      168744 : mkqfis(long a, long b, long c)
     834             : {
     835      168744 :   retmkqfi(stoi(a), stoi(b), stoi(c));
     836             : }
     837             : 
     838             : /**
     839             :  * SECTION: Fixed-length dot-product-like functions on Fl's with
     840             :  * precomputed inverse.
     841             :  */
     842             : 
     843             : /* Computes x0y1 + y0x1 (mod p); assumes p < 2^63. */
     844             : INLINE ulong
     845    47742082 : Fl_addmul2(
     846             :   ulong x0, ulong x1, ulong y0, ulong y1,
     847             :   ulong p, ulong pi)
     848             : {
     849    47742082 :   return Fl_addmulmul_pre(x0,y1,y0,x1,p,pi);
     850             : }
     851             : 
     852             : /* Computes x0y2 + x1y1 + x2y0 (mod p); assumes p < 2^62. */
     853             : INLINE ulong
     854     8755894 : Fl_addmul3(
     855             :   ulong x0, ulong x1, ulong x2, ulong y0, ulong y1, ulong y2,
     856             :   ulong p, ulong pi)
     857             : {
     858             :   ulong l0, l1, h0, h1;
     859             :   LOCAL_OVERFLOW;
     860             :   LOCAL_HIREMAINDER;
     861     8755894 :   l0 = mulll(x0, y2); h0 = hiremainder;
     862     8755894 :   l1 = mulll(x1, y1); h1 = hiremainder;
     863     8755894 :   l1 = addll(l0, l1); h1 = addllx(h0, h1);
     864     8755894 :   l0 = mulll(x2, y0); h0 = hiremainder;
     865     8755894 :   l1 = addll(l0, l1); h1 = addllx(h0, h1);
     866     8755894 :   return remll_pre(h1, l1, p, pi);
     867             : }
     868             : 
     869             : /* Computes x0y3 + x1y2 + x2y1 + x3y0 (mod p); assumes p < 2^62. */
     870             : INLINE ulong
     871     3961128 : Fl_addmul4(
     872             :   ulong x0, ulong x1, ulong x2, ulong x3,
     873             :   ulong y0, ulong y1, ulong y2, ulong y3,
     874             :   ulong p, ulong pi)
     875             : {
     876             :   ulong l0, l1, h0, h1;
     877             :   LOCAL_OVERFLOW;
     878             :   LOCAL_HIREMAINDER;
     879     3961128 :   l0 = mulll(x0, y3); h0 = hiremainder;
     880     3961128 :   l1 = mulll(x1, y2); h1 = hiremainder;
     881     3961128 :   l1 = addll(l0, l1); h1 = addllx(h0, h1);
     882     3961128 :   l0 = mulll(x2, y1); h0 = hiremainder;
     883     3961128 :   l1 = addll(l0, l1); h1 = addllx(h0, h1);
     884     3961128 :   l0 = mulll(x3, y0); h0 = hiremainder;
     885     3961128 :   l1 = addll(l0, l1); h1 = addllx(h0, h1);
     886     3961128 :   return remll_pre(h1, l1, p, pi);
     887             : }
     888             : 
     889             : /* Computes x0y4 + x1y3 + x2y2 + x3y1 + x4y0 (mod p); assumes p < 2^62. */
     890             : INLINE ulong
     891    19727636 : Fl_addmul5(
     892             :   ulong x0, ulong x1, ulong x2, ulong x3, ulong x4,
     893             :   ulong y0, ulong y1, ulong y2, ulong y3, ulong y4,
     894             :   ulong p, ulong pi)
     895             : {
     896             :   ulong l0, l1, h0, h1;
     897             :   LOCAL_OVERFLOW;
     898             :   LOCAL_HIREMAINDER;
     899    19727636 :   l0 = mulll(x0, y4); h0 = hiremainder;
     900    19727636 :   l1 = mulll(x1, y3); h1 = hiremainder;
     901    19727636 :   l1 = addll(l0, l1); h1 = addllx(h0, h1);
     902    19727636 :   l0 = mulll(x2, y2); h0 = hiremainder;
     903    19727636 :   l1 = addll(l0, l1); h1 = addllx(h0, h1);
     904    19727636 :   l0 = mulll(x3, y1); h0 = hiremainder;
     905    19727636 :   l1 = addll(l0, l1); h1 = addllx(h0, h1);
     906    19727636 :   l0 = mulll(x4, y0); h0 = hiremainder;
     907    19727636 :   l1 = addll(l0, l1); h1 = addllx(h0, h1);
     908    19727636 :   return remll_pre(h1, l1, p, pi);
     909             : }
     910             : 
     911             : /* A polmodular database for a given class invariant consists of a t_VEC whose
     912             :  * L-th entry is 0 or a GEN pointing to Phi_L.  This function produces a pair
     913             :  * of databases corresponding to the j-invariant and inv */
     914             : GEN
     915       15225 : polmodular_db_init(long inv)
     916             : {
     917             :   GEN res, jdb, fdb;
     918             :   enum { DEFAULT_MODPOL_DB_LEN = 32 };
     919             : 
     920       15225 :   res = cgetg_block(2 + 1, t_VEC);
     921       15225 :   jdb = zerovec_block(DEFAULT_MODPOL_DB_LEN);
     922       15225 :   gel(res, 1) = jdb;
     923       15225 :   if (inv != INV_J) {
     924         700 :     fdb = zerovec_block(DEFAULT_MODPOL_DB_LEN);
     925         700 :     gel(res, 2) = fdb;
     926             :   } else {
     927       14525 :     gel(res, 2) = gen_0;
     928             :   }
     929       15225 :   return res;
     930             : }
     931             : 
     932             : void
     933       20738 : polmodular_db_add_level(GEN *DB, long L, long inv)
     934             : {
     935             :   long max_L;
     936             :   GEN db;
     937             : 
     938       20738 :   if (inv == INV_J) {
     939       17463 :     db = gel(*DB, 1);
     940             :   } else {
     941        3275 :     db = gel(*DB, 2);
     942        3275 :     if ( isintzero(db)) pari_err_BUG("polmodular_db_add_level");
     943             :   }
     944             : 
     945       20738 :   max_L = lg(db) - 1;
     946       20738 :   if (L > max_L) {
     947             :     GEN newdb;
     948          43 :     long i, newlen = 2 * L;
     949             : 
     950          43 :     newdb = cgetg_block(newlen + 1, t_VEC);
     951          43 :     for (i = 1; i <= max_L; ++i) gel(newdb, i) = gel(db, i);
     952          43 :     for (     ; i <= newlen; ++i) gel(newdb, i) = gen_0;
     953          43 :     killblock(db);
     954          43 :     gel(*DB, (inv == INV_J)? 1: 2) = db = newdb;
     955             :   }
     956       20738 :   if (typ(gel(db, L)) == t_INT) {
     957        6150 :     pari_sp av = avma;
     958        6150 :     GEN x = polmodular0_ZM(L, inv, NULL, NULL, 0, DB); /* may set db[L] */
     959        6150 :     GEN y = gel(db, L);
     960        6150 :     gel(db, L) = gclone(x);
     961        6150 :     if (typ(y) != t_INT) gunclone(y);
     962        6150 :     avma = av;
     963             :   }
     964       20738 : }
     965             : 
     966             : void
     967        3798 : polmodular_db_add_levels(GEN *db, long *levels, long k, long inv)
     968             : {
     969             :   long i;
     970        3798 :   for (i = 0; i < k; ++i) polmodular_db_add_level(db, levels[i], inv);
     971        3798 : }
     972             : 
     973             : GEN
     974      301586 : polmodular_db_for_inv(GEN db, long inv)
     975      301586 : { return (inv == INV_J)? gel(db,1): gel(db,2); }
     976             : 
     977             : /* TODO: Also cache modpoly mod p for most recent p (avoid repeated
     978             :  * reductions in, for example, polclass.c:oneroot_of_classpoly(). */
     979             : GEN
     980      439666 : polmodular_db_getp(GEN db, long L, ulong p)
     981             : {
     982      439666 :   GEN f = gel(db, L);
     983      439666 :   if (isintzero(f)) pari_err_BUG("polmodular_db_getp");
     984      439667 :   return ZM_to_Flm(f, p);
     985             : }
     986             : 
     987             : /**
     988             :  * SECTION: Table of discriminants to use.
     989             :  */
     990             : typedef struct {
     991             :   long GENcode0;  /* used when serializing the struct to a t_VECSMALL */
     992             :   long inv;      /* invariant */
     993             :   long L;        /* modpoly level */
     994             :   long D0;       /* fundamental discriminant */
     995             :   long D1;       /* chosen discriminant */
     996             :   long L0;       /* first generator norm */
     997             :   long L1;       /* second generator norm */
     998             :   long n1;       /* order of L0 in cl(D1) */
     999             :   long n2;       /* order of L0 in cl(D2) where D2 = L^2 D1 */
    1000             :   long dl1;      /* m such that L0^m = L in cl(D1) */
    1001             :   long dl2_0;    /* These two are (m, n) such that L0^m L1^n = form of norm L^2 in D2 */
    1002             :   long dl2_1;    /* This n is always 1 or 0. */
    1003             :   /* this part is not serialized */
    1004             :   long nprimes;  /* number of primes needed for D1 */
    1005             :   long cost;     /* cost to enumerate  subgroup of cl(L^2D): subgroup size is n2 if L1=0, 2*n2 o.w. */
    1006             :   long bits;
    1007             :   ulong *primes;
    1008             :   ulong *traces;
    1009             : } disc_info;
    1010             : 
    1011             : #define MODPOLY_MAX_DCNT    64
    1012             : 
    1013             : /* Flag for last parameter of discriminant_with_classno_at_least.
    1014             :  * Warning: ignoring the sparse factor makes everything slower by
    1015             :  * something like (sparse factor)^3. */
    1016             : #define USE_SPARSE_FACTOR 0
    1017             : #define IGNORE_SPARSE_FACTOR 1
    1018             : 
    1019             : static long
    1020             : discriminant_with_classno_at_least(disc_info Ds[MODPOLY_MAX_DCNT], long L,
    1021             :   long inv, long ignore_sparse);
    1022             : 
    1023             : /**
    1024             :  * SECTION: Hard-coded evaluation functions for modular polynomials of
    1025             :  * small level.
    1026             :  */
    1027             : 
    1028             : /* Based on phi2_eval_ff() in Sutherland's classpoly programme.
    1029             :  * Calculates Phi_2(X, j) (mod p) with 6M+7A (4 reductions, not
    1030             :  * counting those for Phi_2) */
    1031             : INLINE GEN
    1032    22459239 : Flm_Fl_phi2_evalx(GEN phi2, ulong j, ulong p, ulong pi)
    1033             : {
    1034    22459239 :   GEN res = cgetg(6, t_VECSMALL);
    1035             :   ulong j2, t1;
    1036             : 
    1037    22438707 :   res[1] = 0; /* variable name */
    1038             : 
    1039    22438707 :   j2 = Fl_sqr_pre(j, p, pi);
    1040    22495397 :   t1 = Fl_add(j, coeff(phi2, 3, 1), p);
    1041    22460168 :   t1 = Fl_addmul2(j, j2, t1, coeff(phi2, 2, 1), p, pi);
    1042    22505577 :   res[2] = Fl_add(t1, coeff(phi2, 1, 1), p);
    1043             : 
    1044    22477306 :   t1 = Fl_addmul2(j, j2, coeff(phi2, 3, 2), coeff(phi2, 2, 2), p, pi);
    1045    22513342 :   res[3] = Fl_add(t1, coeff(phi2, 2, 1), p);
    1046             : 
    1047    22489306 :   t1 = Fl_mul_pre(j, coeff(phi2, 3, 2), p, pi);
    1048    22495216 :   t1 = Fl_add(t1, coeff(phi2, 3, 1), p);
    1049    22480192 :   res[4] = Fl_sub(t1, j2, p);
    1050             : 
    1051    22466792 :   res[5] = 1;
    1052    22466792 :   return res;
    1053             : }
    1054             : 
    1055             : /* Based on phi3_eval_ff() in Sutherland's classpoly programme.
    1056             :  * Calculates Phi_3(X, j) (mod p) with 13M+13A (6 reductions, not
    1057             :  * counting those for Phi_3) */
    1058             : INLINE GEN
    1059     2921634 : Flm_Fl_phi3_evalx(GEN phi3, ulong j, ulong p, ulong pi)
    1060             : {
    1061     2921634 :   GEN res = cgetg(7, t_VECSMALL);
    1062             :   ulong j2, j3, t1;
    1063             : 
    1064     2919768 :   res[1] = 0; /* variable name */
    1065             : 
    1066     2919768 :   j2 = Fl_sqr_pre(j, p, pi);
    1067     2923866 :   j3 = Fl_mul_pre(j, j2, p, pi);
    1068             : 
    1069     2923978 :   t1 = Fl_add(j, coeff(phi3, 4, 1), p);
    1070     8767500 :   res[2] = Fl_addmul3(j, j2, j3, t1,
    1071     5845000 :                       coeff(phi3, 3, 1), coeff(phi3, 2, 1), p, pi);
    1072             : 
    1073     5853798 :   t1 = Fl_addmul3(j, j2, j3, coeff(phi3, 4, 2),
    1074     5853798 :                   coeff(phi3, 3, 2), coeff(phi3, 2, 2), p, pi);
    1075     2927209 :   res[3] = Fl_add(t1, coeff(phi3, 2, 1), p);
    1076             : 
    1077     5851590 :   t1 = Fl_addmul3(j, j2, j3, coeff(phi3, 4, 3),
    1078     5851590 :                   coeff(phi3, 3, 3), coeff(phi3, 3, 2), p, pi);
    1079     2927305 :   res[4] = Fl_add(t1, coeff(phi3, 3, 1), p);
    1080             : 
    1081     2925811 :   t1 = Fl_addmul2(j, j2, coeff(phi3, 4, 3), coeff(phi3, 4, 2), p, pi);
    1082     2926885 :   t1 = Fl_add(t1, coeff(phi3, 4, 1), p);
    1083     2925407 :   res[5] = Fl_sub(t1, j3, p);
    1084             : 
    1085     2924224 :   res[6] = 1;
    1086     2924224 :   return res;
    1087             : }
    1088             : 
    1089             : /* Based on phi5_eval_ff() in Sutherland's classpoly programme.
    1090             :  * Calculates Phi_5(X, j) (mod p) with 33M+31A (10 reductions, not
    1091             :  * counting those for Phi_5) */
    1092             : INLINE GEN
    1093     3955762 : Flm_Fl_phi5_evalx(GEN phi5, ulong j, ulong p, ulong pi)
    1094             : {
    1095     3955762 :   GEN res = cgetg(9, t_VECSMALL);
    1096             :   ulong j2, j3, j4, j5, t1;
    1097             : 
    1098     3952977 :   res[1] = 0; /* variable name */
    1099             : 
    1100     3952977 :   j2 = Fl_sqr_pre(j, p, pi);
    1101     3959075 :   j3 = Fl_mul_pre(j, j2, p, pi);
    1102     3959119 :   j4 = Fl_sqr_pre(j2, p, pi);
    1103     3959119 :   j5 = Fl_mul_pre(j, j4, p, pi);
    1104             : 
    1105     3959528 :   t1 = Fl_add(j, coeff(phi5, 6, 1), p);
    1106    15829980 :   t1 = Fl_addmul5(j, j2, j3, j4, j5, t1,
    1107     7914990 :                   coeff(phi5, 5, 1), coeff(phi5, 4, 1),
    1108     7914990 :                   coeff(phi5, 3, 1), coeff(phi5, 2, 1),
    1109             :                   p, pi);
    1110     3962661 :   res[2] = Fl_add(t1, coeff(phi5, 1, 1), p);
    1111             : 
    1112    19803080 :   t1 = Fl_addmul5(j, j2, j3, j4, j5,
    1113     7921232 :                   coeff(phi5, 6, 2), coeff(phi5, 5, 2),
    1114    11881848 :                   coeff(phi5, 4, 2), coeff(phi5, 3, 2), coeff(phi5, 2, 2),
    1115             :                   p, pi);
    1116     3963368 :   res[3] = Fl_add(t1, coeff(phi5, 2, 1), p);
    1117             : 
    1118    19806920 :   t1 = Fl_addmul5(j, j2, j3, j4, j5,
    1119     7922768 :                   coeff(phi5, 6, 3), coeff(phi5, 5, 3),
    1120    11884152 :                   coeff(phi5, 4, 3), coeff(phi5, 3, 3), coeff(phi5, 3, 2),
    1121             :                   p, pi);
    1122     3963805 :   res[4] = Fl_add(t1, coeff(phi5, 3, 1), p);
    1123             : 
    1124    19807730 :   t1 = Fl_addmul5(j, j2, j3, j4, j5,
    1125     7923092 :                   coeff(phi5, 6, 4), coeff(phi5, 5, 4),
    1126    11884638 :                   coeff(phi5, 4, 4), coeff(phi5, 4, 3), coeff(phi5, 4, 2),
    1127             :                   p, pi);
    1128     3963827 :   res[5] = Fl_add(t1, coeff(phi5, 4, 1), p);
    1129             : 
    1130    19809505 :   t1 = Fl_addmul5(j, j2, j3, j4, j5,
    1131     7923802 :                   coeff(phi5, 6, 5), coeff(phi5, 5, 5),
    1132    11885703 :                   coeff(phi5, 5, 4), coeff(phi5, 5, 3), coeff(phi5, 5, 2),
    1133             :                   p, pi);
    1134     3964012 :   res[6] = Fl_add(t1, coeff(phi5, 5, 1), p);
    1135             : 
    1136    15848672 :   t1 = Fl_addmul4(j, j2, j3, j4,
    1137     7924336 :                   coeff(phi5, 6, 5), coeff(phi5, 6, 4),
    1138     7924336 :                   coeff(phi5, 6, 3), coeff(phi5, 6, 2),
    1139             :                   p, pi);
    1140     3964332 :   t1 = Fl_add(t1, coeff(phi5, 6, 1), p);
    1141     3962339 :   res[7] = Fl_sub(t1, j5, p);
    1142             : 
    1143     3960880 :   res[8] = 1;
    1144     3960880 :   return res;
    1145             : }
    1146             : 
    1147             : GEN
    1148    36316350 : Flm_Fl_polmodular_evalx(GEN phi, long L, ulong j, ulong p, ulong pi)
    1149             : {
    1150    36316350 :   switch (L) {
    1151    22463731 :     case 2: return Flm_Fl_phi2_evalx(phi, j, p, pi);
    1152     2922173 :     case 3: return Flm_Fl_phi3_evalx(phi, j, p, pi);
    1153     3956428 :     case 5: return Flm_Fl_phi5_evalx(phi, j, p, pi);
    1154             :     default: { /* not GC clean, but gerepileupto-safe */
    1155     6974018 :       GEN j_powers = Fl_powers_pre(j, L + 1, p, pi);
    1156     6959185 :       return Flm_Flc_mul_pre_Flx(phi, j_powers, p, pi, 0);
    1157             :     }
    1158             :   }
    1159             : }
    1160             : 
    1161             : /**
    1162             :  * SECTION: Velu's formula for the codmain curve in the case of small
    1163             :  * prime base field.
    1164             :  */
    1165             : 
    1166             : INLINE ulong
    1167     1481274 : Fl_mul4(ulong x, ulong p)
    1168     1481274 : { return Fl_double(Fl_double(x, p), p); }
    1169             : 
    1170             : INLINE ulong
    1171       77206 : Fl_mul5(ulong x, ulong p)
    1172       77206 : { return Fl_add(x, Fl_mul4(x, p), p); }
    1173             : 
    1174             : INLINE ulong
    1175      740747 : Fl_mul8(ulong x, ulong p)
    1176      740747 : { return Fl_double(Fl_mul4(x, p), p); }
    1177             : 
    1178             : INLINE ulong
    1179      663591 : Fl_mul6(ulong x, ulong p)
    1180      663591 : { return Fl_sub(Fl_mul8(x, p), Fl_double(x, p), p); }
    1181             : 
    1182             : INLINE ulong
    1183       77204 : Fl_mul7(ulong x, ulong p)
    1184       77204 : { return Fl_sub(Fl_mul8(x, p), x, p); }
    1185             : 
    1186             : /* Given an elliptic curve E = [a4, a6] over F_p and a non-zero point
    1187             :  * pt on E, return the quotient E' = E/<P> = [a4_img, a6_img] */
    1188             : static void
    1189       77207 : Fle_quotient_from_kernel_generator(
    1190             :   ulong *a4_img, ulong *a6_img, ulong a4, ulong a6, GEN pt, ulong p, ulong pi)
    1191             : {
    1192       77207 :   pari_sp av = avma;
    1193       77207 :   ulong t = 0, w = 0;
    1194             :   GEN Q;
    1195             :   ulong xQ, yQ, tQ, uQ;
    1196             : 
    1197       77207 :   Q = gcopy(pt);
    1198             :   /* Note that, as L is odd, say L = 2n + 1, we necessarily have
    1199             :    * [(L - 1)/2]P = [n]P = [n - L]P = -[n + 1]P = -[(L + 1)/2]P.  This is
    1200             :    * what the condition Q[1] != xQ tests, so the loop will execute n times. */
    1201             :   do {
    1202      663526 :     xQ = uel(Q, 1);
    1203      663526 :     yQ = uel(Q, 2);
    1204             :     /* tQ = 6 xQ^2 + b2 xQ + b4
    1205             :      *    = 6 xQ^2 + 2 a4 (since b2 = 0 and b4 = 2 a4) */
    1206      663526 :     tQ = Fl_add(Fl_mul6(Fl_sqr_pre(xQ, p, pi), p), Fl_double(a4, p), p);
    1207      663503 :     uQ = Fl_add(Fl_mul4(Fl_sqr_pre(yQ, p, pi), p),
    1208             :                 Fl_mul_pre(tQ, xQ, p, pi), p);
    1209             : 
    1210      663506 :     t = Fl_add(t, tQ, p);
    1211      663469 :     w = Fl_add(w, uQ, p);
    1212      663432 :     Q = gerepileupto(av, Fle_add(pt, Q, a4, p));
    1213      663527 :   } while (uel(Q, 1) != xQ);
    1214             : 
    1215       77206 :   avma = av;
    1216             :   /* a4_img = a4 - 5 * t */
    1217       77206 :   *a4_img = Fl_sub(a4, Fl_mul5(t, p), p);
    1218             :   /* a6_img = a6 - b2 * t - 7 * w = a6 - 7 * w (since a1 = a2 = 0 ==> b2 = 0) */
    1219       77204 :   *a6_img = Fl_sub(a6, Fl_mul7(w, p), p);
    1220       77202 : }
    1221             : 
    1222             : /**
    1223             :  * SECTION: Calculation of modular polynomials.
    1224             :  */
    1225             : 
    1226             : /* Given an elliptic curve [a4, a6] over FF_p, try to find a
    1227             :  * non-trivial L-torsion point on the curve by considering n times a
    1228             :  * random point; val controls the maximum L-valuation expected of n
    1229             :  * times a random point */
    1230             : static GEN
    1231      112738 : find_L_tors_point(
    1232             :   ulong *ival,
    1233             :   ulong a4, ulong a6, ulong p, ulong pi,
    1234             :   ulong n, ulong L, ulong val)
    1235             : {
    1236      112738 :   pari_sp av = avma;
    1237             :   ulong i;
    1238             :   GEN P, Q;
    1239             :   do {
    1240      113751 :     Q = random_Flj_pre(a4, a6, p, pi);
    1241      113751 :     P = Flj_mulu_pre(Q, n, a4, p, pi);
    1242      113741 :   } while (P[3] == 0);
    1243             : 
    1244      218761 :   for (i = 0; i < val; ++i) {
    1245      183216 :     Q = Flj_mulu_pre(P, L, a4, p, pi);
    1246      183239 :     if (Q[3] == 0) break;
    1247      106033 :     P = Q;
    1248             :   }
    1249      112751 :   if (ival) *ival = i;
    1250      112751 :   return gerepilecopy(av, P);
    1251             : }
    1252             : 
    1253             : static GEN
    1254       70162 : select_curve_with_L_tors_point(
    1255             :   ulong *a4, ulong *a6,
    1256             :   ulong L, ulong j, ulong n, ulong card, ulong val,
    1257             :   norm_eqn_t ne)
    1258             : {
    1259       70162 :   pari_sp av = avma;
    1260             :   ulong A4, A4t, A6, A6t;
    1261       70162 :   ulong p = ne->p, pi = ne->pi;
    1262             :   GEN P;
    1263       70162 :   if (card % L != 0) {
    1264           0 :     pari_err_BUG("select_curve_with_L_tors_point: "
    1265             :                  "Cardinality not divisible by L");
    1266             :   }
    1267             : 
    1268       70162 :   Fl_ellj_to_a4a6(j, p, &A4, &A6);
    1269       70164 :   Fl_elltwist_disc(A4, A6, ne->T, p, &A4t, &A6t);
    1270             : 
    1271             :   /* Either E = [a4, a6] or its twist has cardinality divisible by L
    1272             :    * because of the choice of p and t earlier on.  We find out which
    1273             :    * by attempting to find a point of order L on each.  See bot p16 of
    1274             :    * Sutherland 2012. */
    1275       35540 :   while (1) {
    1276             :     ulong i;
    1277      105701 :     P = find_L_tors_point(&i, A4, A6, p, pi, n, L, val);
    1278      105710 :     if (i < val)
    1279       70170 :       break;
    1280       35540 :     avma = av;
    1281       35540 :     lswap(A4, A4t);
    1282       35540 :     lswap(A6, A6t);
    1283             :   }
    1284       70170 :   *a4 = A4;
    1285      140339 :   *a6 = A6; return gerepilecopy(av, P);
    1286             : }
    1287             : 
    1288             : /* Return 1 if the L-Sylow subgroup of the curve [a4, a6] (mod p) is
    1289             :  * cyclic, return 0 if it is not cyclic with "high" probability (I
    1290             :  * guess around 1/L^3 chance it is still cyclic when we return 0).
    1291             :  *
    1292             :  * Based on Sutherland's velu.c:velu_verify_Sylow_cyclic() in classpoly-1.0.1 */
    1293             : INLINE long
    1294       39325 : verify_L_sylow_is_cyclic(long e, ulong a4, ulong a6, ulong p, ulong pi)
    1295             : {
    1296             :   /* Number of times to try to find a point with maximal order in the
    1297             :    * L-Sylow subgroup. */
    1298             :   enum { N_RETRIES = 3 };
    1299       39325 :   pari_sp av = avma;
    1300       39325 :   long i, res = 0;
    1301             :   GEN P;
    1302       63847 :   for (i = 0; i < N_RETRIES; ++i) {
    1303       56809 :     P = random_Flj_pre(a4, a6, p, pi);
    1304       56807 :     P = Flj_mulu_pre(P, e, a4, p, pi);
    1305       56807 :     if (P[3] != 0) { res = 1; break; }
    1306             :   }
    1307       39323 :   avma = av; return res;
    1308             : }
    1309             : 
    1310             : static ulong
    1311       70169 : find_noniso_L_isogenous_curve(
    1312             :   ulong L, ulong n,
    1313             :   norm_eqn_t ne, long e, ulong val, ulong a4, ulong a6, GEN init_pt, long verify)
    1314             : {
    1315             :   pari_sp ltop, av;
    1316       70169 :   ulong p = ne->p, pi = ne->pi, j_res = 0;
    1317       70169 :   GEN pt = init_pt;
    1318       70169 :   ltop = av = avma;
    1319        7038 :   while (1) {
    1320             :     /* c. Use Velu to calculate L-isogenous curve E' = E/<P> */
    1321             :     ulong a4_img, a6_img;
    1322       77207 :     ulong z2 = Fl_sqr_pre(pt[3], p, pi);
    1323       77210 :     pt = mkvecsmall2(Fl_div(pt[1], z2, p),
    1324       77209 :                      Fl_div(pt[2], Fl_mul_pre(z2, pt[3], p, pi), p));
    1325       77207 :     Fle_quotient_from_kernel_generator(&a4_img, &a6_img,
    1326             :                                        a4, a6, pt, p, pi);
    1327             : 
    1328             :     /* d. If j(E') = j_res has a different endo ring to j(E), then
    1329             :      *    return j(E').  Otherwise, go to b. */
    1330       77202 :     if (!verify || verify_L_sylow_is_cyclic(e, a4_img, a6_img, p, pi)) {
    1331       70161 :       j_res = Fl_ellj_pre(a4_img, a6_img, p, pi);
    1332       70171 :       break;
    1333             :     }
    1334             : 
    1335             :     /* b. Generate random point P on E of order L */
    1336        7038 :     avma = av;
    1337        7038 :     pt = find_L_tors_point(NULL, a4, a6, p, pi, n, L, val);
    1338             :   }
    1339      140342 :   avma = ltop; return j_res;
    1340             : }
    1341             : 
    1342             : /* Given a prime L and a j-invariant j (mod p), return the j-invariant
    1343             :  * of a curve which has a different endomorphism ring to j and is
    1344             :  * L-isogenous to j */
    1345             : INLINE ulong
    1346       70158 : compute_L_isogenous_curve(
    1347             :   ulong L, ulong n, norm_eqn_t ne,
    1348             :   ulong j, ulong card, ulong val, long verify)
    1349             : {
    1350             :   ulong a4, a6;
    1351             :   long e;
    1352             :   GEN pt;
    1353             : 
    1354       70158 :   if (ne->p < 5 || j == 0 || j == 1728 % ne->p)
    1355           0 :     pari_err_BUG("compute_L_isogenous_curve");
    1356       70159 :   pt = select_curve_with_L_tors_point(&a4, &a6, L, j, n, card, val, ne);
    1357       70169 :   e = card / L;
    1358       70169 :   if (e * L != card) pari_err_BUG("compute_L_isogenous_curve");
    1359             : 
    1360       70169 :   return find_noniso_L_isogenous_curve(L, n, ne, e, val, a4, a6, pt, verify);
    1361             : }
    1362             : 
    1363             : INLINE GEN
    1364       32285 : get_Lsqr_cycle(const disc_info *dinfo)
    1365             : {
    1366       32285 :   long i, n1 = dinfo->n1, L = dinfo->L;
    1367       32285 :   GEN cyc = cgetg(L, t_VECSMALL);
    1368       32286 :   cyc[1] = 0;
    1369       32286 :   for (i = 2; i <= L / 2; ++i) cyc[i] = cyc[i - 1] + n1;
    1370       32286 :   if ( ! dinfo->L1) {
    1371       12804 :     for ( ; i < L; ++i) cyc[i] = cyc[i - 1] + n1;
    1372             :   } else {
    1373       19482 :     cyc[L - 1] = 2 * dinfo->n2 - n1 / 2;
    1374       19482 :     for (i = L - 2; i > L / 2; --i) cyc[i] = cyc[i + 1] - n1;
    1375             :   }
    1376       32286 :   return cyc;
    1377             : }
    1378             : 
    1379             : INLINE void
    1380      438039 : update_Lsqr_cycle(GEN cyc, const disc_info *dinfo)
    1381             : {
    1382      438039 :   long i, L = dinfo->L;
    1383      438039 :   for (i = 1; i < L; ++i) ++cyc[i];
    1384      438039 :   if (dinfo->L1 && cyc[L - 1] == 2 * dinfo->n2) {
    1385       17787 :     long n1 = dinfo->n1;
    1386       17787 :     for (i = L / 2 + 1; i < L; ++i) cyc[i] -= n1;
    1387             :   }
    1388      438039 : }
    1389             : 
    1390             : static ulong
    1391       32269 : oneroot_of_classpoly(
    1392             :   GEN hilb, GEN factu, norm_eqn_t ne, GEN jdb)
    1393             : {
    1394       32269 :   pari_sp av = avma;
    1395       32269 :   ulong j0, p = ne->p, pi = ne->pi;
    1396       32269 :   long i, nfactors = lg(gel(factu, 1)) - 1;
    1397       32269 :   GEN hilbp = ZX_to_Flx(hilb, p);
    1398             : 
    1399             :   /* TODO: Work out how to use hilb with better invariant */
    1400       32251 :   j0 = Flx_oneroot_split(hilbp, p);
    1401       32285 :   if (j0 == p) {
    1402           0 :     pari_err_BUG("oneroot_of_classpoly: "
    1403             :                  "Didn't find a root of the class polynomial");
    1404             :   }
    1405       33347 :   for (i = 1; i <= nfactors; ++i) {
    1406        1062 :     long L = gel(factu, 1)[i];
    1407        1062 :     long val = gel(factu, 2)[i];
    1408        1062 :     GEN phi = polmodular_db_getp(jdb, L, p);
    1409        1062 :     val += z_lval(ne->v, L);
    1410        1062 :     j0 = descend_volcano(phi, j0, p, pi, 0, L, val, val);
    1411        1062 :     avma = av;
    1412             :   }
    1413       32285 :   avma = av; return j0;
    1414             : }
    1415             : 
    1416             : /* TODO: Precompute the classgp_pcp_t structs and link them to dinfo */
    1417             : INLINE void
    1418       32289 : make_pcp_surface(const disc_info *dinfo, classgp_pcp_t G)
    1419             : {
    1420       32289 :   long k = 1, datalen = 3 * k;
    1421             : 
    1422       32289 :   memset(G, 0, sizeof *G);
    1423             : 
    1424       32289 :   G->_data = cgetg(datalen + 1, t_VECSMALL);
    1425       32288 :   G->L = G->_data + 1;
    1426       32288 :   G->n = G->L + k;
    1427       32288 :   G->o = G->L + k;
    1428             : 
    1429       32288 :   G->k = k;
    1430       32288 :   G->h = G->enum_cnt = dinfo->n1;
    1431       32288 :   G->L[0] = dinfo->L0;
    1432       32288 :   G->n[0] = dinfo->n1;
    1433       32288 :   G->o[0] = dinfo->n1;
    1434       32288 : }
    1435             : 
    1436             : INLINE void
    1437       32290 : make_pcp_floor(const disc_info *dinfo, classgp_pcp_t G)
    1438             : {
    1439       32290 :   long k = dinfo->L1 ? 2 : 1, datalen = 3 * k;
    1440             : 
    1441       32290 :   memset(G, 0, sizeof *G);
    1442       32290 :   G->_data = cgetg(datalen + 1, t_VECSMALL);
    1443       32290 :   G->L = G->_data + 1;
    1444       32290 :   G->n = G->L + k;
    1445       32290 :   G->o = G->L + k;
    1446             : 
    1447       32290 :   G->k = k;
    1448       32290 :   G->h = G->enum_cnt = dinfo->n2 * k;
    1449       32290 :   G->L[0] = dinfo->L0;
    1450       32290 :   G->n[0] = dinfo->n2;
    1451       32290 :   G->o[0] = dinfo->n2;
    1452       32290 :   if (dinfo->L1) {
    1453       19486 :     G->L[1] = dinfo->L1;
    1454       19486 :     G->n[1] = 2;
    1455       19486 :     G->o[1] = 2;
    1456             :   }
    1457       32290 : }
    1458             : 
    1459             : INLINE GEN
    1460       32289 : enum_volcano_surface(const disc_info *dinfo, norm_eqn_t ne, ulong j0, GEN fdb)
    1461             : {
    1462       32289 :   pari_sp av = avma;
    1463             :   classgp_pcp_t G;
    1464       32289 :   make_pcp_surface(dinfo, G);
    1465       32289 :   return gerepileupto(av, enum_roots(j0, ne, fdb, G));
    1466             : }
    1467             : 
    1468             : INLINE GEN
    1469       32290 : enum_volcano_floor(long L, norm_eqn_t ne, ulong j0_pr, GEN fdb, const disc_info *dinfo)
    1470             : {
    1471       32290 :   pari_sp av = avma;
    1472             :   /* L^2 D is the discriminant for the order R = Z + L OO. */
    1473       32290 :   long DR = L * L * ne->D;
    1474       32290 :   long R_cond = L * ne->u; /* conductor(DR); */
    1475       32290 :   long w = R_cond * ne->v;
    1476             :   /* TODO: Calculate these once and for all in polmodular0_ZM(). */
    1477             :   classgp_pcp_t G;
    1478             :   norm_eqn_t eqn;
    1479       32290 :   memcpy(eqn, ne, sizeof *ne);
    1480       32290 :   eqn->D = DR;
    1481       32290 :   eqn->u = R_cond;
    1482       32290 :   eqn->v = w;
    1483       32290 :   make_pcp_floor(dinfo, G);
    1484       32290 :   return gerepileupto(av, enum_roots(j0_pr, eqn, fdb, G));
    1485             : }
    1486             : 
    1487             : INLINE void
    1488       15910 : carray_reverse_inplace(long *arr, long n)
    1489             : {
    1490       15910 :   long lim = n>>1, i;
    1491       15910 :   --n;
    1492       15910 :   for (i = 0; i < lim; i++) lswap(arr[i], arr[n - i]);
    1493       15910 : }
    1494             : 
    1495             : INLINE void
    1496      470312 : append_neighbours(GEN rts, GEN surface_js, long njs, long L, long m, long i)
    1497             : {
    1498      470312 :   long r_idx = (((i - 1) + m) % njs) + 1; /* (i + m) % njs */
    1499      470312 :   long l_idx = smodss((i - 1) - m, njs) + 1; /* (i - m) % njs */
    1500      470302 :   rts[L] = surface_js[l_idx];
    1501      470302 :   rts[L + 1] = surface_js[r_idx];
    1502      470302 : }
    1503             : 
    1504             : INLINE GEN
    1505       34455 : roots_to_coeffs(GEN rts, ulong p, long L)
    1506             : {
    1507       34455 :   long i, k, lrts= lg(rts);
    1508       34455 :   GEN M = cgetg(L+2+1, t_MAT);
    1509      760058 :   for (i = 1; i <= L+2; ++i)
    1510      725574 :     gel(M, i) = cgetg(lrts, t_VECSMALL);
    1511      529915 :   for (i = 1; i < lrts; ++i) {
    1512      495466 :     pari_sp av = avma;
    1513      495466 :     GEN modpol = Flv_roots_to_pol(gel(rts, i), p, 0);
    1514      495431 :     for (k = 1; k <= L + 2; ++k) mael(M, k, i) = modpol[k + 1];
    1515      495431 :     avma = av;
    1516             :   }
    1517       34449 :   return M;
    1518             : }
    1519             : 
    1520             : /* NB: Assumes indices are offset at 0, not at 1 like in GENs;
    1521             :  * i.e. indices[i] will pick out v[indices[i] + 1] from v. */
    1522             : INLINE void
    1523      470334 : vecsmall_pick(GEN res, GEN v, GEN indices)
    1524             : {
    1525             :   long i;
    1526      470334 :   for (i = 1; i < lg(indices); ++i) res[i] = v[indices[i] + 1];
    1527      470334 : }
    1528             : 
    1529             : /* First element of surface_js must lie above the first element of floor_js.
    1530             :  * Reverse surface_js if it is not oriented in the same direction as floor_js */
    1531             : INLINE GEN
    1532       32284 : root_matrix(long L, const disc_info *dinfo, long njinvs, GEN surface_js,
    1533             :   GEN floor_js, ulong n, ulong card, ulong val, norm_eqn_t ne)
    1534             : {
    1535             :   pari_sp av;
    1536       32284 :   long i, m = dinfo->dl1, njs = lg(surface_js) - 1, inv = dinfo->inv, rev;
    1537       32284 :   GEN rt_mat = zero_Flm_copy(L + 1, njinvs), rts, cyc;
    1538       32285 :   av = avma;
    1539             : 
    1540       32285 :   i = 1;
    1541       32285 :   cyc = get_Lsqr_cycle(dinfo);
    1542       32286 :   rts = gel(rt_mat, i);
    1543       32286 :   vecsmall_pick(rts, floor_js, cyc);
    1544       32287 :   append_neighbours(rts, surface_js, njs, L, m, i);
    1545             : 
    1546       32285 :   i = 2;
    1547       32285 :   update_Lsqr_cycle(cyc, dinfo);
    1548       32287 :   rts = gel(rt_mat, i);
    1549       32287 :   vecsmall_pick(rts, floor_js, cyc);
    1550             : 
    1551             :   /* Fix orientation if necessary */
    1552       32289 :   if (modinv_is_double_eta(inv)) {
    1553             :     /* TODO: There is potential for refactoring between this,
    1554             :      * double_eta_initial_js and modfn_preimage. */
    1555        5592 :     pari_sp av0 = avma;
    1556        5592 :     ulong p = ne->p, pi = ne->pi, j;
    1557        5592 :     GEN F = double_eta_Fl(inv, p);
    1558        5592 :     pari_sp av = avma;
    1559        5592 :     ulong r1 = double_eta_power(inv, uel(rts, 1), p, pi);
    1560        5592 :     GEN r, f = Flx_double_eta_jpoly(F, r1, p, pi);
    1561        5592 :     if ((j = Flx_oneroot(f, p)) == p) pari_err_BUG("root_matrix");
    1562        5592 :     j = compute_L_isogenous_curve(L, n, ne, j, card, val, 0);
    1563        5592 :     avma = av;
    1564        5592 :     r1 = double_eta_power(inv, uel(surface_js, i), p, pi);
    1565        5592 :     f = Flx_double_eta_jpoly(F, r1, p, pi);
    1566        5592 :     r = Flx_roots(f, p);
    1567        5592 :     if (lg(r) != 3) pari_err_BUG("root_matrix");
    1568        5592 :     rev = (j != uel(r, 1)) && (j != uel(r, 2));
    1569        5592 :     avma = av0;
    1570             :   } else {
    1571             :     ulong j1pr, j1;
    1572       26694 :     j1pr = modfn_preimage(uel(rts, 1), ne, dinfo->inv);
    1573       26694 :     j1 = compute_L_isogenous_curve(L, n, ne, j1pr, card, val, 0);
    1574       26698 :     rev = j1 != modfn_preimage(uel(surface_js, i), ne, dinfo->inv);
    1575             :   }
    1576       32290 :   if (rev)
    1577       15910 :     carray_reverse_inplace(surface_js + 2, njs - 1);
    1578       32290 :   append_neighbours(rts, surface_js, njs, L, m, i);
    1579             : 
    1580      438035 :   for (i = 3; i <= njinvs; ++i) {
    1581      405746 :     update_Lsqr_cycle(cyc, dinfo);
    1582      405745 :     rts = gel(rt_mat, i);
    1583      405745 :     vecsmall_pick(rts, floor_js, cyc);
    1584      405764 :     append_neighbours(rts, surface_js, njs, L, m, i);
    1585             :   }
    1586       32289 :   avma = av; return rt_mat;
    1587             : }
    1588             : 
    1589             : INLINE void
    1590       34784 : interpolate_coeffs(GEN phi_modp, ulong p, GEN j_invs, GEN coeff_mat)
    1591             : {
    1592       34784 :   pari_sp av = avma;
    1593             :   long i;
    1594       34784 :   GEN pols = Flv_Flm_polint(j_invs, coeff_mat, p, 0);
    1595      762649 :   for (i = 1; i < lg(pols); ++i) {
    1596      727858 :     GEN pol = gel(pols, i);
    1597      727858 :     long k, maxk = lg(pol);
    1598      727858 :     for (k = 2; k < maxk; ++k) coeff(phi_modp, k - 1, i) = pol[k];
    1599             :   }
    1600       34791 :   avma = av;
    1601       34791 : }
    1602             : 
    1603             : INLINE long
    1604      341330 : Flv_lastnonzero(GEN v)
    1605             : {
    1606             :   long i;
    1607    24763480 :   for (i = lg(v) - 1; i > 0; --i)
    1608    24762721 :     if (v[i]) break;
    1609      341330 :   return i;
    1610             : }
    1611             : 
    1612             : /* Assuming the matrix of coefficients in phi corresponds to polynomials
    1613             :  * phi_k^* satisfying Y^c phi_k^*(Y^s) for c in {0, 1, ..., s} satisfying
    1614             :  * c + Lk = L + 1 (mod s), change phi so that the coefficients are for the
    1615             :  * polynomials Y^c phi_k^*(Y^s) (s is the sparsity factor) */
    1616             : INLINE void
    1617       10739 : inflate_polys(GEN phi, long L, long s)
    1618             : {
    1619       10739 :   long k, deg = L + 1;
    1620             :   long maxr;
    1621       10739 :   maxr = nbrows(phi);
    1622      362817 :   for (k = 0; k <= deg; ) {
    1623      341339 :     long i, c = smodss(L * (1 - k) + 1, s);
    1624             :     /* TODO: We actually know that the last non-zero element of gel(phi, k)
    1625             :      * can't be later than index n+1, where n is about (L + 1)/s. */
    1626      341324 :     ++k;
    1627     5530599 :     for (i = Flv_lastnonzero(gel(phi, k)); i > 0; --i) {
    1628     5189275 :       long r = c + (i - 1) * s + 1;
    1629     5189275 :       if (r > maxr) { coeff(phi, i, k) = 0; continue; }
    1630     5127505 :       if (r != i) {
    1631     5020228 :         coeff(phi, r, k) = coeff(phi, i, k);
    1632     5020228 :         coeff(phi, i, k) = 0;
    1633             :       }
    1634             :     }
    1635             :   }
    1636       10739 : }
    1637             : 
    1638             : INLINE void
    1639       39606 : Flv_powu_inplace_pre(GEN v, ulong n, ulong p, ulong pi)
    1640             : {
    1641             :   long i;
    1642       39606 :   for (i = 1; i < lg(v); ++i) v[i] = Fl_powu_pre(v[i], n, p, pi);
    1643       39609 : }
    1644             : 
    1645             : INLINE void
    1646       10739 : normalise_coeffs(GEN coeffs, GEN js, long L, long s, ulong p, ulong pi)
    1647             : {
    1648       10739 :   pari_sp av = avma;
    1649             :   long k;
    1650             :   GEN pows, modinv_js;
    1651             : 
    1652             :   /* NB: In fact it would be correct to return the coefficients "as is" when
    1653             :    * s = 1, but we make that an error anyway since this function should never
    1654             :    * be called with s = 1. */
    1655       10739 :   if (s <= 1) pari_err_BUG("normalise_coeffs");
    1656             : 
    1657             :   /* pows[i + 1] contains 1 / js[i + 1]^i for i = 0, ..., s - 1. */
    1658       10739 :   pows = cgetg(s + 1, t_VEC);
    1659       10739 :   gel(pows, 1) = const_vecsmall(lg(js) - 1, 1);
    1660       10739 :   modinv_js = Flv_inv_pre(js, p, pi);
    1661       10742 :   gel(pows, 2) = modinv_js;
    1662       37318 :   for (k = 3; k <= s; ++k) {
    1663       26577 :     gel(pows, k) = gcopy(modinv_js);
    1664       26577 :     Flv_powu_inplace_pre(gel(pows, k), k - 1, p, pi);
    1665             :   }
    1666             : 
    1667             :   /* For each column of coefficients coeffs[k] = [a0 .. an],
    1668             :    *   replace ai by ai / js[i]^c.
    1669             :    * Said in another way, normalise each row i of coeffs by
    1670             :    * dividing through by js[i - 1]^c (where c depends on i). */
    1671      352044 :   for (k = 1; k < lg(coeffs); ++k) {
    1672      341305 :     long i, c = smodss(L * (1 - (k - 1)) + 1, s);
    1673      341242 :     GEN col = gel(coeffs, k), C = gel(pows, c + 1);
    1674     5908525 :     for (i = 1; i < lg(col); ++i)
    1675     5567222 :       col[i] = Fl_mul_pre(col[i], C[i], p, pi);
    1676             :   }
    1677       10739 :   avma = av;
    1678       10739 : }
    1679             : 
    1680             : INLINE void
    1681        5592 : double_eta_initial_js(
    1682             :   ulong *x0, ulong *x0pr, ulong j0, ulong j0pr, norm_eqn_t ne,
    1683             :   long inv, ulong L, ulong n, ulong card, ulong val)
    1684             : {
    1685        5592 :   pari_sp av0 = avma;
    1686        5592 :   ulong p = ne->p, pi = ne->pi, s2 = ne->s2;
    1687        5592 :   GEN F = double_eta_Fl(inv, p);
    1688        5591 :   pari_sp av = avma;
    1689             :   ulong j1pr, j1, r, t;
    1690             :   GEN f, g;
    1691             : 
    1692        5591 :   *x0pr = modinv_double_eta_from_j(F, inv, j0pr, p, pi, s2);
    1693        5592 :   t = double_eta_power(inv, *x0pr, p, pi);
    1694        5592 :   f = Flx_div_by_X_x(Flx_double_eta_jpoly(F, t, p, pi), j0pr, p, &r);
    1695        5592 :   if (r) pari_err_BUG("double_eta_initial_js");
    1696        5592 :   j1pr = Flx_deg1_root(f, p);
    1697        5592 :   avma = av;
    1698             : 
    1699        5592 :   j1 = compute_L_isogenous_curve(L, n, ne, j1pr, card, val, 0);
    1700        5592 :   f = Flx_double_eta_xpoly(F, j0, p, pi);
    1701        5592 :   g = Flx_double_eta_xpoly(F, j1, p, pi);
    1702             :   /* x0 is the unique common root of f and g */
    1703        5592 :   *x0 = Flx_deg1_root(Flx_gcd(f, g, p), p);
    1704        5592 :   avma = av0;
    1705             : 
    1706        5592 :   if ( ! double_eta_root(inv, x0, *x0, p, pi, s2))
    1707           0 :     pari_err_BUG("double_eta_initial_js");
    1708        5591 : }
    1709             : 
    1710             : /* This is Sutherland 2012, Algorithm 2.1, p16. */
    1711             : static GEN
    1712       32278 : polmodular_split_p_Flm(ulong L, GEN hilb, GEN factu, norm_eqn_t ne, GEN db,
    1713             :   const disc_info *dinfo)
    1714             : {
    1715             :   ulong j0, j0_rt, j0pr, j0pr_rt;
    1716       32278 :   ulong n, card, val, p = ne->p, pi = ne->pi;
    1717       32278 :   long s = modinv_sparse_factor(dinfo->inv);
    1718       32281 :   long nj_selected = ceil((L + 1)/(double)s) + 1;
    1719             :   GEN surface_js, floor_js, rts;
    1720             :   GEN phi_modp;
    1721             :   GEN jdb, fdb;
    1722       32281 :   long switched_signs = 0;
    1723             : 
    1724       32281 :   jdb = polmodular_db_for_inv(db, INV_J);
    1725       32275 :   fdb = polmodular_db_for_inv(db, dinfo->inv);
    1726             : 
    1727             :   /* Precomputation */
    1728       32279 :   card = p + 1 - ne->t;
    1729       32279 :   val = u_lvalrem(card, L, &n); /* n = card / L^{v_L(card)} */
    1730             : 
    1731       32274 :   j0 = oneroot_of_classpoly(hilb, factu, ne, jdb);
    1732       32285 :   j0pr = compute_L_isogenous_curve(L, n, ne, j0, card, val, 1);
    1733       32290 :   if (modinv_is_double_eta(dinfo->inv)) {
    1734        5592 :     double_eta_initial_js(&j0_rt, &j0pr_rt, j0, j0pr, ne, dinfo->inv,
    1735             :         L, n, card, val);
    1736             :   } else {
    1737       26698 :     j0_rt = modfn_root(j0, ne, dinfo->inv);
    1738       26698 :     j0pr_rt = modfn_root(j0pr, ne, dinfo->inv);
    1739             :   }
    1740       32289 :   surface_js = enum_volcano_surface(dinfo, ne, j0_rt, fdb);
    1741       32290 :   floor_js = enum_volcano_floor(L, ne, j0pr_rt, fdb, dinfo);
    1742       32285 :   rts = root_matrix(L, dinfo, nj_selected, surface_js, floor_js,
    1743             :                     n, card, val, ne);
    1744        2166 :   do {
    1745       34455 :     pari_sp btop = avma;
    1746             :     long i;
    1747             :     GEN coeffs, surf;
    1748             : 
    1749       34455 :     coeffs = roots_to_coeffs(rts, p, L);
    1750       34446 :     surf = vecsmall_shorten(surface_js, nj_selected);
    1751       34450 :     if (s > 1) {
    1752       10739 :       normalise_coeffs(coeffs, surf, L, s, p, pi);
    1753       10739 :       Flv_powu_inplace_pre(surf, s, p, pi);
    1754             :     }
    1755       34450 :     phi_modp = zero_Flm_copy(L + 2, L + 2);
    1756       34452 :     interpolate_coeffs(phi_modp, p, surf, coeffs);
    1757       34450 :     if (s > 1) inflate_polys(phi_modp, L, s);
    1758             : 
    1759             :     /* TODO: Calculate just this coefficient of X^L Y^L, so we can do this
    1760             :      * test, then calculate the other coefficients; at the moment we are
    1761             :      * sometimes doing all the roots-to-coeffs, normalisation and interpolation
    1762             :      * work twice. */
    1763       34453 :     if (ucoeff(phi_modp, L + 1, L + 1) == p - 1) break;
    1764             : 
    1765        2166 :     if (switched_signs) pari_err_BUG("polmodular_split_p_Flm");
    1766             : 
    1767        2166 :     avma = btop;
    1768       27438 :     for (i = 1; i < lg(rts); ++i) {
    1769       25272 :       surface_js[i] = Fl_neg(surface_js[i], p);
    1770       25272 :       coeff(rts, L, i) = Fl_neg(coeff(rts, L, i), p);
    1771       25272 :       coeff(rts, L + 1, i) = Fl_neg(coeff(rts, L + 1, i), p);
    1772             :     }
    1773        2166 :     switched_signs = 1;
    1774             :   } while (1);
    1775       32287 :   dbg_printf(4)("  Phi_%lu(X, Y) (mod %lu) = %Ps\n", L, p, phi_modp);
    1776             : 
    1777       32287 :   return phi_modp;
    1778             : }
    1779             : 
    1780             : INLINE void
    1781       32261 : norm_eqn_update(norm_eqn_t ne, GEN vne, ulong t, ulong p, long L)
    1782             : {
    1783             :   long res;
    1784             :   ulong vL_sqr, vL;
    1785             : 
    1786       32261 :   ne->D = vne[1];
    1787       32261 :   ne->u = vne[2];
    1788       32261 :   ne->t = t;
    1789       32261 :   ne->p = p;
    1790       32261 :   ne->pi = get_Fl_red(p);
    1791       32267 :   ne->s2 = Fl_2gener_pre(p, ne->pi);
    1792             : 
    1793       32260 :   vL_sqr = (4 * p - t * t) / -ne->D;
    1794       32260 :   res = uissquareall(vL_sqr, &vL);
    1795       32280 :   if (!res || vL % L) pari_err_BUG("norm_eqn_update");
    1796       32280 :   ne->v = vL;
    1797             : 
    1798             :   /* select twisting parameter */
    1799       64411 :   do ne->T = random_Fl(p); while (krouu(ne->T, p) != -1);
    1800       32280 : }
    1801             : 
    1802             : INLINE void
    1803        2407 : Flv_deriv_pre_inplace(GEN v, long deg, ulong p, ulong pi)
    1804             : {
    1805        2407 :   long i, ln = lg(v), d = deg % p;
    1806        2407 :   for (i = ln - 1; i > 1; --i, --d) v[i] = Fl_mul_pre(v[i - 1], d, p, pi);
    1807        2408 :   v[1] = 0;
    1808        2408 : }
    1809             : 
    1810             : INLINE GEN
    1811        2408 : eval_modpoly_modp(GEN Tp, GEN j_powers, norm_eqn_t ne, int compute_derivs)
    1812             : {
    1813        2408 :   ulong p = ne->p, pi = ne->pi;
    1814        2408 :   long L = lg(j_powers) - 3;
    1815        2408 :   GEN j_pows_p = ZV_to_Flv(j_powers, p);
    1816        2408 :   GEN tmp = cgetg(2 + 2 * compute_derivs, t_VEC);
    1817             :   /* We wrap the result in this t_VEC Tp to trick the
    1818             :    * ZM_*_CRT() functions into thinking it's a matrix. */
    1819        2408 :   gel(tmp, 1) = Flm_Flc_mul_pre(Tp, j_pows_p, p, pi);
    1820        2408 :   if (compute_derivs) {
    1821        1204 :     Flv_deriv_pre_inplace(j_pows_p, L + 1, p, pi);
    1822        1204 :     gel(tmp, 2) = Flm_Flc_mul_pre(Tp, j_pows_p, p, pi);
    1823        1204 :     Flv_deriv_pre_inplace(j_pows_p, L + 1, p, pi);
    1824        1204 :     gel(tmp, 3) = Flm_Flc_mul_pre(Tp, j_pows_p, p, pi);
    1825             :   }
    1826        2408 :   return tmp;
    1827             : }
    1828             : 
    1829             : /* Parallel interface */
    1830             : GEN
    1831       32260 : polmodular_worker(ulong p, ulong t, ulong L,
    1832             :                   GEN hilb, GEN factu, GEN vne, GEN vinfo,
    1833             :                   long derivs, GEN j_powers, GEN fdb)
    1834             : {
    1835       32260 :   pari_sp av = avma;
    1836             :   norm_eqn_t ne;
    1837             :   GEN Tp;
    1838       32260 :   norm_eqn_update(ne, vne, t, p, L);
    1839       32280 :   Tp = polmodular_split_p_Flm(L, hilb, factu, ne, fdb, (const disc_info*)vinfo);
    1840       32284 :   if (!isintzero(j_powers)) Tp = eval_modpoly_modp(Tp, j_powers, ne, derivs);
    1841       32284 :   return gerepileupto(av, Tp);
    1842             : }
    1843             : 
    1844             : static GEN
    1845       18200 : sympol_to_ZM(GEN phi, long L)
    1846             : {
    1847       18200 :   pari_sp av = avma;
    1848       18200 :   GEN res = zeromatcopy(L + 2, L + 2);
    1849       18200 :   long i, j, c = 1;
    1850       79079 :   for (i = 1; i <= L + 1; ++i)
    1851      200662 :     for (j = 1; j <= i; ++j, ++c)
    1852      139783 :       gcoeff(res, i, j) = gcoeff(res, j, i) = gel(phi, c);
    1853       18200 :   gcoeff(res, L + 2, 1) = gcoeff(res, 1, L + 2) = gen_1;
    1854       18200 :   return gerepilecopy(av, res);
    1855             : }
    1856             : 
    1857             : static GEN polmodular_small_ZM(long L, long inv, GEN *db);
    1858             : 
    1859             : INLINE long
    1860       21011 : modinv_max_internal_level(long inv)
    1861             : {
    1862       21011 :   switch (inv) {
    1863       18446 :     case INV_J: return 5;
    1864         364 :     case INV_G2: return 2;
    1865             :     case INV_F:
    1866             :     case INV_F2:
    1867             :     case INV_F4:
    1868         443 :     case INV_F8: return 5;
    1869             :     case INV_W2W5:
    1870         252 :     case INV_W2W5E2: return 7;
    1871             :     case INV_W2W3:
    1872             :     case INV_W2W3E2:
    1873             :     case INV_W3W3:
    1874         441 :     case INV_W3W7:  return 5;
    1875          63 :     case INV_W3W3E2:return 2;
    1876             :     case INV_F3:
    1877             :     case INV_W2W7:
    1878             :     case INV_W2W7E2:
    1879         694 :     case INV_W2W13: return 3;
    1880             :     case INV_W3W5:
    1881             :     case INV_W5W7:
    1882         308 :     case INV_W3W13: return 2;
    1883             :   }
    1884             :   pari_err_BUG("modinv_max_internal_level"); return LONG_MAX;/*LCOV_EXCL_LINE*/
    1885             : }
    1886             : 
    1887             : GEN
    1888       20906 : polmodular0_ZM(long L, long inv, GEN J, GEN Q, int compute_derivs, GEN *db)
    1889             : {
    1890       20906 :   pari_sp ltop = avma;
    1891       20906 :   long k, d, Dcnt, nprimes = 0;
    1892             :   GEN modpoly, plist;
    1893             :   disc_info Ds[MODPOLY_MAX_DCNT];
    1894       20906 :   long lvl = modinv_level(inv);
    1895       20906 :   if (ugcd(L, lvl) != 1)
    1896           7 :     pari_err_DOMAIN("polmodular0_ZM", "invariant",
    1897             :                     "incompatible with", stoi(L), stoi(lvl));
    1898             : 
    1899       20899 :   dbg_printf(1)("Calculating modular polynomial of level %lu for invariant %d\n", L, inv);
    1900       20899 :   if (L <= modinv_max_internal_level(inv)) return polmodular_small_ZM(L,inv,db);
    1901             : 
    1902        2559 :   Dcnt = discriminant_with_classno_at_least(Ds, L, inv, USE_SPARSE_FACTOR);
    1903        2559 :   for (d = 0; d < Dcnt; d++) nprimes += Ds[d].nprimes;
    1904        2559 :   modpoly = cgetg(nprimes+1, t_VEC);
    1905        2559 :   plist = cgetg(nprimes+1, t_VECSMALL);
    1906        5136 :   for (d = 0, k = 1; d < Dcnt; d++)
    1907             :   {
    1908        2577 :     disc_info *dinfo = &Ds[d];
    1909             :     struct pari_mt pt;
    1910        2577 :     const long D = dinfo->D1, DK = dinfo->D0;
    1911        2577 :     const ulong cond = usqrt(D / DK);
    1912        2577 :     long i, pending = 0;
    1913             :     GEN worker, j_powers, factu, hilb;
    1914             : 
    1915        2577 :     polmodular_db_add_level(db, dinfo->L0, inv);
    1916        2577 :     if (dinfo->L1) polmodular_db_add_level(db, dinfo->L1, inv);
    1917        2577 :     factu = factoru(cond);
    1918        2577 :     dbg_printf(1)("Selected discriminant D = %ld = %ld^2 * %ld.\n", D,cond,DK);
    1919             : 
    1920        2577 :     hilb = polclass0(DK, INV_J, 0, db);
    1921        2577 :     if (cond > 1)
    1922          38 :       polmodular_db_add_levels(db, zv_to_longptr(gel(factu,1)), lg(gel(factu,1))-1, INV_J);
    1923             : 
    1924        2577 :     dbg_printf(1)("D = %ld, L0 = %lu, L1 = %lu, ", dinfo->D1, dinfo->L0, dinfo->L1);
    1925        2577 :     dbg_printf(1)("n1 = %lu, n2 = %lu, dl1 = %lu, dl2_0 = %lu, dl2_1 = %lu\n",
    1926             :           dinfo->n1, dinfo->n2, dinfo->dl1, dinfo->dl2_0, dinfo->dl2_1);
    1927        2577 :     dbg_printf(0)("Calculating modular polynomial of level %lu:", L);
    1928             : 
    1929        2577 :     j_powers = gen_0;
    1930        2577 :     if (J) {
    1931          56 :       compute_derivs = !!compute_derivs;
    1932          56 :       j_powers = Fp_powers(J, L+1, Q);
    1933             :     }
    1934        2577 :     worker = strtoclosure("_polmodular_worker", 8, utoi(L), hilb, factu,
    1935             :         mkvecsmall2(D, cond), (GEN)dinfo, stoi(compute_derivs), j_powers, *db);
    1936        2577 :     mt_queue_start_lim(&pt, worker, dinfo->nprimes);
    1937       38840 :     for (i = 0; i < dinfo->nprimes || pending; i++)
    1938             :     {
    1939             :       GEN done;
    1940             :       long workid;
    1941       36263 :       ulong p = dinfo->primes[i], t = dinfo->traces[i];
    1942       36263 :       mt_queue_submit(&pt, p, i < dinfo->nprimes? mkvec2(utoi(p), utoi(t)): NULL);
    1943       36263 :       done = mt_queue_get(&pt, &workid, &pending);
    1944       36263 :       if (done)
    1945             :       {
    1946       32290 :         gel(modpoly, k) = done;
    1947       32290 :         plist[k] = workid; k++;
    1948       32290 :         dbg_printf(0)(" %ld%%", k*100/nprimes);
    1949             :       }
    1950             :     }
    1951        2577 :     dbg_printf(0)("\n");
    1952        2577 :     mt_queue_end(&pt);
    1953        2577 :     killblock((GEN)dinfo->primes);
    1954             :   }
    1955        2559 :   modpoly = nmV_chinese_center(modpoly, plist, NULL);
    1956        2559 :   if (J) modpoly = FpM_red(modpoly, Q);
    1957        2559 :   return gerepileupto(ltop, modpoly);
    1958             : }
    1959             : 
    1960             : GEN
    1961       14574 : polmodular_ZM(long L, long inv)
    1962             : {
    1963             :   GEN db, Phi;
    1964             : 
    1965       14574 :   if (L < 2)
    1966           7 :     pari_err_DOMAIN("polmodular_ZM", "L", "<", gen_2, stoi(L));
    1967             : 
    1968             :   /* TODO: Handle non-prime L.  This is Algorithm 1.1 and Corollary
    1969             :    * 3.4 in Sutherland, "Class polynomials for nonholomorphic modular
    1970             :    * functions". */
    1971       14567 :   if ( ! uisprime(L))
    1972           7 :     pari_err_IMPL("composite level");
    1973             : 
    1974       14560 :   db = polmodular_db_init(inv);
    1975       14560 :   Phi = polmodular0_ZM(L, inv, NULL, NULL, 0, &db);
    1976       14553 :   gunclone_deep(db);
    1977             : 
    1978       14553 :   return Phi;
    1979             : }
    1980             : 
    1981             : GEN
    1982       14497 : polmodular_ZXX(long L, long inv, long vx, long vy)
    1983             : {
    1984       14497 :   pari_sp av = avma;
    1985       14497 :   GEN phi = polmodular_ZM(L, inv);
    1986             : 
    1987       14476 :   if (vx < 0) vx = 0;
    1988       14476 :   if (vy < 0) vy = 1;
    1989       14476 :   if (varncmp(vx, vy) >= 0)
    1990          14 :     pari_err_PRIORITY("polmodular_ZXX", pol_x(vx), "<=", vy);
    1991       14462 :   return gerepilecopy(av, RgM_to_RgXX(phi, vx, vy));
    1992             : }
    1993             : 
    1994             : INLINE GEN
    1995          56 : FpV_deriv(GEN v, long deg, GEN P)
    1996             : {
    1997          56 :   long i, ln = lg(v);
    1998          56 :   GEN dv = cgetg(ln, t_VEC);
    1999         392 :   for (i = ln - 1; i > 1; --i, --deg)
    2000         336 :     gel(dv, i) = Fp_mulu(gel(v, i - 1), deg, P);
    2001          56 :   gel(dv, 1) = gen_0;
    2002          56 :   return dv;
    2003             : }
    2004             : 
    2005             : GEN
    2006         112 : Fp_polmodular_evalx(
    2007             :   long L, long inv, GEN J, GEN P, long v, int compute_derivs)
    2008             : {
    2009         112 :   pari_sp av = avma;
    2010             :   GEN db, phi;
    2011             : 
    2012         112 :   if (L <= modinv_max_internal_level(inv)) {
    2013             :     GEN tmp;
    2014          56 :     GEN phi = RgM_to_FpM(polmodular_ZM(L, inv), P);
    2015          56 :     GEN j_powers = Fp_powers(J, L + 1, P);
    2016          56 :     GEN modpol = RgV_to_RgX(FpM_FpC_mul(phi, j_powers, P), v);
    2017          56 :     if (compute_derivs) {
    2018          28 :       tmp = cgetg(4, t_VEC);
    2019          28 :       gel(tmp, 1) = modpol;
    2020          28 :       j_powers = FpV_deriv(j_powers, L + 1, P);
    2021          28 :       gel(tmp, 2) = RgV_to_RgX(FpM_FpC_mul(phi, j_powers, P), v);
    2022          28 :       j_powers = FpV_deriv(j_powers, L + 1, P);
    2023          28 :       gel(tmp, 3) = RgV_to_RgX(FpM_FpC_mul(phi, j_powers, P), v);
    2024             :     } else
    2025          28 :       tmp = modpol;
    2026          56 :     return gerepilecopy(av, tmp);
    2027             :   }
    2028             : 
    2029          56 :   db = polmodular_db_init(inv);
    2030          56 :   phi = polmodular0_ZM(L, inv, J, P, compute_derivs, &db);
    2031          56 :   phi = RgM_to_RgXV(phi, v);
    2032          56 :   gunclone_deep(db);
    2033          56 :   return gerepilecopy(av, compute_derivs ? phi : gel(phi, 1));
    2034             : }
    2035             : 
    2036             : GEN
    2037         609 : polmodular(long L, long inv, GEN x, long v, long compute_derivs)
    2038             : {
    2039         609 :   pari_sp av = avma;
    2040             :   long tx;
    2041         609 :   GEN J = NULL, P = NULL, res = NULL, one = NULL;
    2042             : 
    2043         609 :   check_modinv(inv);
    2044         602 :   if ( ! x || gequalX(x)) {
    2045         476 :     long xv = 0;
    2046         476 :     if (x) xv = varn(x);
    2047         476 :     if (compute_derivs) pari_err_FLAG("polmodular");
    2048         469 :     return polmodular_ZXX(L, inv, xv, v);
    2049             :   }
    2050             : 
    2051         126 :   tx = typ(x);
    2052         126 :   if (tx == t_INTMOD) {
    2053          56 :     J = gel(x, 2);
    2054          56 :     P = gel(x, 1);
    2055          56 :     one = mkintmod(gen_1, P);
    2056          70 :   } else if (tx == t_FFELT) {
    2057          63 :     J = FF_to_FpXQ_i(x);
    2058          63 :     if (degpol(J) > 0)
    2059           7 :       pari_err_DOMAIN("polmodular", "x", "not in prime subfield ", gen_0, x);
    2060          56 :     J = constant_coeff(J);
    2061          56 :     P = FF_p_i(x);
    2062          56 :     one = p_to_FF(P, 0);
    2063             :   } else {
    2064           7 :     pari_err_TYPE("polmodular", x);
    2065             :   }
    2066             : 
    2067         112 :   if (v < 0) v = 1;
    2068         112 :   res = Fp_polmodular_evalx(L, inv, J, P, v, compute_derivs);
    2069         112 :   res = gmul(res, one);
    2070         112 :   return gerepileupto(av, res);
    2071             : }
    2072             : 
    2073             : /**
    2074             :  * SECTION: Modular polynomials of level <= MAX_INTERNAL_MODPOLY_LEVEL.
    2075             :  */
    2076             : 
    2077             : /* These functions return a vector of unique coefficients of classical
    2078             :  * modular polynomials \Phi_L(X, Y) of small level L.  The number of
    2079             :  * such coefficients is (L + 1)(L + 2)/2 since \Phi is symmetric.
    2080             :  * (Note that we omit the common coefficient of X^{L + 1} and Y^{L + 1} since
    2081             :  * it is always 1.)  See sympol_to_ZM() for how to interpret the coefficients,
    2082             :  * and use that function to get the corresponding full (desymmetrised) matrix
    2083             :  * of coefficients. */
    2084             : 
    2085             : /*  Phi2, the modular polynomial of level 2:
    2086             :  *
    2087             :  *  X^3
    2088             :  *  + X^2 * (-Y^2 + 1488*Y - 162000)
    2089             :  *  + X * (1488*Y^2 + 40773375*Y + 8748000000)
    2090             :  *  + Y^3 - 162000*Y^2 + 8748000000*Y - 157464000000000
    2091             :  *
    2092             :  *  [[3, 0, 1],
    2093             :  *   [2, 2, -1],
    2094             :  *   [2, 1, 1488],
    2095             :  *   [2, 0, -162000],
    2096             :  *   [1, 1, 40773375],
    2097             :  *   [1, 0, 8748000000],
    2098             :  *   [0, 0, -157464000000000]], */
    2099             : static GEN
    2100       15022 : phi2_ZV(void)
    2101             : {
    2102       15022 :   GEN phi2 = cgetg(7, t_VEC);
    2103       15022 :   gel(phi2, 1) = uu32toi(36662, 1908994048);
    2104       15022 :   setsigne(gel(phi2, 1), -1);
    2105       15022 :   gel(phi2, 2) = uu32toi(2, 158065408);
    2106       15022 :   gel(phi2, 3) = stoi(40773375);
    2107       15022 :   gel(phi2, 4) = stoi(-162000);
    2108       15022 :   gel(phi2, 5) = stoi(1488);
    2109       15022 :   gel(phi2, 6) = gen_m1;
    2110       15022 :   return phi2;
    2111             : }
    2112             : 
    2113             : /* L = 3
    2114             :  *
    2115             :  * [4, 0, 1],
    2116             :  * [3, 3, -1],
    2117             :  * [3, 2, 2232],
    2118             :  * [3, 1, -1069956],
    2119             :  * [3, 0, 36864000],
    2120             :  * [2, 2, 2587918086],
    2121             :  * [2, 1, 8900222976000],
    2122             :  * [2, 0, 452984832000000],
    2123             :  * [1, 1, -770845966336000000],
    2124             :  * [1, 0, 1855425871872000000000]
    2125             :  * [0, 0, 0]
    2126             :  *
    2127             :  * X^4
    2128             :  * + X^3 (-Y^3 + 2232*Y^2 - 1069956*Y + 36864000)
    2129             :  * + X^2 (2232*Y^3 + 2587918086*Y^2 + 8900222976000*Y + 452984832000000)
    2130             :  * + X (-1069956*Y^3 + 8900222976000*Y^2 - 770845966336000000*Y + 1855425871872000000000)
    2131             :  * + Y^4 + 36864000*Y^3 + 452984832000000*Y^2 + 1855425871872000000000*Y
    2132             :  *
    2133             :  * 1855425871872000000000 == 2^32 * (100 * 2^32 + 2503270400) */
    2134             : static GEN
    2135        1134 : phi3_ZV(void)
    2136             : {
    2137        1134 :   GEN phi3 = cgetg(11, t_VEC);
    2138        1134 :   pari_sp av = avma;
    2139        1134 :   gel(phi3, 1) = gen_0;
    2140        1134 :   gel(phi3, 2) = gerepileupto(av, shifti(uu32toi(100, 2503270400UL), 32));
    2141        1134 :   gel(phi3, 3) = uu32toi(179476562, 2147483648UL);
    2142        1134 :   setsigne(gel(phi3, 3), -1);
    2143        1134 :   gel(phi3, 4) = uu32toi(105468, 3221225472UL);
    2144        1134 :   gel(phi3, 5) = uu32toi(2072, 1050738688);
    2145        1134 :   gel(phi3, 6) = utoi(2587918086UL);
    2146        1134 :   gel(phi3, 7) = stoi(36864000);
    2147        1134 :   gel(phi3, 8) = stoi(-1069956);
    2148        1134 :   gel(phi3, 9) = stoi(2232);
    2149        1134 :   gel(phi3, 10) = gen_m1;
    2150        1134 :   return phi3;
    2151             : }
    2152             : 
    2153             : static GEN
    2154        1113 : phi5_ZV(void)
    2155             : {
    2156        1113 :   GEN phi5 = cgetg(22, t_VEC);
    2157        1113 :   gel(phi5, 1) = mkintn(5, 0x18c2cc9cUL, 0x484382b2UL, 0xdc000000UL, 0x0UL, 0x0UL);
    2158        1113 :   gel(phi5, 2) = mkintn(5, 0x2638fUL, 0x2ff02690UL, 0x68026000UL, 0x0UL, 0x0UL);
    2159        1113 :   gel(phi5, 3) = mkintn(5, 0x308UL, 0xac9d9a4UL, 0xe0fdab12UL, 0xc0000000UL, 0x0UL);
    2160        1113 :   setsigne(gel(phi5, 3), -1);
    2161        1113 :   gel(phi5, 4) = mkintn(5, 0x13UL, 0xaae09f9dUL, 0x1b5ef872UL, 0x30000000UL, 0x0UL);
    2162        1113 :   gel(phi5, 5) = mkintn(4, 0x1b802fa9UL, 0x77ba0653UL, 0xd2f78000UL, 0x0UL);
    2163        1113 :   gel(phi5, 6) = mkintn(4, 0xfbfdUL, 0x278e4756UL, 0xdf08a7c4UL, 0x40000000UL);
    2164        1113 :   gel(phi5, 7) = mkintn(4, 0x35f922UL, 0x62ccea6fUL, 0x153d0000UL, 0x0UL);
    2165        1113 :   gel(phi5, 8) = mkintn(4, 0x97dUL, 0x29203fafUL, 0xc3036909UL, 0x80000000UL);
    2166        1113 :   setsigne(gel(phi5, 8), -1);
    2167        1113 :   gel(phi5, 9) = mkintn(3, 0x56e9e892UL, 0xd7781867UL, 0xf2ea0000UL);
    2168        1113 :   gel(phi5, 10) = mkintn(3, 0x5d6dUL, 0xe0a58f4eUL, 0x9ee68c14UL);
    2169        1113 :   setsigne(gel(phi5, 10), -1);
    2170        1113 :   gel(phi5, 11) = mkintn(3, 0x1100dUL, 0x85cea769UL, 0x40000000UL);
    2171        1113 :   gel(phi5, 12) = mkintn(3, 0x1b38UL, 0x43cf461fUL, 0x3a900000UL);
    2172        1113 :   gel(phi5, 13) = mkintn(3, 0x14UL, 0xc45a616eUL, 0x4801680fUL);
    2173        1113 :   gel(phi5, 14) = uu32toi(0x17f4350UL, 0x493ca3e0UL);
    2174        1113 :   gel(phi5, 15) = uu32toi(0x183UL, 0xe54ce1f8UL);
    2175        1113 :   gel(phi5, 16) = uu32toi(0x1c9UL, 0x18860000UL);
    2176        1113 :   gel(phi5, 17) = uu32toi(0x39UL, 0x6f7a2206UL);
    2177        1113 :   setsigne(gel(phi5, 17), -1);
    2178        1113 :   gel(phi5, 18) = stoi(2028551200);
    2179        1113 :   gel(phi5, 19) = stoi(-4550940);
    2180        1113 :   gel(phi5, 20) = stoi(3720);
    2181        1113 :   gel(phi5, 21) = gen_m1;
    2182        1113 :   return phi5;
    2183             : }
    2184             : 
    2185             : static GEN
    2186         182 : phi5_f_ZV(void)
    2187             : {
    2188         182 :   GEN phi = zerovec(21);
    2189         182 :   gel(phi, 3) = stoi(4);
    2190         182 :   gel(phi, 21) = gen_m1;
    2191         182 :   return phi;
    2192             : }
    2193             : 
    2194             : static GEN
    2195          21 : phi3_f3_ZV(void)
    2196             : {
    2197          21 :   GEN phi = zerovec(10);
    2198          21 :   gel(phi, 3) = stoi(8);
    2199          21 :   gel(phi, 10) = gen_m1;
    2200          21 :   return phi;
    2201             : }
    2202             : 
    2203             : static GEN
    2204         119 : phi2_g2_ZV(void)
    2205             : {
    2206         119 :   GEN phi = zerovec(6);
    2207         119 :   gel(phi, 1) = stoi(-54000);
    2208         119 :   gel(phi, 3) = stoi(495);
    2209         119 :   gel(phi, 6) = gen_m1;
    2210         119 :   return phi;
    2211             : }
    2212             : 
    2213             : static GEN
    2214          56 : phi5_w2w3_ZV(void)
    2215             : {
    2216          56 :   GEN phi = zerovec(21);
    2217          56 :   gel(phi, 3) = gen_m1;
    2218          56 :   gel(phi, 10) = stoi(5);
    2219          56 :   gel(phi, 21) = gen_m1;
    2220          56 :   return phi;
    2221             : }
    2222             : 
    2223             : static GEN
    2224         112 : phi7_w2w5_ZV(void)
    2225             : {
    2226         112 :   GEN phi = zerovec(36);
    2227         112 :   gel(phi, 3) = gen_m1;
    2228         112 :   gel(phi, 15) = stoi(56);
    2229         112 :   gel(phi, 19) = stoi(42);
    2230         112 :   gel(phi, 24) = stoi(21);
    2231         112 :   gel(phi, 30) = stoi(7);
    2232         112 :   gel(phi, 36) = gen_m1;
    2233         112 :   return phi;
    2234             : }
    2235             : 
    2236             : static GEN
    2237          56 : phi5_w3w3_ZV(void)
    2238             : {
    2239          56 :   GEN phi = zerovec(21);
    2240          56 :   gel(phi, 3) = stoi(9);
    2241          56 :   gel(phi, 6) = stoi(-15);
    2242          56 :   gel(phi, 15) = stoi(5);
    2243          56 :   gel(phi, 21) = gen_m1;
    2244          56 :   return phi;
    2245             : }
    2246             : 
    2247             : static GEN
    2248         168 : phi3_w2w7_ZV(void)
    2249             : {
    2250         168 :   GEN phi = zerovec(10);
    2251         168 :   gel(phi, 3) = gen_m1;
    2252         168 :   gel(phi, 6) = stoi(3);
    2253         168 :   gel(phi, 10) = gen_m1;
    2254         168 :   return phi;
    2255             : }
    2256             : 
    2257             : static GEN
    2258          35 : phi2_w3w5_ZV(void)
    2259             : {
    2260          35 :   GEN phi = zerovec(6);
    2261          35 :   gel(phi, 3) = gen_1;
    2262          35 :   gel(phi, 6) = gen_m1;
    2263          35 :   return phi;
    2264             : }
    2265             : 
    2266             : static GEN
    2267          42 : phi5_w3w7_ZV(void)
    2268             : {
    2269          42 :   GEN phi = zerovec(21);
    2270          42 :   gel(phi, 3) = gen_m1;
    2271          42 :   gel(phi, 6) = stoi(10);
    2272          42 :   gel(phi, 8) = stoi(5);
    2273          42 :   gel(phi, 10) = stoi(35);
    2274          42 :   gel(phi, 13) = stoi(20);
    2275          42 :   gel(phi, 15) = stoi(10);
    2276          42 :   gel(phi, 17) = stoi(5);
    2277          42 :   gel(phi, 19) = stoi(5);
    2278          42 :   gel(phi, 21) = gen_m1;
    2279          42 :   return phi;
    2280             : }
    2281             : 
    2282             : static GEN
    2283          49 : phi3_w2w13_ZV(void)
    2284             : {
    2285          49 :   GEN phi = zerovec(10);
    2286          49 :   gel(phi, 3) = gen_m1;
    2287          49 :   gel(phi, 6) = stoi(3);
    2288          49 :   gel(phi, 8) = stoi(3);
    2289          49 :   gel(phi, 10) = gen_m1;
    2290          49 :   return phi;
    2291             : }
    2292             : 
    2293             : static GEN
    2294          21 : phi2_w3w3e2_ZV(void)
    2295             : {
    2296          21 :   GEN phi = zerovec(6);
    2297          21 :   gel(phi, 3) = stoi(3);
    2298          21 :   gel(phi, 6) = gen_m1;
    2299          21 :   return phi;
    2300             : }
    2301             : 
    2302             : static GEN
    2303          56 : phi2_w5w7_ZV(void)
    2304             : {
    2305          56 :   GEN phi = zerovec(6);
    2306          56 :   gel(phi, 3) = gen_1;
    2307          56 :   gel(phi, 5) = gen_2;
    2308          56 :   gel(phi, 6) = gen_m1;
    2309          56 :   return phi;
    2310             : }
    2311             : 
    2312             : static GEN
    2313          14 : phi2_w3w13_ZV(void)
    2314             : {
    2315          14 :   GEN phi = zerovec(6);
    2316          14 :   gel(phi, 3) = gen_m1;
    2317          14 :   gel(phi, 5) = gen_2;
    2318          14 :   gel(phi, 6) = gen_m1;
    2319          14 :   return phi;
    2320             : }
    2321             : 
    2322             : INLINE long
    2323         140 : modinv_parent(long inv)
    2324             : {
    2325         140 :   switch (inv) {
    2326             :     case INV_F2:
    2327             :     case INV_F4:
    2328          42 :     case INV_F8:     return INV_F;
    2329          14 :     case INV_W2W3E2: return INV_W2W3;
    2330          28 :     case INV_W2W5E2: return INV_W2W5;
    2331          56 :     case INV_W2W7E2: return INV_W2W7;
    2332           0 :     case INV_W3W3E2: return INV_W3W3;
    2333             :     default: pari_err_BUG("modinv_parent"); return -1;/*LCOV_EXCL_LINE*/
    2334             :   }
    2335             : }
    2336             : 
    2337             : /* TODO: Think of a better name than "parent power"; sheesh. */
    2338             : INLINE long
    2339         140 : modinv_parent_power(long inv)
    2340             : {
    2341         140 :   switch (inv) {
    2342          14 :     case INV_F4: return 4;
    2343          14 :     case INV_F8: return 8;
    2344             :     case INV_F2:
    2345             :     case INV_W2W3E2:
    2346             :     case INV_W2W5E2:
    2347             :     case INV_W2W7E2:
    2348         112 :     case INV_W3W3E2: return 2;
    2349             :     default: pari_err_BUG("modinv_parent_power"); return -1;/*LCOV_EXCL_LINE*/
    2350             :   }
    2351             : }
    2352             : 
    2353             : static GEN
    2354         140 : polmodular0_powerup_ZM(long L, long inv, GEN *db)
    2355             : {
    2356         140 :   pari_sp ltop = avma, av;
    2357             :   long s, D, nprimes, N;
    2358             :   GEN mp, pol, P, H;
    2359         140 :   long parent = modinv_parent(inv);
    2360         140 :   long e = modinv_parent_power(inv);
    2361             :   disc_info Ds[MODPOLY_MAX_DCNT];
    2362             :   /* FIXME: We throw away the table of fundamental discriminants here. */
    2363         140 :   long nDs = discriminant_with_classno_at_least(Ds, L, inv, IGNORE_SPARSE_FACTOR);
    2364         140 :   if (nDs != 1) pari_err_BUG("polmodular0_powerup_ZM");
    2365         140 :   D = Ds[0].D1;
    2366         140 :   nprimes = Ds[0].nprimes + 1;
    2367         140 :   mp = polmodular0_ZM(L, parent, NULL, NULL, 0, db);
    2368         140 :   H = polclass0(D, parent, 0, db);
    2369             : 
    2370         140 :   N = L + 2;
    2371         140 :   if (degpol(H) < N) pari_err_BUG("polmodular0_powerup_ZM");
    2372             : 
    2373         140 :   av = avma;
    2374         140 :   pol = ZM_init_CRT(zero_Flm_copy(N, L + 2), 1);
    2375         140 :   P = gen_1;
    2376         476 :   for (s = 1; s < nprimes; ++s) {
    2377             :     pari_sp av1, av2;
    2378         336 :     ulong p = Ds[0].primes[s-1], pi = get_Fl_red(p);
    2379             :     long i;
    2380             :     GEN Hrts, js, Hp, Phip, coeff_mat, phi_modp;
    2381             : 
    2382         336 :     phi_modp = zero_Flm_copy(N, L + 2);
    2383         336 :     av1 = avma;
    2384         336 :     Hp = ZX_to_Flx(H, p);
    2385         336 :     Hrts = Flx_roots(Hp, p);
    2386         336 :     if (lg(Hrts)-1 < N) pari_err_BUG("polmodular0_powerup_ZM");
    2387         336 :     js = cgetg(N + 1, t_VECSMALL);
    2388        2632 :     for (i = 1; i <= N; ++i)
    2389        2296 :       uel(js, i) = Fl_powu_pre(uel(Hrts, i), e, p, pi);
    2390             : 
    2391         336 :     Phip = ZM_to_Flm(mp, p);
    2392         336 :     coeff_mat = zero_Flm_copy(N, L + 2);
    2393         336 :     av2 = avma;
    2394        2632 :     for (i = 1; i <= N; ++i) {
    2395             :       long k;
    2396             :       GEN phi_at_ji, mprts;
    2397             : 
    2398        2296 :       phi_at_ji = Flm_Fl_polmodular_evalx(Phip, L, uel(Hrts, i), p, pi);
    2399        2296 :       mprts = Flx_roots(phi_at_ji, p);
    2400        2296 :       if (lg(mprts) != L + 2) pari_err_BUG("polmodular0_powerup_ZM");
    2401             : 
    2402        2296 :       Flv_powu_inplace_pre(mprts, e, p, pi);
    2403        2296 :       phi_at_ji = Flv_roots_to_pol(mprts, p, 0);
    2404             : 
    2405       18760 :       for (k = 1; k <= L + 2; ++k)
    2406       16464 :         ucoeff(coeff_mat, i, k) = uel(phi_at_ji, k + 1);
    2407        2296 :       avma = av2;
    2408             :     }
    2409             : 
    2410         336 :     interpolate_coeffs(phi_modp, p, js, coeff_mat);
    2411         336 :     avma = av1;
    2412             : 
    2413         336 :     (void) ZM_incremental_CRT(&pol, phi_modp, &P, p);
    2414         336 :     if (gc_needed(av, 2)) gerepileall(av, 2, &pol, &P);
    2415             :   }
    2416         140 :   killblock((GEN)Ds[0].primes); return gerepileupto(ltop, pol);
    2417             : }
    2418             : 
    2419             : /* Returns the modular polynomial with the smallest level for the given
    2420             :  * invariant, except if inv is INV_J, in which case return the modular
    2421             :  * polynomial of level L in {2,3,5}.  NULL is returned if the modular
    2422             :  * polynomial can be calculated using polmodular0_powerup_ZM. */
    2423             : INLINE GEN
    2424       18340 : internal_db(long L, long inv)
    2425             : {
    2426       18340 :   switch (inv) {
    2427       17269 :   case INV_J: switch (L) {
    2428       15022 :     case 2: return phi2_ZV();
    2429        1134 :     case 3: return phi3_ZV();
    2430        1113 :     case 5: return phi5_ZV();
    2431           0 :     default: break;
    2432             :   }
    2433         182 :   case INV_F: return phi5_f_ZV();
    2434          14 :   case INV_F2: return NULL;
    2435          21 :   case INV_F3: return phi3_f3_ZV();
    2436          14 :   case INV_F4: return NULL;
    2437         119 :   case INV_G2: return phi2_g2_ZV();
    2438          56 :   case INV_W2W3: return phi5_w2w3_ZV();
    2439          14 :   case INV_F8: return NULL;
    2440          56 :   case INV_W3W3: return phi5_w3w3_ZV();
    2441         112 :   case INV_W2W5: return phi7_w2w5_ZV();
    2442         168 :   case INV_W2W7: return phi3_w2w7_ZV();
    2443          35 :   case INV_W3W5: return phi2_w3w5_ZV();
    2444          42 :   case INV_W3W7: return phi5_w3w7_ZV();
    2445          14 :   case INV_W2W3E2: return NULL;
    2446          28 :   case INV_W2W5E2: return NULL;
    2447          49 :   case INV_W2W13: return phi3_w2w13_ZV();
    2448          56 :   case INV_W2W7E2: return NULL;
    2449          21 :   case INV_W3W3E2: return phi2_w3w3e2_ZV();
    2450          56 :   case INV_W5W7: return phi2_w5w7_ZV();
    2451          14 :   case INV_W3W13: return phi2_w3w13_ZV();
    2452             :   }
    2453           0 :   pari_err_BUG("internal_db");
    2454           0 :   return NULL;
    2455             : }
    2456             : 
    2457             : /* NB: Should only be called if L <= modinv_max_internal_level(inv) */
    2458             : static GEN
    2459       18340 : polmodular_small_ZM(long L, long inv, GEN *db)
    2460             : {
    2461       18340 :   GEN f = internal_db(L, inv);
    2462       18340 :   if (!f) return polmodular0_powerup_ZM(L, inv, db);
    2463       18200 :   return sympol_to_ZM(f, L);
    2464             : }
    2465             : 
    2466             : /* Each function phi_w?w?_j() returns a vector V containing two
    2467             :  * vectors u and v, and a scalar k, which together represent the
    2468             :  * bivariate polnomial
    2469             :  *
    2470             :  *   phi(X, Y) = \sum_i u[i] X^i + Y \sum_i v[i] X^i + Y^2 X^k
    2471             :  */
    2472             : static GEN
    2473        1078 : phi_w2w3_j(void)
    2474             : {
    2475             :   GEN phi, phi0, phi1;
    2476        1078 :   phi = cgetg(4, t_VEC);
    2477             : 
    2478        1078 :   phi0 = cgetg(14, t_VEC);
    2479        1078 :   gel(phi0, 1) = gen_1;
    2480        1078 :   gel(phi0, 2) = utoineg(0x3cUL);
    2481        1078 :   gel(phi0, 3) = utoi(0x702UL);
    2482        1078 :   gel(phi0, 4) = utoineg(0x797cUL);
    2483        1078 :   gel(phi0, 5) = utoi(0x5046fUL);
    2484        1078 :   gel(phi0, 6) = utoineg(0x1be0b8UL);
    2485        1078 :   gel(phi0, 7) = utoi(0x28ef9cUL);
    2486        1078 :   gel(phi0, 8) = utoi(0x15e2968UL);
    2487        1078 :   gel(phi0, 9) = utoi(0x1b8136fUL);
    2488        1078 :   gel(phi0, 10) = utoi(0xa67674UL);
    2489        1078 :   gel(phi0, 11) = utoi(0x23982UL);
    2490        1078 :   gel(phi0, 12) = utoi(0x294UL);
    2491        1078 :   gel(phi0, 13) = gen_1;
    2492             : 
    2493        1078 :   phi1 = cgetg(13, t_VEC);
    2494        1078 :   gel(phi1, 1) = gen_0;
    2495        1078 :   gel(phi1, 2) = gen_0;
    2496        1078 :   gel(phi1, 3) = gen_m1;
    2497        1078 :   gel(phi1, 4) = utoi(0x23UL);
    2498        1078 :   gel(phi1, 5) = utoineg(0xaeUL);
    2499        1078 :   gel(phi1, 6) = utoineg(0x5b8UL);
    2500        1078 :   gel(phi1, 7) = utoi(0x12d7UL);
    2501        1078 :   gel(phi1, 8) = utoineg(0x7c86UL);
    2502        1078 :   gel(phi1, 9) = utoi(0x37c8UL);
    2503        1078 :   gel(phi1, 10) = utoineg(0x69cUL);
    2504        1078 :   gel(phi1, 11) = utoi(0x48UL);
    2505        1078 :   gel(phi1, 12) = gen_m1;
    2506             : 
    2507        1078 :   gel(phi, 1) = phi0;
    2508        1078 :   gel(phi, 2) = phi1;
    2509        1078 :   gel(phi, 3) = utoi(5); return phi;
    2510             : }
    2511             : 
    2512             : static GEN
    2513        3219 : phi_w3w3_j(void)
    2514             : {
    2515             :   GEN phi, phi0, phi1;
    2516        3219 :   phi = cgetg(4, t_VEC);
    2517             : 
    2518        3219 :   phi0 = cgetg(14, t_VEC);
    2519        3219 :   gel(phi0, 1) = utoi(0x2d9UL);
    2520        3219 :   gel(phi0, 2) = utoi(0x4fbcUL);
    2521        3219 :   gel(phi0, 3) = utoi(0x5828aUL);
    2522        3219 :   gel(phi0, 4) = utoi(0x3a7a3cUL);
    2523        3219 :   gel(phi0, 5) = utoi(0x1bd8edfUL);
    2524        3219 :   gel(phi0, 6) = utoi(0x8348838UL);
    2525        3219 :   gel(phi0, 7) = utoi(0x1983f8acUL);
    2526        3219 :   gel(phi0, 8) = utoi(0x14e4e098UL);
    2527        3219 :   gel(phi0, 9) = utoi(0x69ed1a7UL);
    2528        3219 :   gel(phi0, 10) = utoi(0xc3828cUL);
    2529        3219 :   gel(phi0, 11) = utoi(0x2696aUL);
    2530        3219 :   gel(phi0, 12) = utoi(0x2acUL);
    2531        3219 :   gel(phi0, 13) = gen_1;
    2532             : 
    2533        3219 :   phi1 = cgetg(13, t_VEC);
    2534        3219 :   gel(phi1, 1) = gen_0;
    2535        3219 :   gel(phi1, 2) = utoineg(0x1bUL);
    2536        3219 :   gel(phi1, 3) = utoineg(0x5d6UL);
    2537        3219 :   gel(phi1, 4) = utoineg(0x1c7bUL);
    2538        3219 :   gel(phi1, 5) = utoi(0x7980UL);
    2539        3219 :   gel(phi1, 6) = utoi(0x12168UL);
    2540        3219 :   gel(phi1, 7) = utoineg(0x3528UL);
    2541        3219 :   gel(phi1, 8) = utoineg(0x6174UL);
    2542        3219 :   gel(phi1, 9) = utoi(0x2208UL);
    2543        3219 :   gel(phi1, 10) = utoineg(0x41dUL);
    2544        3219 :   gel(phi1, 11) = utoi(0x36UL);
    2545        3219 :   gel(phi1, 12) = gen_m1;
    2546             : 
    2547        3219 :   gel(phi, 1) = phi0;
    2548        3219 :   gel(phi, 2) = phi1;
    2549        3219 :   gel(phi, 3) = gen_2; return phi;
    2550             : }
    2551             : 
    2552             : static GEN
    2553        3234 : phi_w2w5_j(void)
    2554             : {
    2555             :   GEN phi, phi0, phi1;
    2556        3234 :   phi = cgetg(4, t_VEC);
    2557             : 
    2558        3234 :   phi0 = cgetg(20, t_VEC);
    2559        3234 :   gel(phi0, 1) = gen_1;
    2560        3234 :   gel(phi0, 2) = utoineg(0x2aUL);
    2561        3234 :   gel(phi0, 3) = utoi(0x549UL);
    2562        3234 :   gel(phi0, 4) = utoineg(0x6530UL);
    2563        3234 :   gel(phi0, 5) = utoi(0x60504UL);
    2564        3234 :   gel(phi0, 6) = utoineg(0x3cbbc8UL);
    2565        3234 :   gel(phi0, 7) = utoi(0x1d1ee74UL);
    2566        3234 :   gel(phi0, 8) = utoineg(0x7ef9ab0UL);
    2567        3234 :   gel(phi0, 9) = utoi(0x12b888beUL);
    2568        3234 :   gel(phi0, 10) = utoineg(0x15fa174cUL);
    2569        3234 :   gel(phi0, 11) = utoi(0x615d9feUL);
    2570        3234 :   gel(phi0, 12) = utoi(0xbeca070UL);
    2571        3234 :   gel(phi0, 13) = utoineg(0x88de74cUL);
    2572        3234 :   gel(phi0, 14) = utoineg(0x2b3a268UL);
    2573        3234 :   gel(phi0, 15) = utoi(0x24b3244UL);
    2574        3234 :   gel(phi0, 16) = utoi(0xb56270UL);
    2575        3234 :   gel(phi0, 17) = utoi(0x25989UL);
    2576        3234 :   gel(phi0, 18) = utoi(0x2a6UL);
    2577        3234 :   gel(phi0, 19) = gen_1;
    2578             : 
    2579        3234 :   phi1 = cgetg(19, t_VEC);
    2580        3234 :   gel(phi1, 1) = gen_0;
    2581        3234 :   gel(phi1, 2) = gen_0;
    2582        3234 :   gel(phi1, 3) = gen_m1;
    2583        3234 :   gel(phi1, 4) = utoi(0x1eUL);
    2584        3234 :   gel(phi1, 5) = utoineg(0xffUL);
    2585        3234 :   gel(phi1, 6) = utoi(0x243UL);
    2586        3234 :   gel(phi1, 7) = utoineg(0xf3UL);
    2587        3234 :   gel(phi1, 8) = utoineg(0x5c4UL);
    2588        3234 :   gel(phi1, 9) = utoi(0x107bUL);
    2589        3234 :   gel(phi1, 10) = utoineg(0x11b2fUL);
    2590        3234 :   gel(phi1, 11) = utoi(0x48fa8UL);
    2591        3234 :   gel(phi1, 12) = utoineg(0x6ff7cUL);
    2592        3234 :   gel(phi1, 13) = utoi(0x4bf48UL);
    2593        3234 :   gel(phi1, 14) = utoineg(0x187efUL);
    2594        3234 :   gel(phi1, 15) = utoi(0x404cUL);
    2595        3234 :   gel(phi1, 16) = utoineg(0x582UL);
    2596        3234 :   gel(phi1, 17) = utoi(0x3cUL);
    2597        3234 :   gel(phi1, 18) = gen_m1;
    2598             : 
    2599        3234 :   gel(phi, 1) = phi0;
    2600        3234 :   gel(phi, 2) = phi1;
    2601        3234 :   gel(phi, 3) = utoi(7); return phi;
    2602             : }
    2603             : 
    2604             : static GEN
    2605        6087 : phi_w2w7_j(void)
    2606             : {
    2607             :   GEN phi, phi0, phi1;
    2608        6087 :   phi = cgetg(4, t_VEC);
    2609             : 
    2610        6087 :   phi0 = cgetg(26, t_VEC);
    2611        6087 :   gel(phi0, 1) = gen_1;
    2612        6087 :   gel(phi0, 2) = utoineg(0x24UL);
    2613        6087 :   gel(phi0, 3) = utoi(0x4ceUL);
    2614        6087 :   gel(phi0, 4) = utoineg(0x5d60UL);
    2615        6087 :   gel(phi0, 5) = utoi(0x62b05UL);
    2616        6087 :   gel(phi0, 6) = utoineg(0x47be78UL);
    2617        6087 :   gel(phi0, 7) = utoi(0x2a3880aUL);
    2618        6087 :   gel(phi0, 8) = utoineg(0x114bccf4UL);
    2619        6087 :   gel(phi0, 9) = utoi(0x4b95e79aUL);
    2620        6087 :   gel(phi0, 10) = utoineg(0xe2cfee1cUL);
    2621        6087 :   gel(phi0, 11) = uu32toi(0x1UL, 0xe43d1126UL);
    2622        6087 :   gel(phi0, 12) = uu32toineg(0x2UL, 0xf04dc6f8UL);
    2623        6087 :   gel(phi0, 13) = uu32toi(0x3UL, 0x5384987dUL);
    2624        6087 :   gel(phi0, 14) = uu32toineg(0x2UL, 0xa5ccbe18UL);
    2625        6087 :   gel(phi0, 15) = uu32toi(0x1UL, 0x4c52c8a6UL);
    2626        6087 :   gel(phi0, 16) = utoineg(0x2643fdecUL);
    2627        6087 :   gel(phi0, 17) = utoineg(0x49f5ab66UL);
    2628        6087 :   gel(phi0, 18) = utoi(0x33074d3cUL);
    2629        6087 :   gel(phi0, 19) = utoineg(0x6a3e376UL);
    2630        6087 :   gel(phi0, 20) = utoineg(0x675aa58UL);
    2631        6087 :   gel(phi0, 21) = utoi(0x2674005UL);
    2632        6087 :   gel(phi0, 22) = utoi(0xba5be0UL);
    2633        6087 :   gel(phi0, 23) = utoi(0x2644eUL);
    2634        6087 :   gel(phi0, 24) = utoi(0x2acUL);
    2635        6087 :   gel(phi0, 25) = gen_1;
    2636             : 
    2637        6087 :   phi1 = cgetg(25, t_VEC);
    2638        6087 :   gel(phi1, 1) = gen_0;
    2639        6087 :   gel(phi1, 2) = gen_0;
    2640        6087 :   gel(phi1, 3) = gen_m1;
    2641        6087 :   gel(phi1, 4) = utoi(0x1cUL);
    2642        6087 :   gel(phi1, 5) = utoineg(0x10aUL);
    2643        6087 :   gel(phi1, 6) = utoi(0x3f0UL);
    2644        6087 :   gel(phi1, 7) = utoineg(0x5d3UL);
    2645        6087 :   gel(phi1, 8) = utoi(0x3efUL);
    2646        6087 :   gel(phi1, 9) = utoineg(0x102UL);
    2647        6087 :   gel(phi1, 10) = utoineg(0x5c8UL);
    2648        6087 :   gel(phi1, 11) = utoi(0x102fUL);
    2649        6087 :   gel(phi1, 12) = utoineg(0x13f8aUL);
    2650        6087 :   gel(phi1, 13) = utoi(0x86538UL);
    2651        6087 :   gel(phi1, 14) = utoineg(0x1bbd10UL);
    2652        6087 :   gel(phi1, 15) = utoi(0x3614e8UL);
    2653        6087 :   gel(phi1, 16) = utoineg(0x42f793UL);
    2654        6087 :   gel(phi1, 17) = utoi(0x364698UL);
    2655        6087 :   gel(phi1, 18) = utoineg(0x1c7a10UL);
    2656        6087 :   gel(phi1, 19) = utoi(0x97cc8UL);
    2657        6087 :   gel(phi1, 20) = utoineg(0x1fc8aUL);
    2658        6087 :   gel(phi1, 21) = utoi(0x4210UL);
    2659        6087 :   gel(phi1, 22) = utoineg(0x524UL);
    2660        6087 :   gel(phi1, 23) = utoi(0x38UL);
    2661        6087 :   gel(phi1, 24) = gen_m1;
    2662             : 
    2663        6087 :   gel(phi, 1) = phi0;
    2664        6087 :   gel(phi, 2) = phi1;
    2665        6087 :   gel(phi, 3) = utoi(9); return phi;
    2666             : }
    2667             : 
    2668             : static GEN
    2669        2933 : phi_w2w13_j(void)
    2670             : {
    2671             :   GEN phi, phi0, phi1;
    2672        2933 :   phi = cgetg(4, t_VEC);
    2673             : 
    2674        2933 :   phi0 = cgetg(44, t_VEC);
    2675        2933 :   gel(phi0, 1) = gen_1;
    2676        2933 :   gel(phi0, 2) = utoineg(0x1eUL);
    2677        2933 :   gel(phi0, 3) = utoi(0x45fUL);
    2678        2933 :   gel(phi0, 4) = utoineg(0x5590UL);
    2679        2933 :   gel(phi0, 5) = utoi(0x64407UL);
    2680        2933 :   gel(phi0, 6) = utoineg(0x53a792UL);
    2681        2933 :   gel(phi0, 7) = utoi(0x3b21af3UL);
    2682        2933 :   gel(phi0, 8) = utoineg(0x20d056d0UL);
    2683        2933 :   gel(phi0, 9) = utoi(0xe02db4a6UL);
    2684        2933 :   gel(phi0, 10) = uu32toineg(0x4UL, 0xb23400b0UL);
    2685        2933 :   gel(phi0, 11) = uu32toi(0x14UL, 0x57fbb906UL);
    2686        2933 :   gel(phi0, 12) = uu32toineg(0x49UL, 0xcf80c00UL);
    2687        2933 :   gel(phi0, 13) = uu32toi(0xdeUL, 0x84ff421UL);
    2688        2933 :   gel(phi0, 14) = uu32toineg(0x244UL, 0xc500c156UL);
    2689        2933 :   gel(phi0, 15) = uu32toi(0x52cUL, 0x79162979UL);
    2690        2933 :   gel(phi0, 16) = uu32toineg(0xa64UL, 0x8edc5650UL);
    2691        2933 :   gel(phi0, 17) = uu32toi(0x1289UL, 0x4225bb41UL);
    2692        2933 :   gel(phi0, 18) = uu32toineg(0x1d89UL, 0x2a15229aUL);
    2693        2933 :   gel(phi0, 19) = uu32toi(0x2a3eUL, 0x4539f1ebUL);
    2694        2933 :   gel(phi0, 20) = uu32toineg(0x366aUL, 0xa5ea1130UL);
    2695        2933 :   gel(phi0, 21) = uu32toi(0x3f47UL, 0xa19fecb4UL);
    2696        2933 :   gel(phi0, 22) = uu32toineg(0x4282UL, 0x91a3c4a0UL);
    2697        2933 :   gel(phi0, 23) = uu32toi(0x3f30UL, 0xbaa305b4UL);
    2698        2933 :   gel(phi0, 24) = uu32toineg(0x3635UL, 0xd11c2530UL);
    2699        2933 :   gel(phi0, 25) = uu32toi(0x29e2UL, 0x89df27ebUL);
    2700        2933 :   gel(phi0, 26) = uu32toineg(0x1d03UL, 0x6509d48aUL);
    2701        2933 :   gel(phi0, 27) = uu32toi(0x11e2UL, 0x272cc601UL);
    2702        2933 :   gel(phi0, 28) = uu32toineg(0x9b0UL, 0xacd58ff0UL);
    2703        2933 :   gel(phi0, 29) = uu32toi(0x485UL, 0x608d7db9UL);
    2704        2933 :   gel(phi0, 30) = uu32toineg(0x1bfUL, 0xa941546UL);
    2705        2933 :   gel(phi0, 31) = uu32toi(0x82UL, 0x56e48b21UL);
    2706        2933 :   gel(phi0, 32) = uu32toineg(0x13UL, 0xc36b2340UL);
    2707        2933 :   gel(phi0, 33) = uu32toineg(0x5UL, 0x6637257aUL);
    2708        2933 :   gel(phi0, 34) = uu32toi(0x5UL, 0x40f70bd0UL);
    2709        2933 :   gel(phi0, 35) = uu32toineg(0x1UL, 0xf70842daUL);
    2710        2933 :   gel(phi0, 36) = utoi(0x53eea5f0UL);
    2711        2933 :   gel(phi0, 37) = utoi(0xda17bf3UL);
    2712        2933 :   gel(phi0, 38) = utoineg(0xaf246c2UL);
    2713        2933 :   gel(phi0, 39) = utoi(0x278f847UL);
    2714        2933 :   gel(phi0, 40) = utoi(0xbf5550UL);
    2715        2933 :   gel(phi0, 41) = utoi(0x26f1fUL);
    2716        2933 :   gel(phi0, 42) = utoi(0x2b2UL);
    2717        2933 :   gel(phi0, 43) = gen_1;
    2718             : 
    2719        2933 :   phi1 = cgetg(43, t_VEC);
    2720        2933 :   gel(phi1, 1) = gen_0;
    2721        2933 :   gel(phi1, 2) = gen_0;
    2722        2933 :   gel(phi1, 3) = gen_m1;
    2723        2933 :   gel(phi1, 4) = utoi(0x1aUL);
    2724        2933 :   gel(phi1, 5) = utoineg(0x111UL);
    2725        2933 :   gel(phi1, 6) = utoi(0x5e4UL);
    2726        2933 :   gel(phi1, 7) = utoineg(0x1318UL);
    2727        2933 :   gel(phi1, 8) = utoi(0x2804UL);
    2728        2933 :   gel(phi1, 9) = utoineg(0x3cd6UL);
    2729        2933 :   gel(phi1, 10) = utoi(0x467cUL);
    2730        2933 :   gel(phi1, 11) = utoineg(0x3cd6UL);
    2731        2933 :   gel(phi1, 12) = utoi(0x2804UL);
    2732        2933 :   gel(phi1, 13) = utoineg(0x1318UL);
    2733        2933 :   gel(phi1, 14) = utoi(0x5e3UL);
    2734        2933 :   gel(phi1, 15) = utoineg(0x10dUL);
    2735        2933 :   gel(phi1, 16) = utoineg(0x5ccUL);
    2736        2933 :   gel(phi1, 17) = utoi(0x100bUL);
    2737        2933 :   gel(phi1, 18) = utoineg(0x160e1UL);
    2738        2933 :   gel(phi1, 19) = utoi(0xd2cb0UL);
    2739        2933 :   gel(phi1, 20) = utoineg(0x4c85fcUL);
    2740        2933 :   gel(phi1, 21) = utoi(0x137cb98UL);
    2741        2933 :   gel(phi1, 22) = utoineg(0x3c75568UL);
    2742        2933 :   gel(phi1, 23) = utoi(0x95c69c8UL);
    2743        2933 :   gel(phi1, 24) = utoineg(0x131557bcUL);
    2744        2933 :   gel(phi1, 25) = utoi(0x20aacfd0UL);
    2745        2933 :   gel(phi1, 26) = utoineg(0x2f9164e6UL);
    2746        2933 :   gel(phi1, 27) = utoi(0x3b6a5e40UL);
    2747        2933 :   gel(phi1, 28) = utoineg(0x3ff54344UL);
    2748        2933 :   gel(phi1, 29) = utoi(0x3b6a9140UL);
    2749        2933 :   gel(phi1, 30) = utoineg(0x2f927fa6UL);
    2750        2933 :   gel(phi1, 31) = utoi(0x20ae6450UL);
    2751        2933 :   gel(phi1, 32) = utoineg(0x131cd87cUL);
    2752        2933 :   gel(phi1, 33) = utoi(0x967d1e8UL);
    2753        2933 :   gel(phi1, 34) = utoineg(0x3d48ca8UL);
    2754        2933 :   gel(phi1, 35) = utoi(0x14333b8UL);
    2755        2933 :   gel(phi1, 36) = utoineg(0x5406bcUL);
    2756        2933 :   gel(phi1, 37) = utoi(0x10c130UL);
    2757        2933 :   gel(phi1, 38) = utoineg(0x27ba1UL);
    2758        2933 :   gel(phi1, 39) = utoi(0x433cUL);
    2759        2933 :   gel(phi1, 40) = utoineg(0x4c6UL);
    2760        2933 :   gel(phi1, 41) = utoi(0x34UL);
    2761        2933 :   gel(phi1, 42) = gen_m1;
    2762             : 
    2763        2933 :   gel(phi, 1) = phi0;
    2764        2933 :   gel(phi, 2) = phi1;
    2765        2933 :   gel(phi, 3) = utoi(15); return phi;
    2766             : }
    2767             : 
    2768             : static GEN
    2769        1177 : phi_w3w5_j(void)
    2770             : {
    2771             :   GEN phi, phi0, phi1;
    2772        1177 :   phi = cgetg(4, t_VEC);
    2773             : 
    2774        1177 :   phi0 = cgetg(26, t_VEC);
    2775        1177 :   gel(phi0, 1) = gen_1;
    2776        1177 :   gel(phi0, 2) = utoi(0x18UL);
    2777        1177 :   gel(phi0, 3) = utoi(0xb4UL);
    2778        1177 :   gel(phi0, 4) = utoineg(0x178UL);
    2779        1177 :   gel(phi0, 5) = utoineg(0x2d7eUL);
    2780        1177 :   gel(phi0, 6) = utoineg(0x89b8UL);
    2781        1177 :   gel(phi0, 7) = utoi(0x35c24UL);
    2782        1177 :   gel(phi0, 8) = utoi(0x128a18UL);
    2783        1177 :   gel(phi0, 9) = utoineg(0x12a911UL);
    2784        1177 :   gel(phi0, 10) = utoineg(0xcc0190UL);
    2785        1177 :   gel(phi0, 11) = utoi(0x94368UL);
    2786        1177 :   gel(phi0, 12) = utoi(0x1439d0UL);
    2787        1177 :   gel(phi0, 13) = utoi(0x96f931cUL);
    2788        1177 :   gel(phi0, 14) = utoineg(0x1f59ff0UL);
    2789        1177 :   gel(phi0, 15) = utoi(0x20e7e8UL);
    2790        1177 :   gel(phi0, 16) = utoineg(0x25fdf150UL);
    2791        1177 :   gel(phi0, 17) = utoineg(0x7091511UL);
    2792        1177 :   gel(phi0, 18) = utoi(0x1ef52f8UL);
    2793        1177 :   gel(phi0, 19) = utoi(0x341f2de4UL);
    2794        1177 :   gel(phi0, 20) = utoi(0x25d72c28UL);
    2795        1177 :   gel(phi0, 21) = utoi(0x95d2082UL);
    2796        1177 :   gel(phi0, 22) = utoi(0xd2d828UL);
    2797        1177 :   gel(phi0, 23) = utoi(0x281f4UL);
    2798        1177 :   gel(phi0, 24) = utoi(0x2b8UL);
    2799        1177 :   gel(phi0, 25) = gen_1;
    2800             : 
    2801        1177 :   phi1 = cgetg(25, t_VEC);
    2802        1177 :   gel(phi1, 1) = gen_0;
    2803        1177 :   gel(phi1, 2) = gen_0;
    2804        1177 :   gel(phi1, 3) = gen_0;
    2805        1177 :   gel(phi1, 4) = gen_1;
    2806        1177 :   gel(phi1, 5) = utoi(0xfUL);
    2807        1177 :   gel(phi1, 6) = utoi(0x2eUL);
    2808        1177 :   gel(phi1, 7) = utoineg(0x1fUL);
    2809        1177 :   gel(phi1, 8) = utoineg(0x2dUL);
    2810        1177 :   gel(phi1, 9) = utoineg(0x5caUL);
    2811        1177 :   gel(phi1, 10) = utoineg(0x358UL);
    2812        1177 :   gel(phi1, 11) = utoi(0x2f1cUL);
    2813        1177 :   gel(phi1, 12) = utoi(0xd8eaUL);
    2814        1177 :   gel(phi1, 13) = utoineg(0x38c70UL);
    2815        1177 :   gel(phi1, 14) = utoineg(0x1a964UL);
    2816        1177 :   gel(phi1, 15) = utoi(0x93512UL);
    2817        1177 :   gel(phi1, 16) = utoineg(0x58f2UL);
    2818        1177 :   gel(phi1, 17) = utoineg(0x5af1eUL);
    2819        1177 :   gel(phi1, 18) = utoi(0x1afb8UL);
    2820        1177 :   gel(phi1, 19) = utoi(0xc084UL);
    2821        1177 :   gel(phi1, 20) = utoineg(0x7fcbUL);
    2822        1177 :   gel(phi1, 21) = utoi(0x1c89UL);
    2823        1177 :   gel(phi1, 22) = utoineg(0x32aUL);
    2824        1177 :   gel(phi1, 23) = utoi(0x2dUL);
    2825        1177 :   gel(phi1, 24) = gen_m1;
    2826             : 
    2827        1177 :   gel(phi, 1) = phi0;
    2828        1177 :   gel(phi, 2) = phi1;
    2829        1177 :   gel(phi, 3) = utoi(8); return phi;
    2830             : }
    2831             : 
    2832             : static GEN
    2833        2412 : phi_w3w7_j(void)
    2834             : {
    2835             :   GEN phi, phi0, phi1;
    2836        2412 :   phi = cgetg(4, t_VEC);
    2837             : 
    2838        2412 :   phi0 = cgetg(34, t_VEC);
    2839        2412 :   gel(phi0, 1) = gen_1;
    2840        2412 :   gel(phi0, 2) = utoineg(0x14UL);
    2841        2412 :   gel(phi0, 3) = utoi(0x82UL);
    2842        2412 :   gel(phi0, 4) = utoi(0x1f8UL);
    2843        2412 :   gel(phi0, 5) = utoineg(0x2a45UL);
    2844        2412 :   gel(phi0, 6) = utoi(0x9300UL);
    2845        2412 :   gel(phi0, 7) = utoi(0x32abeUL);
    2846        2412 :   gel(phi0, 8) = utoineg(0x19c91cUL);
    2847        2412 :   gel(phi0, 9) = utoi(0xc1ba9UL);
    2848        2412 :   gel(phi0, 10) = utoi(0x1788f68UL);
    2849        2412 :   gel(phi0, 11) = utoineg(0x2b1989cUL);
    2850        2412 :   gel(phi0, 12) = utoineg(0x7a92408UL);
    2851        2412 :   gel(phi0, 13) = utoi(0x1238d56eUL);
    2852        2412 :   gel(phi0, 14) = utoi(0x13dd66a0UL);
    2853        2412 :   gel(phi0, 15) = utoineg(0x2dbedca8UL);
    2854        2412 :   gel(phi0, 16) = utoineg(0x34282eb8UL);
    2855        2412 :   gel(phi0, 17) = utoi(0x2c2a54d2UL);
    2856        2412 :   gel(phi0, 18) = utoi(0x98db81a8UL);
    2857        2412 :   gel(phi0, 19) = utoineg(0x4088be8UL);
    2858        2412 :   gel(phi0, 20) = utoineg(0xe424a220UL);
    2859        2412 :   gel(phi0, 21) = utoineg(0x67bbb232UL);
    2860        2412 :   gel(phi0, 22) = utoi(0x7dd8bb98UL);
    2861        2412 :   gel(phi0, 23) = uu32toi(0x1UL, 0xcaff744UL);
    2862        2412 :   gel(phi0, 24) = utoineg(0x1d46a378UL);
    2863        2412 :   gel(phi0, 25) = utoineg(0x82fa50f7UL);
    2864        2412 :   gel(phi0, 26) = utoineg(0x700ef38cUL);
    2865        2412 :   gel(phi0, 27) = utoi(0x20aa202eUL);
    2866        2412 :   gel(phi0, 28) = utoi(0x299b3440UL);
    2867        2412 :   gel(phi0, 29) = utoi(0xa476c4bUL);
    2868        2412 :   gel(phi0, 30) = utoi(0xd80558UL);
    2869        2412 :   gel(phi0, 31) = utoi(0x28a32UL);
    2870        2412 :   gel(phi0, 32) = utoi(0x2bcUL);
    2871        2412 :   gel(phi0, 33) = gen_1;
    2872             : 
    2873        2412 :   phi1 = cgetg(33, t_VEC);
    2874        2412 :   gel(phi1, 1) = gen_0;
    2875        2412 :   gel(phi1, 2) = gen_0;
    2876        2412 :   gel(phi1, 3) = gen_0;
    2877        2412 :   gel(phi1, 4) = gen_m1;
    2878        2412 :   gel(phi1, 5) = utoi(0xeUL);
    2879        2412 :   gel(phi1, 6) = utoineg(0x31UL);
    2880        2412 :   gel(phi1, 7) = utoineg(0xeUL);
    2881        2412 :   gel(phi1, 8) = utoi(0x99UL);
    2882        2412 :   gel(phi1, 9) = utoineg(0x8UL);
    2883        2412 :   gel(phi1, 10) = utoineg(0x2eUL);
    2884        2412 :   gel(phi1, 11) = utoineg(0x5ccUL);
    2885        2412 :   gel(phi1, 12) = utoi(0x308UL);
    2886        2412 :   gel(phi1, 13) = utoi(0x2904UL);
    2887        2412 :   gel(phi1, 14) = utoineg(0x15700UL);
    2888        2412 :   gel(phi1, 15) = utoineg(0x2b9ecUL);
    2889        2412 :   gel(phi1, 16) = utoi(0xf0966UL);
    2890        2412 :   gel(phi1, 17) = utoi(0xb3cc8UL);
    2891        2412 :   gel(phi1, 18) = utoineg(0x38241cUL);
    2892        2412 :   gel(phi1, 19) = utoineg(0x8604cUL);
    2893        2412 :   gel(phi1, 20) = utoi(0x578a64UL);
    2894        2412 :   gel(phi1, 21) = utoineg(0x11a798UL);
    2895        2412 :   gel(phi1, 22) = utoineg(0x39c85eUL);
    2896        2412 :   gel(phi1, 23) = utoi(0x1a5084UL);
    2897        2412 :   gel(phi1, 24) = utoi(0xcdeb4UL);
    2898        2412 :   gel(phi1, 25) = utoineg(0xb0364UL);
    2899        2412 :   gel(phi1, 26) = utoi(0x129d4UL);
    2900        2412 :   gel(phi1, 27) = utoi(0x126fcUL);
    2901        2412 :   gel(phi1, 28) = utoineg(0x8649UL);
    2902        2412 :   gel(phi1, 29) = utoi(0x1aa2UL);
    2903        2412 :   gel(phi1, 30) = utoineg(0x2dfUL);
    2904        2412 :   gel(phi1, 31) = utoi(0x2aUL);
    2905        2412 :   gel(phi1, 32) = gen_m1;
    2906             : 
    2907        2412 :   gel(phi, 1) = phi0;
    2908        2412 :   gel(phi, 2) = phi1;
    2909        2412 :   gel(phi, 3) = utoi(10); return phi;
    2910             : }
    2911             : 
    2912             : static GEN
    2913         210 : phi_w3w13_j(void)
    2914             : {
    2915             :   GEN phi, phi0, phi1;
    2916         210 :   phi = cgetg(4, t_VEC);
    2917             : 
    2918         210 :   phi0 = cgetg(58, t_VEC);
    2919         210 :   gel(phi0, 1) = gen_1;
    2920         210 :   gel(phi0, 2) = utoineg(0x10UL);
    2921         210 :   gel(phi0, 3) = utoi(0x58UL);
    2922         210 :   gel(phi0, 4) = utoi(0x258UL);
    2923         210 :   gel(phi0, 5) = utoineg(0x270cUL);
    2924         210 :   gel(phi0, 6) = utoi(0x9c00UL);
    2925         210 :   gel(phi0, 7) = utoi(0x2b40cUL);
    2926         210 :   gel(phi0, 8) = utoineg(0x20e250UL);
    2927         210 :   gel(phi0, 9) = utoi(0x4f46baUL);
    2928         210 :   gel(phi0, 10) = utoi(0x1869448UL);
    2929         210 :   gel(phi0, 11) = utoineg(0xa49ab68UL);
    2930         210 :   gel(phi0, 12) = utoi(0x96c7630UL);
    2931         210 :   gel(phi0, 13) = utoi(0x4f7e0af6UL);
    2932         210 :   gel(phi0, 14) = utoineg(0xea093590UL);
    2933         210 :   gel(phi0, 15) = utoineg(0x6735bc50UL);
    2934         210 :   gel(phi0, 16) = uu32toi(0x5UL, 0x971a2e08UL);
    2935         210 :   gel(phi0, 17) = uu32toineg(0x6UL, 0x29c9d965UL);
    2936         210 :   gel(phi0, 18) = uu32toineg(0xdUL, 0xeb9aa360UL);
    2937         210 :   gel(phi0, 19) = uu32toi(0x26UL, 0xe9c0584UL);
    2938         210 :   gel(phi0, 20) = uu32toineg(0x1UL, 0xb0cadce8UL);
    2939         210 :   gel(phi0, 21) = uu32toineg(0x62UL, 0x73586014UL);
    2940         210 :   gel(phi0, 22) = uu32toi(0x66UL, 0xaf672e38UL);
    2941         210 :   gel(phi0, 23) = uu32toi(0x6bUL, 0x93c28cdcUL);
    2942         210 :   gel(phi0, 24) = uu32toineg(0x11eUL, 0x4f633080UL);
    2943         210 :   gel(phi0, 25) = uu32toi(0x3cUL, 0xcc42461bUL);
    2944         210 :   gel(phi0, 26) = uu32toi(0x17bUL, 0xdec0a78UL);
    2945         210 :   gel(phi0, 27) = uu32toineg(0x166UL, 0x910d8bd0UL);
    2946         210 :   gel(phi0, 28) = uu32toineg(0xd4UL, 0x47873030UL);
    2947         210 :   gel(phi0, 29) = uu32toi(0x204UL, 0x811828baUL);
    2948         210 :   gel(phi0, 30) = uu32toineg(0x50UL, 0x5d713960UL);
    2949         210 :   gel(phi0, 31) = uu32toineg(0x198UL, 0xa27e42b0UL);
    2950         210 :   gel(phi0, 32) = uu32toi(0xe1UL, 0x25685138UL);
    2951         210 :   gel(phi0, 33) = uu32toi(0xe3UL, 0xaa5774bbUL);
    2952         210 :   gel(phi0, 34) = uu32toineg(0xcfUL, 0x392a9a00UL);
    2953         210 :   gel(phi0, 35) = uu32toineg(0x81UL, 0xfb334d04UL);
    2954         210 :   gel(phi0, 36) = uu32toi(0xabUL, 0x59594a68UL);
    2955         210 :   gel(phi0, 37) = uu32toi(0x42UL, 0x356993acUL);
    2956         210 :   gel(phi0, 38) = uu32toineg(0x86UL, 0x307ba678UL);
    2957         210 :   gel(phi0, 39) = uu32toineg(0xbUL, 0x7a9e59dcUL);
    2958         210 :   gel(phi0, 40) = uu32toi(0x4cUL, 0x27935f20UL);
    2959         210 :   gel(phi0, 41) = uu32toineg(0x2UL, 0xe0ac9045UL);
    2960         210 :   gel(phi0, 42) = uu32toineg(0x24UL, 0x14495758UL);
    2961         210 :   gel(phi0, 43) = utoi(0x20973410UL);
    2962         210 :   gel(phi0, 44) = uu32toi(0x13UL, 0x99ff4e00UL);
    2963         210 :   gel(phi0, 45) = uu32toineg(0x1UL, 0xa710d34aUL);
    2964         210 :   gel(phi0, 46) = uu32toineg(0x7UL, 0xfe5405c0UL);
    2965         210 :   gel(phi0, 47) = uu32toi(0x1UL, 0xcdee0f8UL);
    2966         210 :   gel(phi0, 48) = uu32toi(0x2UL, 0x660c92a8UL);
    2967         210 :   gel(phi0, 49) = utoi(0x3f13a35aUL);
    2968         210 :   gel(phi0, 50) = utoineg(0xe4eb4ba0UL);
    2969         210 :   gel(phi0, 51) = utoineg(0x6420f4UL);
    2970         210 :   gel(phi0, 52) = utoi(0x2c624370UL);
    2971         210 :   gel(phi0, 53) = utoi(0xb31b814UL);
    2972         210 :   gel(phi0, 54) = utoi(0xdd3ad8UL);
    2973         210 :   gel(phi0, 55) = utoi(0x29278UL);
    2974         210 :   gel(phi0, 56) = utoi(0x2c0UL);
    2975         210 :   gel(phi0, 57) = gen_1;
    2976             : 
    2977         210 :   phi1 = cgetg(57, t_VEC);
    2978         210 :   gel(phi1, 1) = gen_0;
    2979         210 :   gel(phi1, 2) = gen_0;
    2980         210 :   gel(phi1, 3) = gen_0;
    2981         210 :   gel(phi1, 4) = gen_m1;
    2982         210 :   gel(phi1, 5) = utoi(0xdUL);
    2983         210 :   gel(phi1, 6) = utoineg(0x34UL);
    2984         210 :   gel(phi1, 7) = utoi(0x1aUL);
    2985         210 :   gel(phi1, 8) = utoi(0xf7UL);
    2986         210 :   gel(phi1, 9) = utoineg(0x16cUL);
    2987         210 :   gel(phi1, 10) = utoineg(0xddUL);
    2988         210 :   gel(phi1, 11) = utoi(0x28aUL);
    2989         210 :   gel(phi1, 12) = utoineg(0xddUL);
    2990         210 :   gel(phi1, 13) = utoineg(0x16cUL);
    2991         210 :   gel(phi1, 14) = utoi(0xf6UL);
    2992         210 :   gel(phi1, 15) = utoi(0x1dUL);
    2993         210 :   gel(phi1, 16) = utoineg(0x31UL);
    2994         210 :   gel(phi1, 17) = utoineg(0x5ceUL);
    2995         210 :   gel(phi1, 18) = utoi(0x2e4UL);
    2996         210 :   gel(phi1, 19) = utoi(0x252cUL);
    2997         210 :   gel(phi1, 20) = utoineg(0x1b34cUL);
    2998         210 :   gel(phi1, 21) = utoi(0xaf80UL);
    2999         210 :   gel(phi1, 22) = utoi(0x1cc5f9UL);
    3000         210 :   gel(phi1, 23) = utoineg(0x3e1aa5UL);
    3001         210 :   gel(phi1, 24) = utoineg(0x86d17aUL);
    3002         210 :   gel(phi1, 25) = utoi(0x2427264UL);
    3003         210 :   gel(phi1, 26) = utoineg(0x691c1fUL);
    3004         210 :   gel(phi1, 27) = utoineg(0x862ad4eUL);
    3005         210 :   gel(phi1, 28) = utoi(0xab21e1fUL);
    3006         210 :   gel(phi1, 29) = utoi(0xbc19ddcUL);
    3007         210 :   gel(phi1, 30) = utoineg(0x24331db8UL);
    3008         210 :   gel(phi1, 31) = utoi(0x972c105UL);
    3009         210 :   gel(phi1, 32) = utoi(0x363d7107UL);
    3010         210 :   gel(phi1, 33) = utoineg(0x39696450UL);
    3011         210 :   gel(phi1, 34) = utoineg(0x1bce7c48UL);
    3012         210 :   gel(phi1, 35) = utoi(0x552ecba0UL);
    3013         210 :   gel(phi1, 36) = utoineg(0x1c7771b8UL);
    3014         210 :   gel(phi1, 37) = utoineg(0x393029b8UL);
    3015         210 :   gel(phi1, 38) = utoi(0x3755be97UL);
    3016         210 :   gel(phi1, 39) = utoi(0x83402a9UL);
    3017         210 :   gel(phi1, 40) = utoineg(0x24d5be62UL);
    3018         210 :   gel(phi1, 41) = utoi(0xdb6d90aUL);
    3019         210 :   gel(phi1, 42) = utoi(0xa0ef177UL);
    3020         210 :   gel(phi1, 43) = utoineg(0x99ff162UL);
    3021         210 :   gel(phi1, 44) = utoi(0xb09e27UL);
    3022         210 :   gel(phi1, 45) = utoi(0x26a7adcUL);
    3023         210 :   gel(phi1, 46) = utoineg(0x116e2fcUL);
    3024         210 :   gel(phi1, 47) = utoineg(0x1383b5UL);
    3025         210 :   gel(phi1, 48) = utoi(0x35a9e7UL);
    3026         210 :   gel(phi1, 49) = utoineg(0x1082a0UL);
    3027         210 :   gel(phi1, 50) = utoineg(0x4696UL);
    3028         210 :   gel(phi1, 51) = utoi(0x19f98UL);
    3029         210 :   gel(phi1, 52) = utoineg(0x8bb3UL);
    3030         210 :   gel(phi1, 53) = utoi(0x18bbUL);
    3031         210 :   gel(phi1, 54) = utoineg(0x297UL);
    3032         210 :   gel(phi1, 55) = utoi(0x27UL);
    3033         210 :   gel(phi1, 56) = gen_m1;
    3034             : 
    3035         210 :   gel(phi, 1) = phi0;
    3036         210 :   gel(phi, 2) = phi1;
    3037         210 :   gel(phi, 3) = utoi(16); return phi;
    3038             : }
    3039             : 
    3040             : static GEN
    3041        2902 : phi_w5w7_j(void)
    3042             : {
    3043             :   GEN phi, phi0, phi1;
    3044        2902 :   phi = cgetg(4, t_VEC);
    3045             : 
    3046        2902 :   phi0 = cgetg(50, t_VEC);
    3047        2902 :   gel(phi0, 1) = gen_1;
    3048        2902 :   gel(phi0, 2) = utoi(0xcUL);
    3049        2902 :   gel(phi0, 3) = utoi(0x2aUL);
    3050        2902 :   gel(phi0, 4) = utoi(0x10UL);
    3051        2902 :   gel(phi0, 5) = utoineg(0x69UL);
    3052        2902 :   gel(phi0, 6) = utoineg(0x318UL);
    3053        2902 :   gel(phi0, 7) = utoineg(0x148aUL);
    3054        2902 :   gel(phi0, 8) = utoineg(0x17c4UL);
    3055        2902 :   gel(phi0, 9) = utoi(0x1a73UL);
    3056        2902 :   gel(phi0, 10) = gen_0;
    3057        2902 :   gel(phi0, 11) = utoi(0x338a0UL);
    3058        2902 :   gel(phi0, 12) = utoi(0x61698UL);
    3059        2902 :   gel(phi0, 13) = utoineg(0x96e8UL);
    3060        2902 :   gel(phi0, 14) = utoi(0x140910UL);
    3061        2902 :   gel(phi0, 15) = utoineg(0x45f6b4UL);
    3062        2902 :   gel(phi0, 16) = utoineg(0x309f50UL);
    3063        2902 :   gel(phi0, 17) = utoineg(0xef9f8bUL);
    3064        2902 :   gel(phi0, 18) = utoineg(0x283167cUL);
    3065        2902 :   gel(phi0, 19) = utoi(0x625e20aUL);
    3066        2902 :   gel(phi0, 20) = utoineg(0x16186350UL);
    3067        2902 :   gel(phi0, 21) = utoi(0x46861281UL);
    3068        2902 :   gel(phi0, 22) = utoineg(0x754b96a0UL);
    3069        2902 :   gel(phi0, 23) = uu32toi(0x1UL, 0x421ca02aUL);
    3070        2902 :   gel(phi0, 24) = uu32toineg(0x2UL, 0xdb76a5cUL);
    3071        2902 :   gel(phi0, 25) = uu32toi(0x4UL, 0xf6afd8eUL);
    3072        2902 :   gel(phi0, 26) = uu32toineg(0x6UL, 0xaafd3cb4UL);
    3073        2902 :   gel(phi0, 27) = uu32toi(0x8UL, 0xda2539caUL);
    3074        2902 :   gel(phi0, 28) = uu32toineg(0xfUL, 0x84343790UL);
    3075        2902 :   gel(phi0, 29) = uu32toi(0xfUL, 0x914ff421UL);
    3076        2902 :   gel(phi0, 30) = uu32toineg(0x19UL, 0x3c123950UL);
    3077        2902 :   gel(phi0, 31) = uu32toi(0x15UL, 0x381f722aUL);
    3078        2902 :   gel(phi0, 32) = uu32toineg(0x15UL, 0xe01c0c24UL);
    3079        2902 :   gel(phi0, 33) = uu32toi(0x19UL, 0x3360b375UL);
    3080        2902 :   gel(phi0, 34) = utoineg(0x59fda9c0UL);
    3081        2902 :   gel(phi0, 35) = uu32toi(0x20UL, 0xff55024cUL);
    3082        2902 :   gel(phi0, 36) = uu32toi(0x16UL, 0xcc600800UL);
    3083        2902 :   gel(phi0, 37) = uu32toi(0x24UL, 0x1879c898UL);
    3084        2902 :   gel(phi0, 38) = uu32toi(0x1cUL, 0x37f97498UL);
    3085        2902 :   gel(phi0, 39) = uu32toi(0x19UL, 0x39ec4b60UL);
    3086        2902 :   gel(phi0, 40) = uu32toi(0x10UL, 0x52c660d0UL);
    3087        2902 :   gel(phi0, 41) = uu32toi(0x9UL, 0xcab00333UL);
    3088        2902 :   gel(phi0, 42) = uu32toi(0x4UL, 0x7fe69be4UL);
    3089        2902 :   gel(phi0, 43) = uu32toi(0x1UL, 0xa0c6f116UL);
    3090        2902 :   gel(phi0, 44) = utoi(0x69244638UL);
    3091        2902 :   gel(phi0, 45) = utoi(0xed560f7UL);
    3092        2902 :   gel(phi0, 46) = utoi(0xe7b660UL);
    3093        2902 :   gel(phi0, 47) = utoi(0x29d8aUL);
    3094        2902 :   gel(phi0, 48) = utoi(0x2c4UL);
    3095        2902 :   gel(phi0, 49) = gen_1;
    3096             : 
    3097        2902 :   phi1 = cgetg(49, t_VEC);
    3098        2902 :   gel(phi1, 1) = gen_0;
    3099        2902 :   gel(phi1, 2) = gen_0;
    3100        2902 :   gel(phi1, 3) = gen_0;
    3101        2902 :   gel(phi1, 4) = gen_0;
    3102        2902 :   gel(phi1, 5) = gen_0;
    3103        2902 :   gel(phi1, 6) = gen_1;
    3104        2902 :   gel(phi1, 7) = utoi(0x7UL);
    3105        2902 :   gel(phi1, 8) = utoi(0x8UL);
    3106        2902 :   gel(phi1, 9) = utoineg(0x9UL);
    3107        2902 :   gel(phi1, 10) = gen_0;
    3108        2902 :   gel(phi1, 11) = utoineg(0x13UL);
    3109        2902 :   gel(phi1, 12) = utoineg(0x7UL);
    3110        2902 :   gel(phi1, 13) = utoineg(0x5ceUL);
    3111        2902 :   gel(phi1, 14) = utoineg(0xb0UL);
    3112        2901 :   gel(phi1, 15) = utoi(0x460UL);
    3113        2902 :   gel(phi1, 16) = utoineg(0x194bUL);
    3114        2902 :   gel(phi1, 17) = utoi(0x87c3UL);
    3115        2902 :   gel(phi1, 18) = utoi(0x3cdeUL);
    3116        2902 :   gel(phi1, 19) = utoineg(0xd683UL);
    3117        2902 :   gel(phi1, 20) = utoi(0x6099bUL);
    3118        2902 :   gel(phi1, 21) = utoineg(0x111ea8UL);
    3119        2902 :   gel(phi1, 22) = utoi(0xfa113UL);
    3120        2901 :   gel(phi1, 23) = utoineg(0x1a6561UL);
    3121        2902 :   gel(phi1, 24) = utoineg(0x1e997UL);
    3122        2902 :   gel(phi1, 25) = utoi(0x214e54UL);
    3123        2902 :   gel(phi1, 26) = utoineg(0x29c3f4UL);
    3124        2902 :   gel(phi1, 27) = utoi(0x67e102UL);
    3125        2902 :   gel(phi1, 28) = utoineg(0x227eaaUL);
    3126        2902 :   gel(phi1, 29) = utoi(0x191d10UL);
    3127        2902 :   gel(phi1, 30) = utoi(0x1a9cd5UL);
    3128        2901 :   gel(phi1, 31) = utoineg(0x58386fUL);
    3129        2902 :   gel(phi1, 32) = utoi(0x2e49f6UL);
    3130        2902 :   gel(phi1, 33) = utoineg(0x31194bUL);
    3131        2902 :   gel(phi1, 34) = utoi(0x9e07aUL);
    3132        2902 :   gel(phi1, 35) = utoi(0x260d59UL);
    3133        2902 :   gel(phi1, 36) = utoineg(0x189921UL);
    3134        2902 :   gel(phi1, 37) = utoi(0xeca4aUL);
    3135        2902 :   gel(phi1, 38) = utoineg(0xa3d9cUL);
    3136        2902 :   gel(phi1, 39) = utoineg(0x426daUL);
    3137        2902 :   gel(phi1, 40) = utoi(0x91875UL);
    3138        2902 :   gel(phi1, 41) = utoineg(0x3b55bUL);
    3139        2902 :   gel(phi1, 42) = utoineg(0x56f4UL);
    3140        2902 :   gel(phi1, 43) = utoi(0xcd1bUL);
    3141        2902 :   gel(phi1, 44) = utoineg(0x5159UL);
    3142        2902 :   gel(phi1, 45) = utoi(0x10f4UL);
    3143        2902 :   gel(phi1, 46) = utoineg(0x20dUL);
    3144        2902 :   gel(phi1, 47) = utoi(0x23UL);
    3145        2902 :   gel(phi1, 48) = gen_m1;
    3146             : 
    3147        2902 :   gel(phi, 1) = phi0;
    3148        2902 :   gel(phi, 2) = phi1;
    3149        2902 :   gel(phi, 3) = utoi(12); return phi;
    3150             : }
    3151             : 
    3152             : GEN
    3153       23252 : double_eta_raw(long inv)
    3154             : {
    3155       23252 :   switch (inv) {
    3156             :     case INV_W2W3:
    3157        1078 :     case INV_W2W3E2: return phi_w2w3_j();
    3158             :     case INV_W3W3:
    3159        3219 :     case INV_W3W3E2: return phi_w3w3_j();
    3160             :     case INV_W2W5:
    3161        3234 :     case INV_W2W5E2: return phi_w2w5_j();
    3162             :     case INV_W2W7:
    3163        6087 :     case INV_W2W7E2: return phi_w2w7_j();
    3164        1177 :     case INV_W3W5:   return phi_w3w5_j();
    3165        2412 :     case INV_W3W7:   return phi_w3w7_j();
    3166        2933 :     case INV_W2W13:  return phi_w2w13_j();
    3167         210 :     case INV_W3W13:  return phi_w3w13_j();
    3168        2902 :     case INV_W5W7:   return phi_w5w7_j();
    3169             :     default: pari_err_BUG("double_eta_raw"); return NULL;/*LCOV_EXCL_LINE*/
    3170             :   }
    3171             : }
    3172             : 
    3173             : /**
    3174             :  * SECTION: Select discriminant for given modpoly level.
    3175             :  */
    3176             : 
    3177             : /* require an L1, useful for multi-threading */
    3178             : #define MODPOLY_USE_L1    1
    3179             : /* no bound on L1 other than the fixed bound MAX_L1 - needed to
    3180             :  * handle small L for certain invariants (but not for j) */
    3181             : #define MODPOLY_NO_MAX_L1 2
    3182             : /* don't use any auxilliary primes - needed to handle small L for
    3183             :  * certain invariants (but not for j) */
    3184             : #define MODPOLY_NO_AUX_L  4
    3185             : #define MODPOLY_IGNORE_SPARSE_FACTOR 8
    3186             : 
    3187             : INLINE double
    3188        2699 : modpoly_height_bound(long L, long inv)
    3189             : {
    3190             :   double nbits, nbits2;
    3191             :   double c;
    3192             :   long hf;
    3193             : 
    3194             :   /* proven bound (in bits), derived from: 6l*log(l)+16*l+13*sqrt(l)*log(l) */
    3195        2699 :   nbits = 6.0*L*log2(L)+16/M_LN2*L+8.0*sqrt((double)L)*log2(L);
    3196             :   /* alternative proven bound (in bits), derived from: 6l*log(l)+17*l */
    3197        2699 :   nbits2 = 6.0*L*log2(L)+17/M_LN2*L;
    3198        2699 :   if ( nbits2 < nbits ) nbits = nbits2;
    3199        2699 :   hf = modinv_height_factor(inv);
    3200        2699 :   if (hf > 1) {
    3201             :    /* IMPORTANT: when dividing by the height factor, we only want to reduce
    3202             :    terms related to the bound on j (the roots of Phi_l(X,y)), not terms arising
    3203             :    from binomial coefficients. These arise in lemmas 2 and 3 of the height
    3204             :    bound paper, terms of (log 2)*L and 2.085*(L+1) which we convert here to
    3205             :    binary logs */
    3206             :     /* Massive overestimate: if you care about speed, determine a good height
    3207             :      * bound empirically as done for INV_F below */
    3208        1634 :     nbits2 = nbits - 4.01*L -3.0;
    3209        1634 :     nbits = nbits2/hf + 4.01*L + 3.0;
    3210             :   }
    3211        2699 :   if (inv == INV_F) {
    3212         149 :     if (L < 30) c = 45;
    3213          35 :     else if (L < 100) c = 36;
    3214          21 :     else if (L < 300) c = 32;
    3215           7 :     else if (L < 600) c = 26;
    3216           0 :     else if (L < 1200) c = 24;
    3217           0 :     else if (L < 2400) c = 22;
    3218           0 :     else c = 20;
    3219         149 :     nbits = (6.0*L*log2(L) + c*L)/hf;
    3220             :   }
    3221        2699 :   return nbits;
    3222             : }
    3223             : 
    3224             : /* small enough to write the factorization of a smooth in a BIL bit integer */
    3225             : #define SMOOTH_PRIMES  ((BITS_IN_LONG >> 1) - 1)
    3226             : 
    3227             : #define MAX_ATKIN 255
    3228             : 
    3229             : /* Must have primes at least up to MAX_ATKIN */
    3230             : static const long PRIMES[] = {
    3231             :     0,   2,   3,   5,   7,  11,  13,  17,  19,  23,
    3232             :    29,  31,  37,  41,  43,  47,  53,  59,  61,  67,
    3233             :    71,  73,  79,  83,  89,  97, 101, 103, 107, 109,
    3234             :   113, 127, 131, 137, 139, 149, 151, 157, 163, 167,
    3235             :   173, 179, 181, 191, 193, 197, 199, 211, 223, 227,
    3236             :   229, 233, 239, 241, 251, 257, 263, 269, 271, 277
    3237             : };
    3238             : 
    3239             : #define MAX_L1      255
    3240             : 
    3241             : typedef struct D_entry_struct {
    3242             :   ulong m;
    3243             :   long D, h;
    3244             : } D_entry;
    3245             : 
    3246             : /* Returns a form that generates the classes of norm p^2 in cl(p^2D)
    3247             :  * (i.e. one with order p-1), where p is an odd prime that splits in D
    3248             :  * and does not divide its conductor (but this is not verified) */
    3249             : INLINE GEN
    3250       61211 : qform_primeform2(long p, long D)
    3251             : {
    3252       61211 :   pari_sp ltop = avma, av;
    3253             :   GEN M;
    3254             :   register long k;
    3255             : 
    3256       61211 :   M = factor_pn_1(stoi(p), 1);
    3257       61211 :   av = avma;
    3258      124469 :   for (k = D & 1; k <= p; k += 2)
    3259             :   {
    3260             :     GEN Q;
    3261      124469 :     long ord, a, b, c = (k * k - D) / 4;
    3262             : 
    3263      124469 :     if (!(c % p)) continue;
    3264      107533 :     a = p * p;
    3265      107533 :     b = k * p;
    3266      107533 :     Q = redimag(mkqfis(a, b, c));
    3267             :     /* TODO: How do we know that Q has order dividing p - 1? If we don't, then
    3268             :      * the call to gen_order should be replaced with a call to something with
    3269             :      * fastorder semantics (i.e. return 0 if ord(Q) \ndiv M). */
    3270      107533 :     ord = itos(qfi_order(Q, M));
    3271      107533 :     if (ord == p - 1) {
    3272             :       /* TODO: This check that gen_order returned the correct result should be
    3273             :        * removed when gen_order is replaced with fastorder semantics. */
    3274       61211 :       GEN tst = gpowgs(Q, p - 1);
    3275       61211 :       if (qfb_equal1(tst)) { avma = ltop; return mkqfis(a, b, c); }
    3276           0 :         break;
    3277             :     }
    3278       46322 :     avma = av;
    3279             :   }
    3280           0 :   avma = ltop; return NULL;
    3281             : }
    3282             : 
    3283             : /* Let n = #cl(D); return x such that [L0]^x = [L] in cl(D), or -1 if x was
    3284             :  * not found */
    3285             : INLINE long
    3286      167831 : primeform_discrete_log(long L0, long L, long n, long D)
    3287             : {
    3288      167831 :   pari_sp av = avma;
    3289      167831 :   GEN X, Q, R, DD = stoi(D);
    3290      167831 :   Q = primeform_u(DD, L0);
    3291      167831 :   R = primeform_u(DD, L);
    3292      167831 :   X = qfi_Shanks(R, Q, n);
    3293      167831 :   avma = av; return X? itos(X): -1;
    3294             : }
    3295             : 
    3296             : /* Return the norm of a class group generator appropriate for a discriminant
    3297             :  * that will be used to calculate the modular polynomial of level L and
    3298             :  * invariant inv.  Don't consider norms less than initial_L0 */
    3299             : static long
    3300        2699 : select_L0(long L, long inv, long initial_L0)
    3301             : {
    3302        2699 :   long L0, modinv_N = modinv_level(inv);
    3303             : 
    3304        2699 :   if (modinv_N % L == 0) pari_err_BUG("select_L0");
    3305             : 
    3306             :   /* TODO: Clean up these anomolous L0 choices */
    3307             : 
    3308             :   /* I've no idea why the discriminant-finding code fails with L0=5
    3309             :    * when L=19 and L=29, nor why L0=7 and L0=11 don't work for L=19
    3310             :    * either, nor why this happens for the otherwise unrelated
    3311             :    * invariants Weber-f and (2,3) double-eta. */
    3312        2699 :   if (inv == INV_W3W3E2 && L == 5) return 2;
    3313             : 
    3314        2685 :   if (inv == INV_F || inv == INV_F2 || inv == INV_F4 || inv == INV_F8
    3315        2424 :       || inv == INV_W2W3 || inv == INV_W2W3E2
    3316        2361 :       || inv == INV_W3W3 /* || inv == INV_W3W3E2 */) {
    3317         408 :     if (L == 19) return 13;
    3318         365 :     else if (L == 29 || L == 5) return 7;
    3319         302 :     return 5;
    3320             :   }
    3321        2277 :   if ((inv == INV_W2W5 || inv == INV_W2W5E2)
    3322         140 :       && (L == 7 || L == 19)) return 13;
    3323        2235 :   if ((inv == INV_W2W7 || inv == INV_W2W7E2)
    3324         337 :       && L == 11) return 13;
    3325        2207 :   if (inv == INV_W3W5) {
    3326          63 :     if (L == 7) return 13;
    3327          56 :     else if (L == 17) return 7;
    3328             :   }
    3329        2200 :   if (inv == INV_W3W7) {
    3330         140 :     if (L == 29 || L == 101) return 11;
    3331         112 :     if (L == 11 || L == 19) return 13;
    3332             :   }
    3333        2144 :   if (inv == INV_W5W7 && L == 17) return 3;
    3334             : 
    3335             :   /* L0 = smallest small prime different from L that doesn't divide modinv_N */
    3336        5192 :   for (L0 = unextprime(initial_L0 + 1);
    3337        3013 :        L0 == L || modinv_N % L0 == 0;
    3338         946 :        L0 = unextprime(L0 + 1))
    3339             :     ;
    3340        2123 :   return L0;
    3341             : }
    3342             : 
    3343             : /* Return the order of [L]^n in cl(D), where #cl(D) = ord. */
    3344             : INLINE long
    3345      905402 : primeform_exp_order(long L, long n, long D, long ord)
    3346             : {
    3347      905402 :   pari_sp av = avma;
    3348      905402 :   GEN Q = gpowgs(primeform_u(stoi(D), L), n);
    3349      905402 :   long m = itos(qfi_order(Q, Z_factor(stoi(ord))));
    3350      905402 :   avma = av; return m;
    3351             : }
    3352             : 
    3353             : /* If an ideal of norm modinv_deg is equivalent to an ideal of norm L0, we
    3354             :  * have an orientation ambiguity that we need to avoid. Note that we need to
    3355             :  * check all the possibilities (up to 8), but we can cheaply check inverses
    3356             :  * (so at most 2) */
    3357             : static long
    3358       32999 : orientation_ambiguity(long D1, long L0, long modinv_p1, long modinv_p2, long modinv_N)
    3359             : {
    3360       32999 :   pari_sp av = avma;
    3361       32999 :   long ambiguity = 0;
    3362       32999 :   GEN D = stoi(D1), Q1 = primeform_u(D, modinv_p1), Q2 = NULL;
    3363             : 
    3364       32999 :   if (modinv_p2 > 1)
    3365             :   {
    3366       32999 :     if (modinv_p1 == modinv_p2) Q1 = gsqr(Q1);
    3367             :     else
    3368             :     {
    3369       27659 :       GEN P2 = primeform_u(D, modinv_p2);
    3370       27659 :       GEN Q = gsqr(P2), R = gsqr(Q1);
    3371             :       /* check that p1^2 != p2^{+/-2}, since this leads to
    3372             :        * ambiguities when converting j's to f's */
    3373       27659 :       if (equalii(gel(Q,1), gel(R,1)) && absequalii(gel(Q,2), gel(R,2)))
    3374             :       {
    3375           0 :         dbg_printf(3)("Bad D=%ld, a^2=b^2 problem between modinv_p1=%ld and modinv_p2=%ld\n",
    3376             :                       D1, modinv_p1, modinv_p2);
    3377           0 :         ambiguity = 1;
    3378             :       }
    3379             :       else
    3380             :       { /* generate both p1*p2 and p1*p2^{-1} */
    3381       27659 :         Q2 = gmul(Q1, P2);
    3382       27659 :         P2 = ginv(P2);
    3383       27659 :         Q1 = gmul(Q1, P2);
    3384             :       }
    3385             :     }
    3386             :   }
    3387       32999 :   if (!ambiguity)
    3388             :   {
    3389       32999 :     GEN P = gsqr(primeform_u(D, L0));
    3390       32999 :     if (equalii(gel(P,1), gel(Q1,1))
    3391       31983 :         || (modinv_p2 && modinv_p1 != modinv_p2
    3392       26657 :                       && equalii(gel(P,1), gel(Q2,1)))) {
    3393        1459 :       dbg_printf(3)("Bad D=%ld, a=b^{+/-2} problem between modinv_N=%ld and L0=%ld\n",
    3394             :                     D1, modinv_N, L0);
    3395        1459 :       ambiguity = 1;
    3396             :     }
    3397             :   }
    3398       32999 :   avma = av; return ambiguity;
    3399             : }
    3400             : 
    3401             : static long
    3402      652927 : check_generators(
    3403             :   long *n1_, long *m_,
    3404             :   long D, long h, long n, long subgrp_sz, long L0, long L1)
    3405             : {
    3406      652927 :   long n1, m = primeform_exp_order(L0, n, D, h);
    3407      652927 :   if (m_) *m_ = m;
    3408      652927 :   n1 = n * m;
    3409      652927 :   if (!n1) pari_err_BUG("check_generators");
    3410      652927 :   *n1_ = n1;
    3411      652927 :   if (n1 < subgrp_sz/2 || ( ! L1 && n1 < subgrp_sz))  {
    3412       24121 :     dbg_printf(3)("Bad D1=%ld with n1=%ld, h1=%ld, L1=%ld: "
    3413             :                   "L0 and L1 don't span subgroup of size d in cl(D1)\n",
    3414             :                   D, n, h, L1);
    3415       24121 :     return 0;
    3416             :   }
    3417      628806 :   if (n1 < subgrp_sz && ! (n1 & 1)) {
    3418             :     int res;
    3419             :     /* check whether L1 is generated by L0, use the fact that it has order 2 */
    3420       13598 :     pari_sp av = avma;
    3421       13598 :     GEN D1 = stoi(D);
    3422       13598 :     GEN Q = gpowgs(primeform_u(D1, L0), n1 / 2);
    3423       13598 :     res = gequal(Q, redimag(primeform_u(D1, L1)));
    3424       13598 :     avma = av;
    3425       13598 :     if (res) {
    3426        3983 :       dbg_printf(3)("Bad D1=%ld, with n1=%ld, h1=%ld, L1=%ld: "
    3427             :                     "L1 generated by L0 in cl(D1)\n", D, n, h, L1);
    3428        3983 :       return 0;
    3429             :     }
    3430             :   }
    3431      624823 :   return 1;
    3432             : }
    3433             : 
    3434             : /* Calculate solutions (p, t) to the norm equation
    3435             :  *   4 p = t^2 - v^2 L^2 D   (*)
    3436             :  * corresponding to the descriminant described by Dinfo.
    3437             :  *
    3438             :  * INPUT:
    3439             :  * - max: length of primes and traces
    3440             :  * - xprimes: p to exclude from primes (if they arise)
    3441             :  * - xcnt: length of xprimes
    3442             :  * - minbits: sum of log2(p) must be larger than this
    3443             :  * - Dinfo: discriminant, invariant and L for which we seek solutions to (*)
    3444             :  *
    3445             :  * OUTPUT:
    3446             :  * - primes: array of p in (*)
    3447             :  * - traces: array of t in (*)
    3448             :  * - totbits: sum of log2(p) for p in primes.
    3449             :  *
    3450             :  * RETURN:
    3451             :  * - the number of primes and traces found (these are always the same).
    3452             :  *
    3453             :  * NOTE: primes and traces are both NULL or both non-NULL.
    3454             :  * xprimes can be zero, in which case it is treated as empty. */
    3455             : static long
    3456        9829 : modpoly_pickD_primes(
    3457             :   ulong *primes, ulong *traces, long max, ulong *xprimes, long xcnt,
    3458             :   long *totbits, long minbits, disc_info *Dinfo)
    3459             : {
    3460             :   double bits;
    3461             :   long D, m, n, vcnt, pfilter, one_prime, inv;
    3462             :   ulong maxp;
    3463             :   register ulong a1, a2, v, t, p, a1_start, a1_delta, L0, L1, L, absD;
    3464        9829 :   ulong FF_BITS = BITS_IN_LONG - 2; /* BITS_IN_LONG - NAIL_BITS */
    3465             : 
    3466        9829 :   D = Dinfo->D1; absD = -D;
    3467        9829 :   L0 = Dinfo->L0;
    3468        9829 :   L1 = Dinfo->L1;
    3469        9829 :   L = Dinfo->L;
    3470        9829 :   inv = Dinfo->inv;
    3471             : 
    3472             :   /* make sure pfilter and D don't preclude the possibility of p=(t^2-v^2D)/4 being prime */
    3473        9829 :   pfilter = modinv_pfilter(inv);
    3474        9829 :   if ((pfilter & IQ_FILTER_1MOD3) && ! (D % 3)) return 0;
    3475        9780 :   if ((pfilter & IQ_FILTER_1MOD4) && ! (D & 0xF)) return 0;
    3476             : 
    3477             :   /* Naively estimate the number of primes satisfying 4p=t^2-L^2D with
    3478             :    * t=2 mod L and pfilter. This is roughly
    3479             :    * #{t: t^2 < max p and t=2 mod L} / pi(max p) * filter_density,
    3480             :    * where filter_density is 1, 2, or 4 depending on pfilter.  If this quantity
    3481             :    * is already more than twice the number of bits we need, assume that,
    3482             :    * barring some obstruction, we should have no problem getting enough primes.
    3483             :    * In this case we just verify we can get one prime (which should always be
    3484             :    * true, assuming we chose D properly). */
    3485        9780 :   one_prime = 0;
    3486        9780 :   *totbits = 0;
    3487        9780 :   if (max <= 1 && ! one_prime) {
    3488        7061 :     p = ((pfilter & IQ_FILTER_1MOD3) ? 2 : 1) * ((pfilter & IQ_FILTER_1MOD4) ? 2 : 1);
    3489       14122 :     one_prime = (1UL << ((FF_BITS+1)/2)) * (log2(L*L*(-D))-1)
    3490        7061 :         > p*L*minbits*FF_BITS*M_LN2;
    3491        7061 :     if (one_prime) *totbits = minbits+1;   /* lie */
    3492             :   }
    3493             : 
    3494        9780 :   m = n = 0;
    3495        9780 :   bits = 0.0;
    3496        9780 :   maxp = 0;
    3497       24395 :   for (v = 1; v < 100 && bits < minbits; v++) {
    3498             :     /* Don't allow v dividing the conductor. */
    3499       21626 :     if (ugcd(absD, v) != 1) continue;
    3500             :     /* Avoid v dividing the level. */
    3501       21528 :     if (v > 2 && modinv_is_double_eta(inv) && ugcd(modinv_level(inv), v) != 1)
    3502         966 :       continue;
    3503             :     /* can't get odd p with D=1 mod 8 unless v is even */
    3504       20562 :     if ((v & 1) && (D & 7) == 1) continue;
    3505             :     /* disallow 4 | v for L0=2 (removing this restriction is costly) */
    3506       10149 :     if (L0 == 2 && !(v & 3)) continue;
    3507             :     /* can't get p=3mod4 if v^2D is 0 mod 16 */
    3508       10038 :     if ((pfilter & IQ_FILTER_1MOD4) && !((v*v*D) & 0xF)) continue;
    3509        9955 :     if ((pfilter & IQ_FILTER_1MOD3) && !(v%3) ) continue;
    3510             :     /* avoid L0-volcanos with non-zero height */
    3511        9897 :     if (L0 != 2 && ! (v % L0)) continue;
    3512             :     /* ditto for L1 */
    3513        9876 :     if (L1 && !(v % L1)) continue;
    3514        9876 :     vcnt = 0;
    3515        9876 :     if ((v*v*absD)/4 > (1L<<FF_BITS)/(L*L)) break;
    3516        9830 :     if (both_odd(v,D)) {
    3517           0 :       a1_start = 1;
    3518           0 :       a1_delta = 2;
    3519             :     } else {
    3520        9830 :       a1_start = ((v*v*D) & 7)? 2: 0;
    3521        9830 :       a1_delta = 4;
    3522             :     }
    3523      460493 :     for (a1 = a1_start; bits < minbits; a1 += a1_delta) {
    3524      457738 :       a2 = (a1*a1 + v*v*absD) >> 2;
    3525      457738 :       if (!(a2 % L)) continue;
    3526      387157 :       t = a1*L + 2;
    3527      387157 :       p = a2*L*L + t - 1;
    3528             :       /* double check calculation just in case of overflow or other weirdness */
    3529      387157 :       if (!odd(p) || t*t + v*v*L*L*absD != 4*p)
    3530           0 :         pari_err_BUG("modpoly_pickD_primes");
    3531      387157 :       if (p > (1UL<<FF_BITS)) break;
    3532      387047 :       if (xprimes) {
    3533      307145 :         while (m < xcnt && xprimes[m] < p) m++;
    3534      307145 :         if (m < xcnt && p == xprimes[m]) {
    3535           0 :           dbg_printf(1)("skipping duplicate prime %ld\n", p);
    3536           0 :           continue;
    3537             :         }
    3538             :       }
    3539      387047 :       if (!uisprime(p) || !modinv_good_prime(inv, p)) continue;
    3540       42296 :       if (primes) {
    3541       32650 :         if (n >= max) goto done;
    3542             :         /* TODO: Implement test to filter primes that lead to
    3543             :          * L-valuation != 2 */
    3544       32650 :         primes[n] = p;
    3545       32650 :         traces[n] = t;
    3546             :       }
    3547       42296 :       n++;
    3548       42296 :       vcnt++;
    3549       42296 :       bits += log2(p);
    3550       42296 :       if (p > maxp) maxp = p;
    3551       42296 :       if (one_prime) goto done;
    3552             :     }
    3553        2865 :     if (vcnt)
    3554        2862 :       dbg_printf(3)("%ld primes with v=%ld, maxp=%ld (%.2f bits)\n",
    3555             :                  vcnt, v, maxp, log2(maxp));
    3556             :   }
    3557             : done:
    3558        9780 :   if (!n) {
    3559           9 :     dbg_printf(3)("check_primes failed completely for D=%ld\n", D);
    3560           9 :     return 0;
    3561             :   }
    3562        9771 :   dbg_printf(3)("D=%ld: Found %ld primes totalling %0.2f of %ld bits\n",
    3563             :              D, n, bits, minbits);
    3564        9771 :   if (!*totbits) *totbits = (long)bits;
    3565        9771 :   return n;
    3566             : }
    3567             : 
    3568             : #define MAX_VOLCANO_FLOOR_SIZE 100000000
    3569             : 
    3570             : static long
    3571        2701 : calc_primes_for_discriminants(disc_info Ds[], long Dcnt, long L, long minbits)
    3572             : {
    3573        2701 :   pari_sp av = avma;
    3574             :   long i, j, k, m, n, D1, pcnt, totbits;
    3575             :   ulong *primes, *Dprimes, *Dtraces;
    3576             : 
    3577             :   /* D1 is the discriminant with smallest absolute value among those we found */
    3578        2701 :   D1 = Ds[0].D1;
    3579        7052 :   for (i = 1; i < Dcnt; i++)
    3580        4351 :     if (Ds[i].D1 > D1) D1 = Ds[i].D1;
    3581             : 
    3582             :   /* n is an upper bound on the number of primes we might get. */
    3583        2701 :   n = ceil(minbits / (log2(L * L * (-D1)) - 2)) + 1;
    3584        2701 :   primes = (ulong *) stack_malloc(n * sizeof(*primes));
    3585        2701 :   Dprimes = (ulong *) stack_malloc(n * sizeof(*Dprimes));
    3586        2701 :   Dtraces = (ulong *) stack_malloc(n * sizeof(*Dtraces));
    3587        2719 :   for (i = 0, totbits = 0, pcnt = 0; i < Dcnt && totbits < minbits; i++)
    3588             :   {
    3589        5438 :     long np = modpoly_pickD_primes(Dprimes, Dtraces, n, primes, pcnt,
    3590        5438 :                                    &Ds[i].bits, minbits - totbits, Ds + i);
    3591        2719 :     ulong *T = (ulong *)newblock(2*np);
    3592        2719 :     Ds[i].nprimes = np;
    3593        2719 :     Ds[i].primes = T;    memcpy(T   , Dprimes, np * sizeof(*Dprimes));
    3594        2719 :     Ds[i].traces = T+np; memcpy(T+np, Dtraces, np * sizeof(*Dtraces));
    3595             : 
    3596        2719 :     totbits += Ds[i].bits;
    3597        2719 :     pcnt += np;
    3598             : 
    3599        2719 :     if (totbits >= minbits || i == Dcnt - 1) { Dcnt = i + 1; break; }
    3600             :     /* merge lists */
    3601         552 :     for (j = pcnt - np - 1, k = np - 1, m = pcnt - 1; m >= 0; m--) {
    3602         534 :       if (k >= 0) {
    3603         509 :         if (j >= 0 && primes[j] > Dprimes[k])
    3604         301 :           primes[m] = primes[j--];
    3605             :         else
    3606         208 :           primes[m] = Dprimes[k--];
    3607             :       } else {
    3608          25 :         primes[m] = primes[j--];
    3609             :       }
    3610             :     }
    3611             :   }
    3612        2701 :   if (totbits < minbits) {
    3613           2 :     dbg_printf(1)("Only obtained %ld of %ld bits using %ld discriminants\n",
    3614             :                   totbits, minbits, Dcnt);
    3615           2 :     for (i = 0; i < Dcnt; i++) killblock((GEN)Ds[i].primes);
    3616           2 :     Dcnt = 0;
    3617             :   }
    3618        2701 :   avma = av; return Dcnt;
    3619             : }
    3620             : 
    3621             : /* Select discriminant(s) to use when calculating the modular
    3622             :  * polynomial of level L and invariant inv.
    3623             :  *
    3624             :  * INPUT:
    3625             :  * - L: level of modular polynomial (must be odd)
    3626             :  * - inv: invariant of modular polynomial
    3627             :  * - L0: result of select_L0(L, inv)
    3628             :  * - minbits: height of modular polynomial
    3629             :  * - flags: see below
    3630             :  * - tab: result of scanD0(L0)
    3631             :  * - tablen: length of tab
    3632             :  *
    3633             :  * OUTPUT:
    3634             :  * - Ds: the selected discriminant(s)
    3635             :  *
    3636             :  * RETURN:
    3637             :  * - the number of Ds found
    3638             :  *
    3639             :  * The flags parameter is constructed by ORing zero or more of the
    3640             :  * following values:
    3641             :  * - MODPOLY_USE_L1: force use of second class group generator
    3642             :  * - MODPOLY_NO_AUX_L: don't use auxillary class group elements
    3643             :  * - MODPOLY_IGNORE_SPARSE_FACTOR: obtain D for which h(D) > L + 1
    3644             :  *   rather than h(D) > (L + 1)/s */
    3645             : static long
    3646        2701 : modpoly_pickD(disc_info Ds[MODPOLY_MAX_DCNT], long L, long inv,
    3647             :   long L0, long max_L1, long minbits, long flags, D_entry *tab, long tablen)
    3648             : {
    3649        2701 :   pari_sp ltop = avma, btop;
    3650             :   disc_info Dinfo;
    3651             :   pari_timer T;
    3652             :   long modinv_p1, modinv_p2; /* const after next line */
    3653        2701 :   const long modinv_deg = modinv_degree(&modinv_p1, &modinv_p2, inv);
    3654        2701 :   const long pfilter = modinv_pfilter(inv), modinv_N = modinv_level(inv);
    3655             :   long i, k, use_L1, Dcnt, D0_i, d, cost, enum_cost, best_cost, totbits;
    3656        2701 :   const double L_bits = log2(L);
    3657             : 
    3658        2701 :   if (!odd(L)) pari_err_BUG("modpoly_pickD");
    3659             : 
    3660        2701 :   timer_start(&T);
    3661        2701 :   if (flags & MODPOLY_IGNORE_SPARSE_FACTOR) d = L+2;
    3662        2561 :   else d = ceildivuu(L+1, modinv_sparse_factor(inv)) + 1;
    3663             : 
    3664             :   /* Now set level to 0 unless we will need to compute N-isogenies */
    3665        2701 :   dbg_printf(1)("Using L0=%ld for L=%ld, d=%ld, modinv_N=%ld, modinv_deg=%ld\n",
    3666             :                 L0, L, d, modinv_N, modinv_deg);
    3667             : 
    3668             :   /* We use L1 if (L0|L) == 1 or if we are forced to by flags. */
    3669        2701 :   use_L1 = (kross(L0,L) > 0 || (flags & MODPOLY_USE_L1));
    3670             : 
    3671        2701 :   Dcnt = best_cost = totbits = 0;
    3672        2701 :   dbg_printf(3)("use_L1=%ld\n", use_L1);
    3673        2701 :   dbg_printf(3)("minbits = %ld\n", minbits);
    3674             : 
    3675             :   /* Iterate over the fundamental discriminants for L0 */
    3676     1657002 :   for (D0_i = 0; D0_i < tablen; D0_i++)
    3677             :   {
    3678     1654301 :     D_entry D0_entry = tab[D0_i];
    3679     1654301 :     long m, n0, h0, deg, L1, H_cost, twofactor, D0 = D0_entry.D;
    3680             :     double D0_bits;
    3681     3128222 :     if (! modinv_good_disc(inv, D0)) continue;
    3682     1045380 :     dbg_printf(3)("D0=%ld\n", D0);
    3683             :     /* don't allow either modinv_p1 or modinv_p2 to ramify */
    3684     1045380 :     if (kross(D0, L) < 1
    3685      472046 :         || (modinv_p1 > 1 && kross(D0, modinv_p1) < 1)
    3686      467137 :         || (modinv_p2 > 1 && kross(D0, modinv_p2) < 1)) {
    3687      588465 :       dbg_printf(3)("Bad D0=%ld due to non-split L or ramified level\n", D0);
    3688      588465 :       continue;
    3689             :     }
    3690      456915 :     deg = D0_entry.h; /* class poly degree */
    3691      456915 :     h0 = ((D0_entry.m & 2) ? 2*deg : deg); /* class number */
    3692             :     /* (D0_entry.m & 1) is 1 if ord(L0) < h0 (hence = h0/2),
    3693             :      *                  is 0 if ord(L0) = h0 */
    3694      456915 :     n0 = h0 / ((D0_entry.m & 1) + 1); /* = ord(L0) */
    3695             : 
    3696             :     /* Look for L1: for each smooth prime p */
    3697      456915 :     L1 = 0;
    3698    11137618 :     for (i = 1 ; i <= SMOOTH_PRIMES; i++)
    3699             :     {
    3700    10770340 :       long p = PRIMES[i];
    3701    10770340 :       if (p <= L0) continue;
    3702             :       /* If 1 + (D0 | p) = 1, i.e. p | D0 */
    3703    10139843 :       if (((D0_entry.m >> (2*i)) & 3) == 1) {
    3704             :         /* XXX: Why (p | L) = -1?  Presumably so (L^2 v^2 D0 | p) = -1? */
    3705      327072 :         if (p <= max_L1 && modinv_N % p && kross(p,L) < 0) { L1 = p; break; }
    3706             :       }
    3707             :     }
    3708      456915 :     if (i > SMOOTH_PRIMES && (n0 < h0 || use_L1))
    3709             :     { /* Didn't find suitable L1 though we need one */
    3710      204440 :       dbg_printf(3)("Bad D0=%ld because there is no good L1\n", D0);
    3711      204440 :       continue;
    3712             :     }
    3713      252475 :     dbg_printf(3)("Good D0=%ld with L1=%ld, n0=%ld, h0=%ld, d=%ld\n",
    3714             :                   D0, L1, n0, h0, d);
    3715             : 
    3716             :     /* We're finished if we have sufficiently many discriminants that satisfy
    3717             :      * the cost requirement */
    3718      252475 :     if (totbits > minbits && best_cost && h0*(L-1) > 3*best_cost) break;
    3719             : 
    3720      252475 :     D0_bits = log2(-D0);
    3721             :     /* If L^2 D0 is too big to fit in a BIL bit integer, skip D0. */
    3722      252475 :     if (D0_bits + 2 * L_bits > (BITS_IN_LONG - 1)) continue;
    3723             : 
    3724             :     /* m is the order of L0^n0 in L^2 D0? */
    3725      252475 :     m = primeform_exp_order(L0, n0, L * L * D0, n0 * (L-1));
    3726      252475 :     if (m < (L-1)/2) {
    3727       72095 :       dbg_printf(3)("Bad D0=%ld because %ld is less than (L-1)/2=%ld\n",
    3728           0 :                     D0, m, (L - 1)/2);
    3729       72095 :       continue;
    3730             :     }
    3731             :     /* Heuristic.  Doesn't end up contributing much. */
    3732      180380 :     H_cost = 2 * deg * deg;
    3733             : 
    3734             :     /* 0xc = 0b1100, so D0_entry.m & 0xc == 1 + (D0 | 2) */
    3735      180380 :     if ((D0 & 7) == 5) /* D0 = 5 (mod 8) */
    3736        5337 :       twofactor = ((D0_entry.m & 0xc) ? 1 : 3);
    3737             :     else
    3738      175043 :       twofactor = 0;
    3739             : 
    3740      180380 :     btop = avma;
    3741             :     /* For each small prime... */
    3742      624601 :     for (i = 0; i <= SMOOTH_PRIMES; i++) {
    3743             :       long h1, h2, D1, D2, n1, n2, dl1, dl20, dl21, p, q, j;
    3744             :       double p_bits;
    3745      624496 :       avma = btop;
    3746             :       /* i = 0 corresponds to 1, which we do not want to skip! (i.e. DK = D) */
    3747      624496 :       if (i) {
    3748      887926 :         if (modinv_odd_conductor(inv) && i == 1) continue;
    3749      435327 :         p = PRIMES[i];
    3750             :         /* Don't allow large factors in the conductor. */
    3751      615602 :         if (p > max_L1) break;
    3752      353862 :         if (p == L0 || p == L1 || p == L || p == modinv_p1 || p == modinv_p2)
    3753      122650 :           continue;
    3754      231212 :         p_bits = log2(p);
    3755             :         /* h1 is the class number of D1 = q^2 D0, where q = p^j (j defined in the loop below) */
    3756      231212 :         h1 = h0 * (p - ((D0_entry.m >> (2*i)) & 0x3) + 1);
    3757             :         /* q is the smallest power of p such that h1 >= d ~ "L + 1". */
    3758      231212 :         for (j = 1, q = p; h1 < d; j++, q *= p, h1 *= p)
    3759             :           ;
    3760      231212 :         D1 = q * q * D0;
    3761             :         /* can't have D1 = 0 mod 16 and hope to get any primes congruent to 3 mod 4 */
    3762      231212 :         if ((pfilter & IQ_FILTER_1MOD4) && !(D1 & 0xF)) continue;
    3763             :       } else {
    3764             :         /* i = 0, corresponds to "p = 1". */
    3765      180380 :         h1 = h0;
    3766      180380 :         D1 = D0;
    3767      180380 :         p = q = j = 1;
    3768      180380 :         p_bits = 0;
    3769             :       }
    3770             :       /* include a factor of 4 if D1 is 5 mod 8 */
    3771             :       /* XXX: No idea why he does this. */
    3772      411522 :       if (twofactor && (q & 1)) {
    3773       13052 :         if (modinv_odd_conductor(inv)) continue;
    3774         518 :         D1 *= 4;
    3775         518 :         h1 *= twofactor;
    3776             :       }
    3777             :       /* heuristic early abort; we may miss good D1's, but this saves time */
    3778      398988 :       if (totbits > minbits && best_cost && h1*(L-1) > 2.2*best_cost) continue;
    3779             : 
    3780             :       /* log2(D0 * (p^j)^2 * L^2 * twofactor) > (BIL - 1) -- params too big. */
    3781      776330 :       if (D0_bits + 2*j*p_bits + 2*L_bits
    3782      389752 :           + (twofactor && (q & 1) ? 2.0 : 0.0) > (BITS_IN_LONG-1)) continue;
    3783             : 
    3784      386578 :       if (! check_generators(&n1, NULL, D1, h1, n0, d, L0, L1)) continue;
    3785             : 
    3786      372007 :       if (n1 >= h1) dl1 = -1; /* fill it in later */
    3787      165159 :       else if ((dl1 = primeform_discrete_log(L0, L, n1, D1)) < 0) continue;
    3788      267808 :       dbg_printf(3)("Good D0=%ld, D1=%ld with q=%ld, L1=%ld, n1=%ld, h1=%ld\n",
    3789             :                     D0, D1, q, L1, n1, h1);
    3790      267808 :       if (modinv_deg && orientation_ambiguity(D1, L0, modinv_p1, modinv_p2, modinv_N))
    3791        1459 :         continue;
    3792             : 
    3793      266349 :       D2 = L * L * D1;
    3794      266349 :       h2 = h1 * (L-1);
    3795             :       /* m is the order of L0^n1 in cl(D2) */
    3796      266349 :       if (!check_generators(&n2, &m, D2, h2, n1, d*(L-1), L0, L1)) continue;
    3797             : 
    3798             :       /* This restriction on m is not necessary, but simplifies life later */
    3799      252816 :       if (m < (L-1)/2 || (!L1 && m < L-1)) {
    3800      122944 :         dbg_printf(3)("Bad D2=%ld for D1=%ld, D0=%ld, with n2=%ld, h2=%ld, L1=%ld, "
    3801             :                       "order of L0^n1 in cl(D2) is too small\n", D2, D1, D0, n2, h2, L1);
    3802      122944 :         continue;
    3803             :       }
    3804      129872 :       dl20 = n1;
    3805      129872 :       dl21 = 0;
    3806      129872 :       if (m < L-1) {
    3807       61211 :         GEN Q1 = qform_primeform2(L, D1), Q2, X;
    3808       61211 :         if (!Q1) pari_err_BUG("modpoly_pickD");
    3809       61211 :         Q2 = primeform_u(stoi(D2), L1);
    3810       61211 :         Q2 = gmul(Q1, Q2); /* we know this element has order L-1 */
    3811       61211 :         Q1 = primeform_u(stoi(D2), L0);
    3812       61211 :         k = ((n2 & 1) ? 2*n2 : n2)/(L-1);
    3813       61211 :         Q1 = gpowgs(Q1, k);
    3814       61211 :         X = qfi_Shanks(Q2, Q1, L-1);
    3815       61211 :         if (!X) {
    3816        8909 :           dbg_printf(3)("Bad D2=%ld for D1=%ld, D0=%ld, with n2=%ld, h2=%ld, L1=%ld, "
    3817             :               "form of norm L^2 not generated by L0 and L1\n",
    3818             :               D2, D1, D0, n2, h2, L1);
    3819        8909 :           continue;
    3820             :         }
    3821       52302 :         dl20 = itos(X) * k;
    3822       52302 :         dl21 = 1;
    3823             :       }
    3824      120963 :       if (! (m < L-1 || n2 < d*(L-1)) && n1 >= d && ! use_L1)
    3825       68295 :         L1 = 0;  /* we don't need L1 */
    3826             : 
    3827      120963 :       if (!L1 && use_L1) {
    3828           0 :         dbg_printf(3)("not using D2=%ld for D1=%ld, D0=%ld, with n2=%ld, h2=%ld, L1=%ld, "
    3829             :                    "because we don't need L1 but must use it\n",
    3830             :                    D2, D1, D0, n2, h2, L1);
    3831           0 :         continue;
    3832             :       }
    3833             :       /* don't allow zero dl21 with L1 for the moment, since
    3834             :        * modpoly doesn't handle it - we may change this in the future */
    3835      120963 :       if (L1 && ! dl21) continue;
    3836      120597 :       dbg_printf(3)("Good D0=%ld, D1=%ld, D2=%ld with s=%ld^%ld, L1=%ld, dl2=%ld, n2=%ld, h2=%ld\n",
    3837             :                  D0, D1, D2, p, j, L1, dl20, n2, h2);
    3838             : 
    3839             :       /* This estimate is heuristic and fiddling with the
    3840             :        * parameters 5 and 0.25 can change things quite a bit. */
    3841      120597 :       enum_cost = n2 * (5 * L0 * L0 + 0.25 * L1 * L1);
    3842      120597 :       cost = enum_cost + H_cost;
    3843      120597 :       if (best_cost && cost > 2.2*best_cost) break;
    3844       28428 :       if (best_cost && cost >= 0.99*best_cost) continue;
    3845             : 
    3846        7110 :       Dinfo.GENcode0 = evaltyp(t_VECSMALL)|evallg(13);
    3847        7110 :       Dinfo.inv = inv;
    3848        7110 :       Dinfo.L = L;
    3849        7110 :       Dinfo.D0 = D0;
    3850        7110 :       Dinfo.D1 = D1;
    3851        7110 :       Dinfo.L0 = L0;
    3852        7110 :       Dinfo.L1 = L1;
    3853        7110 :       Dinfo.n1 = n1;
    3854        7110 :       Dinfo.n2 = n2;
    3855        7110 :       Dinfo.dl1 = dl1;
    3856        7110 :       Dinfo.dl2_0 = dl20;
    3857        7110 :       Dinfo.dl2_1 = dl21;
    3858        7110 :       Dinfo.cost = cost;
    3859             : 
    3860        7110 :       if (!modpoly_pickD_primes(NULL, NULL, 0, NULL, 0, &Dinfo.bits, minbits, &Dinfo))
    3861          58 :         continue;
    3862        7052 :       dbg_printf(2)("Best D2=%ld, D1=%ld, D0=%ld with s=%ld^%ld, L1=%ld, "
    3863             :                  "n1=%ld, n2=%ld, cost ratio %.2f, bits=%ld\n",
    3864             :                  D2, D1, D0, p, j, L1, n1, n2,
    3865           0 :                  (double)cost/(d*(L-1)), Dinfo.bits);
    3866             :       /* Insert Dinfo into the Ds array.  Ds is sorted by ascending cost. */
    3867       35473 :       for (j = 0; j < Dcnt; j++)
    3868       32763 :         if (Dinfo.cost < Ds[j].cost) break;
    3869        7052 :       if (n2 > MAX_VOLCANO_FLOOR_SIZE && n2*(L1 ? 2 : 1) > 1.2* (d*(L-1)) ) {
    3870           0 :         dbg_printf(3)("Not using D1=%ld, D2=%ld for space reasons\n", D1, D2);
    3871           0 :         continue;
    3872             :       }
    3873        7052 :       if (j == Dcnt && Dcnt == MODPOLY_MAX_DCNT)
    3874           0 :         continue;
    3875        7052 :       totbits += Dinfo.bits;
    3876        7052 :       if (Dcnt == MODPOLY_MAX_DCNT) totbits -= Ds[Dcnt-1].bits;
    3877        7052 :       if (Dcnt < MODPOLY_MAX_DCNT) Dcnt++;
    3878        7052 :       if (n2 > MAX_VOLCANO_FLOOR_SIZE)
    3879           0 :         dbg_printf(3)("totbits=%ld, minbits=%ld\n", totbits, minbits);
    3880        7052 :       for (k = Dcnt-1; k > j; k--) Ds[k] = Ds[k-1];
    3881        7052 :       Ds[k] = Dinfo;
    3882        7052 :       best_cost = (totbits > minbits)? Ds[Dcnt-1].cost: 0;
    3883             :       /* if we were able to use D1 with s = 1, there is no point in
    3884             :        * using any larger D1 for the same D0 */
    3885        7052 :       if (!i) break;
    3886             :     } /* END FOR over small primes */
    3887             :   } /* END WHILE over D0's */
    3888        2701 :   dbg_printf(2)("  checked %ld of %ld fundamental discriminants to find suitable "
    3889             :                 "discriminant (Dcnt = %ld)\n", D0_i, tablen, Dcnt);
    3890        2701 :   if ( ! Dcnt) {
    3891           0 :     dbg_printf(1)("failed completely for L=%ld\n", L);
    3892           0 :     return 0;
    3893             :   }
    3894             : 
    3895        2701 :   Dcnt = calc_primes_for_discriminants(Ds, Dcnt, L, minbits);
    3896             : 
    3897             :   /* fill in any missing dl1's */
    3898        5418 :   for (i = 0 ; i < Dcnt; i++)
    3899        5389 :     if (Ds[i].dl1 < 0 &&
    3900        2672 :        (Ds[i].dl1 = primeform_discrete_log(L0, L, Ds[i].n1, Ds[i].D1)) < 0)
    3901           0 :         pari_err_BUG("modpoly_pickD");
    3902        2701 :   if (DEBUGLEVEL > 1+3) {
    3903           0 :     err_printf("Selected %ld discriminants using %ld msecs\n", Dcnt, timer_delay(&T));
    3904           0 :     for (i = 0 ; i < Dcnt ; i++)
    3905             :     {
    3906           0 :       GEN H = classno(stoi(Ds[i].D0));
    3907           0 :       long h0 = itos(H);
    3908           0 :       err_printf ("    D0=%ld, h(D0)=%ld, D=%ld, L0=%ld, L1=%ld, "
    3909             :           "cost ratio=%.2f, enum ratio=%.2f,",
    3910           0 :           Ds[i].D0, h0, Ds[i].D1, Ds[i].L0, Ds[i].L1,
    3911           0 :           (double)Ds[i].cost/(d*(L-1)),
    3912           0 :           (double)(Ds[i].n2*(Ds[i].L1 ? 2 : 1))/(d*(L-1)));
    3913           0 :       err_printf (" %ld primes, %ld bits\n", Ds[i].nprimes, Ds[i].bits);
    3914             :     }
    3915             :   }
    3916        2701 :   avma = ltop; return Dcnt;
    3917             : }
    3918             : 
    3919             : static int
    3920    12906347 : _qsort_cmp(const void *a, const void *b)
    3921             : {
    3922    12906347 :   D_entry *x = (D_entry *)a, *y = (D_entry *)b;
    3923             :   long u, v;
    3924             : 
    3925             :   /* u and v are the class numbers of x and y */
    3926    12906347 :   u = x->h * (!!(x->m & 2) + 1);
    3927    12906347 :   v = y->h * (!!(y->m & 2) + 1);
    3928             :   /* Sort by class number */
    3929    12906347 :   if (u < v) return -1;
    3930     8980635 :   if (u > v) return 1;
    3931             :   /* Sort by discriminant (which is < 0, hence the sign reversal) */
    3932     2694509 :   if (x->D > y->D) return -1;
    3933           0 :   if (x->D < y->D) return 1;
    3934           0 :   return 0;
    3935             : }
    3936             : 
    3937             : /* Build a table containing fundamental D, |D| <= maxD whose class groups
    3938             :  * - are cyclic generated by an element of norm L0
    3939             :  * - have class number at most maxh
    3940             :  * The table is ordered using _qsort_cmp above, which ranks the discriminants
    3941             :  * by class number, then by absolute discriminant.
    3942             :  *
    3943             :  * INPUT:
    3944             :  * - maxd: largest allowed discriminant
    3945             :  * - maxh: largest allowed class number
    3946             :  * - L0: norm of class group generator
    3947             :  *
    3948             :  * OUTPUT:
    3949             :  * - tablelen: length of return value
    3950             :  *
    3951             :  * RETURN:
    3952             :  * - array of {D, h(D), kronecker symbols for small p} */
    3953             : static D_entry *
    3954        2701 : scanD0(long *tablelen, long *minD, long maxD, long maxh, long L0)
    3955             : {
    3956             :   pari_sp av;
    3957             :   D_entry *tab;
    3958             :   long d, cnt;
    3959             : 
    3960             :   /* NB: As seen in the loop below, the real class number of D can be */
    3961             :   /* 2*maxh if cl(D) is cyclic. */
    3962        2701 :   if (maxh < 0) pari_err_BUG("scanD0");
    3963             : 
    3964             :   /* Not checked, but L0 should be 2, 3, 5 or 7. */
    3965        2701 :   tab = (D_entry *) stack_malloc((maxD/4)*sizeof(*tab)); /* Overestimate */
    3966             :   /* d = 7, 11, 15, 19, 23, ... */
    3967     6752500 :   for (av = avma, d = *minD, cnt = 0; d <= maxD; d += 4, avma = av)
    3968             :   {
    3969             :     GEN DD, H, fact, ordL, frm;
    3970     6749799 :     long i, j, k, n, h, L1, D = -d;
    3971             :     long *q, *e;
    3972             :     ulong m;
    3973             :     /* Check to see if (D | L0) = 1 */
    3974     6749799 :     if (kross(D, L0) < 1) continue;
    3975             : 
    3976             :     /* [q, e] is the factorisation of d. */
    3977     3054854 :     fact = factoru(d);
    3978     3054854 :     q = zv_to_longptr(gel(fact, 1));
    3979     3054854 :     e = zv_to_longptr(gel(fact, 2));
    3980     3054854 :     k = lg(gel(fact, 1)) - 1;
    3981             : 
    3982             :     /* Check if the discriminant is square-free */
    3983     7907366 :     for (i = 0; i < k; i++)
    3984     5368586 :       if (e[i] > 1) break;
    3985     3054854 :     if (i < k) continue;
    3986             : 
    3987             :     /* L1 initially the first factor of d if small enough, otherwise ignored */
    3988     2538780 :     L1 = (k > 1 && q[0] <= MAX_L1)? q[0]: 0;
    3989             : 
    3990             :     /* restrict to possibly cyclic class groups */
    3991     2538780 :     if (k > 2) continue;
    3992             : 
    3993             :     /* Check if h(D) is too big */
    3994     2054767 :     DD = stoi(D);
    3995     2054767 :     H = classno(DD);
    3996     2054767 :     h = itos(H);
    3997     2054767 :     if (h > 2*maxh || (!L1 && h > maxh)) continue;
    3998             : 
    3999             :     /* Check if ord(q) is not big enough to generate at least half the
    4000             :      * class group (where q is the L0-primeform). */
    4001     1928306 :     frm = primeform_u(DD, L0);
    4002     1928306 :     ordL = qfi_order(redimag(frm), H);
    4003     1928306 :     n = itos(ordL);
    4004     1928306 :     if (n < h/2 || (!L1 && n < h)) continue;
    4005             : 
    4006             :     /* If q is big enough, great!  Otherwise, for each potential L1,
    4007             :      * do a discrete log to see if it is NOT in the subgroup generated
    4008             :      * by L0; stop as soon as such is found. */
    4009     1870916 :     for (j = 0; ; j++) {
    4010     2087531 :       if (n == h || (L1 && !qfi_Shanks(primeform_u(DD, L1), frm, n))) {
    4011     1569168 :         dbg_printf(2)("D0=%ld good with L1=%ld\n", D, L1);
    4012     1569168 :         break;
    4013             :       }
    4014      301748 :       if (!L1) break;
    4015      216615 :       L1 = (j < k && k > 1 && q[j] <= MAX_L1 ? q[j] : 0);
    4016             :     }
    4017             :     /* The first bit of m indicates whether q generates a proper
    4018             :      * subgroup of cl(D) (hence implying that we need L1) or if q
    4019             :      * generates the whole class group. */
    4020     1654301 :     m = (n < h ? 1 : 0);
    4021             :     /* bits i and i+1 of m give the 2-bit number 1 + (D|p) where p is
    4022             :      * the ith prime. */
    4023    49043120 :     for (i = 1 ; i <= SMOOTH_PRIMES; i++)
    4024             :     {
    4025    47388819 :       ulong x  = (ulong) (1 + kross(D, PRIMES[i]));
    4026    47388819 :       m |= x << (2*i);
    4027             :     }
    4028             : 
    4029             :     /* Insert d, h and m into the table */
    4030     1654301 :     tab[cnt].D = D;
    4031     1654301 :     tab[cnt].h = h;
    4032     1654301 :     tab[cnt].m = m;
    4033     1654301 :     cnt++;
    4034             :   }
    4035             : 
    4036             :   /* Sort the table */
    4037        2701 :   qsort(tab, cnt, sizeof(*tab), _qsort_cmp);
    4038        2701 :   *tablelen = cnt; *minD = d; return tab;
    4039             : }
    4040             : 
    4041             : /* Populate Ds with discriminants (and attached data) that can be
    4042             :  * used to calculate the modular polynomial of level L and invariant
    4043             :  * inv.  Return the number of discriminants found. */
    4044             : static long
    4045        2699 : discriminant_with_classno_at_least(disc_info bestD[MODPOLY_MAX_DCNT],
    4046             :   long L, long inv, long ignore_sparse)
    4047             : {
    4048             :   enum { SMALL_L_BOUND = 101 };
    4049        2699 :   long max_max_D = 160000 * (inv ? 2 : 1);
    4050             :   long minD, maxD, maxh, L0, max_L1, minbits, Dcnt, flags, s, d, h, i, tablen;
    4051             :   D_entry *tab;
    4052        2699 :   double eps, cost, best_eps = -1.0, best_cost = -1.0;
    4053             :   disc_info Ds[MODPOLY_MAX_DCNT];
    4054        2699 :   long best_cnt = 0;
    4055             :   pari_timer T;
    4056        2699 :   timer_start(&T);
    4057             : 
    4058        2699 :   s = modinv_sparse_factor(inv);
    4059        2699 :   d = ceildivuu(L+1, s) + 1;
    4060             : 
    4061             :   /* maxD of 10000 allows us to get a satisfactory discriminant in
    4062             :    * under 250ms in most cases. */
    4063        2699 :   maxD = 10000;
    4064             :   /* Allow the class number to overshoot L by 50%.  Must be at least
    4065             :    * 1.1*L, and higher values don't seem to provide much benefit,
    4066             :    * except when L is small, in which case it's necessary to get any
    4067             :    * discriminant at all in some cases. */
    4068        2699 :   maxh = (L / s < SMALL_L_BOUND) ? 10 * L : 1.5 * L;
    4069             : 
    4070        2699 :   flags = ignore_sparse ? MODPOLY_IGNORE_SPARSE_FACTOR : 0;
    4071        2699 :   L0 = select_L0(L, inv, 0);
    4072        2699 :   max_L1 = L / 2 + 2;    /* for L=11 we need L1=7 for j */
    4073        2699 :   minbits = modpoly_height_bound(L, inv);
    4074        2699 :   minD = 7;
    4075             : 
    4076        8097 :   while ( ! best_cnt) {
    4077        5400 :     while (maxD <= max_max_D) {
    4078             :       /* TODO: Find a way to re-use tab when we need multiple modpolys */
    4079        2701 :       tab = scanD0(&tablen, &minD, maxD, maxh, L0);
    4080        2701 :       dbg_printf(1)("Found %ld potential fundamental discriminants\n", tablen);
    4081             : 
    4082        2701 :       Dcnt = modpoly_pickD(Ds, L, inv, L0, max_L1, minbits, flags, tab, tablen);
    4083        2701 :       eps = 0.0;
    4084        2701 :       cost = 0.0;
    4085             : 
    4086        2701 :       if (Dcnt) {
    4087        2699 :         long n1 = 0;
    4088        5416 :         for (i = 0; i < Dcnt; i++) {
    4089        2717 :           n1 = maxss(n1, Ds[i].n1);
    4090        2717 :           cost += Ds[i].cost;
    4091             :         }
    4092        2699 :         eps = (n1 * s - L) / (double)L;
    4093             : 
    4094        2699 :         if (best_cost < 0.0 || cost < best_cost) {
    4095        2699 :           if (best_cnt)
    4096           0 :             for (i = 0; i < best_cnt; i++) killbloc((GEN)bestD[i].primes);
    4097        2699 :           (void) memcpy(bestD, Ds, Dcnt * sizeof(disc_info));
    4098        2699 :           best_cost = cost;
    4099        2699 :           best_cnt = Dcnt;
    4100        2699 :           best_eps = eps;
    4101             :           /* We're satisfied if n1 is within 5% of L. */
    4102        2699 :           if (L / s <= SMALL_L_BOUND || eps < 0.05) break;
    4103             :         } else {
    4104           0 :           for (i = 0; i < Dcnt; i++) killbloc((GEN)Ds[i].primes);
    4105             :         }
    4106             :       } else {
    4107           2 :         if (log2(maxD) > BITS_IN_LONG - 2 * (log2(L) + 2))
    4108             :         {
    4109           0 :           char *err = stack_sprintf("modular polynomial of level %ld and invariant %ld",L,inv);
    4110           0 :           pari_err(e_ARCH, err);
    4111             :         }
    4112             :       }
    4113           2 :       maxD *= 2;
    4114           2 :       minD += 4;
    4115           2 :       dbg_printf(0)("  Doubling discriminant search space (closest: %.1f%%, cost ratio: %.1f)...\n", eps*100, cost/(double)(d*(L-1)));
    4116             :     }
    4117        2699 :     max_max_D *= 2;
    4118             :   }
    4119             : 
    4120        2699 :   if (DEBUGLEVEL > 3) {
    4121           0 :     pari_sp av = avma;
    4122           0 :     err_printf("Found discriminant(s):\n");
    4123           0 :     for (i = 0; i < best_cnt; ++i) {
    4124           0 :       h = itos(classno(stoi(bestD[i].D1)));
    4125           0 :       avma = av;
    4126           0 :       err_printf("  D = %ld, h = %ld, u = %ld, L0 = %ld, L1 = %ld, n1 = %ld, n2 = %ld, cost = %ld\n",
    4127           0 :           bestD[i].D1, h, usqrt(bestD[i].D1 / bestD[i].D0), bestD[i].L0,
    4128           0 :           bestD[i].L1, bestD[i].n1, bestD[i].n2, bestD[i].cost);
    4129             :     }
    4130           0 :     err_printf("(off target by %.1f%%, cost ratio: %.1f)\n",
    4131           0 :                best_eps*100, best_cost/(double)(d*(L-1)));
    4132             :   }
    4133        2699 :   return best_cnt;
    4134             : }

Generated by: LCOV version 1.13