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-bordeaux1.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 16827-3d78da8) Lines: 1128 1214 92.9 %
Date: 2014-09-29 Functions: 132 137 96.4 %
Legend: Lines: hit not hit | Branches: + taken - not taken # not executed Branches: 424 511 83.0 %

           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 modsymbkinit(ulong N, long k, long sign);
      22                 :            : static GEN new_subspace_trivial(GEN W);
      23                 :            : static GEN Cuspidal_subspace_trivial(GEN W0);
      24                 :            : static GEN ZGl2Q_star(GEN v);
      25                 :            : static GEN getMorphism_trivial(GEN WW1, GEN WW2, GEN v);
      26                 :            : static GEN getMorphism(GEN W1, GEN W2, GEN v);
      27                 :            : static GEN voo_act_Gl2Q(GEN g, long k);
      28                 :            : 
      29                 :            : /* Input: P^1(Z/NZ) (formed by create_p1mod)
      30                 :            :    Output: # P^1(Z/NZ) */
      31                 :            : static long
      32                 :        644 : p1_size(GEN p1N) { return lg(gel(p1N,1)) - 1; }
      33                 :            : static ulong
      34                 :   22770811 : p1N_get_N(GEN p1N) { return gel(p1N,3)[2]; }
      35                 :            : static GEN
      36                 :   20788446 : p1N_get_hash(GEN p1N) { return gel(p1N,2); }
      37                 :            : static GEN
      38                 :        203 : p1N_get_fa(GEN p1N) { return gel(p1N,4); }
      39                 :            : static GEN
      40                 :        161 : p1N_get_div(GEN p1N) { return gel(p1N,5); }
      41                 :            : /* modsymb-specific accessors */
      42                 :            : /* W a modsymb or an ellsym */
      43                 :            : static GEN
      44         [ +  + ]:    1969520 : get_modsymb(GEN W) { return lg(W) == 4? gel(W,1): W; }
      45                 :            : static GEN
      46                 :       1337 : modsymb_get_p1N(GEN W) { W = get_modsymb(W); return gel(W,1); }
      47                 :            : static long
      48                 :        329 : modsymb_get_N(GEN W) { return p1N_get_N(modsymb_get_p1N(W)); }
      49                 :            : static GEN
      50                 :         35 : modsymb_get_hashcusps(GEN W) { W = get_modsymb(W); return gel(W,16); }
      51                 :            : static GEN
      52                 :      29953 : modsymb_get_section(GEN W) { W = get_modsymb(W); return gel(W,12); }
      53                 :            : static GEN
      54                 :        462 : modsymb_get_genindex(GEN W) { W = get_modsymb(W); return gel(W,5); }
      55                 :            : static long
      56                 :        294 : modsymb_get_nbgen(GEN W) { return lg(modsymb_get_genindex(W))-1; }
      57                 :            : static long
      58                 :    1936662 : modsymb_get_nbE1(GEN W)
      59                 :            : {
      60                 :            :   GEN W11;
      61                 :    1936662 :   W = get_modsymb(W); W11 = gel(W,11);
      62                 :    1936662 :   return W11[4] - W11[3];
      63                 :            : }
      64                 :            : /* modymbk-specific accessors */
      65                 :            : #if 0
      66                 :            : static long
      67                 :            : modsymbk_get_dim(GEN W) { return gmael(W,3,2)[2]; }
      68                 :            : #endif
      69                 :            : static GEN
      70                 :        336 : modsymbk_get_basis(GEN W) { return gmael(W,3,1); }
      71                 :            : static long
      72                 :        581 : modsymbk_get_weight(GEN W) { return gmael(W,3,2)[1]; }
      73                 :            : static GEN
      74                 :        266 : modsymbk_get_st(GEN W) { return gmael(W,3,3); }
      75                 :            : static GEN
      76                 :        266 : modsymbk_get_link(GEN W) { return gmael(W,3,4); }
      77                 :            : static GEN
      78                 :        266 : modsymbk_get_invphiblock(GEN W) { return gmael(W,3,5); }
      79                 :            : static long
      80                 :         35 : modsymbk_get_sign(GEN W) { return itos(gmael(W,2,1)); }
      81                 :            : static GEN
      82                 :         84 : modsymbk_get_Sigma(GEN W) { return gmael(W,2,2); }
      83                 :            : /* ellsym-specific accessors */
      84                 :            : static GEN
      85                 :    1916656 : ellsym_get_W(GEN E) { GEN W = gel(E,1); return gel(W,1); }
      86                 :            : static GEN
      87                 :    1916656 : ellsym_get_dualell(GEN E) { return gel(E,2); }
      88                 :            : static GEN
      89                 :         63 : ellsym_get_ell(GEN E) { return gel(E,3); }
      90                 :            : 
      91                 :            : /** MODULAR TO SYM **/
      92                 :            : 
      93                 :            : /* q a t_FRAC or t_INT */
      94                 :            : static GEN
      95                 :    1981133 : modulartosym_init(ulong N, GEN q)
      96                 :            : {
      97                 :            :   long l, n;
      98                 :            :   GEN Q;
      99                 :            : 
     100                 :    1981133 :   q = gboundcf(q, 0);
     101                 :    1981133 :   l = lg(q);
     102                 :    1981133 :   Q = cgetg(l, t_VECSMALL);
     103                 :    1981133 :   Q[1] = 1;
     104         [ +  + ]:   20644554 :   for (n=2; n <l; n++) Q[n] = umodiu(gel(q,n), N);
     105         [ +  + ]:   18664730 :   for (n=3; n < l; n++)
     106                 :   16683597 :     Q[n] = Fl_add(Fl_mul(Q[n], Q[n-1], N), Q[n-2], N);
     107                 :    1981133 :   return Q;
     108                 :            : }
     109                 :            : 
     110                 :            : /** INIT MODSYM STRUCTURE, WEIGHT 2 **/
     111                 :            : 
     112                 :            : /* num = [Gamma : Gamma_0(N)] = N * Prod_{p|N} (1+p^-1) */
     113                 :            : static ulong
     114                 :        161 : count_Manin_symbols(ulong N, GEN P)
     115                 :            : {
     116                 :        161 :   long i, l = lg(P);
     117                 :        161 :   ulong num = N;
     118         [ +  + ]:        371 :   for (i = 1; i < l; i++) { ulong p = P[i]; num *= p+1; num /= p; }
     119                 :        161 :   return num;
     120                 :            : }
     121                 :            : /* returns the list of "Manin symbols" (c,d) in (Z/NZ)^2, (c,d,N) = 1
     122                 :            :  * generating H^1(X_0(N), Z) */
     123                 :            : static GEN
     124                 :        161 : generatemsymbols(ulong N, ulong num, GEN divN)
     125                 :            : {
     126                 :        161 :   GEN ret = cgetg(num+1, t_VEC);
     127                 :        161 :   ulong c, d, curn = 0;
     128                 :            :   long i, l;
     129                 :            :   /* generate Manin-symbols in two lists: */
     130                 :            :   /* list 1: (c:1) for 0 <= c < N */
     131         [ +  + ]:      29547 :   for (c = 0; c < N; c++) gel(ret, ++curn) = mkvecsmall2(c, 1);
     132         [ -  + ]:        161 :   if (N == 1) return ret;
     133                 :            :   /* list 2: (c:d) with 1 <= c < N, c | N, 0 <= d < N, gcd(d,N) > 1, gcd(c,d)=1.
     134                 :            :    * Furthermore, d != d0 (mod N/c) with c,d0 already in the list */
     135                 :        161 :   l = lg(divN) - 1;
     136                 :            :   /* c = 1 first */
     137                 :        161 :   gel(ret, ++curn) = mkvecsmall2(1,0);
     138         [ +  + ]:      29225 :   for (d = 2; d < N; d++)
     139         [ +  + ]:      29064 :     if (ugcd(d,N) != 1UL)
     140                 :       8575 :       gel(ret, ++curn) = mkvecsmall2(1,d);
     141                 :            :   /* omit c = 1 (first) and c = N (last) */
     142         [ +  + ]:        504 :   for (i=2; i < l; i++)
     143                 :            :   {
     144                 :            :     ulong Novc, d0;
     145                 :        343 :     c = divN[i];
     146                 :        343 :     Novc = N / c;
     147         [ +  + ]:      17311 :     for (d0 = 2; d0 <= Novc; d0++)
     148                 :            :     {
     149                 :      16968 :       ulong k, d = d0;
     150         [ +  + ]:      16968 :       if (ugcd(d, Novc) == 1UL) continue;
     151         [ +  + ]:      62468 :       for (k = 0; k < c; k++, d += Novc)
     152         [ +  + ]:      56224 :         if (ugcd(c,d) == 1UL)
     153                 :            :         {
     154                 :       2492 :           gel(ret, ++curn) = mkvecsmall2(c,d);
     155                 :       2492 :           break;
     156                 :            :         }
     157                 :            :     }
     158                 :            :   }
     159         [ -  + ]:        161 :   if (curn != num) pari_err_BUG("generatemsymbols [wrong number of symbols]");
     160                 :        161 :   return ret;
     161                 :            : }
     162                 :            : 
     163                 :            : #if OLD_HASH
     164                 :            : static ulong
     165                 :            : hash2(GEN H, long N, long a, long b)
     166                 :            : { return ucoeff(H, smodss(a,N) + 1, smodss(b,N) + 1); }
     167                 :            : /* symbols from generatemsymbols(). Returns H such that
     168                 :            :  * H[ nc mod N, nd mod N ] = index of (c,d) in 'symbols', n < N, (n,N) = 1 */
     169                 :            : static GEN
     170                 :            : inithashmsymbols(ulong N, GEN symbols)
     171                 :            : {
     172                 :            :   GEN H = zero_Flm_copy(N, N);
     173                 :            :   long k, l = lg(symbols);
     174                 :            :   ulong n;
     175                 :            :   for (n = 1; n < N; n++)
     176                 :            :     if (ugcd(n,N) == 1)
     177                 :            :       for (k=1; k < l; k++)
     178                 :            :       {
     179                 :            :         GEN s = gel(symbols, k);
     180                 :            :         ucoeff(H, Fl_mul(s[1],n,N) + 1, Fl_mul(s[2],n,N) + 1) = k;
     181                 :            :       }
     182                 :            :   return H;
     183                 :            : }
     184                 :            : #else
     185                 :            : static GEN
     186                 :        161 : inithashmsymbols(ulong N, GEN symbols)
     187                 :            : {
     188                 :        161 :   GEN H = zerovec(N);
     189                 :        161 :   long k, l = lg(symbols);
     190                 :            :   /* skip the (c:1), 0 <= c < N and (1:0) */
     191         [ +  + ]:      11228 :   for (k=N+2; k < l; k++)
     192                 :            :   {
     193                 :      11067 :     GEN s = gel(symbols, k);
     194                 :      11067 :     ulong c = s[1], d = s[2], Novc = N/c;
     195         [ +  + ]:      11067 :     if (gel(H,c) == gen_0) gel(H,c) = const_vecsmall(Novc+1,0);
     196 [ +  + ][ +  + ]:      11067 :     if (c != 1) { d %= Novc; if (!d) d = Novc; }
     197                 :      11067 :     mael(H, c, d) = k;
     198                 :            :   }
     199                 :        161 :   return H;
     200                 :            : }
     201                 :            : #endif
     202                 :            : 
     203                 :            : /** Helper functions for Sl2(Z) / Gamma_0(N) **/
     204                 :            : /* M a 2x2 ZM in SL2(Z) */
     205                 :            : static GEN
     206                 :       7616 : SL2_inv(GEN M)
     207                 :            : {
     208                 :       7616 :   GEN a=gcoeff(M,1,1), b=gcoeff(M,1,2), c=gcoeff(M,2,1), d=gcoeff(M,2,2);
     209                 :       7616 :   return mkmat2(mkcol2(d, negi(c)), mkcol2(negi(b), a));
     210                 :            : }
     211                 :            : /* M a 2x2 zm in SL2(Z) */
     212                 :            : static GEN
     213                 :      28448 : sl2_inv(GEN M)
     214                 :            : {
     215                 :      28448 :   long a=coeff(M,1,1), b=coeff(M,1,2), c=coeff(M,2,1), d=coeff(M,2,2);
     216                 :      28448 :   return mkmat2(mkvecsmall2(d, -c), mkvecsmall2(-b, a));
     217                 :            : }
     218                 :            : /* Return the zm [a,b; c,d] */
     219                 :            : static GEN
     220                 :     273840 : mat2(long a, long b, long c, long d)
     221                 :     273840 : { return mkmat2(mkvecsmall2(a,c), mkvecsmall2(b,d)); }
     222                 :            : 
     223                 :            : /* Input: a = 2-vector = path = {r/s,x/y}
     224                 :            :  * Output: either [r,x;s,y] or [-r,x;-s,y], whichever has determinant > 0 */
     225                 :            : static GEN
     226                 :     156814 : path_to_matrix(GEN a)
     227                 :            : {
     228                 :     156814 :   GEN v = gel(a,1), w = gel(a,2);
     229                 :     156814 :   long r = v[1], s = v[2], x = w[1], y = w[2];
     230         [ +  + ]:     156814 :   if (cmpii(mulss(r,y), mulss(x,s)) < 0) { r = -r; s = -s; }
     231                 :     156814 :   return mat2(r,x,s,y);
     232                 :            : }
     233                 :            : /* path from c1 to c2 */
     234                 :            : static GEN
     235                 :     115220 : mkpath(GEN c1, GEN c2) { return mat2(c1[1], c2[1], c1[2], c2[2]); }
     236                 :            : static long
     237                 :     116004 : cc(GEN M) { GEN v = gel(M,1); return v[2]; }
     238                 :            : static long
     239                 :     116004 : dd(GEN M) { GEN v = gel(M,2); return v[2]; }
     240                 :            : 
     241                 :            : /*Input: a,b = 2 paths, N = integer
     242                 :            :  *Output: 1 if the a,b are \Gamma_0(N)-equivalent; 0 otherwise */
     243                 :            : static int
     244                 :      13510 : gamma_equiv(GEN a, GEN b, ulong N)
     245                 :            : {
     246                 :      13510 :   pari_sp av = avma;
     247                 :      13510 :   GEN m = path_to_matrix(a);
     248                 :      13510 :   GEN n = path_to_matrix(b);
     249                 :      13510 :   GEN d = subii(mulss(cc(m),dd(n)), mulss(dd(m),cc(n)));
     250                 :      13510 :   ulong res = umodiu(d, N);
     251                 :      13510 :   avma = av; return res == 0;
     252                 :            : }
     253                 :            : /* Input: a,b = 2 paths that are \Gamma_0(N)-equivalent, N = integer
     254                 :            :  * Output: M in \Gamma_0(N) such that Mb=a */
     255                 :            : static GEN
     256                 :       6944 : gamma_equiv_matrix(GEN a, GEN b)
     257                 :            : {
     258                 :       6944 :   GEN m = zm_to_ZM( path_to_matrix(a) );
     259                 :       6944 :   GEN n = zm_to_ZM( path_to_matrix(b) );
     260                 :       6944 :   return ZM_mul(m, SL2_inv(n));
     261                 :            : }
     262                 :            : 
     263                 :            : /*************/
     264                 :            : /* P^1(Z/NZ) */
     265                 :            : /*************/
     266                 :            : 
     267                 :            : /* Input: N = integer
     268                 :            :  * Output: creates P^1(Z/NZ) = [symbols, H, N]
     269                 :            :  *   symbols: list of vectors [x,y] that give a set of representatives
     270                 :            :  *            of P^1(Z/NZ)
     271                 :            :  *   H: an M by M grid whose value at the r,c-th place is the index of the
     272                 :            :  *      "standard representative" equivalent to [r,c] occuring in the first
     273                 :            :  *      list. If gcd(r,c,N) > 1 the grid has value 0. */
     274                 :            : static GEN
     275                 :        161 : create_p1mod(ulong N)
     276                 :            : {
     277                 :        161 :   GEN fa = factoru(N), div = divisorsu(N);
     278                 :        161 :   ulong nsym = count_Manin_symbols(N, gel(fa,1));
     279                 :        161 :   GEN symbols = generatemsymbols(N, nsym, div);
     280                 :        161 :   GEN H = inithashmsymbols(N,symbols);
     281                 :        161 :   return mkvec5(symbols, H, utoipos(N), fa, div);
     282                 :            : }
     283                 :            : 
     284                 :            : /* result is known to be representable as an ulong */
     285                 :            : static ulong
     286                 :     410942 : lcmuu(ulong a, ulong b) { ulong d = ugcd(a,b); return (a/d) * b; }
     287                 :            : /* a != 0 in Z/NZ. Return u in (Z/NZ)^* such that au = gcd(a, N) (mod N)*/
     288                 :            : static ulong
     289                 :    5624668 : Fl_inverse(ulong a, ulong N)
     290                 :            : {
     291                 :            :   pari_sp av;
     292                 :    5624668 :   ulong d, d0, d1, e, u = Fl_invgen(a, N, &d);
     293         [ +  + ]:    5624668 :   if (d == 1) return u;
     294                 :     461895 :   e = N/d;
     295                 :     461895 :   d0 = ucoprime_part(d, e); /* d = d0 d1, d0 coprime to N/d, core(d1) | N/d */
     296         [ +  + ]:     461895 :   if (d0 == 1) return u;
     297                 :     410942 :   av = avma;
     298                 :     410942 :   d1 = d / d0;
     299                 :     410942 :   e = lcmuu(e, d1);
     300                 :     410942 :   u = itou(Z_chinese_coprime(utoipos(u), gen_1,
     301                 :            :                              utoipos(e), utoipos(d0), utoipos(e*d0)));
     302                 :    5624668 :   avma = av; return u;
     303                 :            : }
     304                 :            : /* Let (c : d) in P1(Z/NZ).
     305                 :            :  * If c = 0 return (0:1). If d = 0 return (1:0).
     306                 :            :  * Else replace by (cu : du), where u in (Z/NZ)^* such that C := cu = gcd(c,N).
     307                 :            :  * In create_p1mod(), (c : d) is represented by (C:D) where D = du (mod N/c)
     308                 :            :  * is smallest such that gcd(C,D) = 1. Return (C : du mod N/c), which need
     309                 :            :  * not belong to P1(Z/NZ) ! A second component du mod N/c = 0 is replaced by
     310                 :            :  * N/c in this case to avoid problems with array indices */
     311                 :            : static GEN
     312                 :   20788446 : p1_std_form(long c, long d, ulong N)
     313                 :            : {
     314                 :            :   ulong u;
     315                 :   20788446 :   c = smodss(c, N);
     316                 :   20788446 :   d = smodss(d, N);
     317         [ +  + ]:   20788446 :   if (!c) return mkvecsmall2(0, 1);
     318         [ +  + ]:   18521118 :   if (!d) return mkvecsmall2(1, 0);
     319                 :   18233726 :   u = Fl_invsafe(d, N);
     320         [ +  + ]:   18233726 :   if (u != 0) return mkvecsmall2(Fl_mul(c,u,N), 1); /* (d,N) = 1 */
     321                 :            : 
     322                 :    5619789 :   u = Fl_inverse(c, N);
     323         [ +  + ]:    5619789 :   if (u > 1) { c = Fl_mul(c,u,N); d = Fl_mul(d,u,N); }
     324                 :            :   /* c | N */
     325         [ +  + ]:    5619789 :   if (c != 1) d = d % (N/c);
     326         [ +  + ]:    5619789 :   if (!d) d = N/c;
     327                 :   20788446 :   return mkvecsmall2(c, d);
     328                 :            : }
     329                 :            : 
     330                 :            : /* Input: v = [x,y] = elt of P^1(Z/NZ) = class in Gamma_0(N) \ PSL2(Z)
     331                 :            :  * Output: returns the index of the standard rep equivalent to v */
     332                 :            : static long
     333                 :   20788446 : p1_index(long x, long y, GEN p1N)
     334                 :            : {
     335                 :   20788446 :   ulong N = p1N_get_N(p1N);
     336                 :   20788446 :   GEN H = p1N_get_hash(p1N), c;
     337                 :            : 
     338                 :            : #ifdef OLD_HASH
     339                 :            :   return hash2(p1N_get_hash(p1N), N, x, y);
     340                 :            : #else
     341                 :   20788446 :   c = p1_std_form(x, y, N);
     342                 :   20788446 :   x = c[1];
     343                 :   20788446 :   y = c[2];
     344         [ +  + ]:   20788446 :   if (y == 1) return x+1;
     345         [ +  + ]:    5907181 :   if (y == 0) return N+1;
     346         [ -  + ]:    5619789 :   if (mael(H,x,y) == 0) pari_err_BUG("p1_index");
     347                 :   20788446 :   return mael(H,x,y);
     348                 :            : #endif
     349                 :            : }
     350                 :            : 
     351                 :            : /* Cusps for \Gamma_0(N) */
     352                 :            : 
     353                 :            : /* \sum_{d | N} \phi(gcd(d, N/d)), using multiplicativity. fa = factor(N) */
     354                 :            : static ulong
     355                 :        161 : nbcusp(GEN fa)
     356                 :            : {
     357                 :        161 :   GEN P = gel(fa,1), E = gel(fa,2);
     358                 :        161 :   long i, l = lg(P);
     359                 :        161 :   ulong T = 1;
     360         [ +  + ]:        371 :   for (i = 1; i < l; i++)
     361                 :            :   {
     362                 :        210 :     long e = E[i] >> 1; /* floor(E[i] / 2) */
     363                 :        210 :     ulong p = P[i];
     364         [ +  + ]:        210 :     if (odd(E[i]))
     365                 :        182 :       T *= 2 * upowuu(p, e);
     366                 :            :     else
     367                 :         28 :       T *= (p+1) * upowuu(p, e - 1);
     368                 :            :   }
     369                 :        161 :   return T;
     370                 :            : }
     371                 :            : 
     372                 :            : /* to each cusp in \Gamma_0(N) P1(Q), represented by p/q, we associate a
     373                 :            :  * unique index. Canonical representative: (1:0) or (p:q) with q | N, q < N,
     374                 :            :  * p defined modulo d := gcd(N/q,q), (p,d) = 1.
     375                 :            :  * Return [[N, nbcusps], H, cusps]*/
     376                 :            : static GEN
     377                 :        161 : inithashcusps(GEN p1N)
     378                 :            : {
     379                 :        161 :   ulong N = p1N_get_N(p1N);
     380                 :        161 :   GEN div = p1N_get_div(p1N), H = zerovec(N+1);
     381                 :        161 :   long k, ind, l = lg(div), ncusp = nbcusp(p1N_get_fa(p1N));
     382                 :        161 :   GEN cusps = cgetg(ncusp+1, t_VEC);
     383                 :            : 
     384                 :        161 :   gel(H,1) = mkvecsmall2(0/*empty*/, 1/* first cusp: (1:0) */);
     385                 :        161 :   gel(cusps, 1) = mkvecsmall2(1,0);
     386                 :        161 :   ind = 2;
     387         [ +  + ]:        665 :   for (k=1; k < l-1; k++) /* l-1: remove q = N */
     388                 :            :   {
     389                 :        504 :     ulong p, q = div[k], d = ugcd(q, N/q);
     390                 :        504 :     GEN h = const_vecsmall(d+1,0);
     391                 :        504 :     gel(H,q+1) = h ;
     392         [ +  + ]:       1911 :     for (p = 0; p < d; p++)
     393         [ +  + ]:       1407 :       if (ugcd(p,d) == 1)
     394                 :            :       {
     395                 :        882 :         h[p+1] = ind;
     396                 :        882 :         gel(cusps, ind) = mkvecsmall2(p,q);
     397                 :        882 :         ind++;
     398                 :            :       }
     399                 :            :   }
     400                 :        161 :   return mkvec3(mkvecsmall2(N,ind-1), H, cusps);
     401                 :            : }
     402                 :            : /* c = [p,q], (p,q) = 1, return a canonical representative for
     403                 :            :  * \Gamma_0(N)(p/q) */
     404                 :            : static GEN
     405                 :       4914 : cusp_std_form(GEN c, GEN S)
     406                 :            : {
     407                 :       4914 :   long p, N = gel(S,1)[1], q = smodss(c[2], N);
     408                 :            :   ulong u, d;
     409         [ +  + ]:       4914 :   if (q == 0) return mkvecsmall2(1, 0);
     410                 :       4879 :   p = smodss(c[1], N);
     411                 :       4879 :   u = Fl_inverse(q, N);
     412                 :       4879 :   q = Fl_mul(q,u, N);
     413                 :       4879 :   d = ugcd(q, N/q);
     414                 :       4914 :   return mkvecsmall2(Fl_div(p % d,u % d, d), q);
     415                 :            : }
     416                 :            : /* c = [p,q], (p,q) = 1, return the index of the corresponding cusp.
     417                 :            :  * S from inithashcusps */
     418                 :            : static ulong
     419                 :       4914 : cusp_index(GEN c, GEN S)
     420                 :            : {
     421                 :            :   long p, q;
     422                 :       4914 :   GEN H = gel(S,2);
     423                 :       4914 :   c = cusp_std_form(c, S);
     424                 :       4914 :   p = c[1]; q = c[2];
     425         [ -  + ]:       4914 :   if (!mael(H,q+1,p+1)) pari_err_BUG("cusp_index");
     426                 :       4914 :   return mael(H,q+1,p+1);
     427                 :            : }
     428                 :            : 
     429                 :            : /* M a square invertible ZM, return a ZM iM such that iM M = M iM = d.Id */
     430                 :            : static GEN
     431                 :        147 : ZM_inv_denom(GEN M)
     432                 :            : {
     433                 :        147 :   GEN ciM = ZM_det(M);
     434                 :        147 :   GEN c, iM = Q_primitive_part(ZM_inv(M, ciM), &c);
     435         [ +  + ]:        147 :   if (c) ciM = diviiexact(ciM, c);
     436                 :        147 :   return mkvec2(iM, ciM);
     437                 :            : }
     438                 :            : /* return M^(-1) v, dinv = ZM_inv_denom(M) OR Qevproj_init(M) */
     439                 :            : static GEN
     440                 :        483 : ZC_apply_dinv(GEN dinv, GEN v)
     441                 :            : {
     442                 :            :   GEN x, c, iM;
     443         [ +  + ]:        483 :   if (lg(dinv) == 3)
     444                 :            :   {
     445                 :        217 :     iM = gel(dinv,1);
     446                 :        217 :     c = gel(dinv,2);
     447                 :            :   }
     448                 :            :   else
     449                 :            :   { /* Qevproj_init */
     450                 :        266 :     iM = gel(dinv,2);
     451                 :        266 :     c = gel(dinv,3);
     452                 :        266 :     v = typ(v) == t_MAT? rowpermute(v, gel(dinv,4))
     453         [ -  + ]:        266 :                        : vecpermute(v, gel(dinv,4));
     454                 :            :   }
     455                 :        483 :   x = RgM_RgC_mul(iM, v);
     456         [ +  - ]:        483 :   if (!isint1(c)) x = RgC_Rg_div(x, c);
     457                 :        483 :   return x;
     458                 :            : }
     459                 :            : 
     460                 :            : /* M an n x d ZM of rank d (basis of a Q-subspace), n >= d.
     461                 :            :  * Initialize a projector on M */
     462                 :            : GEN
     463                 :         98 : Qevproj_init(GEN M)
     464                 :            : {
     465                 :         98 :   GEN v = ZM_indexrank(M), perm = gel(v,1);
     466                 :         98 :   GEN MM = rowpermute(M, perm); /* square invertible */
     467                 :         98 :   GEN dinv = ZM_inv_denom(MM), iM = gel(dinv,1), ciM = gel(dinv,2);
     468                 :         98 :   return mkvec4(M, iM, ciM, perm);
     469                 :            : }
     470                 :            : /* T an n x n QM, stabilizing d-dimensional Q-vector space spanned by the
     471                 :            :  * columns of M, pro = Qevproj_init(M). Return d x d matrix of T acting
     472                 :            :  * on M */
     473                 :            : GEN
     474                 :         42 : Qevproj_apply(GEN T, GEN pro)
     475                 :            : {
     476                 :         42 :   GEN M = gel(pro,1), iM = gel(pro,2), ciM = gel(pro,3), perm = gel(pro,4);
     477                 :         42 :   return RgM_Rg_div(RgM_mul(iM, RgM_mul(rowpermute(T,perm), M)), ciM);
     478                 :            : }
     479                 :            : /* Qevproj_apply(T,pro)[,k] */
     480                 :            : GEN
     481                 :         70 : Qevproj_apply_vecei(GEN T, GEN pro, long k)
     482                 :            : {
     483                 :         70 :   GEN M = gel(pro,1), iM = gel(pro,2), ciM = gel(pro,3), perm = gel(pro,4);
     484                 :         70 :   GEN v = RgM_RgC_mul(iM, RgM_RgC_mul(rowpermute(T,perm), gel(M,k)));
     485                 :         70 :   return RgC_Rg_div(v, ciM);
     486                 :            : }
     487                 :            : 
     488                 :            : /* normalize a Q-basis*/
     489                 :            : static GEN
     490                 :        133 : Q_primpart_basis(GEN M)
     491                 :            : {
     492                 :            :   long i, l;
     493                 :        133 :   GEN N = cgetg_copy(M, &l);
     494         [ +  + ]:       3472 :   for (i = 1; i < l; i++) gel(N,i) = Q_primpart(gel(M,i));
     495                 :        133 :   return N;
     496                 :            : }
     497                 :            : static GEN
     498                 :        112 : ZM_ker(GEN M) { return Q_primpart_basis(keri(M)); }
     499                 :            : static GEN
     500                 :         49 : QM_ker(GEN M) { return ZM_ker(Q_primpart(M)); }
     501                 :            : static GEN
     502                 :         21 : QM_image(GEN A)
     503                 :            : {
     504                 :         21 :   A = Q_primpart_basis(A);
     505                 :         21 :   return vecpermute(A, ZM_indeximage(A));
     506                 :            : }
     507                 :            : 
     508                 :            : /* Decompose the subspace H in simple subspaces.
     509                 :            :  * Eg for H = new_subspace */
     510                 :            : GEN
     511                 :         14 : simple_subspaces(GEN W, GEN H)
     512                 :            : {
     513                 :         14 :   ulong p, N = modsymb_get_N(W);
     514                 :            :   long dim, first;
     515                 :            :   forprime_t S;
     516                 :         14 :   GEN T1 = NULL, T2 = NULL, V;
     517         [ -  + ]:         14 :   if (lg(H) == 1) return H;
     518                 :         14 :   (void)u_forprime_init(&S, 2, ULONG_MAX);
     519                 :         14 :   dim = lg(H)-1;
     520                 :         14 :   V = vectrunc_init(dim);
     521                 :         14 :   vectrunc_append(V, Qevproj_init(H));
     522                 :         14 :   first = 1; /* V[1..first-1] contains simple subspaces */
     523         [ +  - ]:         21 :   while ((p = u_forprime_next(&S)))
     524                 :            :   {
     525                 :            :     GEN T;
     526                 :            :     long n, j, lV;
     527         [ -  + ]:         21 :     if (N % p == 0) continue;
     528 [ +  + ][ -  + ]:         21 :     if (T1 && T2)
     529                 :            :     {
     530                 :          0 :       T = RgM_add(T1,T2);
     531                 :          0 :       T2 = NULL;
     532                 :            :     }
     533                 :            :     else
     534                 :            :     {
     535                 :         21 :       T2 = T1;
     536                 :         21 :       T1 = T = modsymbHecke(W, p);
     537                 :            :     }
     538                 :         21 :     lV = lg(V);
     539         [ +  + ]:         42 :     for (j = first; j < lV; j++)
     540                 :            :     {
     541                 :         21 :       pari_sp av = avma;
     542                 :         21 :       GEN Vj = gel(V,j), P = gel(Vj,1);
     543                 :         21 :       GEN TVj = Qevproj_apply(T, Vj); /* c T | V_j */
     544                 :         21 :       GEN ch = QM_charpoly_ZX(TVj), fa = ZX_factor(ch), F, E;
     545                 :            :       long k;
     546                 :         21 :       F = gel(fa, 1);
     547                 :         21 :       E = gel(fa, 2);
     548                 :         21 :       n = lg(F)-1;
     549         [ +  + ]:         21 :       if (n == 1)
     550                 :            :       {
     551         [ +  - ]:         14 :         if (isint1(gel(E,1)))
     552                 :            :         { /* simple subspace */
     553                 :         14 :           swap(gel(V,first), gel(V,j));
     554                 :         14 :           first++;
     555                 :            :         }
     556                 :            :         else
     557                 :          0 :           avma = av;
     558                 :            :       }
     559                 :            :       else
     560                 :            :       { /* can split Vj */
     561                 :            :         GEN pows;
     562                 :          7 :         long D = 1;
     563         [ +  + ]:         42 :         for (k = 1; k <= n; k++)
     564                 :            :         {
     565                 :         35 :           long d = degpol(gel(F,k));
     566         [ +  + ]:         35 :           if (d > D) D = d;
     567                 :            :         }
     568                 :            :         /* remove V[j] */
     569                 :          7 :         swap(gel(V,j), gel(V,lg(V)-1)); setlg(V, lg(V)-1);
     570                 :          7 :         pows = RgM_powers(TVj, minss((long)2*sqrt(D), D));
     571         [ +  + ]:         42 :         for (k = 1; k <= n; k++)
     572                 :            :         {
     573                 :         35 :           GEN f = gel(F,k);
     574                 :         35 :           GEN M = RgX_RgMV_eval(f, pows); /* f(TVj) */
     575                 :         35 :           GEN K = QM_ker(M);
     576                 :         35 :           GEN p = Q_primpart( RgM_mul(P, K) );
     577                 :         35 :           vectrunc_append(V, Qevproj_init(p));
     578 [ +  - ][ +  + ]:         35 :           if (lg(K) == 2 || isint1(gel(E,k)))
     579                 :            :           { /* simple subspace */
     580                 :         28 :             swap(gel(V,first), gel(V, lg(V)-1));
     581                 :         28 :             first++;
     582                 :            :           }
     583                 :            :         }
     584         [ +  - ]:          7 :         if (j < first) j = first;
     585                 :            :       }
     586                 :            :     }
     587         [ +  + ]:         21 :     if (first >= lg(V)) return V;
     588                 :            :   }
     589                 :          0 :   pari_err_BUG("subspaces not found");
     590                 :         14 :   return NULL;
     591                 :            : }
     592                 :            : 
     593                 :            : /* proV = Qevproj_init of a Hecke simple subspace, return [ a_p, p <= B ] */
     594                 :            : GEN
     595                 :          7 : qexpansion(GEN W, GEN proV, long B)
     596                 :            : {
     597                 :          7 :   ulong p, N = modsymb_get_N(W);
     598                 :            :   long i, d;
     599                 :            :   forprime_t S;
     600                 :          7 :   GEN T1 = NULL, T2 = NULL, TV = NULL, ch, v, dTiv, Tiv, dinv, ciM, iM;
     601                 :          7 :   GEN AP = cgetg(uprimepi(B)+1, t_VEC);
     602         [ -  + ]:          7 :   if (B < 2) return AP;
     603                 :          7 :   (void)u_forprime_init(&S, 2, ULONG_MAX);
     604         [ +  - ]:          7 :   while ((p = u_forprime_next(&S)))
     605                 :            :   {
     606                 :            :     GEN T;
     607         [ -  + ]:          7 :     if (N % p == 0) continue;
     608 [ -  + ][ #  # ]:          7 :     if (T1 && T2)
     609                 :            :     {
     610                 :          0 :       T = RgM_add(T1,T2);
     611                 :          0 :       T2 = NULL;
     612                 :            :     }
     613                 :            :     else
     614                 :            :     {
     615                 :          7 :       T2 = T1;
     616                 :          7 :       T1 = T = modsymbHecke(W, p);
     617                 :            :     }
     618                 :          7 :     TV = Qevproj_apply(T, proV); /* T | V */
     619                 :          7 :     ch = QM_charpoly_ZX(TV);
     620         [ +  - ]:          7 :     if (ZX_is_irred(ch)) break;
     621                 :          0 :     ch = NULL;
     622                 :            :   }
     623         [ -  + ]:          7 :   if (!ch) pari_err_BUG("q-Expansion not found");
     624                 :          7 :   d = degpol(ch);
     625                 :          7 :   v = vec_ei(d, 1); /* take v = e_1 */
     626                 :          7 :   Tiv = cgetg(d+1, t_MAT); /* Tiv[i] = T^(i-1)v */
     627                 :          7 :   gel(Tiv, 1) = v;
     628         [ +  + ]:         14 :   for (i = 2; i <= d; i++) gel(Tiv, i) = RgM_RgC_mul(TV, gel(Tiv,i-1));
     629                 :          7 :   Tiv = Q_remove_denom(Tiv, &dTiv);
     630                 :          7 :   dinv = ZM_inv_denom(Tiv);
     631                 :          7 :   iM = gel(dinv,1);
     632                 :          7 :   ciM = gel(dinv,2);
     633         [ -  + ]:          7 :   if (dTiv) ciM = gdiv(ciM, dTiv);
     634                 :          7 :   (void)u_forprime_init(&S, 2, B);
     635                 :          7 :   i = 1;
     636         [ +  + ]:         77 :   while ((p = u_forprime_next(&S)))
     637                 :            :   {
     638                 :         70 :     GEN T = modsymbHecke(W, p), u;
     639                 :         70 :     GEN Tv = Qevproj_apply_vecei(T, proV, 1); /* Tp.v */
     640                 :            :     /* Write Tp.v = \sum u_i T^i v */
     641                 :         70 :     u = RgC_Rg_div(RgM_RgC_mul(iM, Tv), ciM);
     642                 :         70 :     gel(AP, i++) = RgV_to_RgX(u, 0);
     643                 :            :   }
     644                 :          7 :   return mkvec2(ch, AP);
     645                 :            : }
     646                 :            : 
     647                 :            : static GEN
     648                 :         35 : Qevproj_Sigma(GEN W, GEN H)
     649                 :            : {
     650                 :         35 :   long s = modsymbk_get_sign(W);
     651         [ +  + ]:         35 :   if (s)
     652                 :            :   { /* project on +/- component */
     653                 :         21 :     GEN Sigma = modsymbk_get_Sigma(W);
     654                 :         21 :     GEN A = gmul(Sigma, H);
     655         [ +  - ]:         21 :     A = (s > 0)? gadd(A, H): gsub(A, H);
     656                 :            :     /* Im(Sigma + sign) = Ker(Sigma - sign) */
     657                 :         21 :     H = QM_image(A);
     658                 :            :   }
     659                 :         35 :   return H;
     660                 :            : }
     661                 :            : static GEN
     662                 :         21 : new_subspace_trivial(GEN W)
     663                 :            : {
     664                 :         21 :   GEN p1N = modsymb_get_p1N(W), P = gel(p1N_get_fa(p1N), 1);
     665                 :         21 :   ulong N = modsymb_get_N(W);
     666                 :         21 :   long i, nP = lg(P)-1;
     667                 :            :   GEN Snew, K, v;
     668                 :            : 
     669                 :         21 :   K = Cuspidal_subspace_trivial(W);
     670         [ +  + ]:         21 :   if (uisprime(N)) { Snew = K; goto END; }
     671                 :          7 :   v = cgetg(2*nP + 1, t_COL);
     672         [ +  + ]:         21 :   for (i = 1; i <= nP; i++)
     673                 :            :   {
     674                 :         14 :     ulong M = N / P[i];
     675                 :         14 :     GEN T1, Td, Wi = modsymbkinit(M, 2, 0);
     676                 :         14 :     T1 = getMorphism_trivial(W, Wi, mat2(1,0,0,1));
     677                 :         14 :     Td = getMorphism_trivial(W, Wi, mat2(P[i],0,0,1));
     678                 :         14 :     gel(v,2*i-1) = ZM_mul(T1, K); /* multiply by K = restrict to Ker delta */
     679                 :         14 :     gel(v,2*i) = ZM_mul(Td, K);
     680                 :            :   }
     681                 :          7 :   Snew = ZM_mul(K, ZM_ker(matconcat(v)));
     682                 :            : END:
     683                 :         21 :   return Qevproj_Sigma(W, Snew);
     684                 :            : }
     685                 :            : 
     686                 :            : GEN
     687                 :         21 : new_subspace(GEN W)
     688                 :            : {
     689                 :         21 :   GEN p1N = modsymb_get_p1N(W), P = gel(p1N_get_fa(p1N), 1);
     690                 :         21 :   ulong p, N = modsymb_get_N(W);
     691                 :         21 :   long i, nP = lg(P)-1, k = modsymbk_get_weight(W);
     692                 :            :   GEN S, Snew, Sold, pr_S, pr_Sold, v;
     693                 :            :   forprime_t F;
     694         [ +  - ]:         21 :   if (k == 2) return new_subspace_trivial(W);
     695                 :          0 :   S = gel(Cuspidal_subspace(W), 2);
     696         [ #  # ]:          0 :   if (uisprime(N)) { Snew = S; goto END; }
     697                 :          0 :   v = cgetg(2*nP + 1, t_VEC);
     698         [ #  # ]:          0 :   for (i = 1; i <= nP; i++)
     699                 :            :   {
     700                 :          0 :     GEN Wi = modsymbkinit(N/P[i], k, 0);
     701                 :          0 :     gel(v,2*i-1) = getMorphism(Wi, W, mat2(1,0,0,1));
     702                 :          0 :     gel(v,2*i)   = getMorphism(Wi, W, mat2(P[i],0,0,1));
     703                 :            :   }
     704                 :          0 :   v = QM_image(matconcat(v)); /* old forms */
     705                 :            :   /* restrict to cuspidal subspace */
     706                 :          0 :   Sold = Q_primpart_basis(intersect(S, v));
     707                 :          0 :   pr_S = Qevproj_init(S);
     708                 :          0 :   pr_Sold = Qevproj_init(Sold);
     709                 :          0 :   Snew = NULL;
     710                 :          0 :   (void)u_forprime_init(&F, 2, ULONG_MAX);
     711         [ #  # ]:          0 :   while ((p = u_forprime_next(&F)))
     712                 :            :   {
     713                 :          0 :     pari_sp av = avma;
     714                 :            :     GEN T, chS, chSold, chSnew;
     715         [ #  # ]:          0 :     if (N % p == 0) continue;
     716                 :          0 :     T = modsymbHecke(W, p);
     717                 :          0 :     chS = QM_charpoly_ZX(Qevproj_apply(T, pr_S));
     718                 :          0 :     chSold = QM_charpoly_ZX(Qevproj_apply(T, pr_Sold));
     719                 :          0 :     chSnew = RgX_div(chS, chSold); /* = char T | S^new */
     720         [ #  # ]:          0 :     if (!degpol( ZX_gcd(chSnew, chSold) ))
     721                 :            :     {
     722                 :          0 :       GEN M = RgX_RgM_eval(chSnew, T);
     723                 :          0 :       Snew = QM_ker(M);
     724                 :          0 :       break;
     725                 :            :     }
     726                 :          0 :     avma = av;
     727                 :            :   }
     728                 :            : END:
     729                 :         21 :   return Qevproj_Sigma(W, Snew);
     730                 :            : }
     731                 :            : 
     732                 :            : /* Solve the Manin relations for a congruence subgroup \Gamma by constructing
     733                 :            :  * a well-formed fundamental domain for the action of \Gamma on upper half
     734                 :            :  * space. See
     735                 :            :  * Pollack and Stevens, Overconvergent modular symbols and p-adic L-functions
     736                 :            :  * Annales scientifiques de l'ENS 44, fascicule 1 (2011), 1-42
     737                 :            :  * http://math.bu.edu/people/rpollack/Papers/Overconvergent_modular_symbols_and_padic_Lfunctions.pdf
     738                 :            :  *
     739                 :            :  * FIXME: Implemented for \Gamma = \Gamma_0(N) only. */
     740                 :            : 
     741                 :            : #if 0 /* Pollack-Stevens shift their paths so as to solve equations of the
     742                 :            :          form f(z+1) - f(z) = g. We don't (to avoid mistakes) so we will
     743                 :            :          have to solve eqs of the form f(z-1) - f(z) = g */
     744                 :            : /* c = a/b; as a t_VECSMALL [a,b]; return c-1 as a t_VECSMALL */
     745                 :            : static GEN
     746                 :            : Shift_left_cusp(GEN c) { long a=c[1], b=c[2]; return mkvecsmall2(a - b, b); }
     747                 :            : /* c = a/b; as a t_VECSMALL [a,b]; return c+1 as a t_VECSMALL */
     748                 :            : static GEN
     749                 :            : Shift_right_cusp(GEN c) { long a=c[1], b=c[2]; return mkvecsmall2(a + b, b); }
     750                 :            : /*Input: path = [r,s] (thought of as a geodesic between these points)
     751                 :            :  *Output: The path shifted by one to the left, i.e. [r-1,s-1] */
     752                 :            : static GEN
     753                 :            : Shift_left(GEN path)
     754                 :            : {
     755                 :            :   GEN r = gel(path,1), s = gel(path,2);
     756                 :            :   return mkvec2(Shift_left_cusp(r), Shift_left_cusp(s)); }
     757                 :            : /*Input: path = [r,s] (thought of as a geodesic between these points)
     758                 :            :  *Output: The path shifted by one to the right, i.e. [r+1,s+1] */
     759                 :            : GEN
     760                 :            : Shift_right(GEN path)
     761                 :            : {
     762                 :            :   GEN r = gel(path,1), s = gel(path,2);
     763                 :            :   return mkvec2(Shift_right_cusp(r), Shift_right_cusp(s)); }
     764                 :            : #endif
     765                 :            : 
     766                 :            : /* linked lists */
     767                 :            : typedef struct list_t { GEN data; struct list_t *next; } list_t;
     768                 :            : static list_t *
     769                 :      13692 : list_new(GEN x)
     770                 :            : {
     771                 :      13692 :   list_t *L = (list_t*)stack_malloc(sizeof(list_t));
     772                 :      13692 :   L->data = x;
     773                 :      13692 :   L->next = NULL; return L;
     774                 :            : }
     775                 :            : static void
     776                 :      13531 : list_insert(list_t *L, GEN x)
     777                 :            : {
     778                 :      13531 :   list_t *l = list_new(x);
     779                 :      13531 :   l->next = L->next;
     780                 :      13531 :   L->next = l;
     781                 :      13531 : }
     782                 :            : 
     783                 :            : /*Input: N > 1, p1N = P^1(Z/NZ)
     784                 :            :  *Output: a connected fundamental domain for the action of \Gamma_0(N) on
     785                 :            :  *  upper half space.  When \Gamma_0(N) is torsion free, the domain has the
     786                 :            :  *  property that all of its vertices are cusps.  When \Gamma_0(N) has
     787                 :            :  *  three-torsion, 2 extra triangles need to be added.
     788                 :            :  *
     789                 :            :  * The domain is constructed by beginning with the triangle with vertices 0,1
     790                 :            :  * and oo.  Each adjacent triangle is successively tested to see if it contains
     791                 :            :  * points not \Gamma_0(N) equivalent to some point in our region.  If a
     792                 :            :  * triangle contains new points, it is added to the region.  This process is
     793                 :            :  * continued until the region can no longer be extended (and still be a
     794                 :            :  * fundamental domain) by added an adjacent triangle.  The list of cusps
     795                 :            :  * between 0 and 1 are then returned
     796                 :            :  *
     797                 :            :  * Precisely, the function returns a list such that the elements of the list
     798                 :            :  * with odd index are the cusps in increasing order.  The even elements of the
     799                 :            :  * list are either an "x" or a "t".  A "t" represents that there is an element
     800                 :            :  * of order three such that its fixed point is in the triangle directly
     801                 :            :  * adjacent to the our region with vertices given by the cusp before and after
     802                 :            :  * the "t".  The "x" represents that this is not the case. */
     803                 :            : enum { type_X, type_DO /* ? */, type_T };
     804                 :            : static GEN
     805                 :        161 : form_list_of_cusps(ulong N, GEN p1N)
     806                 :            : {
     807                 :        161 :   pari_sp av = avma;
     808                 :        161 :   long i, position, nbC = 2;
     809                 :            :   GEN v, L;
     810                 :            :   list_t *C, *c;
     811                 :            :   /* Let t be the index of a class in PSL2(Z) / \Gamma in our fixed enumeration
     812                 :            :    * v[t] != 0 iff it is the class of z tau^r for z a previous alpha_i
     813                 :            :    * or beta_i.
     814                 :            :    * For \Gamma = \Gamma_0(N), the enumeration is given by p1_index.
     815                 :            :    * We write cl(gamma) = the class of gamma mod \Gamma */
     816                 :        161 :   v = const_vecsmall(p1_size(p1N), 0);
     817                 :        161 :   i = p1_index( 0, 1, p1N); v[i] = 1;
     818                 :        161 :   i = p1_index( 1,-1, p1N); v[i] = 2;
     819                 :        161 :   i = p1_index(-1, 0, p1N); v[i] = 3;
     820                 :            :   /* the value is unused [debugging]: what matters is whether it is != 0 */
     821                 :        161 :   position = 4;
     822                 :            :   /* at this point, Fund = R, v contains the classes of Id, tau, tau^2 */
     823                 :            : 
     824                 :        161 :   C  = list_new(mkvecsmall3(0,1, type_X));
     825                 :        161 :   list_insert(C, mkvecsmall3(1,1,type_DO));
     826                 :            :   /* C is a list of triples[a,b,t], where c = a/b is a cusp, and t is the type
     827                 :            :    * of the path between c and the PREVIOUS cusp in the list, coded as
     828                 :            :    *   type_DO = "?", type_X = "x", type_T = "t"
     829                 :            :    * Initially, C = [0/1,"?",1/1]; */
     830                 :            : 
     831                 :            :   /* loop through the current set of cusps C and check to see if more cusps
     832                 :            :    * should be added */
     833                 :            :   for (;;)
     834                 :            :   {
     835                 :        903 :     int done = 1;
     836         [ +  - ]:      73052 :     for (c = C; c; c = c->next)
     837                 :            :     {
     838                 :            :       GEN cusp1, cusp2, gam;
     839                 :            :       long pos, b1, b2, b;
     840                 :            : 
     841         [ +  + ]:      73052 :       if (!c->next) break;
     842                 :      72149 :       cusp1 = c->data; /* = a1/b1 */
     843                 :      72149 :       cusp2 = (c->next)->data; /* = a2/b2 */
     844         [ +  + ]:      72149 :       if (cusp2[3] != type_DO) continue;
     845                 :            : 
     846                 :            :       /* gam (oo -> 0) = (cusp2 -> cusp1), gam in PSL2(Z) */
     847                 :      26901 :       gam = path_to_matrix(mkpath(cusp2, cusp1)); /* = [a2,a1;b2,b1] */
     848                 :            :       /* we have normalized the cusp representation so that a1 b2 - a2 b1 = 1 */
     849                 :      26901 :       b1 = coeff(gam,2,1); b2 = coeff(gam,2,2);
     850                 :            :       /* gam.1  = (a1 + a2) / (b1 + b2) */
     851                 :      26901 :       b = b1 + b2;
     852                 :            :       /* Determine whether the adjacent triangle *below* (cusp1->cusp2)
     853                 :            :        * should be added */
     854                 :      26901 :       pos = p1_index(b1,b2, p1N); /* did we see cl(gam) before ? */
     855         [ +  + ]:      26901 :       if (v[pos])
     856                 :      13510 :         cusp2[3] = type_X; /* NO */
     857                 :            :       else
     858                 :            :       { /* YES */
     859                 :            :         ulong B1, B2;
     860                 :      13391 :         v[pos] = position;
     861                 :      13391 :         i = p1_index(-(b1+b2), b1, p1N); v[i] = position+1;
     862                 :      13391 :         i = p1_index(b2, -(b1+b2), p1N); v[i] = position+2;
     863                 :            :         /* add cl(gam), cl(gam*TAU), cl(gam*TAU^2) to v */
     864                 :      13391 :         position += 3;
     865                 :            :         /* gam tau gam^(-1) in \Gamma ? */
     866                 :      13391 :         B1 = smodss(b1, N);
     867                 :      13391 :         B2 = smodss(b2, N);
     868         [ +  + ]:      13391 :         if ((Fl_sqr(B2,N) + Fl_sqr(B1,N) + Fl_mul(B1,B2,N)) % N == 0)
     869                 :         21 :           cusp2[3] = type_T;
     870                 :            :         else
     871                 :            :         {
     872                 :      13370 :           long a1 = coeff(gam, 1,1), a2 = coeff(gam, 1,2);
     873                 :      13370 :           long a = a1 + a2; /* gcd(a,b) = 1 */
     874                 :      13370 :           list_insert(c, mkvecsmall3(a,b,type_DO));
     875                 :      13370 :           c = c->next;
     876                 :      13370 :           nbC++;
     877                 :      13370 :           done = 0;
     878                 :            :         }
     879                 :            :       }
     880                 :            :     }
     881         [ +  + ]:        903 :     if (done) break;
     882                 :        742 :   }
     883                 :        161 :   L = cgetg(nbC+1, t_VEC); i = 1;
     884         [ +  + ]:      13853 :   for (c = C; c; c = c->next) gel(L,i++) = c->data;
     885                 :        161 :   return gerepilecopy(av, L);
     886                 :            : }
     887                 :            : 
     888                 :            : /* M in PSL2(Z). Return index of M in P1^(Z/NZ) = Gamma0(N) \ PSL2(Z),
     889                 :            :  * and M0 in Gamma_0(N) such that M = M0 * M', where M' = chosen
     890                 :            :  * section( PSL2(Z) -> P1^(Z/NZ) ). */
     891                 :            : static GEN
     892                 :        742 : Gamma0N_decompose(GEN W, GEN M, long *index)
     893                 :            : {
     894                 :        742 :   GEN p1N = gel(W,1), W3 = gel(W,3), section = modsymb_get_section(W);
     895                 :        742 :   ulong N = p1N_get_N(p1N);
     896                 :        742 :   ulong c = umodiu(gcoeff(M,2,1), N);
     897                 :        742 :   ulong d = umodiu(gcoeff(M,2,2), N);
     898                 :        742 :   long ind = p1_index(c, d, p1N); /* as an elt of P1(Z/NZ) */
     899                 :        742 :   *index = W3[ind]; /* as an elt of F, E2, ... */
     900                 :        742 :   return ZM_zm_mul(M, sl2_inv(gel(section,ind)));
     901                 :            : }
     902                 :            : /* same for a path. Return [[ind], M] */
     903                 :            : static GEN
     904                 :      27706 : path_Gamma0N_decompose(GEN W, GEN path)
     905                 :            : {
     906                 :      27706 :   GEN p1N = gel(W,1);
     907                 :      27706 :   GEN p1index_to_ind = gel(W,3);
     908                 :      27706 :   GEN section = modsymb_get_section(W);
     909                 :      27706 :   GEN M = path_to_matrix(path);
     910                 :      27706 :   long p1index = p1_index(cc(M), dd(M), p1N);
     911                 :      27706 :   long ind = p1index_to_ind[p1index];
     912                 :      27706 :   GEN M0 = ZM_zm_mul(zm_to_ZM(M), sl2_inv(gel(section,p1index)));
     913                 :      27706 :   return mkvec2(mkvecsmall(ind), M0);
     914                 :            : }
     915                 :            : 
     916                 :            : /* Input:  [v1, ..., vn]
     917                 :            :  * Output: [vn, ..., v1]. Shallow function */
     918                 :            : static GEN
     919                 :      40754 : Reverse(GEN v)
     920                 :            : {
     921                 :            :   long i, l;
     922                 :      40754 :   GEN w = cgetg_copy(v, &l);
     923         [ +  + ]:     122262 :   for (i = 1; i < l; i++) gel(w, i) = gel(v, l-i);
     924                 :      40754 :   return w;
     925                 :            : }
     926                 :            : 
     927                 :            : /*Form generators of H_1(X_0(N),{cusps},Z)
     928                 :            : *
     929                 :            : *Input: N = integer > 1, p1N = P^1(Z/NZ)
     930                 :            : *Output: [cusp_list,E,F,T2,T3,E1] where
     931                 :            : *  cusps_list = list of cusps describing fundamental domain of
     932                 :            : *    \Gamma_0(N).
     933                 :            : *  E = list of paths in the boundary of the fundamental domains and oriented
     934                 :            : *    clockwise such that they do not contain a point
     935                 :            : *    fixed by an element of order 2 and they are not an edge of a
     936                 :            : *    triangle containing a fixed point of an element of order 3
     937                 :            : *  F = list of paths in the interior of the domain with each
     938                 :            : *    orientation appearing separately
     939                 :            : * T2 = list of paths in the boundary of domain containing a point fixed
     940                 :            : *    by an element of order 2 (oriented clockwise)
     941                 :            : * T3 = list of paths in the boundard of domain which are the edges of
     942                 :            : *    some triangle containing a fixed point of a matrix of order 3 (both
     943                 :            : *    orientations appear)
     944                 :            : * E1 = a sublist of E such that every path in E is \Gamma_0(N)-equivalent to
     945                 :            : *    either an element of E1 or the flip (reversed orientation) of an element
     946                 :            : *    of E1.
     947                 :            : * (Elements of T2 are \Gamma_0(N)-equivalent to their own flip.)
     948                 :            : *
     949                 :            : * sec = a list from 1..#p1N of matrices describing a section of the map
     950                 :            : *   SL_2(Z) to P^1(Z/NZ) given by [a,b;c,d]-->[c,d].
     951                 :            : *   Given our fixed enumeration of P^1(Z/NZ), the j-th element of the list
     952                 :            : *   represents the image of the j-th element of P^1(Z/NZ) under the section. */
     953                 :            : 
     954                 :            : /* insert path in set T */
     955                 :            : static void
     956                 :      40614 : set_insert(hashtable *T, GEN path)
     957                 :      40614 : { hash_insert(T, path,  (void*)(T->nb + 1)); }
     958                 :            : 
     959                 :            : static GEN
     960                 :       1449 : hash_to_vec(hashtable *h)
     961                 :            : {
     962                 :       1449 :   GEN v = cgetg(h->nb + 1, t_VEC);
     963                 :            :   ulong i;
     964         [ +  + ]:     293622 :   for (i = 0; i < h->len; i++)
     965                 :            :   {
     966                 :     292173 :     hashentry *e = h->table[i];
     967         [ +  + ]:     359604 :     while (e)
     968                 :            :     {
     969                 :      67431 :       GEN key = (GEN)e->key;
     970                 :      67431 :       long index = (long)e->val;
     971                 :      67431 :       gel(v, index) = key;
     972                 :      67431 :       e = e->next;
     973                 :            :     }
     974                 :            :   }
     975                 :       1449 :   return v;
     976                 :            : }
     977                 :            : 
     978                 :            : static long
     979                 :      20664 : path_to_p1_index(GEN path, GEN p1N)
     980                 :            : {
     981                 :      20664 :   GEN M = path_to_matrix(path);
     982                 :      20664 :   return p1_index(cc(M), dd(M), p1N);
     983                 :            : }
     984                 :            : 
     985                 :            : /* Pollack-Stevens sets */
     986                 :            : typedef struct PS_sets_t {
     987                 :            :   hashtable *F, *T2, *T31, *T32, *E1, *E2;
     988                 :            :   GEN E2_in_terms_of_E1, stdE1;
     989                 :            : } PS_sets_t;
     990                 :            : 
     991                 :            : static hashtable *
     992                 :        966 : set_init(long max)
     993                 :        966 : { return hash_create(max, (ulong(*)(void*))&hash_GEN,
     994                 :            :                           (int(*)(void*,void*))&gidentical, 1); }
     995                 :            : static void
     996                 :      13776 : insert_E(GEN path, PS_sets_t *S, GEN p1N)
     997                 :            : {
     998                 :      13776 :   GEN rev = Reverse(path);
     999                 :      13776 :   long std = path_to_p1_index(rev, p1N);
    1000                 :      13776 :   GEN v = gel(S->stdE1, std);
    1001         [ +  + ]:      13776 :   if (v)
    1002                 :            :   { /* [s, p1], where E1[s] = the path p1 \equiv Reverse(path) mod \Gamma */
    1003                 :       6888 :     GEN gamma, p1 = gel(v,2);
    1004                 :       6888 :     long r, s = itos(gel(v,1));
    1005                 :            : 
    1006                 :       6888 :     set_insert(S->E2, path);
    1007                 :       6888 :     r = S->E2->nb;
    1008         [ -  + ]:       6888 :     if (gel(S->E2_in_terms_of_E1, r) != gen_0) pari_err_BUG("insert_E");
    1009                 :            : 
    1010                 :       6888 :     gamma = gamma_equiv_matrix(rev, p1);
    1011                 :            :     /* E2[r] + gamma * E1[s] = 0 */
    1012                 :       6888 :     gel(S->E2_in_terms_of_E1, r) = mkvec2(utoipos(s),
    1013                 :            :                                           to_famat_shallow(gamma,gen_m1));
    1014                 :            :   }
    1015                 :            :   else
    1016                 :            :   {
    1017                 :       6888 :     set_insert(S->E1, path);
    1018                 :       6888 :     std = path_to_p1_index(path, p1N);
    1019                 :       6888 :     gel(S->stdE1, std) = mkvec2(utoipos(S->E1->nb), path);
    1020                 :            :   }
    1021                 :      13776 : }
    1022                 :            : 
    1023                 :            : static GEN
    1024                 :        644 : cusp_infinity() { return mkvecsmall2(1,0); }
    1025                 :            : 
    1026                 :            : static void
    1027                 :        161 : form_E_F_T(ulong N, GEN p1N, GEN *pC, PS_sets_t *S)
    1028                 :            : {
    1029                 :        161 :   GEN C, cusp_list = form_list_of_cusps(N, p1N);
    1030                 :        161 :   long nbgen = lg(cusp_list)-1, nbmanin = p1_size(p1N), r, s, i;
    1031                 :            :   hashtable *F, *T2, *T31, *T32, *E1, *E2;
    1032                 :            : 
    1033                 :        161 :   *pC = C = cgetg(nbgen+1, t_VEC);
    1034         [ +  + ]:      13853 :   for (i = 1; i <= nbgen; i++)
    1035                 :            :   {
    1036                 :      13692 :     GEN c = gel(cusp_list,i);
    1037                 :      13692 :     gel(C,i) = mkvecsmall2(c[1], c[2]);
    1038                 :            :   }
    1039                 :        161 :   S->F  = F  = set_init(nbmanin);
    1040                 :        161 :   S->E1 = E1 = set_init(nbgen);
    1041                 :        161 :   S->E2 = E2 = set_init(nbgen);
    1042                 :        161 :   S->T2 = T2 = set_init(nbgen);
    1043                 :        161 :   S->T31 = T31 = set_init(nbgen);
    1044                 :        161 :   S->T32 = T32 = set_init(nbgen);
    1045                 :            : 
    1046                 :            :   /* T31 represents the three torsion paths going from left to right */
    1047                 :            :   /* T32 represents the three torsion paths going from right to left */
    1048         [ +  + ]:      13692 :   for (r = 1; r < nbgen; r++)
    1049                 :            :   {
    1050                 :      13531 :     GEN c2 = gel(cusp_list,r+1);
    1051         [ +  + ]:      13531 :     if (c2[3] == type_T)
    1052                 :            :     {
    1053                 :         21 :       GEN c1 = gel(cusp_list,r), path = mkpath(c1,c2), path2 = Reverse(path);
    1054                 :         21 :       set_insert(T31, path);
    1055                 :         21 :       set_insert(T32, path2);
    1056                 :            :     }
    1057                 :            :   }
    1058                 :            : 
    1059                 :            :   /* to record relations between E2 and E1 */
    1060                 :        161 :   S->E2_in_terms_of_E1 = zerovec(nbgen);
    1061                 :        161 :   S->stdE1 = const_vec(nbmanin, NULL);
    1062                 :            : 
    1063                 :            :   /* Assumption later: path [oo,0] is E1[1], path [1,oo] is E2[1] */
    1064                 :            :   {
    1065                 :        161 :     GEN oo = cusp_infinity();
    1066                 :        161 :     GEN p1 = mkpath(oo, mkvecsmall2(0,1)); /* [oo, 0] */
    1067                 :        161 :     GEN p2 = mkpath(mkvecsmall2(1,1), oo); /* [1, oo] */
    1068                 :        161 :     insert_E(p1, S, p1N);
    1069                 :        161 :     insert_E(p2, S, p1N);
    1070                 :            :   }
    1071                 :            : 
    1072         [ +  + ]:      13692 :   for (r = 1; r < nbgen; r++)
    1073                 :            :   {
    1074                 :      13531 :     GEN c1 = gel(cusp_list,r);
    1075         [ +  + ]:    3264345 :     for (s = r+1; s <= nbgen; s++)
    1076                 :            :     {
    1077                 :    3250814 :       pari_sp av = avma;
    1078                 :    3250814 :       GEN c2 = gel(cusp_list,s), path;
    1079                 :    3250814 :       GEN d = subii(mulss(c1[1],c2[2]), mulss(c1[2],c2[1]));
    1080                 :    3250814 :       avma = av;
    1081         [ +  + ]:    3250814 :       if (!is_pm1(d)) continue;
    1082                 :            : 
    1083                 :      26901 :       path = mkpath(c1,c2);
    1084         [ +  + ]:      26901 :       if (r+1 == s)
    1085                 :            :       {
    1086                 :      13531 :         GEN w = path;
    1087                 :      13531 :         ulong hash = T31->hash(w); /* T31, T32 use the same hash function */
    1088 [ +  + ][ +  - ]:      13531 :         if (!hash_search2(T31, w, hash) && !hash_search2(T32, w, hash))
    1089                 :            :         {
    1090         [ +  + ]:      13510 :           if (gamma_equiv(path, Reverse(path), N))
    1091                 :         56 :             set_insert(T2, path);
    1092                 :            :           else
    1093                 :      13454 :             insert_E(path, S, p1N);
    1094                 :            :         }
    1095                 :            :       } else {
    1096                 :      13370 :         set_insert(F, mkvec2(path, mkvecsmall2(r,s)));
    1097                 :      13370 :         set_insert(F, mkvec2(Reverse(path), mkvecsmall2(s,r)));
    1098                 :            :       }
    1099                 :            :     }
    1100                 :            :   }
    1101                 :        161 :   setlg(S->E2_in_terms_of_E1, E2->nb+1);
    1102                 :        161 : }
    1103                 :            : 
    1104                 :            : /* v = \sum n_i g_i, return \sum n_i g_i^(-1) */
    1105                 :            : static GEN
    1106                 :        749 : ZGl2Q_star(GEN v)
    1107                 :            : {
    1108                 :            :   long i, l;
    1109                 :            :   GEN w, G;
    1110         [ +  + ]:        749 :   if (typ(v) == t_INT) return v;
    1111                 :        483 :   G = gel(v,1);
    1112                 :        483 :   w = cgetg_copy(G, &l);
    1113         [ +  + ]:       1211 :   for (i = 1; i < l; i++)
    1114                 :            :   {
    1115                 :        728 :     GEN g = gel(G,i);
    1116         [ +  + ]:        728 :     if (typ(g) == t_MAT) g = SL2_inv(g);
    1117                 :        728 :     gel(w,i) = g;
    1118                 :            :   }
    1119                 :        749 :   return ZG_normalize(mkmat2(w, gel(v,2)));
    1120                 :            : }
    1121                 :            : static GEN
    1122                 :        294 : ZGl2QC_star(GEN v)
    1123                 :            : {
    1124                 :            :   long i, l;
    1125                 :        294 :   GEN w = cgetg_copy(v, &l);
    1126         [ +  + ]:        966 :   for (i = 1; i < l; i++) gel(w,i) = ZGl2Q_star(gel(v,i));
    1127                 :        294 :   return w;
    1128                 :            : }
    1129                 :            : 
    1130                 :            : /* Input: h = set of unimodular paths, p1N = P^1(Z/NZ) = Gamma_0(N)\PSL2(Z)
    1131                 :            :  * Output: Each path is converted to a matrix and then an element of P^1(Z/NZ)
    1132                 :            :  * Append the matrix to W[12], append the index that represents
    1133                 :            :  * these elements of P^1 (the classes mod Gamma_0(N) via our fixed
    1134                 :            :  * enumeration to W[2]. */
    1135                 :            : static void
    1136                 :        966 : paths_decompose(GEN W, hashtable *h, int flag)
    1137                 :            : {
    1138                 :        966 :   GEN p1N = modsymb_get_p1N(W), section = modsymb_get_section(W);
    1139                 :        966 :   GEN v = hash_to_vec(h);
    1140                 :        966 :   long i, l = lg(v);
    1141         [ +  + ]:      41580 :   for (i = 1; i < l; i++)
    1142                 :            :   {
    1143                 :      40614 :     GEN e = gel(v,i);
    1144         [ +  + ]:      40614 :     GEN M = path_to_matrix(flag? gel(e,1): e);
    1145                 :      40614 :     long index = p1_index(cc(M), dd(M), p1N);
    1146                 :      40614 :     vecsmalltrunc_append(gel(W,2), index);
    1147                 :      40614 :     gel(section, index) = M;
    1148                 :            :   }
    1149                 :        966 : }
    1150                 :            : static void
    1151                 :        161 : fill_W2_W12(GEN W, PS_sets_t *S)
    1152                 :            : {
    1153                 :        161 :   GEN p1N = gel(W,1);
    1154                 :        161 :   long n = p1_size(p1N);
    1155                 :        161 :   gel(W, 2) = vecsmalltrunc_init(n+1);
    1156                 :        161 :   gel(W,12) = cgetg(n+1, t_VEC);
    1157                 :            :   /* F contains [path, [index cusp1, index cusp2]]. Others contain paths only */
    1158                 :        161 :   paths_decompose(W, S->F, 1);
    1159                 :        161 :   paths_decompose(W, S->E2, 0);
    1160                 :        161 :   paths_decompose(W, S->T32, 0);
    1161                 :        161 :   paths_decompose(W, S->E1, 0);
    1162                 :        161 :   paths_decompose(W, S->T2, 0);
    1163                 :        161 :   paths_decompose(W, S->T31, 0);
    1164                 :        161 : }
    1165                 :            : 
    1166                 :            : /* x t_VECSMALL, corresponds to a map x(i) = j, where 1 <= j <= max for all i
    1167                 :            :  * Return y s.t. y[j] = i or 0 (not in image) */
    1168                 :            : static GEN
    1169                 :        322 : reverse_list(GEN x, long max)
    1170                 :            : {
    1171                 :        322 :   GEN y = const_vecsmall(max, 0);
    1172                 :        322 :   long r, lx = lg(x);
    1173         [ +  + ]:      47901 :   for (r = 1; r < lx; r++) y[ x[r] ] = r;
    1174                 :        322 :   return y;
    1175                 :            : }
    1176                 :            : 
    1177                 :            : /* go from C[a] to C[b]; return the indices of paths
    1178                 :            :  * E.g. if a < b
    1179                 :            :  *   (C[a]->C[a+1], C[a+1]->C[a+2], ... C[b-1]->C[b])
    1180                 :            :  * (else reverse direction)
    1181                 :            :  * = b - a paths */
    1182                 :            : static GEN
    1183                 :      26376 : F_indices(GEN W, long a, long b)
    1184                 :            : {
    1185                 :      26376 :   GEN v = cgetg(labs(b-a) + 1, t_VEC);
    1186                 :      26376 :   long s, k = 1;
    1187         [ +  + ]:      26376 :   if (a < b) {
    1188                 :      13188 :     GEN index_forward = gel(W,13);
    1189         [ +  + ]:     111790 :     for (s = a; s < b; s++) gel(v,k++) = gel(index_forward,s);
    1190                 :            :   } else {
    1191                 :      13188 :     GEN index_backward = gel(W,14);
    1192         [ +  + ]:     111790 :     for (s = a; s > b; s--) gel(v,k++) = gel(index_backward,s);
    1193                 :            :   }
    1194                 :      26376 :   return v;
    1195                 :            : }
    1196                 :            : /* go from C[a] to C[b] via oo; return the indices of paths
    1197                 :            :  * E.g. if a < b
    1198                 :            :  *   (C[a]->C[a-1], ... C[2]->C[1],
    1199                 :            :  *    C[1]->oo, oo-> C[end],
    1200                 :            :  *    C[end]->C[end-1], ... C[b+1]->C[b])
    1201                 :            :  *  a-1 + 2 + end-(b+1)+1 = end - b + a + 1 paths  */
    1202                 :            : static GEN
    1203                 :        364 : F_indices_oo(GEN W, long end, long a, long b)
    1204                 :            : {
    1205                 :        364 :   GEN index_oo = gel(W,15);
    1206                 :        364 :   GEN v = cgetg(end-labs(b-a)+1 + 1, t_VEC);
    1207                 :        364 :   long s, k = 1;
    1208                 :            : 
    1209         [ +  + ]:        364 :   if (a < b) {
    1210                 :        182 :     GEN index_backward = gel(W,14);
    1211         [ -  + ]:        182 :     for (s = a; s > 1; s--) gel(v,k++) = gel(index_backward,s);
    1212                 :        182 :     gel(v,k++) = gel(index_backward,1); /* C[1] -> oo */
    1213                 :        182 :     gel(v,k++) = gel(index_oo,2); /* oo -> C[end] */
    1214         [ +  + ]:       5726 :     for (s = end; s > b; s--) gel(v,k++) = gel(index_backward,s);
    1215                 :            :   } else {
    1216                 :        182 :     GEN index_forward = gel(W,13);
    1217         [ +  + ]:       5726 :     for (s = a; s < end; s++) gel(v,k++) = gel(index_forward,s);
    1218                 :        182 :     gel(v,k++) = gel(index_forward,end); /* C[end] -> oo */
    1219                 :        182 :     gel(v,k++) = gel(index_oo,1); /* oo -> C[1] */
    1220         [ -  + ]:        182 :     for (s = 1; s < b; s++) gel(v,k++) = gel(index_forward,s);
    1221                 :            :   }
    1222                 :        364 :   return v;
    1223                 :            : }
    1224                 :            : /* index of oo -> C[1], oo -> C[end] */
    1225                 :            : static GEN
    1226                 :        161 : indices_oo(GEN W, GEN C)
    1227                 :            : {
    1228                 :        161 :   long end = lg(C)-1;
    1229                 :        161 :   GEN w, v = cgetg(2+1, t_VEC), oo = cusp_infinity();
    1230                 :        161 :   w = mkpath(oo, gel(C,1)); /* oo -> C[1]=0 */
    1231                 :        161 :   gel(v,1) = path_Gamma0N_decompose(W, w);
    1232                 :        161 :   w = mkpath(oo, gel(C,end)); /* oo -> C[end]=1 */
    1233                 :        161 :   gel(v,2) = path_Gamma0N_decompose(W, w);
    1234                 :        161 :   return v;
    1235                 :            : }
    1236                 :            : 
    1237                 :            : /* index of C[1]->C[2], C[2]->C[3], ... C[end-1]->C[end], C[end]->oo
    1238                 :            :  * Recall that C[1] = 0, C[end] = 1 */
    1239                 :            : static GEN
    1240                 :        161 : indices_forward(GEN W, GEN C)
    1241                 :            : {
    1242                 :        161 :   long s, k = 1, end = lg(C)-1;
    1243                 :        161 :   GEN v = cgetg(end+1, t_VEC);
    1244         [ +  + ]:      13853 :   for (s = 1; s <= end; s++)
    1245                 :            :   {
    1246         [ +  + ]:      13692 :     GEN w = mkpath(gel(C,s), s == end? cusp_infinity(): gel(C,s+1));
    1247                 :      13692 :     gel(v,k++) = path_Gamma0N_decompose(W, w);
    1248                 :            :   }
    1249                 :        161 :   return v;
    1250                 :            : }
    1251                 :            : /* index of C[1]->oo, C[2]->C[1], ... C[end]->C[end-1] */
    1252                 :            : static GEN
    1253                 :        161 : indices_backward(GEN W, GEN C)
    1254                 :            : {
    1255                 :        161 :   long s, k = 1, end = lg(C)-1;
    1256                 :        161 :   GEN v = cgetg(end+1, t_VEC);
    1257         [ +  + ]:      13853 :   for (s = 1; s <= end; s++)
    1258                 :            :   {
    1259         [ +  + ]:      13692 :     GEN w = mkpath(gel(C,s), s == 1? cusp_infinity(): gel(C,s-1));
    1260                 :      13692 :     gel(v,k++) = path_Gamma0N_decompose(W, w);
    1261                 :            :   }
    1262                 :        161 :   return v;
    1263                 :            : }
    1264                 :            : 
    1265                 :            : /* N = integer > 1. Returns data describing Delta_0 = Z[P^1(Q)]_0 seen as
    1266                 :            :  * a Gamma_0(N) - module. */
    1267                 :            : static GEN
    1268                 :        161 : modsymbinit_N(ulong N)
    1269                 :            : {
    1270                 :        161 :   GEN p1N = create_p1mod(N);
    1271                 :            :   GEN C, vecF, vecT2, vecT31;
    1272                 :            :   ulong r, s, width;
    1273                 :        161 :   long nball, nbgen, nbp1N = p1_size(p1N);
    1274                 :        161 :   GEN TAU = mkmat2(mkcol2(gen_0,gen_1), mkcol2(gen_m1,gen_m1)); /*[0,-1;1,-1]*/
    1275                 :            :   GEN W, W2, singlerel, annT2, annT31;
    1276                 :            :   GEN F_index;
    1277                 :            :   hashtable *F, *T2, *T31, *T32, *E1, *E2;
    1278                 :            :   PS_sets_t S;
    1279                 :            : 
    1280                 :        161 :   form_E_F_T(N,p1N, &C, &S);
    1281                 :        161 :   E1  = S.E1;
    1282                 :        161 :   E2  = S.E2;
    1283                 :        161 :   T31 = S.T31;
    1284                 :        161 :   T32 = S.T32;
    1285                 :        161 :   F   = S.F;
    1286                 :        161 :   T2  = S.T2;
    1287                 :        161 :   nbgen = lg(C)-1;
    1288                 :            : 
    1289                 :        161 :   W = cgetg(17, t_VEC);
    1290                 :        161 :   gel(W,1) = p1N;
    1291                 :            : 
    1292                 :            :  /* Put our paths in the order: F,E2,T32,E1,T2,T31
    1293                 :            :   * W2[j] associates to the j-th element of this list its index in P1. */
    1294                 :        161 :   fill_W2_W12(W, &S);
    1295                 :        161 :   W2 = gel(W, 2);
    1296                 :        161 :   nball = lg(W2)-1;
    1297                 :        161 :   gel(W,3) = reverse_list(W2, nbp1N);
    1298                 :            : 
    1299                 :        161 :   gel(W,5) = vecslice(gel(W,2), F->nb + E2->nb + T32->nb + 1, nball);
    1300                 :        161 :   gel(W,4) = reverse_list(gel(W,5), nbp1N);
    1301                 :        161 :   gel(W,13) = indices_forward(W, C);
    1302                 :        161 :   gel(W,14) = indices_backward(W, C);
    1303                 :        161 :   gel(W,15) = indices_oo(W, C);
    1304                 :        161 :   gel(W,11) = mkvecsmall5(F->nb,
    1305                 :        161 :                           F->nb + E2->nb,
    1306                 :        161 :                           F->nb + E2->nb + T32->nb,
    1307                 :        161 :                           F->nb + E2->nb + T32->nb + E1->nb,
    1308                 :        161 :                           F->nb + E2->nb + T32->nb + E1->nb + T2->nb);
    1309                 :            : 
    1310                 :            :   /* relations between T32 and T31 [not stored!]
    1311                 :            :    * T32[i] = - T31[i] */
    1312                 :            : 
    1313                 :            :   /* relations of F */
    1314                 :        161 :   width = E1->nb + T2->nb + T31->nb;
    1315                 :            :   /* F_index[r] = [index_1, ..., index_k], where index_i is the p1_index()
    1316                 :            :    * of the elementary unimodular path between 2 consecutive cusps
    1317                 :            :    * [in E1,E2,T2,T31 or T32] */
    1318                 :        161 :   F_index = cgetg(F->nb+1, t_VEC);
    1319                 :        161 :   vecF = hash_to_vec(F);
    1320         [ +  + ]:      26901 :   for (r = 1; r <= F->nb; r++)
    1321                 :            :   {
    1322                 :      26740 :     GEN w = gel(gel(vecF,r), 2);
    1323                 :      26740 :     long a = w[1], b = w[2], d = labs(b - a);
    1324                 :            :     /* c1 = cusp_list[a],  c2 = cusp_list[b], ci != oo */
    1325                 :      53480 :     gel(F_index,r) = (nbgen-d >= d-1)? F_indices(W, a,b)
    1326         [ +  + ]:      26740 :                                      : F_indices_oo(W, lg(C)-1,a,b);
    1327                 :            :   }
    1328                 :            : 
    1329                 :        161 :   singlerel = cgetg(width+1, t_VEC);
    1330                 :            :   /* form the single boundary relation */
    1331         [ +  + ]:       7049 :   for (s = 1; s <= E2->nb; s++)
    1332                 :            :   {
    1333                 :       6888 :     GEN data = gel(S.E2_in_terms_of_E1,s);
    1334                 :       6888 :     long c = itos(gel(data,1));
    1335                 :       6888 :     GEN u = gel(data,2); /* E2[s] = u * E1[c], u = - [gamma] */
    1336                 :       6888 :     GEN gamma = gcoeff(u,1,1);
    1337                 :       6888 :     gel(singlerel, c) = mkmat2(mkcol2(gen_1,gamma),mkcol2(gen_1,gen_m1));
    1338                 :            :   }
    1339         [ +  + ]:        238 :   for (r = E1->nb + 1; r <= width; r++) gel(singlerel, r) = gen_1;
    1340                 :            : 
    1341                 :            :   /* form the 2-torsion relations */
    1342                 :        161 :   annT2 = cgetg(T2->nb+1, t_VEC);
    1343                 :        161 :   vecT2 = hash_to_vec(T2);
    1344         [ +  + ]:        217 :   for (r = 1; r <= T2->nb; r++)
    1345                 :            :   {
    1346                 :         56 :     GEN w = gel(vecT2,r);
    1347                 :         56 :     GEN gamma = gamma_equiv_matrix(Reverse(w), w);
    1348                 :         56 :     gel(annT2, r) = mkmat2(mkcol2(gen_1,gamma), mkcol2(gen_1,gen_1));
    1349                 :            :   }
    1350                 :            : 
    1351                 :            :   /* form the 3-torsion relations */
    1352                 :        161 :   annT31 = cgetg(T31->nb+1, t_VEC);
    1353                 :        161 :   vecT31 = hash_to_vec(T31);
    1354         [ +  + ]:        182 :   for (r = 1; r <= T31->nb; r++)
    1355                 :            :   {
    1356                 :         21 :     GEN M = zm_to_ZM( path_to_matrix( Reverse(gel(vecT31,r)) ) );
    1357                 :         21 :     GEN gamma = ZM_mul(ZM_mul(M, TAU), SL2_inv(M));
    1358                 :         21 :     gel(annT31, r) = mkmat2(mkcol3(gen_1,gamma,ZM_sqr(gamma)),
    1359                 :            :                             mkcol3(gen_1,gen_1,gen_1));
    1360                 :            :   }
    1361                 :        161 :   gel(W,6) = F_index;
    1362                 :        161 :   gel(W,7) = S.E2_in_terms_of_E1;
    1363                 :        161 :   gel(W,8) = annT2;
    1364                 :        161 :   gel(W,9) = annT31;
    1365                 :        161 :   gel(W,10)= singlerel;
    1366                 :        161 :   gel(W,16)= inithashcusps(p1N);
    1367                 :        161 :   return W;
    1368                 :            : }
    1369                 :            : 
    1370                 :            : /* Modular symbols in weight k: Hom_Gamma(Delta, Q[x,y]_{k-2}) */
    1371                 :            : /* A symbol phi is represented by the {phi(g_i)}, {phi(g'_i)}, {phi(g''_i)}
    1372                 :            :  * where the {g_i, g'_i, g''_i} are the Z[\Gamma]-generators of Delta,
    1373                 :            :  * g_i corresponds to E1, g'_i to T2, g''_i to T31.
    1374                 :            :  */
    1375                 :            : 
    1376                 :            : /* FIXME: export. T^1, ..., T^n */
    1377                 :            : static GEN
    1378                 :       1834 : RgX_powers(GEN T, long n)
    1379                 :            : {
    1380                 :       1834 :   GEN v = cgetg(n+1, t_VEC);
    1381                 :            :   long i;
    1382                 :       1834 :   gel(v, 1) = T;
    1383         [ +  + ]:       5292 :   for (i = 1; i < n; i++) gel(v,i+1) = RgX_mul(gel(v,i), T);
    1384                 :       1834 :   return v;
    1385                 :            : }
    1386                 :            : 
    1387                 :            : /* g = [a,b;c,d]. Return (X^{k-2} | g)(X,Y)[X = 1]. */
    1388                 :            : static GEN
    1389                 :          0 : voo_act_Gl2Q(GEN g, long k)
    1390                 :            : {
    1391                 :          0 :   GEN c = gcoeff(g,2,1), d = gcoeff(g,2,2);
    1392                 :          0 :   return RgX_to_RgC(gpowgs(deg1pol_shallow(gneg(c), d, 0), k-2), k-1);
    1393                 :            : }
    1394                 :            : 
    1395                 :            : /* g = [a,b;c,d]. Return (P | g)(X,Y)[X = 1] = P(dX - cY, -b X + aY)[X = 1],
    1396                 :            :  * for P = X^{k-2}, X_^{k-3}Y, ..., Y^{k-2} */
    1397                 :            : GEN
    1398                 :        917 : RgX_act_Gl2Q(GEN g, long k)
    1399                 :            : {
    1400                 :        917 :   GEN a = gcoeff(g,1,1), b = gcoeff(g,1,2);
    1401                 :        917 :   GEN c = gcoeff(g,2,1), d = gcoeff(g,2,2);
    1402                 :        917 :   GEN V1 = RgX_powers(deg1pol_shallow(gneg(c), d, 0), k-2);
    1403                 :        917 :   GEN V2 = RgX_powers(deg1pol_shallow(a, gneg(b), 0), k-2);
    1404                 :        917 :   GEN V = cgetg(k, t_MAT);
    1405                 :            :   long i;
    1406                 :        917 :   gel(V,1)   = RgX_to_RgC(gel(V1, k-2), k-1);
    1407         [ +  + ]:       2646 :   for (i = 1; i < k-2; i++)
    1408                 :            :   {
    1409                 :       1729 :     GEN v1 = gel(V1, k-2-i); /* (d-cY)^(k-2-i) */
    1410                 :       1729 :     GEN v2 = gel(V2, i); /* (-b+aY)^i */
    1411                 :       1729 :     gel(V,i+1) = RgX_to_RgC(RgX_mul(v1,v2), k-1);
    1412                 :            :   }
    1413                 :        917 :   gel(V,k-1) = RgX_to_RgC(gel(V2, k-2), k-1);
    1414                 :        917 :   return V; /* V[i+1] = X^i | g */
    1415                 :            : }
    1416                 :            : /* z in Z[\Gamma], return the matrix of z acting on (X^{k-2},...,Y^{k-2}) */
    1417                 :            : GEN
    1418                 :        749 : RgX_act_ZGl2Q(GEN z, long k)
    1419                 :            : {
    1420                 :            :   long l, j;
    1421                 :        749 :   GEN S = NULL, G, E;
    1422         [ -  + ]:        749 :   if (typ(z) == t_INT) return matid(k-1);
    1423                 :        749 :   G = gel(z,1); l = lg(G);
    1424                 :        749 :   E = gel(z,2);
    1425         [ +  + ]:       1743 :   for (j = 1; j < l; j++)
    1426                 :            :   {
    1427                 :        994 :     GEN M, g = gel(G,j), n = gel(E,j);
    1428         [ +  + ]:        994 :     if (typ(g) == t_INT)
    1429                 :         77 :       M = scalarmat_shallow(n, k-1);
    1430                 :            :     else
    1431                 :            :     {
    1432                 :        917 :       M = RgX_act_Gl2Q(g, k);
    1433         [ +  + ]:        917 :       if (is_pm1(n))
    1434         [ +  + ]:        651 :       { if (signe(n) < 0) M = RgM_neg(M);
    1435                 :            :       } else
    1436                 :        266 :         M = RgM_Rg_mul(M, n);
    1437                 :            :     }
    1438         [ +  + ]:        994 :     S = j == 1? M: RgM_add(S, M);
    1439                 :            :   }
    1440                 :        749 :   return S;
    1441                 :            : }
    1442                 :            : /* Given a vector of elements in Z[G], return it as vector of operators on V
    1443                 :            :  * (given by t_MAT) */
    1444                 :            : static GEN
    1445                 :        294 : ZGl2QC_to_act(GEN v, long k)
    1446                 :            : {
    1447                 :            :   long i, l;
    1448                 :        294 :   GEN w = cgetg_copy(v, &l);
    1449         [ +  + ]:        966 :   for (i = 1; i < l; i++) gel(w,i) = RgX_act_ZGl2Q(gel(v,i), k);
    1450                 :        294 :   return w;
    1451                 :            : }
    1452                 :            : 
    1453                 :            : /* For all V[i] in Z[\Gamma], find the P such that  P . V[i]^* = 0;
    1454                 :            :  * write P in basis X^{k-2}, ..., Y^{k-2} */
    1455                 :            : static GEN
    1456                 :         70 : ZGV_tors(GEN V, long k)
    1457                 :            : {
    1458                 :         70 :   long i, l = lg(V);
    1459                 :         70 :   GEN v = cgetg(l, t_VEC);
    1460         [ +  + ]:         91 :   for (i = 1; i < l; i++)
    1461                 :            :   {
    1462                 :         21 :     GEN a = ZGl2Q_star(gel(V,i));
    1463                 :         21 :     gel(v,i) = ZM_ker(RgX_act_ZGl2Q(a,k));
    1464                 :            :   }
    1465                 :         70 :   return v;
    1466                 :            : }
    1467                 :            : 
    1468                 :            : static long
    1469                 :   90415640 : set_from_index(GEN W11, long i)
    1470                 :            : {
    1471         [ +  + ]:   90415640 :   if (i <= W11[1]) return 1;
    1472         [ +  + ]:   78060325 :   if (i <= W11[2]) return 2;
    1473         [ +  + ]:   48138727 :   if (i <= W11[3]) return 3;
    1474         [ +  + ]:   48095607 :   if (i <= W11[4]) return 4;
    1475         [ +  + ]:    3515253 :   if (i <= W11[5]) return 5;
    1476                 :   90415640 :   return 6;
    1477                 :            : }
    1478                 :            : 
    1479                 :            : static void
    1480                 :       1008 : treat_index(GEN W, GEN M, long index, int negate, GEN v)
    1481                 :            : {
    1482                 :       1008 :   GEN W11 = gel(W,11);
    1483                 :       1008 :   long shift = W11[3]; /* #F + #E2 + T32 */
    1484   [ +  +  +  + ]:       1008 :   switch(set_from_index(W11, index))
    1485                 :            :   {
    1486                 :            :     case 1: /*F*/
    1487                 :            :     {
    1488                 :        133 :       GEN F_index = gel(W,6), ind = gel(F_index, index);
    1489                 :        133 :       long j, l = lg(ind);
    1490         [ +  + ]:        399 :       for (j = 1; j < l; j++)
    1491                 :            :       {
    1492                 :        266 :         GEN IND = gel(ind,j), M0 = gel(IND,2);
    1493                 :        266 :         long index = mael(IND,1,1);
    1494                 :        266 :         treat_index(W, ZM_mul(M,M0), index, negate, v);
    1495                 :            :       }
    1496                 :        133 :       break;
    1497                 :            :     }
    1498                 :            : 
    1499                 :            :     case 2: /*E2, E2[r] + gamma * E1[s] = 0 */
    1500                 :            :     {
    1501                 :        182 :       long r = index - W11[1];
    1502                 :        182 :       GEN E2_in_terms_of_E1= gel(W,7), z = gel(E2_in_terms_of_E1, r);
    1503                 :        182 :       long s = itou(gel(z,1));
    1504                 :            : 
    1505                 :        182 :       index = s;
    1506                 :        182 :       M = G_ZG_mul(M, gel(z,2)); /* M * (-gamma) */
    1507                 :        434 :       gel(v, index) = negate? ZG_sub(gel(v, index), M)
    1508         [ +  + ]:        182 :                             : ZG_add(gel(v, index), M);
    1509                 :        182 :       break;
    1510                 :            :     }
    1511                 :            : 
    1512                 :            :     case 3: /*T32, T32[i] = -T31[i] */
    1513                 :            :     {
    1514                 :         21 :       long T3shift = W11[5] - W11[2]; /* #T32 + #E1 + #T2 */
    1515                 :         21 :       index += T3shift;
    1516                 :         21 :       index -= shift;
    1517         [ +  + ]:         21 :       gel(v, index) = ZG_add(gel(v, index),
    1518                 :            :                              to_famat_shallow(M,negate?gen_1:gen_m1));
    1519                 :         21 :       break;
    1520                 :            :     }
    1521                 :            :     default: /*E1,T2,T31*/
    1522                 :        672 :       index -= shift;
    1523         [ +  + ]:        672 :       gel(v, index) = ZG_add(gel(v, index),
    1524                 :            :                              to_famat_shallow(M,negate?gen_m1:gen_1));
    1525                 :        672 :       break;
    1526                 :            :   }
    1527                 :       1008 : }
    1528                 :            : static void
    1529                 :   90414632 : treat_index_trivial(GEN W, long index, int negate, GEN v)
    1530                 :            : {
    1531                 :   90414632 :   GEN W11 = gel(W,11);
    1532                 :   90414632 :   long shift = W11[3]; /* #F + #E2 + T32 */
    1533   [ +  +  +  +  :   90414632 :   switch(set_from_index(W11, index))
                      - ]
    1534                 :            :   {
    1535                 :            :     case 1: /*F*/
    1536                 :            :     {
    1537                 :   12355182 :       GEN F_index = gel(W,6), ind = gel(F_index, index);
    1538                 :   12355182 :       long j, l = lg(ind);
    1539         [ +  + ]:   80208604 :       for (j = 1; j < l; j++)
    1540                 :            :       {
    1541                 :   67853422 :         GEN IND = gel(ind,j);
    1542                 :   67853422 :         treat_index_trivial(W, mael(IND,1,1), negate, v);
    1543                 :            :       }
    1544                 :   12355182 :       break;
    1545                 :            :     }
    1546                 :            : 
    1547                 :            :     case 2: /*E2, E2[r] + gamma * E1[s] = 0 */
    1548                 :            :     {
    1549                 :   29921416 :       long r = index - W11[1];
    1550                 :   29921416 :       GEN E2_in_terms_of_E1= gel(W,7), z = gel(E2_in_terms_of_E1, r);
    1551                 :   29921416 :       long s = itou(gel(z,1));
    1552                 :   29921416 :       index = s;
    1553                 :   61012007 :       gel(v, index) = negate? addiu(gel(v, index), 1)
    1554         [ +  + ]:   29921416 :                             : subiu(gel(v, index), 1);
    1555                 :   29921416 :       break;
    1556                 :            :     }
    1557                 :            : 
    1558                 :            :     case 3: case 5: case 6: /*T32,T2,T31*/
    1559                 :    3558352 :       break;
    1560                 :            : 
    1561                 :            :     case 4: /*E1*/
    1562                 :   44579682 :       index -= shift;
    1563                 :   93576049 :       gel(v, index) = negate? subiu(gel(v, index), 1)
    1564         [ +  + ]:   44579682 :                             : addiu(gel(v, index), 1);
    1565                 :   44579682 :       break;
    1566                 :            :   }
    1567                 :   90414632 : }
    1568                 :            : static GEN
    1569                 :      64932 : cusp_to_frac(GEN c) { return gdivgs(stoi(c[1]), c[2]); }
    1570                 :            : 
    1571                 :            : /* cusp in P^1(Q): expresses {1/0 -> cusp} as \sum x_i g_i, see modsymblog */
    1572                 :            : static void
    1573                 :        588 : cusplog(GEN v, GEN W, int negate, GEN cusp)
    1574                 :            : {
    1575                 :            :   GEN PQ, P,Q, a,b,c,d;
    1576                 :            :   long i, lx;
    1577                 :            : 
    1578         [ +  + ]:       1043 :   if (!cusp[2]) return;
    1579                 :            : 
    1580                 :        455 :   PQ = ZV_allpnqn( gboundcf(cusp_to_frac(cusp), 0) );
    1581                 :        455 :   P = gel(PQ,1);
    1582                 :        455 :   Q = gel(PQ,2);
    1583                 :        455 :   lx = lg(P);
    1584                 :        455 :   a = gen_1;
    1585                 :        455 :   c = gen_0;
    1586         [ +  + ]:       1197 :   for (i = 1; i < lx; i++, c = d, a = b)
    1587                 :            :   {
    1588                 :            :     long index;
    1589                 :            :     GEN M;
    1590                 :        742 :     b = gel(P,i);
    1591                 :        742 :     d = gel(Q,i);
    1592         [ +  + ]:        742 :     if (!odd(i)) { a = negi(a); c = negi(c); }
    1593                 :        742 :     M = Gamma0N_decompose(W, mkmat2(mkcol2(a,c), mkcol2(b,d)), &index);
    1594                 :        742 :     treat_index(W, M, index, negate, v);
    1595                 :            :   }
    1596                 :            : }
    1597                 :            : /* cusplog for cusp [a,b] in case the action is trivial, q = a/b */
    1598                 :            : static void
    1599                 :    1981133 : cusplog_trivial_frac(GEN v, GEN W, int negate, GEN q)
    1600                 :            : {
    1601                 :    1981133 :   GEN Q, W3 = gel(W,3), p1N = gel(W,1);
    1602                 :    1981133 :   ulong c,d, N = p1N_get_N(p1N);
    1603                 :            :   long i, lx;
    1604                 :            : 
    1605                 :    1981133 :   Q = modulartosym_init(N, q);
    1606                 :    1981133 :   lx = lg(Q);
    1607                 :    1981133 :   c = 0;
    1608         [ +  + ]:   22625687 :   for (i = 1; i < lx; i++, c = d)
    1609                 :            :   {
    1610                 :            :     long index;
    1611                 :   20644554 :     d = Q[i];
    1612 [ +  + ][ +  + ]:   20644554 :     if (c && !odd(i)) c = N - c;
    1613                 :   20644554 :     index = W3[ p1_index(c,d,p1N) ];
    1614                 :   20644554 :     treat_index_trivial(W, index, negate, v);
    1615                 :            :   }
    1616                 :    1981133 : }
    1617                 :            : /* cusplog in case the action is trivial */
    1618                 :            : static void
    1619                 :      66150 : cusplog_trivial(GEN v, GEN W, int negate, GEN cusp)
    1620                 :            : {
    1621         [ +  + ]:      66150 :   if (!cusp[2]) return;
    1622                 :      66150 :   return cusplog_trivial_frac(v,W,negate,cusp_to_frac(cusp));
    1623                 :            : }
    1624                 :            : 
    1625                 :            : /* Expresses path as \sum x_i g_i, where the g_i are our distinguished
    1626                 :            :  * generators and x_i \in Z[\Gamma]. Returns [x_1,...,x_n].
    1627                 :            :  * W from solve_manin_relations */
    1628                 :            : GEN
    1629                 :        294 : modsymblog(GEN W, GEN path)
    1630                 :            : {
    1631                 :        294 :   pari_sp av = avma;
    1632                 :            :   GEN v;
    1633         [ -  + ]:        294 :   if (typ(W) != t_VEC) pari_err_TYPE("modsymblog", W);
    1634                 :        294 :   W = get_modsymb(W);
    1635                 :        294 :   v = zerovec(modsymb_get_nbgen(W));
    1636                 :        294 :   cusplog(v, W, 0, gel(path,2));
    1637                 :        294 :   cusplog(v, W, 1, gel(path,1));
    1638                 :        294 :   return gerepilecopy(av, v);
    1639                 :            : }
    1640                 :            : /* in case the action is trivial: v += modsymblog(path) */
    1641                 :            : static void
    1642                 :      33075 : modsymblog_trivial(GEN v, GEN W, GEN path)
    1643                 :            : {
    1644                 :      33075 :   cusplog_trivial(v, W, 0, gel(path,2));
    1645                 :      33075 :   cusplog_trivial(v, W, 1, gel(path,1));
    1646                 :      33075 : }
    1647                 :            : 
    1648                 :            : /** HECKE OPERATORS **/
    1649                 :            : /* [a,b;c,d] * cusp */
    1650                 :            : static GEN
    1651                 :      66738 : cusp_mul(GEN cusp, long a, long b, long c, long d)
    1652                 :            : {
    1653                 :      66738 :   long x = cusp[1], y = cusp[2];
    1654                 :      66738 :   return mkvecsmall2(a*x+b*y, c*x+d*y);
    1655                 :            : }
    1656                 :            : static GEN
    1657                 :      33369 : path_mul(GEN path, long a, long b, long c, long d)
    1658                 :            : {
    1659                 :      33369 :   GEN c1 = cusp_mul(gel(path,1), a,b,c,d);
    1660                 :      33369 :   GEN c2 = cusp_mul(gel(path,2), a,b,c,d);
    1661                 :      33369 :   return mkpath(c1,c2);
    1662                 :            : }
    1663                 :            : /* f in Gl2(Q), act on path (zm) */
    1664                 :            : static GEN
    1665                 :      33369 : Gl2Q_act_path(GEN f, GEN path)
    1666                 :      33369 : { return path_mul(path, coeff(f,1,1),coeff(f,1,2),coeff(f,2,1),coeff(f,2,2)); }
    1667                 :            : 
    1668                 :            : static GEN
    1669                 :    1936095 : init_act_trivial(GEN W) { return zerocol(modsymb_get_nbE1(W)); }
    1670                 :            : 
    1671                 :            : /* map from WW1 -> WW2, weight 2, trivial action. v a Gl2_Q or a t_VEC
    1672                 :            :  * of Gl2_Q (\sum v[i] in Z[Gl2(Q)]). Return the matrix associated to the
    1673                 :            :  * action of v. */
    1674                 :            : static GEN
    1675                 :        371 : getMorphism_trivial(GEN WW1, GEN WW2, GEN v)
    1676                 :            : {
    1677                 :        371 :   GEN W1 = get_modsymb(WW1), W2 = get_modsymb(WW2);
    1678                 :        371 :   GEN section = modsymb_get_section(W1), gen = gel(W1,5);
    1679                 :        371 :   long i, j, lv, nbE1 = modsymb_get_nbE1(W1);
    1680                 :        371 :   GEN T = cgetg(nbE1+1, t_MAT);
    1681         [ +  + ]:        371 :   if (typ(v) != t_VEC) v = mkvec(v);
    1682                 :        371 :   lv = lg(v);
    1683         [ +  + ]:      19810 :   for (i = 1; i <= nbE1; i++)
    1684                 :            :   {
    1685                 :      19439 :     long e = gen[i]; /* path-index of E1-element */
    1686                 :      19439 :     GEN w = gel(section, e); /* path_to_matrix() */
    1687                 :      19439 :     GEN t = init_act_trivial(W2);
    1688         [ +  + ]:      52514 :     for (j = 1; j < lv; j++)
    1689                 :      33075 :       modsymblog_trivial(t, W2, Gl2Q_act_path(gel(v,j), w));
    1690                 :      19439 :     gel(T,i) = t;
    1691                 :            :   }
    1692                 :        371 :   return T;
    1693                 :            : }
    1694                 :            : 
    1695                 :            : /* f zm/ZM in Gl_2(Q), acts from the left on Delta, which is generated by
    1696                 :            :  * (g_i) as Z[Gamma1]-module, and by (G_i) as Z[Gamma2]-module.
    1697                 :            :  * We have f.G_j = \sum_i \lambda_{i,j} g_i,   \lambda_{i,j} in Z[Gamma1]
    1698                 :            :  *
    1699                 :            :  * For phi in Hom_Gamma1(D,V), g in D, phi | f is in Hom_Gamma2(D,V) and
    1700                 :            :  *  (phi | f)(G_j) = phi(f.G_j) | f
    1701                 :            :  *                 = phi( \sum_i \lambda_{i,j} g_i ) | f
    1702                 :            :  *                 = \sum_i phi(g_i) | (\lambda_{i,j}^* f)
    1703                 :            :  *                 = \sum_i phi(g_i) | \mu_{i,j}
    1704                 :            :  * Return the \mu_{i,j} matrix as operators on V (t_MAT) */
    1705                 :            : static GEN
    1706                 :        133 : init_dual_act(GEN f, GEN W1, GEN W2)
    1707                 :            : {
    1708                 :        133 :   GEN section = modsymb_get_section(W2), gen = modsymb_get_genindex(W2);
    1709                 :        133 :   long j, dim = lg(gen)-1, k = modsymbk_get_weight(W1);
    1710                 :        133 :   GEN T = cgetg(dim+1, t_MAT), F;
    1711         [ -  + ]:        133 :   if (typ(gel(f,1)) == t_VEC)
    1712                 :            :   {
    1713                 :          0 :     F = f;
    1714                 :          0 :     f = ZM_to_zm(F);
    1715                 :            :   }
    1716                 :            :   else
    1717                 :        133 :     F = zm_to_ZM(f);
    1718                 :            :   /* f zm = F ZM */
    1719         [ +  + ]:        427 :   for (j = 1; j <= dim; j++)
    1720                 :            :   {
    1721                 :        294 :     pari_sp av = avma;
    1722                 :        294 :     GEN w = gel(section, gen[j]); /* path_to_matrix( E1/T2/T3 element ) */
    1723                 :        294 :     GEN l = modsymblog(W1, Gl2Q_act_path(f, w)); /* lambda_{i,j} */
    1724                 :        294 :     l = ZGl2QC_star(l); /* lambda_{i,j}^* */
    1725                 :        294 :     l = ZGC_G_mul(l, F); /* mu_{i,j} */
    1726                 :        294 :     gel(T,j) = gerepilecopy(av, ZGl2QC_to_act(l, k)); /* as operators on V */
    1727                 :            :   }
    1728                 :        133 :   return T;
    1729                 :            : }
    1730                 :            : /* phi in Hom_Gamma1(Delta, V), return the matrix whose colums are the
    1731                 :            :  *   \sum_i phi(g_i) | \mu_{i,j} = (phi|f)(G_j),
    1732                 :            :  * see init_dual_act. */
    1733                 :            : static GEN
    1734                 :        511 : dual_act(GEN phi, GEN mu)
    1735                 :            : {
    1736                 :        511 :   long l = lg(mu), a, j;
    1737                 :        511 :   GEN ind = gel(phi,2), pols = gel(phi,3);
    1738                 :        511 :   GEN v = cgetg(l, t_MAT);
    1739         [ +  + ]:       1701 :   for (j = 1; j < l; j++)
    1740                 :            :   {
    1741                 :       1190 :     GEN T = NULL;
    1742         [ +  + ]:       3528 :     for (a = 1; a < lg(ind); a++)
    1743                 :            :     {
    1744                 :       2338 :       long i = ind[a];
    1745                 :       2338 :       GEN t = gel(pols, a); /* phi(G_i) */
    1746                 :       2338 :       t = RgM_RgC_mul(gcoeff(mu,i,j), t);
    1747         [ +  + ]:       2338 :       T = T? RgC_add(T, t): t;
    1748                 :            :     }
    1749                 :       1190 :     gel(v,j) = T;
    1750                 :            :   }
    1751                 :        511 :   return v;
    1752                 :            : }
    1753                 :            : 
    1754                 :            : /* phi in Hom(Delta, V), phi(G_k) = vecT[k]. Write phi as
    1755                 :            :  *   \sum_{i,j} mu_{i,j} phi_{i,j}, mu_{i,j} in Q */
    1756                 :            : static GEN
    1757                 :        266 : getMorphism_basis(GEN W, GEN vecT)
    1758                 :            : {
    1759                 :        266 :   GEN basis = modsymbk_get_basis(W);
    1760                 :        266 :   long i, j, r, lvecT = lg(vecT), dim = lg(basis)-1;
    1761                 :        266 :   GEN st = modsymbk_get_st(W);
    1762                 :        266 :   GEN link = modsymbk_get_link(W);
    1763                 :        266 :   GEN invphiblock = modsymbk_get_invphiblock(W);
    1764                 :        266 :   long s = st[1], t = st[2];
    1765                 :        266 :   GEN R = zerovec(dim), Q, Ls, T0, T1, Ts, mu_st;
    1766         [ +  + ]:        616 :   for (r = 2; r < lvecT; r++)
    1767                 :            :   {
    1768                 :            :     GEN Tr, L;
    1769         [ +  + ]:        350 :     if (r == s) continue;
    1770                 :         84 :     Tr = gel(vecT,r); /* Phi(G_r), r != 1,s */
    1771                 :         84 :     L = gel(link, r);
    1772                 :         84 :     Q = ZC_apply_dinv(gel(invphiblock,r), Tr);
    1773                 :            :     /* write Phi(G_r) as sum_{a,b} mu_{a,b} Phi_{a,b}(G_r) */
    1774         [ +  + ]:        336 :     for (j = 1; j < lg(L); j++) gel(R, L[j]) = gel(Q,j);
    1775                 :            :   }
    1776                 :        266 :   Ls = gel(link, s);
    1777                 :        266 :   T1 = gel(vecT,1); /* Phi(G_1) */
    1778                 :        266 :   gel(R, Ls[t]) = mu_st = gel(T1, 1);
    1779                 :            : 
    1780                 :        266 :   T0 = NULL;
    1781         [ +  + ]:        616 :   for (i = 2; i < lg(link); i++)
    1782                 :            :   {
    1783                 :            :     GEN L;
    1784         [ +  + ]:        350 :     if (i == s) continue;
    1785                 :         84 :     L = gel(link,i);
    1786         [ +  + ]:        336 :     for (j =1 ; j < lg(L); j++)
    1787                 :            :     {
    1788                 :        252 :       long n = L[j]; /* phi_{i,j} = basis[n] */
    1789                 :        252 :       GEN mu_ij = gel(R, n);
    1790                 :        252 :       GEN phi_ij = gel(basis, n), pols = gel(phi_ij,3);
    1791                 :        252 :       GEN z = RgC_Rg_mul(gel(pols, 3), mu_ij);
    1792         [ +  + ]:        252 :       T0 = T0? RgC_add(T0, z): z; /* += mu_{i,j} Phi_{i,j} (G_s) */
    1793                 :            :     }
    1794                 :            :   }
    1795                 :        266 :   Ts = gel(vecT,s); /* Phi(G_s) */
    1796         [ +  + ]:        266 :   if (T0) Ts = RgC_sub(Ts, T0);
    1797                 :            :   /* solve \sum_{j!=t} mu_{s,j} Phi_{s,j}(G_s) = Ts */
    1798                 :        266 :   Q = ZC_apply_dinv(gel(invphiblock,s), Ts);
    1799         [ +  + ]:        840 :   for (j = 1; j < t; j++) gel(R, Ls[j]) = gel(Q,j);
    1800                 :            :   /* avoid mu_{s,t} */
    1801         [ +  + ]:        308 :   for (j = t; j < lg(Q); j++) gel(R, Ls[j+1]) = gel(Q,j);
    1802                 :        266 :   return R;
    1803                 :            : }
    1804                 :            : 
    1805                 :            : #if DEBUG
    1806                 :            : static void
    1807                 :            : checkrelation(GEN W, GEN T)
    1808                 :            : {
    1809                 :            :   GEN s, annT2, annT31, singlerel;
    1810                 :            :   long k = modsymbk_get_weight(W);
    1811                 :            :   long nbE1 = modsymb_get_nbE1(W), nbT2, nbT31;
    1812                 :            :   long i, l;
    1813                 :            :   W = get_modsymb(W);
    1814                 :            :   annT2 = gel(W,8); nbT2 = lg(annT2)-1;
    1815                 :            :   annT31 = gel(W,9);nbT31 = lg(annT31)-1;
    1816                 :            :   singlerel = gel(W,10);
    1817                 :            :   l = lg(singlerel);
    1818                 :            :   s = NULL;
    1819                 :            :   for (i = 1; i <= nbE1; i++)
    1820                 :            :   {
    1821                 :            :     GEN z = gel(singlerel, i);
    1822                 :            :     GEN M = RgX_act_ZGl2Q(ZGl2Q_star(z), k);
    1823                 :            :     GEN a = RgM_RgC_mul(M, gel(T,i));
    1824                 :            :     s = s? gadd(s, a): a;
    1825                 :            :   }
    1826                 :            :   for (; i < l; i++)
    1827                 :            :   {
    1828                 :            :     GEN a = gel(T,i);
    1829                 :            :     s = s? gadd(s, a): a;
    1830                 :            :   }
    1831                 :            :   if (!gcmp0(s)) pari_err_BUG("checkrelation (singlerel)");
    1832                 :            :   for (i = 1; i <= nbT2; i++)
    1833                 :            :   {
    1834                 :            :     GEN z = gel(annT2, i);
    1835                 :            :     GEN M = RgX_act_ZGl2Q(ZGl2Q_star(z), k);
    1836                 :            :     GEN a = RgM_RgC_mul(M, gel(T,i + nbE1));
    1837                 :            :     if (!gcmp0(a)) pari_err_BUG("checkrelation (T2)");
    1838                 :            :   }
    1839                 :            :   for (i = 1; i <= nbT31; i++)
    1840                 :            :   {
    1841                 :            :     GEN z = gel(annT31, i);
    1842                 :            :     GEN M = RgX_act_ZGl2Q(ZGl2Q_star(z), k);
    1843                 :            :     GEN a = RgM_RgC_mul(M, gel(T,i + nbE1 + nbT2));
    1844                 :            :     if (!gcmp0(a)) pari_err_BUG("checkrelation (T31)");
    1845                 :            :   }
    1846                 :            : }
    1847                 :            : /* phi_i(G_j) */
    1848                 :            : static GEN
    1849                 :            : eval_phii_Gj(GEN W, long i, long j)
    1850                 :            : {
    1851                 :            :   GEN basis = modsymbk_get_basis(W), b = gel(basis,i);
    1852                 :            :   GEN ind = gel(b,2), pols = gel(b,3);
    1853                 :            :   long s;
    1854                 :            :   for (s = 1; s < lg(ind); s++)
    1855                 :            :     if (ind[s] == j) return gel(pols,s);
    1856                 :            :   return zerocol(lg(gel(pols,1))-1);
    1857                 :            : }
    1858                 :            : /* check that \sum d_i phi_i(G_j)  = T_j for all j */
    1859                 :            : static void
    1860                 :            : checkdec(GEN W, GEN D, GEN T)
    1861                 :            : {
    1862                 :            :   long i, j;
    1863                 :            :   checkrelation(W,T);
    1864                 :            :   for (j = 1; j < lg(T); j++)
    1865                 :            :   {
    1866                 :            :     GEN S = gen_0;
    1867                 :            :     for (i = 1; i < lg(D); i++)
    1868                 :            :     {
    1869                 :            :       GEN d = gel(D,i);
    1870                 :            :       if (gcmp0(d)) continue;
    1871                 :            :       S = gadd(S, gmul(d, eval_phii_Gj(W, i, j)));
    1872                 :            :     }
    1873                 :            :     /* S = \sum_i d_i phi_i(G_j) */
    1874                 :            :     if (!gequal(S, gel(T,j)))
    1875                 :            :       pari_warn(warner, "checkdec j = %ld\n\tS = %Ps\n\tT = %Ps", j,S,gel(T,j));
    1876                 :            :   }
    1877                 :            : }
    1878                 :            : #endif
    1879                 :            : 
    1880                 :            : /* map op: Hom(W1,V) -> Hom(W2,V), given by \sum v[i], v[i] in Gl2(Q) */
    1881                 :            : static GEN
    1882                 :         70 : getMorphism(GEN W1, GEN W2, GEN v)
    1883                 :            : {
    1884                 :         70 :   pari_sp av0 = avma;
    1885                 :         70 :   GEN basis1 = modsymbk_get_basis(W1);
    1886                 :         70 :   long i, a, lv, dim1 = lg(basis1)-1;
    1887                 :         70 :   GEN M = cgetg(dim1+1, t_MAT), act;
    1888         [ +  + ]:         70 :   if (typ(v) != t_VEC) v = mkvec(v);
    1889                 :         70 :   lv = lg(v); act = cgetg(lv, t_MAT);
    1890         [ +  + ]:        203 :   for (i = 1; i < lv; i++) gel(act,i) = init_dual_act(gel(v,i), W1, W2);
    1891         [ +  + ]:        336 :   for (a = 1; a <= dim1; a++)
    1892                 :            :   {
    1893                 :        266 :     pari_sp av = avma;
    1894                 :        266 :     GEN phi = gel(basis1, a), D, T = NULL;
    1895         [ +  + ]:        777 :     for (i = 1; i < lv; i++)
    1896                 :            :     {
    1897                 :        511 :       GEN t = dual_act(phi, gel(act, i));
    1898         [ +  + ]:        511 :       T = T? gerepileupto(av, RgM_add(T,t)): t;
    1899                 :            :     }
    1900                 :            :     /* T = (phi|op)(G_1,...,G_d2) */
    1901                 :        266 :     D = getMorphism_basis(W2, T);
    1902                 :            : #if DEBUG
    1903                 :            :     checkdec(W2,D,T);
    1904                 :            : #endif
    1905                 :        266 :     gel(M,a) = gerepilecopy(av, D);
    1906                 :            :   }
    1907                 :         70 :   return gerepilecopy(av0, M);
    1908                 :            : }
    1909                 :            : 
    1910                 :            : static GEN
    1911                 :        224 : Tp_matrices(ulong p)
    1912                 :            : {
    1913                 :        224 :   GEN v = cgetg(p+2, t_VEC);
    1914                 :            :   ulong i;
    1915         [ +  + ]:       1351 :   for (i = 1; i <= p; i++) gel(v,i) = mat2(1, i-1, 0, p);
    1916                 :        224 :   gel(v,i) = mat2(p, 0, 0, 1);
    1917                 :        224 :   return v;
    1918                 :            : }
    1919                 :            : static GEN
    1920                 :         28 : Up_matrices(ulong p)
    1921                 :            : {
    1922                 :         28 :   GEN v = cgetg(p+1, t_VEC);
    1923                 :            :   ulong i;
    1924         [ +  + ]:        294 :   for (i = 1; i <= p; i++) gel(v,i) = mat2(1, i-1, 0, p);
    1925                 :         28 :   return v;
    1926                 :            : }
    1927                 :            : GEN
    1928                 :        252 : modsymbHecke(GEN W, ulong p)
    1929                 :            : {
    1930         [ +  + ]:        252 :   GEN v = modsymb_get_N(W) % p? Tp_matrices(p): Up_matrices(p);
    1931                 :        504 :   return modsymbk_get_weight(W) == 2? getMorphism_trivial(W,W,v)
    1932         [ +  + ]:        252 :                                     : getMorphism(W, W, v);
    1933                 :            : }
    1934                 :            : GEN
    1935                 :        161 : modsymbSigma(GEN W)
    1936                 :            : {
    1937                 :        161 :   GEN v = mat2(-1,0,0,1);
    1938                 :        322 :   return modsymbk_get_weight(W) == 2? getMorphism_trivial(W,W,v)
    1939         [ +  + ]:        161 :                                     : getMorphism(W, W, v);
    1940                 :            : }
    1941                 :            : 
    1942                 :            : #if 0
    1943                 :            : /* is \Gamma_0(N) cusp1 = \Gamma_0(N) cusp2 ? */
    1944                 :            : static int
    1945                 :            : iscuspeq(ulong N, GEN cusp1, GEN cusp2)
    1946                 :            : {
    1947                 :            :   long p1, q1, p2, q2, s1, s2, d;
    1948                 :            :   p1 = cusp1[1]; p2 = cusp2[1];
    1949                 :            :   q1 = cusp1[2]; q2 = cusp2[2];
    1950                 :            :   d = Fl_mul(smodss(q1,N),smodss(q2,N), N);
    1951                 :            :   d = ugcd(d, N);
    1952                 :            : 
    1953                 :            :   s1 = q1 > 2? Fl_inv(smodss(p1,q1), q1): 1;
    1954                 :            :   s2 = q2 > 2? Fl_inv(smodss(p2,q2), q2): 1;
    1955                 :            :   return Fl_mul(s1,q2,d) == Fl_mul(s2,q1,d);
    1956                 :            : }
    1957                 :            : #endif
    1958                 :            : 
    1959                 :            : static GEN
    1960                 :         35 : Cuspidal_subspace_trivial(GEN W0)
    1961                 :            : {
    1962                 :         35 :   GEN W = get_modsymb(W0);
    1963                 :         35 :   GEN section = modsymb_get_section(W), gen = modsymb_get_genindex(W);
    1964                 :         35 :   GEN S = modsymb_get_hashcusps(W);
    1965                 :         35 :   long j, nbE1 = modsymb_get_nbE1(W), ncusp = gel(S,1)[2];
    1966                 :         35 :   GEN T = zeromatcopy(ncusp,nbE1);
    1967         [ +  + ]:       2492 :   for (j = 1; j <= nbE1; j++)
    1968                 :            :   {
    1969                 :       2457 :     long e = gen[j]; /* path-index of E1-element */
    1970                 :       2457 :     GEN t = gel(section, e); /* path_to_matrix() */
    1971                 :       2457 :     long i1 = cusp_index(gel(t,1), S);
    1972                 :       2457 :     long i2 = cusp_index(gel(t,2), S);
    1973         [ +  + ]:       2457 :     if (i1 != i2)
    1974                 :            :     {
    1975                 :       1694 :       gcoeff(T, i1, j) = gen_1;
    1976                 :       1694 :       gcoeff(T, i2, j) = gen_m1;
    1977                 :            :     }
    1978                 :            :   }
    1979                 :         35 :   return ZM_ker(T);
    1980                 :            : }
    1981                 :            : 
    1982                 :            : /* return E_c(r) */
    1983                 :            : static GEN
    1984                 :          0 : get_Ec_r(GEN c, long k)
    1985                 :            : {
    1986                 :          0 :   long p = c[1], q = c[2], u, v;
    1987                 :            :   GEN gr;
    1988                 :          0 :   (void)cbezout(p, q, &u, &v);
    1989                 :          0 :   gr = mat2(p, -v, q, u); /* g . (1:0) = (p:q) */
    1990                 :          0 :   return voo_act_Gl2Q(zm_to_ZM(sl2_inv(gr)), k);
    1991                 :            : }
    1992                 :            : /* returns the modular symbol associated to the cusp c := p/q via the rule
    1993                 :            :  * E_c(path from a to b in Delta_0) := E_c(b) - E_c(a), where
    1994                 :            :  * E_c(r) := 0 if r != c mod Gamma
    1995                 :            :  *           v_oo | gamma_r^(-1)
    1996                 :            :  * where v_oo is stable by T = [1,1;0,1] (i.e x^(k-2)) and
    1997                 :            :  * gamma_r . (1:0) = r, for some gamma_r in SL_2(Z) * */
    1998                 :            : GEN
    1999                 :          0 : Eisenstein_symbol(GEN W, GEN c)
    2000                 :            : {
    2001                 :          0 :   GEN section = modsymb_get_section(W), gen = modsymb_get_genindex(W);
    2002                 :          0 :   GEN S = modsymb_get_hashcusps(W);
    2003                 :          0 :   long j, ic = cusp_index(c, S), l = lg(gen), k = modsymbk_get_weight(W);
    2004                 :          0 :   GEN vecT = cgetg(l, t_COL);
    2005         [ #  # ]:          0 :   for (j = 1; j < l; j++)
    2006                 :            :   {
    2007                 :          0 :     GEN vj = NULL, g = gel(section, gen[j]); /* path_to_matrix(generator) */
    2008                 :          0 :     GEN c1 = gel(g,1), c2 = gel(g,2);
    2009                 :          0 :     long i1 = cusp_index(c1, S);
    2010                 :          0 :     long i2 = cusp_index(c2, S);
    2011         [ #  # ]:          0 :     if (i1 == ic) vj = get_Ec_r(c1, k);
    2012         [ #  # ]:          0 :     if (i2 == ic)
    2013                 :            :     {
    2014                 :          0 :       GEN s = get_Ec_r(c2, k);
    2015         [ #  # ]:          0 :       vj = vj? gsub(vj, s): gneg(s);
    2016                 :            :     }
    2017         [ #  # ]:          0 :     if (!vj) vj = zerocol(k-1);
    2018                 :          0 :     gel(vecT, j) = vj;
    2019                 :            :   }
    2020                 :          0 :   return getMorphism_basis(W, vecT);
    2021                 :            : }
    2022                 :            : 
    2023                 :            : static GEN
    2024                 :         14 : EC_subspace_trivial(GEN W)
    2025                 :            : {
    2026                 :         14 :   GEN M, ch, chC, chE, T, TC, C = Cuspidal_subspace_trivial(W);
    2027                 :         14 :   long N = modsymb_get_N(W);
    2028                 :            :   ulong p;
    2029                 :            :   forprime_t S;
    2030                 :         14 :   (void)u_forprime_init(&S, 2, ULONG_MAX);
    2031         [ +  - ]:         14 :   while ((p = u_forprime_next(&S)))
    2032         [ +  - ]:         14 :     if (N % p) break;
    2033                 :         14 :   T = modsymbHecke(W, p);
    2034                 :         14 :   ch = QM_charpoly_ZX(T);
    2035                 :         14 :   TC = Qevproj_apply(T, Qevproj_init(C)); /* T_p | TC */
    2036                 :         14 :   chC = QM_charpoly_ZX(TC);
    2037                 :         14 :   chE = RgX_div(ch, chC); /* charpoly(T_p | E_k), coprime to chC */
    2038                 :         14 :   M = RgX_RgM_eval(chE, T);
    2039                 :         14 :   return mkvec2(Qevproj_Sigma(W, QM_ker(M)), C);
    2040                 :            : }
    2041                 :            : 
    2042                 :            : static GEN
    2043                 :          7 : Eisenstein_subspace_trivial(GEN W)
    2044                 :            : {
    2045                 :          7 :   GEN ES = EC_subspace_trivial(W);
    2046                 :          7 :   return gel(ES,1);
    2047                 :            : }
    2048                 :            : GEN
    2049                 :          7 : Eisenstein_subspace(GEN W)
    2050                 :            : {
    2051                 :            :   GEN S, cusps, M;
    2052                 :            :   long i, l;
    2053         [ +  - ]:          7 :   if (modsymbk_get_weight(W) == 2) return Eisenstein_subspace_trivial(W);
    2054                 :          0 :   S = modsymb_get_hashcusps(W);
    2055                 :          0 :   cusps = gel(S,3);
    2056                 :          0 :   l = lg(cusps);
    2057                 :          0 :   M = cgetg(l, t_MAT);
    2058         [ #  # ]:          0 :   for (i = 1; i < l; i++) gel(M,i) = Eisenstein_symbol(W, gel(cusps,i));
    2059                 :          7 :   return Qevproj_Sigma(W, QM_image(M));
    2060                 :            : }
    2061                 :            : 
    2062                 :            : GEN
    2063                 :          7 : Cuspidal_subspace(GEN W)
    2064                 :            : {
    2065                 :            :   GEN E;
    2066                 :            :   GEN M, T, TE, ch, chS, chE;
    2067                 :            :   forprime_t S;
    2068                 :            :   ulong p, N;
    2069                 :            : 
    2070         [ +  - ]:          7 :   if (modsymbk_get_weight(W) == 2) return EC_subspace_trivial(W);
    2071                 :          0 :   E = Eisenstein_subspace(W);
    2072                 :          0 :   N = modsymb_get_N(W);
    2073                 :          0 :   (void)u_forprime_init(&S, 2, ULONG_MAX);
    2074         [ #  # ]:          0 :   while ((p = u_forprime_next(&S)))
    2075         [ #  # ]:          0 :     if (N % p) break;
    2076                 :          0 :   T = modsymbHecke(W, p);
    2077                 :          0 :   ch = QM_charpoly_ZX(T);
    2078                 :          0 :   TE = Qevproj_apply(T, Qevproj_init(E)); /* T_p | E */
    2079                 :          0 :   chE = QM_charpoly_ZX(TE);
    2080                 :          0 :   chS = RgX_div(ch, chE); /* charpoly(T_p | S_k), coprime to chE */
    2081                 :          0 :   M = RgX_RgM_eval(chS, T);
    2082                 :          7 :   return mkvec2(E, Qevproj_Sigma(W, QM_ker(M)));
    2083                 :            : }
    2084                 :            : 
    2085                 :            : /** INIT ELLSYM STRUCTURE **/
    2086                 :            : /* V a vector of ZM. If all of them have 0 last row, return NULL.
    2087                 :            :  * Otherwise return [m,i,j], where m = V[i][last,j] contains the value
    2088                 :            :  * of smallest absolute value */
    2089                 :            : static GEN
    2090                 :         49 : RgMV_find_non_zero_last_row(long offset, GEN V)
    2091                 :            : {
    2092                 :         49 :   long i, lasti = 0, lastj = 0, lV = lg(V);
    2093                 :         49 :   GEN m = NULL;
    2094         [ +  + ]:        105 :   for (i = 1; i < lV; i++)
    2095                 :            :   {
    2096                 :         56 :     GEN M = gel(V,i);
    2097                 :         56 :     long j, n, l = lg(M);
    2098         [ +  + ]:         56 :     if (l == 1) continue;
    2099                 :         42 :     n = nbrows(M);
    2100         [ +  + ]:        175 :     for (j = 1; j < l; j++)
    2101                 :            :     {
    2102                 :        133 :       GEN a = gcoeff(M, n, j);
    2103 [ +  + ][ +  + ]:        133 :       if (!gcmp0(a) && (!m || absi_cmp(a, m) < 0))
                 [ +  + ]
    2104                 :            :       {
    2105                 :         56 :         m = a; lasti = i; lastj = j;
    2106         [ -  + ]:         56 :         if (is_pm1(m)) goto END;
    2107                 :            :       }
    2108                 :            :     }
    2109                 :            :   }
    2110                 :            : END:
    2111         [ +  + ]:         49 :   if (!m) return NULL;
    2112                 :         49 :   return mkvec2(m, mkvecsmall2(lasti+offset, lastj));
    2113                 :            : }
    2114                 :            : /* invert the d_oo := (\gamma_oo - 1) operator, acting on
    2115                 :            :  * [x^(k-2), ..., y^(k-2)] */
    2116                 :            : static GEN
    2117                 :         35 : Delta_inv(GEN doo, long k)
    2118                 :            : {
    2119                 :         35 :   GEN M = RgX_act_ZGl2Q(doo, k);
    2120                 :         35 :   M = RgM_minor(M, k-1, 1); /* 1st column and last row are 0 */
    2121                 :         35 :   return ZM_inv_denom(M);
    2122                 :            : }
    2123                 :            : /* The ZX P = \sum a_i x^i y^{k-2-i} is given by the ZV [a_0, ..., a_k-2]~,
    2124                 :            :  * return Q and d such that P = doo Q + d y^k-2, where d in Z and Q */
    2125                 :            : static GEN
    2126                 :        133 : doo_decompose(GEN dinv, GEN P, GEN *pd)
    2127                 :            : {
    2128                 :        133 :   long l = lg(P); *pd = gel(P, l-1);
    2129                 :        133 :   P = vecslice(P, 1, l-2);
    2130                 :        133 :   return shallowconcat(gen_0, ZC_apply_dinv(dinv, P));
    2131                 :            : }
    2132                 :            : 
    2133                 :            : static GEN
    2134                 :        133 : get_phi_ij(long i,long j,long n, long s,long t,GEN P_st,GEN Q_st,GEN d_st,
    2135                 :            :            GEN P_ij, GEN lP_ij, GEN dinv)
    2136                 :            : {
    2137                 :            :   GEN ind, pols;
    2138 [ +  + ][ +  + ]:        133 :   if (i == s && j == t)
    2139                 :            :   {
    2140                 :         35 :     ind = mkvecsmall(1);
    2141                 :         35 :     pols = mkvec(scalarcol_shallow(gen_1, lg(P_st)-1)); /* x^{k-2} */
    2142                 :            :   }
    2143                 :            :   else
    2144                 :            :   {
    2145                 :         98 :     GEN d_ij, Q_ij = doo_decompose(dinv, lP_ij, &d_ij);
    2146                 :         98 :     GEN a = ZC_Z_mul(P_ij, d_st);
    2147                 :         98 :     GEN b = ZC_Z_mul(P_st, negi(d_ij));
    2148                 :         98 :     GEN c = RgC_sub(RgC_Rg_mul(Q_ij, d_st), RgC_Rg_mul(Q_st, d_ij));
    2149         [ +  + ]:         98 :     if (i == s) { /* j != t */
    2150                 :         77 :       ind = mkvecsmall2(1, s);
    2151                 :         77 :       pols = mkvec2(c, ZC_add(a, b));
    2152                 :            :     } else {
    2153                 :         21 :       ind = mkvecsmall3(1, i, s);
    2154                 :         98 :       pols = mkvec3(c, a, b); /* image of g_1, g_i, g_s */
    2155                 :            :     }
    2156                 :            :   }
    2157                 :        133 :   return mkvec3(mkvecsmall3(i,j,n), ind, pols);
    2158                 :            : }
    2159                 :            : 
    2160                 :            : static GEN
    2161                 :        161 : add_Sigma(GEN W, long sign)
    2162                 :            : {
    2163                 :        161 :   gel(W, 2) = mkvec2(stoi(sign), modsymbSigma(W));
    2164                 :        161 :   return W;
    2165                 :            : }
    2166                 :            : static GEN
    2167                 :        126 : modsymbkinit_trivial(GEN WN, long sign)
    2168                 :            : {
    2169                 :        126 :   long dim = modsymb_get_nbE1(WN);
    2170                 :        126 :   GEN W = mkvec3(WN, gen_0, mkvec2(gen_0,mkvecsmall2(2, dim)));
    2171                 :        126 :   return add_Sigma(W, sign);
    2172                 :            : }
    2173                 :            : /* sum of #cols of the matrices contained in V */
    2174                 :            : static long
    2175                 :         70 : RgMV_dim(GEN V)
    2176                 :            : {
    2177                 :         70 :   long l = lg(V), d = 0, i;
    2178         [ +  + ]:         91 :   for (i = 1; i < l; i++) d += lg(gel(V,i)) - 1;
    2179                 :         70 :   return d;
    2180                 :            : }
    2181                 :            : static GEN
    2182                 :         35 : modsymbkinit_nontrivial(GEN WN, long k, long sign)
    2183                 :            : {
    2184                 :         35 :   GEN annT2 = gel(WN,8), annT31 = gel(WN,9), singlerel = gel(WN,10);
    2185                 :            :   GEN link, basis, monomials, invphiblock;
    2186                 :         35 :   long nbE1 = modsymb_get_nbE1(WN);
    2187                 :         35 :   GEN dinv = Delta_inv(ZG_neg( ZGl2Q_star(gel(singlerel,1)) ), k);
    2188                 :         35 :   GEN p1 = cgetg(nbE1+1, t_VEC), remove;
    2189                 :         35 :   GEN p2 = ZGV_tors(annT2, k);
    2190                 :         35 :   GEN p3 = ZGV_tors(annT31, k);
    2191                 :         35 :   GEN gentor = shallowconcat(p2, p3);
    2192                 :            :   GEN P_st, lP_st, Q_st, d_st, W;
    2193                 :            :   long n, i, dim, s, t, u;
    2194                 :         35 :   gel(p1, 1) = cgetg(1,t_MAT); /* dummy */
    2195         [ +  + ]:         56 :   for (i = 2; i <= nbE1; i++) /* skip 1st element = (\gamma_oo-1)g_oo */
    2196                 :            :   {
    2197                 :         21 :     GEN z = gel(singlerel, i);
    2198                 :         21 :     gel(p1, i) = RgX_act_ZGl2Q(ZGl2Q_star(z), k);
    2199                 :            :   }
    2200                 :         35 :   remove = RgMV_find_non_zero_last_row(nbE1, gentor);
    2201         [ +  + ]:         35 :   if (!remove) remove = RgMV_find_non_zero_last_row(0, p1);
    2202         [ -  + ]:         35 :   if (!remove) pari_err_BUG("msinit [no y^k-2]");
    2203                 :         35 :   remove = gel(remove,2); /* [s,t] */
    2204                 :         35 :   s = remove[1];
    2205                 :         35 :   t = remove[2];
    2206                 :            :   /* +1 because of = x^(k-2), but -1 because of Manin relation */
    2207                 :         35 :   dim = (k-1)*(nbE1-1) + RgMV_dim(p2) + RgMV_dim(p3);
    2208                 :            :   /* Let (g_1,...,g_d) be the Gamma-generators of Delta, g_1 = g_oo.
    2209                 :            :    * We describe modular symbols by the collection phi(g_1), ..., phi(g_d)
    2210                 :            :    * \in V := Q[x,y]_{k-2}, with right Gamma action.
    2211                 :            :    * For each i = 1, .., d, let V_i \subset V be the Q-vector space of
    2212                 :            :    * allowed values for phi(g_i): with basis (P^{i,j}) given by the monomials
    2213                 :            :    * x^(j-1) y^{k-2-(j-1)}, j = 1 .. k-1
    2214                 :            :    * (g_i in E_1) or the solution of the torsion equations (1 + gamma)P = 0
    2215                 :            :    * (g_i in T2) or (1 + gamma + gamma^2)P = 0 (g_i in T31).
    2216                 :            :    *
    2217                 :            :    * The Manin relation (singlerel) is of the form \sum_i \lambda_i g_i = 0,
    2218                 :            :    * where \lambda_i = 1 if g_i in T2 or T31, and \lambda_i = (1 - \gamma_i)
    2219                 :            :    * for g_i in E1.
    2220                 :            :    *
    2221                 :            :    * If phi \in Hom_Gamma(Delta, V), it is defined by phi(g_i) := P_i in V
    2222                 :            :    * with \sum_i P_i . \lambda_i^*, where (\sum n_i g_i)^* :=
    2223                 :            :    * \sum n_i \gamma_i^(-1).
    2224                 :            :    *
    2225                 :            :    * We single out gamma_1 / g_1 (g_oo in Pollack-Stevens paper) and
    2226                 :            :    * write P_{i,j} \lambda_i^* =  Q_{i,j} (\gamma_1 - 1)^* + d_{i,j} y^{k-2}
    2227                 :            :    * where d_{i,j} is a scalar and Q_{i,j} in V; we normalize Q_{i,j} to
    2228                 :            :    * that the coefficient of x^{k-2} is 0.
    2229                 :            :    *
    2230                 :            :    * There exist (s,t) such that d_{s,t} != 0.
    2231                 :            :    * A Q-basis of the (dual) space of modular symbols is given by the
    2232                 :            :    * functions phi_{i,j}, 2 <= i <= d, 1 <= j <= k-1, mapping
    2233                 :            :    *  g_1 -> d_{s,t} Q_{i,j} - d_{i,j} Q_{s,t} + [(i,j)=(s,t)] x^{k-2}
    2234                 :            :    * If i != s
    2235                 :            :    *   g_i -> d_{s,t} P_{i,j}
    2236                 :            :    *   g_s -> - d_{i,j} P_{s,t}
    2237                 :            :    * If i = s, j != t
    2238                 :            :    *   g_i -> d_{s,t} P_{i,j} - d_{i,j} P_{s,t}
    2239                 :            :    * And everything else to 0. */
    2240                 :         35 :   monomials = idmat(k-1); /* represent the monomials x^{k-2}, ... , y^{k-2} */
    2241         [ +  + ]:         35 :   if (s <= nbE1) /* in E1 */
    2242                 :            :   {
    2243                 :         14 :     P_st = gel(monomials, t);
    2244                 :         14 :     lP_st = gmael(p1, s, t); /* P_{s,t} lambda_s^* */
    2245                 :            :   }
    2246                 :            :   else /* in T2, T31 */
    2247                 :            :   {
    2248                 :         21 :     P_st = gmael(gentor, s - nbE1, t);
    2249                 :         21 :     lP_st = P_st;
    2250                 :            :   }
    2251                 :         35 :   Q_st = doo_decompose(dinv, lP_st, &d_st);
    2252                 :         35 :   basis = cgetg(dim+1, t_VEC);
    2253                 :         35 :   link = cgetg(nbE1 + lg(gentor), t_VEC);
    2254                 :         35 :   gel(link,1) = cgetg(1,t_VECSMALL); /* dummy */
    2255                 :         35 :   n = 1;
    2256         [ +  + ]:         56 :   for (i = 2; i <= nbE1; i++)
    2257                 :            :   {
    2258                 :         21 :     GEN L = cgetg(k, t_VECSMALL);
    2259                 :            :     long j;
    2260                 :            :     /* link[i][j] = n gives correspondance between phi_{i,j} and basis[n] */
    2261                 :         21 :     gel(link,i) = L;
    2262         [ +  + ]:         84 :     for (j = 1; j < k; j++)
    2263                 :            :     {
    2264                 :         63 :       GEN lP_ij = gmael(p1, i, j); /* P_{i,j} lambda_i^* */
    2265                 :         63 :       GEN P_ij = gel(monomials,j);
    2266                 :         63 :       L[j] = n;
    2267                 :         63 :       gel(basis, n) = get_phi_ij(i,j,n, s,t, P_st, Q_st, d_st, P_ij, lP_ij, dinv);
    2268                 :         63 :       n++;
    2269                 :            :     }
    2270                 :            :   }
    2271         [ +  + ]:         56 :   for (u = 1; u < lg(gentor); u++,i++)
    2272                 :            :   {
    2273                 :         21 :     GEN V = gel(gentor,u);
    2274                 :         21 :     long j, lV = lg(V);
    2275                 :         21 :     GEN L = cgetg(lV, t_VECSMALL);
    2276                 :         21 :     gel(link,i) = L;
    2277         [ +  + ]:         91 :     for (j = 1; j < lV; j++)
    2278                 :            :     {
    2279                 :         70 :       GEN lP_ij = gel(V, j); /* P_{i,j} lambda_i^* = P_{i,j} */
    2280                 :         70 :       GEN P_ij = lP_ij;
    2281                 :         70 :       L[j] = n;
    2282                 :         70 :       gel(basis, n) = get_phi_ij(i,j,n, s,t, P_st, Q_st, d_st, P_ij, lP_ij, dinv);
    2283                 :         70 :       n++;
    2284                 :            :     }
    2285                 :            :   }
    2286                 :         35 :   invphiblock = cgetg(lg(link), t_VEC);
    2287                 :         35 :   gel(invphiblock,1) = cgetg(1, t_MAT); /* dummy */
    2288         [ +  + ]:         77 :   for (i = 2; i < lg(link); i++)
    2289                 :            :   {
    2290                 :         42 :     GEN M, inv, B = gel(link,i);
    2291                 :         42 :     long j, lB = lg(B);
    2292         [ +  + ]:         42 :     if (i == s) { B = vecsplice(B, t); lB--; } /* remove phi_st */
    2293                 :         42 :     M = cgetg(lB, t_MAT);
    2294         [ +  + ]:        140 :     for (j = 1; j < lB; j++)
    2295                 :            :     {
    2296                 :         98 :       GEN phi_ij = gel(basis, B[j]), pols = gel(phi_ij,3);
    2297                 :         98 :       gel(M, j) = gel(pols, 2); /* phi_ij(g_i) */
    2298                 :            :     }
    2299 [ +  + ][ +  + ]:         42 :     if (i <= nbE1 && i != s) /* maximal rank k-1 */
    2300                 :          7 :       inv = ZM_inv_denom(M);
    2301                 :            :     else /* i = s (rank k-2) or from torsion: rank k/3 or k/2 */
    2302                 :         35 :       inv = Qevproj_init(M);
    2303                 :         42 :     gel(invphiblock,i) = inv;
    2304                 :            :   }
    2305                 :         35 :   W = mkvec3(WN, gen_0, mkvec5(basis, mkvecsmall2(k, dim), mkvecsmall2(s,t),
    2306                 :            :                                link, invphiblock));
    2307                 :         35 :   return add_Sigma(W, sign);
    2308                 :            : }
    2309                 :            : /* WN = modsymbinit_N(N) */
    2310                 :            : static GEN
    2311                 :        161 : modsymbkinit(ulong N, long k, long sign)
    2312                 :            : {
    2313                 :        161 :   GEN WN = modsymbinit_N(N);
    2314                 :        161 :   return k == 2? modsymbkinit_trivial(WN, sign)
    2315         [ +  + ]:        161 :                : modsymbkinit_nontrivial(WN, k, sign);
    2316                 :            : }
    2317                 :            : GEN
    2318                 :         84 : msinit(GEN N, GEN K, long sign)
    2319                 :            : {
    2320                 :         84 :   pari_sp av = avma;
    2321                 :            :   GEN W;
    2322                 :            :   long k;
    2323         [ -  + ]:         84 :   if (typ(N) != t_INT) pari_err_TYPE("msinit", N);
    2324         [ -  + ]:         84 :   if (typ(K) != t_INT) pari_err_TYPE("msinit", K);
    2325                 :         84 :   k = itos(K);
    2326         [ -  + ]:         84 :   if (k < 2) pari_err_DOMAIN("msinit","k", "<", gen_2,K);
    2327         [ -  + ]:         84 :   if (odd(k)) pari_err_IMPL("msinit [odd weight]");
    2328         [ -  + ]:         84 :   if (signe(N) <= 0) pari_err_DOMAIN("msinit","N", "<=", gen_0,N);
    2329         [ -  + ]:         84 :   if (equali1(N)) pari_err_IMPL("msinit [ N = 1 ]");
    2330                 :         84 :   W = modsymbkinit(itou(N), k, sign);
    2331                 :         84 :   return gerepilecopy(av, W);
    2332                 :            : }
    2333                 :            : 
    2334                 :            : /* c t_FRAC */
    2335                 :            : GEN
    2336                 :    1916656 : Q_xpm(GEN E, GEN c)
    2337                 :            : {
    2338                 :    1916656 :   pari_sp av = avma;
    2339                 :    1916656 :   GEN dualell = ellsym_get_dualell(E), W = ellsym_get_W(E);
    2340                 :    1916656 :   GEN v = init_act_trivial(W);
    2341                 :    1916656 :   GEN oo_0 = gmael(W,15,1);
    2342                 :    1916656 :   long index = gel(oo_0,1)[1];
    2343                 :            : 
    2344                 :    1916656 :   cusplog_trivial_frac(v, W, 0, c); /* oo -> (a:b), c = a/b */
    2345                 :    1916656 :   treat_index_trivial(W, index, 1, v); /* - (oo->0) to correct to 0 -> (a:b) */
    2346                 :    1916656 :   return gerepileuptoint(av, ZV_dotproduct(dualell,v));
    2347                 :            : }
    2348                 :            : 
    2349                 :            : /* <0->oo>: unused */
    2350                 :            : static GEN
    2351                 :          0 : xpmoo(GEN E)
    2352                 :            : {
    2353                 :          0 :   pari_sp av = avma;
    2354                 :          0 :   GEN dualell = ellsym_get_dualell(E), W = ellsym_get_W(E);
    2355                 :          0 :   GEN v = init_act_trivial(W);
    2356                 :          0 :   GEN oo_0 = gmael(W,15,1);
    2357                 :          0 :   long index = gel(oo_0,1)[1];
    2358                 :          0 :   treat_index_trivial(W, index, 1, v); /* - (oo->0) to correct to 0 -> (a:b) */
    2359                 :          0 :   return gerepileuptoint(av, ZV_dotproduct(dualell,v));
    2360                 :            : }
    2361                 :            : 
    2362                 :            : /* E = ellsym struct; image of <0->a/b> */
    2363                 :            : GEN
    2364                 :          0 : xpm(GEN E, GEN a, GEN b)
    2365         [ #  # ]:          0 : { return signe(b)? Q_xpm(E, gdiv(a,b)): xpmoo(E); }
    2366                 :            : 
    2367                 :            : static GEN
    2368                 :         63 : twistcurve(GEN e, GEN D)
    2369                 :            : {
    2370                 :         63 :   GEN D2 = sqri(D);
    2371                 :         63 :   GEN a4 = mulii(mulsi(-27, D2), ell_get_c4(e));
    2372                 :         63 :   GEN a6 = mulii(mulsi(-54, mulii(D, D2)), ell_get_c6(e));
    2373                 :         63 :   return ellinit(mkvec2(a4,a6),NULL,0);
    2374                 :            : }
    2375                 :            : 
    2376                 :            : /* sum_{a < |D|} (D/a)*xpm(E,a/D) */
    2377                 :            : static GEN
    2378                 :        154 : get_X(GEN E, long D)
    2379                 :            : {
    2380                 :        154 :   ulong a, d = (ulong)labs(D);
    2381                 :        154 :   GEN t = gen_0;
    2382                 :        154 :   GEN nc = icopy(gen_1), c = mkfrac(nc, utoipos(d));
    2383         [ +  + ]:        630 :   for (a=1; a < d; a++)
    2384                 :            :   {
    2385                 :        476 :     long s = kross(D,a);
    2386                 :            :     GEN x;
    2387         [ +  + ]:        476 :     if (!s) continue;
    2388                 :        364 :     nc[2] = a; x = Q_xpm(E, c);
    2389         [ +  + ]:        364 :     t = (s > 0)? addii(t, x): subii(t, x);
    2390                 :            :   }
    2391                 :        154 :   return t;
    2392                 :            : }
    2393                 :            : /* quotient of the Neron periods of E^(d) and E, divided by sqrt(d) */
    2394                 :            : static long
    2395                 :         63 : get_alpha_d(GEN E, long d)
    2396                 :            : {
    2397         [ +  + ]:         63 :   if (odd(d)) return 1;
    2398         [ +  + ]:         21 :   if (!mpodd(ell_get_c4(E))) return 2; /* additive reduction at 2 */
    2399                 :            :   /* reduction not additive */
    2400 [ +  - ][ -  + ]:         63 :   return (d % 8 == 0 && !mpodd(ell_get_a1(E)))? 2: 1;
    2401                 :            : }
    2402                 :            : /* write L(E,1) = Q*w1, return the rational Q */
    2403                 :            : static GEN
    2404                 :         63 : get_Q(GEN E)
    2405                 :            : {
    2406                 :            :   GEN L, N, tam, T, n, w1;
    2407                 :            :   long ex, t, t2;
    2408                 :         63 :   E = ellanal_globalred_all(E, NULL, &N, &tam);
    2409                 :         63 :   T = elltors(E); t = itos(gel(T,1)); t2 = t*t;
    2410                 :         63 :   w1 = gel(ellR_omega(E,DEFAULTPREC), 1);
    2411                 :            : 
    2412                 :            :   /* |Sha| = n^2 */
    2413                 :         63 :   L = ellL1(E, 0, DEFAULTPREC);
    2414                 :         63 :   n = sqrtr(divrr(mulru(L, t2), mulri(w1,tam)));
    2415                 :         63 :   n = grndtoi(n, &ex);
    2416         [ -  + ]:         63 :   if (ex > -5) pari_err_BUG("ellsym (can't compute analytic |Sha|)");
    2417                 :         63 :   return gdivgs(mulii(tam,sqri(n)), t2);
    2418                 :            : }
    2419                 :            : 
    2420                 :            : /* return C such that C*L(E,1)_{xpm} = L(E,1) / w1 */
    2421                 :            : static GEN
    2422                 :         63 : ell_get_scale(GEN E, long s)
    2423                 :            : {
    2424                 :         63 :   pari_sp av = avma;
    2425                 :         63 :   GEN Q, X = NULL, e;
    2426                 :            :   long d;
    2427                 :            : 
    2428                 :            :   /* find D = s*d such that twist by D has rank 0 */
    2429         [ +  - ]:        406 :   for (d = 1; d < LONG_MAX; d++,avma=av)
    2430                 :            :   {
    2431         [ -  + ]:        406 :     if (s < 0)
    2432         [ #  # ]:          0 :     { if (!unegisfundamental(d)) continue; }
    2433                 :            :     else
    2434         [ +  + ]:        406 :     { if (!uposisfundamental(d)) continue; }
    2435         [ -  + ]:        154 :     X = get_X(E, s < 0? -d: d);
    2436         [ +  + ]:        154 :     if (signe(X)) break;
    2437                 :            :   }
    2438         [ -  + ]:         63 :   if (d == LONG_MAX) pari_err_BUG("ellsym (no suitable twist)");
    2439         [ -  + ]:         63 :   if (s < 0) d = -d;
    2440                 :         63 :   e = ellsym_get_ell(E);
    2441                 :         63 :   Q = get_Q(twistcurve(e, stoi(d)));
    2442                 :         63 :   return gerepileupto(av, gdiv(gmulsg(get_alpha_d(e,d), Q), X));
    2443                 :            : }
    2444                 :            : 
    2445                 :            : GEN
    2446                 :         63 : ellsym(GEN E, long sign)
    2447                 :            : {
    2448                 :         63 :   pari_sp av = avma;
    2449                 :            :   GEN cond;
    2450                 :            :   GEN W, K, dualell, scale_L;
    2451                 :            :   ulong p, N;
    2452                 :            :   forprime_t T;
    2453                 :            : 
    2454                 :         63 :   E = ellminimalmodel(E, NULL);
    2455                 :         63 :   cond = gel(ellglobalred(E), 1);
    2456                 :         63 :   N = itou(cond);
    2457                 :         63 :   W = modsymbkinit(N, 2, sign);
    2458                 :            : 
    2459         [ -  + ]:         63 :   if (!sign) pari_err(e_MISC, "missing sign in ellsym");
    2460                 :         63 :   K = keri(shallowtrans(gsubgs(modsymbk_get_Sigma(W), sign)));
    2461                 :            :   /* loop for p <= count_Manin_symbols(N) / 6 would be enough */
    2462                 :         63 :   (void)u_forprime_init(&T, 2, ULONG_MAX);
    2463         [ +  - ]:         91 :   while( (p = u_forprime_next(&T)) )
    2464                 :            :   {
    2465                 :            :     GEN Tp, ap, M, K2;
    2466         [ +  + ]:         91 :     if (N % p == 0) continue;
    2467                 :         70 :     Tp = modsymbHecke(W, p);
    2468                 :         70 :     ap = ellap(E, utoipos(p));
    2469                 :         70 :     M = RgM_Rg_add_shallow(Tp, negi(ap));
    2470                 :         70 :     K2 = keri( ZM_mul(shallowtrans(M), K) );
    2471         [ +  - ]:         70 :     if (lg(K2) < lg(K)) K = ZM_mul(K, K2);
    2472         [ +  + ]:         70 :     if (lg(K2)-1 == 1) break;
    2473                 :            :   }
    2474         [ -  + ]:         63 :   if (!p) pari_err_BUG("ellsym: ran out of primes");
    2475                 :            :   /* linear form = 0 on all Im(Tp - ap) and Im(S - sign) */
    2476                 :         63 :   dualell = Q_primpart(gel(K,1)); settyp(dualell, t_VEC);
    2477                 :         63 :   scale_L = ell_get_scale(mkvec3(W, dualell, E), sign);
    2478                 :         63 :   return gerepilecopy(av, mkvec4(W, dualell, E, scale_L));
    2479                 :            : }

Generated by: LCOV version 1.9