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.8.0 lcov report (development 18953-2660f45) Lines: 1981 2051 96.6 %
Date: 2016-05-29 Functions: 210 211 99.5 %
Legend: Lines: hit not hit | Branches: + taken - not taken # not executed Branches: 853 1020 83.6 %

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

Generated by: LCOV version 1.9