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

Generated by: LCOV version 1.11