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

Generated by: LCOV version 1.11