"Fossies" - the Fresh Open Source Software Archive

Member "pari-2.13.1/src/basemath/mftrace.c" (14 Jan 2021, 367208 Bytes) of package /linux/misc/pari-2.13.1.tar.gz:


As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) C and C++ source code syntax highlighting (style: standard) with prefixed line numbers and code folding option. Alternatively you can here view or download the uninterpreted source code file. For more information about "mftrace.c" see the Fossies "Dox" file reference documentation and the latest Fossies "Diffs" side-by-side code changes report: 2.13.0_vs_2.13.1.

    1 /* Copyright (C) 2016  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 /*************************************************************************/
   16 /*                                                                       */
   17 /*              Modular forms package based on trace formulas            */
   18 /*                                                                       */
   19 /*************************************************************************/
   20 #include "pari.h"
   21 #include "paripriv.h"
   22 
   23 enum {
   24   MF_SPLIT = 1,
   25   MF_EISENSPACE,
   26   MF_FRICKE,
   27   MF_MF2INIT,
   28   MF_SPLITN
   29 };
   30 
   31 typedef struct {
   32   GEN vnew, vfull, DATA, VCHIP;
   33   long n, newHIT, newTOTAL, cuspHIT, cuspTOTAL;
   34 } cachenew_t;
   35 
   36 static void init_cachenew(cachenew_t *c, long n, long N, GEN f);
   37 static GEN mfinit_i(GEN NK, long space);
   38 static GEN mfinit_Nkchi(long N, long k, GEN CHI, long space, long flraw);
   39 static GEN mf2init_Nkchi(long N, long k, GEN CHI, long space, long flraw);
   40 static GEN mf2basis(long N, long r, GEN CHI, GEN *pCHI1, long space);
   41 static GEN mfeisensteinbasis(long N, long k, GEN CHI);
   42 static GEN mfeisensteindec(GEN mf, GEN F);
   43 static GEN initwt1newtrace(GEN mf);
   44 static GEN initwt1trace(GEN mf);
   45 static GEN myfactoru(long N);
   46 static GEN mydivisorsu(long N);
   47 static GEN Qab_Czeta(long k, long ord, GEN C, long vt);
   48 static GEN mfcoefs_i(GEN F, long n, long d);
   49 static GEN bhnmat_extend(GEN M, long m,long l, GEN S, cachenew_t *cache);
   50 static GEN initnewtrace(long N, GEN CHI);
   51 static void dbg_cachenew(cachenew_t *C);
   52 static GEN hecke_i(long m, long l, GEN V, GEN F, GEN DATA);
   53 static GEN c_Ek(long n, long d, GEN F);
   54 static GEN RgV_heckef2(long n, long d, GEN V, GEN F, GEN DATA);
   55 static GEN mfcusptrace_i(long N, long k, long n, GEN Dn, GEN TDATA);
   56 static GEN mfnewtracecache(long N, long k, long n, cachenew_t *cache);
   57 static GEN colnewtrace(long n0, long n, long d, long N, long k, cachenew_t *c);
   58 static GEN dihan(GEN bnr, GEN w, GEN k0j, ulong n);
   59 static GEN sigchi(long k, GEN CHI, long n);
   60 static GEN sigchi2(long k, GEN CHI1, GEN CHI2, long n, long ord);
   61 static GEN mflineardivtomat(long N, GEN vF, long n);
   62 static GEN mfdihedralcusp(long N, GEN CHI);
   63 static long mfdihedralcuspdim(long N, GEN CHI);
   64 static GEN mfdihedralnew(long N, GEN CHI);
   65 static GEN mfdihedralall(GEN LIM);
   66 static long mfwt1cuspdim(long N, GEN CHI);
   67 static long mf2dim_Nkchi(long N, long k, GEN CHI, ulong space);
   68 static long mfdim_Nkchi(long N, long k, GEN CHI, long space);
   69 static GEN charLFwtk(long N, long k, GEN CHI, long ord, long t);
   70 static GEN mfeisensteingacx(GEN E,long w,GEN ga,long n,long prec);
   71 static GEN mfgaexpansion(GEN mf, GEN F, GEN gamma, long n, long prec);
   72 static GEN mfEHmat(long n, long r);
   73 static GEN mfEHcoef(long r, long N);
   74 static GEN mftobasis_i(GEN mf, GEN F);
   75 
   76 static GEN
   77 mkgNK(GEN N, GEN k, GEN CHI, GEN P) { return mkvec4(N, k, CHI, P); }
   78 static GEN
   79 mkNK(long N, long k, GEN CHI) { return mkgNK(stoi(N), stoi(k), CHI, pol_x(1)); }
   80 GEN
   81 MF_get_CHI(GEN mf) { return gmael(mf,1,3); }
   82 GEN
   83 MF_get_gN(GEN mf) { return gmael(mf,1,1); }
   84 long
   85 MF_get_N(GEN mf) { return itou(MF_get_gN(mf)); }
   86 GEN
   87 MF_get_gk(GEN mf) { return gmael(mf,1,2); }
   88 long
   89 MF_get_k(GEN mf)
   90 {
   91   GEN gk = MF_get_gk(mf);
   92   if (typ(gk)!=t_INT) pari_err_IMPL("half-integral weight");
   93   return itou(gk);
   94 }
   95 long
   96 MF_get_r(GEN mf)
   97 {
   98   GEN gk = MF_get_gk(mf);
   99   if (typ(gk) == t_INT) pari_err_IMPL("integral weight");
  100   return itou(gel(gk, 1)) >> 1;
  101 }
  102 long
  103 MF_get_space(GEN mf) { return itos(gmael(mf,1,4)); }
  104 GEN
  105 MF_get_E(GEN mf) { return gel(mf,2); }
  106 GEN
  107 MF_get_S(GEN mf) { return gel(mf,3); }
  108 GEN
  109 MF_get_basis(GEN mf) { return shallowconcat(gel(mf,2), gel(mf,3)); }
  110 long
  111 MF_get_dim(GEN mf)
  112 {
  113   switch(MF_get_space(mf))
  114   {
  115     case mf_FULL:
  116       return lg(MF_get_S(mf)) - 1 + lg(MF_get_E(mf))-1;
  117     case mf_EISEN:
  118       return lg(MF_get_E(mf))-1;
  119     default: /* mf_NEW, mf_CUSP, mf_OLD */
  120       return lg(MF_get_S(mf)) - 1;
  121   }
  122 }
  123 GEN
  124 MFnew_get_vj(GEN mf) { return gel(mf,4); }
  125 GEN
  126 MFcusp_get_vMjd(GEN mf) { return gel(mf,4); }
  127 GEN
  128 MF_get_M(GEN mf) { return gmael(mf,5,3); }
  129 GEN
  130 MF_get_Minv(GEN mf) { return gmael(mf,5,2); }
  131 GEN
  132 MF_get_Mindex(GEN mf) { return gmael(mf,5,1); }
  133 
  134 /* ordinary gtocol forgets about initial 0s */
  135 GEN
  136 sertocol(GEN S) { return gtocol0(S, -(lg(S) - 2 + valp(S))); }
  137 /*******************************************************************/
  138 /*     Linear algebra in cyclotomic fields (TODO: export this)     */
  139 /*******************************************************************/
  140 /* return r and split prime p giving projection Q(zeta_n) -> Fp, zeta -> r */
  141 static ulong
  142 QabM_init(long n, ulong *p)
  143 {
  144   ulong pinit = 1000000007;
  145   forprime_t T;
  146   if (n <= 1) { *p = pinit; return 0; }
  147   u_forprime_arith_init(&T, pinit, ULONG_MAX, 1, n);
  148   *p = u_forprime_next(&T);
  149   return Flx_oneroot(ZX_to_Flx(polcyclo(n, 0), *p), *p);
  150 }
  151 static ulong
  152 Qab_to_Fl(GEN P, ulong r, ulong p)
  153 {
  154   ulong t;
  155   GEN den;
  156   P = Q_remove_denom(liftpol_shallow(P), &den);
  157   if (typ(P) == t_POL) { GEN Pp = ZX_to_Flx(P, p); t = Flx_eval(Pp, r, p); }
  158   else t = umodiu(P, p);
  159   if (den) t = Fl_div(t, umodiu(den, p), p);
  160   return t;
  161 }
  162 static GEN
  163 QabC_to_Flc(GEN C, ulong r, ulong p)
  164 {
  165   long i, l = lg(C);
  166   GEN A = cgetg(l, t_VECSMALL);
  167   for (i = 1; i < l; i++) uel(A,i) = Qab_to_Fl(gel(C,i), r, p);
  168   return A;
  169 }
  170 static GEN
  171 QabM_to_Flm(GEN M, ulong r, ulong p)
  172 {
  173   long i, l;
  174   GEN A = cgetg_copy(M, &l);
  175   for (i = 1; i < l; i++)
  176     gel(A, i) = QabC_to_Flc(gel(M, i), r, p);
  177   return A;
  178 }
  179 /* A a t_POL */
  180 static GEN
  181 QabX_to_Flx(GEN A, ulong r, ulong p)
  182 {
  183   long i, l = lg(A);
  184   GEN a = cgetg(l, t_VECSMALL);
  185   a[1] = ((ulong)A[1])&VARNBITS;
  186   for (i = 2; i < l; i++) uel(a,i) = Qab_to_Fl(gel(A,i), r, p);
  187   return Flx_renormalize(a, l);
  188 }
  189 
  190 /* FIXME: remove */
  191 static GEN
  192 ZabM_pseudoinv_i(GEN M, GEN P, long n, GEN *pv, GEN *den, int ratlift)
  193 {
  194   GEN v = ZabM_indexrank(M, P, n);
  195   if (pv) *pv = v;
  196   M = shallowmatextract(M,gel(v,1),gel(v,2));
  197   return ratlift? ZabM_inv_ratlift(M, P, n, den): ZabM_inv(M, P, n, den);
  198 }
  199 
  200 /* M matrix with coeff in Q(\chi)), where Q(\chi) = Q(X)/(P) for
  201  * P = cyclotomic Phi_n. Assume M rational if n <= 2 */
  202 static GEN
  203 QabM_ker(GEN M, GEN P, long n)
  204 {
  205   if (n <= 2) return QM_ker(M);
  206   return ZabM_ker(Q_primpart(liftpol_shallow(M)), P, n);
  207 }
  208 /* pseudo-inverse of M. FIXME: should replace QabM_pseudoinv */
  209 static GEN
  210 QabM_pseudoinv_i(GEN M, GEN P, long n, GEN *pv, GEN *pden)
  211 {
  212   GEN cM, Mi;
  213   if (n <= 2)
  214   {
  215     M = Q_primitive_part(M, &cM);
  216     Mi = ZM_pseudoinv(M, pv, pden); /* M^(-1) = Mi / (cM * den) */
  217   }
  218   else
  219   {
  220     M = Q_primitive_part(liftpol_shallow(M), &cM);
  221     Mi = ZabM_pseudoinv(M, P, n, pv, pden);
  222   }
  223   *pden = mul_content(*pden, cM);
  224   return Mi;
  225 }
  226 /* FIXME: delete */
  227 static GEN
  228 QabM_pseudoinv(GEN M, GEN P, long n, GEN *pv, GEN *pden)
  229 {
  230   GEN Mi = QabM_pseudoinv_i(M, P, n, pv, pden);
  231   return P? gmodulo(Mi, P): Mi;
  232 }
  233 
  234 static GEN
  235 QabM_indexrank(GEN M, GEN P, long n)
  236 {
  237   GEN z;
  238   if (n <= 2)
  239   {
  240     M = vec_Q_primpart(M);
  241     z = ZM_indexrank(M); /* M^(-1) = Mi / (cM * den) */
  242   }
  243   else
  244   {
  245     M = vec_Q_primpart(liftpol_shallow(M));
  246     z = ZabM_indexrank(M, P, n);
  247   }
  248   return z;
  249 }
  250 
  251 /*********************************************************************/
  252 /*                    Simple arithmetic functions                    */
  253 /*********************************************************************/
  254 /* TODO: most of these should be exported and used in ifactor1.c */
  255 /* phi(n) */
  256 static ulong
  257 myeulerphiu(ulong n)
  258 {
  259   pari_sp av;
  260   if (n == 1) return 1;
  261   av = avma; return gc_ulong(av, eulerphiu_fact(myfactoru(n)));
  262 }
  263 static long
  264 mymoebiusu(ulong n)
  265 {
  266   pari_sp av;
  267   if (n == 1) return 1;
  268   av = avma; return gc_long(av, moebiusu_fact(myfactoru(n)));
  269 }
  270 
  271 static long
  272 mynumdivu(long N)
  273 {
  274   pari_sp av;
  275   if (N == 1) return 1;
  276   av = avma; return gc_long(av, numdivu_fact(myfactoru(N)));
  277 }
  278 
  279 /* N\prod_{p|N} (1+1/p) */
  280 static long
  281 mypsiu(ulong N)
  282 {
  283   pari_sp av = avma;
  284   GEN P = gel(myfactoru(N), 1);
  285   long j, l = lg(P), res = N;
  286   for (j = 1; j < l; j++) res += res/P[j];
  287   return gc_long(av,res);
  288 }
  289 /* write n = mf^2. Return m, set f. */
  290 static ulong
  291 mycore(ulong n, long *pf)
  292 {
  293   pari_sp av = avma;
  294   GEN fa = myfactoru(n), P = gel(fa,1), E = gel(fa,2);
  295   long i, l = lg(P), m = 1, f = 1;
  296   for (i = 1; i < l; i++)
  297   {
  298     long j, p = P[i], e = E[i];
  299     if (e & 1) m *= p;
  300     for (j = 2; j <= e; j+=2) f *= p;
  301   }
  302   *pf = f; return gc_long(av,m);
  303 }
  304 
  305 /* fa = factorization of -D > 0, return -D0 > 0 (where D0 is fundamental) */
  306 static long
  307 corediscs_fact(GEN fa)
  308 {
  309   GEN P = gel(fa,1), E = gel(fa,2);
  310   long i, l = lg(P), m = 1;
  311   for (i = 1; i < l; i++)
  312   {
  313     long p = P[i], e = E[i];
  314     if (e & 1) m *= p;
  315   }
  316   if ((m&3L) != 3) m <<= 2;
  317   return m;
  318 }
  319 static long
  320 mubeta(long n)
  321 {
  322   pari_sp av = avma;
  323   GEN E = gel(myfactoru(n), 2);
  324   long i, s = 1, l = lg(E);
  325   for (i = 1; i < l; i++)
  326   {
  327     long e = E[i];
  328     if (e >= 3) return gc_long(av,0);
  329     if (e == 1) s *= -2;
  330   }
  331   return gc_long(av,s);
  332 }
  333 
  334 /* n = n1*n2, n1 = ppo(n, m); return mubeta(n1)*moebiusu(n2).
  335  * N.B. If n from newt_params we, in fact, never return 0 */
  336 static long
  337 mubeta2(long n, long m)
  338 {
  339   pari_sp av = avma;
  340   GEN fa = myfactoru(n), P = gel(fa,1), E = gel(fa,2);
  341   long i, s = 1, l = lg(P);
  342   for (i = 1; i < l; i++)
  343   {
  344     long p = P[i], e = E[i];
  345     if (m % p)
  346     { /* p^e in n1 */
  347       if (e >= 3) return gc_long(av,0);
  348       if (e == 1) s *= -2;
  349     }
  350     else
  351     { /* in n2 */
  352       if (e >= 2) return gc_long(av,0);
  353       s = -s;
  354     }
  355   }
  356   return gc_long(av,s);
  357 }
  358 
  359 /* write N = prod p^{ep} and n = df^2, d squarefree.
  360  * set g  = ppo(gcd(sqfpart(N), f), FC)
  361  *     N2 = prod p^if(e==1 || p|n, ep-1, ep-2) */
  362 static void
  363 newt_params(long N, long n, long FC, long *pg, long *pN2)
  364 {
  365   GEN fa = myfactoru(N), P = gel(fa,1), E = gel(fa,2);
  366   long i, g = 1, N2 = 1, l = lg(P);
  367   for (i = 1; i < l; i++)
  368   {
  369     long p = P[i], e = E[i];
  370     if (e == 1)
  371     { if (FC % p && n % (p*p) == 0) g *= p; }
  372     else
  373       N2 *= upowuu(p,(n % p)? e-2: e-1);
  374   }
  375   *pg = g; *pN2 = N2;
  376 }
  377 /* simplified version of newt_params for n = 1 (newdim) */
  378 static void
  379 newd_params(long N, long *pN2)
  380 {
  381   GEN fa = myfactoru(N), P = gel(fa,1), E = gel(fa,2);
  382   long i, N2 = 1, l = lg(P);
  383   for (i = 1; i < l; i++)
  384   {
  385     long p = P[i], e = E[i];
  386     if (e > 2) N2 *= upowuu(p, e-2);
  387   }
  388   *pN2 = N2;
  389 }
  390 
  391 static long
  392 newd_params2(long N)
  393 {
  394   GEN fa = myfactoru(N), P = gel(fa,1), E = gel(fa,2);
  395   long i, N2 = 1, l = lg(P);
  396   for (i = 1; i < l; i++)
  397   {
  398     long p = P[i], e = E[i];
  399     if (e >= 2) N2 *= upowuu(p, e);
  400   }
  401   return N2;
  402 }
  403 
  404 /*******************************************************************/
  405 /*   Relative trace between cyclotomic fields (TODO: export this)  */
  406 /*******************************************************************/
  407 /* g>=1; return g * prod_{p | g, (p,q) = 1} (1-1/p) */
  408 static long
  409 phipart(long g, long q)
  410 {
  411   if (g > 1)
  412   {
  413     GEN P = gel(myfactoru(g), 1);
  414     long i, l = lg(P);
  415     for (i = 1; i < l; i++) { long p = P[i]; if (q % p) g -= g / p; }
  416   }
  417   return g;
  418 }
  419 /* Set s,v s.t. Trace(zeta_N^k) from Q(zeta_N) to Q(\zeta_N) = s * zeta_M^v
  420  * With k > 0, N = M*d and N, M != 2 mod 4 */
  421 static long
  422 tracerelz(long *pv, long d, long M, long k)
  423 {
  424   long s, g, q, muq;
  425   if (d == 1) { *pv = k; return 1; }
  426   *pv = 0; g = ugcd(k, d); q = d / g;
  427   muq = mymoebiusu(q); if (!muq) return 0;
  428   if (M != 1)
  429   {
  430     long v = Fl_invsafe(q % M, M);
  431     if (!v) return 0;
  432     *pv = (v * (k/g)) % M;
  433   }
  434   s = phipart(g, M*q); if (muq < 0) s = -s;
  435   return s;
  436 }
  437 /* Pi = polcyclo(i), i = m or n. Let Ki = Q(zeta_i), initialize Tr_{Kn/Km} */
  438 GEN
  439 Qab_trace_init(long n, long m, GEN Pn, GEN Pm)
  440 {
  441   long a, i, j, N, M, vt, d, D;
  442   GEN T, G;
  443 
  444   if (m == n || n <= 2) return mkvec(Pm);
  445   vt = varn(Pn);
  446   d = degpol(Pn);
  447   /* if (N != n) zeta_N = zeta_n^2 and zeta_n = - zeta_N^{(N+1)/2} */
  448   N = ((n & 3) == 2)? n >> 1: n;
  449   M = ((m & 3) == 2)? m >> 1: m; /* M | N | n */
  450   a = N / M;
  451   T = const_vec(d, NULL);
  452   D = d / degpol(Pm); /* relative degree */
  453   if (D == 1) G = NULL;
  454   else
  455   { /* zeta_M = zeta_n^A; s_j(zeta_M) = zeta_M <=> j = 1 (mod J) */
  456     long lG, A = (N == n)? a: (a << 1), J = n / ugcd(n, A);
  457     G = coprimes_zv(n);
  458     for (j = lG = 1; j < n; j += J)
  459       if (G[j]) G[lG++] = j;
  460     setlg(G, lG); /* Gal(Q(zeta_n) / Q(zeta_m)) */
  461   }
  462   T = const_vec(d, NULL);
  463   gel(T,1) = utoipos(D); /* Tr 1 */
  464   for (i = 1; i < d; i++)
  465   { /* if n = 2N, zeta_n^i = (-1)^i zeta_N^k */
  466     long s, v, k;
  467     GEN t;
  468 
  469     if (gel(T, i+1)) continue;
  470     k = (N == n)? i: ((odd(i)? i + N: i) >> 1);
  471     if ((s = tracerelz(&v, a, M, k)))
  472     {
  473       if (m != M) v *= 2;/* Tr = s * zeta_m^v */
  474       if (n != N && odd(i)) s = -s;
  475       t = Qab_Czeta(v, m, stoi(s), vt);
  476     }
  477     else
  478       t = gen_0;
  479     /* t = Tr_{Kn/Km} zeta_n^i; fill using Galois action */
  480     if (!G)
  481       gel(T, i + 1) = t;
  482     else
  483       for (j = 1; j <= D; j++)
  484       {
  485         long z = Fl_mul(i,G[j], n);
  486         if (z < d) gel(T, z + 1) = t;
  487       }
  488   }
  489   return mkvec3(Pm, Pn, T);
  490 }
  491 /* x a t_POL modulo Phi_n */
  492 static GEN
  493 tracerel_i(GEN T, GEN x)
  494 {
  495   long k, l = lg(x);
  496   GEN S;
  497   if (l == 2) return gen_0;
  498   S = gmul(gel(T,1), gel(x,2));
  499   for (k = 3; k < l; k++) S = gadd(S, gmul(gel(T,k-1), gel(x,k)));
  500   return S;
  501 }
  502 static GEN
  503 tracerel(GEN a, GEN v, GEN z)
  504 {
  505   a = liftpol_shallow(a);
  506   a = simplify_shallow(z? gmul(z,a): a);
  507   if (typ(a) == t_POL)
  508   {
  509     GEN T = gel(v,3);
  510     long degrel = itou(gel(T,1));
  511     a = tracerel_i(T, RgX_rem(a, gel(v,2)));
  512     if (degrel != 1) a = gdivgs(a, degrel);
  513     if (typ(a) == t_POL) a = RgX_rem(a, gel(v,1));
  514   }
  515   return a;
  516 }
  517 static GEN
  518 tracerel_z(GEN v, long t)
  519 {
  520   GEN Pn = gel(v,2);
  521   return t? pol_xn(t, varn(Pn)): NULL;
  522 }
  523 /* v = Qab_trace_init(n,m); x is a t_VEC of polmodulo Phi_n; Kn = Q(zeta_n)
  524  * [Kn:Km]^(-1) Tr_{Kn/Km} (zeta_n^t * x); 0 <= t < [Kn:Km] */
  525 GEN
  526 Qab_tracerel(GEN v, long t, GEN a)
  527 {
  528   if (lg(v) != 4) return a; /* => t = 0 */
  529   return tracerel(a, v, tracerel_z(v, t));
  530 }
  531 GEN
  532 QabV_tracerel(GEN v, long t, GEN x)
  533 {
  534   GEN z;
  535   if (lg(v) != 4) return x; /* => t = 0 */
  536   z = tracerel_z(v, t);
  537   pari_APPLY_same(tracerel(gel(x,i), v, z));
  538 }
  539 GEN
  540 QabM_tracerel(GEN v, long t, GEN x)
  541 {
  542   if (lg(v) != 4) return x;
  543   pari_APPLY_same(QabV_tracerel(v, t, gel(x,i)));
  544 }
  545 
  546 /* C*zeta_o^k mod X^o - 1 */
  547 static GEN
  548 Qab_Czeta(long k, long o, GEN C, long vt)
  549 {
  550   if (!k) return C;
  551   if (!odd(o))
  552   { /* optimization: reduce max degree by a factor 2 for free */
  553     o >>= 1;
  554     if (k >= o) { k -= o; C = gneg(C); if (!k) return C; }
  555   }
  556   return monomial(C, k, vt);
  557 }
  558 /* zeta_o^k */
  559 static GEN
  560 Qab_zeta(long k, long o, long vt) { return Qab_Czeta(k, o, gen_1, vt); }
  561 
  562 /*              Operations on Dirichlet characters                       */
  563 
  564 /* A Dirichlet character can be given in GP in different formats, but in this
  565  * package, it will be a vector CHI=[G,chi,ord], where G is the (Z/MZ)^* to
  566  * which the character belongs, chi is the character in Conrey format, ord is
  567  * the order */
  568 
  569 static GEN
  570 gmfcharorder(GEN CHI) { return gel(CHI, 3); }
  571 long
  572 mfcharorder(GEN CHI) { return itou(gmfcharorder(CHI)); }
  573 static long
  574 mfcharistrivial(GEN CHI) { return !CHI || mfcharorder(CHI) == 1; }
  575 static GEN
  576 gmfcharmodulus(GEN CHI) { return gmael3(CHI, 1, 1, 1); }
  577 long
  578 mfcharmodulus(GEN CHI) { return itou(gmfcharmodulus(CHI)); }
  579 GEN
  580 mfcharpol(GEN CHI) { return gel(CHI,4); }
  581 
  582 /* vz[i+1] = image of (zeta_o)^i in Fp */
  583 static ulong
  584 Qab_Czeta_Fl(long k, GEN vz, ulong C, ulong p)
  585 {
  586   long o;
  587   if (!k) return C;
  588   o = lg(vz)-2;
  589   if ((k << 1) == o) return Fl_neg(C,p);
  590   return Fl_mul(C, vz[k+1], p);
  591 }
  592 
  593 static long
  594 znchareval_i(GEN CHI, long n, GEN ord)
  595 { return itos(znchareval(gel(CHI,1), gel(CHI,2), stoi(n), ord)); }
  596 
  597 /* n coprime with the modulus of CHI */
  598 static GEN
  599 mfchareval(GEN CHI, long n)
  600 {
  601   GEN Pn, C, go = gmfcharorder(CHI);
  602   long k, o = go[2];
  603   if (o == 1) return gen_1;
  604   k = znchareval_i(CHI, n, go);
  605   Pn = mfcharpol(CHI);
  606   C = Qab_zeta(k, o, varn(Pn));
  607   if (typ(C) != t_POL) return C;
  608   return gmodulo(C, Pn);
  609 }
  610 /* d a multiple of ord(CHI); n coprime with char modulus;
  611  * return x s.t. CHI(n) = \zeta_d^x] */
  612 static long
  613 mfcharevalord(GEN CHI, long n, long d)
  614 {
  615   if (mfcharorder(CHI) == 1) return 0;
  616   return znchareval_i(CHI, n, utoi(d));
  617 }
  618 
  619 /* G a znstar, L a Conrey log: return a 'mfchar' */
  620 static GEN
  621 mfcharGL(GEN G, GEN L)
  622 {
  623   GEN o = zncharorder(G,L);
  624   long ord = itou(o), vt = fetch_user_var("t");
  625   return mkvec4(G, L, o, polcyclo(ord,vt));
  626 }
  627 static GEN
  628 mfchartrivial()
  629 { return mfcharGL(znstar0(gen_1,1), cgetg(1,t_COL)); }
  630 /* convert a generic character into an 'mfchar' */
  631 static GEN
  632 get_mfchar(GEN CHI)
  633 {
  634   GEN G, L;
  635   if (typ(CHI) != t_VEC) CHI = znchar(CHI);
  636   else
  637   {
  638     long l = lg(CHI);
  639     if ((l != 3 && l != 5) || !checkznstar_i(gel(CHI,1)))
  640       pari_err_TYPE("checkNF [chi]", CHI);
  641     if (l == 5) return CHI;
  642   }
  643   G = gel(CHI,1);
  644   L = gel(CHI,2); if (typ(L) != t_COL) L = znconreylog(G,L);
  645   return mfcharGL(G, L);
  646 }
  647 
  648 /* parse [N], [N,k], [N,k,CHI]. If 'joker' is set, allow wildcard for CHI */
  649 static GEN
  650 checkCHI(GEN NK, long N, int joker)
  651 {
  652   GEN CHI;
  653   if (lg(NK) == 3)
  654     CHI = mfchartrivial();
  655   else
  656   {
  657     long i, l;
  658     CHI = gel(NK,3); l = lg(CHI);
  659     if (isintzero(CHI) && joker)
  660       CHI = NULL; /* all character orbits */
  661     else if (isintm1(CHI) && joker > 1)
  662       CHI = gen_m1; /* sum over all character orbits */
  663     else if ((typ(CHI) == t_VEC &&
  664              (l == 1 || l != 3 || !checkznstar_i(gel(CHI,1)))) && joker)
  665     {
  666       CHI = shallowtrans(CHI); /* list of characters */
  667       for (i = 1; i < l; i++) gel(CHI,i) = get_mfchar(gel(CHI,i));
  668     }
  669     else
  670     {
  671       CHI = get_mfchar(CHI); /* single char */
  672       if (N % mfcharmodulus(CHI)) pari_err_TYPE("checkNF [chi]", NK);
  673     }
  674   }
  675   return CHI;
  676 }
  677 /* support half-integral weight */
  678 static void
  679 checkNK2(GEN NK, long *N, long *nk, long *dk, GEN *CHI, int joker)
  680 {
  681   long l = lg(NK);
  682   GEN T;
  683   if (typ(NK) != t_VEC || l < 3 || l > 4) pari_err_TYPE("checkNK", NK);
  684   T = gel(NK,1); if (typ(T) != t_INT) pari_err_TYPE("checkNF [N]", NK);
  685   *N = itos(T); if (*N <= 0) pari_err_TYPE("checkNF [N <= 0]", NK);
  686   T = gel(NK,2);
  687   switch(typ(T))
  688   {
  689     case t_INT:  *nk = itos(T); *dk = 1; break;
  690     case t_FRAC:
  691       *nk = itos(gel(T,1));
  692       *dk = itou(gel(T,2)); if (*dk == 2) break;
  693     default: pari_err_TYPE("checkNF [k]", NK);
  694   }
  695   *CHI = checkCHI(NK, *N, joker);
  696 }
  697 /* don't support half-integral weight */
  698 static void
  699 checkNK(GEN NK, long *N, long *k, GEN *CHI, int joker)
  700 {
  701   long d;
  702   checkNK2(NK, N, k, &d, CHI, joker);
  703   if (d != 1) pari_err_TYPE("checkNF [k]", NK);
  704 }
  705 
  706 static GEN
  707 mfchargalois(long N, int odd, GEN flagorder)
  708 {
  709   GEN G = znstar0(utoi(N), 1), L = chargalois(G, flagorder);
  710   long l = lg(L), i, j;
  711   for (i = j = 1; i < l; i++)
  712   {
  713     GEN chi = znconreyfromchar(G, gel(L,i));
  714     if (zncharisodd(G,chi) == odd) gel(L,j++) = mfcharGL(G,chi);
  715   }
  716   setlg(L, j); return L;
  717 }
  718 /* possible characters for nontrivial S_1(N, chi) */
  719 static GEN
  720 mfwt1chars(long N, GEN vCHI)
  721 {
  722   if (vCHI) return vCHI; /*do not filter, user knows best*/
  723   /* Tate's theorem */
  724   return mfchargalois(N, 1, uisprime(N)? mkvecsmall2(2,4): NULL);
  725 }
  726 static GEN
  727 mfchars(long N, long k, long dk, GEN vCHI)
  728 { return vCHI? vCHI: mfchargalois(N, (dk == 2)? 0: (k & 1), NULL); }
  729 
  730 /* wrappers from mfchar to znchar */
  731 static long
  732 mfcharparity(GEN CHI)
  733 {
  734   if (!CHI) return 1;
  735   return zncharisodd(gel(CHI,1), gel(CHI,2)) ? -1 : 1;
  736 }
  737 /* if CHI is primitive, return CHI itself, not a copy */
  738 static GEN
  739 mfchartoprimitive(GEN CHI, long *pF)
  740 {
  741   pari_sp av;
  742   GEN chi, F;
  743   if (!CHI) { if (pF) *pF = 1; return mfchartrivial(); }
  744   av = avma; F = znconreyconductor(gel(CHI,1), gel(CHI,2), &chi);
  745   if (typ(F) == t_INT) set_avma(av);
  746   else
  747   {
  748     CHI = leafcopy(CHI);
  749     gel(CHI,1) = znstar0(F, 1);
  750     gel(CHI,2) = chi;
  751   }
  752   if (pF) *pF = mfcharmodulus(CHI);
  753   return CHI;
  754 }
  755 static long
  756 mfcharconductor(GEN CHI)
  757 {
  758   pari_sp av = avma;
  759   GEN res = znconreyconductor(gel(CHI,1), gel(CHI,2), NULL);
  760   if (typ(res) == t_VEC) res = gel(res, 1);
  761   return gc_long(av, itos(res));
  762 }
  763 
  764 /*                      Operations on mf closures                    */
  765 static GEN
  766 tagparams(long t, GEN NK) { return mkvec2(mkvecsmall(t), NK); }
  767 static GEN
  768 lfuntag(long t, GEN x) { return mkvec2(mkvecsmall(t), x); }
  769 static GEN
  770 tag0(long t, GEN NK) { retmkvec(tagparams(t,NK)); }
  771 static GEN
  772 tag(long t, GEN NK, GEN x) { retmkvec2(tagparams(t,NK), x); }
  773 static GEN
  774 tag2(long t, GEN NK, GEN x, GEN y) { retmkvec3(tagparams(t,NK), x,y); }
  775 static GEN
  776 tag3(long t, GEN NK, GEN x,GEN y,GEN z) { retmkvec4(tagparams(t,NK), x,y,z); }
  777 static GEN
  778 tag4(long t, GEN NK, GEN x,GEN y,GEN z,GEN a)
  779 { retmkvec5(tagparams(t,NK), x,y,z,a); }
  780 /* is F a "modular form" ? */
  781 int
  782 checkmf_i(GEN F)
  783 { return typ(F) == t_VEC
  784     && lg(F) > 1 && typ(gel(F,1)) == t_VEC
  785     && lg(gel(F,1)) == 3
  786     && typ(gmael(F,1,1)) == t_VECSMALL
  787     && typ(gmael(F,1,2)) == t_VEC; }
  788 long mf_get_type(GEN F) { return gmael(F,1,1)[1]; }
  789 GEN mf_get_gN(GEN F) { return gmael3(F,1,2,1); }
  790 GEN mf_get_gk(GEN F) { return gmael3(F,1,2,2); }
  791 /* k - 1/2, assume k in 1/2 + Z */
  792 long mf_get_r(GEN F) { return itou(gel(mf_get_gk(F),1)) >> 1; }
  793 long mf_get_N(GEN F) { return itou(mf_get_gN(F)); }
  794 long mf_get_k(GEN F)
  795 {
  796   GEN gk = mf_get_gk(F);
  797   if (typ(gk)!=t_INT) pari_err_IMPL("half-integral weight");
  798   return itou(gk);
  799 }
  800 GEN mf_get_CHI(GEN F) { return gmael3(F,1,2,3); }
  801 GEN mf_get_field(GEN F) { return gmael3(F,1,2,4); }
  802 GEN mf_get_NK(GEN F) { return gmael(F,1,2); }
  803 static void
  804 mf_setfield(GEN f, GEN P)
  805 {
  806   gel(f,1) = leafcopy(gel(f,1));
  807   gmael(f,1,2) = leafcopy(gmael(f,1,2));
  808   gmael3(f,1,2,4) = P;
  809 }
  810 
  811 /* UTILITY FUNCTIONS */
  812 GEN
  813 mftocol(GEN F, long lim, long d)
  814 { GEN c = mfcoefs_i(F, lim, d); settyp(c,t_COL); return c; }
  815 GEN
  816 mfvectomat(GEN vF, long lim, long d)
  817 {
  818   long j, l = lg(vF);
  819   GEN M = cgetg(l, t_MAT);
  820   for (j = 1; j < l; j++) gel(M,j) = mftocol(gel(vF,j), lim, d);
  821   return M;
  822 }
  823 
  824 static GEN
  825 RgV_to_ser_full(GEN x) { return RgV_to_ser(x, 0, lg(x)+1); }
  826 /* TODO: delete */
  827 static GEN
  828 mfcoefsser(GEN F, long n) { return RgV_to_ser_full(mfcoefs_i(F,n,1)); }
  829 static GEN
  830 sertovecslice(GEN S, long n)
  831 {
  832   GEN v = gtovec0(S, -(lg(S) - 2 + valp(S)));
  833   long l = lg(v), n2 = n + 2;
  834   if (l < n2) pari_err_BUG("sertovecslice [n too large]");
  835   return (l == n2)? v: vecslice(v, 1, n2-1);
  836 }
  837 
  838 /* a, b two RgV of the same length, multiply as truncated power series */
  839 static GEN
  840 RgV_mul_RgXn(GEN a, GEN b)
  841 {
  842   long n = lg(a)-1;
  843   GEN c;
  844   a = RgV_to_RgX(a,0);
  845   b = RgV_to_RgX(b,0); c = RgXn_mul(a, b, n);
  846   c = RgX_to_RgC(c,n); settyp(c,t_VEC); return c;
  847 }
  848 /* divide as truncated power series */
  849 static GEN
  850 RgV_div_RgXn(GEN a, GEN b)
  851 {
  852   long n = lg(a)-1;
  853   GEN c;
  854   a = RgV_to_RgX(a,0);
  855   b = RgV_to_RgX(b,0); c = RgXn_mul(a, RgXn_inv(b,n), n);
  856   c = RgX_to_RgC(c,n); settyp(c,t_VEC); return c;
  857 }
  858 /* a^b */
  859 static GEN
  860 RgV_pows_RgXn(GEN a, long b)
  861 {
  862   long n = lg(a)-1;
  863   GEN c;
  864   a = RgV_to_RgX(a,0);
  865   if (b < 0) { a = RgXn_inv(a, n); b = -b; }
  866   c = RgXn_powu_i(a,b,n);
  867   c = RgX_to_RgC(c,n); settyp(c,t_VEC); return c;
  868 }
  869 
  870 /* assume lg(V) >= n*d + 2 */
  871 static GEN
  872 c_deflate(long n, long d, GEN v)
  873 {
  874   long i, id, l = n+2;
  875   GEN w;
  876   if (d == 1) return lg(v) == l ? v: vecslice(v, 1, l-1);
  877   w = cgetg(l, typ(v));
  878   for (i = id = 1; i < l; i++, id += d) gel(w, i) = gel(v, id);
  879   return w;
  880 }
  881 
  882 static void
  883 err_cyclo(void)
  884 { pari_err_IMPL("changing cyclotomic fields in mf"); }
  885 /* Q(zeta_a) = Q(zeta_b) ? */
  886 static int
  887 same_cyc(long a, long b)
  888 { return (a == b) || (odd(a) && b == (a<<1)) || (odd(b) && a == (b<<1)); }
  889 /* need to combine elements in Q(CHI1) and Q(CHI2) with result in Q(CHI),
  890  * CHI = CHI1 * CHI2 or CHI / CHI2 times some character of order 2 */
  891 static GEN
  892 chicompat(GEN CHI, GEN CHI1, GEN CHI2)
  893 {
  894   long o1 = mfcharorder(CHI1);
  895   long o2 = mfcharorder(CHI2), O, o;
  896   GEN T1, T2, P, Po;
  897   if (o1 <= 2 && o2 <= 2) return NULL;
  898   o = mfcharorder(CHI);
  899   Po = mfcharpol(CHI);
  900   P = mfcharpol(CHI1);
  901   if (o1 == o2)
  902   {
  903     if (o1 == o) return NULL;
  904     if (!same_cyc(o1,o)) err_cyclo();
  905     return mkvec4(P, gen_1,gen_1, Qab_trace_init(o1, o, P, Po));
  906   }
  907   O = ulcm(o1, o2);
  908   if (!same_cyc(O,o)) err_cyclo();
  909   if (O != o1) P = (O == o2)? mfcharpol(CHI2): polcyclo(O, varn(P));
  910   T1 = o1 <= 2? gen_1: utoipos(O / o1);
  911   T2 = o2 <= 2? gen_1: utoipos(O / o2);
  912   return mkvec4(P, T1, T2, O == o? gen_1: Qab_trace_init(O, o, P, Po));
  913 }
  914 /* *F a vector of cyclotomic numbers */
  915 static void
  916 compatlift(GEN *F, long o, GEN P)
  917 {
  918   long i, l;
  919   GEN f = *F, g = cgetg_copy(f,&l);
  920   for (i = 1; i < l; i++)
  921   {
  922     GEN fi = lift_shallow(gel(f,i));
  923     gel(g,i) = gmodulo(typ(fi)==t_POL? RgX_inflate(fi,o): fi, P);
  924   }
  925   *F = g;
  926 }
  927 static void
  928 chicompatlift(GEN T, GEN *F, GEN *G)
  929 {
  930   long o1 = itou(gel(T,2)), o2 = itou(gel(T,3));
  931   GEN P = gel(T,1);
  932   if (o1 != 1) compatlift(F, o1, P);
  933   if (o2 != 1 && G) compatlift(G, o2, P);
  934 }
  935 static GEN
  936 chicompatfix(GEN T, GEN F)
  937 {
  938   GEN V = gel(T,4);
  939   if (typ(V) == t_VEC) F = gmodulo(QabV_tracerel(V, 0, F), gel(V,1));
  940   return F;
  941 }
  942 
  943 static GEN
  944 c_mul(long n, long d, GEN S)
  945 {
  946   pari_sp av = avma;
  947   long nd = n*d;
  948   GEN F = gel(S,2), G = gel(S,3);
  949   F = mfcoefs_i(F, nd, 1);
  950   G = mfcoefs_i(G, nd, 1);
  951   if (lg(S) == 5) chicompatlift(gel(S,4),&F,&G);
  952   F = c_deflate(n, d, RgV_mul_RgXn(F,G));
  953   if (lg(S) == 5) F = chicompatfix(gel(S,4), F);
  954   return gerepilecopy(av, F);
  955 }
  956 static GEN
  957 c_pow(long n, long d, GEN S)
  958 {
  959   pari_sp av = avma;
  960   long nd = n*d;
  961   GEN F = gel(S,2), a = gel(S,3), f = mfcoefs_i(F,nd,1);
  962   if (lg(S) == 5) chicompatlift(gel(S,4),&F, NULL);
  963   f = RgV_pows_RgXn(f, itos(a));
  964   f = c_deflate(n, d, f);
  965   if (lg(S) == 5) f = chicompatfix(gel(S,4), f);
  966   return gerepilecopy(av, f);
  967 }
  968 
  969 /* F * Theta */
  970 static GEN
  971 mfmultheta(GEN F)
  972 {
  973   if (typ(mf_get_gk(F)) == t_FRAC && mf_get_type(F) == t_MF_DIV)
  974   {
  975     GEN T = gel(F,3); /* hopefully mfTheta() */
  976     if (mf_get_type(T) == t_MF_THETA && mf_get_N(T) == 4) return gel(F,2);
  977   }
  978   return mfmul(F, mfTheta(NULL));
  979 }
  980 
  981 static GEN
  982 c_bracket(long n, long d, GEN S)
  983 {
  984   pari_sp av = avma;
  985   long i, nd = n*d;
  986   GEN F = gel(S,2), G = gel(S,3), tF, tG, C, mpow, res, gk, gl;
  987   GEN VF = mfcoefs_i(F, nd, 1);
  988   GEN VG = mfcoefs_i(G, nd, 1);
  989   ulong j, m = itou(gel(S,4));
  990 
  991   if (!n)
  992   {
  993     if (m > 0) { set_avma(av); return mkvec(gen_0); }
  994     return gerepilecopy(av, mkvec(gmul(gel(VF, 1), gel(VG, 1))));
  995   }
  996   tF = cgetg(nd+2, t_VEC);
  997   tG = cgetg(nd+2, t_VEC);
  998   res = NULL; gk = mf_get_gk(F); gl = mf_get_gk(G);
  999   /* pow[i,j+1] = i^j */
 1000   if (lg(S) == 6) chicompatlift(gel(S,5),&VF,&VG);
 1001   mpow = cgetg(m+2, t_MAT);
 1002   gel(mpow,1) = const_col(nd, gen_1);
 1003   for (j = 1; j <= m; j++)
 1004   {
 1005     GEN c = cgetg(nd+1, t_COL);
 1006     gel(mpow,j+1) = c;
 1007     for (i = 1; i <= nd; i++) gel(c,i) = muliu(gcoeff(mpow,i,j), i);
 1008   }
 1009   C = binomial(gaddgs(gk, m-1), m);
 1010   if (odd(m)) C = gneg(C);
 1011   for (j = 0; j <= m; j++)
 1012   { /* C = (-1)^(m-j) binom(m+l-1, j) binom(m+k-1,m-j) */
 1013     GEN c;
 1014     gel(tF,1) = j == 0? gel(VF,1): gen_0;
 1015     gel(tG,1) = j == m? gel(VG,1): gen_0;
 1016     gel(tF,2) = gel(VF,2); /* assume nd >= 1 */
 1017     gel(tG,2) = gel(VG,2);
 1018     for (i = 2; i <= nd; i++)
 1019     {
 1020       gel(tF, i+1) = gmul(gcoeff(mpow,i,j+1),   gel(VF, i+1));
 1021       gel(tG, i+1) = gmul(gcoeff(mpow,i,m-j+1), gel(VG, i+1));
 1022     }
 1023     c = gmul(C, c_deflate(n, d, RgV_mul_RgXn(tF, tG)));
 1024     res = res? gadd(res, c): c;
 1025     if (j < m)
 1026       C = gdiv(gmul(C, gmulsg(m-j, gaddgs(gl,m-j-1))),
 1027                gmulsg(-(j+1), gaddgs(gk,j)));
 1028   }
 1029   if (lg(S) == 6) res = chicompatfix(gel(S,5), res);
 1030   return gerepileupto(av, res);
 1031 }
 1032 /* linear combination \sum L[j] vecF[j] */
 1033 static GEN
 1034 c_linear(long n, long d, GEN F, GEN L, GEN dL)
 1035 {
 1036   pari_sp av = avma;
 1037   long j, l = lg(L);
 1038   GEN S = NULL;
 1039   for (j = 1; j < l; j++)
 1040   {
 1041     GEN c = gel(L,j);
 1042     if (gequal0(c)) continue;
 1043     c = gmul(c, mfcoefs_i(gel(F,j), n, d));
 1044     S = S? gadd(S,c): c;
 1045   }
 1046   if (!S) return zerovec(n+1);
 1047   if (!is_pm1(dL)) S = gdiv(S, dL);
 1048   return gerepileupto(av, S);
 1049 }
 1050 
 1051 /* B_d(T_j Trace^new) as t_MF_BD(t_MF_HECKE(t_MF_NEWTRACE)) or
 1052  * t_MF_HECKE(t_MF_NEWTRACE)
 1053  * or t_MF_NEWTRACE in level N. Set d and j, return t_MF_NEWTRACE component*/
 1054 static GEN
 1055 bhn_parse(GEN f, long *d, long *j)
 1056 {
 1057   long t = mf_get_type(f);
 1058   *d = *j = 1;
 1059   if (t == t_MF_BD) { *d = itos(gel(f,3)); f = gel(f,2); t = mf_get_type(f); }
 1060   if (t == t_MF_HECKE) { *j = gel(f,2)[1]; f = gel(f,3); }
 1061   return f;
 1062 }
 1063 /* f as above, return the t_MF_NEWTRACE component */
 1064 static GEN
 1065 bhn_newtrace(GEN f)
 1066 {
 1067   long t = mf_get_type(f);
 1068   if (t == t_MF_BD) { f = gel(f,2); t = mf_get_type(f); }
 1069   if (t == t_MF_HECKE) f = gel(f,3);
 1070   return f;
 1071 }
 1072 static int
 1073 ok_bhn_linear(GEN vf)
 1074 {
 1075   long i, N0 = 0, l = lg(vf);
 1076   GEN CHI, gk;
 1077   if (l == 1) return 1;
 1078   gk = mf_get_gk(gel(vf,1));
 1079   CHI = mf_get_CHI(gel(vf,1));
 1080   for (i = 1; i < l; i++)
 1081   {
 1082     GEN f = bhn_newtrace(gel(vf,i));
 1083     long N = mf_get_N(f);
 1084     if (mf_get_type(f) != t_MF_NEWTRACE) return 0;
 1085     if (N < N0) return 0; /* largest level must come last */
 1086     N0 = N;
 1087     if (!gequal(gk,mf_get_gk(f))) return 0; /* same k */
 1088     if (!gequal(gel(mf_get_CHI(f),2), gel(CHI,2))) return 0; /* same CHI */
 1089   }
 1090   return 1;
 1091 }
 1092 
 1093 /* vF not empty, same hypotheses as bhnmat_extend */
 1094 static GEN
 1095 bhnmat_extend_nocache(GEN M, long N, long n, long d, GEN vF)
 1096 {
 1097   cachenew_t cache;
 1098   long l = lg(vF);
 1099   GEN f;
 1100   if (l == 1) return M? M: cgetg(1, t_MAT);
 1101   f = bhn_newtrace(gel(vF,1)); /* N.B. mf_get_N(f) divides N */
 1102   init_cachenew(&cache, n*d, N, f);
 1103   M = bhnmat_extend(M, n, d, vF, &cache);
 1104   dbg_cachenew(&cache); return M;
 1105 }
 1106 /* c_linear of "bhn" mf closures, same hypotheses as bhnmat_extend */
 1107 static GEN
 1108 c_linear_bhn(long n, long d, GEN F)
 1109 {
 1110   pari_sp av;
 1111   GEN M, v, vF = gel(F,2), L = gel(F,3), dL = gel(F,4);
 1112   if (lg(L) == 1) return zerovec(n+1);
 1113   av = avma;
 1114   M = bhnmat_extend_nocache(NULL, mf_get_N(F), n, d, vF);
 1115   v = RgM_RgC_mul(M,L); settyp(v, t_VEC);
 1116   if (!is_pm1(dL)) v = gdiv(v, dL);
 1117   return gerepileupto(av, v);
 1118 }
 1119 
 1120 /* c in K, K := Q[X]/(T) vz = vector of consecutive powers of root z of T
 1121  * attached to an embedding s: K -> C. Return s(c) in C */
 1122 static GEN
 1123 Rg_embed1(GEN c, GEN vz)
 1124 {
 1125   long t = typ(c);
 1126   if (t == t_POLMOD) { c = gel(c,2); t = typ(c); }
 1127   if (t == t_POL) c = RgX_RgV_eval(c, vz);
 1128   return c;
 1129 }
 1130 /* return s(P) in C[X] */
 1131 static GEN
 1132 RgX_embed1(GEN P, GEN vz)
 1133 {
 1134   long i, l;
 1135   GEN Q = cgetg_copy(P, &l);
 1136   Q[1] = P[1];
 1137   for (i = 2; i < l; i++) gel(Q,i) = Rg_embed1(gel(P,i), vz);
 1138   return normalizepol_lg(Q,l); /* normally a no-op */
 1139 }
 1140 /* return s(P) in C^n */
 1141 static GEN
 1142 vecembed1(GEN P, GEN vz)
 1143 {
 1144   long i, l;
 1145   GEN Q = cgetg_copy(P, &l);
 1146   for (i = 1; i < l; i++) gel(Q,i) = Rg_embed1(gel(P,i), vz);
 1147   return Q;
 1148 }
 1149 /* P in L = K[X]/(U), K = Q[t]/T; s an embedding of K -> C attached
 1150  * to a root of T, extended to an embedding of L -> C attached to a root
 1151  * of s(U); vT powers of the root of T, vU powers of the root of s(U).
 1152  * Return s(P) in C^n */
 1153 static GEN
 1154 Rg_embed2(GEN P, long vt, GEN vT, GEN vU)
 1155 {
 1156   long i, l;
 1157   GEN Q;
 1158   P = liftpol_shallow(P);
 1159   if (typ(P) != t_POL) return P;
 1160   if (varn(P) == vt) return Rg_embed1(P, vT);
 1161   /* varn(P) == vx */
 1162   Q = cgetg_copy(P, &l); Q[1] = P[1];
 1163   for (i = 2; i < l; i++) gel(Q,i) = Rg_embed1(gel(P,i), vT);
 1164   return Rg_embed1(Q, vU);
 1165 }
 1166 static GEN
 1167 vecembed2(GEN P, long vt, GEN vT, GEN vU)
 1168 {
 1169   long i, l;
 1170   GEN Q = cgetg_copy(P, &l);
 1171   for (i = 1; i < l; i++) gel(Q,i) = Rg_embed2(gel(P,i), vt, vT, vU);
 1172   return Q;
 1173 }
 1174 static GEN
 1175 RgX_embed2(GEN P, long vt, GEN vT, GEN vU)
 1176 {
 1177   long i, l;
 1178   GEN Q = cgetg_copy(P, &l);
 1179   for (i = 2; i < l; i++) gel(Q,i) = Rg_embed2(gel(P,i), vt, vT, vU);
 1180   Q[1] = P[1]; return normalizepol_lg(Q,l);
 1181 }
 1182 /* embed polynomial f in variable vx [ may be a scalar ], E from getembed */
 1183 static GEN
 1184 RgX_embed(GEN f, long vx, GEN E)
 1185 {
 1186   GEN vT;
 1187   if (typ(f) != t_POL || varn(f) != vx) return mfembed(E, f);
 1188   if (lg(E) == 1) return f;
 1189   vT = gel(E,2);
 1190   if (lg(E) == 3)
 1191     f = RgX_embed1(f, vT);
 1192   else
 1193     f = RgX_embed2(f, varn(gel(E,1)), vT, gel(E,3));
 1194   return f;
 1195 }
 1196 /* embed vector, E from getembed */
 1197 GEN
 1198 mfvecembed(GEN E, GEN v)
 1199 {
 1200   GEN vT;
 1201   if (lg(E) == 1) return v;
 1202   vT = gel(E,2);
 1203   if (lg(E) == 3)
 1204     v = vecembed1(v, vT);
 1205   else
 1206     v = vecembed2(v, varn(gel(E,1)), vT, gel(E,3));
 1207   return v;
 1208 }
 1209 GEN
 1210 mfmatembed(GEN E, GEN f)
 1211 {
 1212   long i, l;
 1213   GEN g;
 1214   if (lg(E) == 1) return f;
 1215   g = cgetg_copy(f, &l);
 1216   for (i = 1; i < l; i++) gel(g,i) = mfvecembed(E, gel(f,i));
 1217   return g;
 1218 }
 1219 /* embed vector of polynomials in var vx */
 1220 static GEN
 1221 RgXV_embed(GEN f, long vx, GEN E)
 1222 {
 1223   long i, l;
 1224   GEN v;
 1225   if (lg(E) == 1) return f;
 1226   v = cgetg_copy(f, &l);
 1227   for (i = 1; i < l; i++) gel(v,i) = RgX_embed(gel(f,i), vx, E);
 1228   return v;
 1229 }
 1230 
 1231 /* embed scalar */
 1232 GEN
 1233 mfembed(GEN E, GEN f)
 1234 {
 1235   GEN vT;
 1236   if (lg(E) == 1) return f;
 1237   vT = gel(E,2);
 1238   if (lg(E) == 3)
 1239     f = Rg_embed1(f, vT);
 1240   else
 1241     f = Rg_embed2(f, varn(gel(E,1)), vT, gel(E,3));
 1242   return f;
 1243 }
 1244 /* vector of the sigma(f), sigma in vE */
 1245 static GEN
 1246 RgX_embedall(GEN f, long vx, GEN vE)
 1247 {
 1248   long i, l = lg(vE);
 1249   GEN v;
 1250   if (l == 2) return RgX_embed(f, vx, gel(vE,1));
 1251   v = cgetg(l, t_VEC);
 1252   for (i = 1; i < l; i++) gel(v,i) = RgX_embed(f, vx, gel(vE,i));
 1253   return v;
 1254 }
 1255 /* matrix whose colums are the sigma(v), sigma in vE */
 1256 static GEN
 1257 RgC_embedall(GEN v, GEN vE)
 1258 {
 1259   long j, l = lg(vE);
 1260   GEN M = cgetg(l, t_MAT);
 1261   for (j = 1; j < l; j++) gel(M,j) = mfvecembed(gel(vE,j), v);
 1262   return M;
 1263 }
 1264 /* vector of the sigma(v), sigma in vE */
 1265 static GEN
 1266 Rg_embedall_i(GEN v, GEN vE)
 1267 {
 1268   long j, l = lg(vE);
 1269   GEN M = cgetg(l, t_VEC);
 1270   for (j = 1; j < l; j++) gel(M,j) = mfembed(gel(vE,j), v);
 1271   return M;
 1272 }
 1273 /* vector of the sigma(v), sigma in vE; if #vE == 1, return v */
 1274 static GEN
 1275 Rg_embedall(GEN v, GEN vE)
 1276 { return (lg(vE) == 2)? mfembed(gel(vE,1), v): Rg_embedall_i(v, vE); }
 1277 
 1278 static GEN
 1279 c_div_i(long n, GEN S)
 1280 {
 1281   GEN F = gel(S,2), G = gel(S,3);
 1282   GEN a0, a0i, H;
 1283   F = mfcoefs_i(F, n, 1);
 1284   G = mfcoefs_i(G, n, 1);
 1285   if (lg(S) == 5) chicompatlift(gel(S,4),&F,&G);
 1286   F = RgV_to_ser_full(F);
 1287   G = RgV_to_ser_full(G);
 1288   a0 = polcoef_i(G, 0, -1); /* != 0 */
 1289   if (gequal1(a0)) a0 = a0i = NULL;
 1290   else
 1291   {
 1292     a0i = ginv(a0);
 1293     G = gmul(ser_unscale(G,a0), a0i);
 1294     F = gmul(ser_unscale(F,a0), a0i);
 1295   }
 1296   H = gdiv(F, G);
 1297   if (a0) H = ser_unscale(H,a0i);
 1298   H = sertovecslice(H, n);
 1299   if (lg(S) == 5) H = chicompatfix(gel(S,4), H);
 1300   return H;
 1301 }
 1302 static GEN
 1303 c_div(long n, long d, GEN S)
 1304 {
 1305   pari_sp av = avma;
 1306   GEN D = (d==1)? c_div_i(n, S): c_deflate(n, d, c_div_i(n*d, S));
 1307   return gerepilecopy(av, D);
 1308 }
 1309 
 1310 static GEN
 1311 c_shift(long n, long d, GEN F, GEN gsh)
 1312 {
 1313   pari_sp av = avma;
 1314   GEN vF;
 1315   long sh = itos(gsh), n1 = n*d + sh;
 1316   if (n1 < 0) return zerovec(n+1);
 1317   vF = mfcoefs_i(F, n1, 1);
 1318   if (sh < 0) vF = shallowconcat(zerovec(-sh), vF);
 1319   else vF = vecslice(vF, sh+1, n1+1);
 1320   return gerepilecopy(av, c_deflate(n, d, vF));
 1321 }
 1322 
 1323 static GEN
 1324 c_deriv(long n, long d, GEN F, GEN gm)
 1325 {
 1326   pari_sp av = avma;
 1327   GEN V = mfcoefs_i(F, n, d), res;
 1328   long i, m = itos(gm);
 1329   if (!m) return V;
 1330   res = cgetg(n+2, t_VEC); gel(res,1) = gen_0;
 1331   if (m < 0)
 1332   { for (i=1; i <= n; i++) gel(res, i+1) = gdiv(gel(V, i+1), powuu(i,-m)); }
 1333   else
 1334   { for (i=1; i <= n; i++) gel(res, i+1) = gmul(gel(V,i+1), powuu(i,m)); }
 1335   return gerepileupto(av, res);
 1336 }
 1337 
 1338 static GEN
 1339 c_derivE2(long n, long d, GEN F, GEN gm)
 1340 {
 1341   pari_sp av = avma;
 1342   GEN VF, VE, res, tmp, gk;
 1343   long i, m = itos(gm), nd;
 1344   if (m == 0) return mfcoefs_i(F, n, d);
 1345   nd = n*d;
 1346   VF = mfcoefs_i(F, nd, 1); VE = mfcoefs_i(mfEk(2), nd, 1);
 1347   gk = mf_get_gk(F);
 1348   if (m == 1)
 1349   {
 1350     res = cgetg(n+2, t_VEC);
 1351     for (i = 0; i <= n; i++) gel(res, i+1) = gmulsg(i, gel(VF, i*d+1));
 1352     tmp = c_deflate(n, d, RgV_mul_RgXn(VF, VE));
 1353     return gerepileupto(av, gsub(res, gmul(gdivgs(gk, 12), tmp)));
 1354   }
 1355   else
 1356   {
 1357     long j;
 1358     for (j = 1; j <= m; j++)
 1359     {
 1360       tmp = RgV_mul_RgXn(VF, VE);
 1361       for (i = 0; i <= nd; i++) gel(VF, i+1) = gmulsg(i, gel(VF, i+1));
 1362       VF = gsub(VF, gmul(gdivgs(gaddgs(gk, 2*(j-1)), 12), tmp));
 1363     }
 1364     return gerepilecopy(av, c_deflate(n, d, VF));
 1365   }
 1366 }
 1367 
 1368 /* Twist by the character (D/.) */
 1369 static GEN
 1370 c_twist(long n, long d, GEN F, GEN D)
 1371 {
 1372   pari_sp av = avma;
 1373   GEN V = mfcoefs_i(F, n, d), res = cgetg(n+2, t_VEC);
 1374   long i;
 1375   for (i = 0; i <= n; i++)
 1376     gel(res, i + 1) = gmulsg(krois(D, i), gel(V, i+1));
 1377   return gerepileupto(av, res);
 1378 }
 1379 
 1380 /* form F given by closure, compute T(n)(F) as closure */
 1381 static GEN
 1382 c_hecke(long m, long l, GEN DATA, GEN F)
 1383 {
 1384   pari_sp av = avma;
 1385   return gerepilecopy(av, hecke_i(m, l, NULL, F, DATA));
 1386 }
 1387 static GEN
 1388 c_const(long n, long d, GEN C)
 1389 {
 1390   GEN V = zerovec(n+1);
 1391   long i, j, l = lg(C);
 1392   if (l > d*n+2) l = d*n+2;
 1393   for (i = j = 1; i < l; i+=d, j++) gel(V, j) = gcopy(gel(C,i));
 1394   return V;
 1395 }
 1396 
 1397 /* m > 0 */
 1398 static GEN
 1399 eta3_ZXn(long m)
 1400 {
 1401   long l = m+2, n, k;
 1402   GEN P = cgetg(l,t_POL);
 1403   P[1] = evalsigne(1)|evalvarn(0);
 1404   for (n = 2; n < l; n++) gel(P,n) = gen_0;
 1405   for (n = k = 0;; n++)
 1406   {
 1407     if (k + n >= m) { setlg(P, k+3); return P; }
 1408     k += n;
 1409     /* now k = n(n+1) / 2 */
 1410     gel(P, k+2) = odd(n)? utoineg(2*n+1): utoipos(2*n+1);
 1411   }
 1412 }
 1413 
 1414 static GEN
 1415 c_delta(long n, long d)
 1416 {
 1417   pari_sp ltop = avma;
 1418   long N = n*d;
 1419   GEN e;
 1420   if (!N) return mkvec(gen_0);
 1421   e = eta3_ZXn(N);
 1422   e = ZXn_sqr(e,N);
 1423   e = ZXn_sqr(e,N);
 1424   e = ZXn_sqr(e,N); /* eta(x)^24 */
 1425   settyp(e, t_VEC);
 1426   gel(e,1) = gen_0; /* Delta(x) = x*eta(x)^24 as a t_VEC */
 1427   return gerepilecopy(ltop, c_deflate(n, d, e));
 1428 }
 1429 
 1430 /* return s(d) such that s|f <=> d | f^2 */
 1431 static long
 1432 mysqrtu(ulong d)
 1433 {
 1434   GEN fa = myfactoru(d), P = gel(fa,1), E = gel(fa,2);
 1435   long l = lg(P), i, s = 1;
 1436   for (i = 1; i < l; i++) s *= upowuu(P[i], (E[i]+1)>>1);
 1437   return s;
 1438 }
 1439 static GEN
 1440 c_theta(long n, long d, GEN psi)
 1441 {
 1442   long lim = usqrt(n*d), F = mfcharmodulus(psi), par = mfcharparity(psi);
 1443   long f, d2 = d == 1? 1: mysqrtu(d);
 1444   GEN V = zerovec(n + 1);
 1445   for (f = d2; f <= lim; f += d2)
 1446     if (ugcd(F, f) == 1)
 1447     {
 1448       pari_sp av = avma;
 1449       GEN c = mfchareval(psi, f);
 1450       gel(V, f*f/d + 1) = gerepileupto(av, par < 0 ? gmulgs(c,2*f) : gmul2n(c,1));
 1451     }
 1452   if (F == 1) gel(V, 1) = gen_1;
 1453   return V; /* no gerepile needed */
 1454 }
 1455 
 1456 static GEN
 1457 c_etaquo(long n, long d, GEN eta, GEN gs)
 1458 {
 1459   pari_sp av = avma;
 1460   long s = itos(gs), nd = n*d, nds = nd - s + 1;
 1461   GEN c;
 1462   if (nds <= 0) return zerovec(n+1);
 1463   c = RgX_to_RgC(eta_product_ZXn(eta, nds), nds); settyp(c, t_VEC);
 1464   if (s > 0) c = shallowconcat(zerovec(s), c);
 1465   return gerepilecopy(av, c_deflate(n, d, c));
 1466 }
 1467 
 1468 static GEN
 1469 c_ell(long n, long d, GEN E)
 1470 {
 1471   pari_sp av = avma;
 1472   GEN v;
 1473   if (d == 1) return concat(gen_0, anell(E, n));
 1474   v = shallowconcat(gen_0, anell(E, n*d));
 1475   return gerepilecopy(av, c_deflate(n, d, v));
 1476 }
 1477 
 1478 static GEN
 1479 c_cusptrace(long n, long d, GEN F)
 1480 {
 1481   pari_sp av = avma;
 1482   GEN D = gel(F,2), res = cgetg(n+2, t_VEC);
 1483   long i, N = mf_get_N(F), k = mf_get_k(F);
 1484   gel(res, 1) = gen_0;
 1485   for (i = 1; i <= n; i++)
 1486     gel(res, i+1) = mfcusptrace_i(N, k, i*d, mydivisorsu(i*d), D);
 1487   return gerepilecopy(av, res);
 1488 }
 1489 
 1490 static GEN
 1491 c_newtrace(long n, long d, GEN F)
 1492 {
 1493   pari_sp av = avma;
 1494   cachenew_t cache;
 1495   long N = mf_get_N(F);
 1496   GEN v;
 1497   init_cachenew(&cache, n*d, N, F);
 1498   v = colnewtrace(0, n, d, N, mf_get_k(F), &cache);
 1499   settyp(v, t_VEC); return gerepilecopy(av, v);
 1500 }
 1501 
 1502 static GEN
 1503 c_Bd(long n, long d, GEN F, GEN A)
 1504 {
 1505   pari_sp av = avma;
 1506   long a = itou(A), ad = ugcd(a,d), aad = a/ad, i, j;
 1507   GEN w, v = mfcoefs_i(F, n/aad, d/ad);
 1508   if (a == 1) return v;
 1509   n++; w = zerovec(n);
 1510   for (i = j = 1; j <= n; i++, j += aad) gel(w,j) = gcopy(gel(v,i));
 1511   return gerepileupto(av, w);
 1512 }
 1513 
 1514 static GEN
 1515 c_dihedral(long n, long d, GEN bnr, GEN w, GEN k0j)
 1516 {
 1517   pari_sp av = avma;
 1518   GEN V = dihan(bnr, w, k0j, n*d);
 1519   GEN Tinit = gel(w,3), Pm = gel(Tinit,1);
 1520   GEN A = c_deflate(n, d, V);
 1521   if (degpol(Pm) == 1 || RgV_is_ZV(A)) return gerepilecopy(av, A);
 1522   return gerepileupto(av, gmodulo(A, Pm));
 1523 }
 1524 
 1525 static GEN
 1526 c_mfEH(long n, long d, GEN F)
 1527 {
 1528   pari_sp av = avma;
 1529   GEN v, M, A;
 1530   long i, r = mf_get_r(F);
 1531   if (n == 1)
 1532     return gerepilecopy(av, mkvec2(mfEHcoef(r,0),mfEHcoef(r,d)));
 1533   /* speedup mfcoef */
 1534   if (r == 1)
 1535   {
 1536     v = cgetg(n+2, t_VEC);
 1537     gel(v,1) = sstoQ(-1,12);
 1538     for (i = 1; i <= n; i++)
 1539     {
 1540       long id = i*d, a = id & 3;
 1541       gel(v,i+1) = (a==1 || a==2)? gen_0: sstoQ(hclassno6u(id), 6);
 1542     }
 1543     return v; /* no gerepile needed */
 1544   }
 1545   M = mfEHmat(n*d+1,r);
 1546   if (d > 1)
 1547   {
 1548     long l = lg(M);
 1549     for (i = 1; i < l; i++) gel(M,i) = c_deflate(n, d, gel(M,i));
 1550   }
 1551   A = gel(F,2); /* [num(B), den(B)] */
 1552   v = RgC_Rg_div(RgM_RgC_mul(M, gel(A,1)), gel(A,2));
 1553   settyp(v,t_VEC); return gerepileupto(av, v);
 1554 }
 1555 
 1556 static GEN
 1557 c_mfeisen(long n, long d, GEN F)
 1558 {
 1559   pari_sp av = avma;
 1560   GEN v, vchi, E0, P, T, CHI, gk = mf_get_gk(F);
 1561   long i, k;
 1562   if (typ(gk) != t_INT) return c_mfEH(n, d, F);
 1563   k = itou(gk);
 1564   vchi = gel(F,2);
 1565   E0 = gel(vchi,1);
 1566   T = gel(vchi,2);
 1567   P = gel(T,1);
 1568   CHI = gel(vchi,3);
 1569   v = cgetg(n+2, t_VEC);
 1570   gel(v, 1) = gcopy(E0); /* E(0) */
 1571   if (lg(vchi) == 5)
 1572   { /* E_k(chi1,chi2) */
 1573     GEN CHI2 = gel(vchi,4), F3 = gel(F,3);
 1574     long ord = F3[1], j = F3[2];
 1575     for (i = 1; i <= n; i++) gel(v, i+1) = sigchi2(k, CHI, CHI2, i*d, ord);
 1576     v = QabV_tracerel(T, j, v);
 1577   }
 1578   else
 1579   { /* E_k(chi) */
 1580     for (i = 1; i <= n; i++) gel(v, i+1) = sigchi(k, CHI, i*d);
 1581   }
 1582   if (degpol(P) != 1 && !RgV_is_QV(v)) return gerepileupto(av, gmodulo(v, P));
 1583   return gerepilecopy(av, v);
 1584 }
 1585 
 1586 /* L(chi_D, 1-k) */
 1587 static GEN
 1588 lfunquadneg_naive(long D, long k)
 1589 {
 1590   GEN B, dS, S = gen_0;
 1591   long r, N = labs(D);
 1592   pari_sp av;
 1593   if (k == 1 && N == 1) return gneg(ghalf);
 1594   /* B = N^k * denom(B) * B(x/N) */
 1595   B = ZX_rescale(Q_remove_denom(bernpol(k, 0), &dS), utoi(N));
 1596   dS = mul_denom(dS, stoi(-N*k));
 1597   av = avma;
 1598   for (r = 0; r < N; r++)
 1599   {
 1600     long c = kross(D, r);
 1601     if (c)
 1602     {
 1603       GEN tmp = poleval(B, utoi(r));
 1604       S = c > 0 ? addii(S, tmp) : subii(S, tmp);
 1605       S = gerepileuptoint(av, S);
 1606     }
 1607   }
 1608   return gdiv(S, dS);
 1609 }
 1610 
 1611 /* Returns vector of coeffs from F[0], F[d], ..., F[d*n] */
 1612 static GEN
 1613 mfcoefs_i(GEN F, long n, long d)
 1614 {
 1615   if (n < 0) return gen_0;
 1616   switch(mf_get_type(F))
 1617   {
 1618     case t_MF_CONST: return c_const(n, d, gel(F,2));
 1619     case t_MF_EISEN: return c_mfeisen(n, d, F);
 1620     case t_MF_Ek: return c_Ek(n, d, F);
 1621     case t_MF_DELTA: return c_delta(n, d);
 1622     case t_MF_THETA: return c_theta(n, d, gel(F,2));
 1623     case t_MF_ETAQUO: return c_etaquo(n, d, gel(F,2), gel(F,3));
 1624     case t_MF_ELL: return c_ell(n, d, gel(F,2));
 1625     case t_MF_MUL: return c_mul(n, d, F);
 1626     case t_MF_POW: return c_pow(n, d, F);
 1627     case t_MF_BRACKET: return c_bracket(n, d, F);
 1628     case t_MF_LINEAR: return c_linear(n, d, gel(F,2), gel(F,3), gel(F,4));
 1629     case t_MF_LINEAR_BHN: return c_linear_bhn(n, d, F);
 1630     case t_MF_DIV: return c_div(n, d, F);
 1631     case t_MF_SHIFT: return c_shift(n, d, gel(F,2), gel(F,3));
 1632     case t_MF_DERIV: return c_deriv(n, d, gel(F,2), gel(F,3));
 1633     case t_MF_DERIVE2: return c_derivE2(n, d, gel(F,2), gel(F,3));
 1634     case t_MF_TWIST: return c_twist(n, d, gel(F,2), gel(F,3));
 1635     case t_MF_HECKE: return c_hecke(n, d, gel(F,2), gel(F,3));
 1636     case t_MF_BD: return c_Bd(n, d, gel(F,2), gel(F,3));
 1637     case t_MF_TRACE: return c_cusptrace(n, d, F);
 1638     case t_MF_NEWTRACE: return c_newtrace(n, d, F);
 1639     case t_MF_DIHEDRAL: return c_dihedral(n, d, gel(F,2), gel(F,3), gel(F,4));
 1640     default: pari_err_TYPE("mfcoefs",F); return NULL;/*LCOV_EXCL_LINE*/
 1641   }
 1642 }
 1643 
 1644 static GEN
 1645 matdeflate(long n, long d, GEN M)
 1646 {
 1647   long i, l;
 1648   GEN A;
 1649   /*  if (d == 1) return M; */
 1650   A = cgetg_copy(M,&l);
 1651   for (i = 1; i < l; i++) gel(A,i) = c_deflate(n,d,gel(M,i));
 1652   return A;
 1653 }
 1654 static int
 1655 space_is_cusp(long space) { return space != mf_FULL && space != mf_EISEN; }
 1656 /* safe with flraw mf */
 1657 static GEN
 1658 mfcoefs_mf(GEN mf, long n, long d)
 1659 {
 1660   GEN MS, ME, E = MF_get_E(mf), S = MF_get_S(mf), M = MF_get_M(mf);
 1661   long lE = lg(E), lS = lg(S), l = lE+lS-1;
 1662 
 1663   if (l == 1) return cgetg(1, t_MAT);
 1664   if (typ(M) == t_MAT && lg(M) != 1 && (n+1)*d < nbrows(M))
 1665     return matdeflate(n, d, M); /*cached; lg = 1 is possible from mfinit */
 1666   ME = (lE == 1)? cgetg(1, t_MAT): mfvectomat(E, n, d);
 1667   if (lS == 1)
 1668     MS = cgetg(1, t_MAT);
 1669   else if (mf_get_type(gel(S,1)) == t_MF_DIV) /*k 1/2-integer or k=1 (exotic)*/
 1670     MS = matdeflate(n,d, mflineardivtomat(MF_get_N(mf), S, n*d));
 1671   else if (MF_get_k(mf) == 1) /* k = 1 (dihedral) */
 1672   {
 1673     GEN M = mfvectomat(gmael(S,1,2), n, d);
 1674     long i;
 1675     MS = cgetg(lS, t_MAT);
 1676     for (i = 1; i < lS; i++)
 1677     {
 1678       GEN f = gel(S,i), dc = gel(f,4), c = RgM_RgC_mul(M, gel(f,3));
 1679       if (!equali1(dc)) c = RgC_Rg_div(c,dc);
 1680       gel(MS,i) = c;
 1681     }
 1682   }
 1683   else /* k >= 2 integer */
 1684     MS = bhnmat_extend_nocache(NULL, MF_get_N(mf), n, d, S);
 1685   return shallowconcat(ME,MS);
 1686 }
 1687 GEN
 1688 mfcoefs(GEN F, long n, long d)
 1689 {
 1690   if (!checkmf_i(F))
 1691   {
 1692     pari_sp av = avma;
 1693     GEN mf = checkMF_i(F); if (!mf) pari_err_TYPE("mfcoefs", F);
 1694     return gerepilecopy(av, mfcoefs_mf(mf,n,d));
 1695   }
 1696   if (d <= 0) pari_err_DOMAIN("mfcoefs", "d", "<=", gen_0, stoi(d));
 1697   if (n < 0) return cgetg(1, t_VEC);
 1698   return mfcoefs_i(F, n, d);
 1699 }
 1700 
 1701 /* assume k >= 0 */
 1702 static GEN
 1703 mfak_i(GEN F, long k)
 1704 {
 1705   if (!k) return gel(mfcoefs_i(F,0,1), 1);
 1706   return gel(mfcoefs_i(F,1,k), 2);
 1707 }
 1708 GEN
 1709 mfcoef(GEN F, long n)
 1710 {
 1711   pari_sp av = avma;
 1712   if (!checkmf_i(F)) pari_err_TYPE("mfcoef",F);
 1713   return n < 0? gen_0: gerepilecopy(av, mfak_i(F, n));
 1714 }
 1715 
 1716 static GEN
 1717 paramconst() { return tagparams(t_MF_CONST, mkNK(1,0,mfchartrivial())); }
 1718 static GEN
 1719 mftrivial(void) { retmkvec2(paramconst(), cgetg(1,t_VEC)); }
 1720 static GEN
 1721 mf1(void) { retmkvec2(paramconst(), mkvec(gen_1)); }
 1722 
 1723 /* induce mfchar CHI to G */
 1724 static GEN
 1725 induce(GEN G, GEN CHI)
 1726 {
 1727   GEN o, chi;
 1728   if (typ(CHI) == t_INT) /* Kronecker */
 1729   {
 1730     chi = znchar_quad(G, CHI);
 1731     o = ZV_equal0(chi)? gen_1: gen_2;
 1732     CHI = mkvec4(G,chi,o,cgetg(1,t_VEC));
 1733   }
 1734   else
 1735   {
 1736     if (mfcharmodulus(CHI) == itos(znstar_get_N(G))) return CHI;
 1737     CHI = leafcopy(CHI);
 1738     chi = zncharinduce(gel(CHI,1), gel(CHI,2), G);
 1739     gel(CHI,1) = G;
 1740     gel(CHI,2) = chi;
 1741   }
 1742   return CHI;
 1743 }
 1744 /* induce mfchar CHI to znstar(G) */
 1745 static GEN
 1746 induceN(long N, GEN CHI)
 1747 {
 1748   if (mfcharmodulus(CHI) != N) CHI = induce(znstar0(utoipos(N),1), CHI);
 1749   return CHI;
 1750 }
 1751 /* *pCHI1 and *pCHI2 are mfchar, induce to common modulus */
 1752 static void
 1753 char2(GEN *pCHI1, GEN *pCHI2)
 1754 {
 1755   GEN CHI1 = *pCHI1, G1 = gel(CHI1,1), N1 = znstar_get_N(G1);
 1756   GEN CHI2 = *pCHI2, G2 = gel(CHI2,1), N2 = znstar_get_N(G2);
 1757   if (!equalii(N1,N2))
 1758   {
 1759     GEN G, d = gcdii(N1,N2);
 1760     if      (equalii(N2,d)) *pCHI2 = induce(G1, CHI2);
 1761     else if (equalii(N1,d)) *pCHI1 = induce(G2, CHI1);
 1762     else
 1763     {
 1764       if (!equali1(d)) N2 = diviiexact(N2,d);
 1765       G = znstar0(mulii(N1,N2), 1);
 1766       *pCHI1 = induce(G, CHI1);
 1767       *pCHI2 = induce(G, CHI2);
 1768     }
 1769   }
 1770 }
 1771 /* mfchar or charinit wrt same modulus; outputs a mfchar */
 1772 static GEN
 1773 mfcharmul_i(GEN CHI1, GEN CHI2)
 1774 {
 1775   GEN G = gel(CHI1,1), chi3 = zncharmul(G, gel(CHI1,2), gel(CHI2,2));
 1776   return mfcharGL(G, chi3);
 1777 }
 1778 /* mfchar or charinit; outputs a mfchar */
 1779 static GEN
 1780 mfcharmul(GEN CHI1, GEN CHI2)
 1781 {
 1782   char2(&CHI1, &CHI2); return mfcharmul_i(CHI1,CHI2);
 1783 }
 1784 /* mfchar or charinit; outputs a mfchar */
 1785 static GEN
 1786 mfcharpow(GEN CHI, GEN n)
 1787 {
 1788   GEN G, chi;
 1789   G = gel(CHI,1); chi = zncharpow(G, gel(CHI,2), n);
 1790   return mfchartoprimitive(mfcharGL(G, chi), NULL);
 1791 }
 1792 /* mfchar or charinit wrt same modulus; outputs a mfchar */
 1793 static GEN
 1794 mfchardiv_i(GEN CHI1, GEN CHI2)
 1795 {
 1796   GEN G = gel(CHI1,1), chi3 = znchardiv(G, gel(CHI1,2), gel(CHI2,2));
 1797   return mfcharGL(G, chi3);
 1798 }
 1799 /* mfchar or charinit; outputs a mfchar */
 1800 static GEN
 1801 mfchardiv(GEN CHI1, GEN CHI2)
 1802 {
 1803   char2(&CHI1, &CHI2); return mfchardiv_i(CHI1,CHI2);
 1804 }
 1805 static GEN
 1806 mfcharconj(GEN CHI)
 1807 {
 1808   CHI = leafcopy(CHI);
 1809   gel(CHI,2) = zncharconj(gel(CHI,1), gel(CHI,2));
 1810   return CHI;
 1811 }
 1812 
 1813 /* CHI mfchar, assume 4 | N. Multiply CHI by \chi_{-4} */
 1814 static GEN
 1815 mfchilift(GEN CHI, long N)
 1816 {
 1817   CHI = induceN(N, CHI);
 1818   return mfcharmul_i(CHI, induce(gel(CHI,1), stoi(-4)));
 1819 }
 1820 /* CHI defined mod N, N4 = N/4;
 1821  * if CHI is defined mod N4 return CHI;
 1822  * else if CHI' = CHI*(-4,.) is defined mod N4, return CHI' (primitive)
 1823  * else error */
 1824 static GEN
 1825 mfcharchiliftprim(GEN CHI, long N4)
 1826 {
 1827   long FC = mfcharconductor(CHI);
 1828   GEN CHIP;
 1829   if (N4 % FC == 0) return CHI;
 1830   CHIP = mfchartoprimitive(mfchilift(CHI, N4 << 2), &FC);
 1831   if (N4 % FC) pari_err_TYPE("mfkohnenbasis [incorrect CHI]", CHI);
 1832   return CHIP;
 1833 }
 1834 /* ensure CHI(-1) = (-1)^k [k integer] or 1 [half-integer], by multiplying
 1835  * by (-4/.) if needed */
 1836 static GEN
 1837 mfchiadjust(GEN CHI, GEN gk, long N)
 1838 {
 1839   long par = mfcharparity(CHI);
 1840   if (typ(gk) == t_INT &&  mpodd(gk)) par = -par;
 1841   return par == 1 ? CHI : mfchilift(CHI, N);
 1842 }
 1843 
 1844 static GEN
 1845 mfsamefield(GEN T, GEN P, GEN Q)
 1846 {
 1847   if (degpol(P) == 1) return Q;
 1848   if (degpol(Q) == 1) return P;
 1849   if (!gequal(P,Q)) pari_err_TYPE("mfsamefield [different fields]",mkvec2(P,Q));
 1850   if (T) err_cyclo();
 1851   return P;
 1852 }
 1853 
 1854 GEN
 1855 mfmul(GEN f, GEN g)
 1856 {
 1857   pari_sp av = avma;
 1858   GEN T, N, K, NK, CHI, CHIf, CHIg;
 1859   if (!checkmf_i(f)) pari_err_TYPE("mfmul",f);
 1860   if (!checkmf_i(g)) pari_err_TYPE("mfmul",g);
 1861   N = lcmii(mf_get_gN(f), mf_get_gN(g));
 1862   K = gadd(mf_get_gk(f), mf_get_gk(g));
 1863   CHIf = mf_get_CHI(f);
 1864   CHIg = mf_get_CHI(g);
 1865   CHI = mfchiadjust(mfcharmul(CHIf,CHIg), K, itos(N));
 1866   T = chicompat(CHI, CHIf, CHIg);
 1867   NK = mkgNK(N, K, CHI, mfsamefield(T, mf_get_field(f), mf_get_field(g)));
 1868   return gerepilecopy(av, T? tag3(t_MF_MUL,NK,f,g,T): tag2(t_MF_MUL,NK,f,g));
 1869 }
 1870 GEN
 1871 mfpow(GEN f, long n)
 1872 {
 1873   pari_sp av = avma;
 1874   GEN T, KK, NK, gn, CHI, CHIf;
 1875   if (!checkmf_i(f)) pari_err_TYPE("mfpow",f);
 1876   if (!n) return mf1();
 1877   if (n == 1) return gcopy(f);
 1878   KK = gmulsg(n,mf_get_gk(f));
 1879   gn = stoi(n);
 1880   CHIf = mf_get_CHI(f);
 1881   CHI = mfchiadjust(mfcharpow(CHIf,gn), KK, mf_get_N(f));
 1882   T = chicompat(CHI, CHIf, CHIf);
 1883   NK = mkgNK(mf_get_gN(f), KK, CHI, mf_get_field(f));
 1884   return gerepilecopy(av, T? tag3(t_MF_POW,NK,f,gn,T): tag2(t_MF_POW,NK,f,gn));
 1885 }
 1886 GEN
 1887 mfbracket(GEN f, GEN g, long m)
 1888 {
 1889   pari_sp av = avma;
 1890   GEN T, N, K, NK, CHI, CHIf, CHIg;
 1891   if (!checkmf_i(f)) pari_err_TYPE("mfbracket",f);
 1892   if (!checkmf_i(g)) pari_err_TYPE("mfbracket",g);
 1893   if (m < 0) pari_err_TYPE("mfbracket [m<0]",stoi(m));
 1894   K = gaddgs(gadd(mf_get_gk(f), mf_get_gk(g)), 2*m);
 1895   if (gsigne(K) < 0) pari_err_IMPL("mfbracket for this form");
 1896   N = lcmii(mf_get_gN(f), mf_get_gN(g));
 1897   CHIf = mf_get_CHI(f);
 1898   CHIg = mf_get_CHI(g);
 1899   CHI = mfcharmul(CHIf, CHIg);
 1900   CHI = mfchiadjust(CHI, K, itou(N));
 1901   T = chicompat(CHI, CHIf, CHIg);
 1902   NK = mkgNK(N, K, CHI, mfsamefield(T, mf_get_field(f), mf_get_field(g)));
 1903   return gerepilecopy(av, T? tag4(t_MF_BRACKET, NK, f, g, utoi(m), T)
 1904                            : tag3(t_MF_BRACKET, NK, f, g, utoi(m)));
 1905 }
 1906 
 1907 /* remove 0 entries in L */
 1908 static int
 1909 mflinear_strip(GEN *pF, GEN *pL)
 1910 {
 1911   pari_sp av = avma;
 1912   GEN F = *pF, L = *pL;
 1913   long i, j, l = lg(L);
 1914   GEN F2 = cgetg(l, t_VEC), L2 = cgetg(l, t_VEC);
 1915   for (i = j = 1; i < l; i++)
 1916   {
 1917     if (gequal0(gel(L,i))) continue;
 1918     gel(F2,j) = gel(F,i);
 1919     gel(L2,j) = gel(L,i); j++;
 1920   }
 1921   if (j == l) set_avma(av);
 1922   else
 1923   {
 1924     setlg(F2,j); *pF = F2;
 1925     setlg(L2,j); *pL = L2;
 1926   }
 1927   return (j > 1);
 1928 }
 1929 static GEN
 1930 taglinear_i(long t, GEN NK, GEN F, GEN L)
 1931 {
 1932   GEN dL;
 1933   L = Q_remove_denom(L, &dL); if (!dL) dL = gen_1;
 1934   return tag3(t, NK, F, L, dL);
 1935 }
 1936 static GEN
 1937 taglinear(GEN NK, GEN F, GEN L)
 1938 {
 1939   long t = ok_bhn_linear(F)? t_MF_LINEAR_BHN: t_MF_LINEAR;
 1940    return taglinear_i(t, NK, F, L);
 1941 }
 1942 /* assume F has parameters NK = [N,K,CHI] */
 1943 static GEN
 1944 mflinear_i(GEN NK, GEN F, GEN L)
 1945 {
 1946   if (!mflinear_strip(&F,&L)) return mftrivial();
 1947   return taglinear(NK, F,L);
 1948 }
 1949 static GEN
 1950 mflinear_bhn(GEN mf, GEN L)
 1951 {
 1952   long i, l;
 1953   GEN P, NK, F = MF_get_S(mf);
 1954   if (!mflinear_strip(&F,&L)) return mftrivial();
 1955   l = lg(L); P = pol_x(1);
 1956   for (i = 1; i < l; i++)
 1957   {
 1958     GEN c = gel(L,i);
 1959     if (typ(c) == t_POLMOD && varn(gel(c,1)) == 1)
 1960       P = mfsamefield(NULL, P, gel(c,1));
 1961   }
 1962   NK = mkgNK(MF_get_gN(mf), MF_get_gk(mf), MF_get_CHI(mf), P);
 1963   return taglinear_i(t_MF_LINEAR_BHN,  NK, F,L);
 1964 }
 1965 
 1966 /* F vector of forms with same weight and character but varying level, return
 1967  * global [N,k,chi,P] */
 1968 static GEN
 1969 vecmfNK(GEN F)
 1970 {
 1971   long i, l = lg(F);
 1972   GEN N, f;
 1973   if (l == 1) return mkNK(1, 0, mfchartrivial());
 1974   f = gel(F,1); N = mf_get_gN(f);
 1975   for (i = 2; i < l; i++) N = lcmii(N, mf_get_gN(gel(F,i)));
 1976   return mkgNK(N, mf_get_gk(f), mf_get_CHI(f), mf_get_field(f));
 1977 }
 1978 /* do not use mflinear: mflineardivtomat rely on F being constant across the
 1979  * basis where mflinear strips the ones matched by 0 coeffs. Assume k and CHI
 1980  * constant, N is allowed to vary. */
 1981 static GEN
 1982 vecmflinear(GEN F, GEN C)
 1983 {
 1984   long i, t, l = lg(C);
 1985   GEN NK, v = cgetg(l, t_VEC);
 1986   if (l == 1) return v;
 1987   t = ok_bhn_linear(F)? t_MF_LINEAR_BHN: t_MF_LINEAR;
 1988   NK = vecmfNK(F);
 1989   for (i = 1; i < l; i++) gel(v,i) = taglinear_i(t, NK, F, gel(C,i));
 1990   return v;
 1991 }
 1992 /* vecmflinear(F,C), then divide everything by E, which has valuation 0 */
 1993 static GEN
 1994 vecmflineardiv0(GEN F, GEN C, GEN E)
 1995 {
 1996   GEN v = vecmflinear(F, C);
 1997   long i, l = lg(v);
 1998   if (l == 1) return v;
 1999   gel(v,1) = mfdiv_val(gel(v,1), E, 0);
 2000   for (i = 2; i < l; i++)
 2001   { /* v[i] /= E */
 2002     GEN f = shallowcopy(gel(v,1));
 2003     gel(f,2) = gel(v,i);
 2004     gel(v,i) = f;
 2005   }
 2006   return v;
 2007 }
 2008 
 2009 /* Non empty linear combination of linear combinations of same
 2010  * F_j=\sum_i \mu_{i,j}G_i so R = \sum_i (\sum_j(\la_j\mu_{i,j})) G_i */
 2011 static GEN
 2012 mflinear_linear(GEN F, GEN L, int strip)
 2013 {
 2014   long l = lg(F), j;
 2015   GEN vF, M = cgetg(l, t_MAT);
 2016   L = shallowcopy(L);
 2017   for (j = 1; j < l; j++)
 2018   {
 2019     GEN f = gel(F,j), c = gel(f,3), d = gel(f,4);
 2020     if (typ(c) == t_VEC) c = shallowtrans(c);
 2021     if (!isint1(d)) gel(L,j) = gdiv(gel(L,j),d);
 2022     gel(M,j) = c;
 2023   }
 2024   vF = gmael(F,1,2); L = RgM_RgC_mul(M,L);
 2025   if (strip && !mflinear_strip(&vF,&L)) return mftrivial();
 2026   return taglinear(vecmfNK(vF), vF, L);
 2027 }
 2028 /* F nonempty vector of forms of the form mfdiv(mflinear(B,v), E) where E
 2029  * does not vanish at oo, or mflinear(B,v). Apply mflinear(F, L) */
 2030 static GEN
 2031 mflineardiv_linear(GEN F, GEN L, int strip)
 2032 {
 2033   long l = lg(F), j;
 2034   GEN v, E, f;
 2035   if (lg(L) != l) pari_err_DIM("mflineardiv_linear");
 2036   f = gel(F,1); /* l > 1 */
 2037   if (mf_get_type(f) != t_MF_DIV) return mflinear_linear(F,L,strip);
 2038   E = gel(f,3);
 2039   v = cgetg(l, t_VEC);
 2040   for (j = 1; j < l; j++) { GEN f = gel(F,j); gel(v,j) = gel(f,2); }
 2041   return mfdiv_val(mflinear_linear(v,L,strip), E, 0);
 2042 }
 2043 static GEN
 2044 vecmflineardiv_linear(GEN F, GEN M)
 2045 {
 2046   long i, l = lg(M);
 2047   GEN v = cgetg(l, t_VEC);
 2048   for (i = 1; i < l; i++) gel(v,i) = mflineardiv_linear(F, gel(M,i), 0);
 2049   return v;
 2050 }
 2051 
 2052 static GEN
 2053 tobasis(GEN mf, GEN F, GEN L)
 2054 {
 2055   if (checkmf_i(L) && mf) return mftobasis(mf, L, 0);
 2056   if (typ(F) != t_VEC) pari_err_TYPE("mflinear",F);
 2057   if (!is_vec_t(typ(L))) pari_err_TYPE("mflinear",L);
 2058   if (lg(L) != lg(F)) pari_err_DIM("mflinear");
 2059   return L;
 2060 }
 2061 GEN
 2062 mflinear(GEN F, GEN L)
 2063 {
 2064   pari_sp av = avma;
 2065   GEN G, NK, P, mf = checkMF_i(F), N = NULL, K = NULL, CHI = NULL;
 2066   long i, l;
 2067   if (mf)
 2068   {
 2069     GEN gk = MF_get_gk(mf);
 2070     F = MF_get_basis(F);
 2071     if (typ(gk) != t_INT)
 2072       return gerepilecopy(av, mflineardiv_linear(F, L, 1));
 2073     if (itou(gk) > 1 && space_is_cusp(MF_get_space(mf)))
 2074     {
 2075       L = tobasis(mf, F, L);
 2076       return gerepilecopy(av, mflinear_bhn(mf, L));
 2077     }
 2078   }
 2079   L = tobasis(mf, F, L);
 2080   if (!mflinear_strip(&F,&L)) return mftrivial();
 2081 
 2082   l = lg(F);
 2083   if (l == 2 && gequal1(gel(L,1))) return gerepilecopy(av, gel(F,1));
 2084   P = pol_x(1);
 2085   for (i = 1; i < l; i++)
 2086   {
 2087     GEN f = gel(F,i), c = gel(L,i), Ni, Ki;
 2088     if (!checkmf_i(f)) pari_err_TYPE("mflinear", f);
 2089     Ni = mf_get_gN(f); N = N? lcmii(N, Ni): Ni;
 2090     Ki = mf_get_gk(f);
 2091     if (!K) K = Ki;
 2092     else if (!gequal(K, Ki))
 2093       pari_err_TYPE("mflinear [different weights]", mkvec2(K,Ki));
 2094     P = mfsamefield(NULL, P, mf_get_field(f));
 2095     if (typ(c) == t_POLMOD && varn(gel(c,1)) == 1)
 2096       P = mfsamefield(NULL, P, gel(c,1));
 2097   }
 2098   G = znstar0(N,1);
 2099   for (i = 1; i < l; i++)
 2100   {
 2101     GEN CHI2 = mf_get_CHI(gel(F,i));
 2102     CHI2 = induce(G, CHI2);
 2103     if (!CHI) CHI = CHI2;
 2104     else if (!gequal(CHI, CHI2))
 2105       pari_err_TYPE("mflinear [different characters]", mkvec2(CHI,CHI2));
 2106   }
 2107   NK = mkgNK(N, K, CHI, P);
 2108   return gerepilecopy(av, taglinear(NK,F,L));
 2109 }
 2110 
 2111 GEN
 2112 mfshift(GEN F, long sh)
 2113 {
 2114   pari_sp av = avma;
 2115   if (!checkmf_i(F)) pari_err_TYPE("mfshift",F);
 2116   return gerepilecopy(av, tag2(t_MF_SHIFT, mf_get_NK(F), F, stoi(sh)));
 2117 }
 2118 static long
 2119 mfval(GEN F)
 2120 {
 2121   pari_sp av = avma;
 2122   long i = 0, n, sb;
 2123   GEN gk, gN;
 2124   if (!checkmf_i(F)) pari_err_TYPE("mfval", F);
 2125   gN = mf_get_gN(F);
 2126   gk = mf_get_gk(F);
 2127   sb = mfsturmNgk(itou(gN), gk);
 2128   for (n = 1; n <= sb;)
 2129   {
 2130     GEN v;
 2131     if (n > 0.5*sb) n = sb+1;
 2132     v = mfcoefs_i(F, n, 1);
 2133     for (; i <= n; i++)
 2134       if (!gequal0(gel(v, i+1))) return gc_long(av,i);
 2135     n <<= 1;
 2136   }
 2137   return gc_long(av,-1);
 2138 }
 2139 
 2140 GEN
 2141 mfdiv_val(GEN f, GEN g, long vg)
 2142 {
 2143   GEN T, N, K, NK, CHI, CHIf, CHIg;
 2144   if (vg) { f = mfshift(f,vg); g = mfshift(g,vg); }
 2145   N = lcmii(mf_get_gN(f), mf_get_gN(g));
 2146   K = gsub(mf_get_gk(f), mf_get_gk(g));
 2147   CHIf = mf_get_CHI(f);
 2148   CHIg = mf_get_CHI(g);
 2149   CHI = mfchiadjust(mfchardiv(CHIf, CHIg), K, itos(N));
 2150   T = chicompat(CHI, CHIf, CHIg);
 2151   NK = mkgNK(N, K, CHI, mfsamefield(T, mf_get_field(f), mf_get_field(g)));
 2152   return T? tag3(t_MF_DIV, NK, f, g, T): tag2(t_MF_DIV, NK, f, g);
 2153 }
 2154 GEN
 2155 mfdiv(GEN F, GEN G)
 2156 {
 2157   pari_sp av = avma;
 2158   long v = mfval(G);
 2159   if (!checkmf_i(F)) pari_err_TYPE("mfdiv", F);
 2160   if (v < 0 || (v && !gequal0(mfcoefs(F, v-1, 1))))
 2161     pari_err_DOMAIN("mfdiv", "ord(G)", ">", strtoGENstr("ord(F)"),
 2162                     mkvec2(F, G));
 2163   return gerepilecopy(av, mfdiv_val(F, G, v));
 2164 }
 2165 GEN
 2166 mfderiv(GEN F, long m)
 2167 {
 2168   pari_sp av = avma;
 2169   GEN NK, gk;
 2170   if (!checkmf_i(F)) pari_err_TYPE("mfderiv",F);
 2171   gk = gaddgs(mf_get_gk(F), 2*m);
 2172   NK = mkgNK(mf_get_gN(F), gk, mf_get_CHI(F), mf_get_field(F));
 2173   return gerepilecopy(av, tag2(t_MF_DERIV, NK, F, stoi(m)));
 2174 }
 2175 GEN
 2176 mfderivE2(GEN F, long m)
 2177 {
 2178   pari_sp av = avma;
 2179   GEN NK, gk;
 2180   if (!checkmf_i(F)) pari_err_TYPE("mfderivE2",F);
 2181   if (m < 0) pari_err_DOMAIN("mfderivE2","m","<",gen_0,stoi(m));
 2182   gk = gaddgs(mf_get_gk(F), 2*m);
 2183   NK = mkgNK(mf_get_gN(F), gk, mf_get_CHI(F), mf_get_field(F));
 2184   return gerepilecopy(av, tag2(t_MF_DERIVE2, NK, F, stoi(m)));
 2185 }
 2186 
 2187 GEN
 2188 mftwist(GEN F, GEN D)
 2189 {
 2190   pari_sp av = avma;
 2191   GEN NK, CHI, NT, Da;
 2192   long q;
 2193   if (!checkmf_i(F)) pari_err_TYPE("mftwist", F);
 2194   if (typ(D) != t_INT) pari_err_TYPE("mftwist", D);
 2195   Da = mpabs_shallow(D);
 2196   CHI = mf_get_CHI(F); q = mfcharconductor(CHI);
 2197   NT = glcm(glcm(mf_get_gN(F), mulsi(q, Da)), sqri(Da));
 2198   NK = mkgNK(NT, mf_get_gk(F), CHI, mf_get_field(F));
 2199   return gerepilecopy(av, tag2(t_MF_TWIST, NK, F, D));
 2200 }
 2201 
 2202 /***************************************************************/
 2203 /*                 Generic cache handling                      */
 2204 /***************************************************************/
 2205 enum { cache_FACT, cache_DIV, cache_H, cache_D, cache_DIH };
 2206 typedef struct {
 2207   const char *name;
 2208   GEN cache;
 2209   ulong minself, maxself;
 2210   void (*init)(long);
 2211   ulong miss, maxmiss;
 2212   long compressed;
 2213 } cache;
 2214 
 2215 static void constfact(long lim);
 2216 static void constdiv(long lim);
 2217 static void consttabh(long lim);
 2218 static void consttabdihedral(long lim);
 2219 static void constcoredisc(long lim);
 2220 static THREAD cache caches[] = {
 2221 { "Factors",  NULL,  50000,    50000, &constfact, 0, 0, 0 },
 2222 { "Divisors", NULL,  50000,    50000, &constdiv, 0, 0, 0 },
 2223 { "H",        NULL, 100000, 10000000, &consttabh, 0, 0, 1 },
 2224 { "CorediscF",NULL, 100000, 10000000, &constcoredisc, 0, 0, 0 },
 2225 { "Dihedral", NULL,   1000,     3000, &consttabdihedral, 0, 0, 0 },
 2226 };
 2227 
 2228 static void
 2229 cache_reset(long id) { caches[id].miss = caches[id].maxmiss = 0; }
 2230 static void
 2231 cache_delete(long id) { guncloneNULL(caches[id].cache); }
 2232 static void
 2233 cache_set(long id, GEN S)
 2234 {
 2235   GEN old = caches[id].cache;
 2236   caches[id].cache = gclone(S);
 2237   guncloneNULL(old);
 2238 }
 2239 
 2240 /* handle a cache miss: store stats, possibly reset table; return value
 2241  * if (now) cached; return NULL on failure. HACK: some caches contain an
 2242  * ulong where the 0 value is impossible, and return it (typecast to GEN) */
 2243 static GEN
 2244 cache_get(long id, ulong D)
 2245 {
 2246   cache *S = &caches[id];
 2247   const ulong d = S->compressed? D>>1: D;
 2248   ulong max, l;
 2249 
 2250   if (!S->cache)
 2251   {
 2252     max = maxuu(minuu(D, S->maxself), S->minself);
 2253     S->init(max);
 2254     l = lg(S->cache);
 2255   }
 2256   else
 2257   {
 2258     l = lg(S->cache);
 2259     if (l <= d)
 2260     {
 2261       if (D > S->maxmiss) S->maxmiss = D;
 2262       if (DEBUGLEVEL >= 3)
 2263         err_printf("miss in cache %s: %lu, max = %lu\n",
 2264                    S->name, D, S->maxmiss);
 2265       if (S->miss++ >= 5 && D < S->maxself)
 2266       {
 2267         max = minuu(S->maxself, (long)(S->maxmiss * 1.2));
 2268         if (max <= S->maxself)
 2269         {
 2270           if (DEBUGLEVEL >= 3)
 2271             err_printf("resetting cache %s to %lu\n", S->name, max);
 2272           S->init(max); l = lg(S->cache);
 2273         }
 2274       }
 2275     }
 2276   }
 2277   return (l <= d)? NULL: gel(S->cache, d);
 2278 }
 2279 static GEN
 2280 cache_report(long id)
 2281 {
 2282   cache *S = &caches[id];
 2283   GEN v = zerocol(5);
 2284   gel(v,1) = strtoGENstr(S->name);
 2285   if (S->cache)
 2286   {
 2287     gel(v,2) = utoi(lg(S->cache)-1);
 2288     gel(v,3) = utoi(S->miss);
 2289     gel(v,4) = utoi(S->maxmiss);
 2290     gel(v,5) = utoi(gsizebyte(S->cache));
 2291   }
 2292   return v;
 2293 }
 2294 GEN
 2295 getcache(void)
 2296 {
 2297   pari_sp av = avma;
 2298   GEN M = cgetg(6, t_MAT);
 2299   gel(M,1) = cache_report(cache_FACT);
 2300   gel(M,2) = cache_report(cache_DIV);
 2301   gel(M,3) = cache_report(cache_H);
 2302   gel(M,4) = cache_report(cache_D);
 2303   gel(M,5) = cache_report(cache_DIH);
 2304   return gerepilecopy(av, shallowtrans(M));
 2305 }
 2306 
 2307 void
 2308 pari_close_mf(void)
 2309 {
 2310   cache_delete(cache_FACT);
 2311   cache_delete(cache_DIV);
 2312   cache_delete(cache_H);
 2313   cache_delete(cache_D);
 2314   cache_delete(cache_DIH);
 2315 }
 2316 
 2317 /*************************************************************************/
 2318 /* a odd, update local cache (recycle memory) */
 2319 static GEN
 2320 update_factor_cache(long a, long lim, long *pb)
 2321 {
 2322   const long step = 16000; /* even; don't increase this: RAM cache thrashing */
 2323   if (a + 2*step > lim)
 2324     *pb = lim; /* fuse last 2 chunks */
 2325   else
 2326     *pb = a + step;
 2327   return vecfactoroddu_i(a, *pb);
 2328 }
 2329 /* assume lim < MAX_LONG/8 */
 2330 static void
 2331 constcoredisc(long lim)
 2332 {
 2333   pari_sp av2, av = avma;
 2334   GEN D = caches[cache_D].cache, CACHE = NULL;
 2335   long cachea, cacheb, N, LIM = !D ? 4 : lg(D)-1;
 2336   if (lim <= 0) lim = 5;
 2337   if (lim <= LIM) return;
 2338   cache_reset(cache_D);
 2339   D = zero_zv(lim);
 2340   av2 = avma;
 2341   cachea = cacheb = 0;
 2342   for (N = 1; N <= lim; N+=2)
 2343   { /* N odd */
 2344     long i, d, d2;
 2345     GEN F;
 2346     if (N > cacheb)
 2347     {
 2348       set_avma(av2); cachea = N;
 2349       CACHE = update_factor_cache(N, lim, &cacheb);
 2350     }
 2351     F = gel(CACHE, ((N-cachea)>>1)+1); /* factoru(N) */
 2352     D[N] = d = corediscs_fact(F); /* = 3 mod 4 or 4 mod 16 */
 2353     d2 = odd(d)? d<<3: d<<1;
 2354     for (i = 1;;)
 2355     {
 2356       if ((N << i) > lim) break;
 2357       D[N<<i] = d2; i++;
 2358       if ((N << i) > lim) break;
 2359       D[N<<i] = d; i++;
 2360     }
 2361   }
 2362   cache_set(cache_D, D);
 2363   set_avma(av);
 2364 }
 2365 
 2366 static void
 2367 constfact(long lim)
 2368 {
 2369   pari_sp av;
 2370   GEN VFACT = caches[cache_FACT].cache;
 2371   long LIM = VFACT? lg(VFACT)-1: 4;
 2372   if (lim <= 0) lim = 5;
 2373   if (lim <= LIM) return;
 2374   cache_reset(cache_FACT); av = avma;
 2375   cache_set(cache_FACT, vecfactoru_i(1,lim)); set_avma(av);
 2376 }
 2377 static void
 2378 constdiv(long lim)
 2379 {
 2380   pari_sp av;
 2381   GEN VFACT, VDIV = caches[cache_DIV].cache;
 2382   long N, LIM = VDIV? lg(VDIV)-1: 4;
 2383   if (lim <= 0) lim = 5;
 2384   if (lim <= LIM) return;
 2385   constfact(lim);
 2386   VFACT = caches[cache_FACT].cache;
 2387   cache_reset(cache_DIV); av = avma;
 2388   VDIV  = cgetg(lim+1, t_VEC);
 2389   for (N = 1; N <= lim; N++) gel(VDIV,N) = divisorsu_fact(gel(VFACT,N));
 2390   cache_set(cache_DIV, VDIV); set_avma(av);
 2391 }
 2392 
 2393 /* n > 1, D = divisors(n); sets L = 2*lambda(n), S = sigma(n) */
 2394 static void
 2395 lamsig(GEN D, long *pL, long *pS)
 2396 {
 2397   pari_sp av = avma;
 2398   long i, l = lg(D), L = 1, S = D[l-1]+1;
 2399   for (i = 2; i < l; i++) /* skip d = 1 */
 2400   {
 2401     long d = D[i], nd = D[l-i]; /* nd = n/d */
 2402     if (d < nd) { L += d; S += d + nd; }
 2403     else
 2404     {
 2405       L <<= 1; if (d == nd) { L += d; S += d; }
 2406       break;
 2407     }
 2408   }
 2409   set_avma(av); *pL = L; *pS = S;
 2410 }
 2411 /* table of 6 * Hurwitz class numbers D <= lim */
 2412 static void
 2413 consttabh(long lim)
 2414 {
 2415   pari_sp av = avma, av2;
 2416   GEN VHDH0, VDIV, CACHE = NULL;
 2417   GEN VHDH = caches[cache_H].cache;
 2418   long r, N, cachea, cacheb, lim0 = VHDH? lg(VHDH)-1: 2, LIM = lim0 << 1;
 2419 
 2420   if (lim <= 0) lim = 5;
 2421   if (lim <= LIM) return;
 2422   cache_reset(cache_H);
 2423   r = lim&3L; if (r) lim += 4-r;
 2424   cache_get(cache_DIV, lim);
 2425   VDIV = caches[cache_DIV].cache;
 2426   VHDH0 = cgetg(lim/2 + 1, t_VECSMALL);
 2427   VHDH0[1] = 2;
 2428   VHDH0[2] = 3;
 2429   for (N = 3; N <= lim0; N++) VHDH0[N] = VHDH[N];
 2430   av2 = avma;
 2431   cachea = cacheb = 0;
 2432   for (N = LIM + 3; N <= lim; N += 4)
 2433   {
 2434     long s = 0, limt = usqrt(N>>2), flsq = 0, ind, t, L, S;
 2435     GEN DN, DN2;
 2436     if (N + 2 >= lg(VDIV))
 2437     { /* use local cache */
 2438       GEN F;
 2439       if (N + 2 > cacheb)
 2440       {
 2441         set_avma(av2); cachea = N;
 2442         CACHE = update_factor_cache(N, lim+2, &cacheb);
 2443       }
 2444       F = gel(CACHE, ((N-cachea)>>1)+1); /* factoru(N) */
 2445       DN = divisorsu_fact(F);
 2446       F = gel(CACHE, ((N-cachea)>>1)+2); /* factoru(N+2) */
 2447       DN2 = divisorsu_fact(F);
 2448     }
 2449     else
 2450     { /* use global cache */
 2451       DN = gel(VDIV,N);
 2452       DN2 = gel(VDIV,N+2);
 2453     }
 2454     ind = N >> 1;
 2455     for (t = 1; t <= limt; t++)
 2456     {
 2457       ind -= (t<<2)-2; /* N/2 - 2t^2 */
 2458       if (ind) s += VHDH0[ind]; else flsq = 1;
 2459     }
 2460     lamsig(DN, &L,&S);
 2461     VHDH0[N >> 1] = 2*S - 3*L - 2*s + flsq;
 2462     s = 0; flsq = 0; limt = (usqrt(N+2) - 1) >> 1;
 2463     ind = (N+1) >> 1;
 2464     for (t = 1; t <= limt; t++)
 2465     {
 2466       ind -= t<<2; /* (N+1)/2 - 2t(t+1) */
 2467       if (ind) s += VHDH0[ind]; else flsq = 1;
 2468     }
 2469     lamsig(DN2, &L,&S);
 2470     VHDH0[(N+1) >> 1] = S - 3*(L >> 1) - s - flsq;
 2471   }
 2472   cache_set(cache_H, VHDH0); set_avma(av);
 2473 }
 2474 
 2475 /*************************************************************************/
 2476 /* Core functions using factorizations, divisors of class numbers caches */
 2477 /* TODO: myfactoru and factorization cache should be exported */
 2478 static GEN
 2479 myfactoru(long N)
 2480 {
 2481   GEN z = cache_get(cache_FACT, N);
 2482   return z? gcopy(z): factoru(N);
 2483 }
 2484 static GEN
 2485 mydivisorsu(long N)
 2486 {
 2487   GEN z = cache_get(cache_DIV, N);
 2488   return z? leafcopy(z): divisorsu(N);
 2489 }
 2490 /* write -n = Df^2, D < 0 fundamental discriminant. Return D, set f. */
 2491 static long
 2492 mycoredisc2neg(ulong n, long *pf)
 2493 {
 2494   ulong m, D = (ulong)cache_get(cache_D, n);
 2495   if (D) { *pf = usqrt(n/D); return -(long)D; }
 2496   m = mycore(n, pf);
 2497   if ((m&3) != 3) { m <<= 2; *pf >>= 1; }
 2498   return (long)-m;
 2499 }
 2500 /* write n = Df^2, D > 0 fundamental discriminant. Return D, set f. */
 2501 static long
 2502 mycoredisc2pos(ulong n, long *pf)
 2503 {
 2504   ulong m = mycore(n, pf);
 2505   if ((m&3) != 1) { m <<= 2; *pf >>= 1; }
 2506   return (long)m;
 2507 }
 2508 
 2509 /* 1+p+...+p^e, e >= 1 */
 2510 static ulong
 2511 usumpow(ulong p, long e)
 2512 {
 2513   ulong q = 1+p;
 2514   long i;
 2515   for (i = 1; i < e; i++) q = p*q + 1;
 2516   return q;
 2517 }
 2518 /* Hurwitz(D0 F^2)/ Hurwitz(D0)
 2519  * = \sum_{f|F}  f \prod_{p|f} (1-kro(D0/p)/p)
 2520  * = \prod_{p^e || F} (1 + (p^e-1) / (p-1) * (p-kro(D0/p))) */
 2521 static long
 2522 get_sh(long F, long D0)
 2523 {
 2524   GEN fa = myfactoru(F), P = gel(fa,1), E = gel(fa,2);
 2525   long i, l = lg(P), t = 1;
 2526   for (i = 1; i < l; i++)
 2527   {
 2528     long p = P[i], e = E[i], s = kross(D0,p);
 2529     if (e == 1) { t *= 1 + p - s; continue; }
 2530     if (s == 1) { t *= upowuu(p,e); continue; }
 2531     t *= 1 + usumpow(p,e-1)*(p-s);
 2532   }
 2533   return t;
 2534 }
 2535 /* d > 0, d = 0,3 (mod 4). Return 6*hclassno(d); -d must be fundamental
 2536  * Faster than quadclassunit up to 5*10^5 or so */
 2537 static ulong
 2538 hclassno6u_count(ulong d)
 2539 {
 2540   ulong a, b, b2, h = 0;
 2541   int f = 0;
 2542 
 2543   if (d > 500000)
 2544     return 6 * itou(gel(quadclassunit0(utoineg(d), 0, NULL, 0), 1));
 2545 
 2546   /* this part would work with -d non fundamental */
 2547   b = d&1; b2 = (1+d)>>2;
 2548   if (!b)
 2549   {
 2550     for (a=1; a*a<b2; a++)
 2551       if (b2%a == 0) h++;
 2552     f = (a*a==b2); b=2; b2=(4+d)>>2;
 2553   }
 2554   while (b2*3 < d)
 2555   {
 2556     if (b2%b == 0) h++;
 2557     for (a=b+1; a*a < b2; a++)
 2558       if (b2%a == 0) h += 2;
 2559     if (a*a == b2) h++;
 2560     b += 2; b2 = (b*b+d)>>2;
 2561   }
 2562   if (b2*3 == d) return 6*h+2;
 2563   if (f) return 6*h+3;
 2564   return 6*h;
 2565 }
 2566 /* D > 0; 6 * hclassno(D), using D = D0*F^2 */
 2567 static long
 2568 hclassno6u_2(ulong D, long D0, long F)
 2569 {
 2570   long h;
 2571   if (F == 1) h = hclassno6u_count(D);
 2572   else
 2573   { /* second chance */
 2574     h = (ulong)cache_get(cache_H, -D0);
 2575     if (!h) h = hclassno6u_count(-D0);
 2576     h *= get_sh(F,D0);
 2577   }
 2578   return h;
 2579 }
 2580 /* D > 0; 6 * hclassno(D) (6*Hurwitz). Beware, cached value for D (=0,3 mod 4)
 2581  * is stored at D>>1 */
 2582 ulong
 2583 hclassno6u(ulong D)
 2584 {
 2585   ulong z = (ulong)cache_get(cache_H, D);
 2586   long D0, F;
 2587   if (z) return z;
 2588   D0 = mycoredisc2neg(D, &F);
 2589   return hclassno6u_2(D,D0,F);
 2590 }
 2591 /* same, where the decomposition D = D0*F^2 is already known */
 2592 static ulong
 2593 hclassno6u_i(ulong D, long D0, long F)
 2594 {
 2595   ulong z = (ulong)cache_get(cache_H, D);
 2596   if (z) return z;
 2597   return hclassno6u_2(D,D0,F);
 2598 }
 2599 
 2600 #if 0
 2601 /* D > 0, return h(-D) [ordinary class number].
 2602  * Assume consttabh(D or more) was previously called */
 2603 static long
 2604 hfromH(long D)
 2605 {
 2606   pari_sp ltop = avma;
 2607   GEN m, d, fa = myfactoru(D), P = gel(fa,1), E = gel(fa,2);
 2608   GEN VH = caches[cache_H].cache;
 2609   long i, nd, S, l = lg(P);
 2610 
 2611   /* n = d[i] loops through squarefree divisors of f, where f^2 = largest square
 2612    * divisor of N = |D|; m[i] = moebius(n) */
 2613   nd = 1 << (l-1);
 2614   d = cgetg(nd+1, t_VECSMALL);
 2615   m = cgetg(nd+1, t_VECSMALL);
 2616   d[1] = 1; S = VH[D >> 1]; /* 6 hclassno(-D) */
 2617   m[1] = 1; nd = 1;
 2618   i = 1;
 2619   if (P[1] == 2 && E[1] <= 3) /* need D/n^2 to be a discriminant */
 2620   { if (odd(E[1]) || (E[1] == 2 && (D & 15) == 4)) i = 2; }
 2621   for (; i<l; i++)
 2622   {
 2623     long j, p = P[i];
 2624     if (E[i] == 1) continue;
 2625     for (j=1; j<=nd; j++)
 2626     {
 2627       long n, s, hn;
 2628       d[nd+j] = n = d[j] * p;
 2629       m[nd+j] = s = - m[j]; /* moebius(n) */
 2630       hn = VH[(D/(n*n)) >> 1]; /* 6 hclassno(-D/n^2) */
 2631       if (s > 0) S += hn; else S -= hn;
 2632     }
 2633     nd <<= 1;
 2634   }
 2635   return gc_long(ltop, S/6);
 2636 }
 2637 #endif
 2638 /* D < -4 fundamental, h(D), ordinary class number */
 2639 static long
 2640 myh(long D)
 2641 {
 2642   ulong z = (ulong)cache_get(cache_H, -D);
 2643   if (z) return z/6; /* should be hfromH(-D) if D nonfundamental */
 2644   return itou(quadclassno(stoi(D)));
 2645 }
 2646 
 2647 /*************************************************************************/
 2648 /*                          TRACE FORMULAS                               */
 2649 /* CHIP primitive, initialize for t_POLMOD output */
 2650 static GEN
 2651 mfcharinit(GEN CHIP)
 2652 {
 2653   long n, o, l, vt, N = mfcharmodulus(CHIP);
 2654   GEN c, v, V, G, Pn;
 2655   if (N == 1) return mkvec2(mkvec(gen_1), pol_x(0));
 2656   G = gel(CHIP,1);
 2657   v = ncharvecexpo(G, znconrey_normalized(G, gel(CHIP,2)));
 2658   l = lg(v); V = cgetg(l, t_VEC);
 2659   o = mfcharorder(CHIP);
 2660   Pn = mfcharpol(CHIP); vt = varn(Pn);
 2661   if (o <= 2)
 2662   {
 2663     for (n = 1; n < l; n++)
 2664     {
 2665       if (v[n] < 0) c = gen_0; else c = v[n]? gen_m1: gen_1;
 2666       gel(V,n) = c;
 2667     }
 2668   }
 2669   else
 2670   {
 2671     for (n = 1; n < l; n++)
 2672     {
 2673       if (v[n] < 0) c = gen_0;
 2674       else
 2675       {
 2676         c = Qab_zeta(v[n], o, vt);
 2677         if (typ(c) == t_POL && lg(c) >= lg(Pn)) c = RgX_rem(c, Pn);
 2678       }
 2679       gel(V,n) = c;
 2680     }
 2681   }
 2682   return mkvec2(V, Pn);
 2683 }
 2684 static GEN
 2685 vchip_lift(GEN VCHI, long x, GEN C)
 2686 {
 2687   GEN V = gel(VCHI,1);
 2688   long F = lg(V)-1;
 2689   if (F == 1) return C;
 2690   x %= F;
 2691   if (!x) return C;
 2692   if (x <= 0) x += F;
 2693   return gmul(C, gel(V, x));
 2694 }
 2695 static long
 2696 vchip_FC(GEN VCHI) { return lg(gel(VCHI,1))-1; }
 2697 static GEN
 2698 vchip_mod(GEN VCHI, GEN S)
 2699 { return (typ(S) == t_POL)? RgX_rem(S, gel(VCHI,2)): S; }
 2700 static GEN
 2701 vchip_polmod(GEN VCHI, GEN S)
 2702 { return (typ(S) == t_POL)? mkpolmod(S, gel(VCHI,2)): S; }
 2703 
 2704 /* ceil(m/d) */
 2705 static long
 2706 ceildiv(long m, long d)
 2707 {
 2708   long q;
 2709   if (!m) return 0;
 2710   q = m/d; return m%d? q+1: q;
 2711 }
 2712 
 2713 /* contribution of scalar matrices in dimension formula */
 2714 static GEN
 2715 A1(long N, long k)
 2716 { return sstoQ(mypsiu(N)*(k-1), 12); }
 2717 static long
 2718 ceilA1(long N, long k)
 2719 { return ceildiv(mypsiu(N) * (k-1), 12); }
 2720 
 2721 /* sturm bound, slightly larger than dimension */
 2722 long
 2723 mfsturmNk(long N, long k) { return (mypsiu(N) * k) / 12; }
 2724 long
 2725 mfsturmNgk(long N, GEN k)
 2726 {
 2727   long n,d; Qtoss(k,&n,&d);
 2728   return 1 + (mypsiu(N)*n)/(d == 1? 12: 24);
 2729 }
 2730 static long
 2731 mfsturmmf(GEN F) { return mfsturmNgk(mf_get_N(F), mf_get_gk(F)); }
 2732 
 2733 /* List of all solutions of x^2 + x + 1 = 0 modulo N, x modulo N */
 2734 static GEN
 2735 sqrtm3modN(long N)
 2736 {
 2737   pari_sp av;
 2738   GEN fa, P, E, B, mB, A, Q, T, R, v, gen_m3;
 2739   long l, i, n, ct, fl3 = 0, Ninit;
 2740   if (!odd(N) || (N%9) == 0) return cgetg(1,t_VECSMALL);
 2741   Ninit = N;
 2742   if ((N%3) == 0) { N /= 3; fl3 = 1; }
 2743   fa = myfactoru(N); P = gel(fa, 1); E = gel(fa, 2);
 2744   l = lg(P);
 2745   for (i = 1; i < l; i++)
 2746     if ((P[i]%3) == 2) return cgetg(1,t_VECSMALL);
 2747   A = cgetg(l, t_VECSMALL);
 2748   B = cgetg(l, t_VECSMALL);
 2749   mB= cgetg(l, t_VECSMALL);
 2750   Q = cgetg(l, t_VECSMALL); gen_m3 = utoineg(3);
 2751   for (i = 1; i < l; i++)
 2752   {
 2753     long p = P[i], e = E[i];
 2754     Q[i] = upowuu(p,e);
 2755     B[i] = itou( Zp_sqrt(gen_m3, utoipos(p), e) );
 2756     mB[i]= Q[i] - B[i];
 2757   }
 2758   ct = 1 << (l-1);
 2759   T = ZV_producttree(Q);
 2760   R = ZV_chinesetree(Q,T);
 2761   v = cgetg(ct+1, t_VECSMALL);
 2762   av = avma;
 2763   for (n = 1; n <= ct; n++)
 2764   {
 2765     long m = n-1, r;
 2766     for (i = 1; i < l; i++)
 2767     {
 2768       A[i] = (m&1L)? mB[i]: B[i];
 2769       m >>= 1;
 2770     }
 2771     r = itou( ZV_chinese_tree(A, Q, T, R) );
 2772     if (fl3) while (r%3) r += N;
 2773     set_avma(av); v[n] = odd(r) ? (r-1) >> 1 : (r+Ninit-1) >> 1;
 2774   }
 2775   return v;
 2776 }
 2777 
 2778 /* number of elliptic points of order 3 in X0(N) */
 2779 static long
 2780 nu3(long N)
 2781 {
 2782   long i, l;
 2783   GEN P;
 2784   if (!odd(N) || (N%9) == 0) return 0;
 2785   if ((N%3) == 0) N /= 3;
 2786   P = gel(myfactoru(N), 1); l = lg(P);
 2787   for (i = 1; i < l; i++) if ((P[i]%3) == 2) return 0;
 2788   return 1L<<(l-1);
 2789 }
 2790 /* number of elliptic points of order 2 in X0(N) */
 2791 static long
 2792 nu2(long N)
 2793 {
 2794   long i, l;
 2795   GEN P;
 2796   if ((N&3L) == 0) return 0;
 2797   if (!odd(N)) N >>= 1;
 2798   P = gel(myfactoru(N), 1); l = lg(P);
 2799   for (i = 1; i < l; i++) if ((P[i]&3L) == 3) return 0;
 2800   return 1L<<(l-1);
 2801 }
 2802 
 2803 /* contribution of elliptic matrices of order 3 in dimension formula
 2804  * Only depends on CHIP the primitive char attached to CHI */
 2805 static GEN
 2806 A21(long N, long k, GEN CHI)
 2807 {
 2808   GEN res, G, chi, o;
 2809   long a21, i, limx, S;
 2810   if ((N&1L) == 0) return gen_0;
 2811   a21 = k%3 - 1;
 2812   if (!a21) return gen_0;
 2813   if (N <= 3) return sstoQ(a21, 3);
 2814   if (!CHI) return sstoQ(nu3(N) * a21, 3);
 2815   res = sqrtm3modN(N); limx = (N - 1) >> 1;
 2816   G = gel(CHI,1); chi = gel(CHI,2);
 2817   o = gmfcharorder(CHI);
 2818   for (S = 0, i = 1; i < lg(res); i++)
 2819   { /* (x,N) = 1; S += chi(x) + chi(x^2) */
 2820     long x = res[i];
 2821     if (x <= limx)
 2822     { /* CHI(x)=e(c/o), 3rd-root of 1 */
 2823       GEN c = znchareval(G, chi, utoi(x), o);
 2824       if (!signe(c)) S += 2; else S--;
 2825     }
 2826   }
 2827   return sstoQ(a21 * S, 3);
 2828 }
 2829 
 2830 /* List of all square roots of -1 modulo N */
 2831 static GEN
 2832 sqrtm1modN(long N)
 2833 {
 2834   pari_sp av;
 2835   GEN fa, P, E, B, mB, A, Q, T, R, v;
 2836   long l, i, n, ct, fleven = 0;
 2837   if ((N&3L) == 0) return cgetg(1,t_VECSMALL);
 2838   if ((N&1L) == 0) { N >>= 1; fleven = 1; }
 2839   fa = myfactoru(N); P = gel(fa,1); E = gel(fa,2);
 2840   l = lg(P);
 2841   for (i = 1; i < l; i++)
 2842     if ((P[i]&3L) == 3) return cgetg(1,t_VECSMALL);
 2843   A = cgetg(l, t_VECSMALL);
 2844   B = cgetg(l, t_VECSMALL);
 2845   mB= cgetg(l, t_VECSMALL);
 2846   Q = cgetg(l, t_VECSMALL);
 2847   for (i = 1; i < l; i++)
 2848   {
 2849     long p = P[i], e = E[i];
 2850     Q[i] = upowuu(p,e);
 2851     B[i] = itou( Zp_sqrt(gen_m1, utoipos(p), e) );
 2852     mB[i]= Q[i] - B[i];
 2853   }
 2854   ct = 1 << (l-1);
 2855   T = ZV_producttree(Q);
 2856   R = ZV_chinesetree(Q,T);
 2857   v = cgetg(ct+1, t_VECSMALL);
 2858   av = avma;
 2859   for (n = 1; n <= ct; n++)
 2860   {
 2861     long m = n-1, r;
 2862     for (i = 1; i < l; i++)
 2863     {
 2864       A[i] = (m&1L)? mB[i]: B[i];
 2865       m >>= 1;
 2866     }
 2867     r = itou( ZV_chinese_tree(A, Q, T, R) );
 2868     if (fleven && !odd(r)) r += N;
 2869     set_avma(av); v[n] = r;
 2870   }
 2871   return v;
 2872 }
 2873 
 2874 /* contribution of elliptic matrices of order 4 in dimension formula.
 2875  * Only depends on CHIP the primitive char attached to CHI */
 2876 static GEN
 2877 A22(long N, long k, GEN CHI)
 2878 {
 2879   GEN G, chi, o, res;
 2880   long S, a22, i, limx, o2;
 2881   if ((N&3L) == 0) return gen_0;
 2882   a22 = (k & 3L) - 1; /* (k % 4) - 1 */
 2883   if (!a22) return gen_0;
 2884   if (N <= 2) return sstoQ(a22, 4);
 2885   if (!CHI) return sstoQ(nu2(N)*a22, 4);
 2886   if (mfcharparity(CHI) == -1) return gen_0;
 2887   res = sqrtm1modN(N); limx = (N - 1) >> 1;
 2888   G = gel(CHI,1); chi = gel(CHI,2);
 2889   o = gmfcharorder(CHI);
 2890   o2 = itou(o)>>1;
 2891   for (S = 0, i = 1; i < lg(res); i++)
 2892   { /* (x,N) = 1, S += real(chi(x)) */
 2893     long x = res[i];
 2894     if (x <= limx)
 2895     { /* CHI(x)=e(c/o), 4th-root of 1 */
 2896       long c = itou( znchareval(G, chi, utoi(x), o) );
 2897       if (!c) S++; else if (c == o2) S--;
 2898     }
 2899   }
 2900   return sstoQ(a22 * S, 2);
 2901 }
 2902 
 2903 /* sumdiv(N,d,eulerphi(gcd(d,N/d))) */
 2904 static long
 2905 nuinf(long N)
 2906 {
 2907   GEN fa = myfactoru(N), P = gel(fa,1), E = gel(fa,2);
 2908   long i, t = 1, l = lg(P);
 2909   for (i=1; i<l; i++)
 2910   {
 2911     long p = P[i], e = E[i];
 2912     if (odd(e))
 2913       t *= upowuu(p,e>>1) << 1;
 2914     else
 2915       t *= upowuu(p,(e>>1)-1) * (p+1);
 2916   }
 2917   return t;
 2918 }
 2919 
 2920 /* contribution of hyperbolic matrices in dimension formula */
 2921 static GEN
 2922 A3(long N, long FC)
 2923 {
 2924   long i, S, NF, l;
 2925   GEN D;
 2926   if (FC == 1) return sstoQ(nuinf(N),2);
 2927   D = mydivisorsu(N); l = lg(D);
 2928   S = 0; NF = N/FC;
 2929   for (i = 1; i < l; i++)
 2930   {
 2931     long g = ugcd(D[i], D[l-i]);
 2932     if (NF%g == 0) S += myeulerphiu(g);
 2933   }
 2934   return sstoQ(S, 2);
 2935 }
 2936 
 2937 /* special contribution in weight 2 in dimension formula */
 2938 static long
 2939 A4(long k, long FC)
 2940 { return (k==2 && FC==1)? 1: 0; }
 2941 /* gcd(x,N) */
 2942 static long
 2943 myugcd(GEN GCD, ulong x)
 2944 {
 2945   ulong N = lg(GCD)-1;
 2946   if (x >= N) x %= N;
 2947   return GCD[x+1];
 2948 }
 2949 /* 1_{gcd(x,N) = 1} * chi(x), return NULL if 0 */
 2950 static GEN
 2951 mychicgcd(GEN GCD, GEN VCHI, long x)
 2952 {
 2953   long N = lg(GCD)-1;
 2954   if (N == 1) return gen_1;
 2955   x = umodsu(x, N);
 2956   if (GCD[x+1] != 1) return NULL;
 2957   x %= vchip_FC(VCHI); if (!x) return gen_1;
 2958   return gel(gel(VCHI,1), x);
 2959 }
 2960 
 2961 /* contribution of scalar matrices to trace formula */
 2962 static GEN
 2963 TA1(long N, long k, GEN VCHI, GEN GCD, long n)
 2964 {
 2965   GEN S;
 2966   ulong m;
 2967   if (!uissquareall(n, &m)) return gen_0;
 2968   if (m == 1) return A1(N,k); /* common */
 2969   S = mychicgcd(GCD, VCHI, m);
 2970   return S? gmul(gmul(powuu(m, k-2), A1(N,k)), S): gen_0;
 2971 }
 2972 
 2973 /* All square roots modulo 4N, x modulo 2N, precomputed to accelerate TA2 */
 2974 static GEN
 2975 mksqr(long N)
 2976 {
 2977   pari_sp av = avma;
 2978   long x, N2 = N << 1, N4 = N << 2;
 2979   GEN v = const_vec(N2, cgetg(1, t_VECSMALL));
 2980   gel(v, N2) = mkvecsmall(0); /* x = 0 */
 2981   for (x = 1; x <= N; x++)
 2982   {
 2983     long r = (((x*x - 1)%N4) >> 1) + 1;
 2984     gel(v,r) = vecsmall_append(gel(v,r), x);
 2985   }
 2986   return gerepilecopy(av, v);
 2987 }
 2988 
 2989 static GEN
 2990 mkgcd(long N)
 2991 {
 2992   GEN GCD, d;
 2993   long i, N2;
 2994   if (N == 1) return mkvecsmall(N);
 2995   GCD = cgetg(N + 1, t_VECSMALL);
 2996   d = GCD+1; /* GCD[i+1] = d[i] = gcd(i,N) = gcd(N-i,N), i = 0..N-1 */
 2997   d[0] = N; d[1] = d[N-1] = 1; N2 = N>>1;
 2998   for (i = 2; i <= N2; i++) d[i] = d[N-i] = ugcd(N, i);
 2999   return GCD;
 3000 }
 3001 
 3002 /* Table of \sum_{x^2-tx+n=0 mod Ng}chi(x) for all g dividing gcd(N,F),
 3003  * F^2 largest such that (t^2-4n)/F^2=0 or 1 mod 4; t >= 0 */
 3004 static GEN
 3005 mutglistall(long t, long N, long NF, GEN VCHI, long n, GEN MUP, GEN li, GEN GCD)
 3006 {
 3007   long i, lx = lg(li);
 3008   GEN DNF = mydivisorsu(NF), v = zerovec(NF);
 3009   long j, g, lDNF = lg(DNF);
 3010   for (i = 1; i < lx; i++)
 3011   {
 3012     long x = (li[i] + t) >> 1, y, lD;
 3013     GEN D, c = mychicgcd(GCD, VCHI, x);
 3014     if (li[i] && li[i] != N)
 3015     {
 3016       GEN c2 = mychicgcd(GCD, VCHI, t - x);
 3017       if (c2) c = c? gadd(c, c2): c2;
 3018     }
 3019     if (!c) continue;
 3020     y = (x*(x - t) + n) / N; /* exact division */
 3021     D = mydivisorsu(ugcd(labs(y), NF)); lD = lg(D);
 3022     for (j=1; j < lD; j++) { g = D[j]; gel(v,g) = gadd(gel(v,g), c); }
 3023   }
 3024   /* j = 1 corresponds to g = 1, and MUP[1] = 1 */
 3025   for (j=2; j < lDNF; j++) { g = DNF[j]; gel(v,g) = gmulsg(MUP[g], gel(v,g)); }
 3026   return v;
 3027 }
 3028 
 3029 /* special case (N,F) = 1: easier */
 3030 static GEN
 3031 mutg1(long t, long N, GEN VCHI, GEN li, GEN GCD)
 3032 { /* (N,F) = 1 */
 3033   GEN S = NULL;
 3034   long i, lx = lg(li);
 3035   for (i = 1; i < lx; i++)
 3036   {
 3037     long x = (li[i] + t) >> 1;
 3038     GEN c = mychicgcd(GCD, VCHI, x);
 3039     if (c) S = S? gadd(S, c): c;
 3040     if (li[i] && li[i] != N)
 3041     {
 3042       c = mychicgcd(GCD, VCHI, t - x);
 3043       if (c) S = S? gadd(S, c): c;
 3044     }
 3045     if (S && !signe(S)) S = NULL; /* strive hard to add gen_0 */
 3046   }
 3047   return S; /* single value */
 3048 }
 3049 
 3050 /* Gegenbauer pol; n > 2, P = \sum_{0<=j<=n/2} (-1)^j (n-j)!/j!(n-2*j)! X^j */
 3051 static GEN
 3052 mfrhopol(long n)
 3053 {
 3054 #ifdef LONG_IS_64BIT
 3055   const long M = 2642249;
 3056 #else
 3057   const long M = 1629;
 3058 #endif
 3059   long j, d = n >> 1; /* >= 1 */
 3060   GEN P = cgetg(d + 3, t_POL);
 3061 
 3062   if (n > M) pari_err_IMPL("mfrhopol for large weight"); /* avoid overflow */
 3063   P[1] = evalvarn(0)|evalsigne(1);
 3064   gel(P,2) = gen_1;
 3065   gel(P,3) = utoineg(n-1); /* j = 1 */
 3066   if (d > 1) gel(P,4) = utoipos(((n-3)*(n-2)) >> 1); /* j = 2 */
 3067   if (d > 2) gel(P,5) = utoineg(((n-5)*(n-4)*(n-3)) / 6); /* j = 3 */
 3068   for (j = 4; j <= d; j++)
 3069     gel(P,j+2) = divis(mulis(gel(P,j+1), (n-2*j+1)*(n-2*j+2)), (n-j+1)*(-j));
 3070   return P;
 3071 }
 3072 
 3073 /* polrecip(Q)(t2), assume Q(0) = 1 */
 3074 static GEN
 3075 ZXrecip_u_eval(GEN Q, ulong t2)
 3076 {
 3077   GEN T = addiu(gel(Q,3), t2);
 3078   long l = lg(Q), j;
 3079   for (j = 4; j < l; j++) T = addii(gel(Q,j), mului(t2, T));
 3080   return T;
 3081 }
 3082 /* return sh * sqrt(n)^nu * G_nu(t/(2*sqrt(n))) for t != 0
 3083  * else (sh/2) * sqrt(n)^nu * G_nu(0) [ implies nu is even ]
 3084  * G_nu(z) = \sum_{0<=j<=nu/2} (-1)^j (nu-j)!/j!(nu-2*j)! * (2z)^(nu-2*j)) */
 3085 static GEN
 3086 mfrhopowsimp(GEN Q, GEN sh, long nu, long t, long t2, long n)
 3087 {
 3088   GEN T;
 3089   switch (nu)
 3090   {
 3091     case 0: return t? sh: gmul2n(sh,-1);
 3092     case 1: return gmulsg(t, sh);
 3093     case 2: return t? gmulsg(t2 - n, sh): gmul(gmul2n(stoi(-n), -1), sh);
 3094     case 3: return gmul(mulss(t, t2 - 2*n), sh);
 3095     default:
 3096       if (!t) return gmul(gmul2n(gel(Q, lg(Q) - 1), -1), sh);
 3097       T = ZXrecip_u_eval(Q, t2); if (odd(nu)) T = mulsi(t, T);
 3098       return gmul(T, sh);
 3099   }
 3100 }
 3101 
 3102 /* contribution of elliptic matrices to trace formula */
 3103 static GEN
 3104 TA2(long N, long k, GEN VCHI, long n, GEN SQRTS, GEN MUP, GEN GCD)
 3105 {
 3106   const long n4 = n << 2, N4 = N << 2, nu = k - 2;
 3107   const long st = (!odd(N) && odd(n)) ? 2 : 1;
 3108   long limt, t;
 3109   GEN S, Q;
 3110 
 3111   limt = usqrt(n4);
 3112   if (limt*limt == n4) limt--;
 3113   Q = nu > 3 ? ZX_z_unscale(mfrhopol(nu), n) : NULL;
 3114   S = gen_0;
 3115   for (t = odd(k)? st: 0; t <= limt; t += st) /* t^2 < 4n */
 3116   {
 3117     pari_sp av = avma;
 3118     long t2 = t*t, D = n4 - t2, F, D0, NF;
 3119     GEN sh, li;
 3120 
 3121     li = gel(SQRTS, (umodsu(-D - 1, N4) >> 1) + 1);
 3122     if (lg(li) == 1) continue;
 3123     D0 = mycoredisc2neg(D, &F);
 3124     NF = myugcd(GCD, F);
 3125     if (NF == 1)
 3126     { /* (N,F) = 1 => single value in mutglistall */
 3127       GEN mut = mutg1(t, N, VCHI, li, GCD);
 3128       if (!mut) { set_avma(av); continue; }
 3129       sh = gmul(sstoQ(hclassno6u_i(D,D0,F),6), mut);
 3130     }
 3131     else
 3132     {
 3133       GEN v = mutglistall(t, N, NF, VCHI, n, MUP, li, GCD);
 3134       GEN DF = mydivisorsu(F);
 3135       long i, lDF = lg(DF);
 3136       sh = gen_0;
 3137       for (i = 1; i < lDF; i++)
 3138       {
 3139         long Ff, f = DF[i], g = myugcd(GCD, f);
 3140         GEN mut = gel(v, g);
 3141         if (gequal0(mut)) continue;
 3142         Ff = DF[lDF-i]; /* F/f */
 3143         if (Ff == 1) sh = gadd(sh, mut);
 3144         else
 3145         {
 3146           GEN P = gel(myfactoru(Ff), 1);
 3147           long j, lP = lg(P);
 3148           for (j = 1; j < lP; j++) { long p = P[j]; Ff -= kross(D0, p)*Ff/p; }
 3149           sh = gadd(sh, gmulsg(Ff, mut));
 3150         }
 3151       }
 3152       if (gequal0(sh)) { set_avma(av); continue; }
 3153       if (D0 == -3) sh = gdivgs(sh, 3);
 3154       else if (D0 == -4) sh = gdivgs(sh, 2);
 3155       else sh = gmulgs(sh, myh(D0));
 3156     }
 3157     S = gerepileupto(av, gadd(S, mfrhopowsimp(Q,sh,nu,t,t2,n)));
 3158   }
 3159   return S;
 3160 }
 3161 
 3162 /* compute global auxiliary data for TA3 */
 3163 static GEN
 3164 mkbez(long N, long FC)
 3165 {
 3166   long ct, i, NF = N/FC;
 3167   GEN w, D = mydivisorsu(N);
 3168   long l = lg(D);
 3169 
 3170   w = cgetg(l, t_VEC);
 3171   for (i = ct = 1; i < l; i++)
 3172   {
 3173     long u, v, h, c = D[i], Nc = D[l-i];
 3174     if (c > Nc) break;
 3175     h = cbezout(c, Nc, &u, &v);
 3176     if (h == 1) /* shortcut */
 3177       gel(w, ct++) = mkvecsmall4(1,u*c,1,i);
 3178     else if (!(NF%h))
 3179       gel(w, ct++) = mkvecsmall4(h,u*(c/h),myeulerphiu(h),i);
 3180   }
 3181   setlg(w,ct); stackdummy((pari_sp)(w+ct),(pari_sp)(w+l));
 3182   return w;
 3183 }
 3184 
 3185 /* contribution of hyperbolic matrices to trace formula, d * nd = n,
 3186  * DN = divisorsu(N) */
 3187 static GEN
 3188 auxsum(GEN VCHI, GEN GCD, long d, long nd, GEN DN, GEN BEZ)
 3189 {
 3190   GEN S = gen_0;
 3191   long ct, g = nd - d, lDN = lg(DN), lBEZ = lg(BEZ);
 3192   for (ct = 1; ct < lBEZ; ct++)
 3193   {
 3194     GEN y, B = gel(BEZ, ct);
 3195     long ic, c, Nc, uch, h = B[1];
 3196     if (g%h) continue;
 3197     uch = B[2];
 3198     ic  = B[4];
 3199     c = DN[ic];
 3200     Nc= DN[lDN - ic]; /* Nc = N/c */
 3201     if (ugcd(Nc, nd) == 1)
 3202       y = mychicgcd(GCD, VCHI, d + uch*g); /* 0 if (c,d) > 1 */
 3203     else
 3204       y = NULL;
 3205     if (c != Nc && ugcd(Nc, d) == 1)
 3206     {
 3207       GEN y2 = mychicgcd(GCD, VCHI, nd - uch*g); /* 0 if (c,nd) > 1 */
 3208       if (y2) y = y? gadd(y, y2): y2;
 3209     }
 3210     if (y) S = gadd(S, gmulsg(B[3], y));
 3211   }
 3212   return S;
 3213 }
 3214 
 3215 static GEN
 3216 TA3(long N, long k, GEN VCHI, GEN GCD, GEN Dn, GEN BEZ)
 3217 {
 3218   GEN S = gen_0, DN = mydivisorsu(N);
 3219   long i, l = lg(Dn);
 3220   for (i = 1; i < l; i++)
 3221   {
 3222     long d = Dn[i], nd = Dn[l-i]; /* = n/d */
 3223     GEN t, u;
 3224     if (d > nd) break;
 3225     t = auxsum(VCHI, GCD, d, nd, DN, BEZ);
 3226     if (isintzero(t)) continue;
 3227     u = powuu(d,k-1); if (d == nd) u = gmul2n(u,-1);
 3228     S = gadd(S, gmul(u,t));
 3229   }
 3230   return S;
 3231 }
 3232 
 3233 /* special contribution in weight 2 in trace formula */
 3234 static long
 3235 TA4(long k, GEN VCHIP, GEN Dn, GEN GCD)
 3236 {
 3237   long i, l, S;
 3238   if (k != 2 || vchip_FC(VCHIP) != 1) return 0;
 3239   l = lg(Dn); S = 0;
 3240   for (i = 1; i < l; i++)
 3241   {
 3242     long d = Dn[i]; /* gcd(N,n/d) == 1? */
 3243     if (myugcd(GCD, Dn[l-i]) == 1) S += d;
 3244   }
 3245   return S;
 3246 }
 3247 
 3248 /* precomputation of products occurring im mutg, again to accelerate TA2 */
 3249 static GEN
 3250 mkmup(long N)
 3251 {
 3252   GEN fa = myfactoru(N), P = gel(fa,1), D = divisorsu_fact(fa);
 3253   long i, lP = lg(P), lD = lg(D);
 3254   GEN MUP = zero_zv(N);
 3255   MUP[1] = 1;
 3256   for (i = 2; i < lD; i++)
 3257   {
 3258     long j, g = D[i], Ng = D[lD-i]; /*  N/g */
 3259     for (j = 1; j < lP; j++) { long p = P[j]; if (Ng%p) g += g/p; }
 3260     MUP[D[i]] = g;
 3261   }
 3262   return MUP;
 3263 }
 3264 
 3265 /* quadratic nonresidues mod p; p odd prime, p^2 fits in a long */
 3266 static GEN
 3267 non_residues(long p)
 3268 {
 3269   long i, j, p2 = p >> 1;
 3270   GEN v = cgetg(p2+1, t_VECSMALL), w = const_vecsmall(p-1, 1);
 3271   for (i = 2; i <= p2; i++) w[(i*i) % p] = 0; /* no need to check 1 */
 3272   for (i = 2, j = 1; i < p; i++) if (w[i]) v[j++] = i;
 3273   return v;
 3274 }
 3275 
 3276 /* CHIP primitive. Return t_VECSMALL v of length q such that
 3277  * Tr^new_{N,CHIP}(n) = 0 whenever v[(n%q) + 1] is nonzero */
 3278 static GEN
 3279 mfnewzerodata(long N, GEN CHIP)
 3280 {
 3281   GEN V, M, L, faN = myfactoru(N), PN = gel(faN,1), EN = gel(faN,2);
 3282   GEN G = gel(CHIP,1), chi = gel(CHIP,2);
 3283   GEN fa = znstar_get_faN(G), P = ZV_to_zv(gel(fa,1)), E = gel(fa,2);
 3284   long i, mod, j = 1, l = lg(PN);
 3285 
 3286   M = cgetg(l, t_VECSMALL); M[1] = 0;
 3287   V = cgetg(l, t_VEC);
 3288   /* Tr^new(n) = 0 if (n mod M[i]) in V[i]  */
 3289   if ((N & 3) == 0)
 3290   {
 3291     long e = EN[1];
 3292     long c = (lg(P) > 1 && P[1] == 2)? E[1]: 0; /* c = v_2(FC) */
 3293     /* e >= 2 */
 3294     if (c == e-1) return NULL; /* Tr^new = 0 */
 3295     if (c == e)
 3296     {
 3297       if (e == 2)
 3298       { /* sc: -4 */
 3299         gel(V,1) = mkvecsmall(3);
 3300         M[1] = 4;
 3301       }
 3302       else if (e == 3)
 3303       { /* sc: -8 (CHI_2(-1)=-1<=>chi[1]=1) and 8 (CHI_2(-1)=1 <=> chi[1]=0) */
 3304         long t = signe(gel(chi,1))? 7: 3;
 3305         gel(V,1) = mkvecsmall2(5, t);
 3306         M[1] = 8;
 3307       }
 3308     }
 3309     else if (e == 5 && c == 3)
 3310     { /* sc: -8 (CHI_2(-1)=-1<=>chi[1]=1) and 8 (CHI_2(-1)=1 <=> chi[1]=0) */
 3311       long t = signe(gel(chi,1))? 7: 3;
 3312       gel(V,1) = mkvecsmalln(6, 2L,4L,5L,6L,8L,t);
 3313       M[1] = 8;
 3314     }
 3315     else if ((e == 4 && c == 2) || (e == 5 && c <= 2) || (e == 6 && c <= 2)
 3316          || (e >= 7 && c == e - 3))
 3317     { /* sc: 4 */
 3318       gel(V,1) = mkvecsmall3(0,2,3);
 3319       M[1] = 4;
 3320     }
 3321     else if ((e <= 4 && c == 0) || (e >= 5 && c == e - 2))
 3322     { /* sc: 2 */
 3323       gel(V,1) = mkvecsmall(0);
 3324       M[1] = 2;
 3325     }
 3326     else if ((e == 6 && c == 3) || (e >= 7 && c <= e - 4))
 3327     { /* sc: -2 */
 3328       gel(V,1) = mkvecsmalln(7, 0L,2L,3L,4L,5L,6L,7L);
 3329       M[1] = 8;
 3330     }
 3331   }
 3332   j = M[1]? 2: 1;
 3333   for (i = odd(N)? 1: 2; i < l; i++) /* skip p=2, done above */
 3334   {
 3335     long p = PN[i], e = EN[i];
 3336     long z = zv_search(P, p), c = z? E[z]: 0; /* c = v_p(FC) */
 3337     if ((e <= 2 && c == 1 && itos(gel(chi,z)) == (p>>1)) /* ord(CHI_p)=2 */
 3338         || (e >= 3 && c <= e - 2))
 3339     { /* sc: -p */
 3340       GEN v = non_residues(p);
 3341       if (e != 1) v = vecsmall_prepend(v, 0);
 3342       gel(V,j) = v;
 3343       M[j] = p; j++;
 3344     }
 3345     else if (e >= 2 && c < e)
 3346     { /* sc: p */
 3347       gel(V,j) = mkvecsmall(0);
 3348       M[j] = p; j++;
 3349     }
 3350   }
 3351   if (j == 1) return cgetg(1, t_VECSMALL);
 3352   setlg(V,j); setlg(M,j); mod = zv_prod(M);
 3353   L = zero_zv(mod);
 3354   for (i = 1; i < j; i++)
 3355   {
 3356     GEN v = gel(V,i);
 3357     long s, m = M[i], lv = lg(v);
 3358     for (s = 1; s < lv; s++)
 3359     {
 3360       long a = v[s] + 1;
 3361       do { L[a] = 1; a += m; } while (a <= mod);
 3362     }
 3363   }
 3364   return L;
 3365 }
 3366 /* v=mfnewzerodata(N,CHI); returns TRUE if newtrace(n) must be zero,
 3367  * (but newtrace(n) may still be zero if we return FALSE) */
 3368 static long
 3369 mfnewchkzero(GEN v, long n) { long q = lg(v)-1; return q && v[(n%q) + 1]; }
 3370 
 3371 /* if (!VCHIP): from mftraceform_cusp;
 3372  * else from initnewtrace and CHI is known to be primitive */
 3373 static GEN
 3374 inittrace(long N, GEN CHI, GEN VCHIP)
 3375 {
 3376   long FC;
 3377   if (VCHIP)
 3378     FC = mfcharmodulus(CHI);
 3379   else
 3380     VCHIP = mfcharinit(mfchartoprimitive(CHI, &FC));
 3381   return mkvecn(5, mksqr(N), mkmup(N), mkgcd(N), VCHIP, mkbez(N, FC));
 3382 }
 3383 
 3384 /* p > 2 prime; return a sorted t_VECSMALL of primes s.t Tr^new(p) = 0 for all
 3385  * weights > 2 */
 3386 static GEN
 3387 inittrconj(long N, long FC)
 3388 {
 3389   GEN fa, P, E, v;
 3390   long i, k, l;
 3391 
 3392   if (FC != 1) return cgetg(1,t_VECSMALL);
 3393 
 3394   fa = myfactoru(N >> vals(N));
 3395   P = gel(fa,1); l = lg(P);
 3396   E = gel(fa,2);
 3397   v = cgetg(l, t_VECSMALL);
 3398   for (i = k = 1; i < l; i++)
 3399   {
 3400     long j, p = P[i]; /* > 2 */
 3401     for (j = 1; j < l; j++)
 3402       if (j != i && E[j] == 1 && kross(-p, P[j]) == 1) v[k++] = p;
 3403   }
 3404   setlg(v,k); return v;
 3405 }
 3406 
 3407 /* assume CHIP primitive, f(CHIP) | N; NZ = mfnewzerodata(N,CHIP) */
 3408 static GEN
 3409 initnewtrace_i(long N, GEN CHIP, GEN NZ)
 3410 {
 3411   GEN T = const_vec(N, cgetg(1,t_VEC)), D, VCHIP;
 3412   long FC = mfcharmodulus(CHIP), N1, N2, i, l;
 3413 
 3414   if (!NZ) NZ = mkvecsmall(1); /*Tr^new = 0; initialize data nevertheless*/
 3415   VCHIP = mfcharinit(CHIP);
 3416   N1 = N/FC; newd_params(N1, &N2);
 3417   D = mydivisorsu(N1/N2); l = lg(D);
 3418   N2 *= FC;
 3419   for (i = 1; i < l; i++)
 3420   {
 3421     long M = D[i]*N2;
 3422     gel(T,M) = inittrace(M, CHIP, VCHIP);
 3423   }
 3424   gel(T,N) = shallowconcat(gel(T,N), mkvec2(NZ, inittrconj(N,FC)));
 3425   return T;
 3426 }
 3427 /* don't initialize if Tr^new = 0, return NULL */
 3428 static GEN
 3429 initnewtrace(long N, GEN CHI)
 3430 {
 3431   GEN CHIP = mfchartoprimitive(CHI, NULL), NZ = mfnewzerodata(N,CHIP);
 3432   return NZ? initnewtrace_i(N, CHIP, NZ): NULL;
 3433 }
 3434 
 3435 /* (-1)^k */
 3436 static long
 3437 m1pk(long k) { return odd(k)? -1 : 1; }
 3438 static long
 3439 badchar(long N, long k, GEN CHI)
 3440 { return mfcharparity(CHI) != m1pk(k) || (CHI && N % mfcharconductor(CHI)); }
 3441 
 3442 /* dimension of space of cusp forms S_k(\G_0(N),CHI)
 3443  * Only depends on CHIP the primitive char attached to CHI */
 3444 long
 3445 mfcuspdim(long N, long k, GEN CHI)
 3446 {
 3447   pari_sp av = avma;
 3448   long FC;
 3449   GEN s;
 3450   if (k <= 0) return 0;
 3451   if (k == 1) return mfwt1cuspdim(N, CHI);
 3452   FC = CHI? mfcharconductor(CHI): 1;
 3453   if (FC == 1) CHI = NULL;
 3454   s = gsub(A1(N, k), gadd(A21(N, k, CHI), A22(N, k, CHI)));
 3455   s = gadd(s, gsubsg(A4(k, FC), A3(N, FC)));
 3456   return gc_long(av, itos(s));
 3457 }
 3458 
 3459 /* dimension of whole space M_k(\G_0(N),CHI)
 3460  * Only depends on CHIP the primitive char attached to CHI; assumes !badchar */
 3461 long
 3462 mffulldim(long N, long k, GEN CHI)
 3463 {
 3464   pari_sp av = avma;
 3465   long FC = CHI? mfcharconductor(CHI): 1;
 3466   GEN s;
 3467   if (k <= 0) return (k == 0 && FC == 1)? 1: 0;
 3468   if (k == 1) return gc_long(av, itos(A3(N, FC)) + mfwt1cuspdim(N, CHI));
 3469   if (FC == 1) CHI = NULL;
 3470   s = gsub(A1(N, k), gadd(A21(N, k, CHI), A22(N, k, CHI)));
 3471   s = gadd(s, A3(N, FC));
 3472   return gc_long(av, itos(s));
 3473 }
 3474 
 3475 /* Dimension of the space of Eisenstein series */
 3476 long
 3477 mfeisensteindim(long N, long k, GEN CHI)
 3478 {
 3479   pari_sp av = avma;
 3480   long s, FC = CHI? mfcharconductor(CHI): 1;
 3481   if (k <= 0) return (k == 0 && FC == 1)? 1: 0;
 3482   s = itos(gmul2n(A3(N, FC), 1));
 3483   if (k > 1) s -= A4(k, FC); else s >>= 1;
 3484   return gc_long(av,s);
 3485 }
 3486 
 3487 enum { _SQRTS = 1, _MUP, _GCD, _VCHIP, _BEZ, _NEWLZ, _TRCONJ };
 3488 /* Trace of T(n) on space of cuspforms; only depends on CHIP the primitive char
 3489  * attached to CHI */
 3490 static GEN
 3491 mfcusptrace_i(long N, long k, long n, GEN Dn, GEN S)
 3492 {
 3493   pari_sp av = avma;
 3494   GEN a, b, VCHIP, GCD;
 3495   long t;
 3496   if (!n) return gen_0;
 3497   VCHIP = gel(S,_VCHIP);
 3498   GCD = gel(S,_GCD);
 3499   t = TA4(k, VCHIP, Dn, GCD);
 3500   a = TA1(N, k, VCHIP, GCD, n); if (t) a = gaddgs(a,t);
 3501   b = TA2(N, k, VCHIP, n, gel(S,_SQRTS), gel(S,_MUP), GCD);
 3502   b = gadd(b, TA3(N, k, VCHIP, GCD, Dn, gel(S,_BEZ)));
 3503   b = gsub(a,b);
 3504   if (typ(b) != t_POL) return gerepileupto(av, b);
 3505   return gerepilecopy(av, vchip_polmod(VCHIP, b));
 3506 }
 3507 
 3508 static GEN
 3509 mfcusptracecache(long N, long k, long n, GEN Dn, GEN S, cachenew_t *cache)
 3510 {
 3511   GEN C = NULL, T = gel(cache->vfull,N);
 3512   long lcache = lg(T);
 3513   if (n < lcache) C = gel(T, n);
 3514   if (C) cache->cuspHIT++; else C = mfcusptrace_i(N, k, n, Dn, S);
 3515   cache->cuspTOTAL++;
 3516   if (n < lcache) gel(T,n) = C;
 3517   return C;
 3518 }
 3519 
 3520 /* return the divisors of n, known to be among the elements of D */
 3521 static GEN
 3522 div_restrict(GEN D, ulong n)
 3523 {
 3524   long i, j, l;
 3525   GEN v, VDIV = caches[cache_DIV].cache;
 3526   if (lg(VDIV) > n) return gel(VDIV,n);
 3527   l = lg(D);
 3528   v = cgetg(l, t_VECSMALL);
 3529   for (i = j = 1; i < l; i++)
 3530   {
 3531     ulong d = D[i];
 3532     if (n % d == 0) v[j++] = d;
 3533   }
 3534   setlg(v,j); return v;
 3535 }
 3536 
 3537 /* for some prime divisors of N, Tr^new(p) = 0 */
 3538 static int
 3539 trconj(GEN T, long N, long n)
 3540 { return (lg(T) > 1 && N % n == 0 && zv_search(T, n)); }
 3541 
 3542 /* n > 0; trace formula on new space */
 3543 static GEN
 3544 mfnewtrace_i(long N, long k, long n, cachenew_t *cache)
 3545 {
 3546   GEN VCHIP, s, Dn, DN1, SN, S = cache->DATA;
 3547   long FC, N1, N2, N1N2, g, i, j, lDN1;
 3548 
 3549   if (!S) return gen_0;
 3550   SN = gel(S,N);
 3551   if (mfnewchkzero(gel(SN,_NEWLZ), n)) return gen_0;
 3552   if (k > 2 && trconj(gel(SN,_TRCONJ), N, n)) return gen_0;
 3553   VCHIP = gel(SN, _VCHIP); FC = vchip_FC(VCHIP);
 3554   N1 = N/FC; newt_params(N1, n, FC, &g, &N2);
 3555   N1N2 = N1/N2;
 3556   DN1 = mydivisorsu(N1N2); lDN1 = lg(DN1);
 3557   N2 *= FC;
 3558   Dn = mydivisorsu(n); /* this one is probably out of cache */
 3559   s = gmulsg(mubeta2(N1N2,n), mfcusptracecache(N2, k, n, Dn, gel(S,N2), cache));
 3560   for (i = 2; i < lDN1; i++)
 3561   { /* skip M1 = 1, done above */
 3562     long M1 = DN1[i], N1M1 = DN1[lDN1-i];
 3563     GEN Dg = mydivisorsu(ugcd(M1, g));
 3564     M1 *= N2;
 3565     s = gadd(s, gmulsg(mubeta2(N1M1,n),
 3566                        mfcusptracecache(M1, k, n, Dn, gel(S,M1), cache)));
 3567     for (j = 2; j < lg(Dg); j++) /* skip d = 1, done above */
 3568     {
 3569       long d = Dg[j], ndd = n/(d*d), M = M1/d;
 3570       GEN z = mulsi(mubeta2(N1M1,ndd), powuu(d,k-1)), C = vchip_lift(VCHIP,d,z);
 3571       GEN Dndd = div_restrict(Dn, ndd);
 3572       s = gadd(s, gmul(C, mfcusptracecache(M, k, ndd, Dndd, gel(S,M), cache)));
 3573     }
 3574     s = vchip_mod(VCHIP, s);
 3575   }
 3576   return vchip_polmod(VCHIP, s);
 3577 }
 3578 
 3579 /* mfcuspdim(N,k,CHI) - mfnewdim(N,k,CHI); CHIP primitive (for efficiency) */
 3580 static long
 3581 mfolddim_i(long N, long k, GEN CHIP)
 3582 {
 3583   long S, i, l, FC = mfcharmodulus(CHIP), N1 = N/FC, N2;
 3584   GEN D;
 3585   newd_params(N1, &N2); /* will ensure mubeta != 0 */
 3586   D = mydivisorsu(N1/N2); l = lg(D);
 3587   N2 *= FC; S = 0;
 3588   for (i = 2; i < l; i++)
 3589   {
 3590     long M = D[l-i]*N2, d = mfcuspdim(M, k, CHIP);
 3591     if (d) S -= mubeta(D[i]) * d;
 3592   }
 3593   return S;
 3594 }
 3595 long
 3596 mfolddim(long N, long k, GEN CHI)
 3597 {
 3598   pari_sp av = avma;
 3599   GEN CHIP = mfchartoprimitive(CHI, NULL);
 3600   return gc_long(av, mfolddim_i(N, k, CHIP));
 3601 }
 3602 /* Only depends on CHIP the primitive char attached to CHI; assumes !badchar */
 3603 long
 3604 mfnewdim(long N, long k, GEN CHI)
 3605 {
 3606   pari_sp av;
 3607   long S;
 3608   GEN CHIP = mfchartoprimitive(CHI, NULL);
 3609   S = mfcuspdim(N, k, CHIP); if (!S) return 0;
 3610   av = avma; return gc_long(av, S - mfolddim_i(N, k, CHIP));
 3611 }
 3612 
 3613 /* trace form, given as closure */
 3614 static GEN
 3615 mftraceform_new(long N, long k, GEN CHI)
 3616 {
 3617   GEN T;
 3618   if (k == 1) return initwt1newtrace(mfinit_Nkchi(N, 1, CHI, mf_CUSP, 0));
 3619   T = initnewtrace(N,CHI); if (!T) return mftrivial();
 3620   return tag(t_MF_NEWTRACE, mkNK(N,k,CHI), T);
 3621 }
 3622 static GEN
 3623 mftraceform_cusp(long N, long k, GEN CHI)
 3624 {
 3625   if (k == 1) return initwt1trace(mfinit_Nkchi(N, 1, CHI, mf_CUSP, 0));
 3626   return tag(t_MF_TRACE, mkNK(N,k,CHI), inittrace(N,CHI,NULL));
 3627 }
 3628 static GEN
 3629 mftraceform_i(GEN NK, long space)
 3630 {
 3631   GEN CHI;
 3632   long N, k;
 3633   checkNK(NK, &N, &k, &CHI, 0);
 3634   if (!mfdim_Nkchi(N, k, CHI, space)) return mftrivial();
 3635   switch(space)
 3636   {
 3637     case mf_NEW: return mftraceform_new(N, k, CHI);
 3638     case mf_CUSP:return mftraceform_cusp(N, k, CHI);
 3639   }
 3640   pari_err_DOMAIN("mftraceform", "space", "=", utoi(space), NK);
 3641   return NULL;/*LCOV_EXCL_LINE*/
 3642 }
 3643 GEN
 3644 mftraceform(GEN NK, long space)
 3645 { pari_sp av = avma; return gerepilecopy(av, mftraceform_i(NK,space)); }
 3646 
 3647 static GEN
 3648 hecke_data(long N, long n)
 3649 { return mkvecsmall3(n, u_ppo(n, N), N); }
 3650 /* 1/2-integral weight */
 3651 static GEN
 3652 heckef2_data(long N, long n)
 3653 {
 3654   ulong f, fN, fN2;
 3655   if (!uissquareall(n, &f)) return NULL;
 3656   fN = u_ppo(f, N); fN2 = fN*fN;
 3657   return mkvec2(myfactoru(fN), mkvecsmall4(n, N, fN2, n/fN2));
 3658 }
 3659 /* N = mf_get_N(F) or a multiple */
 3660 static GEN
 3661 mfhecke_i(long n, long N, GEN F)
 3662 {
 3663   if (n == 1) return F;
 3664   return tag2(t_MF_HECKE, mf_get_NK(F), hecke_data(N,n), F);
 3665 }
 3666 
 3667 GEN
 3668 mfhecke(GEN mf, GEN F, long n)
 3669 {
 3670   pari_sp av = avma;
 3671   GEN NK, CHI, gk, DATA;
 3672   long N, nk, dk;
 3673   mf = checkMF(mf);
 3674   if (!checkmf_i(F)) pari_err_TYPE("mfhecke",F);
 3675   if (n <= 0) pari_err_TYPE("mfhecke [n <= 0]", stoi(n));
 3676   if (n == 1) return gcopy(F);
 3677   gk = mf_get_gk(F);
 3678   Qtoss(gk,&nk,&dk);
 3679   CHI = mf_get_CHI(F);
 3680   N = MF_get_N(mf);
 3681   if (dk == 2)
 3682   {
 3683     DATA = heckef2_data(N,n);
 3684     if (!DATA) return mftrivial();
 3685   }
 3686   else
 3687     DATA = hecke_data(N,n);
 3688   NK = mkgNK(lcmii(stoi(N), mf_get_gN(F)), gk, CHI, mf_get_field(F));
 3689   return gerepilecopy(av, tag2(t_MF_HECKE, NK, DATA, F));
 3690 }
 3691 
 3692 /* form F given by closure, compute B(d)(F) as closure (q -> q^d) */
 3693 static GEN
 3694 mfbd_i(GEN F, long d)
 3695 {
 3696   GEN D, NK, gk, CHI;
 3697   if (d == 1) return F;
 3698   if (d <= 0) pari_err_TYPE("mfbd [d <= 0]", stoi(d));
 3699   if (mf_get_type(F) != t_MF_BD) D = utoi(d);
 3700   else { D = mului(d, gel(F,3)); F = gel(F,2); }
 3701   gk = mf_get_gk(F); CHI = mf_get_CHI(F);
 3702   if (typ(gk) != t_INT) CHI = mfcharmul(CHI, get_mfchar(utoi(d << 2)));
 3703   NK = mkgNK(muliu(mf_get_gN(F), d), gk, CHI, mf_get_field(F));
 3704   return tag2(t_MF_BD, NK, F, D);
 3705 }
 3706 GEN
 3707 mfbd(GEN F, long d)
 3708 {
 3709   pari_sp av = avma;
 3710   if (!checkmf_i(F)) pari_err_TYPE("mfbd",F);
 3711   return gerepilecopy(av, mfbd_i(F, d));
 3712 }
 3713 
 3714 /* A[i+1] = a(t*i^2) */
 3715 static GEN
 3716 RgV_shimura(GEN A, long n, long t, long N, long r, GEN CHI)
 3717 {
 3718   GEN R, a0, Pn = mfcharpol(CHI);
 3719   long m, st, ord = mfcharorder(CHI), vt = varn(Pn), Nt = t == 1? N: ulcm(N,t);
 3720 
 3721   R = cgetg(n + 2, t_VEC);
 3722   st = odd(r)? -t: t;
 3723   a0 = gel(A, 1);
 3724   if (!gequal0(a0))
 3725   {
 3726     long o = mfcharorder(CHI);
 3727     if (st != 1 && odd(o)) o <<= 1;
 3728     a0 = gmul(a0, charLFwtk(Nt, r, CHI, o, st));
 3729   }
 3730   gel(R, 1) = a0;
 3731   for (m = 1; m <= n; m++)
 3732   {
 3733     GEN Dm = mydivisorsu(u_ppo(m, Nt)), S = gel(A, m*m + 1);
 3734     long i, l = lg(Dm);
 3735     for (i = 2; i < l; i++)
 3736     { /* (e,Nt) = 1; skip i = 1: e = 1, done above */
 3737       long e = Dm[i], me = m / e, a = mfcharevalord(CHI, e, ord);
 3738       GEN c, C = powuu(e, r - 1);
 3739       if (kross(st, e) == -1) C = negi(C);
 3740       c = Qab_Czeta(a, ord, C, vt);
 3741       S = gadd(S, gmul(c, gel(A, me*me + 1)));
 3742     }
 3743     gel(R, m+1) = S;
 3744   }
 3745   return degpol(Pn) > 1? gmodulo(R, Pn): R;
 3746 }
 3747 
 3748 static long
 3749 mfisinkohnen(GEN mf, GEN F)
 3750 {
 3751   GEN v, gk = MF_get_gk(mf), CHI = MF_get_CHI(mf);
 3752   long i, eps, N4 = MF_get_N(mf) >> 2, sb = mfsturmNgk(N4 << 4, gk) + 1;
 3753   eps = N4 % mfcharconductor(CHI)? -1 : 1;
 3754   if (odd(MF_get_r(mf))) eps = -eps;
 3755   v = mfcoefs(F, sb, 1);
 3756   for (i = 2;     i <= sb; i+=4) if (!gequal0(gel(v,i+1))) return 0;
 3757   for (i = 2+eps; i <= sb; i+=4) if (!gequal0(gel(v,i+1))) return 0;
 3758   return 1;
 3759 }
 3760 
 3761 static long
 3762 mfshimura_space_cusp(GEN mf)
 3763 {
 3764   long N4;
 3765   if (MF_get_r(mf) == 1 && (N4 = MF_get_N(mf) >> 2) >= 4)
 3766   {
 3767     GEN E = gel(myfactoru(N4), 2);
 3768     long ma = vecsmall_max(E);
 3769     if (ma > 2 || (ma == 2 && !mfcharistrivial(MF_get_CHI(mf)))) return 0;
 3770   }
 3771   return 1;
 3772 }
 3773 
 3774 /* D is either a discriminant (not necessarily fundamental) with
 3775    sign(D)=(-1)^{k-1/2}*eps, or a positive squarefree integer t, which is then
 3776    transformed into a fundamental discriminant of the correct sign. */
 3777 GEN
 3778 mfshimura(GEN mf, GEN F, long t)
 3779 {
 3780   pari_sp av = avma;
 3781   GEN G, res, mf2, CHI;
 3782   long sb, M, r, N, space = mf_FULL;
 3783 
 3784   if (!checkmf_i(F)) pari_err_TYPE("mfshimura",F);
 3785   mf = checkMF(mf);
 3786   r = MF_get_r(mf);
 3787   if (r <= 0) pari_err_DOMAIN("mfshimura", "weight", "<=", ghalf, mf_get_gk(F));
 3788   if (t <= 0 || !uissquarefree(t)) pari_err_TYPE("mfshimura [t]", stoi(t));
 3789   N = MF_get_N(mf); M = N >> 1;
 3790   if (mfiscuspidal(mf,F))
 3791   {
 3792     if (mfshimura_space_cusp(mf)) space = mf_CUSP;
 3793     if (mfisinkohnen(mf,F)) M = N >> 2;
 3794   }
 3795   CHI = MF_get_CHI(mf);
 3796   mf2 = mfinit_Nkchi(M, r << 1, mfcharpow(CHI, gen_2), space, 0);
 3797   sb = mfsturm(mf2);
 3798   G = RgV_shimura(mfcoefs_i(F, sb*sb, t), sb, t, N, r, CHI);
 3799   res = mftobasis_i(mf2, G);
 3800   /* not mflinear(mf2,): we want lowest possible level */
 3801   G = mflinear(MF_get_basis(mf2), res);
 3802   return gerepilecopy(av, mkvec3(mf2, G, res));
 3803 }
 3804 
 3805 /* W ZabM (ZM if n = 1), a t_INT or NULL, b t_INT, ZXQ mod P or NULL.
 3806  * Write a/b = A/d with d t_INT and A Zab return [W,d,A,P] */
 3807 static GEN
 3808 mkMinv(GEN W, GEN a, GEN b, GEN P)
 3809 {
 3810   GEN A = (b && typ(b) == t_POL)? Q_remove_denom(QXQ_inv(b,P), &b): NULL;
 3811   if (a && b)
 3812   {
 3813     a = Qdivii(a,b);
 3814     if (typ(a) == t_INT) b = gen_1; else { b = gel(a,2); a = gel(a,1); }
 3815     if (is_pm1(a)) a = NULL;
 3816   }
 3817   if (a) A = A? ZX_Z_mul(A,a): a; else if (!A) A = gen_1;
 3818   if (!b) b = gen_1;
 3819   if (!P) P = gen_0;
 3820   return mkvec4(W,b,A,P);
 3821 }
 3822 /* M square invertible QabM, return [M',d], M*M' = d*Id */
 3823 static GEN
 3824 QabM_Minv(GEN M, GEN P, long n)
 3825 {
 3826   GEN dW, W, dM;
 3827   M = Q_remove_denom(M, &dM);
 3828   W = P? ZabM_inv(liftpol_shallow(M), P, n, &dW): ZM_inv(M, &dW);
 3829   return mkMinv(W, dM, dW, P);
 3830 }
 3831 /* Simplified form of mfclean, after a QabM_indexrank: M a ZabM with full
 3832  * column rank and z = indexrank(M) is known */
 3833 static GEN
 3834 mfclean2(GEN M, GEN z, GEN P, long n)
 3835 {
 3836   GEN d, Minv, y = gel(z,1), W = rowpermute(M, y);
 3837   W = P? ZabM_inv(liftpol_shallow(W), P, n, &d): ZM_inv(W, &d);
 3838   M = rowslice(M, 1, y[lg(y)-1]);
 3839   Minv = mkMinv(W, NULL, d, P);
 3840   return mkvec3(y, Minv, M);
 3841 }
 3842 /* M QabM, lg(M)>1 and [y,z] its rank profile. Let Minv be the inverse of the
 3843  * invertible square matrix in mkMinv format. Return [y,Minv, M[..y[#y],]]
 3844  * P cyclotomic polynomial of order n > 2 or NULL */
 3845 static GEN
 3846 mfclean(GEN M, GEN P, long n, int ratlift)
 3847 {
 3848   GEN W, v, y, z, d, Minv, dM, MdM = Q_remove_denom(M, &dM);
 3849   if (n <= 2)
 3850     W = ZM_pseudoinv(MdM, &v, &d);
 3851   else
 3852     W = ZabM_pseudoinv_i(liftpol_shallow(MdM), P, n, &v, &d, ratlift);
 3853   y = gel(v,1);
 3854   z = gel(v,2);
 3855   if (lg(z) != lg(MdM)) M = vecpermute(M,z);
 3856   M = rowslice(M, 1, y[lg(y)-1]);
 3857   Minv = mkMinv(W, dM, d, P);
 3858   return mkvec3(y, Minv, M);
 3859 }
 3860 /* call mfclean using only CHI */
 3861 static GEN
 3862 mfcleanCHI(GEN M, GEN CHI, int ratlift)
 3863 {
 3864   long n = mfcharorder(CHI);
 3865   GEN P = (n <= 2)? NULL: mfcharpol(CHI);
 3866   return mfclean(M, P, n, ratlift);
 3867 }
 3868 
 3869 /* DATA component of a t_MF_NEWTRACE. Was it stripped to save memory ? */
 3870 static int
 3871 newtrace_stripped(GEN DATA)
 3872 { return DATA && (lg(DATA) == 5 && typ(gel(DATA,3)) == t_INT); }
 3873 /* f a t_MF_NEWTRACE */
 3874 static GEN
 3875 newtrace_DATA(long N, GEN f)
 3876 {
 3877   GEN DATA = gel(f,2);
 3878   return newtrace_stripped(DATA)? initnewtrace(N, DATA): DATA;
 3879 }
 3880 /* reset cachenew for new level incorporating new DATA, tf a t_MF_NEWTRACE
 3881  * (+ possibly initialize 'full' for new allowed levels) */
 3882 static void
 3883 reset_cachenew(cachenew_t *cache, long N, GEN tf)
 3884 {
 3885   long i, n, l;
 3886   GEN v, DATA = newtrace_DATA(N,tf);
 3887   cache->DATA = DATA;
 3888   if (!DATA) return;
 3889   n = cache->n;
 3890   v = cache->vfull; l = N+1; /* = lg(DATA) */
 3891   for (i = 1; i < l; i++)
 3892     if (typ(gel(v,i)) == t_INT && lg(gel(DATA,i)) != 1)
 3893       gel(v,i) = const_vec(n, NULL);
 3894   cache->VCHIP = gel(gel(DATA,N),_VCHIP);
 3895 }
 3896 /* initialize a cache of newtrace / cusptrace up to index n and level | N;
 3897  * DATA may be NULL (<=> Tr^new = 0). tf a t_MF_NEWTRACE */
 3898 static void
 3899 init_cachenew(cachenew_t *cache, long n, long N, GEN tf)
 3900 {
 3901   long i, l = N+1; /* = lg(tf.DATA) when DATA != NULL */
 3902   GEN v;
 3903   cache->n = n;
 3904   cache->vnew = v = cgetg(l, t_VEC);
 3905   for (i = 1; i < l; i++) gel(v,i) = (N % i)? gen_0: const_vec(n, NULL);
 3906   cache->newHIT = cache->newTOTAL = cache->cuspHIT = cache->cuspTOTAL = 0;
 3907   cache->vfull = v = zerovec(N);
 3908   reset_cachenew(cache, N, tf);
 3909 }
 3910 static void
 3911 dbg_cachenew(cachenew_t *C)
 3912 {
 3913   if (DEBUGLEVEL >= 2 && C)
 3914     err_printf("newtrace cache hits: new = %ld/%ld, cusp = %ld/%ld\n",
 3915                     C->newHIT, C->newTOTAL, C->cuspHIT, C->cuspTOTAL);
 3916 }
 3917 
 3918 /* newtrace_{N,k}(d*i), i = n0, ..., n */
 3919 static GEN
 3920 colnewtrace(long n0, long n, long d, long N, long k, cachenew_t *cache)
 3921 {
 3922   GEN v = cgetg(n-n0+2, t_COL);
 3923   long i;
 3924   for (i = n0; i <= n; i++) gel(v, i-n0+1) = mfnewtracecache(N, k, i*d, cache);
 3925   return v;
 3926 }
 3927 /* T_n(l*m0, l*(m0+1), ..., l*m) F, F = t_MF_NEWTRACE [N,k],DATA, cache
 3928  * contains DATA != NULL as well as cached values of F */
 3929 static GEN
 3930 heckenewtrace(long m0, long m, long l, long N, long NBIG, long k, long n, cachenew_t *cache)
 3931 {
 3932   long lD, a, k1, nl = n*l;
 3933   GEN D, V, v = colnewtrace(m0, m, nl, N, k, cache); /* d=1 */
 3934   GEN VCHIP;
 3935   if (n == 1) return v;
 3936   VCHIP = cache->VCHIP;
 3937   D = mydivisorsu(u_ppo(n, NBIG)); lD = lg(D);
 3938   k1 = k - 1;
 3939   for (a = 2; a < lD; a++)
 3940   { /* d > 1, (d,NBIG) = 1 */
 3941     long i, j, d = D[a], c = ugcd(l, d), dl = d/c, m0d = ceildiv(m0, dl);
 3942     GEN C = vchip_lift(VCHIP, d, powuu(d, k1));
 3943     /* m0=0: i = 1 => skip F(0) = 0 */
 3944     if (!m0) { i = 1; j = dl; } else { i = 0; j = m0d*dl; }
 3945     V = colnewtrace(m0d, m/dl, nl/(d*c), N, k, cache);
 3946     /* C = chi(d) d^(k-1) */
 3947     for (; j <= m; i++, j += dl)
 3948       gel(v,j-m0+1) = gadd(gel(v,j-m0+1), vchip_mod(VCHIP, gmul(C,gel(V,i+1))));
 3949   }
 3950   return v;
 3951 }
 3952 
 3953 /* Given v = an[i], return an[d*i] */
 3954 static GEN
 3955 anextract(GEN v, long n, long d)
 3956 {
 3957   GEN w = cgetg(n+2, t_VEC);
 3958   long i;
 3959   for (i = 0; i <= n; i++) gel(w, i+1) = gel(v, i*d+1);
 3960   return w;
 3961 }
 3962 /* T_n(F)(0, l, ..., l*m) */
 3963 static GEN
 3964 hecke_i(long m, long l, GEN V, GEN F, GEN DATA)
 3965 {
 3966   long k, n, nNBIG, NBIG, lD, M, a, t, nl;
 3967   GEN D, v, CHI;
 3968   if (typ(DATA) == t_VEC)
 3969   { /* 1/2-integral k */
 3970     if (!V) { GEN S = gel(DATA,2); V = mfcoefs_i(F, m*l*S[3], S[4]); }
 3971     return RgV_heckef2(m, l, V, F, DATA);
 3972   }
 3973   k = mf_get_k(F);
 3974   n = DATA[1]; nl = n*l;
 3975   nNBIG = DATA[2];
 3976   NBIG = DATA[3];
 3977   if (nNBIG == 1) return V? V: mfcoefs_i(F,m,nl);
 3978   if (!V && mf_get_type(F) == t_MF_NEWTRACE)
 3979   { /* inline F to allow cache, T_n at level NBIG acting on Tr^new(N,k,CHI) */
 3980     cachenew_t cache;
 3981     long N = mf_get_N(F);
 3982     init_cachenew(&cache, m*nl, N, F);
 3983     v = heckenewtrace(0, m, l, N, NBIG, k, n, &cache);
 3984     dbg_cachenew(&cache);
 3985     settyp(v, t_VEC); return v;
 3986   }
 3987   CHI = mf_get_CHI(F);
 3988   D = mydivisorsu(nNBIG); lD = lg(D);
 3989   M = m + 1;
 3990   t = nNBIG * ugcd(nNBIG, l);
 3991   if (!V) V = mfcoefs_i(F, m * t, nl / t); /* usually nl = t */
 3992   v = anextract(V, m, t); /* mfcoefs(F, m, nl); d = 1 */
 3993   for (a = 2; a < lD; a++)
 3994   { /* d > 1, (d, NBIG) = 1 */
 3995     long d = D[a], c = ugcd(l, d), dl = d/c, i, idl;
 3996     GEN C = gmul(mfchareval(CHI, d), powuu(d, k-1));
 3997     GEN w = anextract(V, m/dl, t/(d*c)); /* mfcoefs(F, m/dl, nl/(d*c)) */
 3998     for (i = idl = 1; idl <= M; i++, idl += dl)
 3999       gel(v,idl) = gadd(gel(v,idl), gmul(C, gel(w,i)));
 4000   }
 4001   return v;
 4002 }
 4003 
 4004 static GEN
 4005 mkmf(GEN x1, GEN x2, GEN x3, GEN x4, GEN x5)
 4006 {
 4007   GEN MF = obj_init(5, MF_SPLITN);
 4008   gel(MF,1) = x1;
 4009   gel(MF,2) = x2;
 4010   gel(MF,3) = x3;
 4011   gel(MF,4) = x4;
 4012   gel(MF,5) = x5; return MF;
 4013 }
 4014 
 4015 /* return an integer b such that p | b => T_p^k Tr^new = 0, for all k > 0 */
 4016 static long
 4017 get_badj(long N, long FC)
 4018 {
 4019   GEN fa = myfactoru(N), P = gel(fa,1), E = gel(fa,2);
 4020   long i, b = 1, l = lg(P);
 4021   for (i = 1; i < l; i++)
 4022     if (E[i] > 1 && u_lval(FC, P[i]) < E[i]) b *= P[i];
 4023   return b;
 4024 }
 4025 /* in place, assume perm strictly increasing */
 4026 static void
 4027 vecpermute_inplace(GEN v, GEN perm)
 4028 {
 4029   long i, l = lg(perm);
 4030   for (i = 1; i < l; i++) gel(v,i) = gel(v,perm[i]);
 4031 }
 4032 
 4033 /* Find basis of newspace using closures; assume k >= 2 and !badchar.
 4034  * Return NULL if space is empty, else
 4035  * [mf1, list of closures T(j)traceform, list of corresponding j, matrix] */
 4036 static GEN
 4037 mfnewinit(long N, long k, GEN CHI, cachenew_t *cache, long init)
 4038 {
 4039   GEN S, vj, M, CHIP, mf1, listj, P, tf;
 4040   long j, ct, ctlj, dim, jin, SB, sb, two, ord, FC, badj;
 4041 
 4042   dim = mfnewdim(N, k, CHI);
 4043   if (!dim && !init) return NULL;
 4044   sb = mfsturmNk(N, k);
 4045   CHIP = mfchartoprimitive(CHI, &FC);
 4046   /* remove newtrace data from S to save space in output: negligible slowdown */
 4047   tf = tag(t_MF_NEWTRACE, mkNK(N,k,CHIP), CHIP);
 4048   badj = get_badj(N, FC);
 4049   /* try sbsmall first: Sturm bound not sharp for new space */
 4050   SB = ceilA1(N, k);
 4051   listj = cgetg(2*sb + 3, t_VECSMALL);
 4052   for (j = ctlj = 1; ctlj < 2*sb + 3; j++)
 4053     if (ugcd(j, badj) == 1) listj[ctlj++] = j;
 4054   if (init)
 4055   {
 4056     init_cachenew(cache, (SB+1)*listj[dim+1], N, tf);
 4057     if (init == -1 || !dim) return NULL; /* old space or dim = 0 */
 4058   }
 4059   else
 4060     reset_cachenew(cache, N, tf);
 4061   /* cache.DATA is not NULL */
 4062   ord = mfcharorder(CHIP);
 4063   P = ord <= 2? NULL: mfcharpol(CHIP);
 4064   vj = cgetg(dim+1, t_VECSMALL);
 4065   M = cgetg(dim+1, t_MAT);
 4066   for (two = 1, ct = 0, jin = 1; two <= 2; two++)
 4067   {
 4068     long a, jlim = jin + sb;
 4069     for (a = jin; a <= jlim; a++)
 4070     {
 4071       GEN z, vecz;
 4072       ct++; vj[ct] = listj[a];
 4073       gel(M, ct) = heckenewtrace(0, SB, 1, N, N, k, vj[ct], cache);
 4074       if (ct < dim) continue;
 4075 
 4076       z = QabM_indexrank(M, P, ord);
 4077       vecz = gel(z, 2); ct = lg(vecz) - 1;
 4078       if (ct == dim) { M = mkvec3(z, gen_0, M); break; } /*maximal rank, done*/
 4079       vecpermute_inplace(M, vecz);
 4080       vecpermute_inplace(vj, vecz);
 4081     }
 4082     if (a <= jlim) break;
 4083     /* sbsmall was not sufficient, use Sturm bound: must extend M */
 4084     for (j = 1; j <= ct; j++)
 4085     {
 4086       GEN t = heckenewtrace(SB + 1, sb, 1, N, N, k, vj[j], cache);
 4087       gel(M,j) = shallowconcat(gel(M, j), t);
 4088     }
 4089     jin = jlim + 1; SB = sb;
 4090   }
 4091   S = cgetg(dim + 1, t_VEC);
 4092   for (j = 1; j <= dim; j++) gel(S, j) = mfhecke_i(vj[j], N, tf);
 4093   dbg_cachenew(cache);
 4094   mf1 =