Code coverage tests

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

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

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

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

Generated by: LCOV version 1.16