Code coverage tests

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

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

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

LCOV - code coverage report
Current view: top level - basemath - modsym.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.10.0 lcov report (development 21343-6216058) Lines: 2147 2250 95.4 %
Date: 2017-11-19 06:21:17 Functions: 230 232 99.1 %
Legend: Lines: hit not hit

          Line data    Source code
       1             : /* Copyright (C) 2011  The PARI group.
       2             : 
       3             : This file is part of the PARI/GP package.
       4             : 
       5             : PARI/GP is free software; you can redistribute it and/or modify it under the
       6             : terms of the GNU General Public License as published by the Free Software
       7             : Foundation. It is distributed in the hope that it will be useful, but WITHOUT
       8             : ANY WARRANTY WHATSOEVER.
       9             : 
      10             : Check the License for details. You should have received a copy of it, along
      11             : with the package; see the file 'COPYING'. If not, write to the Free Software
      12             : Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */
      13             : 
      14             : #include "pari.h"
      15             : #include "paripriv.h"
      16             : 
      17             : /* Adapted from shp_package/moments by Robert Pollack
      18             :  * http://www.math.mcgill.ca/darmon/programs/shp/shp.html */
      19             : static GEN mskinit(ulong N, long k, long sign);
      20             : static GEN mshecke_i(GEN W, ulong p);
      21             : static GEN ZSl2_star(GEN v);
      22             : static GEN getMorphism(GEN W1, GEN W2, GEN v);
      23             : static GEN voo_act_Gl2Q(GEN g, long k);
      24             : 
      25             : /* Input: P^1(Z/NZ) (formed by create_p1mod)
      26             :    Output: # P^1(Z/NZ) */
      27             : static long
      28        4872 : p1_size(GEN p1N) { return lg(gel(p1N,1)) - 1; }
      29             : static ulong
      30     5085696 : p1N_get_N(GEN p1N) { return gel(p1N,3)[2]; }
      31             : static GEN
      32     2095345 : p1N_get_hash(GEN p1N) { return gel(p1N,2); }
      33             : static GEN
      34        1323 : p1N_get_fa(GEN p1N) { return gel(p1N,4); }
      35             : static GEN
      36        1218 : p1N_get_div(GEN p1N) { return gel(p1N,5); }
      37             : static GEN
      38     2002630 : p1N_get_invsafe(GEN p1N) { return gel(p1N,6); }
      39             : static GEN
      40      543676 : p1N_get_inverse(GEN p1N) { return gel(p1N,7); }
      41             : /* ms-specific accessors */
      42             : /* W = msinit or msfromell */
      43             : static GEN
      44     1058827 : get_ms(GEN W) { return lg(W) == 4? gel(W,1): W; }
      45             : static GEN
      46       11669 : ms_get_p1N(GEN W) { W = get_ms(W); return gel(W,1); }
      47             : static long
      48        4256 : ms_get_N(GEN W) { return p1N_get_N(ms_get_p1N(W)); }
      49             : static GEN
      50        1659 : ms_get_hashcusps(GEN W) { W = get_ms(W); return gel(W,16); }
      51             : static GEN
      52      584766 : ms_get_section(GEN W) { W = get_ms(W); return gel(W,12); }
      53             : static GEN
      54      164808 : ms_get_genindex(GEN W) { W = get_ms(W); return gel(W,5); }
      55             : static long
      56      157920 : ms_get_nbgen(GEN W) { return lg(ms_get_genindex(W))-1; }
      57             : static long
      58      132202 : ms_get_nbE1(GEN W)
      59             : {
      60             :   GEN W11;
      61      132202 :   W = get_ms(W); W11 = gel(W,11);
      62      132202 :   return W11[4] - W11[3];
      63             : }
      64             : /* msk-specific accessors */
      65             : static long
      66           7 : msk_get_dim(GEN W) { return gmael(W,3,2)[2]; }
      67             : static GEN
      68       61341 : msk_get_basis(GEN W) { return gmael(W,3,1); }
      69             : static long
      70       11669 : msk_get_weight(GEN W) { return gmael(W,3,2)[1]; }
      71             : static GEN
      72       58373 : msk_get_st(GEN W) { return gmael(W,3,3); }
      73             : static GEN
      74       58373 : msk_get_link(GEN W) { return gmael(W,3,4); }
      75             : static GEN
      76       58373 : msk_get_invphiblock(GEN W) { return gmael(W,3,5); }
      77             : static long
      78        4592 : msk_get_sign(GEN W)
      79             : {
      80        4592 :   GEN t = gel(W,2);
      81        4592 :   return typ(t)==t_INT? 0: itos(gel(t,1));
      82             : }
      83             : static GEN
      84         553 : msk_get_star(GEN W) { return gmael(W,2,2); }
      85             : static GEN
      86        3472 : msk_get_starproj(GEN W) { return gmael(W,2,3); }
      87             : 
      88             : long
      89           7 : msgetlevel(GEN W) { checkms(W); return ms_get_N(W); }
      90             : long
      91           7 : msgetweight(GEN W) { checkms(W); return msk_get_weight(W); }
      92             : long
      93          21 : msgetsign(GEN W) { checkms(W); return msk_get_sign(W); }
      94             : 
      95             : void
      96        5761 : checkms(GEN W)
      97             : {
      98        5761 :   if (typ(W) != t_VEC || lg(W) != 4)
      99           0 :     pari_err_TYPE("checkms [please apply msinit]", W);
     100        5761 : }
     101             : 
     102             : /** MODULAR TO SYM **/
     103             : 
     104             : /* q a t_FRAC or t_INT */
     105             : static GEN
     106        6552 : Q_log_init(ulong N, GEN q)
     107             : {
     108             :   long l, n;
     109             :   GEN Q;
     110             : 
     111        6552 :   q = gboundcf(q, 0);
     112        6552 :   l = lg(q);
     113        6552 :   Q = cgetg(l, t_VECSMALL);
     114        6552 :   Q[1] = 1;
     115        6552 :   for (n=2; n <l; n++) Q[n] = umodiu(gel(q,n), N);
     116       15547 :   for (n=3; n < l; n++)
     117        8995 :     Q[n] = Fl_add(Fl_mul(Q[n], Q[n-1], N), Q[n-2], N);
     118        6552 :   return Q;
     119             : }
     120             : 
     121             : /** INIT MODSYM STRUCTURE, WEIGHT 2 **/
     122             : 
     123             : /* num = [Gamma : Gamma_0(N)] = N * Prod_{p|N} (1+p^-1) */
     124             : static ulong
     125        1218 : count_Manin_symbols(ulong N, GEN P)
     126             : {
     127        1218 :   long i, l = lg(P);
     128        1218 :   ulong num = N;
     129        1218 :   for (i = 1; i < l; i++) { ulong p = P[i]; num *= p+1; num /= p; }
     130        1218 :   return num;
     131             : }
     132             : /* returns the list of "Manin symbols" (c,d) in (Z/NZ)^2, (c,d,N) = 1
     133             :  * generating H^1(X_0(N), Z) */
     134             : static GEN
     135        1218 : generatemsymbols(ulong N, ulong num, GEN divN)
     136             : {
     137        1218 :   GEN ret = cgetg(num+1, t_VEC);
     138        1218 :   ulong c, d, curn = 0;
     139             :   long i, l;
     140             :   /* generate Manin-symbols in two lists: */
     141             :   /* list 1: (c:1) for 0 <= c < N */
     142        1218 :   for (c = 0; c < N; c++) gel(ret, ++curn) = mkvecsmall2(c, 1);
     143        1218 :   if (N == 1) return ret;
     144             :   /* list 2: (c:d) with 1 <= c < N, c | N, 0 <= d < N, gcd(d,N) > 1, gcd(c,d)=1.
     145             :    * Furthermore, d != d0 (mod N/c) with c,d0 already in the list */
     146        1218 :   l = lg(divN) - 1;
     147             :   /* c = 1 first */
     148        1218 :   gel(ret, ++curn) = mkvecsmall2(1,0);
     149      122304 :   for (d = 2; d < N; d++)
     150      121086 :     if (ugcd(d,N) != 1UL)
     151       42987 :       gel(ret, ++curn) = mkvecsmall2(1,d);
     152             :   /* omit c = 1 (first) and c = N (last) */
     153        3346 :   for (i=2; i < l; i++)
     154             :   {
     155             :     ulong Novc, d0;
     156        2128 :     c = divN[i];
     157        2128 :     Novc = N / c;
     158       62664 :     for (d0 = 2; d0 <= Novc; d0++)
     159             :     {
     160       60536 :       ulong k, d = d0;
     161       60536 :       if (ugcd(d, Novc) == 1UL) continue;
     162       88067 :       for (k = 0; k < c; k++, d += Novc)
     163       79576 :         if (ugcd(c,d) == 1UL)
     164             :         {
     165       11186 :           gel(ret, ++curn) = mkvecsmall2(c,d);
     166       11186 :           break;
     167             :         }
     168             :     }
     169             :   }
     170        1218 :   if (curn != num) pari_err_BUG("generatemsymbols [wrong number of symbols]");
     171        1218 :   return ret;
     172             : }
     173             : 
     174             : static GEN
     175        1218 : inithashmsymbols(ulong N, GEN symbols)
     176             : {
     177        1218 :   GEN H = zerovec(N);
     178        1218 :   long k, l = lg(symbols);
     179             :   /* skip the (c:1), 0 <= c < N and (1:0) */
     180       55391 :   for (k=N+2; k < l; k++)
     181             :   {
     182       54173 :     GEN s = gel(symbols, k);
     183       54173 :     ulong c = s[1], d = s[2], Novc = N/c;
     184       54173 :     if (gel(H,c) == gen_0) gel(H,c) = const_vecsmall(Novc+1,0);
     185       54173 :     if (c != 1) { d %= Novc; if (!d) d = Novc; }
     186       54173 :     mael(H, c, d) = k;
     187             :   }
     188        1218 :   return H;
     189             : }
     190             : 
     191             : /** Helper functions for Sl2(Z) / Gamma_0(N) **/
     192             : /* [a,b;c,d] */
     193             : static GEN
     194     1301510 : mkmat22(GEN a, GEN b, GEN c, GEN d) { retmkmat2(mkcol2(a,c),mkcol2(b,d)); }
     195             : /* M a 2x2 ZM in SL2(Z) */
     196             : static GEN
     197     1046493 : SL2_inv(GEN M)
     198             : {
     199     1046493 :   GEN a=gcoeff(M,1,1), b=gcoeff(M,1,2), c=gcoeff(M,2,1), d=gcoeff(M,2,2);
     200     1046493 :   return mkmat22(d,negi(b), negi(c),a);
     201             : }
     202             : /* M a 2x2 mat2 in SL2(Z) */
     203             : static GEN
     204      573146 : sl2_inv(GEN M)
     205             : {
     206      573146 :   long a=coeff(M,1,1), b=coeff(M,1,2), c=coeff(M,2,1), d=coeff(M,2,2);
     207      573146 :   return mkvec2(mkvecsmall2(d, -c), mkvecsmall2(-b, a));
     208             : }
     209             : /* Return the mat2 [a,b; c,d], not a zm to avoid GP problems */
     210             : static GEN
     211     1027845 : mat2(long a, long b, long c, long d)
     212     1027845 : { return mkvec2(mkvecsmall2(a,c), mkvecsmall2(b,d)); }
     213             : static GEN
     214      278306 : mat2_to_ZM(GEN M)
     215             : {
     216      278306 :   GEN A = gel(M,1), B = gel(M,2);
     217      278306 :   retmkmat2(mkcol2s(A[1],A[2]), mkcol2s(B[1],B[2]));
     218             : }
     219             : 
     220             : /* Input: a = 2-vector = path = {r/s,x/y}
     221             :  * Output: either [r,x;s,y] or [-r,x;-s,y], whichever has determinant > 0 */
     222             : static GEN
     223       62559 : path_to_ZM(GEN a)
     224             : {
     225       62559 :   GEN v = gel(a,1), w = gel(a,2);
     226       62559 :   long r = v[1], s = v[2], x = w[1], y = w[2];
     227       62559 :   if (cmpii(mulss(r,y), mulss(x,s)) < 0) { r = -r; s = -s; }
     228       62559 :   return mkmat22(stoi(r),stoi(x),stoi(s),stoi(y));
     229             : }
     230             : static GEN
     231      630091 : path_to_zm(GEN a)
     232             : {
     233      630091 :   GEN v = gel(a,1), w = gel(a,2);
     234      630091 :   long r = v[1], s = v[2], x = w[1], y = w[2];
     235      630091 :   if (cmpii(mulss(r,y), mulss(x,s)) < 0) { r = -r; s = -s; }
     236      630091 :   return mat2(r,x,s,y);
     237             : }
     238             : /* path from c1 to c2 */
     239             : static GEN
     240      362285 : mkpath(GEN c1, GEN c2) { return mat2(c1[1], c2[1], c1[2], c2[2]); }
     241             : static long
     242      512309 : cc(GEN M) { GEN v = gel(M,1); return v[2]; }
     243             : static long
     244      512309 : dd(GEN M) { GEN v = gel(M,2); return v[2]; }
     245             : 
     246             : /*Input: a,b = 2 paths, N = integer
     247             :  *Output: 1 if the a,b are \Gamma_0(N)-equivalent; 0 otherwise */
     248             : static int
     249       59087 : gamma_equiv(GEN a, GEN b, ulong N)
     250             : {
     251       59087 :   pari_sp av = avma;
     252       59087 :   GEN m = path_to_zm(a);
     253       59087 :   GEN n = path_to_zm(b);
     254       59087 :   GEN d = subii(mulss(cc(m),dd(n)), mulss(dd(m),cc(n)));
     255       59087 :   ulong res = umodiu(d, N);
     256       59087 :   avma = av; return res == 0;
     257             : }
     258             : /* Input: a,b = 2 paths that are \Gamma_0(N)-equivalent, N = integer
     259             :  * Output: M in \Gamma_0(N) such that Mb=a */
     260             : static GEN
     261       31073 : gamma_equiv_matrix(GEN a, GEN b)
     262             : {
     263       31073 :   GEN m = path_to_ZM(a);
     264       31073 :   GEN n = path_to_ZM(b);
     265       31073 :   return ZM_mul(m, SL2_inv(n));
     266             : }
     267             : 
     268             : /*************/
     269             : /* P^1(Z/NZ) */
     270             : /*************/
     271             : /* a != 0 in Z/NZ. Return v in (Z/NZ)^* such that av = gcd(a, N) (mod N)*/
     272             : static ulong
     273      322427 : Fl_inverse(ulong a, ulong N) { ulong g; return Fl_invgen(a,N,&g); }
     274             : 
     275             : /* Input: N = integer
     276             :  * Output: creates P^1(Z/NZ) = [symbols, H, N]
     277             :  *   symbols: list of vectors [x,y] that give a set of representatives
     278             :  *            of P^1(Z/NZ)
     279             :  *   H: an M by M grid whose value at the r,c-th place is the index of the
     280             :  *      "standard representative" equivalent to [r,c] occuring in the first
     281             :  *      list. If gcd(r,c,N) > 1 the grid has value 0. */
     282             : static GEN
     283        1218 : create_p1mod(ulong N)
     284             : {
     285        1218 :   GEN fa = factoru(N), div = divisorsu_fact(fa);
     286        1218 :   ulong i, nsym = count_Manin_symbols(N, gel(fa,1));
     287        1218 :   GEN symbols = generatemsymbols(N, nsym, div);
     288        1218 :   GEN H = inithashmsymbols(N,symbols);
     289        1218 :   GEN invsafe = cgetg(N, t_VECSMALL), inverse = cgetg(N, t_VECSMALL);
     290      123522 :   for (i = 1; i < N; i++)
     291             :   {
     292      122304 :     invsafe[i] = Fl_invsafe(i,N);
     293      122304 :     inverse[i] = Fl_inverse(i,N);
     294             :   }
     295        1218 :   return mkvecn(7, symbols, H, utoipos(N), fa, div, invsafe, inverse);
     296             : }
     297             : 
     298             : /* Let (c : d) in P1(Z/NZ).
     299             :  * If c = 0 return (0:1). If d = 0 return (1:0).
     300             :  * Else replace by (cu : du), where u in (Z/NZ)^* such that C := cu = gcd(c,N).
     301             :  * In create_p1mod(), (c : d) is represented by (C:D) where D = du (mod N/c)
     302             :  * is smallest such that gcd(C,D) = 1. Return (C : du mod N/c), which need
     303             :  * not belong to P1(Z/NZ) ! A second component du mod N/c = 0 is replaced by
     304             :  * N/c in this case to avoid problems with array indices */
     305             : static void
     306     2095345 : p1_std_form(long *pc, long *pd, GEN p1N)
     307             : {
     308     2095345 :   ulong N = p1N_get_N(p1N);
     309             :   ulong u;
     310     2095345 :   *pc = smodss(*pc, N); if (!*pc) { *pd = 1; return; }
     311     2030917 :   *pd = smodss(*pd, N); if (!*pd) { *pc = 1; return; }
     312     2002630 :   u = p1N_get_invsafe(p1N)[*pd];
     313     2002630 :   if (u) { *pc = Fl_mul(*pc,u,N); *pd = 1; return; } /* (d,N) = 1 */
     314             : 
     315      543676 :   u = p1N_get_inverse(p1N)[*pc];
     316      543676 :   if (u > 1) { *pc = Fl_mul(*pc,u,N); *pd = Fl_mul(*pd,u,N); }
     317             :   /* c | N */
     318      543676 :   if (*pc != 1) *pd %= (N / *pc);
     319      543676 :   if (!*pd) *pd = N / *pc;
     320             : }
     321             : 
     322             : /* Input: v = [x,y] = elt of P^1(Z/NZ) = class in Gamma_0(N) \ PSL2(Z)
     323             :  * Output: returns the index of the standard rep equivalent to v */
     324             : static long
     325     2095345 : p1_index(long x, long y, GEN p1N)
     326             : {
     327     2095345 :   ulong N = p1N_get_N(p1N);
     328     2095345 :   GEN H = p1N_get_hash(p1N);
     329             : 
     330     2095345 :   p1_std_form(&x, &y, p1N);
     331     2095345 :   if (y == 1) return x+1;
     332      571963 :   if (y == 0) return N+1;
     333      543676 :   if (mael(H,x,y) == 0) pari_err_BUG("p1_index");
     334      543676 :   return mael(H,x,y);
     335             : }
     336             : 
     337             : /* Cusps for \Gamma_0(N) */
     338             : 
     339             : /* \sum_{d | N} \phi(gcd(d, N/d)), using multiplicativity. fa = factor(N) */
     340             : ulong
     341        1239 : mfnumcuspsu_fact(GEN fa)
     342             : {
     343        1239 :   GEN P = gel(fa,1), E = gel(fa,2);
     344        1239 :   long i, l = lg(P);
     345        1239 :   ulong T = 1;
     346        3122 :   for (i = 1; i < l; i++)
     347             :   {
     348        1883 :     long e = E[i], e2 = e >> 1; /* floor(E[i] / 2) */
     349        1883 :     ulong p = P[i];
     350        1883 :     if (odd(e))
     351        1680 :       T *= 2 * upowuu(p, e2);
     352             :     else
     353         203 :       T *= (p+1) * upowuu(p, e2-1);
     354             :   }
     355        1239 :   return T;
     356             : }
     357             : ulong
     358           7 : mfnumcuspsu(ulong n)
     359             : {
     360           7 :   pari_sp av = avma;
     361           7 :   ulong t = mfnumcuspsu_fact( factoru(n) );
     362           7 :   avma = av; return t;
     363             : }
     364             : /* \sum_{d | N} \phi(gcd(d, N/d)), using multiplicativity. fa = factor(N) */
     365             : GEN
     366          14 : mfnumcusps_fact(GEN fa)
     367             : {
     368          14 :   GEN P = gel(fa,1), E = gel(fa,2), T = gen_1;
     369          14 :   long i, l = lg(P);
     370          35 :   for (i = 1; i < l; i++)
     371             :   {
     372          21 :     GEN p = gel(P,i), c;
     373          21 :     long e = itos(gel(E,i)), e2 = e >> 1; /* floor(E[i] / 2) */
     374          21 :     if (odd(e))
     375           0 :       c = shifti(powiu(p, e2), 1);
     376             :     else
     377          21 :       c = mulii(addiu(p,1), powiu(p, e2-1));
     378          21 :     T = T? mulii(T, c): c;
     379             :   }
     380          14 :   return T? T: gen_1;
     381             : }
     382             : GEN
     383          21 : mfnumcusps(GEN n)
     384             : {
     385          21 :   pari_sp av = avma;
     386          21 :   GEN F = check_arith_pos(n,"mfnumcusps");
     387          21 :   if (!F)
     388             :   {
     389          14 :     if (lgefint(n) == 3) return utoi( mfnumcuspsu(n[2]) );
     390           7 :     F = absZ_factor(n);
     391             :   }
     392          14 :   return gerepileuptoint(av, mfnumcusps_fact(F));
     393             : }
     394             : 
     395             : 
     396             : /* to each cusp in \Gamma_0(N) P1(Q), represented by p/q, we associate a
     397             :  * unique index. Canonical representative: (1:0) or (p:q) with q | N, q < N,
     398             :  * p defined modulo d := gcd(N/q,q), (p,d) = 1.
     399             :  * Return [[N, nbcusps], H, cusps]*/
     400             : static GEN
     401        1218 : inithashcusps(GEN p1N)
     402             : {
     403        1218 :   ulong N = p1N_get_N(p1N);
     404        1218 :   GEN div = p1N_get_div(p1N), H = zerovec(N+1);
     405        1218 :   long k, ind, l = lg(div), ncusp = mfnumcuspsu_fact(p1N_get_fa(p1N));
     406        1218 :   GEN cusps = cgetg(ncusp+1, t_VEC);
     407             : 
     408        1218 :   gel(H,1) = mkvecsmall2(0/*empty*/, 1/* first cusp: (1:0) */);
     409        1218 :   gel(cusps, 1) = mkvecsmall2(1,0);
     410        1218 :   ind = 2;
     411        4564 :   for (k=1; k < l-1; k++) /* l-1: remove q = N */
     412             :   {
     413        3346 :     ulong p, q = div[k], d = ugcd(q, N/q);
     414        3346 :     GEN h = const_vecsmall(d+1,0);
     415        3346 :     gel(H,q+1) = h ;
     416        8841 :     for (p = 0; p < d; p++)
     417        5495 :       if (ugcd(p,d) == 1)
     418             :       {
     419        4284 :         h[p+1] = ind;
     420        4284 :         gel(cusps, ind) = mkvecsmall2(p,q);
     421        4284 :         ind++;
     422             :       }
     423             :   }
     424        1218 :   return mkvec3(mkvecsmall2(N,ind-1), H, cusps);
     425             : }
     426             : /* c = [p,q], (p,q) = 1, return a canonical representative for
     427             :  * \Gamma_0(N)(p/q) */
     428             : static GEN
     429      201789 : cusp_std_form(GEN c, GEN S)
     430             : {
     431      201789 :   long p, N = gel(S,1)[1], q = smodss(c[2], N);
     432             :   ulong u, d;
     433      201789 :   if (q == 0) return mkvecsmall2(1, 0);
     434      200123 :   p = smodss(c[1], N);
     435      200123 :   u = Fl_inverse(q, N);
     436      200123 :   q = Fl_mul(q,u, N);
     437      200123 :   d = ugcd(q, N/q);
     438      200123 :   return mkvecsmall2(Fl_div(p % d,u % d, d), q);
     439             : }
     440             : /* c = [p,q], (p,q) = 1, return the index of the corresponding cusp.
     441             :  * S from inithashcusps */
     442             : static ulong
     443      201789 : cusp_index(GEN c, GEN S)
     444             : {
     445             :   long p, q;
     446      201789 :   GEN H = gel(S,2);
     447      201789 :   c = cusp_std_form(c, S);
     448      201789 :   p = c[1]; q = c[2];
     449      201789 :   if (!mael(H,q+1,p+1)) pari_err_BUG("cusp_index");
     450      201789 :   return mael(H,q+1,p+1);
     451             : }
     452             : 
     453             : /* M a square invertible ZM, return a ZM iM such that iM M = M iM = d.Id */
     454             : static GEN
     455        2821 : ZM_inv_denom(GEN M)
     456             : {
     457        2821 :   GEN diM, iM = ZM_inv_ratlift(M, &diM);
     458        2821 :   return mkvec2(iM, diM);
     459             : }
     460             : /* return M^(-1) v, dinv = ZM_inv_denom(M) OR Qevproj_init(M) */
     461             : static GEN
     462      735504 : ZC_apply_dinv(GEN dinv, GEN v)
     463             : {
     464             :   GEN x, c, iM;
     465      735504 :   if (lg(dinv) == 3)
     466             :   {
     467      658693 :     iM = gel(dinv,1);
     468      658693 :     c = gel(dinv,2);
     469             :   }
     470             :   else
     471             :   { /* Qevproj_init */
     472       76811 :     iM = gel(dinv,2);
     473       76811 :     c = gel(dinv,3);
     474      153622 :     v = typ(v) == t_MAT? rowpermute(v, gel(dinv,4))
     475       76811 :                        : vecpermute(v, gel(dinv,4));
     476             :   }
     477      735504 :   x = RgM_RgC_mul(iM, v);
     478      735504 :   if (!isint1(c)) x = RgC_Rg_div(x, c);
     479      735504 :   return x;
     480             : }
     481             : 
     482             : /* M an n x d ZM of rank d (basis of a Q-subspace), n >= d.
     483             :  * Initialize a projector on M */
     484             : GEN
     485        4165 : Qevproj_init(GEN M)
     486             : {
     487             :   GEN v, perm, MM, iM, diM;
     488        4165 :   v = ZM_indexrank(M); perm = gel(v,1);
     489        4165 :   MM = rowpermute(M, perm); /* square invertible */
     490        4165 :   iM = ZM_inv_ratlift(MM, &diM);
     491        4165 :   return mkvec4(M, iM, diM, perm);
     492             : }
     493             : static int
     494         161 : is_Qevproj(GEN x)
     495         161 : { return typ(x) == t_VEC && lg(x) == 5 && typ(gel(x,1)) == t_MAT; }
     496             : 
     497             : /* same with typechecks */
     498             : static GEN
     499         700 : Qevproj_init0(GEN M)
     500             : {
     501         700 :   switch(typ(M))
     502             :   {
     503             :     case t_VEC:
     504         651 :       if (lg(M) == 5) return M;
     505           0 :       break;
     506             :     case t_COL:
     507          42 :       M = mkmat(M);/*fall through*/
     508             :     case t_MAT:
     509          49 :       M = Q_primpart(M);
     510          49 :       RgM_check_ZM(M,"Qevproj_init");
     511          49 :       return Qevproj_init(M);
     512             :   }
     513           0 :   pari_err_TYPE("Qevproj_init",M);
     514           0 :   return NULL;
     515             : }
     516             : 
     517             : /* T an n x n QM, stabilizing d-dimensional Q-vector space spanned by the
     518             :  * columns of M, pro = Qevproj_init(M). Return d x d matrix of T acting
     519             :  * on M */
     520             : GEN
     521        3045 : Qevproj_apply(GEN T, GEN pro)
     522             : {
     523        3045 :   GEN M = gel(pro,1), iM = gel(pro,2), ciM = gel(pro,3), perm = gel(pro,4);
     524        3045 :   return RgM_Rg_div(RgM_mul(iM, RgM_mul(rowpermute(T,perm), M)), ciM);
     525             : }
     526             : /* Qevproj_apply(T,pro)[,k] */
     527             : GEN
     528         777 : Qevproj_apply_vecei(GEN T, GEN pro, long k)
     529             : {
     530         777 :   GEN M = gel(pro,1), iM = gel(pro,2), ciM = gel(pro,3), perm = gel(pro,4);
     531         777 :   GEN v = RgM_RgC_mul(iM, RgM_RgC_mul(rowpermute(T,perm), gel(M,k)));
     532         777 :   return RgC_Rg_div(v, ciM);
     533             : }
     534             : 
     535             : static GEN
     536        1162 : QM_ker_r(GEN M) { return ZM_ker(Q_primpart(M)); }
     537             : static GEN
     538         840 : QM_image(GEN A)
     539             : {
     540         840 :   A = vec_Q_primpart(A);
     541         840 :   return vecpermute(A, ZM_indeximage(A));
     542             : }
     543             : 
     544             : static int
     545         420 : cmp_dim(void *E, GEN a, GEN b)
     546             : {
     547             :   long k;
     548             :   (void)E;
     549         420 :   a = gel(a,1);
     550         420 :   b = gel(b,1); k = lg(a)-lg(b);
     551         420 :   return k? ((k > 0)? 1: -1): 0;
     552             : }
     553             : 
     554             : /* FIXME: could use ZX_roots for deglim = 1 */
     555             : static GEN
     556         322 : ZX_factor_limit(GEN T, long deglim, long *pl)
     557             : {
     558         322 :   GEN fa = ZX_factor(T), P, E;
     559             :   long i, l;
     560         322 :   P = gel(fa,1); *pl = l = lg(P);
     561         322 :   if (deglim <= 0) return fa;
     562         224 :   E = gel(fa,2);
     563         567 :   for (i = 1; i < l; i++)
     564         406 :     if (degpol(gel(P,i)) > deglim) break;
     565         224 :   setlg(P,i);
     566         224 :   setlg(E,i); return fa;
     567             : }
     568             : 
     569             : /* Decompose the subspace H (Qevproj format) in simple subspaces.
     570             :  * Eg for H = msnew */
     571             : static GEN
     572         252 : mssplit_i(GEN W, GEN H, long deglim)
     573             : {
     574         252 :   ulong p, N = ms_get_N(W);
     575             :   long first, dim;
     576             :   forprime_t S;
     577         252 :   GEN T1 = NULL, T2 = NULL, V;
     578         252 :   dim = lg(gel(H,1))-1;
     579         252 :   V = vectrunc_init(dim+1);
     580         252 :   if (!dim) return V;
     581         245 :   (void)u_forprime_init(&S, 2, ULONG_MAX);
     582         245 :   vectrunc_append(V, H);
     583         245 :   first = 1; /* V[1..first-1] contains simple subspaces */
     584         616 :   while ((p = u_forprime_next(&S)))
     585             :   {
     586             :     GEN T;
     587             :     long j, lV;
     588         371 :     if (N % p == 0) continue;
     589         315 :     if (T1 && T2) {
     590          21 :       T = RgM_add(T1,T2);
     591          21 :       T2 = NULL;
     592             :     } else {
     593         294 :       T2 = T1;
     594         294 :       T1 = T = mshecke(W, p, NULL);
     595             :     }
     596         315 :     lV = lg(V);
     597         637 :     for (j = first; j < lV; j++)
     598             :     {
     599         322 :       pari_sp av = avma;
     600             :       long lP;
     601         322 :       GEN Vj = gel(V,j), P = gel(Vj,1);
     602         322 :       GEN TVj = Qevproj_apply(T, Vj); /* c T | V_j */
     603         322 :       GEN ch = QM_charpoly_ZX(TVj), fa = ZX_factor_limit(ch,deglim, &lP);
     604         322 :       GEN F = gel(fa, 1), E = gel(fa, 2);
     605         322 :       long k, lF = lg(F);
     606         322 :       if (lF == 2 && lP == 2)
     607             :       {
     608         322 :         if (isint1(gel(E,1)))
     609             :         { /* simple subspace */
     610         161 :           swap(gel(V,first), gel(V,j));
     611         161 :           first++;
     612             :         }
     613             :         else
     614           0 :           avma = av;
     615             :       }
     616         161 :       else if (lF == 1) /* discard V[j] */
     617           7 :       { swap(gel(V,j), gel(V,lg(V)-1)); setlg(V, lg(V)-1); }
     618             :       else
     619             :       { /* can split Vj */
     620             :         GEN pows;
     621         154 :         long D = 1;
     622         616 :         for (k = 1; k < lF; k++)
     623             :         {
     624         462 :           long d = degpol(gel(F,k));
     625         462 :           if (d > D) D = d;
     626             :         }
     627             :         /* remove V[j] */
     628         154 :         swap(gel(V,j), gel(V,lg(V)-1)); setlg(V, lg(V)-1);
     629         154 :         pows = RgM_powers(TVj, minss((long)2*sqrt((double)D), D));
     630         616 :         for (k = 1; k < lF; k++)
     631             :         {
     632         462 :           GEN f = gel(F,k);
     633         462 :           GEN K = QM_ker_r( RgX_RgMV_eval(f, pows)) ; /* Ker f(TVj) */
     634         462 :           GEN p = vec_Q_primpart( RgM_mul(P, K) );
     635         462 :           vectrunc_append(V, Qevproj_init(p));
     636         462 :           if (lg(K) == 2 || isint1(gel(E,k)))
     637             :           { /* simple subspace */
     638         385 :             swap(gel(V,first), gel(V, lg(V)-1));
     639         385 :             first++;
     640             :           }
     641             :         }
     642         154 :         if (j < first) j = first;
     643             :       }
     644             :     }
     645         315 :     if (first >= lg(V)) {
     646         245 :       gen_sort_inplace(V, NULL, cmp_dim, NULL);
     647         245 :       return V;
     648             :     }
     649             :   }
     650           0 :   pari_err_BUG("subspaces not found");
     651           0 :   return NULL;
     652             : }
     653             : GEN
     654         252 : mssplit(GEN W, GEN H, long deglim)
     655             : {
     656         252 :   pari_sp av = avma;
     657         252 :   checkms(W);
     658         252 :   if (!msk_get_sign(W))
     659           0 :     pari_err_DOMAIN("mssplit","abs(sign)","!=",gen_1,gen_0);
     660         252 :   H = Qevproj_init0(H);
     661         252 :   return gerepilecopy(av, mssplit_i(W,H,deglim));
     662             : }
     663             : 
     664             : /* proV = Qevproj_init of a Hecke simple subspace, return [ a_n, n <= B ] */
     665             : static GEN
     666         238 : msqexpansion_i(GEN W, GEN proV, ulong B)
     667             : {
     668         238 :   ulong p, N = ms_get_N(W), sqrtB;
     669         238 :   long i, d, k = msk_get_weight(W);
     670             :   forprime_t S;
     671         238 :   GEN T1=NULL, T2=NULL, TV=NULL, ch=NULL, v, dTiv, Tiv, diM, iM, L;
     672         238 :   switch(B)
     673             :   {
     674           0 :     case 0: return cgetg(1,t_VEC);
     675           0 :     case 1: return mkvec(gen_1);
     676             :   }
     677         238 :   (void)u_forprime_init(&S, 2, ULONG_MAX);
     678         588 :   while ((p = u_forprime_next(&S)))
     679             :   {
     680             :     GEN T;
     681         350 :     if (N % p == 0) continue;
     682         259 :     if (T1 && T2)
     683             :     {
     684           0 :       T = RgM_add(T1,T2);
     685           0 :       T2 = NULL;
     686             :     }
     687             :     else
     688             :     {
     689         259 :       T2 = T1;
     690         259 :       T1 = T = mshecke(W, p, NULL);
     691             :     }
     692         259 :     TV = Qevproj_apply(T, proV); /* T | V */
     693         259 :     ch = QM_charpoly_ZX(TV);
     694         259 :     if (ZX_is_irred(ch)) break;
     695          21 :     ch = NULL;
     696             :   }
     697         238 :   if (!ch) pari_err_BUG("q-Expansion not found");
     698             :   /* T generates the Hecke algebra (acting on V) */
     699         238 :   d = degpol(ch);
     700         238 :   v = vec_ei(d, 1); /* take v = e_1 */
     701         238 :   Tiv = cgetg(d+1, t_MAT); /* Tiv[i] = T^(i-1)v */
     702         238 :   gel(Tiv, 1) = v;
     703         238 :   for (i = 2; i <= d; i++) gel(Tiv, i) = RgM_RgC_mul(TV, gel(Tiv,i-1));
     704         238 :   Tiv = Q_remove_denom(Tiv, &dTiv);
     705         238 :   iM = ZM_inv_ratlift(Tiv, &diM);
     706         238 :   if (dTiv) diM = gdiv(diM, dTiv);
     707         238 :   L = const_vec(B,NULL);
     708         238 :   sqrtB = usqrt(B);
     709         238 :   gel(L,1) = d > 1? mkpolmod(gen_1,ch): gen_1;
     710        2359 :   for (p = 2; p <= B; p++)
     711             :   {
     712        2121 :     pari_sp av = avma;
     713             :     GEN T, u, Tv, ap, P;
     714             :     ulong m;
     715        2121 :     if (gel(L,p)) continue;  /* p not prime */
     716         777 :     T = mshecke(W, p, NULL);
     717         777 :     Tv = Qevproj_apply_vecei(T, proV, 1); /* Tp.v */
     718             :     /* Write Tp.v = \sum u_i T^i v */
     719         777 :     u = RgC_Rg_div(RgM_RgC_mul(iM, Tv), diM);
     720         777 :     ap = gerepilecopy(av, RgV_to_RgX(u, 0));
     721         777 :     if (d > 1)
     722         399 :       ap = mkpolmod(ap,ch);
     723             :     else
     724         378 :       ap = simplify_shallow(ap);
     725         777 :     gel(L,p) = ap;
     726         777 :     if (!(N % p))
     727             :     { /* p divides the level */
     728         147 :       ulong C = B/p;
     729         546 :       for (m=1; m<=C; m++)
     730         399 :         if (gel(L,m)) gel(L,m*p) = gmul(gel(L,m), ap);
     731         147 :       continue;
     732             :     }
     733         630 :     P = powuu(p,k-1);
     734         630 :     if (p <= sqrtB) {
     735         105 :       ulong pj, oldpj = 1;
     736         490 :       for (pj = p; pj <= B; oldpj=pj, pj *= p)
     737             :       {
     738         385 :         GEN apj = (pj==p)? ap
     739         385 :                          : gsub(gmul(ap,gel(L,oldpj)), gmul(P,gel(L,oldpj/p)));
     740         385 :         gel(L,pj) = apj;
     741        2989 :         for (m = B/pj; m > 1; m--)
     742        2604 :           if (gel(L,m) && m%p) gel(L,m*pj) = gmul(gel(L,m), apj);
     743             :       }
     744             :     } else {
     745         525 :       gel(L,p) = ap;
     746        1043 :       for (m = B/p; m > 1; m--)
     747         518 :         if (gel(L,m)) gel(L,m*p) = gmul(gel(L,m), ap);
     748             :     }
     749             :   }
     750         238 :   return L;
     751             : }
     752             : GEN
     753         238 : msqexpansion(GEN W, GEN proV, ulong B)
     754             : {
     755         238 :   pari_sp av = avma;
     756         238 :   checkms(W);
     757         238 :   proV = Qevproj_init0(proV);
     758         238 :   return gerepilecopy(av, msqexpansion_i(W,proV,B));
     759             : }
     760             : 
     761             : static GEN
     762         602 : Qevproj_apply2(GEN T, GEN pro1, GEN pro2)
     763             : {
     764         602 :   GEN M = gel(pro1,1), iM = gel(pro2,2), ciM = gel(pro2,3), perm = gel(pro2,4);
     765         602 :   return RgM_Rg_div(RgM_mul(iM, RgM_mul(rowpermute(T,perm), M)), ciM);
     766             : }
     767             : static GEN
     768         259 : Qevproj_apply0(GEN T, GEN pro)
     769             : {
     770         259 :   GEN iM = gel(pro,2), perm = gel(pro,4);
     771         259 :   return vec_Q_primpart(ZM_mul(iM, rowpermute(T,perm)));
     772             : }
     773             : /* T a ZC or ZM */
     774             : GEN
     775        1057 : Qevproj_down(GEN T, GEN pro)
     776             : {
     777        1057 :   GEN iM = gel(pro,2), ciM = gel(pro,3), perm = gel(pro,4);
     778        1057 :   if (typ(T) == t_COL)
     779        1057 :     return RgC_Rg_div(ZM_ZC_mul(iM, vecpermute(T,perm)), ciM);
     780             :   else
     781           0 :     return RgM_Rg_div(ZM_mul(iM, rowpermute(T,perm)), ciM);
     782             : }
     783             : 
     784             : static GEN
     785         294 : Qevproj_star(GEN W, GEN H)
     786             : {
     787         294 :   long s = msk_get_sign(W);
     788         294 :   if (s)
     789             :   { /* project on +/- component */
     790         259 :     GEN A = RgM_mul(msk_get_star(W), H);
     791         259 :     A = (s > 0)? gadd(A, H): gsub(A, H);
     792             :     /* Im(star + sign) = Ker(star - sign) */
     793         259 :     H = QM_image(A);
     794         259 :     H = Qevproj_apply0(H, msk_get_starproj(W));
     795             :   }
     796         294 :   return H;
     797             : }
     798             : 
     799             : static GEN
     800        2387 : Tp_matrices(ulong p)
     801             : {
     802        2387 :   GEN v = cgetg(p+2, t_VEC);
     803             :   ulong i;
     804        2387 :   for (i = 1; i <= p; i++) gel(v,i) = mat2(1, i-1, 0, p);
     805        2387 :   gel(v,i) = mat2(p, 0, 0, 1);
     806        2387 :   return v;
     807             : }
     808             : static GEN
     809         924 : Up_matrices(ulong p)
     810             : {
     811         924 :   GEN v = cgetg(p+1, t_VEC);
     812             :   ulong i;
     813         924 :   for (i = 1; i <= p; i++) gel(v,i) = mat2(1, i-1, 0, p);
     814         924 :   return v;
     815             : }
     816             : 
     817             : /* M = N/p. Classes of Gamma_0(M) / Gamma_O(N) when p | M */
     818             : static GEN
     819         168 : NP_matrices(ulong M, ulong p)
     820             : {
     821         168 :   GEN v = cgetg(p+1, t_VEC);
     822             :   ulong i;
     823         168 :   for (i = 1; i <= p; i++) gel(v,i) = mat2(1, 0, (i-1)*M, 1);
     824         168 :   return v;
     825             : }
     826             : /* M = N/p. Extra class of Gamma_0(M) / Gamma_O(N) when p \nmid M */
     827             : static GEN
     828          84 : NP_matrix_extra(ulong M, ulong p)
     829             : {
     830          84 :   long w,z, d = cbezout(p, -M, &w, &z);
     831          84 :   if (d != 1) return NULL;
     832          84 :   return mat2(w,z,M,p);
     833             : }
     834             : static GEN
     835          98 : WQ_matrix(long N, long Q)
     836             : {
     837          98 :   long w,z, d = cbezout(Q, N/Q, &w, &z);
     838          98 :   if (d != 1) return NULL;
     839          98 :   return mat2(Q,1,-N*z,Q*w);
     840             : }
     841             : 
     842             : GEN
     843         266 : msnew(GEN W)
     844             : {
     845         266 :   pari_sp av = avma;
     846         266 :   GEN S = mscuspidal(W, 0);
     847         266 :   ulong N = ms_get_N(W);
     848         266 :   long s = msk_get_sign(W);
     849         266 :   if (!uisprime(N))
     850             :   {
     851         105 :     GEN p1N = ms_get_p1N(W), P = gel(p1N_get_fa(p1N), 1);
     852         105 :     long i, nP = lg(P)-1, k = msk_get_weight(W);
     853         105 :     GEN v = cgetg(2*nP + 1, t_COL);
     854         105 :     S = gel(S,1); /* Q basis */
     855         273 :     for (i = 1; i <= nP; i++)
     856             :     {
     857         168 :       pari_sp av = avma, av2;
     858         168 :       long M = N/P[i];
     859         168 :       GEN T1,Td, Wi = mskinit(M, k, s);
     860         168 :       GEN v1 = NP_matrices(M, P[i]);
     861         168 :       GEN vd = Up_matrices(P[i]);
     862             :       /* p^2 \nmid N */
     863         168 :       if (M % P[i])
     864             :       {
     865          84 :         v1 = shallowconcat(v1, mkvec(NP_matrix_extra(M,P[i])));
     866          84 :         vd = shallowconcat(vd, mkvec(WQ_matrix(N,P[i])));
     867             :       }
     868         168 :       T1 = getMorphism(W, Wi, v1);
     869         168 :       Td = getMorphism(W, Wi, vd);
     870         168 :       if (s)
     871             :       {
     872         154 :         T1 = Qevproj_apply2(T1, msk_get_starproj(W), msk_get_starproj(Wi));
     873         154 :         Td = Qevproj_apply2(Td, msk_get_starproj(W), msk_get_starproj(Wi));
     874             :       }
     875         168 :       av2 = avma;
     876         168 :       T1 = RgM_mul(T1,S);
     877         168 :       Td = RgM_mul(Td,S);  /* multiply by S = restrict to mscusp */
     878         168 :       gerepileallsp(av, av2, 2, &T1, &Td);
     879         168 :       gel(v,2*i-1) = T1;
     880         168 :       gel(v,2*i)   = Td;
     881             :     }
     882         105 :     S = ZM_mul(S, QM_ker_r(matconcat(v))); /* Snew */
     883         105 :     S = Qevproj_init(vec_Q_primpart(S));
     884             :   }
     885         266 :   return gerepilecopy(av, S);
     886             : }
     887             : 
     888             : /* Solve the Manin relations for a congruence subgroup \Gamma by constructing
     889             :  * a well-formed fundamental domain for the action of \Gamma on upper half
     890             :  * space. See
     891             :  * Pollack and Stevens, Overconvergent modular symbols and p-adic L-functions
     892             :  * Annales scientifiques de l'ENS 44, fascicule 1 (2011), 1-42
     893             :  * http://math.bu.edu/people/rpollack/Papers/Overconvergent_modular_symbols_and_padic_Lfunctions.pdf
     894             :  *
     895             :  * FIXME: Implemented for \Gamma = \Gamma_0(N) only. */
     896             : 
     897             : #if 0 /* Pollack-Stevens shift their paths so as to solve equations of the
     898             :          form f(z+1) - f(z) = g. We don't (to avoid mistakes) so we will
     899             :          have to solve eqs of the form f(z-1) - f(z) = g */
     900             : /* c = a/b; as a t_VECSMALL [a,b]; return c-1 as a t_VECSMALL */
     901             : static GEN
     902             : Shift_left_cusp(GEN c) { long a=c[1], b=c[2]; return mkvecsmall2(a - b, b); }
     903             : /* c = a/b; as a t_VECSMALL [a,b]; return c+1 as a t_VECSMALL */
     904             : static GEN
     905             : Shift_right_cusp(GEN c) { long a=c[1], b=c[2]; return mkvecsmall2(a + b, b); }
     906             : /*Input: path = [r,s] (thought of as a geodesic between these points)
     907             :  *Output: The path shifted by one to the left, i.e. [r-1,s-1] */
     908             : static GEN
     909             : Shift_left(GEN path)
     910             : {
     911             :   GEN r = gel(path,1), s = gel(path,2);
     912             :   return mkvec2(Shift_left_cusp(r), Shift_left_cusp(s)); }
     913             : /*Input: path = [r,s] (thought of as a geodesic between these points)
     914             :  *Output: The path shifted by one to the right, i.e. [r+1,s+1] */
     915             : GEN
     916             : Shift_right(GEN path)
     917             : {
     918             :   GEN r = gel(path,1), s = gel(path,2);
     919             :   return mkvec2(Shift_right_cusp(r), Shift_right_cusp(s)); }
     920             : #endif
     921             : 
     922             : /* linked lists */
     923             : typedef struct list_t { GEN data; struct list_t *next; } list_t;
     924             : static list_t *
     925       60718 : list_new(GEN x)
     926             : {
     927       60718 :   list_t *L = (list_t*)stack_malloc(sizeof(list_t));
     928       60718 :   L->data = x;
     929       60718 :   L->next = NULL; return L;
     930             : }
     931             : static void
     932       59500 : list_insert(list_t *L, GEN x)
     933             : {
     934       59500 :   list_t *l = list_new(x);
     935       59500 :   l->next = L->next;
     936       59500 :   L->next = l;
     937       59500 : }
     938             : 
     939             : /*Input: N > 1, p1N = P^1(Z/NZ)
     940             :  *Output: a connected fundamental domain for the action of \Gamma_0(N) on
     941             :  *  upper half space.  When \Gamma_0(N) is torsion free, the domain has the
     942             :  *  property that all of its vertices are cusps.  When \Gamma_0(N) has
     943             :  *  three-torsion, 2 extra triangles need to be added.
     944             :  *
     945             :  * The domain is constructed by beginning with the triangle with vertices 0,1
     946             :  * and oo.  Each adjacent triangle is successively tested to see if it contains
     947             :  * points not \Gamma_0(N) equivalent to some point in our region.  If a
     948             :  * triangle contains new points, it is added to the region.  This process is
     949             :  * continued until the region can no longer be extended (and still be a
     950             :  * fundamental domain) by added an adjacent triangle.  The list of cusps
     951             :  * between 0 and 1 are then returned
     952             :  *
     953             :  * Precisely, the function returns a list such that the elements of the list
     954             :  * with odd index are the cusps in increasing order.  The even elements of the
     955             :  * list are either an "x" or a "t".  A "t" represents that there is an element
     956             :  * of order three such that its fixed point is in the triangle directly
     957             :  * adjacent to the our region with vertices given by the cusp before and after
     958             :  * the "t".  The "x" represents that this is not the case. */
     959             : enum { type_X, type_DO /* ? */, type_T };
     960             : static GEN
     961        1218 : form_list_of_cusps(ulong N, GEN p1N)
     962             : {
     963        1218 :   pari_sp av = avma;
     964        1218 :   long i, position, nbC = 2;
     965             :   GEN v, L;
     966             :   list_t *C, *c;
     967             :   /* Let t be the index of a class in PSL2(Z) / \Gamma in our fixed enumeration
     968             :    * v[t] != 0 iff it is the class of z tau^r for z a previous alpha_i
     969             :    * or beta_i.
     970             :    * For \Gamma = \Gamma_0(N), the enumeration is given by p1_index.
     971             :    * We write cl(gamma) = the class of gamma mod \Gamma */
     972        1218 :   v = const_vecsmall(p1_size(p1N), 0);
     973        1218 :   i = p1_index( 0, 1, p1N); v[i] = 1;
     974        1218 :   i = p1_index( 1,-1, p1N); v[i] = 2;
     975        1218 :   i = p1_index(-1, 0, p1N); v[i] = 3;
     976             :   /* the value is unused [debugging]: what matters is whether it is != 0 */
     977        1218 :   position = 4;
     978             :   /* at this point, Fund = R, v contains the classes of Id, tau, tau^2 */
     979             : 
     980        1218 :   C  = list_new(mkvecsmall3(0,1, type_X));
     981        1218 :   list_insert(C, mkvecsmall3(1,1,type_DO));
     982             :   /* C is a list of triples[a,b,t], where c = a/b is a cusp, and t is the type
     983             :    * of the path between c and the PREVIOUS cusp in the list, coded as
     984             :    *   type_DO = "?", type_X = "x", type_T = "t"
     985             :    * Initially, C = [0/1,"?",1/1]; */
     986             : 
     987             :   /* loop through the current set of cusps C and check to see if more cusps
     988             :    * should be added */
     989             :   for (;;)
     990             :   {
     991        6727 :     int done = 1;
     992      296884 :     for (c = C; c; c = c->next)
     993             :     {
     994             :       GEN cusp1, cusp2, gam;
     995             :       long pos, b1, b2, b;
     996             : 
     997      296884 :       if (!c->next) break;
     998      290157 :       cusp1 = c->data; /* = a1/b1 */
     999      290157 :       cusp2 = (c->next)->data; /* = a2/b2 */
    1000      290157 :       if (cusp2[3] != type_DO) continue;
    1001             : 
    1002             :       /* gam (oo -> 0) = (cusp2 -> cusp1), gam in PSL2(Z) */
    1003      117782 :       gam = path_to_zm(mkpath(cusp2, cusp1)); /* = [a2,a1;b2,b1] */
    1004             :       /* we have normalized the cusp representation so that a1 b2 - a2 b1 = 1 */
    1005      117782 :       b1 = coeff(gam,2,1); b2 = coeff(gam,2,2);
    1006             :       /* gam.1  = (a1 + a2) / (b1 + b2) */
    1007      117782 :       b = b1 + b2;
    1008             :       /* Determine whether the adjacent triangle *below* (cusp1->cusp2)
    1009             :        * should be added */
    1010      117782 :       pos = p1_index(b1,b2, p1N); /* did we see cl(gam) before ? */
    1011      117782 :       if (v[pos])
    1012       59087 :         cusp2[3] = type_X; /* NO */
    1013             :       else
    1014             :       { /* YES */
    1015             :         ulong B1, B2;
    1016       58695 :         v[pos] = position;
    1017       58695 :         i = p1_index(-(b1+b2), b1, p1N); v[i] = position+1;
    1018       58695 :         i = p1_index(b2, -(b1+b2), p1N); v[i] = position+2;
    1019             :         /* add cl(gam), cl(gam*TAU), cl(gam*TAU^2) to v */
    1020       58695 :         position += 3;
    1021             :         /* gam tau gam^(-1) in \Gamma ? */
    1022       58695 :         B1 = smodss(b1, N);
    1023       58695 :         B2 = smodss(b2, N);
    1024       58695 :         if ((Fl_sqr(B2,N) + Fl_sqr(B1,N) + Fl_mul(B1,B2,N)) % N == 0)
    1025         413 :           cusp2[3] = type_T;
    1026             :         else
    1027             :         {
    1028       58282 :           long a1 = coeff(gam, 1,1), a2 = coeff(gam, 1,2);
    1029       58282 :           long a = a1 + a2; /* gcd(a,b) = 1 */
    1030       58282 :           list_insert(c, mkvecsmall3(a,b,type_DO));
    1031       58282 :           c = c->next;
    1032       58282 :           nbC++;
    1033       58282 :           done = 0;
    1034             :         }
    1035             :       }
    1036             :     }
    1037        6727 :     if (done) break;
    1038        5509 :   }
    1039        1218 :   L = cgetg(nbC+1, t_VEC); i = 1;
    1040        1218 :   for (c = C; c; c = c->next) gel(L,i++) = c->data;
    1041        1218 :   return gerepilecopy(av, L);
    1042             : }
    1043             : 
    1044             : /* M in PSL2(Z). Return index of M in P1^(Z/NZ) = Gamma0(N) \ PSL2(Z),
    1045             :  * and M0 in Gamma_0(N) such that M = M0 * M', where M' = chosen
    1046             :  * section( PSL2(Z) -> P1^(Z/NZ) ). */
    1047             : static GEN
    1048      446698 : Gamma0N_decompose(GEN W, GEN M, long *index)
    1049             : {
    1050      446698 :   GEN p1N = gel(W,1), W3 = gel(W,3), section = ms_get_section(W);
    1051             :   GEN A;
    1052      446698 :   ulong N = p1N_get_N(p1N);
    1053      446698 :   ulong c = umodiu(gcoeff(M,2,1), N);
    1054      446698 :   ulong d = umodiu(gcoeff(M,2,2), N);
    1055      446698 :   long s, ind = p1_index(c, d, p1N); /* as an elt of P1(Z/NZ) */
    1056      446698 :   *index = W3[ind]; /* as an elt of F, E2, ... */
    1057      446698 :   M = ZM_zm_mul(M, sl2_inv(gel(section,ind)));
    1058             :   /* normalize mod +/-Id */
    1059      446698 :   A = gcoeff(M,1,1);
    1060      446698 :   s = signe(A);
    1061      446698 :   if (s < 0)
    1062      221459 :     M = ZM_neg(M);
    1063      225239 :   else if (!s)
    1064             :   {
    1065           0 :     GEN C = gcoeff(M,2,1);
    1066           0 :     if (signe(C) < 0) M = ZM_neg(M);
    1067             :   }
    1068      446698 :   return M;
    1069             : }
    1070             : /* same for a path. Return [[ind], M] */
    1071             : static GEN
    1072      123872 : path_Gamma0N_decompose(GEN W, GEN path)
    1073             : {
    1074      123872 :   GEN p1N = gel(W,1);
    1075      123872 :   GEN p1index_to_ind = gel(W,3);
    1076      123872 :   GEN section = ms_get_section(W);
    1077      123872 :   GEN M = path_to_zm(path);
    1078      123872 :   long p1index = p1_index(cc(M), dd(M), p1N);
    1079      123872 :   long ind = p1index_to_ind[p1index];
    1080      123872 :   GEN M0 = ZM_zm_mul(mat2_to_ZM(M), sl2_inv(gel(section,p1index)));
    1081      123872 :   return mkvec2(mkvecsmall(ind), M0);
    1082             : }
    1083             : 
    1084             : /*Form generators of H_1(X_0(N),{cusps},Z)
    1085             : *
    1086             : *Input: N = integer > 1, p1N = P^1(Z/NZ)
    1087             : *Output: [cusp_list,E,F,T2,T3,E1] where
    1088             : *  cusps_list = list of cusps describing fundamental domain of
    1089             : *    \Gamma_0(N).
    1090             : *  E = list of paths in the boundary of the fundamental domains and oriented
    1091             : *    clockwise such that they do not contain a point
    1092             : *    fixed by an element of order 2 and they are not an edge of a
    1093             : *    triangle containing a fixed point of an element of order 3
    1094             : *  F = list of paths in the interior of the domain with each
    1095             : *    orientation appearing separately
    1096             : * T2 = list of paths in the boundary of domain containing a point fixed
    1097             : *    by an element of order 2 (oriented clockwise)
    1098             : * T3 = list of paths in the boundard of domain which are the edges of
    1099             : *    some triangle containing a fixed point of a matrix of order 3 (both
    1100             : *    orientations appear)
    1101             : * E1 = a sublist of E such that every path in E is \Gamma_0(N)-equivalent to
    1102             : *    either an element of E1 or the flip (reversed orientation) of an element
    1103             : *    of E1.
    1104             : * (Elements of T2 are \Gamma_0(N)-equivalent to their own flip.)
    1105             : *
    1106             : * sec = a list from 1..#p1N of matrices describing a section of the map
    1107             : *   SL_2(Z) to P^1(Z/NZ) given by [a,b;c,d]-->[c,d].
    1108             : *   Given our fixed enumeration of P^1(Z/NZ), the j-th element of the list
    1109             : *   represents the image of the j-th element of P^1(Z/NZ) under the section. */
    1110             : 
    1111             : /* insert path in set T */
    1112             : static void
    1113      178913 : set_insert(hashtable *T, GEN path)
    1114      178913 : { hash_insert(T, path,  (void*)(T->nb + 1)); }
    1115             : 
    1116             : static GEN
    1117       10962 : hash_to_vec(hashtable *h)
    1118             : {
    1119       10962 :   GEN v = cgetg(h->nb + 1, t_VEC);
    1120             :   ulong i;
    1121     1484756 :   for (i = 0; i < h->len; i++)
    1122             :   {
    1123     1473794 :     hashentry *e = h->table[i];
    1124     3244101 :     while (e)
    1125             :     {
    1126      296513 :       GEN key = (GEN)e->key;
    1127      296513 :       long index = (long)e->val;
    1128      296513 :       gel(v, index) = key;
    1129      296513 :       e = e->next;
    1130             :     }
    1131             :   }
    1132       10962 :   return v;
    1133             : }
    1134             : 
    1135             : static long
    1136       91350 : path_to_p1_index(GEN path, GEN p1N)
    1137             : {
    1138       91350 :   GEN M = path_to_zm(path);
    1139       91350 :   return p1_index(cc(M), dd(M), p1N);
    1140             : }
    1141             : 
    1142             : /* Pollack-Stevens sets */
    1143             : typedef struct PS_sets_t {
    1144             :   hashtable *F, *T2, *T31, *T32, *E1, *E2;
    1145             :   GEN E2_in_terms_of_E1, stdE1;
    1146             : } PS_sets_t;
    1147             : 
    1148             : static hashtable *
    1149       10787 : set_init(long max)
    1150       10787 : { return hash_create(max, (ulong(*)(void*))&hash_GEN,
    1151             :                           (int(*)(void*,void*))&gidentical, 1); }
    1152             : static void
    1153       60900 : insert_E(GEN path, PS_sets_t *S, GEN p1N)
    1154             : {
    1155       60900 :   GEN rev = vecreverse(path);
    1156       60900 :   long std = path_to_p1_index(rev, p1N);
    1157       60900 :   GEN v = gel(S->stdE1, std);
    1158       60900 :   if (v)
    1159             :   { /* [s, p1], where E1[s] = the path p1 \equiv vecreverse(path) mod \Gamma */
    1160       30450 :     GEN gamma, p1 = gel(v,2);
    1161       30450 :     long r, s = itos(gel(v,1));
    1162             : 
    1163       30450 :     set_insert(S->E2, path);
    1164       30450 :     r = S->E2->nb;
    1165       30450 :     if (gel(S->E2_in_terms_of_E1, r) != gen_0) pari_err_BUG("insert_E");
    1166             : 
    1167       30450 :     gamma = gamma_equiv_matrix(rev, p1);
    1168             :     /* E2[r] + gamma * E1[s] = 0 */
    1169       30450 :     gel(S->E2_in_terms_of_E1, r) = mkvec2(utoipos(s),
    1170             :                                           to_famat_shallow(gamma,gen_m1));
    1171             :   }
    1172             :   else
    1173             :   {
    1174       30450 :     set_insert(S->E1, path);
    1175       30450 :     std = path_to_p1_index(path, p1N);
    1176       30450 :     gel(S->stdE1, std) = mkvec2(utoipos(S->E1->nb), path);
    1177             :   }
    1178       60900 : }
    1179             : 
    1180             : static GEN
    1181        4872 : cusp_infinity(void) { return mkvecsmall2(1,0); }
    1182             : 
    1183             : static void
    1184        1218 : form_E_F_T(ulong N, GEN p1N, GEN *pC, PS_sets_t *S)
    1185             : {
    1186        1218 :   GEN C, cusp_list = form_list_of_cusps(N, p1N);
    1187        1218 :   long nbgen = lg(cusp_list)-1, nbmanin = p1_size(p1N), r, s, i;
    1188             :   hashtable *F, *T2, *T31, *T32, *E1, *E2;
    1189             : 
    1190        1218 :   *pC = C = cgetg(nbgen+1, t_VEC);
    1191       61936 :   for (i = 1; i <= nbgen; i++)
    1192             :   {
    1193       60718 :     GEN c = gel(cusp_list,i);
    1194       60718 :     gel(C,i) = mkvecsmall2(c[1], c[2]);
    1195             :   }
    1196        1218 :   S->F  = F  = set_init(nbmanin);
    1197        1218 :   S->E1 = E1 = set_init(nbgen);
    1198        1218 :   S->E2 = E2 = set_init(nbgen);
    1199        1218 :   S->T2 = T2 = set_init(nbgen);
    1200        1218 :   S->T31 = T31 = set_init(nbgen);
    1201        1218 :   S->T32 = T32 = set_init(nbgen);
    1202             : 
    1203             :   /* T31 represents the three torsion paths going from left to right */
    1204             :   /* T32 represents the three torsion paths going from right to left */
    1205       60718 :   for (r = 1; r < nbgen; r++)
    1206             :   {
    1207       59500 :     GEN c2 = gel(cusp_list,r+1);
    1208       59500 :     if (c2[3] == type_T)
    1209             :     {
    1210         413 :       GEN c1 = gel(cusp_list,r), path = mkpath(c1,c2), path2 = vecreverse(path);
    1211         413 :       set_insert(T31, path);
    1212         413 :       set_insert(T32, path2);
    1213             :     }
    1214             :   }
    1215             : 
    1216             :   /* to record relations between E2 and E1 */
    1217        1218 :   S->E2_in_terms_of_E1 = zerovec(nbgen);
    1218        1218 :   S->stdE1 = const_vec(nbmanin, NULL);
    1219             : 
    1220             :   /* Assumption later: path [oo,0] is E1[1], path [1,oo] is E2[1] */
    1221             :   {
    1222        1218 :     GEN oo = cusp_infinity();
    1223        1218 :     GEN p1 = mkpath(oo, mkvecsmall2(0,1)); /* [oo, 0] */
    1224        1218 :     GEN p2 = mkpath(mkvecsmall2(1,1), oo); /* [1, oo] */
    1225        1218 :     insert_E(p1, S, p1N);
    1226        1218 :     insert_E(p2, S, p1N);
    1227             :   }
    1228             : 
    1229       60718 :   for (r = 1; r < nbgen; r++)
    1230             :   {
    1231       59500 :     GEN c1 = gel(cusp_list,r);
    1232    14195188 :     for (s = r+1; s <= nbgen; s++)
    1233             :     {
    1234    14135688 :       pari_sp av = avma;
    1235    14135688 :       GEN c2 = gel(cusp_list,s), path;
    1236    14135688 :       GEN d = subii(mulss(c1[1],c2[2]), mulss(c1[2],c2[1]));
    1237    14135688 :       avma = av;
    1238    14135688 :       if (!is_pm1(d)) continue;
    1239             : 
    1240      117782 :       path = mkpath(c1,c2);
    1241      117782 :       if (r+1 == s)
    1242             :       {
    1243       59500 :         GEN w = path;
    1244       59500 :         ulong hash = T31->hash(w); /* T31, T32 use the same hash function */
    1245       59500 :         if (!hash_search2(T31, w, hash) && !hash_search2(T32, w, hash))
    1246             :         {
    1247       59087 :           if (gamma_equiv(path, vecreverse(path), N))
    1248         623 :             set_insert(T2, path);
    1249             :           else
    1250       58464 :             insert_E(path, S, p1N);
    1251             :         }
    1252             :       } else {
    1253       58282 :         set_insert(F, mkvec2(path, mkvecsmall2(r,s)));
    1254       58282 :         set_insert(F, mkvec2(vecreverse(path), mkvecsmall2(s,r)));
    1255             :       }
    1256             :     }
    1257             :   }
    1258        1218 :   setlg(S->E2_in_terms_of_E1, E2->nb+1);
    1259        1218 : }
    1260             : 
    1261             : /* v = \sum n_i g_i, g_i in Sl(2,Z), return \sum n_i g_i^(-1) */
    1262             : static GEN
    1263      742483 : ZSl2_star(GEN v)
    1264             : {
    1265             :   long i, l;
    1266             :   GEN w, G;
    1267      742483 :   if (typ(v) == t_INT) return v;
    1268      742483 :   G = gel(v,1);
    1269      742483 :   w = cgetg_copy(G, &l);
    1270     1761081 :   for (i = 1; i < l; i++)
    1271             :   {
    1272     1018598 :     GEN g = gel(G,i);
    1273     1018598 :     if (typ(g) == t_MAT) g = SL2_inv(g);
    1274     1018598 :     gel(w,i) = g;
    1275             :   }
    1276      742483 :   return ZG_normalize(mkmat2(w, gel(v,2)));
    1277             : }
    1278             : static void
    1279      156968 : ZSl2C_star_inplace(GEN v)
    1280             : {
    1281      156968 :   long i, l = lg(v);
    1282      156968 :   for (i = 1; i < l; i++) gel(v,i) = ZSl2_star(gel(v,i));
    1283      156968 : }
    1284             : 
    1285             : /* Input: h = set of unimodular paths, p1N = P^1(Z/NZ) = Gamma_0(N)\PSL2(Z)
    1286             :  * Output: Each path is converted to a matrix and then an element of P^1(Z/NZ)
    1287             :  * Append the matrix to W[12], append the index that represents
    1288             :  * these elements of P^1 (the classes mod Gamma_0(N) via our fixed
    1289             :  * enumeration to W[2]. */
    1290             : static void
    1291        7308 : paths_decompose(GEN W, hashtable *h, int flag)
    1292             : {
    1293        7308 :   GEN p1N = ms_get_p1N(W), section = ms_get_section(W);
    1294        7308 :   GEN v = hash_to_vec(h);
    1295        7308 :   long i, l = lg(v);
    1296      186221 :   for (i = 1; i < l; i++)
    1297             :   {
    1298      178913 :     GEN e = gel(v,i);
    1299      178913 :     GEN M = path_to_zm(flag? gel(e,1): e);
    1300      178913 :     long index = p1_index(cc(M), dd(M), p1N);
    1301      178913 :     vecsmalltrunc_append(gel(W,2), index);
    1302      178913 :     gel(section, index) = M;
    1303             :   }
    1304        7308 : }
    1305             : static void
    1306        1218 : fill_W2_W12(GEN W, PS_sets_t *S)
    1307             : {
    1308        1218 :   GEN p1N = gel(W,1);
    1309        1218 :   long n = p1_size(p1N);
    1310        1218 :   gel(W, 2) = vecsmalltrunc_init(n+1);
    1311        1218 :   gel(W,12) = cgetg(n+1, t_VEC);
    1312             :   /* F contains [path, [index cusp1, index cusp2]]. Others contain paths only */
    1313        1218 :   paths_decompose(W, S->F, 1);
    1314        1218 :   paths_decompose(W, S->E2, 0);
    1315        1218 :   paths_decompose(W, S->T32, 0);
    1316        1218 :   paths_decompose(W, S->E1, 0);
    1317        1218 :   paths_decompose(W, S->T2, 0);
    1318        1218 :   paths_decompose(W, S->T31, 0);
    1319        1218 : }
    1320             : 
    1321             : /* x t_VECSMALL, corresponds to a map x(i) = j, where 1 <= j <= max for all i
    1322             :  * Return y s.t. y[j] = i or 0 (not in image) */
    1323             : static GEN
    1324        2436 : reverse_list(GEN x, long max)
    1325             : {
    1326        2436 :   GEN y = const_vecsmall(max, 0);
    1327        2436 :   long r, lx = lg(x);
    1328        2436 :   for (r = 1; r < lx; r++) y[ x[r] ] = r;
    1329        2436 :   return y;
    1330             : }
    1331             : 
    1332             : /* go from C[a] to C[b]; return the indices of paths
    1333             :  * E.g. if a < b
    1334             :  *   (C[a]->C[a+1], C[a+1]->C[a+2], ... C[b-1]->C[b])
    1335             :  * (else reverse direction)
    1336             :  * = b - a paths */
    1337             : static GEN
    1338      113876 : F_indices(GEN W, long a, long b)
    1339             : {
    1340      113876 :   GEN v = cgetg(labs(b-a) + 1, t_VEC);
    1341      113876 :   long s, k = 1;
    1342      113876 :   if (a < b) {
    1343       56938 :     GEN index_forward = gel(W,13);
    1344       56938 :     for (s = a; s < b; s++) gel(v,k++) = gel(index_forward,s);
    1345             :   } else {
    1346       56938 :     GEN index_backward = gel(W,14);
    1347       56938 :     for (s = a; s > b; s--) gel(v,k++) = gel(index_backward,s);
    1348             :   }
    1349      113876 :   return v;
    1350             : }
    1351             : /* go from C[a] to C[b] via oo; return the indices of paths
    1352             :  * E.g. if a < b
    1353             :  *   (C[a]->C[a-1], ... C[2]->C[1],
    1354             :  *    C[1]->oo, oo-> C[end],
    1355             :  *    C[end]->C[end-1], ... C[b+1]->C[b])
    1356             :  *  a-1 + 2 + end-(b+1)+1 = end - b + a + 1 paths  */
    1357             : static GEN
    1358        2688 : F_indices_oo(GEN W, long end, long a, long b)
    1359             : {
    1360        2688 :   GEN index_oo = gel(W,15);
    1361        2688 :   GEN v = cgetg(end-labs(b-a)+1 + 1, t_VEC);
    1362        2688 :   long s, k = 1;
    1363             : 
    1364        2688 :   if (a < b) {
    1365        1344 :     GEN index_backward = gel(W,14);
    1366        1344 :     for (s = a; s > 1; s--) gel(v,k++) = gel(index_backward,s);
    1367        1344 :     gel(v,k++) = gel(index_backward,1); /* C[1] -> oo */
    1368        1344 :     gel(v,k++) = gel(index_oo,2); /* oo -> C[end] */
    1369        1344 :     for (s = end; s > b; s--) gel(v,k++) = gel(index_backward,s);
    1370             :   } else {
    1371        1344 :     GEN index_forward = gel(W,13);
    1372        1344 :     for (s = a; s < end; s++) gel(v,k++) = gel(index_forward,s);
    1373        1344 :     gel(v,k++) = gel(index_forward,end); /* C[end] -> oo */
    1374        1344 :     gel(v,k++) = gel(index_oo,1); /* oo -> C[1] */
    1375        1344 :     for (s = 1; s < b; s++) gel(v,k++) = gel(index_forward,s);
    1376             :   }
    1377        2688 :   return v;
    1378             : }
    1379             : /* index of oo -> C[1], oo -> C[end] */
    1380             : static GEN
    1381        1218 : indices_oo(GEN W, GEN C)
    1382             : {
    1383        1218 :   long end = lg(C)-1;
    1384        1218 :   GEN w, v = cgetg(2+1, t_VEC), oo = cusp_infinity();
    1385        1218 :   w = mkpath(oo, gel(C,1)); /* oo -> C[1]=0 */
    1386        1218 :   gel(v,1) = path_Gamma0N_decompose(W, w);
    1387        1218 :   w = mkpath(oo, gel(C,end)); /* oo -> C[end]=1 */
    1388        1218 :   gel(v,2) = path_Gamma0N_decompose(W, w);
    1389        1218 :   return v;
    1390             : }
    1391             : 
    1392             : /* index of C[1]->C[2], C[2]->C[3], ... C[end-1]->C[end], C[end]->oo
    1393             :  * Recall that C[1] = 0, C[end] = 1 */
    1394             : static GEN
    1395        1218 : indices_forward(GEN W, GEN C)
    1396             : {
    1397        1218 :   long s, k = 1, end = lg(C)-1;
    1398        1218 :   GEN v = cgetg(end+1, t_VEC);
    1399       61936 :   for (s = 1; s <= end; s++)
    1400             :   {
    1401       60718 :     GEN w = mkpath(gel(C,s), s == end? cusp_infinity(): gel(C,s+1));
    1402       60718 :     gel(v,k++) = path_Gamma0N_decompose(W, w);
    1403             :   }
    1404        1218 :   return v;
    1405             : }
    1406             : /* index of C[1]->oo, C[2]->C[1], ... C[end]->C[end-1] */
    1407             : static GEN
    1408        1218 : indices_backward(GEN W, GEN C)
    1409             : {
    1410        1218 :   long s, k = 1, end = lg(C)-1;
    1411        1218 :   GEN v = cgetg(end+1, t_VEC);
    1412       61936 :   for (s = 1; s <= end; s++)
    1413             :   {
    1414       60718 :     GEN w = mkpath(gel(C,s), s == 1? cusp_infinity(): gel(C,s-1));
    1415       60718 :     gel(v,k++) = path_Gamma0N_decompose(W, w);
    1416             :   }
    1417        1218 :   return v;
    1418             : }
    1419             : 
    1420             : /* N = integer > 1. Returns data describing Delta_0 = Z[P^1(Q)]_0 seen as
    1421             :  * a Gamma_0(N) - module. */
    1422             : static GEN
    1423        1218 : msinit_N(ulong N)
    1424             : {
    1425        1218 :   GEN p1N = create_p1mod(N);
    1426             :   GEN C, vecF, vecT2, vecT31;
    1427             :   ulong r, s, width;
    1428        1218 :   long nball, nbgen, nbp1N = p1_size(p1N);
    1429        1218 :   GEN TAU = mkmat22(gen_0,gen_m1, gen_1,gen_m1); /*[0,-1;1,-1]*/
    1430             :   GEN W, W2, singlerel, annT2, annT31;
    1431             :   GEN F_index;
    1432             :   hashtable *F, *T2, *T31, *T32, *E1, *E2;
    1433             :   PS_sets_t S;
    1434             : 
    1435        1218 :   form_E_F_T(N,p1N, &C, &S);
    1436        1218 :   E1  = S.E1;
    1437        1218 :   E2  = S.E2;
    1438        1218 :   T31 = S.T31;
    1439        1218 :   T32 = S.T32;
    1440        1218 :   F   = S.F;
    1441        1218 :   T2  = S.T2;
    1442        1218 :   nbgen = lg(C)-1;
    1443             : 
    1444        1218 :   W = cgetg(17, t_VEC);
    1445        1218 :   gel(W,1) = p1N;
    1446             : 
    1447             :  /* Put our paths in the order: F,E2,T32,E1,T2,T31
    1448             :   * W2[j] associates to the j-th element of this list its index in P1. */
    1449        1218 :   fill_W2_W12(W, &S);
    1450        1218 :   W2 = gel(W, 2);
    1451        1218 :   nball = lg(W2)-1;
    1452        1218 :   gel(W,3) = reverse_list(W2, nbp1N);
    1453             : 
    1454        1218 :   gel(W,5) = vecslice(gel(W,2), F->nb + E2->nb + T32->nb + 1, nball);
    1455        1218 :   gel(W,4) = reverse_list(gel(W,5), nbp1N);
    1456        1218 :   gel(W,13) = indices_forward(W, C);
    1457        1218 :   gel(W,14) = indices_backward(W, C);
    1458        1218 :   gel(W,15) = indices_oo(W, C);
    1459        6090 :   gel(W,11) = mkvecsmall5(F->nb,
    1460        1218 :                           F->nb + E2->nb,
    1461        1218 :                           F->nb + E2->nb + T32->nb,
    1462        1218 :                           F->nb + E2->nb + T32->nb + E1->nb,
    1463        1218 :                           F->nb + E2->nb + T32->nb + E1->nb + T2->nb);
    1464             : 
    1465             :   /* relations between T32 and T31 [not stored!]
    1466             :    * T32[i] = - T31[i] */
    1467             : 
    1468             :   /* relations of F */
    1469        1218 :   width = E1->nb + T2->nb + T31->nb;
    1470             :   /* F_index[r] = [index_1, ..., index_k], where index_i is the p1_index()
    1471             :    * of the elementary unimodular path between 2 consecutive cusps
    1472             :    * [in E1,E2,T2,T31 or T32] */
    1473        1218 :   F_index = cgetg(F->nb+1, t_VEC);
    1474        1218 :   vecF = hash_to_vec(F);
    1475      117782 :   for (r = 1; r <= F->nb; r++)
    1476             :   {
    1477      116564 :     GEN w = gel(gel(vecF,r), 2);
    1478      116564 :     long a = w[1], b = w[2], d = labs(b - a);
    1479             :     /* c1 = cusp_list[a],  c2 = cusp_list[b], ci != oo */
    1480      233128 :     gel(F_index,r) = (nbgen-d >= d-1)? F_indices(W, a,b)
    1481      116564 :                                      : F_indices_oo(W, lg(C)-1,a,b);
    1482             :   }
    1483             : 
    1484        1218 :   singlerel = cgetg(width+1, t_VEC);
    1485             :   /* form the single boundary relation */
    1486       31668 :   for (s = 1; s <= E2->nb; s++)
    1487             :   {
    1488       30450 :     GEN data = gel(S.E2_in_terms_of_E1,s);
    1489       30450 :     long c = itos(gel(data,1));
    1490       30450 :     GEN u = gel(data,2); /* E2[s] = u * E1[c], u = - [gamma] */
    1491       30450 :     GEN gamma = gcoeff(u,1,1);
    1492       30450 :     gel(singlerel, c) = mkmat22(gen_1,gen_1, gamma,gen_m1);
    1493             :   }
    1494        1218 :   for (r = E1->nb + 1; r <= width; r++) gel(singlerel, r) = gen_1;
    1495             : 
    1496             :   /* form the 2-torsion relations */
    1497        1218 :   annT2 = cgetg(T2->nb+1, t_VEC);
    1498        1218 :   vecT2 = hash_to_vec(T2);
    1499        1841 :   for (r = 1; r <= T2->nb; r++)
    1500             :   {
    1501         623 :     GEN w = gel(vecT2,r);
    1502         623 :     GEN gamma = gamma_equiv_matrix(vecreverse(w), w);
    1503         623 :     gel(annT2, r) = mkmat22(gen_1,gen_1, gamma,gen_1);
    1504             :   }
    1505             : 
    1506             :   /* form the 3-torsion relations */
    1507        1218 :   annT31 = cgetg(T31->nb+1, t_VEC);
    1508        1218 :   vecT31 = hash_to_vec(T31);
    1509        1631 :   for (r = 1; r <= T31->nb; r++)
    1510             :   {
    1511         413 :     GEN M = path_to_ZM( vecreverse(gel(vecT31,r)) );
    1512         413 :     GEN gamma = ZM_mul(ZM_mul(M, TAU), SL2_inv(M));
    1513         413 :     gel(annT31, r) = mkmat2(mkcol3(gen_1,gamma,ZM_sqr(gamma)),
    1514             :                             mkcol3(gen_1,gen_1,gen_1));
    1515             :   }
    1516        1218 :   gel(W,6) = F_index;
    1517        1218 :   gel(W,7) = S.E2_in_terms_of_E1;
    1518        1218 :   gel(W,8) = annT2;
    1519        1218 :   gel(W,9) = annT31;
    1520        1218 :   gel(W,10)= singlerel;
    1521        1218 :   gel(W,16)= inithashcusps(p1N);
    1522        1218 :   return W;
    1523             : }
    1524             : static GEN
    1525          98 : cusp_to_P1Q(GEN c) { return c[2]? gdivgs(stoi(c[1]), c[2]): mkoo(); }
    1526             : GEN
    1527          14 : mspathgens(GEN W)
    1528             : {
    1529          14 :   pari_sp av = avma;
    1530             :   long i,j, l, nbE1, nbT2, nbT31;
    1531             :   GEN R, r, g, section, gen, annT2, annT31, singlerel;
    1532          14 :   checkms(W); W = get_ms(W);
    1533          14 :   section = ms_get_section(W);
    1534          14 :   gen = ms_get_genindex(W);
    1535          14 :   l = lg(gen);
    1536          14 :   g = cgetg(l,t_VEC);
    1537          63 :   for (i=1; i<l; i++)
    1538             :   {
    1539          49 :     GEN p = gel(section,gen[i]);
    1540          49 :     gel(g,i) = mkvec2(cusp_to_P1Q(gel(p,1)), cusp_to_P1Q(gel(p,2)));
    1541             :   }
    1542          14 :   nbE1 = ms_get_nbE1(W);
    1543          14 :   annT2 = gel(W,8); nbT2 = lg(annT2)-1;
    1544          14 :   annT31 = gel(W,9);nbT31 = lg(annT31)-1;
    1545          14 :   singlerel = gel(W,10);
    1546          14 :   R = cgetg(nbT2+nbT31+2, t_VEC);
    1547          14 :   l = lg(singlerel);
    1548          14 :   r = cgetg(l, t_VEC);
    1549          42 :   for (i = 1; i <= nbE1; i++)
    1550          28 :     gel(r,i) = mkvec2(gel(singlerel, i), stoi(i));
    1551          35 :   for (; i < l; i++)
    1552          21 :     gel(r,i) = mkvec2(gen_1, stoi(i));
    1553          14 :   gel(R,1) = r; j = 2;
    1554          35 :   for (i = 1; i <= nbT2; i++,j++)
    1555          21 :     gel(R,j) = mkvec( mkvec2(gel(annT2,i), stoi(i + nbE1)) );
    1556          14 :   for (i = 1; i <= nbT31; i++,j++)
    1557           0 :     gel(R,j) = mkvec( mkvec2(gel(annT31,i), stoi(i + nbE1 + nbT2)) );
    1558          14 :   return gerepilecopy(av, mkvec2(g,R));
    1559             : }
    1560             : 
    1561             : /* Modular symbols in weight k: Hom_Gamma(Delta, Q[x,y]_{k-2}) */
    1562             : /* A symbol phi is represented by the {phi(g_i)}, {phi(g'_i)}, {phi(g''_i)}
    1563             :  * where the {g_i, g'_i, g''_i} are the Z[\Gamma]-generators of Delta,
    1564             :  * g_i corresponds to E1, g'_i to T2, g''_i to T31.
    1565             :  */
    1566             : 
    1567             : /* FIXME: export. T^1, ..., T^n */
    1568             : static GEN
    1569      510706 : RgX_powers(GEN T, long n)
    1570             : {
    1571      510706 :   GEN v = cgetg(n+1, t_VEC);
    1572             :   long i;
    1573      510706 :   gel(v, 1) = T;
    1574      510706 :   for (i = 1; i < n; i++) gel(v,i+1) = RgX_mul(gel(v,i), T);
    1575      510706 :   return v;
    1576             : }
    1577             : 
    1578             : /* g = [a,b;c,d] a mat2. Return (X^{k-2} | g)(X,Y)[X = 1]. */
    1579             : static GEN
    1580        2576 : voo_act_Gl2Q(GEN g, long k)
    1581             : {
    1582        2576 :   GEN mc = stoi(-coeff(g,2,1)), d = stoi(coeff(g,2,2));
    1583        2576 :   return RgX_to_RgC(gpowgs(deg1pol_shallow(mc, d, 0), k-2), k-1);
    1584             : }
    1585             : 
    1586             : struct m_act {
    1587             :   long dim, k, p;
    1588             :   GEN q;
    1589             : };
    1590             : 
    1591             : /* g = [a,b;c,d]. Return (P | g)(X,Y)[X = 1] = P(dX - cY, -b X + aY)[X = 1],
    1592             :  * for P = X^{k-2}, X_^{k-3}Y, ..., Y^{k-2} */
    1593             : GEN
    1594      255353 : RgX_act_Gl2Q(GEN g, long k)
    1595             : {
    1596             :   GEN a,b,c,d, V1,V2,V;
    1597             :   long i;
    1598      255353 :   if (k == 2) return matid(1);
    1599      255353 :   a = gcoeff(g,1,1); b = gcoeff(g,1,2);
    1600      255353 :   c = gcoeff(g,2,1); d = gcoeff(g,2,2);
    1601      255353 :   V1 = RgX_powers(deg1pol_shallow(gneg(c), d, 0), k-2); /* d - c Y */
    1602      255353 :   V2 = RgX_powers(deg1pol_shallow(a, gneg(b), 0), k-2); /*-b + a Y */
    1603      255353 :   V = cgetg(k, t_MAT);
    1604      255353 :   gel(V,1)   = RgX_to_RgC(gel(V1, k-2), k-1);
    1605      617722 :   for (i = 1; i < k-2; i++)
    1606             :   {
    1607      362369 :     GEN v1 = gel(V1, k-2-i); /* (d-cY)^(k-2-i) */
    1608      362369 :     GEN v2 = gel(V2, i); /* (-b+aY)^i */
    1609      362369 :     gel(V,i+1) = RgX_to_RgC(RgX_mul(v1,v2), k-1);
    1610             :   }
    1611      255353 :   gel(V,k-1) = RgX_to_RgC(gel(V2, k-2), k-1);
    1612      255353 :   return V; /* V[i+1] = X^i | g */
    1613             : }
    1614             : /* z in Z[Gl2(Q)], return the matrix of z acting on V */
    1615             : static GEN
    1616      521409 : act_ZGl2Q(GEN z, struct m_act *T, GEN(*act)(struct m_act*,GEN), hashtable *H)
    1617             : {
    1618      521409 :   GEN S = NULL, G, E;
    1619             :   pari_sp av;
    1620             :   long l, j;
    1621             :   /* paranoia: should'n t occur */
    1622      521409 :   if (typ(z) == t_INT) return scalarmat_shallow(z, T->dim);
    1623      521409 :   G = gel(z,1); l = lg(G);
    1624      521409 :   E = gel(z,2);
    1625      521409 :   if (H)
    1626             :   { /* First pass, identify matrices in Sl_2 to convert to operators;
    1627             :      * insert operators in hashtable. This allows GC in 2nd pass */
    1628     1528905 :     for (j = 1; j < l; j++)
    1629             :     {
    1630     1011129 :       GEN g = gel(G,j);
    1631     1011129 :       if (typ(g) != t_INT)
    1632             :       {
    1633     1011129 :         ulong hash = H->hash(g);
    1634     1011129 :         hashentry *e = hash_search2(H,g,hash);
    1635     1011129 :         if (!e) hash_insert2(H,g,act(T,g),hash);
    1636             :       }
    1637             :     }
    1638             :   }
    1639      521409 :   av = avma;
    1640     1540007 :   for (j = 1; j < l; j++)
    1641             :   {
    1642     1018598 :     GEN M, g = gel(G,j), n = gel(E,j);
    1643     1018598 :     if (typ(g) == t_INT) /* = 1 */
    1644        3591 :       M = n; /* n*Id_dim */
    1645             :     else
    1646             :     {
    1647     1015007 :       if (H)
    1648     1011129 :         M = (GEN)hash_search(H,g)->val; /*search succeeds because of 1st pass*/
    1649             :       else
    1650        3878 :         M = act(T,g);
    1651     1015007 :       if (is_pm1(n))
    1652     1009337 :       { if (signe(n) < 0) M = RgM_neg(M); }
    1653             :       else
    1654        5670 :         M = RgM_Rg_mul(M, n);
    1655             :     }
    1656     1018598 :     if (!S) { S = M; continue; }
    1657      497189 :     S = gadd(S, M);
    1658      497189 :     if (gc_needed(av,1))
    1659             :     {
    1660           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"act_ZGl2Q, j = %ld",j);
    1661           0 :       S = gerepileupto(av, S);
    1662             :     }
    1663             :   }
    1664      521409 :   return gerepilecopy(av, S);
    1665             : }
    1666             : static GEN
    1667      255353 : _RgX_act_Gl2Q(struct m_act *S, GEN z) { return RgX_act_Gl2Q(z, S->k); }
    1668             : /* acting on (X^{k-2},...,Y^{k-2}) */
    1669             : GEN
    1670        3619 : RgX_act_ZGl2Q(GEN z, long k)
    1671             : {
    1672             :   struct m_act T;
    1673        3619 :   T.k = k;
    1674        3619 :   T.dim = k-1;
    1675        3619 :   return act_ZGl2Q(z, &T, _RgX_act_Gl2Q, NULL);
    1676             : }
    1677             : 
    1678             : /* Given a sparse vector of elements in Z[G], convert it to a (sparse) vector
    1679             :  * of operators on V (given by t_MAT) */
    1680             : static void
    1681       39116 : ZGl2QC_to_act(struct m_act *S, GEN(*act)(struct m_act*,GEN), GEN v, hashtable *H)
    1682             : {
    1683       39116 :   GEN val = gel(v,2);
    1684       39116 :   long i, l = lg(val);
    1685       39116 :   for (i = 1; i < l; i++) gel(val,i) = act_ZGl2Q(gel(val,i), S, act, H);
    1686       39116 : }
    1687             : 
    1688             : /* For all V[i] in Z[\Gamma], find the P such that  P . V[i]^* = 0;
    1689             :  * write P in basis X^{k-2}, ..., Y^{k-2} */
    1690             : static GEN
    1691        1106 : ZGV_tors(GEN V, long k)
    1692             : {
    1693        1106 :   long i, l = lg(V);
    1694        1106 :   GEN v = cgetg(l, t_VEC);
    1695        1554 :   for (i = 1; i < l; i++)
    1696             :   {
    1697         448 :     GEN a = ZSl2_star(gel(V,i));
    1698         448 :     gel(v,i) = ZM_ker(RgX_act_ZGl2Q(a,k));
    1699             :   }
    1700        1106 :   return v;
    1701             : }
    1702             : 
    1703             : static long
    1704     6560057 : set_from_index(GEN W11, long i)
    1705             : {
    1706     6560057 :   if (i <= W11[1]) return 1;
    1707     5689432 :   if (i <= W11[2]) return 2;
    1708     3007557 :   if (i <= W11[3]) return 3;
    1709     3002762 :   if (i <= W11[4]) return 4;
    1710       21203 :   if (i <= W11[5]) return 5;
    1711        4550 :   return 6;
    1712             : }
    1713             : 
    1714             : /* det M = 1 */
    1715             : static void
    1716     1401323 : treat_index(GEN W, GEN M, long index, GEN v)
    1717             : {
    1718     1401323 :   GEN W11 = gel(W,11);
    1719     1401323 :   long shift = W11[3]; /* #F + #E2 + T32 */
    1720     1401323 :   switch(set_from_index(W11, index))
    1721             :   {
    1722             :     case 1: /*F*/
    1723             :     {
    1724      230349 :       GEN F_index = gel(W,6), ind = gel(F_index, index);
    1725      230349 :       long j, l = lg(ind);
    1726     1184974 :       for (j = 1; j < l; j++)
    1727             :       {
    1728      954625 :         GEN IND = gel(ind,j), M0 = gel(IND,2);
    1729      954625 :         long index = mael(IND,1,1);
    1730      954625 :         treat_index(W, ZM_mul(M,M0), index, v);
    1731             :       }
    1732      230349 :       break;
    1733             :     }
    1734             : 
    1735             :     case 2: /*E2, E2[r] + gamma * E1[s] = 0 */
    1736             :     {
    1737      532420 :       long r = index - W11[1];
    1738      532420 :       GEN E2_in_terms_of_E1= gel(W,7), z = gel(E2_in_terms_of_E1, r);
    1739      532420 :       long s = itou(gel(z,1));
    1740             : 
    1741      532420 :       index = s;
    1742      532420 :       M = G_ZG_mul(M, gel(z,2)); /* M * (-gamma) */
    1743      532420 :       gel(v, index) = ZG_add(gel(v, index), M);
    1744      532420 :       break;
    1745             :     }
    1746             : 
    1747             :     case 3: /*T32, T32[i] = -T31[i] */
    1748             :     {
    1749        3675 :       long T3shift = W11[5] - W11[2]; /* #T32 + #E1 + #T2 */
    1750        3675 :       index += T3shift;
    1751        3675 :       index -= shift;
    1752        3675 :       gel(v, index) = ZG_add(gel(v, index), to_famat_shallow(M,gen_m1));
    1753        3675 :       break;
    1754             :     }
    1755             :     default: /*E1,T2,T31*/
    1756      634879 :       index -= shift;
    1757      634879 :       gel(v, index) = ZG_add(gel(v, index), to_famat_shallow(M,gen_1));
    1758      634879 :       break;
    1759             :   }
    1760     1401323 : }
    1761             : static void
    1762     5158734 : treat_index_trivial(GEN v, GEN W, long index)
    1763             : {
    1764     5158734 :   GEN W11 = gel(W,11);
    1765     5158734 :   long shift = W11[3]; /* #F + #E2 + T32 */
    1766     5158734 :   switch(set_from_index(W11, index))
    1767             :   {
    1768             :     case 1: /*F*/
    1769             :     {
    1770      640276 :       GEN F_index = gel(W,6), ind = gel(F_index, index);
    1771      640276 :       long j, l = lg(ind);
    1772     4783324 :       for (j = 1; j < l; j++)
    1773             :       {
    1774     4143048 :         GEN IND = gel(ind,j);
    1775     4143048 :         treat_index_trivial(v, W, mael(IND,1,1));
    1776             :       }
    1777      640276 :       break;
    1778             :     }
    1779             : 
    1780             :     case 2: /*E2, E2[r] + gamma * E1[s] = 0 */
    1781             :     {
    1782     2149455 :       long r = index - W11[1];
    1783     2149455 :       GEN E2_in_terms_of_E1 = gel(W,7), z = gel(E2_in_terms_of_E1, r);
    1784     2149455 :       long s = itou(gel(z,1));
    1785     2149455 :       v[s]--;
    1786     2149455 :       break;
    1787             :     }
    1788             : 
    1789             :     case 3: case 5: case 6: /*T32,T2,T31*/
    1790        9989 :       break;
    1791             : 
    1792             :     case 4: /*E1*/
    1793     2359014 :       v[index-shift]++;
    1794     2359014 :       break;
    1795             :   }
    1796     5158734 : }
    1797             : 
    1798             : static GEN
    1799      157052 : M2_log(GEN W, GEN M)
    1800             : {
    1801      157052 :   GEN a = gcoeff(M,1,1), b = gcoeff(M,1,2);
    1802      157052 :   GEN c = gcoeff(M,2,1), d = gcoeff(M,2,2);
    1803             :   GEN  u, v, D, V;
    1804             :   long index, s;
    1805             : 
    1806      157052 :   W = get_ms(W);
    1807      157052 :   V = zerovec(ms_get_nbgen(W));
    1808             : 
    1809      157052 :   D = subii(mulii(a,d), mulii(b,c));
    1810      157052 :   s = signe(D);
    1811      157052 :   if (!s) return V;
    1812      157045 :   if (is_pm1(D))
    1813             :   { /* shortcut, no need to apply Manin's trick */
    1814       55314 :     if (s < 0) {
    1815        3612 :       b = negi(b);
    1816        3612 :       d = negi(d);
    1817             :     }
    1818       55314 :     M = Gamma0N_decompose(W, mkmat22(a,b, c,d), &index);
    1819       55314 :     treat_index(W, M, index, V);
    1820             :   }
    1821             :   else
    1822             :   {
    1823             :     GEN U, B, P, Q, PQ, C1,C2;
    1824             :     long i, l;
    1825      101731 :     (void)bezout(a,c,&u,&v);
    1826      101731 :     B = addii(mulii(b,u), mulii(d,v));
    1827             :     /* [u,v;-c,a] [a,b; c,d] = [1,B; 0,D], i.e. M = U [1,B;0,D] */
    1828      101731 :     U = mkmat22(a,negi(v), c,u);
    1829             : 
    1830             :     /* {1/0 -> B/D} as \sum g_i, g_i unimodular paths */
    1831      101731 :     PQ = ZV_allpnqn( gboundcf(gdiv(B,D), 0) );
    1832      101731 :     P = gel(PQ,1); l = lg(P);
    1833      101731 :     Q = gel(PQ,2);
    1834      101731 :     C1 = gel(U,1);
    1835      493115 :     for (i = 1; i < l; i++, C1 = C2)
    1836             :     {
    1837             :       GEN M;
    1838      391384 :       C2 = ZM_ZC_mul(U, mkcol2(gel(P,i), gel(Q,i)));
    1839      391384 :       if (!odd(i)) C1 = ZC_neg(C1);
    1840      391384 :       M = Gamma0N_decompose(W, mkmat2(C1,C2), &index);
    1841      391384 :       treat_index(W, M, index, V);
    1842             :     }
    1843             :   }
    1844      157045 :   return V;
    1845             : }
    1846             : 
    1847             : /* express +oo->q=a/b in terms of the Z[G]-generators, trivial action */
    1848             : static void
    1849        6552 : Q_log_trivial(GEN v, GEN W, GEN q)
    1850             : {
    1851        6552 :   GEN Q, W3 = gel(W,3), p1N = gel(W,1);
    1852        6552 :   ulong c,d, N = p1N_get_N(p1N);
    1853             :   long i, lx;
    1854             : 
    1855        6552 :   Q = Q_log_init(N, q);
    1856        6552 :   lx = lg(Q);
    1857        6552 :   c = 0;
    1858       28350 :   for (i = 1; i < lx; i++, c = d)
    1859             :   {
    1860             :     long index;
    1861       21798 :     d = Q[i];
    1862       21798 :     if (c && !odd(i)) c = N - c;
    1863       21798 :     index = W3[ p1_index(c,d,p1N) ];
    1864       21798 :     treat_index_trivial(v, W, index);
    1865             :   }
    1866        6552 : }
    1867             : static void
    1868      436282 : M2_log_trivial(GEN V, GEN W, GEN M)
    1869             : {
    1870      436282 :   GEN p1N = gel(W,1), W3 = gel(W,3);
    1871      436282 :   ulong N = p1N_get_N(p1N);
    1872      436282 :   GEN a = gcoeff(M,1,1), b = gcoeff(M,1,2);
    1873      436282 :   GEN c = gcoeff(M,2,1), d = gcoeff(M,2,2);
    1874             :   GEN  u, v, D;
    1875             :   long index, s;
    1876             : 
    1877      436282 :   D = subii(mulii(a,d), mulii(b,c));
    1878      436282 :   s = signe(D);
    1879      441189 :   if (!s) return;
    1880      436282 :   if (is_pm1(D))
    1881             :   { /* shortcut, not need to apply Manin's trick */
    1882      180243 :     if (s < 0) d = negi(d);
    1883      180243 :     index = W3[ p1_index(umodiu(c,N),umodiu(d,N),p1N) ];
    1884      180243 :     treat_index_trivial(V, W, index);
    1885             :   }
    1886             :   else
    1887             :   {
    1888             :     GEN U, B, P, Q, PQ;
    1889             :     long i, l;
    1890      256039 :     if (!signe(c)) { Q_log_trivial(V,W,gdiv(b,d)); return; }
    1891      251132 :     (void)bezout(a,c,&u,&v);
    1892      251132 :     B = addii(mulii(b,u), mulii(d,v));
    1893             :     /* [u,v;-c,a] [a,b; c,d] = [1,B; 0,D], i.e. M = U [1,B;0,D] */
    1894      251132 :     U = mkvec2(c, u);
    1895             : 
    1896             :     /* {1/0 -> B/D} as \sum g_i, g_i unimodular paths */
    1897      251132 :     PQ = ZV_allpnqn( gboundcf(gdiv(B,D), 0) );
    1898      251132 :     P = gel(PQ,1); l = lg(P);
    1899      251132 :     Q = gel(PQ,2);
    1900     1064777 :     for (i = 1; i < l; i++, c = d)
    1901             :     {
    1902      813645 :       d = addii(mulii(gel(U,1),gel(P,i)), mulii(gel(U,2),gel(Q,i)));
    1903      813645 :       if (!odd(i)) c = negi(c);
    1904      813645 :       index = W3[ p1_index(umodiu(c,N),umodiu(d,N),p1N) ];
    1905      813645 :       treat_index_trivial(V, W, index);
    1906             :     }
    1907             :   }
    1908             : }
    1909             : 
    1910             : static GEN
    1911         238 : cusp_to_ZC(GEN c)
    1912             : {
    1913         238 :   switch(typ(c))
    1914             :   {
    1915             :     case t_INFINITY:
    1916          28 :       return mkcol2(gen_1,gen_0);
    1917             :     case t_INT:
    1918          84 :       return mkcol2(c,gen_1);
    1919             :     case t_FRAC:
    1920         126 :       return mkcol2(gel(c,1),gel(c,2));
    1921             :     case t_VECSMALL:
    1922           0 :       return mkcol2(stoi(c[1]), stoi(c[2]));
    1923             :     default:
    1924           0 :       pari_err_TYPE("mspathlog",c);
    1925           0 :       return NULL;
    1926             :   }
    1927             : }
    1928             : static GEN
    1929         119 : path2_to_M2(GEN p)
    1930         119 : { return mkmat2(cusp_to_ZC(gel(p,1)), cusp_to_ZC(gel(p,2))); }
    1931             : static GEN
    1932         133 : path_to_M2(GEN p)
    1933             : {
    1934         133 :   if (lg(p) != 3) pari_err_TYPE("mspathlog",p);
    1935         126 :   switch(typ(p))
    1936             :   {
    1937             :     case t_MAT:
    1938           7 :       RgM_check_ZM(p,"mspathlog");
    1939           7 :       break;
    1940             :     case t_VEC:
    1941         119 :       p = path2_to_M2(p);
    1942         119 :       break;
    1943           0 :     default: pari_err_TYPE("mspathlog",p);
    1944             :   }
    1945         126 :   return p;
    1946             : }
    1947             : /* Expresses path p as \sum x_i g_i, where the g_i are our distinguished
    1948             :  * generators and x_i \in Z[\Gamma]. Returns [x_1,...,x_n] */
    1949             : GEN
    1950          98 : mspathlog(GEN W, GEN p)
    1951             : {
    1952          98 :   pari_sp av = avma;
    1953          98 :   checkms(W);
    1954          98 :   return gerepilecopy(av, M2_log(W, path_to_M2(p)));
    1955             : }
    1956             : 
    1957             : /** HECKE OPERATORS **/
    1958             : /* [a,b;c,d] * cusp */
    1959             : static GEN
    1960     1181376 : cusp_mul(long a, long b, long c, long d, GEN cusp)
    1961             : {
    1962     1181376 :   long x = cusp[1], y = cusp[2];
    1963     1181376 :   long A = a*x+b*y, B = c*x+d*y, u = cgcd(A,B);
    1964     1181376 :   if (u != 1) { A /= u; B /= u; }
    1965     1181376 :   return mkcol2s(A, B);
    1966             : }
    1967             : /* f in Gl2(Q), act on path (zm), return path_to_M2(f.path) */
    1968             : static GEN
    1969      590688 : Gl2Q_act_path(GEN f, GEN path)
    1970             : {
    1971      590688 :   long a = coeff(f,1,1), b = coeff(f,1,2);
    1972      590688 :   long c = coeff(f,2,1), d = coeff(f,2,2);
    1973      590688 :   GEN c1 = cusp_mul(a,b,c,d, gel(path,1));
    1974      590688 :   GEN c2 = cusp_mul(a,b,c,d, gel(path,2));
    1975      590688 :   return mkmat2(c1,c2);
    1976             : }
    1977             : 
    1978             : static GEN
    1979      127477 : init_act_trivial(GEN W) { return const_vecsmall(ms_get_nbE1(W), 0); }
    1980             : static GEN
    1981          35 : mspathlog_trivial(GEN W, GEN p)
    1982             : {
    1983             :   GEN v;
    1984          35 :   W = get_ms(W);
    1985          35 :   v = init_act_trivial(W);
    1986          35 :   M2_log_trivial(v, W, path_to_M2(p));
    1987          28 :   return v;
    1988             : }
    1989             : 
    1990             : /* map from W1=Hom(Delta_0(N1),Q) -> W2=Hom(Delta_0(N2),Q), weight 2,
    1991             :  * trivial action. v a Gl2_Q or a t_VEC of Gl2_Q (\sum v[i] in Z[Gl2(Q)]).
    1992             :  * Return the matrix attached to the action of v. */
    1993             : static GEN
    1994        2478 : getMorphism_trivial(GEN WW1, GEN WW2, GEN v)
    1995             : {
    1996        2478 :   GEN W1 = get_ms(WW1), W2 = get_ms(WW2);
    1997        2478 :   GEN section = ms_get_section(W2), gen = ms_get_genindex(W2);
    1998        2478 :   long j, lv, d2 = ms_get_nbE1(W2);
    1999        2478 :   GEN T = cgetg(d2+1, t_MAT);
    2000        2478 :   lv = lg(v);
    2001      128275 :   for (j = 1; j <= d2; j++)
    2002             :   {
    2003      125797 :     GEN w = gel(section, gen[j]);
    2004      125797 :     GEN t = init_act_trivial(W1);
    2005             :     long l;
    2006      125797 :     for (l = 1; l < lv; l++) M2_log_trivial(t, W1, Gl2Q_act_path(gel(v,l), w));
    2007      125797 :     gel(T,j) = t;
    2008             :   }
    2009        2478 :   return shallowtrans(zm_to_ZM(T));
    2010             : }
    2011             : 
    2012             : static GEN
    2013      156968 : RgV_sparse(GEN v, GEN *pind)
    2014             : {
    2015             :   long i, l, k;
    2016      156968 :   GEN w = cgetg_copy(v,&l), ind = cgetg(l, t_VECSMALL);
    2017    16914646 :   for (i = k = 1; i < l; i++)
    2018             :   {
    2019    16757678 :     GEN c = gel(v,i);
    2020    16757678 :     if (typ(c) == t_INT) continue;
    2021      738864 :     gel(w,k) = c; ind[k] = i; k++;
    2022             :   }
    2023      156968 :   setlg(w,k); setlg(ind,k);
    2024      156968 :   *pind = ind; return w;
    2025             : }
    2026             : 
    2027             : static hashtable *
    2028        3479 : Gl2act_cache(long dim) { return set_init(dim*10); }
    2029             : 
    2030             : /* f zm/ZM in Gl_2(Q), acts from the left on Delta, which is generated by
    2031             :  * (g_i) as Z[Gamma1]-module, and by (G_i) as Z[Gamma2]-module.
    2032             :  * We have f.G_j = \sum_i \lambda_{i,j} g_i,   \lambda_{i,j} in Z[Gamma1]
    2033             :  * For phi in Hom_Gamma1(D,V), g in D, phi | f is in Hom_Gamma2(D,V) and
    2034             :  *  (phi | f)(G_j) = phi(f.G_j) | f
    2035             :  *                 = phi( \sum_i \lambda_{i,j} g_i ) | f
    2036             :  *                 = \sum_i phi(g_i) | (\lambda_{i,j}^* f)
    2037             :  *                 = \sum_i phi(g_i) | \mu_{i,j}(f)
    2038             :  * More generally
    2039             :  *  (\sum_k (phi |v_k))(G_j) = \sum_i phi(g_i) | \Mu_{i,j}
    2040             :  * with \Mu_{i,j} = \sum_k \mu{i,j}(v_k)
    2041             :  * Return the \Mu_{i,j} matrix as vector of sparse columns of operators on V */
    2042             : static GEN
    2043        3031 : init_dual_act(GEN v, GEN W1, GEN W2, struct m_act *S,
    2044             :               GEN(*act)(struct m_act *,GEN))
    2045             : {
    2046        3031 :   GEN section = ms_get_section(W2), gen = ms_get_genindex(W2);
    2047             :   /* HACK: the actions we consider in dimension 1 are trivial and in
    2048             :    * characteristic != 2, 3 => torsion generators are 0
    2049             :    * [satisfy e.g. (1+gamma).g = 0 => \phi(g) | 1+gamma  = 0 => \phi(g) = 0 */
    2050        3031 :   long j, lv = lg(v), dim = S->dim == 1? ms_get_nbE1(W2): lg(gen)-1;
    2051        3031 :   GEN T = cgetg(dim+1, t_VEC);
    2052        3031 :   hashtable *H = Gl2act_cache(dim);
    2053             : 
    2054       39613 :   for (j = 1; j <= dim; j++)
    2055             :   {
    2056       36582 :     pari_sp av = avma;
    2057       36582 :     GEN w = gel(section, gen[j]); /* path_to_zm( E1/T2/T3 element ) */
    2058       36582 :     GEN t = NULL;
    2059             :     long k;
    2060      191016 :     for (k = 1; k < lv; k++)
    2061             :     {
    2062      154434 :       GEN ind, L, F, tk, f = gel(v,k);
    2063      154434 :       if (typ(gel(f,1)) == t_VECSMALL) F = mat2_to_ZM(f);
    2064           0 :       else { F = f; f = ZM_to_zm(F); }
    2065             :       /* f zm = F ZM */
    2066      154434 :       L = M2_log(W1, Gl2Q_act_path(f,w)); /* L[i] = lambda_{i,j} */
    2067      154434 :       L = RgV_sparse(L,&ind);
    2068      154434 :       ZSl2C_star_inplace(L); /* L[i] = lambda_{i,j}^* */
    2069      154434 :       if (!ZM_isidentity(F)) ZGC_G_mul_inplace(L, F);
    2070      154434 :       tk = mkvec2(ind,L); /* L[i] = mu_{i,j}(v[k]) */
    2071      154434 :       t = t? ZGCs_add(t, tk): tk;
    2072             :     }
    2073       36582 :     gel(T,j) = gerepilecopy(av, t);
    2074             :   }
    2075        3031 :   for(j = 1; j <= dim; j++) ZGl2QC_to_act(S, act, gel(T,j), H);
    2076        3031 :   return T;
    2077             : }
    2078             : 
    2079             : /* modular symbol given by phi[j] = \phi(G_j)
    2080             :  * \sum L[i]*phi[i], L a sparse column of operators */
    2081             : static GEN
    2082      348124 : dense_act_col(GEN col, GEN phi)
    2083             : {
    2084      348124 :   GEN s = NULL, colind = gel(col,1), colval = gel(col,2);
    2085      348124 :   long i, l = lg(colind), lphi = lg(phi);
    2086     5472061 :   for (i = 1; i < l; i++)
    2087             :   {
    2088     5126023 :     long a = colind[i];
    2089             :     GEN t;
    2090     5126023 :     if (a >= lphi) break; /* happens if k=2: torsion generator t omitted */
    2091     5123937 :     t = gel(phi, a); /* phi(G_a) */
    2092     5123937 :     t = RgM_RgC_mul(gel(colval,i), t);
    2093     5123937 :     s = s? RgC_add(s, t): t;
    2094             :   }
    2095      348124 :   return s;
    2096             : }
    2097             : /* modular symbol given by \phi( G[ind[j]] ) = val[j]
    2098             :  * \sum L[i]*phi[i], L a sparse column of operators */
    2099             : static GEN
    2100      771176 : sparse_act_col(GEN col, GEN phi)
    2101             : {
    2102      771176 :   GEN s = NULL, colind = gel(col,1), colval = gel(col,2);
    2103      771176 :   GEN ind = gel(phi,2), val = gel(phi,3);
    2104      771176 :   long a, l = lg(ind);
    2105     3003833 :   for (a = 1; a < l; a++)
    2106             :   {
    2107     2232657 :     GEN t = gel(val, a); /* phi(G_i) */
    2108     2232657 :     long i = zv_search(colind, ind[a]);
    2109     2232657 :     if (!i) continue;
    2110      531804 :     t = RgM_RgC_mul(gel(colval,i), t);
    2111      531804 :     s = s? RgC_add(s, t): t;
    2112             :   }
    2113      771176 :   return s;
    2114             : }
    2115             : static int
    2116       67900 : phi_sparse(GEN phi) { return typ(gel(phi,1)) == t_VECSMALL; }
    2117             : /* phi in Hom_Gamma1(Delta, V), return the matrix whose colums are the
    2118             :  *   \sum_i phi(g_i) | \mu_{i,j} = (phi|f)(G_j),
    2119             :  * see init_dual_act. */
    2120             : static GEN
    2121       67900 : dual_act(long dimV, GEN act, GEN phi)
    2122             : {
    2123       67900 :   long l = lg(act), j;
    2124       67900 :   GEN v = cgetg(l, t_MAT);
    2125       67900 :   GEN (*ACT)(GEN,GEN) = phi_sparse(phi)? sparse_act_col: dense_act_col;
    2126     1184288 :   for (j = 1; j < l; j++)
    2127             :   {
    2128     1116388 :     pari_sp av = avma;
    2129     1116388 :     GEN s = ACT(gel(act,j), phi);
    2130     1116388 :     gel(v,j) = s? gerepileupto(av,s): zerocol(dimV);
    2131             :   }
    2132       67900 :   return v;
    2133             : }
    2134             : 
    2135             : /* \phi in Hom(Delta, V), \phi(G_k) = phi[k]. Write \phi as
    2136             :  *   \sum_{i,j} mu_{i,j} phi_{i,j}, mu_{i,j} in Q */
    2137             : static GEN
    2138       58373 : getMorphism_basis(GEN W, GEN phi)
    2139             : {
    2140       58373 :   GEN basis = msk_get_basis(W);
    2141       58373 :   long i, j, r, lvecT = lg(phi), dim = lg(basis)-1;
    2142       58373 :   GEN st = msk_get_st(W);
    2143       58373 :   GEN link = msk_get_link(W);
    2144       58373 :   GEN invphiblock = msk_get_invphiblock(W);
    2145       58373 :   long s = st[1], t = st[2];
    2146       58373 :   GEN R = zerocol(dim), Q, Ls, T0, T1, Ts;
    2147      781788 :   for (r = 2; r < lvecT; r++)
    2148             :   {
    2149             :     GEN Tr, L;
    2150      723415 :     if (r == s) continue;
    2151      665042 :     Tr = gel(phi,r); /* Phi(G_r), r != 1,s */
    2152      665042 :     L = gel(link, r);
    2153      665042 :     Q = ZC_apply_dinv(gel(invphiblock,r), Tr);
    2154             :     /* write Phi(G_r) as sum_{a,b} mu_{a,b} Phi_{a,b}(G_r) */
    2155      665042 :     for (j = 1; j < lg(L); j++) gel(R, L[j]) = gel(Q,j);
    2156             :   }
    2157       58373 :   Ls = gel(link, s);
    2158       58373 :   T1 = gel(phi,1); /* Phi(G_1) */
    2159       58373 :   gel(R, Ls[t]) = gel(T1, 1);
    2160             : 
    2161       58373 :   T0 = NULL;
    2162      781788 :   for (i = 2; i < lg(link); i++)
    2163             :   {
    2164             :     GEN L;
    2165      723415 :     if (i == s) continue;
    2166      665042 :     L = gel(link,i);
    2167     3485286 :     for (j =1 ; j < lg(L); j++)
    2168             :     {
    2169     2820244 :       long n = L[j]; /* phi_{i,j} = basis[n] */
    2170     2820244 :       GEN mu_ij = gel(R, n);
    2171     2820244 :       GEN phi_ij = gel(basis, n), pols = gel(phi_ij,3);
    2172     2820244 :       GEN z = RgC_Rg_mul(gel(pols, 3), mu_ij);
    2173     2820244 :       T0 = T0? RgC_add(T0, z): z; /* += mu_{i,j} Phi_{i,j} (G_s) */
    2174             :     }
    2175             :   }
    2176       58373 :   Ts = gel(phi,s); /* Phi(G_s) */
    2177       58373 :   if (T0) Ts = RgC_sub(Ts, T0);
    2178             :   /* solve \sum_{j!=t} mu_{s,j} Phi_{s,j}(G_s) = Ts */
    2179       58373 :   Q = ZC_apply_dinv(gel(invphiblock,s), Ts);
    2180       58373 :   for (j = 1; j < t; j++) gel(R, Ls[j]) = gel(Q,j);
    2181             :   /* avoid mu_{s,t} */
    2182       58373 :   for (j = t; j < lg(Q); j++) gel(R, Ls[j+1]) = gel(Q,j);
    2183       58373 :   return R;
    2184             : }
    2185             : 
    2186             : /* a = s(g_i) for some modular symbol s; b in Z[G]
    2187             :  * return s(b.g_i) = b^* . s(g_i) */
    2188             : static GEN
    2189         147 : ZGl2Q_act_s(GEN b, GEN a, long k)
    2190             : {
    2191         147 :   if (typ(b) == t_INT)
    2192             :   {
    2193          56 :     if (!signe(b)) return gen_0;
    2194          14 :     switch(typ(a))
    2195             :     {
    2196             :       case t_POL:
    2197          14 :         a = RgX_to_RgC(a, k-1); /*fall through*/
    2198             :       case t_COL:
    2199          14 :         a = RgC_Rg_mul(a,b);
    2200          14 :         break;
    2201           0 :       default: a = scalarcol_shallow(b,k-1);
    2202             :     }
    2203             :   }
    2204             :   else
    2205             :   {
    2206          91 :     b = RgX_act_ZGl2Q(ZSl2_star(b), k);
    2207          91 :     switch(typ(a))
    2208             :     {
    2209             :       case t_POL:
    2210          63 :         a = RgX_to_RgC(a, k-1); /*fall through*/
    2211             :       case t_COL:
    2212          91 :         a = RgM_RgC_mul(b,a);
    2213          91 :         break;
    2214           0 :       default: a = RgC_Rg_mul(gel(b,1),a);
    2215             :     }
    2216             :   }
    2217         105 :   return a;
    2218             : }
    2219             : 
    2220             : static int
    2221          21 : checksymbol(GEN W, GEN s)
    2222             : {
    2223             :   GEN t, annT2, annT31, singlerel;
    2224             :   long i, k, l, nbE1, nbT2, nbT31;
    2225          21 :   k = msk_get_weight(W);
    2226          21 :   W = get_ms(W);
    2227          21 :   nbE1 = ms_get_nbE1(W);
    2228          21 :   singlerel = gel(W,10);
    2229          21 :   l = lg(singlerel);
    2230          21 :   if (k == 2)
    2231             :   {
    2232           0 :     for (i = nbE1+1; i < l; i++)
    2233           0 :       if (!gequal0(gel(s,i))) return 0;
    2234           0 :     return 1;
    2235             :   }
    2236          21 :   annT2 = gel(W,8); nbT2 = lg(annT2)-1;
    2237          21 :   annT31 = gel(W,9);nbT31 = lg(annT31)-1;
    2238          21 :   t = NULL;
    2239          84 :   for (i = 1; i < l; i++)
    2240             :   {
    2241          63 :     GEN a = gel(s,i);
    2242          63 :     a = ZGl2Q_act_s(gel(singlerel,i), a, k);
    2243          63 :     t = t? gadd(t, a): a;
    2244             :   }
    2245          21 :   if (!gequal0(t)) return 0;
    2246          14 :   for (i = 1; i <= nbT2; i++)
    2247             :   {
    2248           0 :     GEN a = gel(s,i + nbE1);
    2249           0 :     a = ZGl2Q_act_s(gel(annT2,i), a, k);
    2250           0 :     if (!gequal0(a)) return 0;
    2251             :   }
    2252          28 :   for (i = 1; i <= nbT31; i++)
    2253             :   {
    2254          14 :     GEN a = gel(s,i + nbE1 + nbT2);
    2255          14 :     a = ZGl2Q_act_s(gel(annT31,i), a, k);
    2256          14 :     if (!gequal0(a)) return 0;
    2257             :   }
    2258          14 :   return 1;
    2259             : }
    2260             : long
    2261          28 : msissymbol(GEN W, GEN s)
    2262             : {
    2263             :   long k, nbgen;
    2264          28 :   checkms(W);
    2265          28 :   k = msk_get_weight(W);
    2266          28 :   nbgen = ms_get_nbgen(W);
    2267          28 :   switch(typ(s))
    2268             :   {
    2269             :     case t_VEC: /* values s(g_i) */
    2270          21 :       if (lg(s)-1 != nbgen) return 0;
    2271          21 :       break;
    2272             :     case t_COL:
    2273           7 :       if (msk_get_sign(W))
    2274             :       {
    2275           0 :         GEN star = gel(msk_get_starproj(W), 1);
    2276           0 :         if (lg(star) == lg(s)) return 1;
    2277             :       }
    2278           7 :       if (k == 2) /* on the dual basis of (g_i) */
    2279             :       {
    2280           0 :         if (lg(s)-1 != nbgen) return 0;
    2281             :       }
    2282             :       else
    2283             :       {
    2284           7 :         GEN basis = msk_get_basis(W);
    2285           7 :         return (lg(s) == lg(basis));
    2286             :       }
    2287           0 :       break;
    2288           0 :     default: return 0;
    2289             :   }
    2290          21 :   return checksymbol(W,s);
    2291             : }
    2292             : #if DEBUG
    2293             : /* phi is a sparse symbol from msk_get_basis, return phi(G_j) */
    2294             : static GEN
    2295             : phi_Gj(GEN W, GEN phi, long j)
    2296             : {
    2297             :   GEN ind = gel(phi,2), pols = gel(phi,3);
    2298             :   long i = vecsmall_isin(ind,j);
    2299             :   return i? gel(pols,i): NULL;
    2300             : }
    2301             : /* check that \sum d_i phi_i(G_j)  = T_j for all j */
    2302             : static void
    2303             : checkdec(GEN W, GEN D, GEN T)
    2304             : {
    2305             :   GEN B = msk_get_basis(W);
    2306             :   long i, j;
    2307             :   if (!checksymbol(W,T)) pari_err_BUG("checkdec");
    2308             :   for (j = 1; j < lg(T); j++)
    2309             :   {
    2310             :     GEN S = gen_0;
    2311             :     for (i = 1; i < lg(D); i++)
    2312             :     {
    2313             :       GEN d = gel(D,i), v = phi_Gj(W, gel(B,i), j);
    2314             :       if (!v || gequal0(d)) continue;
    2315             :       S = gadd(S, gmul(d, v));
    2316             :     }
    2317             :     /* S = \sum_i d_i phi_i(G_j) */
    2318             :     if (!gequal(S, gel(T,j)))
    2319             :       pari_warn(warner, "checkdec j = %ld\n\tS = %Ps\n\tT = %Ps", j,S,gel(T,j));
    2320             :   }
    2321             : }
    2322             : #endif
    2323             : 
    2324             : /* map op: W1 = Hom(Delta_0(N1),V) -> W2 = Hom(Delta_0(N2),V), given by
    2325             :  * \sum v[i], v[i] in Gl2(Q) */
    2326             : static GEN
    2327        5054 : getMorphism(GEN W1, GEN W2, GEN v)
    2328             : {
    2329             :   struct m_act S;
    2330             :   GEN B1, M, act;
    2331        5054 :   long a, l, k = msk_get_weight(W1);
    2332        5054 :   if (k == 2) return getMorphism_trivial(W1,W2,v);
    2333        2576 :   S.k = k;
    2334        2576 :   S.dim = k-1;
    2335        2576 :   act = init_dual_act(v,W1,W2,&S, _RgX_act_Gl2Q);
    2336        2576 :   B1 = msk_get_basis(W1);
    2337        2576 :   l = lg(B1); M = cgetg(l, t_MAT);
    2338       60018 :   for (a = 1; a < l; a++)
    2339             :   {
    2340       57442 :     pari_sp av = avma;
    2341       57442 :     GEN phi = dual_act(S.dim, act, gel(B1,a));
    2342       57442 :     GEN D = getMorphism_basis(W2, phi);
    2343             : #if DEBUG
    2344             :     checkdec(W2,D,T);
    2345             : #endif
    2346       57442 :     gel(M,a) = gerepilecopy(av, D);
    2347             :   }
    2348        2576 :   return M;
    2349             : }
    2350             : static GEN
    2351        3934 : msendo(GEN W, GEN v) { return getMorphism(W, W, v); }
    2352             : 
    2353             : static GEN
    2354        2422 : endo_project(GEN W, GEN e, GEN H)
    2355             : {
    2356        2422 :   if (msk_get_sign(W)) e = Qevproj_apply(e, msk_get_starproj(W));
    2357        2422 :   if (H) e = Qevproj_apply(e, Qevproj_init0(H));
    2358        2422 :   return e;
    2359             : }
    2360             : static GEN
    2361        2688 : mshecke_i(GEN W, ulong p)
    2362             : {
    2363        2688 :   GEN v = ms_get_N(W) % p? Tp_matrices(p): Up_matrices(p);
    2364        2688 :   return msendo(W,v);
    2365             : }
    2366             : GEN
    2367        2387 : mshecke(GEN W, long p, GEN H)
    2368             : {
    2369        2387 :   pari_sp av = avma;
    2370             :   GEN T;
    2371        2387 :   checkms(W);
    2372        2387 :   if (p <= 1) pari_err_PRIME("mshecke",stoi(p));
    2373        2387 :   T = mshecke_i(W,p);
    2374        2387 :   T = endo_project(W,T,H);
    2375        2387 :   return gerepilecopy(av, T);
    2376             : }
    2377             : 
    2378             : static GEN
    2379          35 : msatkinlehner_i(GEN W, long Q)
    2380             : {
    2381          35 :   long N = ms_get_N(W);
    2382             :   GEN v;
    2383          35 :   if (Q == 1) return matid(msk_get_dim(W));
    2384          28 :   if (Q == N) return msendo(W, mkvec(mat2(0,1,-N,0)));
    2385          21 :   if (N % Q) pari_err_DOMAIN("msatkinlehner","N % Q","!=",gen_0,stoi(Q));
    2386          14 :   v = WQ_matrix(N, Q);
    2387          14 :   if (!v) pari_err_DOMAIN("msatkinlehner","gcd(Q,N/Q)","!=",gen_1,stoi(Q));
    2388          14 :   return msendo(W,mkvec(v));
    2389             : }
    2390             : GEN
    2391          35 : msatkinlehner(GEN W, long Q, GEN H)
    2392             : {
    2393          35 :   pari_sp av = avma;
    2394             :   GEN w;
    2395             :   long k;
    2396          35 :   checkms(W);
    2397          35 :   k = msk_get_weight(W);
    2398          35 :   if (Q <= 0) pari_err_DOMAIN("msatkinlehner","Q","<=",gen_0,stoi(Q));
    2399          35 :   w = msatkinlehner_i(W,Q);
    2400          28 :   w = endo_project(W,w,H);
    2401          28 :   if (k > 2 && Q != 1) w = RgM_Rg_div(w, powuu(Q,(k-2)>>1));
    2402          28 :   return gerepilecopy(av, w);
    2403             : }
    2404             : 
    2405             : static GEN
    2406        1225 : msstar_i(GEN W) { return msendo(W, mkvec(mat2(-1,0,0,1))); }
    2407             : GEN
    2408           7 : msstar(GEN W, GEN H)
    2409             : {
    2410           7 :   pari_sp av = avma;
    2411             :   GEN s;
    2412           7 :   checkms(W);
    2413           7 :   s = msstar_i(W);
    2414           7 :   s = endo_project(W,s,H);
    2415           7 :   return gerepilecopy(av, s);
    2416             : }
    2417             : 
    2418             : #if 0
    2419             : /* is \Gamma_0(N) cusp1 = \Gamma_0(N) cusp2 ? */
    2420             : static int
    2421             : iscuspeq(ulong N, GEN cusp1, GEN cusp2)
    2422             : {
    2423             :   long p1, q1, p2, q2, s1, s2, d;
    2424             :   p1 = cusp1[1]; p2 = cusp2[1];
    2425             :   q1 = cusp1[2]; q2 = cusp2[2];
    2426             :   d = Fl_mul(smodss(q1,N),smodss(q2,N), N);
    2427             :   d = ugcd(d, N);
    2428             : 
    2429             :   s1 = q1 > 2? Fl_inv(smodss(p1,q1), q1): 1;
    2430             :   s2 = q2 > 2? Fl_inv(smodss(p2,q2), q2): 1;
    2431             :   return Fl_mul(s1,q2,d) == Fl_mul(s2,q1,d);
    2432             : }
    2433             : #endif
    2434             : 
    2435             : /* return E_c(r) */
    2436             : static GEN
    2437        2576 : get_Ec_r(GEN c, long k)
    2438             : {
    2439        2576 :   long p = c[1], q = c[2], u, v;
    2440             :   GEN gr;
    2441        2576 :   (void)cbezout(p, q, &u, &v);
    2442        2576 :   gr = mat2(p, -v, q, u); /* g . (1:0) = (p:q) */
    2443        2576 :   return voo_act_Gl2Q(sl2_inv(gr), k);
    2444             : }
    2445             : /* returns the modular symbol attached to the cusp c := p/q via the rule
    2446             :  * E_c(path from a to b in Delta_0) := E_c(b) - E_c(a), where
    2447             :  * E_c(r) := 0 if r != c mod Gamma
    2448             :  *           v_oo | gamma_r^(-1)
    2449             :  * where v_oo is stable by T = [1,1;0,1] (i.e x^(k-2)) and
    2450             :  * gamma_r . (1:0) = r, for some gamma_r in SL_2(Z) * */
    2451             : static GEN
    2452         434 : msfromcusp_trivial(GEN W, GEN c)
    2453             : {
    2454         434 :   GEN section = ms_get_section(W), gen = ms_get_genindex(W);
    2455         434 :   GEN S = ms_get_hashcusps(W);
    2456         434 :   long j, ic = cusp_index(c, S), l = ms_get_nbE1(W)+1;
    2457         434 :   GEN phi = cgetg(l, t_COL);
    2458       90034 :   for (j = 1; j < l; j++)
    2459             :   {
    2460       89600 :     GEN vj, g = gel(section, gen[j]); /* path_to_zm(generator) */
    2461       89600 :     GEN c1 = gel(g,1), c2 = gel(g,2);
    2462       89600 :     long i1 = cusp_index(c1, S);
    2463       89600 :     long i2 = cusp_index(c2, S);
    2464       89600 :     if (i1 == ic)
    2465        3199 :       vj = (i2 == ic)?  gen_0: gen_1;
    2466             :     else
    2467       86401 :       vj = (i2 == ic)? gen_m1: gen_0;
    2468       89600 :     gel(phi, j) = vj;
    2469             :   }
    2470         434 :   return phi;
    2471             : }
    2472             : static GEN
    2473        1365 : msfromcusp_i(GEN W, GEN c)
    2474             : {
    2475             :   GEN section, gen, S, phi;
    2476        1365 :   long j, ic, l, k = msk_get_weight(W);
    2477        1365 :   if (k == 2) return msfromcusp_trivial(W, c);
    2478         931 :   k = msk_get_weight(W);
    2479         931 :   section = ms_get_section(W);
    2480         931 :   gen = ms_get_genindex(W);
    2481         931 :   S = ms_get_hashcusps(W);
    2482         931 :   ic = cusp_index(c, S);
    2483         931 :   l = lg(gen);
    2484         931 :   phi = cgetg(l, t_COL);
    2485       11543 :   for (j = 1; j < l; j++)
    2486             :   {
    2487       10612 :     GEN vj = NULL, g = gel(section, gen[j]); /* path_to_zm(generator) */
    2488       10612 :     GEN c1 = gel(g,1), c2 = gel(g,2);
    2489       10612 :     long i1 = cusp_index(c1, S);
    2490       10612 :     long i2 = cusp_index(c2, S);
    2491       10612 :     if (i1 == ic) vj = get_Ec_r(c1, k);
    2492       10612 :     if (i2 == ic)
    2493             :     {
    2494        1288 :       GEN s = get_Ec_r(c2, k);
    2495        1288 :       vj = vj? gsub(vj, s): gneg(s);
    2496             :     }
    2497       10612 :     if (!vj) vj = zerocol(k-1);
    2498       10612 :     gel(phi, j) = vj;
    2499             :   }
    2500         931 :   return getMorphism_basis(W, phi);
    2501             : }
    2502             : GEN
    2503          21 : msfromcusp(GEN W, GEN c)
    2504             : {
    2505          21 :   pari_sp av = avma;
    2506             :   long N;
    2507          21 :   checkms(W);
    2508          21 :   N = ms_get_N(W);
    2509          21 :   switch(typ(c))
    2510             :   {
    2511             :     case t_INFINITY:
    2512           7 :       c = mkvecsmall2(1,0);
    2513           7 :       break;
    2514             :     case t_INT:
    2515           7 :       c = mkvecsmall2(smodis(c,N), 1);
    2516           7 :       break;
    2517             :     case t_FRAC:
    2518           7 :       c = mkvecsmall2(smodis(gel(c,1),N), smodis(gel(c,2),N));
    2519           7 :       break;
    2520             :     default:
    2521           0 :       pari_err_TYPE("msfromcusp",c);
    2522             :   }
    2523          21 :   return gerepilecopy(av, msfromcusp_i(W,c));
    2524             : }
    2525             : 
    2526             : static GEN
    2527         294 : mseisenstein_i(GEN W)
    2528             : {
    2529         294 :   GEN M, S = ms_get_hashcusps(W), cusps = gel(S,3);
    2530         294 :   long i, l = lg(cusps);
    2531         294 :   if (msk_get_weight(W)==2) l--;
    2532         294 :   M = cgetg(l, t_MAT);
    2533         294 :   for (i = 1; i < l; i++) gel(M,i) = msfromcusp_i(W, gel(cusps,i));
    2534         294 :   return Qevproj_star(W, QM_image(M));
    2535             : }
    2536             : GEN
    2537         294 : mseisenstein(GEN W)
    2538             : {
    2539         294 :   pari_sp av = avma;
    2540         294 :   checkms(W);
    2541         294 :   return gerepilecopy(av, Qevproj_init(mseisenstein_i(W)));
    2542             : }
    2543             : 
    2544             : /* upper bound for log_2 |charpoly(T_p|S)|, where S is a cuspidal subspace of
    2545             :  * dimension d, k is the weight */
    2546             : #if 0
    2547             : static long
    2548             : TpS_char_bound(ulong p, long k, long d)
    2549             : { /* |eigenvalue| <= 2 p^(k-1)/2 */
    2550             :   return d * (2 + (log2((double)p)*(k-1))/2);
    2551             : }
    2552             : #endif
    2553             : static long
    2554         287 : TpE_char_bound(ulong p, long k, long d)
    2555             : { /* |eigenvalue| <= 2 p^(k-1) */
    2556         287 :   return d * (2 + log2((double)p)*(k-1));
    2557             : }
    2558             : 
    2559             : GEN
    2560         287 : mscuspidal(GEN W, long flag)
    2561             : {
    2562         287 :   pari_sp av = avma;
    2563             :   GEN S, E, M, T, TE, chE;
    2564             :   long bit;
    2565             :   forprime_t F;
    2566             :   ulong p, N;
    2567             :   pari_timer ti;
    2568             : 
    2569         287 :   E = mseisenstein(W);
    2570         287 :   N = ms_get_N(W);
    2571         287 :   (void)u_forprime_init(&F, 2, ULONG_MAX);
    2572         287 :   while ((p = u_forprime_next(&F)))
    2573         399 :     if (N % p) break;
    2574         287 :   if (DEBUGLEVEL) timer_start(&ti);
    2575         287 :   T = mshecke(W, p, NULL);
    2576         287 :   if (DEBUGLEVEL) timer_printf(&ti,"Tp, p = %ld", p);
    2577         287 :   TE = Qevproj_apply(T, E); /* T_p | E */
    2578         287 :   if (DEBUGLEVEL) timer_printf(&ti,"Qevproj_init(E)");
    2579         287 :   bit = TpE_char_bound(p, msk_get_weight(W), lg(TE)-1);
    2580         287 :   chE = QM_charpoly_ZX_bound(TE, bit);
    2581         287 :   chE = ZX_radical(chE);
    2582         287 :   M = RgX_RgM_eval(chE, T);
    2583         287 :   S = Qevproj_init(QM_image(M));
    2584         287 :   return gerepilecopy(av, flag? mkvec2(S,E): S);
    2585             : }
    2586             : 
    2587             : /** INIT ELLSYM STRUCTURE **/
    2588             : /* V a vector of ZM. If all of them have 0 last row, return NULL.
    2589             :  * Otherwise return [m,i,j], where m = V[i][last,j] contains the value
    2590             :  * of smallest absolute value */
    2591             : static GEN
    2592         812 : RgMV_find_non_zero_last_row(long offset, GEN V)
    2593             : {
    2594         812 :   long i, lasti = 0, lastj = 0, lV = lg(V);
    2595         812 :   GEN m = NULL;
    2596        3668 :   for (i = 1; i < lV; i++)
    2597             :   {
    2598        2856 :     GEN M = gel(V,i);
    2599        2856 :     long j, n, l = lg(M);
    2600        2856 :     if (l == 1) continue;
    2601        2597 :     n = nbrows(M);
    2602       12908 :     for (j = 1; j < l; j++)
    2603             :     {
    2604       10311 :       GEN a = gcoeff(M, n, j);
    2605       10311 :       if (!gequal0(a) && (!m || abscmpii(a, m) < 0))
    2606             :       {
    2607        1414 :         m = a; lasti = i; lastj = j;
    2608        1414 :         if (is_pm1(m)) goto END;
    2609             :       }
    2610             :     }
    2611             :   }
    2612             : END:
    2613         812 :   if (!m) return NULL;
    2614         553 :   return mkvec2(m, mkvecsmall2(lasti+offset, lastj));
    2615             : }
    2616             : /* invert the d_oo := (\gamma_oo - 1) operator, acting on
    2617             :  * [x^(k-2), ..., y^(k-2)] */
    2618             : static GEN
    2619         553 : Delta_inv(GEN doo, long k)
    2620             : {
    2621         553 :   GEN M = RgX_act_ZGl2Q(doo, k);
    2622         553 :   M = RgM_minor(M, k-1, 1); /* 1st column and last row are 0 */
    2623         553 :   return ZM_inv_denom(M);
    2624             : }
    2625             : /* The ZX P = \sum a_i x^i y^{k-2-i} is given by the ZV [a_0, ..., a_k-2]~,
    2626             :  * return Q and d such that P = doo Q + d y^k-2, where d in Z and Q */
    2627             : static GEN
    2628       12089 : doo_decompose(GEN dinv, GEN P, GEN *pd)
    2629             : {
    2630       12089 :   long l = lg(P); *pd = gel(P, l-1);
    2631       12089 :   P = vecslice(P, 1, l-2);
    2632       12089 :   return shallowconcat(gen_0, ZC_apply_dinv(dinv, P));
    2633             : }
    2634             : 
    2635             : static GEN
    2636       12089 : get_phi_ij(long i,long j,long n, long s,long t,GEN P_st,GEN Q_st,GEN d_st,
    2637             :            GEN P_ij, GEN lP_ij, GEN dinv)
    2638             : {
    2639             :   GEN ind, pols;
    2640       12089 :   if (i == s && j == t)
    2641             :   {
    2642         553 :     ind = mkvecsmall(1);
    2643         553 :     pols = mkvec(scalarcol_shallow(gen_1, lg(P_st)-1)); /* x^{k-2} */
    2644             :   }
    2645             :   else
    2646             :   {
    2647       11536 :     GEN d_ij, Q_ij = doo_decompose(dinv, lP_ij, &d_ij);
    2648       11536 :     GEN a = ZC_Z_mul(P_ij, d_st);
    2649       11536 :     GEN b = ZC_Z_mul(P_st, negi(d_ij));
    2650       11536 :     GEN c = RgC_sub(RgC_Rg_mul(Q_ij, d_st), RgC_Rg_mul(Q_st, d_ij));
    2651       11536 :     if (i == s) { /* j != t */
    2652        1526 :       ind = mkvecsmall2(1, s);
    2653        1526 :       pols = mkvec2(c, ZC_add(a, b));
    2654             :     } else {
    2655       10010 :       ind = mkvecsmall3(1, i, s);
    2656       10010 :       pols = mkvec3(c, a, b); /* image of g_1, g_i, g_s */
    2657             :     }
    2658       11536 :     pols = Q_primpart(pols);
    2659             :   }
    2660       12089 :   return mkvec3(mkvecsmall3(i,j,n), ind, pols);
    2661             : }
    2662             : 
    2663             : static GEN
    2664         665 : mskinit_trivial(GEN WN)
    2665             : {
    2666         665 :   long dim = ms_get_nbE1(WN);
    2667         665 :   return mkvec3(WN, gen_0, mkvec2(gen_0,mkvecsmall2(2, dim)));
    2668             : }
    2669             : /* sum of #cols of the matrices contained in V */
    2670             : static long
    2671        1106 : RgMV_dim(GEN V)
    2672             : {
    2673        1106 :   long l = lg(V), d = 0, i;
    2674        1106 :   for (i = 1; i < l; i++) d += lg(gel(V,i)) - 1;
    2675        1106 :   return d;
    2676             : }
    2677             : static GEN
    2678         553 : mskinit_nontrivial(GEN WN, long k)
    2679             : {
    2680         553 :   GEN annT2 = gel(WN,8), annT31 = gel(WN,9), singlerel = gel(WN,10);
    2681             :   GEN link, basis, monomials, invphiblock;
    2682         553 :   long nbE1 = ms_get_nbE1(WN);
    2683         553 :   GEN dinv = Delta_inv(ZG_neg( ZSl2_star(gel(singlerel,1)) ), k);
    2684         553 :   GEN p1 = cgetg(nbE1+1, t_VEC), remove;
    2685         553 :   GEN p2 = ZGV_tors(annT2, k);
    2686         553 :   GEN p3 = ZGV_tors(annT31, k);
    2687         553 :   GEN gentor = shallowconcat(p2, p3);
    2688             :   GEN P_st, lP_st, Q_st, d_st;
    2689             :   long n, i, dim, s, t, u;
    2690         553 :   gel(p1, 1) = cgetg(1,t_MAT); /* dummy */
    2691        3080 :   for (i = 2; i <= nbE1; i++) /* skip 1st element = (\gamma_oo-1)g_oo */
    2692             :   {
    2693        2527 :     GEN z = gel(singlerel, i);
    2694        2527 :     gel(p1, i) = RgX_act_ZGl2Q(ZSl2_star(z), k);
    2695             :   }
    2696         553 :   remove = RgMV_find_non_zero_last_row(nbE1, gentor);
    2697         553 :   if (!remove) remove = RgMV_find_non_zero_last_row(0, p1);
    2698         553 :   if (!remove) pari_err_BUG("msinit [no y^k-2]");
    2699         553 :   remove = gel(remove,2); /* [s,t] */
    2700         553 :   s = remove[1];
    2701         553 :   t = remove[2];
    2702             :   /* +1 because of = x^(k-2), but -1 because of Manin relation */
    2703         553 :   dim = (k-1)*(nbE1-1) + RgMV_dim(p2) + RgMV_dim(p3);
    2704             :   /* Let (g_1,...,g_d) be the Gamma-generators of Delta, g_1 = g_oo.
    2705             :    * We describe modular symbols by the collection phi(g_1), ..., phi(g_d)
    2706             :    * \in V := Q[x,y]_{k-2}, with right Gamma action.
    2707             :    * For each i = 1, .., d, let V_i \subset V be the Q-vector space of
    2708             :    * allowed values for phi(g_i): with basis (P^{i,j}) given by the monomials
    2709             :    * x^(j-1) y^{k-2-(j-1)}, j = 1 .. k-1
    2710             :    * (g_i in E_1) or the solution of the torsion equations (1 + gamma)P = 0
    2711             :    * (g_i in T2) or (1 + gamma + gamma^2)P = 0 (g_i in T31). All such P
    2712             :    * are chosen in Z[x,y] with Q_content 1.
    2713             :    *
    2714             :    * The Manin relation (singlerel) is of the form \sum_i \lambda_i g_i = 0,
    2715             :    * where \lambda_i = 1 if g_i in T2 or T31, and \lambda_i = (1 - \gamma_i)
    2716             :    * for g_i in E1.
    2717             :    *
    2718             :    * If phi \in Hom_Gamma(Delta, V), it is defined by phi(g_i) := P_i in V
    2719             :    * with \sum_i P_i . \lambda_i^* = 0, where (\sum n_i g_i)^* :=
    2720             :    * \sum n_i \gamma_i^(-1).
    2721             :    *
    2722             :    * We single out gamma_1 / g_1 (g_oo in Pollack-Stevens paper) and
    2723             :    * write P_{i,j} \lambda_i^* =  Q_{i,j} (\gamma_1 - 1)^* + d_{i,j} y^{k-2}
    2724             :    * where d_{i,j} is a scalar and Q_{i,j} in V; we normalize Q_{i,j} to
    2725             :    * that the coefficient of x^{k-2} is 0.
    2726             :    *
    2727             :    * There exist (s,t) such that d_{s,t} != 0.
    2728             :    * A Q-basis of the (dual) space of modular symbols is given by the
    2729             :    * functions phi_{i,j}, 2 <= i <= d, 1 <= j <= k-1, mapping
    2730             :    *  g_1 -> d_{s,t} Q_{i,j} - d_{i,j} Q_{s,t} + [(i,j)=(s,t)] x^{k-2}
    2731             :    * If i != s
    2732             :    *   g_i -> d_{s,t} P_{i,j}
    2733             :    *   g_s -> - d_{i,j} P_{s,t}
    2734             :    * If i = s, j != t
    2735             :    *   g_i -> d_{s,t} P_{i,j} - d_{i,j} P_{s,t}
    2736             :    * And everything else to 0. Again we normalize the phi_{i,j} such that
    2737             :    * their image has content 1. */
    2738         553 :   monomials = matid(k-1); /* represent the monomials x^{k-2}, ... , y^{k-2} */
    2739         553 :   if (s <= nbE1) /* in E1 */
    2740             :   {
    2741         259 :     P_st = gel(monomials, t);
    2742         259 :     lP_st = gmael(p1, s, t); /* P_{s,t} lambda_s^* */
    2743             :   }
    2744             :   else /* in T2, T31 */
    2745             :   {
    2746         294 :     P_st = gmael(gentor, s - nbE1, t);
    2747         294 :     lP_st = P_st;
    2748             :   }
    2749         553 :   Q_st = doo_decompose(dinv, lP_st, &d_st);
    2750         553 :   basis = cgetg(dim+1, t_VEC);
    2751         553 :   link = cgetg(nbE1 + lg(gentor), t_VEC);
    2752         553 :   gel(link,1) = cgetg(1,t_VECSMALL); /* dummy */
    2753         553 :   n = 1;
    2754        3080 :   for (i = 2; i <= nbE1; i++)
    2755             :   {
    2756        2527 :     GEN L = cgetg(k, t_VECSMALL);
    2757             :     long j;
    2758             :     /* link[i][j] = n gives correspondance between phi_{i,j} and basis[n] */
    2759        2527 :     gel(link,i) = L;
    2760       13160 :     for (j = 1; j < k; j++)
    2761             :     {
    2762       10633 :       GEN lP_ij = gmael(p1, i, j); /* P_{i,j} lambda_i^* */
    2763       10633 :       GEN P_ij = gel(monomials,j);
    2764       10633 :       L[j] = n;
    2765       10633 :       gel(basis, n) = get_phi_ij(i,j,n, s,t, P_st, Q_st, d_st, P_ij, lP_ij, dinv);
    2766       10633 :       n++;
    2767             :     }
    2768             :   }
    2769        1001 :   for (u = 1; u < lg(gentor); u++,i++)
    2770             :   {
    2771         448 :     GEN V = gel(gentor,u);
    2772         448 :     long j, lV = lg(V);
    2773         448 :     GEN L = cgetg(lV, t_VECSMALL);
    2774         448 :     gel(link,i) = L;
    2775        1904 :     for (j = 1; j < lV; j++)
    2776             :     {
    2777        1456 :       GEN lP_ij = gel(V, j); /* P_{i,j} lambda_i^* = P_{i,j} */
    2778        1456 :       GEN P_ij = lP_ij;
    2779        1456 :       L[j] = n;
    2780        1456 :       gel(basis, n) = get_phi_ij(i,j,n, s,t, P_st, Q_st, d_st, P_ij, lP_ij, dinv);
    2781        1456 :       n++;
    2782             :     }
    2783             :   }
    2784         553 :   invphiblock = cgetg(lg(link), t_VEC);
    2785         553 :   gel(invphiblock,1) = cgetg(1, t_MAT); /* dummy */
    2786        3528 :   for (i = 2; i < lg(link); i++)
    2787             :   {
    2788        2975 :     GEN M, inv, B = gel(link,i);
    2789        2975 :     long j, lB = lg(B);
    2790        2975 :     if (i == s) { B = vecsplice(B, t); lB--; } /* remove phi_st */
    2791        2975 :     M = cgetg(lB, t_MAT);
    2792       14511 :     for (j = 1; j < lB; j++)
    2793             :     {
    2794       11536 :       GEN phi_ij = gel(basis, B[j]), pols = gel(phi_ij,3);
    2795       11536 :       gel(M, j) = gel(pols, 2); /* phi_ij(g_i) */
    2796             :     }
    2797        2975 :     if (i <= nbE1 && i != s) /* maximal rank k-1 */
    2798        2268 :       inv = ZM_inv_denom(M);
    2799             :     else /* i = s (rank k-2) or from torsion: rank k/3 or k/2 */
    2800         707 :       inv = Qevproj_init(M);
    2801        2975 :     gel(invphiblock,i) = inv;
    2802             :   }
    2803         553 :   return mkvec3(WN, gen_0, mkvec5(basis, mkvecsmall2(k, dim), mkvecsmall2(s,t),
    2804             :                                   link, invphiblock));
    2805             : }
    2806             : static GEN
    2807        1218 : add_star(GEN W, long sign)
    2808             : {
    2809        1218 :   GEN s = msstar_i(W);
    2810        1218 :   GEN K = sign? QM_ker_r(gsubgs(s, sign)): cgetg(1,t_MAT);
    2811        1218 :   gel(W,2) = mkvec3(stoi(sign), s, Qevproj_init(K));
    2812        1218 :   return W;
    2813             : }
    2814             : /* WN = msinit_N(N) */
    2815             : static GEN
    2816        1218 : mskinit(ulong N, long k, long sign)
    2817             : {
    2818        1218 :   GEN WN = msinit_N(N);
    2819        1218 :   GEN W = k == 2? mskinit_trivial(WN)
    2820        1218 :                 : mskinit_nontrivial(WN, k);
    2821        1218 :   return add_star(W, sign);
    2822             : }
    2823             : GEN
    2824         364 : msinit(GEN N, GEN K, long sign)
    2825             : {
    2826         364 :   pari_sp av = avma;
    2827             :   GEN W;
    2828             :   long k;
    2829         364 :   if (typ(N) != t_INT) pari_err_TYPE("msinit", N);
    2830         364 :   if (typ(K) != t_INT) pari_err_TYPE("msinit", K);
    2831         364 :   k = itos(K);
    2832         364 :   if (k < 2) pari_err_DOMAIN("msinit","k", "<", gen_2,K);
    2833         364 :   if (odd(k)) pari_err_IMPL("msinit [odd weight]");
    2834         364 :   if (signe(N) <= 0) pari_err_DOMAIN("msinit","N", "<=", gen_0,N);
    2835         364 :   if (equali1(N)) pari_err_IMPL("msinit [ N = 1 ]");
    2836         364 :   W = mskinit(itou(N), k, sign);
    2837         364 :   return gerepilecopy(av, W);
    2838             : }
    2839             : 
    2840             : /* W = msinit, xpm integral modular symbol of weight 2, c t_FRAC
    2841             :  * Return image of <oo->c> */
    2842             : static GEN
    2843        1645 : Q_xpm(GEN W, GEN xpm, GEN c)
    2844             : {
    2845        1645 :   pari_sp av = avma;
    2846             :   GEN v;
    2847        1645 :   W = get_ms(W);
    2848        1645 :   v = init_act_trivial(W);
    2849        1645 :   Q_log_trivial(v, W, c); /* oo -> (a:b), c = a/b */
    2850        1645 :   return gerepileuptoint(av, ZV_zc_mul(xpm, v));
    2851             : }
    2852             : 
    2853             : static GEN
    2854          35 : eval_single(GEN s, long k, GEN B, long v)
    2855             : {
    2856             :   long i, l;
    2857          35 :   GEN A = cgetg_copy(s,&l);
    2858          35 :   for (i=1; i<l; i++) gel(A,i) = ZGl2Q_act_s(gel(B,i), gel(s,i), k);
    2859          35 :   A = RgV_sum(A);
    2860          35 :   if (is_vec_t(typ(A))) A = RgV_to_RgX(A, v);
    2861          35 :   return A;
    2862             : }
    2863             : /* Evaluate symbol s on mspathlog B (= sum p_i g_i, p_i in Z[G]). Allow
    2864             :  * s = t_MAT [ collection of symbols, return a vector ]*/
    2865             : static GEN
    2866          63 : mseval_by_values(GEN W, GEN s, GEN p, long v)
    2867             : {
    2868          63 :   long i, l, k = msk_get_weight(W);
    2869             :   GEN A;
    2870          63 :   if (k == 2)
    2871             :   { /* trivial represention: don't bother with Z[G] */
    2872          35 :     GEN B = mspathlog_trivial(W,p);
    2873          28 :     if (typ(s) != t_MAT) return RgV_zc_mul(s,B);
    2874           0 :     l = lg(s); A = cgetg(l, t_VEC);
    2875           0 :     for (i = 1; i < l; i++) gel(A,i) = RgV_zc_mul(gel(s,i), B);
    2876             :   }
    2877             :   else
    2878             :   {
    2879          28 :     GEN B = mspathlog(W,p);
    2880          28 :     if (typ(s) != t_MAT) return eval_single(s, k, B, v);
    2881           7 :     l = lg(s); A = cgetg(l, t_VEC);
    2882           7 :     for (i = 1; i < l; i++) gel(A,i) = eval_single(gel(s,i), k, B, v);
    2883             :   }
    2884           7 :   return A;
    2885             : }
    2886             : 
    2887             : /* express symbol on the basis phi_{i,j} */
    2888             : static GEN
    2889         385 : symtophi(GEN W, GEN s)
    2890             : {
    2891         385 :   GEN e, basis = msk_get_basis(W);
    2892         385 :   long i, l = lg(basis);
    2893         385 :   if (lg(s) != l) pari_err_TYPE("mseval",s);
    2894         385 :   e = const_vec(ms_get_nbgen(W), gen_0);
    2895       12565 :   for (i=1; i<l; i++)
    2896             :   {
    2897       12180 :     GEN phi, ind, pols, c = gel(s,i);
    2898             :     long j, m;
    2899       12180 :     if (gequal0(c)) continue;
    2900       12019 :     phi = gel(basis,i);
    2901       12019 :     ind = gel(phi,2); m = lg(ind);
    2902       12019 :     pols = gel(phi,3);
    2903       46249 :     for (j=1; j<m; j++)
    2904             :     {
    2905       34230 :       long t = ind[j];
    2906       34230 :       gel(e,t) = gadd(gel(e,t), gmul(c, gel(pols,j)));
    2907             :     }
    2908             :   }
    2909         385 :   return e;
    2910             : }
    2911             : /* evaluate symbol s on path p */
    2912             : GEN
    2913         952 : mseval(GEN W, GEN s, GEN p)
    2914             : {
    2915         952 :   pari_sp av = avma;
    2916         952 :   long i, k, l, v = 0;
    2917         952 :   checkms(W);
    2918         952 :   k = msk_get_weight(W);
    2919         952 :   switch(typ(s))
    2920             :   {
    2921             :     case t_VEC: /* values s(g_i) */
    2922           7 :       if (lg(s)-1 != ms_get_nbgen(W)) pari_err_TYPE("mseval",s);
    2923           7 :       if (!p) return gcopy(s);
    2924           0 :       v = gvar(s);
    2925           0 :       break;
    2926             :     case t_COL:
    2927         931 :       if (msk_get_sign(W))
    2928             :       {
    2929         336 :         GEN star = gel(msk_get_starproj(W), 1);
    2930         336 :         if (lg(star) == lg(s)) s = RgM_RgC_mul(star, s);
    2931             :       }
    2932         931 :       if (k == 2) /* on the dual basis of (g_i) */
    2933             :       {
    2934         560 :         if (lg(s)-1 != ms_get_nbE1(W)) pari_err_TYPE("mseval",s);
    2935         560 :         if (!p) return gtrans(s);
    2936             :       }
    2937             :       else
    2938         371 :         s = symtophi(W,s);
    2939         406 :       break;
    2940             :     case t_MAT:
    2941          14 :       if (!p) pari_err_TYPE("mseval",s);
    2942          14 :       l = lg(s);
    2943          14 :       if (l == 1) return cgetg(1, t_VEC);
    2944           7 :       if (msk_get_sign(W))
    2945             :       {
    2946           0 :         GEN star = gel(msk_get_starproj(W), 1);
    2947           0 :         if (lg(star) == lgcols(s)) s = RgM_mul(star, s);
    2948             :       }
    2949           7 :       if (k == 2)
    2950           0 :       { if (nbrows(s) != ms_get_nbE1(W)) pari_err_TYPE("mseval",s); }
    2951             :       else
    2952             :       {
    2953           7 :         GEN t = cgetg(l, t_MAT);
    2954           7 :         for (i = 1; i < l; i++) gel(t,i) = symtophi(W,gel(s,i));
    2955           7 :         s = t;
    2956             :       }
    2957           7 :       break;
    2958           0 :     default: pari_err_TYPE("mseval",s);
    2959             :   }
    2960         413 :   if (p)
    2961          63 :     s = mseval_by_values(W, s, p, v);
    2962             :   else
    2963             :   {
    2964         350 :     l = lg(s);
    2965        3577 :     for (i = 1; i < l; i++)
    2966             :     {
    2967        3227 :       GEN c = gel(s,i);
    2968        3227 :       if (!isintzero(c)) gel(s,i) = RgV_to_RgX(gel(s,i), v);
    2969             :     }
    2970             :   }
    2971         406 :   return gerepilecopy(av, s);
    2972             : }
    2973             : 
    2974             : /* sum_{a <= |D|} (D/a)*xpm(E,a/|D|) */
    2975             : static GEN
    2976         539 : get_X(GEN W, GEN xpm, long D)
    2977             : {
    2978         539 :   ulong a, d = (ulong)labs(D);
    2979         539 :   GEN t = gen_0;
    2980             :   GEN nc, c;
    2981         539 :   if (d == 1) return Q_xpm(W, xpm, gen_0);
    2982         238 :   nc = icopy(gen_1);
    2983         238 :   c = mkfrac(nc, utoipos(d));
    2984        2072 :   for (a=1; a < d; a++)
    2985             :   {
    2986        1834 :     long s = kross(D,a);
    2987             :     GEN x;
    2988        1834 :     if (!s) continue;
    2989        1344 :     nc[2] = a; x = Q_xpm(W, xpm, c);
    2990        1344 :     t = (s > 0)? addii(t, x): subii(t, x);
    2991             :   }
    2992         238 :   return t;
    2993             : }
    2994             : static long
    2995         385 : torsion_order(GEN E) { GEN T = elltors(E); return itos(gel(T,1)); }
    2996             : /* E of rank 0, minimal model; write L(E,1) = Q*w1(E) != 0 and return the
    2997             :  * rational Q; tam = product of all Tamagawa (incl. c_oo(E)). */
    2998             : static GEN
    2999         385 : get_Q(GEN E, GEN tam)
    3000             : {
    3001         385 :   GEN L, sha, w1 = gel(ellR_omega(E,DEFAULTPREC), 1);
    3002         385 :   long ex, t = torsion_order(E), t2 = t*t;
    3003             : 
    3004         385 :   L = ellL1(E, 0, DEFAULTPREC);
    3005         385 :   sha = divrr(mulru(L, t2), mulri(w1,tam)); /* integral = |Sha| by BSD */
    3006         385 :   sha = sqri( grndtoi(sqrtr(sha), &ex) ); /* |Sha| is a square */
    3007         385 :   if (ex > -5) pari_err_BUG("msfromell (can't compute analytic |Sha|)");
    3008         385 :   return gdivgs(mulii(tam,sha), t2);
    3009             : }
    3010             : 
    3011             : /* E given by a minimal model; D != 0. Compare Euler factor of L(E,(D/.),1)
    3012             :  * with L(E^D,1). Return
    3013             :  *   \prod_{p|D} (p-a_p(E)+eps_{E}(p)) / p,
    3014             :  * where eps(p) = 0 if p | N_E and 1 otherwise */
    3015             : static GEN
    3016         161 : get_Euler(GEN E, GEN D)
    3017             : {
    3018         161 :   GEN a = gen_1, b = gen_1, P = gel(absZ_factor(D), 1);
    3019         161 :   long i, l = lg(P);
    3020         336 :   for (i = 1; i < l; i++)
    3021             :   {
    3022         175 :     GEN p = gel(P,i);
    3023         175 :     a = mulii(a, ellcard(E, p));
    3024         175 :     b = mulii(b, p);
    3025             :   }
    3026         161 :   return gdiv(a, b);
    3027             : }
    3028             : 
    3029             : /* E given by a minimal model, xpm in the sign(D) part with the same
    3030             :  * eigenvalues as E (unique up to multiplication with a rational).
    3031             :  * Let X(D) = \sum_{a <= |D|} (D/a) * xpm(E, a/|D|)
    3032             :  * Return the rational correction factor A such that
    3033             :  *   A * X(D) = L(E, (D/.), 1) / \Omega(E^D)
    3034             :  * for fundamental D (such that E^D has rank 0 otherwise both sides vanish). */
    3035             : static GEN
    3036         539 : ell_get_scale_d(GEN E, GEN W, GEN xpm, long D)
    3037             : {
    3038         539 :   GEN gD, cb, N, Q, tam, u, Ed, X = get_X(W, xpm, D);
    3039             : 
    3040         539 :   if (!signe(X)) return NULL;
    3041         385 :   if (D == 1)
    3042             :   {
    3043         224 :     gD = NULL;
    3044         224 :     Ed = E;
    3045             :   }
    3046             :   else
    3047             :   {
    3048         161 :     gD = stoi(D);
    3049         161 :     Ed = ellinit(elltwist(E, gD), NULL, DEFAULTPREC);
    3050             :   }
    3051         385 :   Ed = ellanal_globalred_all(Ed, &cb, &N, &tam);
    3052         385 :   Q =  get_Q(Ed, tam);
    3053         385 :   if (cb)
    3054             :   { /* \tilde{u} in Pal's "Periods of quadratic twists of elliptic curves" */
    3055         182 :     u = gel(cb,1); /* Omega(E^D_min) = u * Omega(E^D) */
    3056         182 :     if (abscmpiu(Q_denom(u), 2) > 0) pari_err_BUG("msfromell [ell_get_scale]");
    3057         182 :     Q = gmul(Q,u);
    3058             :   }
    3059             :   /* L(E^D,1) = Q * w1(E^D_min) */
    3060         385 :   if (gD) Q = gmul(Q, get_Euler(Ed, gD));
    3061         385 :   if (D != 1) obj_free(Ed);
    3062             :   /* L(E^D,1) / Omega(E^D) = Q. Divide by X to get A */
    3063         385 :   return gdiv(Q, X);
    3064             : }
    3065             : 
    3066             : /* Let W = msinit(conductor(E), 2), xpm an integral modular symbol with the same
    3067             :  * eigenvalues as L_E. There exist a unique C such that
    3068             :  *   C*L(E,(D/.),1)_{xpm} = L(E,(D/.),1) / w1(E_D) != 0, for all D fundamental,
    3069             :  * sign(D) = s, and such that E_D has rank 0. Return the normalized symbol
    3070             :  * C * xpm */
    3071             : static GEN
    3072         385 : ell_get_scale(GEN E, GEN W, GEN xpm, long s)
    3073             : {
    3074             :   long d;
    3075             :   /* find D = s*d such that twist by D has rank 0 */
    3076        1106 :   for (d = 1; d < LONG_MAX; d++)
    3077             :   {
    3078        1106 :     pari_sp av = avma;
    3079             :     GEN C;
    3080        1106 :     long D = s > 0? d: -d;
    3081        1106 :     if (!sisfundamental(D)) continue;
    3082         539 :     C = ell_get_scale_d(E, W, xpm, D);
    3083         539 :     if (C) return RgC_Rg_mul(xpm, C);
    3084         154 :     avma = av;
    3085             :   }
    3086           0 :   pari_err_BUG("msfromell (no suitable twist)");
    3087           0 :   return NULL;
    3088             : }
    3089             : 
    3090             : /* v != 0 */
    3091             : static GEN
    3092         343 : Flc_normalize(GEN v, ulong p)
    3093             : {
    3094         343 :   long i, l = lg(v);
    3095         525 :   for (i = 1; i < l; i++)
    3096         525 :     if (v[i])
    3097             :     {
    3098         343 :       if (v[i] != 1) v = Flv_Fl_div(v, v[i], p);
    3099         343 :       return v;
    3100             :     }
    3101           0 :   return NULL;
    3102             : }
    3103             : 
    3104             : /* K \cap Ker M  [F_l vector spaces]. K = NULL means full space */
    3105             : static GEN
    3106         301 : msfromell_ker(GEN K, GEN M, ulong l)
    3107             : {
    3108         301 :   GEN B, Ml = ZM_to_Flm(M, l);
    3109         301 :   if (K) Ml = Flm_mul(Ml, K, l);
    3110         301 :   B = Flm_ker(Ml, l);
    3111         301 :   if (!K) K = B;
    3112           7 :   else if (lg(B) < lg(K))
    3113           7 :     K = Flm_mul(K, B, l);
    3114         301 :   return K;
    3115             : }
    3116             : /* K = \cap_p Ker(T_p - a_p), 2-dimensional. Set *xl to the 1-dimensional
    3117             :  * Fl-basis  such that star . xl = sign . xl if sign != 0 and
    3118             :  * star * xl[1] = xl[1]; star * xl[2] = -xl[2] if sign = 0 */
    3119             : static void
    3120         294 : msfromell_l(GEN *pxl, GEN K, GEN star, long sign, ulong l)
    3121             : {
    3122         294 :   GEN s = ZM_to_Flm(star, l);
    3123         294 :   GEN a = gel(K,1), Sa = Flm_Flc_mul(s,a,l);
    3124         294 :   GEN b = gel(K,2);
    3125         294 :   GEN t = Flv_add(a,Sa,l), xp, xm;
    3126         294 :   if (zv_equal0(t))
    3127             :   {
    3128          14 :     xm = a;
    3129          14 :     xp = Flv_add(b,Flm_Flc_mul(s,b,l), l);
    3130             :   }
    3131             :   else
    3132             :   {
    3133         280 :     xp = t; t = Flv_sub(a, Sa, l);
    3134         280 :     xm = zv_equal0(t)? Flv_sub(b, Flm_Flc_mul(s,b,l), l): t;
    3135             :   }
    3136             :   /* xp = 0 on Im(S - 1), xm = 0 on Im(S + 1) */
    3137         294 :   if (sign > 0)
    3138         231 :     *pxl = mkmat(Flc_normalize(xp, l));
    3139          63 :   else if (sign < 0)
    3140          14 :     *pxl = mkmat(Flc_normalize(xm, l));
    3141             :   else
    3142          49 :     *pxl = mkmat2(Flc_normalize(xp, l), Flc_normalize(xm, l));
    3143         294 : }
    3144             : static GEN
    3145         294 : msfromell_ratlift(GEN x, GEN q)
    3146             : {
    3147         294 :   GEN B = sqrti(shifti(q,-1));
    3148         294 :   GEN r = FpM_ratlift(x, q, B, B, NULL);
    3149         294 :   if (r) r = Q_primpart(r);
    3150         294 :   return r;
    3151             : }
    3152             : static int
    3153         294 : msfromell_check(GEN x, GEN vT, GEN star, long sign)
    3154             : {
    3155             :   long i, l;
    3156             :   GEN sx;
    3157         294 :   if (!x) return 0;
    3158         294 :   l = lg(vT);
    3159         595 :   for (i = 1; i < l; i++)
    3160             :   {
    3161         301 :     GEN T = gel(vT,i);
    3162         301 :     if (!gequal0(ZM_mul(T, x))) return 0; /* fail */
    3163             :   }
    3164         294 :   sx = ZM_mul(star,x);
    3165         294 :   if (sign)
    3166         245 :     return ZV_equal(gel(sx,1), sign > 0? gel(x,1): ZC_neg(gel(x,1)));
    3167             :   else
    3168          49 :     return ZV_equal(gel(sx,1),gel(x,1)) && ZV_equal(gel(sx,2),ZC_neg(gel(x,2)));
    3169             : }
    3170             : static GEN
    3171         315 : msfromell_scale(GEN E, GEN W, long sign, GEN x)
    3172             : {
    3173         315 :   if (sign)
    3174         245 :     x = ell_get_scale(E, W, gel(x,1), sign);
    3175             :   else
    3176             :   {
    3177          70 :     GEN p = ell_get_scale(E, W, gel(x,1), 1);
    3178          70 :     GEN m = ell_get_scale(E, W, gel(x,2),-1);
    3179          70 :     x = mkvec2(p,m);
    3180             :   }
    3181         315 :   return x;
    3182             : }
    3183             : GEN
    3184         294 : msfromell(GEN E0, long sign)
    3185             : {
    3186         294 :   pari_sp av = avma;
    3187         294 :   GEN E, cond, W, x = NULL, K = NULL, star, q, vT, xl, xr;
    3188             :   long lE, single;
    3189             :   ulong p, l, N;
    3190             :   forprime_t S, Sl;
    3191             : 
    3192         294 :   if (typ(E0) != t_VEC) pari_err_TYPE("msfromell",E0);
    3193         294 :   lE = lg(E0);
    3194         294 :   if (lE == 1) return cgetg(1,t_VEC);
    3195         294 :   single = (typ(gel(E0,1)) != t_VEC);
    3196         294 :   E = single ? E0: gel(E0,1);
    3197             : 
    3198         294 :   E = ellminimalmodel(E, NULL);
    3199         294 :   cond = ellQ_get_N(E);
    3200         294 :   N = itou(cond);
    3201         294 :   W = mskinit(N, 2, 0);
    3202         294 :   star = msk_get_star(W);
    3203         294 :   init_modular_small(&Sl);
    3204             :   /* loop for p <= count_Manin_symbols(N) / 6 would be enough */
    3205         294 :   (void)u_forprime_init(&S, 2, ULONG_MAX);
    3206         294 :   vT = cgetg(1, t_VEC);
    3207         294 :   l = u_forprime_next(&Sl);
    3208         294 :   while( (p = u_forprime_next(&S)) )
    3209             :   {
    3210             :     GEN M;
    3211         343 :     if (N % p == 0) continue;
    3212         301 :     M = RgM_Rg_sub_shallow(mshecke_i(W, p), ellap(E, utoipos(p)));
    3213         301 :     vT = shallowconcat(vT, mkvec(M)); /* for certification at the end */
    3214         301 :     K = msfromell_ker(K, M, l);
    3215         301 :     if (lg(K) == 3) break;
    3216             :   }
    3217         294 :   if (!p) pari_err_BUG("msfromell: ran out of primes");
    3218             : 
    3219             :   /* mod one l should be enough */
    3220         294 :   msfromell_l(&xl, K, star, sign, l);
    3221         294 :   x = ZM_init_CRT(xl, l);
    3222         294 :   q = utoipos(l);
    3223         294 :   xr = msfromell_ratlift(x, q);
    3224             :   /* paranoia */
    3225         588 :   while (!msfromell_check(xr, vT, star, sign) && (l = u_forprime_next(&Sl)) )
    3226             :   {
    3227           0 :     GEN K = NULL;
    3228           0 :     long i, lvT = lg(vT);
    3229           0 :     for (i = 1; i < lvT; i++)
    3230             :     {
    3231           0 :       K = msfromell_ker(K, gel(vT,i), l);
    3232           0 :       if (lg(K) == 3) break;
    3233             :     }
    3234           0 :     if (i >= lvT) { x = NULL; continue; }
    3235           0 :     msfromell_l(&xl, K, star, sign, l);
    3236           0 :     ZM_incremental_CRT(&x, xl, &q, l);
    3237           0 :     xr = msfromell_ratlift(x, q);
    3238             :   }
    3239             :   /* linear form = 0 on all Im(Tp - ap) and Im(S - sign) if sign != 0 */
    3240             : 
    3241         294 :   if (single)
    3242         287 :     x = msfromell_scale(E, W, sign, xr);
    3243             :   else
    3244             :   {
    3245           7 :     GEN v = cgetg(lE, t_VEC);
    3246             :     long i;
    3247           7 :     for (i=1; i < lE; i++) gel(v,i) = msfromell_scale(gel(E0,i), W, sign, xr);
    3248           7 :     x = v;
    3249             :   }
    3250         294 :   return gerepilecopy(av, mkvec2(W, x));
    3251             : }
    3252             : 
    3253             : GEN
    3254          14 : msfromhecke(GEN W, GEN v, GEN H)
    3255             : {
    3256          14 :   pari_sp av = avma;
    3257          14 :   long i, l = lg(v);
    3258          14 :   GEN K = NULL;
    3259          14 :   checkms(W);
    3260          14 :   if (typ(v) != t_VEC) pari_err_TYPE("msfromhecke",v);
    3261          35 :   for (i = 1; i < l; i++)
    3262             :   {
    3263          21 :     GEN K2, T, p, P, c = gel(v,i);
    3264          21 :     if (typ(c) != t_VEC || lg(c) != 3) pari_err_TYPE("msfromhecke",v);
    3265          21 :     p = gel(c,1);
    3266          21 :     if (typ(p) != t_INT) pari_err_TYPE("msfromhecke",v);
    3267          21 :     P = gel(c,2);
    3268          21 :     switch(typ(P))
    3269             :     {
    3270             :       case t_INT:
    3271          14 :         P = deg1pol_shallow(gen_1, negi(P), 0);
    3272          14 :         break;
    3273             :       case t_POL:
    3274           7 :         if (RgX_is_ZX(P)) break;
    3275             :       default:
    3276           0 :         pari_err_TYPE("msfromhecke",v);
    3277             :     };
    3278          21 :     T = mshecke(W, itos(p), H);
    3279          21 :     T = Q_primpart(RgX_RgM_eval(P, T));
    3280          21 :     if (K) T = ZM_mul(T,K);
    3281          21 :     K2 = ZM_ker(T);
    3282          21 :     if (!K) K = K2;
    3283           7 :     else if (lg(K2) < lg(K)) K = ZM_mul(K,K2);
    3284             :   }
    3285          14 :   return gerepilecopy(av, K);
    3286             : }
    3287             : 
    3288             : /* OVERCONVERGENT MODULAR SYMBOLS */
    3289             : 
    3290             : static GEN
    3291        2765 : mspadic_get_Wp(GEN W) { return gel(W,1); }
    3292             : static GEN
    3293         455 : mspadic_get_Tp(GEN W) { return gel(W,2); }
    3294             : static GEN
    3295         455 : mspadic_get_bin(GEN W) { return gel(W,3); }
    3296             : static GEN
    3297         448 : mspadic_get_actUp(GEN W) { return gel(W,4); }
    3298             : static GEN
    3299         448 : mspadic_get_q(GEN W) { return gel(W,5); }
    3300             : static long
    3301        1372 : mspadic_get_p(GEN W) { return gel(W,6)[1]; }
    3302             : static long
    3303        1148 : mspadic_get_n(GEN W) { return gel(W,6)[2]; }
    3304             : static long
    3305         161 : mspadic_get_flag(GEN W) { return gel(W,6)[3]; }
    3306             : static GEN
    3307         455 : mspadic_get_M(GEN W) { return gel(W,7); }
    3308             : static GEN
    3309         455 : mspadic_get_C(GEN W) { return gel(W,8); }
    3310             : static long
    3311         917 : mspadic_get_weight(GEN W) { return msk_get_weight(mspadic_get_Wp(W)); }
    3312             : 
    3313             : void
    3314         924 : checkmspadic(GEN W)
    3315             : {
    3316         924 :   if (typ(W) != t_VEC || lg(W) != 9) pari_err_TYPE("checkmspadic",W);
    3317         924 :   checkms(mspadic_get_Wp(W));
    3318         924 : }
    3319             : 
    3320             : /* f in M_2(Z) \cap GL_2(Q), p \nmid a [ and for the result to mean anything
    3321             :  * p | c, but not needed here]. Return the matrix M in M_D(Z), D = M+k-1
    3322             :  * such that, if v = \int x^i d mu, i < D, is a vector of D moments of mu,
    3323             :  * then M * v is the vector of moments of mu | f  mod p^D */
    3324             : static GEN
    3325      252917 : moments_act(struct m_act *S, GEN f)
    3326             : {
    3327      252917 :   pari_sp av = avma;
    3328      252917 :   long j, k = S->k, D = S->dim;
    3329      252917 :   GEN a = gcoeff(f,1,1), b = gcoeff(f,1,2);
    3330      252917 :   GEN c = gcoeff(f,2,1), d = gcoeff(f,2,2);
    3331      252917 :   GEN u,z,C, q = S->q, mat = cgetg(D+1, t_MAT);
    3332             : 
    3333      252917 :   a = modii(a,q);
    3334      252917 :   z = FpX_powu(deg1pol(c,a,0), k-2, q); /* (a+cx)^(k-2) */
    3335             :   /* u := (b+dx) / (a+cx) mod (q,x^D) = (b/a +d/a*x) / (1 - (-c/a)*x) */
    3336      252917 :   if (!equali1(a))
    3337             :   {
    3338      248696 :     GEN ai = Fp_inv(a,q);
    3339      248696 :     b = Fp_mul(b,ai,q);
    3340      248696 :     c = Fp_mul(c,ai,q);
    3341      248696 :     d = Fp_mul(d,ai,q);
    3342             :   }
    3343      252917 :   u = cgetg(D+2,t_POL); u[1] = evalsigne(1)|evalvarn(0);
    3344      252917 :   gel(u, 2) = gen_1;
    3345      252917 :   gel(u, 3) = C = Fp_neg(c,q);
    3346      252917 :   for (j = 4; j < D+2; j++) gel(u,j) = Fp_mul(gel(u,j-1), C, q);
    3347      252917 :   u = FpX_red(RgXn_mul(deg1pol(d,b,0), u, D), q);
    3348     2120272 :   for (j = 1; j <= D; j++)
    3349             :   {
    3350     1867355 :     gel(mat,j) = RgX_to_RgC(z, D); /* (a+cx)^(k-2) * ((b+dx)/(a+cx))^(j-1) */
    3351     1867355 :     if (j != D) z = FpX_red(RgXn_mul(z, u, D), q);
    3352             :   }
    3353      252917 :   return gerepilecopy(av, shallowtrans(mat));
    3354             : }
    3355             : 
    3356             : static GEN
    3357         455 : init_moments_act(GEN W, long p, long n, GEN q, GEN v)
    3358             : {
    3359             :   struct m_act S;
    3360         455 :   long k = msk_get_weight(W);
    3361         455 :   S.p = p;
    3362         455 :   S.k = k;
    3363         455 :   S.q = q;
    3364         455 :   S.dim = n+k-1;
    3365         455 :   return init_dual_act(v,W,W,&S, moments_act);
    3366             : }
    3367             : 
    3368             : static void
    3369        6552 : clean_tail(GEN phi, long c, GEN q)
    3370             : {
    3371        6552 :   long a, l = lg(phi);
    3372      208418 :   for (a = 1; a < l; a++)
    3373             :   {
    3374      201866 :     GEN P = FpV_red(gel(phi,a), q); /* phi(G_a) = vector of moments */
    3375      201866 :     long j, lP = lg(P);
    3376      201866 :     for (j = c; j < lP; j++) gel(P,j) = gen_0; /* reset garbage to 0 */
    3377      201866 :     gel(phi,a) = P;
    3378             :   }
    3379        6552 : }
    3380             : /* concat z to all phi[i] */
    3381             : static GEN
    3382         602 : concat2(GEN phi, GEN z)
    3383             : {
    3384             :   long i, l;
    3385         602 :   GEN v = cgetg_copy(phi,&l);
    3386         602 :   for (i = 1; i < l; i++) gel(v,i) = shallowconcat(gel(phi,i), z);
    3387         602 :   return v;
    3388             : }
    3389             : static GEN
    3390         602 : red_mod_FilM(GEN phi, ulong p, long k, long flag)
    3391             : {
    3392             :   long a, l;
    3393         602 :   GEN den = gen_1, v = cgetg_copy(phi, &l);
    3394         602 :   if (flag)
    3395             :   {
    3396         343 :     phi = Q_remove_denom(phi, &den);
    3397         343 :     if (!den) { den = gen_1; flag = 0; }
    3398             :   }
    3399       28630 :   for (a = 1; a < l; a++)
    3400             :   {
    3401       28028 :     GEN P = gel(phi,a), q = den;
    3402             :     long j;
    3403      201866 :     for (j = lg(P)-1; j >= k+1; j--)
    3404             :     {
    3405      173838 :       q = muliu(q,p);
    3406      173838 :       gel(P,j) = modii(gel(P,j),q);
    3407             :     }
    3408       28028 :     q = muliu(q,p);
    3409       91196 :     for (     ; j >= 1; j--)
    3410       63168 :       gel(P,j) = modii(gel(P,j),q);
    3411       28028 :     gel(v,a) = P;
    3412             :   }
    3413         602 :   if (flag) v = gdiv(v, den);
    3414         602 :   return v;
    3415             : }
    3416             : 
    3417             : /* denom(C) | p^(2(k-1) - v_p(ap)) */
    3418             : static GEN
    3419         154 : oms_dim2(GEN W, GEN phi, GEN C, GEN ap)
    3420             : {
    3421         154 :   long t, i, k = mspadic_get_weight(W);
    3422         154 :   long p = mspadic_get_p(W), n = mspadic_get_n(W);
    3423         154 :   GEN phi1 = gel(phi,1), phi2 = gel(phi,2);
    3424         154 :   GEN v, q = mspadic_get_q(W);
    3425         154 :   GEN act = mspadic_get_actUp(W);
    3426             : 
    3427         154 :   t = signe(ap)? Z_lval(ap,p) : k-1;
    3428         154 :   phi1 = concat2(phi1, zerovec(n));
    3429         154 :   phi2 = concat2(phi2, zerovec(n));
    3430        2107 :   for (i = 1; i <= n; i++)
    3431             :   {
    3432        1953 :     phi1 = dual_act(k-1, act, phi1);
    3433        1953 :     phi1 = dual_act(k-1, act, phi1);
    3434        1953 :     clean_tail(phi1, k + i*t, q);
    3435             : 
    3436        1953 :     phi2 = dual_act(k-1, act, phi2);
    3437        1953 :     phi2 = dual_act(k-1, act, phi2);
    3438        1953 :     clean_tail(phi2, k + i*t, q);
    3439             :   }
    3440         154 :   C = gpowgs(C,n);
    3441         154 :   v = RgM_RgC_mul(C, mkcol2(phi1,phi2));
    3442         154 :   phi1 = red_mod_FilM(gel(v,1), p, k, 1);
    3443         154 :   phi2 = red_mod_FilM(gel(v,2), p, k, 1);
    3444         154 :   return mkvec2(phi1,phi2);
    3445             : }
    3446             : 
    3447             : /* flag = 0 iff alpha is a p-unit */
    3448             : static GEN
    3449         294 : oms_dim1(GEN W, GEN phi, GEN alpha, long flag)
    3450             : {
    3451         294 :   long i, k = mspadic_get_weight(W);
    3452         294 :   long p = mspadic_get_p(W), n = mspadic_get_n(W);
    3453         294 :   GEN q = mspadic_get_q(W);
    3454         294 :   GEN act = mspadic_get_actUp(W);
    3455         294 :   phi = concat2(phi, zerovec(n));
    3456        2940 :   for (i = 1; i <= n; i++)
    3457             :   {
    3458        2646 :     phi = dual_act(k-1, act, phi);
    3459        2646 :     clean_tail(phi, k + i, q);
    3460             :   }
    3461         294 :   phi = gmul(lift_shallow(gpowgs(alpha,n)), phi);
    3462         294 :   phi = red_mod_FilM(phi, p, k, flag);
    3463         294 :   return mkvec(phi);
    3464             : }
    3465             : 
    3466             : /* lift polynomial P in RgX[X,Y]_{k-2} to a distribution \mu such that
    3467             :  * \int (Y - X z)^(k-2) d\mu(z) = P(X,Y)
    3468             :  * Return the t_VEC of k-1 first moments of \mu: \int z^i d\mu(z), 0<= i < k-1.
    3469             :  *   \sum_j (-1)^(k-2-j) binomial(k-2,j) Y^j \int z^(k-2-j) d\mu(z) = P(1,Y)
    3470             :  * Input is P(1,Y), bin = vecbinomial(k-2): bin[j] = binomial(k-2,j-1) */
    3471             : static GEN
    3472       37667 : RgX_to_moments(GEN P, GEN bin)
    3473             : {
    3474       37667 :   long j, k = lg(bin);
    3475             :   GEN Pd, Bd;
    3476       37667 :   if (typ(P) != t_POL) P = scalarpol(P,0);
    3477       37667 :   P = RgX_to_RgC(P, k-1); /* deg <= k-2 */
    3478       37667 :   settyp(P, t_VEC);
    3479       37667 :   Pd = P+1;  /* Pd[i] = coeff(P,i) */
    3480       37667 :   Bd = bin+1;/* Bd[i] = binomial(k-2,i) */
    3481       45290 :   for (j = 1; j < k-2; j++)
    3482             :   {
    3483        7623 :     GEN c = gel(Pd,j);
    3484        7623 :     if (odd(j)) c = gneg(c);
    3485        7623 :     gel(Pd,j) = gdiv(c, gel(Bd,j));
    3486             :   }
    3487       37667 :   return vecreverse(P);
    3488             : }
    3489             : static GEN
    3490         847 : RgXC_to_moments(GEN v, GEN bin)
    3491             : {
    3492             :   long i, l;
    3493         847 :   GEN w = cgetg_copy(v,&l);
    3494         847 :   for (i=1; i<l; i++) gel(w,i) = RgX_to_moments(gel(v,i),bin);
    3495         847 :   return w;
    3496             : }
    3497             : 
    3498             : /* W an mspadic, assume O[2] is integral, den is the cancelled denominator
    3499             :  * or NULL, L = log(path) */
    3500             : static GEN
    3501        2534 : omseval_int(struct m_act *S, GEN PHI, GEN L, hashtable *H)
    3502             : {
    3503             :   long a, lphi;
    3504        2534 :   GEN ind, v = cgetg_copy(PHI, &lphi);
    3505             : 
    3506        2534 :   L = RgV_sparse(L,&ind);
    3507        2534 :   ZSl2C_star_inplace(L); /* lambda_{i,j}^* */
    3508        2534 :   L = mkvec2(ind,L);
    3509        2534 :   ZGl2QC_to_act(S, moments_act, L, H); /* as operators on V */
    3510        5446 :   for (a = 1; a < lphi; a++)
    3511             :   {
    3512        2912 :     GEN T = dense_act_col(L, gel(PHI,a));
    3513        2912 :     if (T) T = FpC_red(T,S->q); else T = zerocol(S->dim);
    3514        2912 :     gel(v,a) = T;
    3515             :   }
    3516        2534 :   return v;
    3517             : }
    3518             : 
    3519             : GEN
    3520          14 : msomseval(GEN W, GEN phi, GEN path)
    3521             : {
    3522             :   struct m_act S;
    3523          14 :   pari_sp av = avma;
    3524             :   GEN v, Wp;
    3525             :   long n, vden;
    3526          14 :   checkmspadic(W);
    3527          14 :   if (typ(phi) != t_COL || lg(phi) != 4)  pari_err_TYPE("msomseval",phi);
    3528          14 :   vden = itos(gel(phi,2));
    3529          14 :   phi = gel(phi,1);
    3530          14 :   n = mspadic_get_n(W);
    3531          14 :   Wp= mspadic_get_Wp(W);
    3532          14 :   S.k = mspadic_get_weight(W);
    3533          14 :   S.p = mspadic_get_p(W);
    3534          14 :   S.q = powuu(S.p, n+vden);
    3535          14 :   S.dim = n + S.k - 1;
    3536          14 :   v = omseval_int(&S, phi, mspathlog(Wp,path), NULL);
    3537          14 :   return gerepilecopy(av, v);
    3538             : }
    3539             : /* W = msinit(N,k,...); if flag < 0 or flag >= k-1, allow all symbols;
    3540             :  * else commit to v_p(a_p) <= flag (ordinary if flag = 0)*/
    3541             : GEN
    3542         462 : mspadicinit(GEN W, long p, long n, long flag)
    3543             : {
    3544         462 :   pari_sp av = avma;
    3545             :   long a, N, k;
    3546             :   GEN P, C, M, bin, Wp, Tp, q, pn, actUp, teich, pas;
    3547             : 
    3548         462 :   checkms(W);
    3549         462 :   N = ms_get_N(W);
    3550         462 :   k = msk_get_weight(W);
    3551         462 :   if (flag < 0) flag = 1; /* worst case */
    3552         343 :   else if (flag >= k) flag = k-1;
    3553             : 
    3554         462 :   bin = vecbinomial(k-2);
    3555         462 :   Tp = mshecke(W, p, NULL);
    3556         462 :   if (N % p == 0)
    3557             :   {
    3558          70 :     if ((N/p) % p == 0) pari_err_IMPL("mspadicinit when p^2 | N");
    3559             :     /* a_p != 0 */
    3560          63 :     Wp = W;
    3561          63 :     M = gen_0;
    3562          63 :     flag = (k-2) / 2; /* exact valuation */
    3563             :     /* will multiply by matrix with denominator p^(k-2)/2 in mspadicint.
    3564             :      * Except if p = 2 (multiply by alpha^2) */
    3565          63 :     if (p == 2) n += k-2; else n += (k-2)/2;
    3566          63 :     pn = powuu(p,n);
    3567             :     /* For accuracy mod p^n, oms_dim1 require p^(k/2*n) */
    3568          63 :     q = powiu(pn, k/2);
    3569             :   }
    3570             :   else
    3571             :   { /* p-stabilize */
    3572         392 :     long s = msk_get_sign(W);
    3573             :     GEN M1, M2;
    3574             : 
    3575         392 :     Wp = mskinit(N*p, k, s);
    3576         392 :     M1 = getMorphism(W, Wp, mkvec(mat2(1,0,0,1)));
    3577         392 :     M2 = getMorphism(W, Wp, mkvec(mat2(p,0,0,1)));
    3578         392 :     if (s)
    3579             :     {
    3580         147 :       GEN SW = msk_get_starproj(W), SWp = msk_get_starproj(Wp);
    3581         147 :       M1 = Qevproj_apply2(M1, SW, SWp);
    3582         147 :       M2 = Qevproj_apply2(M2, SW, SWp);
    3583             :     }
    3584         392 :     M = mkvec2(M1,M2);
    3585         392 :     n += Z_lval(Q_denom(M), p); /*den. introduced by p-stabilization*/
    3586             :     /* in supersingular case: will multiply by matrix with denominator p^k
    3587             :      * in mspadicint. Except if p = 2 (multiply by alpha^2) */
    3588         392 :     if (flag) { if (p == 2) n += 2*k-2; else n += k; }
    3589         392 :     pn = powuu(p,n);
    3590             :     /* For accuracy mod p^n, supersingular require p^((2k-1-v_p(a_p))*n) */
    3591         392 :     if (flag) /* k-1 also takes care of a_p = 0. Worst case v_p(a_p) = flag */
    3592         231 :       q = powiu(pn, 2*k-1 - flag);
    3593             :     else
    3594         161 :       q = pn;
    3595             :   }
    3596         455 :   actUp = init_moments_act(Wp, p, n, q, Up_matrices(p));
    3597             : 
    3598         455 :   if (p == 2) C = gen_0;
    3599             :   else
    3600             :   {
    3601         399 :     pas = matpascal(n);
    3602         399 :     teich = teichmullerinit(p, n+1);
    3603         399 :     P = gpowers(utoipos(p), n);
    3604         399 :     C = cgetg(p, t_VEC);
    3605        1911 :     for (a = 1; a < p; a++)
    3606             :     { /* powb[j+1] = ((a - w(a)) / p)^j mod p^n */
    3607        1512 :       GEN powb = Fp_powers(diviuexact(subui(a, gel(teich,a)), p), n, pn);
    3608        1512 :       GEN Ca = cgetg(n+2, t_VEC);
    3609        1512 :       long j, r, ai = Fl_inv(a, p); /* a^(-1) */
    3610        1512 :       gel(C,a) = Ca;
    3611       18018 :       for (j = 0; j <= n; j++)
    3612             :       {
    3613       16506 :         GEN Caj = cgetg(j+2, t_VEC);
    3614       16506 :         GEN atij = gel(teich, Fl_powu(ai,j,p));/* w(a)^(-j) = w(a^(-j) mod p) */
    3615       16506 :         gel(Ca,j+1) = Caj;
    3616      133294 :         for (r = 0; r <= j; r++)
    3617             :         {
    3618      116788 :           GEN c = Fp_mul(gcoeff(pas,j+1,r+1), gel(powb, j-r+1), pn);
    3619      116788 :           c = Fp_mul(c,atij,pn); /* binomial(j,r)*b^(j-r)*w(a)^(-j) mod p^n */
    3620      116788 :           gel(Caj,r+1) = mulii(c, gel(P,j+1)); /* p^j * c mod p^(n+j) */
    3621             :         }
    3622             :       }
    3623             :     }
    3624             :   }
    3625         455 :   return gerepilecopy(av, mkvecn(8, Wp,Tp, bin, actUp, q,
    3626             :                                  mkvecsmall3(p,n,flag), M, C));
    3627             : }
    3628             : 
    3629             : #if 0
    3630             : /* assume phi an ordinary OMS */
    3631             : static GEN
    3632             : omsactgl2(GEN W, GEN phi, GEN M)
    3633             : {
    3634             :   GEN q, Wp, act;
    3635             :   long p, k, n;
    3636             :   checkmspadic(W);
    3637             :   Wp = mspadic_get_Wp(W);
    3638             :   p = mspadic_get_p(W);
    3639             :   k = mspadic_get_weight(W);
    3640             :   n = mspadic_get_n(W);
    3641             :   q = mspadic_get_q(W);
    3642             :   act = init_moments_act(Wp, p, n, q, M);
    3643             :   phi = gel(phi,1);
    3644             :   return dual_act(k-1, act, gel(phi,1));
    3645             : }
    3646             : #endif
    3647             : 
    3648             : static GEN
    3649         455 : eigenvalue(GEN T, GEN x)
    3650             : {
    3651         455 :   long i, l = lg(x);
    3652         581 :   for (i = 1; i < l; i++)
    3653         581 :     if (!isintzero(gel(x,i))) break;
    3654         455 :   if (i == l) pari_err_DOMAIN("mstooms", "phi", "=", gen_0, x);
    3655         455 :   return gdiv(RgMrow_RgC_mul(T,x,i), gel(x,i));
    3656             : }
    3657             : 
    3658             : /* p coprime to ap, return unit root of x^2 - ap*x + p^(k-1), accuracy p^n */
    3659             : static GEN
    3660         231 : ms_unit_eigenvalue(GEN ap, long k, GEN p, long n)
    3661             : {
    3662         231 :   GEN sqrtD, D = subii(sqri(ap), shifti(powiu(p,k-1),2));
    3663         231 :   if (absequaliu(p,2))
    3664             :   {
    3665           7 :     n++; sqrtD = Zp_sqrt(D, p, n);
    3666           7 :     if (mod4(sqrtD) != mod4(ap)) sqrtD = negi(sqrtD);
    3667             :   }
    3668             :   else
    3669         224 :     sqrtD = Zp_sqrtlift(D, ap, p, n);
    3670             :   /* sqrtD = ap (mod p) */
    3671         231 :   return gmul2n(gadd(ap, cvtop(sqrtD,p,n)), -1);
    3672             : }
    3673             : 
    3674             : /* W = msinit(N,k,...); phi = T_p/U_p - eigensymbol */
    3675             : GEN
    3676         455 : mstooms(GEN W, GEN phi)
    3677             : {
    3678         455 :   pari_sp av = avma;
    3679             :   GEN Wp, bin, Tp, c, alpha, ap, phi0, M;
    3680             :   long k, p, vden;
    3681             : 
    3682         455 :   checkmspadic(W);
    3683         455 :   if (typ(phi) != t_COL)
    3684             :   {
    3685         161 :     if (!is_Qevproj(phi)) pari_err_TYPE("mstooms",phi);
    3686         161 :     phi = gel(phi,1);
    3687         161 :     if (lg(phi) != 2) pari_err_TYPE("mstooms [dim_Q (eigenspace) > 1]",phi);
    3688         161 :     phi = gel(phi,1);
    3689             :   }
    3690             : 
    3691         455 :   Wp = mspadic_get_Wp(W);
    3692         455 :   Tp = mspadic_get_Tp(W);
    3693         455 :   bin = mspadic_get_bin(W);
    3694         455 :   k = msk_get_weight(Wp);
    3695         455 :   p = mspadic_get_p(W);
    3696         455 :   M = mspadic_get_M(W);
    3697             : 
    3698         455 :   phi = Q_remove_denom(phi, &c);
    3699         455 :   ap = eigenvalue(Tp, phi);
    3700         455 :   vden = c? Z_lvalrem(c, p, &c): 0;
    3701             : 
    3702         455 :   if (typ(M) == t_INT)
    3703             :   { /* p | N */
    3704             :     GEN c1;
    3705          63 :     alpha = ap;
    3706          63 :     alpha = ginv(alpha);
    3707          63 :     phi0 = mseval(Wp, phi, NULL);
    3708          63 :     phi0 = RgXC_to_moments(phi0, bin);
    3709          63 :     phi0 = Q_remove_denom(phi0, &c1);
    3710          63 :     if (c1) { vden += Z_lvalrem(c1, p, &c1); c = mul_denom(c,c1); }
    3711          63 :     if (umodiu(ap,p)) /* p \nmid a_p */
    3712          28 :       phi = oms_dim1(W, phi0, alpha, 0);
    3713             :     else
    3714             :     {
    3715          35 :       phi = oms_dim1(W, phi0, alpha, 1);
    3716          35 :       phi = Q_remove_denom(phi, &c1);
    3717          35 :       if (c1) { vden += Z_lvalrem(c1, p, &c1); c = mul_denom(c,c1); }
    3718             :     }
    3719             :   }
    3720             :   else
    3721             :   { /* p-stabilize */
    3722             :     GEN M1, M2, phi1, phi2, c1;
    3723         392 :     if (typ(M) != t_VEC || lg(M) != 3) pari_err_TYPE("mstooms",W);
    3724         392 :     M1 = gel(M,1);
    3725         392 :     M2 = gel(M,2);
    3726             : 
    3727         392 :     phi1 = RgM_RgC_mul(M1, phi);
    3728         392 :     phi2 = RgM_RgC_mul(M2, phi);
    3729         392 :     phi1 = mseval(Wp, phi1, NULL);
    3730         392 :     phi2 = mseval(Wp, phi2, NULL);
    3731             : 
    3732         392 :     phi1 = RgXC_to_moments(phi1, bin);
    3733         392 :     phi2 = RgXC_to_moments(phi2, bin);
    3734         392 :     phi = Q_remove_denom(mkvec2(phi1,phi2), &c1);
    3735         392 :     phi1 = gel(phi,1);
    3736         392 :     phi2 = gel(phi,2);
    3737         392 :     if (c1) { vden += Z_lvalrem(c1, p, &c1); c = mul_denom(c,c1); }
    3738             :     /* all polynomials multiplied by c p^vden */
    3739         392 :     if (umodiu(ap, p))
    3740             :     {
    3741         231 :       alpha = ms_unit_eigenvalue(ap, k, utoipos(p), mspadic_get_n(W));
    3742         231 :       alpha = ginv(alpha);
    3743         231 :       phi0 = gsub(phi1, gmul(lift_shallow(alpha),phi2));
    3744         231 :       phi = oms_dim1(W, phi0, alpha, 0);
    3745             :     }
    3746             :     else
    3747             :     { /* p | ap, alpha = [a_p, -1; p^(k-1), 0] */
    3748         161 :       long flag = mspadic_get_flag(W);
    3749         161 :       if (!flag || (signe(ap) && Z_lval(ap,p) < flag))
    3750           7 :         pari_err_TYPE("mstooms [v_p(ap) > mspadicinit flag]", phi);
    3751         154 :       alpha = mkmat22(ap,gen_m1, powuu(p, k-1),gen_0);
    3752         154 :       alpha = ginv(alpha);
    3753         154 :       phi = oms_dim2(W, mkvec2(phi1,phi2), gsqr(alpha), ap);
    3754         154 :       phi = Q_remove_denom(phi, &c1);
    3755         154 :       if (c1) { vden += Z_lvalrem(c1, p, &c1); c = mul_denom(c,c1); }
    3756             :     }
    3757             :   }
    3758         448 :   if (vden) c = mul_denom(c, powuu(p,vden));
    3759         448 :   if (p == 2) alpha = gsqr(alpha);
    3760         448 :   if (c) alpha = gdiv(alpha,c);
    3761         448 :   if (typ(alpha) == t_MAT)
    3762             :   { /* express in basis (omega,-p phi(omega)) */
    3763         154 :     gcoeff(alpha,2,1) = gdivgs(gcoeff(alpha,2,1), -p);
    3764         154 :     gcoeff(alpha,2,2) = gdivgs(gcoeff(alpha,2,2), -p);
    3765             :     /* at the end of mspadicint we shall multiply result by [1,0;0,-1/p]*alpha
    3766             :      * vden + k is the denominator of this matrix */
    3767             :   }
    3768             :   /* phi is integral-valued */
    3769         448 :   return gerepilecopy(av, mkcol3(phi, stoi(vden), alpha));
    3770             : }
    3771             : 
    3772             : /* HACK: the v[j] have different lengths */
    3773             : static GEN
    3774        1778 : FpVV_dotproduct(GEN v, GEN w, GEN p)
    3775             : {
    3776        1778 :   long j, l = lg(v);
    3777        1778 :   GEN T = cgetg(l, t_VEC);
    3778        1778 :   for (j = 1; j < l; j++) gel(T,j) = FpV_dotproduct(gel(v,j),w,p);
    3779        1778 :   return T;
    3780             : }
    3781             : 
    3782             : /* \int (-4z)^j given \int z^j */
    3783             : static GEN
    3784          98 : twistmoment_minus(GEN v)
    3785             : {
    3786             :   long i, l;
    3787          98 :   GEN w = cgetg_copy(v, &l);
    3788        2009 :   for (i = 1; i < l; i++)
    3789             :   {
    3790        1911 :     GEN c = gel(v,i);
    3791        1911 :     if (i > 1) c = gmul2n(c, (i-1)<<1);
    3792        1911 :     gel(w,i) = odd(i)? c: gneg(c);
    3793             :   }
    3794          98 :   return w;
    3795             : }
    3796             : /* \int (4z)^j given \int z^j */
    3797             : static GEN
    3798          98 : twistmoment_plus(GEN v)
    3799             : {
    3800             :   long i, l;
    3801          98 :   GEN w = cgetg_copy(v, &l);
    3802        2009 :   for (i = 1; i < l; i++)
    3803             :   {
    3804        1911 :     GEN c = gel(v,i);
    3805        1911 :     if (i > 1) c = gmul2n(c, (i-1)<<1);
    3806        1911 :     gel(w,i) = c;
    3807             :   }
    3808          98 :   return w;
    3809             : }
    3810             : /* W an mspadic, phi eigensymbol, p \nmid D. Return C(x) mod FilM */
    3811             : GEN
    3812         455 : mspadicmoments(GEN W, GEN PHI, long D)
    3813             : {
    3814         455 :   pari_sp av = avma;
    3815         455 :   long la, ia, b, lphi, aD = labs(D), pp, p, k, n, vden;
    3816             :   GEN Wp, Dact, Dk, v, C, gp, pn, phi;
    3817             :   struct m_act S;
    3818             :   hashtable *H;
    3819             : 
    3820         455 :   checkmspadic(W);
    3821         455 :   Wp = mspadic_get_Wp(W);
    3822         455 :   p = mspadic_get_p(W);
    3823         455 :   k = mspadic_get_weight(W);
    3824         455 :   n = mspadic_get_n(W);
    3825         455 :   C = mspadic_get_C(W);
    3826         455 :   if (typ(PHI) != t_COL || lg(PHI) != 4 || typ(gel(PHI,1)) != t_VEC)
    3827         448 :     PHI = mstooms(W, PHI);
    3828         448 :   vden = itos( gel(PHI,2) );
    3829         448 :   phi = gel(PHI,1);
    3830         448 :   if (p == 2)
    3831          56 :   { la = 3; pp = 4; }
    3832             :   else
    3833         392 :   { la = p; pp = p; }
    3834         448 :   v = cgetg_copy(phi, &lphi);
    3835         448 :   for (b = 1; b < lphi; b++) gel(v,b) = cgetg(la, t_VEC);
    3836         448 :   pn = powuu(p, n + vden);
    3837         448 :   gp = utoipos(p);
    3838             : 
    3839         448 :   S.p = p;
    3840         448 :   S.k = k;
    3841         448 :   S.q = pn;
    3842         448 :   S.dim = n+k-1;
    3843             : 
    3844         448 :   Dact = NULL;
    3845         448 :   Dk = NULL;
    3846         448 :   if (D != 1)
    3847             :   {
    3848          56 :     GEN gaD = utoi(aD);
    3849          56 :     if (!sisfundamental(D)) pari_err_TYPE("mspadicmoments", stoi(D));
    3850          56 :     if (D % p == 0) pari_err_DOMAIN("mspadicmoments", "p","|", stoi(D), gp);
    3851          56 :     Dact = cgetg(aD, t_VEC);
    3852         504 :     for (b = 1; b < aD; b++)
    3853             :     {
    3854         448 :       GEN z = NULL;
    3855         448 :       if (ugcd(b,aD) == 1)
    3856         448 :         z = moments_act(&S, mkmat22(gaD,utoipos(b), gen_0,gaD));
    3857         448 :       gel(Dact,b) = z;
    3858             :     }
    3859          56 :     if (k != 2) Dk = Fp_pows(stoi(D), 2-k, pn);
    3860             :   }
    3861             : 
    3862         448 :   H = Gl2act_cache(ms_get_nbgen(Wp));
    3863             : 
    3864        2058 :   for (ia = 1; ia < la; ia++)
    3865             :   {
    3866             :     GEN path, vca;
    3867        1610 :     long i, a = ia;
    3868        1610 :     if (p == 2 && a == 2) a = -1;
    3869        1610 :     if (Dact) /* twist by D */
    3870             :     {
    3871             :       long c;
    3872         182 :       vca = const_vec(lphi-1,NULL);
    3873        1274 :       for (b = 1; b < aD; b++)
    3874             :       {
    3875        1092 :         long s = kross(D, b);
    3876             :         GEN z, T;
    3877        1092 :         if (!s) continue;
    3878        1092 :         z = addii(mulss(a, aD), muluu(pp, b));
    3879             :         /* oo -> a/pp + pp/|D|*/
    3880        1092 :         path = mkmat22(gen_1,z, gen_0,muluu(pp, aD));
    3881        1092 :         T = omseval_int(&S, phi, M2_log(Wp,path), H);
    3882        2184 :         for (c = 1; c < lphi; c++)
    3883             :         {
    3884        1092 :           z = FpM_FpC_mul(gel(Dact,b), gel(T,c), pn);
    3885        1092 :           if (s < 0) ZV_neg_inplace(z);
    3886        1092 :           gel(vca, c) = gel(vca,c)? ZC_add(gel(vca,c), z): z;
    3887             :         }
    3888             :       }
    3889         252 :       if (Dk) for(c = 1; c < lphi; c++)
    3890          70 :         gel(vca,c) = FpC_Fp_mul(gel(vca,c), Dk, pn);
    3891             :     }
    3892             :     else
    3893             :     {
    3894        1428 :       path = mkmat22(gen_1,stoi(a), gen_0, utoipos(pp));
    3895        1428 :       vca = omseval_int(&S, phi, M2_log(Wp,path), H);
    3896             :     }
    3897        1610 :     if (p != 2)
    3898             :     {
    3899        1498 :       GEN Ca = gel(C,a);
    3900        3276 :       for (i = 1; i < lphi; i++)
    3901        1778 :         gmael(v,i,a) = FpVV_dotproduct(Ca, gel(vca,i), pn);
    3902             :     }
    3903             :     else
    3904             :     {
    3905         112 :       if (ia == 1) /* \tilde{a} = 1 */
    3906          56 :       { for (i = 1; i < lphi; i++) gel(vca,i) = twistmoment_plus(gel(vca,i)); }
    3907             :       else /* \tilde{a} = -1 */
    3908          56 :       { for (i = 1; i < lphi; i++) gel(vca,i) = twistmoment_minus(gel(vca,i)); }
    3909         112 :       for (i = 1; i < lphi; i++) gmael(v,i,ia) = gel(vca,i);
    3910             :     }
    3911             :   }
    3912         448 :   return gerepilecopy(av, mkvec3(v, gel(PHI,3), mkvecsmall4(p,n+vden,n,D)));
    3913             : }
    3914             : static void
    3915        1855 : checkoms(GEN v)
    3916             : {
    3917        1855 :   if (typ(v) != t_VEC || lg(v) != 4 || typ(gel(v,1)) != t_VEC
    3918        1855 :       || typ(gel(v,3))!=t_VECSMALL)
    3919           0 :     pari_err_TYPE("checkoms [apply mspadicmoments]", v);
    3920        1855 : }
    3921             : static long
    3922        4158 : oms_get_p(GEN oms) { return gel(oms,3)[1]; }
    3923             : static long
    3924        4060 : oms_get_n(GEN oms) { return gel(oms,3)[2]; }
    3925             : static long
    3926        2401 : oms_get_n0(GEN oms) { return gel(oms,3)[3]; }
    3927             : static long
    3928        1855 : oms_get_D(GEN oms) { return gel(oms,3)[4]; }
    3929             : static int
    3930          98 : oms_is_supersingular(GEN oms) { GEN v = gel(oms,1); return lg(v) == 3; }
    3931             : 
    3932             : /* sum(j = 1, n, (-1)^(j+1)/j * x^j) */
    3933             : static GEN
    3934         749 : log1x(long n)
    3935             : {
    3936         749 :   long i, l = n+3;
    3937         749 :   GEN v = cgetg(l, t_POL);
    3938         749 :   v[1] = evalvarn(0)|evalsigne(1); gel(v,2) = gen_0;
    3939         749 :   for (i = 3; i < l; i++) gel(v,i) = ginv(stoi(odd(i)? i-2: 2-i));
    3940         749 :   return v;
    3941             : }
    3942             : 
    3943             : /* S = (1+x)^zk log(1+x)^logj (mod x^(n+1)) */
    3944             : static GEN
    3945        1757 : xlog1x(long n, long zk, long logj, long *pteich)
    3946             : {
    3947        1757 :   GEN S = logj? RgXn_powu_i(log1x(n), logj, n+1): NULL;
    3948        1757 :   if (zk)
    3949             :   {
    3950        1183 :     GEN L = deg1pol_shallow(gen_1, gen_1, 0); /* x+1 */
    3951        1183 :     *pteich += zk;
    3952        1183 :     if (zk < 0) { L = RgXn_inv(L,n+1); zk = -zk; }
    3953        1183 :     if (zk != 1) L = RgXn_powu_i(L, zk, n+1);
    3954        1183 :     S = S? RgXn_mul(S, L, n+1): L;
    3955             :   }
    3956        1757 :   return S;
    3957             : }
    3958             : 
    3959             : /* oms from mspadicmoments; integrate teichmuller^i * S(x) [S = NULL: 1]*/
    3960             : static GEN
    3961        2303 : mspadicint(GEN oms, long teichi, GEN S)
    3962             : {
    3963        2303 :   pari_sp av = avma;
    3964        2303 :   long p = oms_get_p(oms), n = oms_get_n(oms), n0 = oms_get_n0(oms);
    3965        2303 :   GEN vT = gel(oms,1), alpha = gel(oms,2), gp = utoipos(p);
    3966        2303 :   long loss = S? Z_lval(Q_denom(S), p): 0;
    3967        2303 :   long nfinal = minss(n-loss, n0);
    3968        2303 :   long i, la, l = lg(vT);
    3969        2303 :   GEN res = cgetg(l, t_COL), teich = NULL;
    3970             : 
    3971        2303 :   if (S) S = RgX_to_RgC(S,lg(gmael(vT,1,1))-1);
    3972        2303 :   if (p == 2)
    3973             :   {
    3974         448 :     la = 3; /* corresponds to [1,-1] */
    3975         448 :     teichi &= 1;
    3976             :   }
    3977             :   else
    3978             :   {
    3979        1855 :     la = p; /* corresponds to [1,2,...,p-1] */
    3980        1855 :     teichi = smodss(teichi, p-1);
    3981        1855 :     if (teichi) teich = teichmullerinit(p, n);
    3982             :   }
    3983        5320 :   for (i=1; i<l; i++)
    3984             :   {
    3985        3017 :     pari_sp av2 = avma;
    3986        3017 :     GEN s = gen_0, T = gel(vT,i);
    3987             :     long ia;
    3988       13895 :     for (ia = 1; ia < la; ia++)
    3989             :     { /* Ta[j+1] correct mod p^n */
    3990       10878 :       GEN Ta = gel(T,ia), v = S? RgV_dotproduct(Ta, S): gel(Ta,1);
    3991       10878 :       if (teichi && ia != 1)
    3992             :       {
    3993        3843 :         if (p != 2)
    3994        3626 :           v = gmul(v, gel(teich, Fl_powu(ia,teichi,p)));
    3995             :         else
    3996         217 :           if (teichi) v = gneg(v);
    3997             :       }
    3998       10878 :       s = gadd(s, v);
    3999             :     }
    4000        3017 :     s = gadd(s, zeropadic(gp,nfinal));
    4001        3017 :     gel(res,i) = gerepileupto(av2, s);
    4002             :   }
    4003        2303 :   return gerepileupto(av, gmul(alpha, res));
    4004             : }
    4005             : /* integrate P = polynomial in log(x); vlog[j+1] = mspadicint(0,log(1+x)^j) */
    4006             : static GEN
    4007         539 : mspadicint_RgXlog(GEN P, GEN vlog)
    4008             : {
    4009         539 :   long i, d = degpol(P);
    4010         539 :   GEN s = gmul(gel(P,2), gel(vlog,1));
    4011         539 :   for (i = 1; i <= d; i++) s = gadd(s, gmul(gel(P,i+2), gel(vlog,i+1)));
    4012         539 :   return s;
    4013             : };
    4014             : 
    4015             : /* oms from mspadicmoments */
    4016             : GEN
    4017          98 : mspadicseries(GEN oms, long teichi)
    4018             : {
    4019          98 :   pari_sp av = avma;
    4020             :   GEN S, L, X, vlog, s, s2, u, logu, bin;
    4021             :   long j, p, m, n, step, stop;
    4022          98 :   checkoms(oms);
    4023          98 :   n = oms_get_n0(oms);
    4024          98 :   if (n < 1)
    4025             :   {
    4026           0 :     s = zeroser(0,0);
    4027           0 :     if (oms_is_supersingular(oms)) s = mkvec2(s,s);
    4028           0 :     return gerepilecopy(av, s);
    4029             :   }
    4030          98 :   p = oms_get_p(oms);
    4031          98 :   vlog = cgetg(n+1, t_VEC);
    4032          98 :   step = p == 2? 2: 1;
    4033          98 :   stop = 0;
    4034          98 :   S = NULL;
    4035          98 :   L = log1x(n);
    4036         644 :   for (j = 0; j < n; j++)
    4037             :   {
    4038         616 :     if (j) stop += step + u_lval(j,p); /* = step*j + v_p(j!) */
    4039         616 :     if (stop >= n) break;
    4040             :     /* S = log(1+x)^j */
    4041         546 :     gel(vlog,j+1) = mspadicint(oms,teichi,S);
    4042         546 :     S = S? RgXn_mul(S, L, n+1): L;
    4043             :   }
    4044          98 :   m = j;
    4045          98 :   u = utoipos(p == 2? 5: 1+p);
    4046          98 :   logu = glog(cvtop(u, utoipos(p), 4*m), 0);
    4047          98 :   X = gdiv(pol_x(0), logu);
    4048          98 :   s = cgetg(m+1, t_VEC);
    4049          98 :   s2 = oms_is_supersingular(oms)? cgetg(m+1, t_VEC): NULL;
    4050          98 :   bin = pol_1(0);
    4051         539 :   for (j = 0; j < m; j++)
    4052             :   { /* bin = binomial(x/log(1+p+O(p^(4*n))), j) mod x^m */
    4053         539 :     GEN a, v = mspadicint_RgXlog(bin, vlog);
    4054         539 :     int done = 1;
    4055         539 :     gel(s,j+1) = a = gel(v,1);
    4056         539 :     if (!gequal0(a) || valp(a) > 0) done = 0; else setlg(s,j+1);
    4057         539 :     if (s2)
    4058             :     {
    4059         119 :       gel(s2,j+1) = a = gel(v,2);
    4060         119 :       if (!gequal0(a) || valp(a) > 0) done = 0; else setlg(s2,j+1);
    4061             :     }
    4062         539 :     if (done || j == m-1) break;
    4063         441 :     bin = RgXn_mul(bin, gdivgs(gsubgs(X, j), j+1), m);
    4064             :   }
    4065          98 :   s = gtoser(s,0,lg(s)-1);
    4066          98 :   if (s2) { s2 = gtoser(s2,0,lg(s2)-1); s = mkvec2(s, s2); }
    4067          98 :   if (kross(oms_get_D(oms), p) >= 0) return gerepilecopy(av, s);
    4068           7 :   return gerepileupto(av, gneg(s));
    4069             : }
    4070             : static void
    4071        1820 : parse_chi(GEN s, GEN *s1, GEN *s2)
    4072             : {
    4073        1820 :   if (!s) *s1 = *s2 = gen_0;
    4074        1687 :   else switch(typ(s))
    4075             :   {
    4076        1183 :     case t_INT: *s1 = *s2 = s; break;
    4077             :     case t_VEC:
    4078         504 :       if (lg(s) == 3)
    4079             :       {
    4080         504 :         *s1 = gel(s,1);
    4081         504 :         *s2 = gel(s,2);
    4082         504 :         if (typ(*s1) == t_INT && typ(*s2) == t_INT) break;
    4083             :       }
    4084           0 :     default: pari_err_TYPE("mspadicL",s);
    4085           0 :              *s1 = *s2 = NULL;
    4086             :   }
    4087        1820 : }
    4088             : /* oms from mspadicmoments
    4089             :  * r-th derivative of L(f,chi^s,psi) in direction <chi>
    4090             :    - s \in Z_p \times \Z/(p-1)\Z, s-> chi^s=<\chi>^s_1 omega^s_2)
    4091             :    - Z -> Z_p \times \Z/(p-1)\Z par s-> (s, s mod p-1).
    4092             :  */
    4093             : GEN
    4094        1757 : mspadicL(GEN oms, GEN s, long r)
    4095             : {
    4096        1757 :   pari_sp av = avma;
    4097             :   GEN s1, s2, z, S;
    4098             :   long p, n, teich;
    4099        1757 :   checkoms(oms);
    4100        1757 :   p = oms_get_p(oms);
    4101        1757 :   n = oms_get_n(oms);
    4102        1757 :   parse_chi(s, &s1,&s2);
    4103        1757 :   teich = umodiu(subii(s2,s1), p==2? 2: p-1);
    4104        1757 :   S = xlog1x(n, itos(s1), r, &teich);
    4105        1757 :   z = mspadicint(oms, teich, S);
    4106        1757 :   if (lg(z) == 2) z = gel(z,1);
    4107        1757 :   if (kross(oms_get_D(oms), p) < 0) z = gneg(z);
    4108        1757 :   return gerepilecopy(av, z);
    4109             : }
    4110             : 
    4111             : GEN
    4112          63 : ellpadicL(GEN E, GEN pp, long n, GEN s, long r, GEN DD)
    4113             : {
    4114          63 :   pari_sp av = avma;
    4115             :   GEN L, W, Wp, xpm, NE, s1,s2, oms, den;
    4116             :   long sign, D;
    4117             :   ulong p;
    4118             : 
    4119          63 :   if (DD && !Z_isfundamental(DD))
    4120           0 :     pari_err_DOMAIN("ellpadicL", "isfundamental(D)", "=", gen_0, DD);
    4121          63 :   if (typ(pp) != t_INT) pari_err_TYPE("ellpadicL",pp);
    4122          63 :   if (cmpis(pp,2) < 0) pari_err_PRIME("ellpadicL",pp);
    4123          63 :   if (n <= 0) pari_err_DOMAIN("ellpadicL","precision","<=",gen_0,stoi(n));
    4124          63 :   if (r < 0) pari_err_DOMAIN("ellpadicL","r","<",gen_0,stoi(r));
    4125          63 :   parse_chi(s, &s1,&s2);
    4126          63 :   if (!DD) { sign = 1; D = 1; }
    4127             :   else
    4128             :   {
    4129           0 :     sign = signe(DD); D = itos(DD);
    4130           0 :     if (!sign) pari_err_DOMAIN("ellpadicL", "D", "=", gen_0, DD);
    4131             :   }
    4132          63 :   if (mpodd(s2)) sign = -sign;
    4133          63 :   W = msfromell(E, sign);
    4134          63 :   xpm = gel(W,2);
    4135          63 :   W = gel(W,1);
    4136             : 
    4137          63 :   p = itou(pp);
    4138          63 :   NE = ellQ_get_N(E);
    4139          63 :   if (dvdii(NE, sqri(pp))) pari_err_IMPL("additive reduction in ellpadicL");
    4140             : 
    4141          63 :   xpm = Q_remove_denom(xpm,&den);
    4142          63 :   if (!den) den = gen_1;
    4143          63 :   n += Z_lval(den, p);
    4144             : 
    4145          63 :   Wp = mspadicinit(W, p, n, umodiu(ellap(E,pp),p)? 0: 1);
    4146          63 :   oms = mspadicmoments(Wp, xpm, D);
    4147          63 :   L = mspadicL(oms, s, r);
    4148          63 :   return gerepileupto(av, gdiv(L,den));
    4149             : }
    4150             : 
    4151             : /* D coprime to p; Euler factor for the twisted L-function L(E,(D|.)),
    4152             :  * small difference with L(E_D) */
    4153             : static GEN
    4154           0 : ellpadicLeul(GEN E, GEN ED, GEN ND, GEN p, long n, GEN D, long r)
    4155             : {
    4156           0 :   GEN Z = ellpadicL(E, p, n, 0, r, D);
    4157           0 :   GEN F, U, apD = ellap(ED,p);
    4158           0 :   if (typ(Z) == t_COL)
    4159             :   { /* p | a_p(E_D), frobenius on E_D */
    4160           0 :     F = mkmat22(gen_0, negi(p), gen_1, apD);
    4161           0 :     U = RgM_RgC_mul(gpowgs(gsubsg(1, gdiv(F,p)), -2), Z);
    4162           0 :     settyp(U, t_VEC);
    4163             :   }
    4164             :   else
    4165             :   {
    4166           0 :     U = Z;
    4167           0 :     if (dvdii(ND,p)) /* assume a_p(E_D) = -1 */
    4168           0 :       U = gdivgs(U, 2);
    4169             :     else
    4170             :     {
    4171           0 :       GEN a = ms_unit_eigenvalue(apD, 2, p, n);
    4172           0 :       U = gmul(U, gpowgs(gsubsg(1, ginv(a)), -2));
    4173             :     }
    4174             :   }
    4175           0 :   return U;
    4176             : }
    4177             : 
    4178             : GEN
    4179           0 : ellpadicbsd(GEN E, GEN p, long n, GEN D)
    4180             : {
    4181           0 :   pari_sp av = avma;
    4182             :   GEN ED, tam, Lstar, N, C;
    4183             :   long r, vN;
    4184           0 :   checkell(E);
    4185           0 :   if (ell_get_type(E) != t_ELL_Q) pari_err_TYPE("ellpadicbsd",E);
    4186           0 :   if (D) {
    4187           0 :     if (typ(D) != t_INT) pari_err_TYPE("ellpadicbsd",D);
    4188           0 :     if (equali1(D)) D = NULL;
    4189             :   }
    4190           0 :   if (typ(p) != t_INT) pari_err_TYPE("ellpadicbsd",D);
    4191           0 :   if (n <= 0) pari_err_DOMAIN("ellpadicbsd","precision","<=",gen_0,stoi(n));
    4192           0 :   ED = D? ellinit(elltwist(E,D), gen_1, 0): E;
    4193           0 :   ED = ellanal_globalred_all(ED, NULL, &N, &tam);
    4194           0 :   r = itos( gel(ellanalyticrank_bitprec(ED, NULL, 32), 1) );
    4195             :   /* additive reduction ? */
    4196           0 :   vN = Z_pval(N, p);
    4197           0 :   if (vN >= 2) pari_err_DOMAIN("ellpadicbsd","v_p(N)", ">", gen_1, stoi(vN));
    4198           0 :   if (vN == 1 && equali1(ellap(ED,p)))
    4199           0 :     pari_err_IMPL("ellpadicbsd in the multiplicative reduction case");
    4200             :   /* TODO: should be something like
    4201             :    *   Lstar = ellpadicLeul(E,p,n,r+1,D)/(r+1)!/Linvariant(ED,p,n);
    4202             :    */
    4203           0 :   Lstar = ellpadicLeul(E, ED, N, p, n, D, r);
    4204           0 :   C = mulii(tam,mpfact(r));
    4205           0 :   if (D) C = gmul(C, get_Euler(ED, D));
    4206           0 :   C = gdiv(sqru(torsion_order(ED)), C);
    4207           0 :   if (D) obj_free(ED);
    4208           0 :   return gerepileupto(av, gmul(Lstar, C));
    4209             : }

Generated by: LCOV version 1.11