/* Copyright (C) 2011 The PARI group. This file is part of the PARI/GP package. PARI/GP is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation. It is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY WHATSOEVER. Check the License for details. You should have received a copy of it, along with the package; see the file 'COPYING'. If not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ #include "pari.h" #include "paripriv.h" /* Adapted from shp_package/moments by Robert Pollack * http://www.math.mcgill.ca/darmon/programs/shp/shp.html */ static GEN mskinit(ulong N, long k, long sign); static GEN mshecke_i(GEN W, ulong p); static GEN ZSl2_star(GEN v); static GEN getMorphism(GEN W1, GEN W2, GEN v); static GEN voo_act_Gl2Q(GEN g, long k); /* Input: P^1(Z/NZ) (formed by create_p1mod) Output: # P^1(Z/NZ) */ static long p1_size(GEN p1N) { return lg(gel(p1N,1)) - 1; } static ulong p1N_get_N(GEN p1N) { return gel(p1N,3)[2]; } static GEN p1N_get_hash(GEN p1N) { return gel(p1N,2); } static GEN p1N_get_fa(GEN p1N) { return gel(p1N,4); } static GEN p1N_get_div(GEN p1N) { return gel(p1N,5); } static GEN p1N_get_invsafe(GEN p1N) { return gel(p1N,6); } static GEN p1N_get_inverse(GEN p1N) { return gel(p1N,7); } /* ms-specific accessors */ /* W = msinit, return the output of msinit_N */ static GEN get_msN(GEN W) { return lg(W) == 4? gel(W,1): W; } static GEN msN_get_p1N(GEN W) { return gel(W,1); } static GEN msN_get_genindex(GEN W) { return gel(W,5); } static GEN msN_get_E2fromE1(GEN W) { return gel(W,7); } static GEN msN_get_annT2(GEN W) { return gel(W,8); } static GEN msN_get_annT31(GEN W) { return gel(W,9); } static GEN msN_get_singlerel(GEN W) { return gel(W,10); } static GEN msN_get_section(GEN W) { return gel(W,12); } static GEN ms_get_p1N(GEN W) { return msN_get_p1N(get_msN(W)); } static long ms_get_N(GEN W) { return p1N_get_N(ms_get_p1N(W)); } static GEN ms_get_hashcusps(GEN W) { W = get_msN(W); return gel(W,16); } static GEN ms_get_section(GEN W) { return msN_get_section(get_msN(W)); } static GEN ms_get_genindex(GEN W) { return msN_get_genindex(get_msN(W)); } static long ms_get_nbgen(GEN W) { return lg(ms_get_genindex(W))-1; } static long ms_get_nbE1(GEN W) { GEN W11; W = get_msN(W); W11 = gel(W,11); return W11[4] - W11[3]; } /* msk-specific accessors */ static long msk_get_dim(GEN W) { return gmael(W,3,2)[2]; } static GEN msk_get_basis(GEN W) { return gmael(W,3,1); } static long msk_get_weight(GEN W) { return gmael(W,3,2)[1]; } static long msk_get_sign(GEN W) { GEN t = gel(W,2); return typ(t)==t_INT? 0: itos(gel(t,1)); } static GEN msk_get_star(GEN W) { return gmael(W,2,2); } static GEN msk_get_starproj(GEN W) { return gmael(W,2,3); } static int is_Qevproj(GEN x) { return typ(x) == t_VEC && lg(x) == 5 && typ(gel(x,1)) == t_MAT; } long msdim(GEN W) { if (is_Qevproj(W)) return lg(gel(W,1)) - 1; checkms(W); if (!msk_get_sign(W)) return msk_get_dim(W); return lg(gel(msk_get_starproj(W), 1)) - 1; } long msgetlevel(GEN W) { checkms(W); return ms_get_N(W); } long msgetweight(GEN W) { checkms(W); return msk_get_weight(W); } long msgetsign(GEN W) { checkms(W); return msk_get_sign(W); } void checkms(GEN W) { if (typ(W) != t_VEC || lg(W) != 4) pari_err_TYPE("checkms [please apply msinit]", W); } /** MODULAR TO SYM **/ /* q a t_FRAC or t_INT */ static GEN Q_log_init(ulong N, GEN q) { long l, n; GEN Q; q = gboundcf(q, 0); l = lg(q); Q = cgetg(l, t_VECSMALL); Q[1] = 1; for (n=2; n 1, gcd(c,d)=1. * Furthermore, d != d0 (mod N/c) with c,d0 already in the list */ l = lg(divN) - 1; /* c = 1 first */ gel(ret, ++curn) = mkvecsmall2(1,0); for (d = 2; d < N; d++) if (ugcd(d,N) != 1UL) gel(ret, ++curn) = mkvecsmall2(1,d); /* omit c = 1 (first) and c = N (last) */ for (i=2; i < l; i++) { ulong Novc, d0; c = divN[i]; Novc = N / c; for (d0 = 2; d0 <= Novc; d0++) { ulong k, d = d0; if (ugcd(d, Novc) == 1UL) continue; for (k = 0; k < c; k++, d += Novc) if (ugcd(c,d) == 1UL) { gel(ret, ++curn) = mkvecsmall2(c,d); break; } } } if (curn != num) pari_err_BUG("generatemsymbols [wrong number of symbols]"); return ret; } static GEN inithashmsymbols(ulong N, GEN symbols) { GEN H = zerovec(N); long k, l = lg(symbols); /* skip the (c:1), 0 <= c < N and (1:0) */ for (k=N+2; k < l; k++) { GEN s = gel(symbols, k); ulong c = s[1], d = s[2], Novc = N/c; if (gel(H,c) == gen_0) gel(H,c) = const_vecsmall(Novc+1,0); if (c != 1) { d %= Novc; if (!d) d = Novc; } mael(H, c, d) = k; } return H; } /** Helper functions for Sl2(Z) / Gamma_0(N) **/ /* M a 2x2 ZM in SL2(Z) */ static GEN SL2_inv(GEN M) { GEN a = gcoeff(M,1,1), b = gcoeff(M,1,2); GEN c = gcoeff(M,2,1), d = gcoeff(M,2,2); return mkmat22(d,negi(b), negi(c),a); } /* SL2_inv(M)[2] */ static GEN SL2_inv2(GEN M) { GEN a = gcoeff(M,1,1), b = gcoeff(M,1,2); return mkcol2(negi(b),a); } /* M a 2x2 mat2 in SL2(Z) */ static GEN sl2_inv(GEN M) { long a=coeff(M,1,1), b=coeff(M,1,2), c=coeff(M,2,1), d=coeff(M,2,2); return mkvec2(mkvecsmall2(d, -c), mkvecsmall2(-b, a)); } /* Return the mat2 [a,b; c,d], not a zm to avoid GP problems */ static GEN mat2(long a, long b, long c, long d) { return mkvec2(mkvecsmall2(a,c), mkvecsmall2(b,d)); } static GEN mat2_to_ZM(GEN M) { GEN A = gel(M,1), B = gel(M,2); retmkmat2(mkcol2s(A[1],A[2]), mkcol2s(B[1],B[2])); } /* Input: a = 2-vector = path = {r/s,x/y} * Output: either [r,x;s,y] or [-r,x;-s,y], whichever has determinant > 0 */ static GEN path_to_ZM(GEN a) { GEN v = gel(a,1), w = gel(a,2); long r = v[1], s = v[2], x = w[1], y = w[2]; if (cmpii(mulss(r,y), mulss(x,s)) < 0) { r = -r; s = -s; } return mkmat22(stoi(r),stoi(x),stoi(s),stoi(y)); } static GEN path_to_zm(GEN a) { GEN v = gel(a,1), w = gel(a,2); long r = v[1], s = v[2], x = w[1], y = w[2]; if (cmpii(mulss(r,y), mulss(x,s)) < 0) { r = -r; s = -s; } return mat2(r,x,s,y); } /* path from c1 to c2 */ static GEN mkpath(GEN c1, GEN c2) { return mat2(c1[1], c2[1], c1[2], c2[2]); } static long cc(GEN M) { GEN v = gel(M,1); return v[2]; } static long dd(GEN M) { GEN v = gel(M,2); return v[2]; } /*Input: a,b = 2 paths, N = integer *Output: 1 if the a,b are \Gamma_0(N)-equivalent; 0 otherwise */ static int gamma_equiv(GEN a, GEN b, ulong N) { pari_sp av = avma; GEN m = path_to_zm(a); GEN n = path_to_zm(b); GEN d = subii(mulss(cc(m),dd(n)), mulss(dd(m),cc(n))); ulong res = umodiu(d, N); avma = av; return res == 0; } /* Input: a,b = 2 paths that are \Gamma_0(N)-equivalent, N = integer * Output: M in \Gamma_0(N) such that Mb=a */ static GEN gamma_equiv_matrix(GEN a, GEN b) { GEN m = path_to_ZM(a); GEN n = path_to_ZM(b); return ZM_mul(m, SL2_inv(n)); } /*************/ /* P^1(Z/NZ) */ /*************/ /* a != 0 in Z/NZ. Return v in (Z/NZ)^* such that av = gcd(a, N) (mod N)*/ static ulong Fl_inverse(ulong a, ulong N) { ulong g; return Fl_invgen(a,N,&g); } /* Input: N = integer * Output: creates P^1(Z/NZ) = [symbols, H, N] * symbols: list of vectors [x,y] that give a set of representatives * of P^1(Z/NZ) * H: an M by M grid whose value at the r,c-th place is the index of the * "standard representative" equivalent to [r,c] occuring in the first * list. If gcd(r,c,N) > 1 the grid has value 0. */ static GEN create_p1mod(ulong N) { GEN fa = factoru(N), div = divisorsu_fact(fa); ulong i, nsym = count_Manin_symbols(N, gel(fa,1)); GEN symbols = generatemsymbols(N, nsym, div); GEN H = inithashmsymbols(N,symbols); GEN invsafe = cgetg(N, t_VECSMALL), inverse = cgetg(N, t_VECSMALL); for (i = 1; i < N; i++) { invsafe[i] = Fl_invsafe(i,N); inverse[i] = Fl_inverse(i,N); } return mkvecn(7, symbols, H, utoipos(N), fa, div, invsafe, inverse); } /* Let (c : d) in P1(Z/NZ). * If c = 0 return (0:1). If d = 0 return (1:0). * Else replace by (cu : du), where u in (Z/NZ)^* such that C := cu = gcd(c,N). * In create_p1mod(), (c : d) is represented by (C:D) where D = du (mod N/c) * is smallest such that gcd(C,D) = 1. Return (C : du mod N/c), which need * not belong to P1(Z/NZ) ! A second component du mod N/c = 0 is replaced by * N/c in this case to avoid problems with array indices */ static void p1_std_form(long *pc, long *pd, GEN p1N) { ulong N = p1N_get_N(p1N); ulong u; *pc = smodss(*pc, N); if (!*pc) { *pd = 1; return; } *pd = smodss(*pd, N); if (!*pd) { *pc = 1; return; } u = p1N_get_invsafe(p1N)[*pd]; if (u) { *pc = Fl_mul(*pc,u,N); *pd = 1; return; } /* (d,N) = 1 */ u = p1N_get_inverse(p1N)[*pc]; if (u > 1) { *pc = Fl_mul(*pc,u,N); *pd = Fl_mul(*pd,u,N); } /* c | N */ if (*pc != 1) *pd %= (N / *pc); if (!*pd) *pd = N / *pc; } /* Input: v = [x,y] = elt of P^1(Z/NZ) = class in Gamma_0(N) \ PSL2(Z) * Output: returns the index of the standard rep equivalent to v */ static long p1_index(long x, long y, GEN p1N) { ulong N = p1N_get_N(p1N); GEN H = p1N_get_hash(p1N); p1_std_form(&x, &y, p1N); if (y == 1) return x+1; if (y == 0) return N+1; if (mael(H,x,y) == 0) pari_err_BUG("p1_index"); return mael(H,x,y); } /* Cusps for \Gamma_0(N) */ /* \sum_{d | N} \phi(gcd(d, N/d)), using multiplicativity. fa = factor(N) */ ulong mfnumcuspsu_fact(GEN fa) { GEN P = gel(fa,1), E = gel(fa,2); long i, l = lg(P); ulong T = 1; for (i = 1; i < l; i++) { long e = E[i], e2 = e >> 1; /* floor(E[i] / 2) */ ulong p = P[i]; if (odd(e)) T *= 2 * upowuu(p, e2); else T *= (p+1) * upowuu(p, e2-1); } return T; } ulong mfnumcuspsu(ulong n) { pari_sp av = avma; ulong t = mfnumcuspsu_fact( factoru(n) ); avma = av; return t; } /* \sum_{d | N} \phi(gcd(d, N/d)), using multiplicativity. fa = factor(N) */ GEN mfnumcusps_fact(GEN fa) { GEN P = gel(fa,1), E = gel(fa,2), T = gen_1; long i, l = lg(P); for (i = 1; i < l; i++) { GEN p = gel(P,i), c; long e = itos(gel(E,i)), e2 = e >> 1; /* floor(E[i] / 2) */ if (odd(e)) c = shifti(powiu(p, e2), 1); else c = mulii(addiu(p,1), powiu(p, e2-1)); T = T? mulii(T, c): c; } return T? T: gen_1; } GEN mfnumcusps(GEN n) { pari_sp av = avma; GEN F = check_arith_pos(n,"mfnumcusps"); if (!F) { if (lgefint(n) == 3) return utoi( mfnumcuspsu(n[2]) ); F = absZ_factor(n); } return gerepileuptoint(av, mfnumcusps_fact(F)); } /* to each cusp in \Gamma_0(N) P1(Q), represented by p/q, we associate a * unique index. Canonical representative: (1:0) or (p:q) with q | N, q < N, * p defined modulo d := gcd(N/q,q), (p,d) = 1. * Return [[N, nbcusps], H, cusps]*/ static GEN inithashcusps(GEN p1N) { ulong N = p1N_get_N(p1N); GEN div = p1N_get_div(p1N), H = zerovec(N+1); long k, ind, l = lg(div), ncusp = mfnumcuspsu_fact(p1N_get_fa(p1N)); GEN cusps = cgetg(ncusp+1, t_VEC); gel(H,1) = mkvecsmall2(0/*empty*/, 1/* first cusp: (1:0) */); gel(cusps, 1) = mkvecsmall2(1,0); ind = 2; for (k=1; k < l-1; k++) /* l-1: remove q = N */ { ulong p, q = div[k], d = ugcd(q, N/q); GEN h = const_vecsmall(d+1,0); gel(H,q+1) = h ; for (p = 0; p < d; p++) if (ugcd(p,d) == 1) { h[p+1] = ind; gel(cusps, ind) = mkvecsmall2(p,q); ind++; } } return mkvec3(mkvecsmall2(N,ind-1), H, cusps); } /* c = [p,q], (p,q) = 1, return a canonical representative for * \Gamma_0(N)(p/q) */ static GEN cusp_std_form(GEN c, GEN S) { long p, N = gel(S,1)[1], q = smodss(c[2], N); ulong u, d; if (q == 0) return mkvecsmall2(1, 0); p = smodss(c[1], N); u = Fl_inverse(q, N); q = Fl_mul(q,u, N); d = ugcd(q, N/q); return mkvecsmall2(Fl_div(p % d,u % d, d), q); } /* c = [p,q], (p,q) = 1, return the index of the corresponding cusp. * S from inithashcusps */ static ulong cusp_index(GEN c, GEN S) { long p, q; GEN H = gel(S,2); c = cusp_std_form(c, S); p = c[1]; q = c[2]; if (!mael(H,q+1,p+1)) pari_err_BUG("cusp_index"); return mael(H,q+1,p+1); } /* M a square invertible ZM, return a ZM iM such that iM M = M iM = d.Id */ static GEN ZM_inv_denom(GEN M) { GEN diM, iM = ZM_inv(M, &diM); return mkvec2(iM, diM); } /* return M^(-1) v, dinv = ZM_inv_denom(M) OR Qevproj_init(M) */ static GEN ZC_apply_dinv(GEN dinv, GEN v) { GEN x, c, iM; if (lg(dinv) == 3) { iM = gel(dinv,1); c = gel(dinv,2); } else { /* Qevproj_init */ iM = gel(dinv,2); c = gel(dinv,3); v = typ(v) == t_MAT? rowpermute(v, gel(dinv,4)) : vecpermute(v, gel(dinv,4)); } x = RgM_RgC_mul(iM, v); if (!isint1(c)) x = RgC_Rg_div(x, c); return x; } /* M an n x d ZM of rank d (basis of a Q-subspace), n >= d. * Initialize a projector on M */ GEN Qevproj_init(GEN M) { GEN v, perm, MM, iM, diM; v = ZM_indexrank(M); perm = gel(v,1); MM = rowpermute(M, perm); /* square invertible */ iM = ZM_inv(MM, &diM); return mkvec4(M, iM, diM, perm); } /* same with typechecks */ static GEN Qevproj_init0(GEN M) { switch(typ(M)) { case t_VEC: if (lg(M) == 5) return M; break; case t_COL: M = mkmat(M);/*fall through*/ case t_MAT: M = Q_primpart(M); RgM_check_ZM(M,"Qevproj_init"); return Qevproj_init(M); } pari_err_TYPE("Qevproj_init",M); return NULL; } /* T an n x n QM, pro = Qevproj_init(M), pro2 = Qevproj_init(M2); TM \subset M2. * Express these column vectors on M2's basis */ static GEN Qevproj_apply2(GEN T, GEN pro, GEN pro2) { GEN M = gel(pro,1), iM = gel(pro2,2), ciM = gel(pro2,3), perm = gel(pro2,4); return RgM_Rg_div(RgM_mul(iM, RgM_mul(rowpermute(T,perm), M)), ciM); } /* T an n x n QM, stabilizing d-dimensional Q-vector space spanned by the * d columns of M, pro = Qevproj_init(M). Return dxd matrix of T acting on M */ GEN Qevproj_apply(GEN T, GEN pro) { return Qevproj_apply2(T, pro, pro); } /* Qevproj_apply(T,pro)[,k] */ GEN Qevproj_apply_vecei(GEN T, GEN pro, long k) { GEN M = gel(pro,1), iM = gel(pro,2), ciM = gel(pro,3), perm = gel(pro,4); GEN v = RgM_RgC_mul(iM, RgM_RgC_mul(rowpermute(T,perm), gel(M,k))); return RgC_Rg_div(v, ciM); } static GEN QM_ker_r(GEN M) { return ZM_ker(Q_primpart(M)); } static GEN QM_image(GEN A) { A = vec_Q_primpart(A); return vecpermute(A, ZM_indeximage(A)); } static int cmp_dim(void *E, GEN a, GEN b) { long k; (void)E; a = gel(a,1); b = gel(b,1); k = lg(a)-lg(b); return k? ((k > 0)? 1: -1): 0; } /* FIXME: could use ZX_roots for deglim = 1 */ static GEN ZX_factor_limit(GEN T, long deglim, long *pl) { GEN fa = ZX_factor(T), P, E; long i, l; P = gel(fa,1); *pl = l = lg(P); if (deglim <= 0) return fa; E = gel(fa,2); for (i = 1; i < l; i++) if (degpol(gel(P,i)) > deglim) break; setlg(P,i); setlg(E,i); return fa; } /* Decompose the subspace H (Qevproj format) in simple subspaces. * Eg for H = msnew */ static GEN mssplit_i(GEN W, GEN H, long deglim) { ulong p, N = ms_get_N(W); long first, dim; forprime_t S; GEN T1 = NULL, T2 = NULL, V; dim = lg(gel(H,1))-1; V = vectrunc_init(dim+1); if (!dim) return V; (void)u_forprime_init(&S, 2, ULONG_MAX); vectrunc_append(V, H); first = 1; /* V[1..first-1] contains simple subspaces */ while ((p = u_forprime_next(&S))) { GEN T; long j, lV; if (N % p == 0) continue; if (T1 && T2) { T = RgM_add(T1,T2); T2 = NULL; } else { T2 = T1; T1 = T = mshecke(W, p, NULL); } lV = lg(V); for (j = first; j < lV; j++) { pari_sp av = avma; long lP; GEN Vj = gel(V,j), P = gel(Vj,1); GEN TVj = Qevproj_apply(T, Vj); /* c T | V_j */ GEN ch = QM_charpoly_ZX(TVj), fa = ZX_factor_limit(ch,deglim, &lP); GEN F = gel(fa, 1), E = gel(fa, 2); long k, lF = lg(F); if (lF == 2 && lP == 2) { if (isint1(gel(E,1))) { /* simple subspace */ swap(gel(V,first), gel(V,j)); first++; } else avma = av; } else if (lF == 1) /* discard V[j] */ { swap(gel(V,j), gel(V,lg(V)-1)); setlg(V, lg(V)-1); } else { /* can split Vj */ GEN pows; long D = 1; for (k = 1; k < lF; k++) { long d = degpol(gel(F,k)); if (d > D) D = d; } /* remove V[j] */ swap(gel(V,j), gel(V,lg(V)-1)); setlg(V, lg(V)-1); pows = RgM_powers(TVj, minss((long)2*sqrt((double)D), D)); for (k = 1; k < lF; k++) { GEN f = gel(F,k); GEN K = QM_ker_r( RgX_RgMV_eval(f, pows)) ; /* Ker f(TVj) */ GEN p = vec_Q_primpart( RgM_mul(P, K) ); vectrunc_append(V, Qevproj_init(p)); if (lg(K) == 2 || isint1(gel(E,k))) { /* simple subspace */ swap(gel(V,first), gel(V, lg(V)-1)); first++; } } if (j < first) j = first; } } if (first >= lg(V)) { gen_sort_inplace(V, NULL, cmp_dim, NULL); return V; } } pari_err_BUG("subspaces not found"); return NULL; } GEN mssplit(GEN W, GEN H, long deglim) { pari_sp av = avma; checkms(W); if (!msk_get_sign(W)) pari_err_DOMAIN("mssplit","abs(sign)","!=",gen_1,gen_0); if (!H) H = msnew(W); H = Qevproj_init0(H); return gerepilecopy(av, mssplit_i(W,H,deglim)); } /* proV = Qevproj_init of a Hecke simple subspace, return [ a_n, n <= B ] */ static GEN msqexpansion_i(GEN W, GEN proV, ulong B) { ulong p, N = ms_get_N(W), sqrtB; long i, d, k = msk_get_weight(W); forprime_t S; GEN T1=NULL, T2=NULL, TV=NULL, ch=NULL, v, dTiv, Tiv, diM, iM, L; switch(B) { case 0: return cgetg(1,t_VEC); case 1: return mkvec(gen_1); } (void)u_forprime_init(&S, 2, ULONG_MAX); while ((p = u_forprime_next(&S))) { GEN T; if (N % p == 0) continue; if (T1 && T2) { T = RgM_add(T1,T2); T2 = NULL; } else { T2 = T1; T1 = T = mshecke(W, p, NULL); } TV = Qevproj_apply(T, proV); /* T | V */ ch = QM_charpoly_ZX(TV); if (ZX_is_irred(ch)) break; ch = NULL; } if (!ch) pari_err_BUG("q-Expansion not found"); /* T generates the Hecke algebra (acting on V) */ d = degpol(ch); v = vec_ei(d, 1); /* take v = e_1 */ Tiv = cgetg(d+1, t_MAT); /* Tiv[i] = T^(i-1)v */ gel(Tiv, 1) = v; for (i = 2; i <= d; i++) gel(Tiv, i) = RgM_RgC_mul(TV, gel(Tiv,i-1)); Tiv = Q_remove_denom(Tiv, &dTiv); iM = ZM_inv(Tiv, &diM); if (dTiv) diM = gdiv(diM, dTiv); L = const_vec(B,NULL); sqrtB = usqrt(B); gel(L,1) = d > 1? mkpolmod(gen_1,ch): gen_1; for (p = 2; p <= B; p++) { pari_sp av = avma; GEN T, u, Tv, ap, P; ulong m; if (gel(L,p)) continue; /* p not prime */ T = mshecke(W, p, NULL); Tv = Qevproj_apply_vecei(T, proV, 1); /* Tp.v */ /* Write Tp.v = \sum u_i T^i v */ u = RgC_Rg_div(RgM_RgC_mul(iM, Tv), diM); ap = gerepilecopy(av, RgV_to_RgX(u, 0)); if (d > 1) ap = mkpolmod(ap,ch); else ap = simplify_shallow(ap); gel(L,p) = ap; if (!(N % p)) { /* p divides the level */ ulong C = B/p; for (m=1; m<=C; m++) if (gel(L,m)) gel(L,m*p) = gmul(gel(L,m), ap); continue; } P = powuu(p,k-1); if (p <= sqrtB) { ulong pj, oldpj = 1; for (pj = p; pj <= B; oldpj=pj, pj *= p) { GEN apj = (pj==p)? ap : gsub(gmul(ap,gel(L,oldpj)), gmul(P,gel(L,oldpj/p))); gel(L,pj) = apj; for (m = B/pj; m > 1; m--) if (gel(L,m) && m%p) gel(L,m*pj) = gmul(gel(L,m), apj); } } else { gel(L,p) = ap; for (m = B/p; m > 1; m--) if (gel(L,m)) gel(L,m*p) = gmul(gel(L,m), ap); } } return L; } GEN msqexpansion(GEN W, GEN proV, ulong B) { pari_sp av = avma; checkms(W); proV = Qevproj_init0(proV); return gerepilecopy(av, msqexpansion_i(W,proV,B)); } static GEN Qevproj_apply0(GEN T, GEN pro) { GEN iM = gel(pro,2), perm = gel(pro,4); return vec_Q_primpart(ZM_mul(iM, rowpermute(T,perm))); } /* T a ZC or ZM */ GEN Qevproj_down(GEN T, GEN pro) { GEN iM = gel(pro,2), ciM = gel(pro,3), perm = gel(pro,4); if (typ(T) == t_COL) return RgC_Rg_div(ZM_ZC_mul(iM, vecpermute(T,perm)), ciM); else return RgM_Rg_div(ZM_mul(iM, rowpermute(T,perm)), ciM); } static GEN Qevproj_star(GEN W, GEN H) { long s = msk_get_sign(W); if (s) { /* project on +/- component */ GEN A = RgM_mul(msk_get_star(W), H); A = (s > 0)? gadd(A, H): gsub(A, H); /* Im(star + sign) = Ker(star - sign) */ H = QM_image(A); H = Qevproj_apply0(H, msk_get_starproj(W)); } return H; } static GEN Tp_matrices(ulong p) { GEN v = cgetg(p+2, t_VEC); ulong i; for (i = 1; i <= p; i++) gel(v,i) = mat2(1, i-1, 0, p); gel(v,i) = mat2(p, 0, 0, 1); return v; } static GEN Up_matrices(ulong p) { GEN v = cgetg(p+1, t_VEC); ulong i; for (i = 1; i <= p; i++) gel(v,i) = mat2(1, i-1, 0, p); return v; } /* M = N/p. Classes of Gamma_0(M) / Gamma_O(N) when p | M */ static GEN NP_matrices(ulong M, ulong p) { GEN v = cgetg(p+1, t_VEC); ulong i; for (i = 1; i <= p; i++) gel(v,i) = mat2(1, 0, (i-1)*M, 1); return v; } /* M = N/p. Extra class of Gamma_0(M) / Gamma_O(N) when p \nmid M */ static GEN NP_matrix_extra(ulong M, ulong p) { long w,z, d = cbezout(p, -M, &w, &z); if (d != 1) return NULL; return mat2(w,z,M,p); } static GEN WQ_matrix(long N, long Q) { long w,z, d = cbezout(Q, N/Q, &w, &z); if (d != 1) return NULL; return mat2(Q,1,-N*z,Q*w); } GEN msnew(GEN W) { pari_sp av = avma; GEN S = mscuspidal(W, 0); ulong N = ms_get_N(W); long s = msk_get_sign(W), k = msk_get_weight(W); if (N > 1 && (!uisprime(N) || (k == 12 || k > 14))) { GEN p1N = ms_get_p1N(W), P = gel(p1N_get_fa(p1N), 1); long i, nP = lg(P)-1; GEN v = cgetg(2*nP + 1, t_COL); S = gel(S,1); /* Q basis */ for (i = 1; i <= nP; i++) { pari_sp av = avma, av2; long M = N/P[i]; GEN T1,Td, Wi = mskinit(M, k, s); GEN v1 = NP_matrices(M, P[i]); GEN vd = Up_matrices(P[i]); /* p^2 \nmid N */ if (M % P[i]) { v1 = shallowconcat(v1, mkvec(NP_matrix_extra(M,P[i]))); vd = shallowconcat(vd, mkvec(WQ_matrix(N,P[i]))); } T1 = getMorphism(W, Wi, v1); Td = getMorphism(W, Wi, vd); if (s) { T1 = Qevproj_apply2(T1, msk_get_starproj(W), msk_get_starproj(Wi)); Td = Qevproj_apply2(Td, msk_get_starproj(W), msk_get_starproj(Wi)); } av2 = avma; T1 = RgM_mul(T1,S); Td = RgM_mul(Td,S); /* multiply by S = restrict to mscusp */ gerepileallsp(av, av2, 2, &T1, &Td); gel(v,2*i-1) = T1; gel(v,2*i) = Td; } S = ZM_mul(S, QM_ker_r(matconcat(v))); /* Snew */ S = Qevproj_init(vec_Q_primpart(S)); } return gerepilecopy(av, S); } /* Solve the Manin relations for a congruence subgroup \Gamma by constructing * a well-formed fundamental domain for the action of \Gamma on upper half * space. See * Pollack and Stevens, Overconvergent modular symbols and p-adic L-functions * Annales scientifiques de l'ENS 44, fascicule 1 (2011), 1-42 * http://math.bu.edu/people/rpollack/Papers/Overconvergent_modular_symbols_and_padic_Lfunctions.pdf * * FIXME: Implemented for \Gamma = \Gamma_0(N) only. */ #if 0 /* Pollack-Stevens shift their paths so as to solve equations of the form f(z+1) - f(z) = g. We don't (to avoid mistakes) so we will have to solve eqs of the form f(z-1) - f(z) = g */ /* c = a/b; as a t_VECSMALL [a,b]; return c-1 as a t_VECSMALL */ static GEN Shift_left_cusp(GEN c) { long a=c[1], b=c[2]; return mkvecsmall2(a - b, b); } /* c = a/b; as a t_VECSMALL [a,b]; return c+1 as a t_VECSMALL */ static GEN Shift_right_cusp(GEN c) { long a=c[1], b=c[2]; return mkvecsmall2(a + b, b); } /*Input: path = [r,s] (thought of as a geodesic between these points) *Output: The path shifted by one to the left, i.e. [r-1,s-1] */ static GEN Shift_left(GEN path) { GEN r = gel(path,1), s = gel(path,2); return mkvec2(Shift_left_cusp(r), Shift_left_cusp(s)); } /*Input: path = [r,s] (thought of as a geodesic between these points) *Output: The path shifted by one to the right, i.e. [r+1,s+1] */ GEN Shift_right(GEN path) { GEN r = gel(path,1), s = gel(path,2); return mkvec2(Shift_right_cusp(r), Shift_right_cusp(s)); } #endif /* linked lists */ typedef struct list_t { GEN data; struct list_t *next; } list_t; static list_t * list_new(GEN x) { list_t *L = (list_t*)stack_malloc(sizeof(list_t)); L->data = x; L->next = NULL; return L; } static void list_insert(list_t *L, GEN x) { list_t *l = list_new(x); l->next = L->next; L->next = l; } /*Input: N > 1, p1N = P^1(Z/NZ) *Output: a connected fundamental domain for the action of \Gamma_0(N) on * upper half space. When \Gamma_0(N) is torsion free, the domain has the * property that all of its vertices are cusps. When \Gamma_0(N) has * three-torsion, 2 extra triangles need to be added. * * The domain is constructed by beginning with the triangle with vertices 0,1 * and oo. Each adjacent triangle is successively tested to see if it contains * points not \Gamma_0(N) equivalent to some point in our region. If a * triangle contains new points, it is added to the region. This process is * continued until the region can no longer be extended (and still be a * fundamental domain) by added an adjacent triangle. The list of cusps * between 0 and 1 are then returned * * Precisely, the function returns a list such that the elements of the list * with odd index are the cusps in increasing order. The even elements of the * list are either an "x" or a "t". A "t" represents that there is an element * of order three such that its fixed point is in the triangle directly * adjacent to the our region with vertices given by the cusp before and after * the "t". The "x" represents that this is not the case. */ enum { type_X, type_DO /* ? */, type_T }; static GEN form_list_of_cusps(ulong N, GEN p1N) { pari_sp av = avma; long i, position, nbC = 2; GEN v, L; list_t *C, *c; /* Let t be the index of a class in PSL2(Z) / \Gamma in our fixed enumeration * v[t] != 0 iff it is the class of z tau^r for z a previous alpha_i * or beta_i. * For \Gamma = \Gamma_0(N), the enumeration is given by p1_index. * We write cl(gamma) = the class of gamma mod \Gamma */ v = const_vecsmall(p1_size(p1N), 0); i = p1_index( 0, 1, p1N); v[i] = 1; i = p1_index( 1,-1, p1N); v[i] = 2; i = p1_index(-1, 0, p1N); v[i] = 3; /* the value is unused [debugging]: what matters is whether it is != 0 */ position = 4; /* at this point, Fund = R, v contains the classes of Id, tau, tau^2 */ C = list_new(mkvecsmall3(0,1, type_X)); list_insert(C, mkvecsmall3(1,1,type_DO)); /* C is a list of triples[a,b,t], where c = a/b is a cusp, and t is the type * of the path between c and the PREVIOUS cusp in the list, coded as * type_DO = "?", type_X = "x", type_T = "t" * Initially, C = [0/1,"?",1/1]; */ /* loop through the current set of cusps C and check to see if more cusps * should be added */ for (;;) { int done = 1; for (c = C; c; c = c->next) { GEN cusp1, cusp2, gam; long pos, b1, b2, b; if (!c->next) break; cusp1 = c->data; /* = a1/b1 */ cusp2 = (c->next)->data; /* = a2/b2 */ if (cusp2[3] != type_DO) continue; /* gam (oo -> 0) = (cusp2 -> cusp1), gam in PSL2(Z) */ gam = path_to_zm(mkpath(cusp2, cusp1)); /* = [a2,a1;b2,b1] */ /* we have normalized the cusp representation so that a1 b2 - a2 b1 = 1 */ b1 = coeff(gam,2,1); b2 = coeff(gam,2,2); /* gam.1 = (a1 + a2) / (b1 + b2) */ b = b1 + b2; /* Determine whether the adjacent triangle *below* (cusp1->cusp2) * should be added */ pos = p1_index(b1,b2, p1N); /* did we see cl(gam) before ? */ if (v[pos]) cusp2[3] = type_X; /* NO */ else { /* YES */ ulong B1, B2; v[pos] = position; i = p1_index(-(b1+b2), b1, p1N); v[i] = position+1; i = p1_index(b2, -(b1+b2), p1N); v[i] = position+2; /* add cl(gam), cl(gam*TAU), cl(gam*TAU^2) to v */ position += 3; /* gam tau gam^(-1) in \Gamma ? */ B1 = smodss(b1, N); B2 = smodss(b2, N); if ((Fl_sqr(B2,N) + Fl_sqr(B1,N) + Fl_mul(B1,B2,N)) % N == 0) cusp2[3] = type_T; else { long a1 = coeff(gam, 1,1), a2 = coeff(gam, 1,2); long a = a1 + a2; /* gcd(a,b) = 1 */ list_insert(c, mkvecsmall3(a,b,type_DO)); c = c->next; nbC++; done = 0; } } } if (done) break; } L = cgetg(nbC+1, t_VEC); i = 1; for (c = C; c; c = c->next) gel(L,i++) = c->data; return gerepilecopy(av, L); } /* W an msN. M in PSL2(Z). Return index of M in P1^(Z/NZ) = Gamma0(N) \ PSL2(Z), * and M0 in Gamma_0(N) such that M = M0 * M', where M' = chosen * section( PSL2(Z) -> P1^(Z/NZ) ). */ static GEN Gamma0N_decompose(GEN W, GEN M, long *index) { GEN p1N = msN_get_p1N(W), W3 = gel(W,3), section = msN_get_section(W); GEN A; ulong N = p1N_get_N(p1N); ulong c = umodiu(gcoeff(M,2,1), N); ulong d = umodiu(gcoeff(M,2,2), N); long s, ind = p1_index(c, d, p1N); /* as an elt of P1(Z/NZ) */ *index = W3[ind]; /* as an elt of F, E2, ... */ M = ZM_zm_mul(M, sl2_inv(gel(section,ind))); /* normalize mod +/-Id */ A = gcoeff(M,1,1); s = signe(A); if (s < 0) M = ZM_neg(M); else if (!s) { GEN C = gcoeff(M,2,1); if (signe(C) < 0) M = ZM_neg(M); } return M; } /* W an msN; as above for a path. Return [[ind], M] */ static GEN path_Gamma0N_decompose(GEN W, GEN path) { GEN p1N = msN_get_p1N(W); GEN p1index_to_ind = gel(W,3); GEN section = msN_get_section(W); GEN M = path_to_zm(path); long p1index = p1_index(cc(M), dd(M), p1N); long ind = p1index_to_ind[p1index]; GEN M0 = ZM_zm_mul(mat2_to_ZM(M), sl2_inv(gel(section,p1index))); return mkvec2(mkvecsmall(ind), M0); } /*Form generators of H_1(X_0(N),{cusps},Z) * *Input: N = integer > 1, p1N = P^1(Z/NZ) *Output: [cusp_list,E,F,T2,T3,E1] where * cusps_list = list of cusps describing fundamental domain of * \Gamma_0(N). * E = list of paths in the boundary of the fundamental domains and oriented * clockwise such that they do not contain a point * fixed by an element of order 2 and they are not an edge of a * triangle containing a fixed point of an element of order 3 * F = list of paths in the interior of the domain with each * orientation appearing separately * T2 = list of paths in the boundary of domain containing a point fixed * by an element of order 2 (oriented clockwise) * T3 = list of paths in the boundard of domain which are the edges of * some triangle containing a fixed point of a matrix of order 3 (both * orientations appear) * E1 = a sublist of E such that every path in E is \Gamma_0(N)-equivalent to * either an element of E1 or the flip (reversed orientation) of an element * of E1. * (Elements of T2 are \Gamma_0(N)-equivalent to their own flip.) * * sec = a list from 1..#p1N of matrices describing a section of the map * SL_2(Z) to P^1(Z/NZ) given by [a,b;c,d]-->[c,d]. * Given our fixed enumeration of P^1(Z/NZ), the j-th element of the list * represents the image of the j-th element of P^1(Z/NZ) under the section. */ /* insert path in set T */ static void set_insert(hashtable *T, GEN path) { hash_insert(T, path, (void*)(T->nb + 1)); } static GEN hash_to_vec(hashtable *h) { GEN v = cgetg(h->nb + 1, t_VEC); ulong i; for (i = 0; i < h->len; i++) { hashentry *e = h->table[i]; while (e) { GEN key = (GEN)e->key; long index = (long)e->val; gel(v, index) = key; e = e->next; } } return v; } static long path_to_p1_index(GEN path, GEN p1N) { GEN M = path_to_zm(path); return p1_index(cc(M), dd(M), p1N); } /* Pollack-Stevens sets */ typedef struct PS_sets_t { hashtable *F, *T2, *T31, *T32, *E1, *E2; GEN E2fromE1, stdE1; } PS_sets_t; static hashtable * set_init(long max) { return hash_create(max, (ulong(*)(void*))&hash_GEN, (int(*)(void*,void*))&gidentical, 1); } /* T = E2fromE1[i] = [c, gamma] */ static ulong E2fromE1_c(GEN T) { return itou(gel(T,1)); } static GEN E2fromE1_Zgamma(GEN T) { return gel(T,2); } static GEN E2fromE1_gamma(GEN T) { return gcoeff(gel(T,2),1,1); } static void insert_E(GEN path, PS_sets_t *S, GEN p1N) { GEN rev = vecreverse(path); long std = path_to_p1_index(rev, p1N); GEN v = gel(S->stdE1, std); if (v) { /* [s, p1], where E1[s] is the path p1 = vecreverse(path) mod \Gamma */ GEN gamma, p1 = gel(v,2); long r, s = itou(gel(v,1)); set_insert(S->E2, path); r = S->E2->nb; if (gel(S->E2fromE1, r) != gen_0) pari_err_BUG("insert_E"); gamma = gamma_equiv_matrix(rev, p1); /* reverse(E2[r]) = gamma * E1[s] */ gel(S->E2fromE1, r) = mkvec2(utoipos(s), to_famat_shallow(gamma,gen_m1)); } else { set_insert(S->E1, path); std = path_to_p1_index(path, p1N); gel(S->stdE1, std) = mkvec2(utoipos(S->E1->nb), path); } } static GEN cusp_infinity(void) { return mkvecsmall2(1,0); } static void form_E_F_T(ulong N, GEN p1N, GEN *pC, PS_sets_t *S) { GEN C, cusp_list = form_list_of_cusps(N, p1N); long nbgen = lg(cusp_list)-1, nbmanin = p1_size(p1N), r, s, i; hashtable *F, *T2, *T31, *T32, *E1, *E2; *pC = C = cgetg(nbgen+1, t_VEC); for (i = 1; i <= nbgen; i++) { GEN c = gel(cusp_list,i); gel(C,i) = mkvecsmall2(c[1], c[2]); } S->F = F = set_init(nbmanin); S->E1 = E1 = set_init(nbgen); S->E2 = E2 = set_init(nbgen); S->T2 = T2 = set_init(nbgen); S->T31 = T31 = set_init(nbgen); S->T32 = T32 = set_init(nbgen); /* T31 represents the three torsion paths going from left to right */ /* T32 represents the three torsion paths going from right to left */ for (r = 1; r < nbgen; r++) { GEN c2 = gel(cusp_list,r+1); if (c2[3] == type_T) { GEN c1 = gel(cusp_list,r), path = mkpath(c1,c2), path2 = vecreverse(path); set_insert(T31, path); set_insert(T32, path2); } } /* to record relations between E2 and E1 */ S->E2fromE1 = zerovec(nbgen); S->stdE1 = const_vec(nbmanin, NULL); /* Assumption later: path [oo,0] is E1[1], path [1,oo] is E2[1] */ { GEN oo = cusp_infinity(); GEN p1 = mkpath(oo, mkvecsmall2(0,1)); /* [oo, 0] */ GEN p2 = mkpath(mkvecsmall2(1,1), oo); /* [1, oo] */ insert_E(p1, S, p1N); insert_E(p2, S, p1N); } for (r = 1; r < nbgen; r++) { GEN c1 = gel(cusp_list,r); for (s = r+1; s <= nbgen; s++) { pari_sp av = avma; GEN c2 = gel(cusp_list,s), path; GEN d = subii(mulss(c1[1],c2[2]), mulss(c1[2],c2[1])); avma = av; if (!is_pm1(d)) continue; path = mkpath(c1,c2); if (r+1 == s) { GEN w = path; ulong hash = T31->hash(w); /* T31, T32 use the same hash function */ if (!hash_search2(T31, w, hash) && !hash_search2(T32, w, hash)) { if (gamma_equiv(path, vecreverse(path), N)) set_insert(T2, path); else insert_E(path, S, p1N); } } else { set_insert(F, mkvec2(path, mkvecsmall2(r,s))); set_insert(F, mkvec2(vecreverse(path), mkvecsmall2(s,r))); } } } setlg(S->E2fromE1, E2->nb+1); } /* v = \sum n_i g_i, g_i in Sl(2,Z), return \sum n_i g_i^(-1) */ static GEN ZSl2_star(GEN v) { long i, l; GEN w, G; if (typ(v) == t_INT) return v; G = gel(v,1); w = cgetg_copy(G, &l); for (i = 1; i < l; i++) { GEN g = gel(G,i); if (typ(g) == t_MAT) g = SL2_inv(g); gel(w,i) = g; } return ZG_normalize(mkmat2(w, gel(v,2))); } /* Input: h = set of unimodular paths, p1N = P^1(Z/NZ) = Gamma_0(N)\PSL2(Z) * Output: Each path is converted to a matrix and then an element of P^1(Z/NZ) * Append the matrix to W[12], append the index that represents * these elements of P^1 (the classes mod Gamma_0(N) via our fixed * enumeration to W[2]. */ static void paths_decompose(GEN W, hashtable *h, int flag) { GEN p1N = ms_get_p1N(W), section = ms_get_section(W); GEN v = hash_to_vec(h); long i, l = lg(v); for (i = 1; i < l; i++) { GEN e = gel(v,i); GEN M = path_to_zm(flag? gel(e,1): e); long index = p1_index(cc(M), dd(M), p1N); vecsmalltrunc_append(gel(W,2), index); gel(section, index) = M; } } static void fill_W2_W12(GEN W, PS_sets_t *S) { GEN p1N = msN_get_p1N(W); long n = p1_size(p1N); gel(W, 2) = vecsmalltrunc_init(n+1); gel(W,12) = cgetg(n+1, t_VEC); /* F contains [path, [index cusp1, index cusp2]]. Others contain paths only */ paths_decompose(W, S->F, 1); paths_decompose(W, S->E2, 0); paths_decompose(W, S->T32, 0); paths_decompose(W, S->E1, 0); paths_decompose(W, S->T2, 0); paths_decompose(W, S->T31, 0); } /* x t_VECSMALL, corresponds to a map x(i) = j, where 1 <= j <= max for all i * Return y s.t. y[j] = i or 0 (not in image) */ static GEN reverse_list(GEN x, long max) { GEN y = const_vecsmall(max, 0); long r, lx = lg(x); for (r = 1; r < lx; r++) y[ x[r] ] = r; return y; } /* go from C[a] to C[b]; return the indices of paths * E.g. if a < b * (C[a]->C[a+1], C[a+1]->C[a+2], ... C[b-1]->C[b]) * (else reverse direction) * = b - a paths */ static GEN F_indices(GEN W, long a, long b) { GEN v = cgetg(labs(b-a) + 1, t_VEC); long s, k = 1; if (a < b) { GEN index_forward = gel(W,13); for (s = a; s < b; s++) gel(v,k++) = gel(index_forward,s); } else { GEN index_backward = gel(W,14); for (s = a; s > b; s--) gel(v,k++) = gel(index_backward,s); } return v; } /* go from C[a] to C[b] via oo; return the indices of paths * E.g. if a < b * (C[a]->C[a-1], ... C[2]->C[1], * C[1]->oo, oo-> C[end], * C[end]->C[end-1], ... C[b+1]->C[b]) * a-1 + 2 + end-(b+1)+1 = end - b + a + 1 paths */ static GEN F_indices_oo(GEN W, long end, long a, long b) { GEN index_oo = gel(W,15); GEN v = cgetg(end-labs(b-a)+1 + 1, t_VEC); long s, k = 1; if (a < b) { GEN index_backward = gel(W,14); for (s = a; s > 1; s--) gel(v,k++) = gel(index_backward,s); gel(v,k++) = gel(index_backward,1); /* C[1] -> oo */ gel(v,k++) = gel(index_oo,2); /* oo -> C[end] */ for (s = end; s > b; s--) gel(v,k++) = gel(index_backward,s); } else { GEN index_forward = gel(W,13); for (s = a; s < end; s++) gel(v,k++) = gel(index_forward,s); gel(v,k++) = gel(index_forward,end); /* C[end] -> oo */ gel(v,k++) = gel(index_oo,1); /* oo -> C[1] */ for (s = 1; s < b; s++) gel(v,k++) = gel(index_forward,s); } return v; } /* index of oo -> C[1], oo -> C[end] */ static GEN indices_oo(GEN W, GEN C) { long end = lg(C)-1; GEN w, v = cgetg(2+1, t_VEC), oo = cusp_infinity(); w = mkpath(oo, gel(C,1)); /* oo -> C[1]=0 */ gel(v,1) = path_Gamma0N_decompose(W, w); w = mkpath(oo, gel(C,end)); /* oo -> C[end]=1 */ gel(v,2) = path_Gamma0N_decompose(W, w); return v; } /* index of C[1]->C[2], C[2]->C[3], ... C[end-1]->C[end], C[end]->oo * Recall that C[1] = 0, C[end] = 1 */ static GEN indices_forward(GEN W, GEN C) { long s, k = 1, end = lg(C)-1; GEN v = cgetg(end+1, t_VEC); for (s = 1; s <= end; s++) { GEN w = mkpath(gel(C,s), s == end? cusp_infinity(): gel(C,s+1)); gel(v,k++) = path_Gamma0N_decompose(W, w); } return v; } /* index of C[1]->oo, C[2]->C[1], ... C[end]->C[end-1] */ static GEN indices_backward(GEN W, GEN C) { long s, k = 1, end = lg(C)-1; GEN v = cgetg(end+1, t_VEC); for (s = 1; s <= end; s++) { GEN w = mkpath(gel(C,s), s == 1? cusp_infinity(): gel(C,s-1)); gel(v,k++) = path_Gamma0N_decompose(W, w); } return v; } /*[0,-1;1,-1]*/ static GEN mkTAU() { return mkmat22(gen_0,gen_m1, gen_1,gen_m1); } /* S */ static GEN mkS() { return mkmat22(gen_0,gen_1, gen_m1,gen_0); } /* N = integer > 1. Returns data describing Delta_0 = Z[P^1(Q)]_0 seen as * a Gamma_0(N) - module. */ static GEN msinit_N(ulong N) { GEN p1N, C, vecF, vecT2, vecT31, TAU, W, W2, singlerel, annT2, annT31; GEN F_index; ulong r, s, width; long nball, nbgen, nbp1N; hashtable *F, *T2, *T31, *T32, *E1, *E2; PS_sets_t S; W = zerovec(16); gel(W,1) = p1N = create_p1mod(N); gel(W,16)= inithashcusps(p1N); TAU = mkTAU(); if (N == 1) { gel(W,5) = mkvecsmall(1); /* cheat because sets are not disjoint if N=1 */ gel(W,11) = mkvecsmall5(0, 0, 1, 1, 2); gel(W,12) = mkvec(mat2(1,0,0,1)); gel(W,8) = mkvec( mkmat22(gen_1,gen_1, mkS(),gen_1) ); gel(W,9) = mkvec( mkmat2(mkcol3(gen_1,TAU,ZM_sqr(TAU)), mkcol3(gen_1,gen_1,gen_1)) ); return W; } nbp1N = p1_size(p1N); form_E_F_T(N,p1N, &C, &S); E1 = S.E1; E2 = S.E2; T31 = S.T31; T32 = S.T32; F = S.F; T2 = S.T2; nbgen = lg(C)-1; /* Put our paths in the order: F,E2,T32,E1,T2,T31 * W2[j] associates to the j-th element of this list its index in P1. */ fill_W2_W12(W, &S); W2 = gel(W, 2); nball = lg(W2)-1; gel(W,3) = reverse_list(W2, nbp1N); gel(W,5) = vecslice(gel(W,2), F->nb + E2->nb + T32->nb + 1, nball); gel(W,4) = reverse_list(gel(W,5), nbp1N); gel(W,13) = indices_forward(W, C); gel(W,14) = indices_backward(W, C); gel(W,15) = indices_oo(W, C); gel(W,11) = mkvecsmall5(F->nb, F->nb + E2->nb, F->nb + E2->nb + T32->nb, F->nb + E2->nb + T32->nb + E1->nb, F->nb + E2->nb + T32->nb + E1->nb + T2->nb); /* relations between T32 and T31 [not stored!] * T32[i] = - T31[i] */ /* relations of F */ width = E1->nb + T2->nb + T31->nb; /* F_index[r] = [index_1, ..., index_k], where index_i is the p1_index() * of the elementary unimodular path between 2 consecutive cusps * [in E1,E2,T2,T31 or T32] */ F_index = cgetg(F->nb+1, t_VEC); vecF = hash_to_vec(F); for (r = 1; r <= F->nb; r++) { GEN w = gel(gel(vecF,r), 2); long a = w[1], b = w[2], d = labs(b - a); /* c1 = cusp_list[a], c2 = cusp_list[b], ci != oo */ gel(F_index,r) = (nbgen-d >= d-1)? F_indices(W, a,b) : F_indices_oo(W, lg(C)-1,a,b); } singlerel = cgetg(width+1, t_VEC); /* form the single boundary relation */ for (s = 1; s <= E2->nb; s++) { /* reverse(E2[s]) = gamma * E1[c] */ GEN T = gel(S.E2fromE1,s), gamma = E2fromE1_gamma(T); gel(singlerel, E2fromE1_c(T)) = mkmat22(gen_1,gen_1, gamma,gen_m1); } for (r = E1->nb + 1; r <= width; r++) gel(singlerel, r) = gen_1; /* form the 2-torsion relations */ annT2 = cgetg(T2->nb+1, t_VEC); vecT2 = hash_to_vec(T2); for (r = 1; r <= T2->nb; r++) { GEN w = gel(vecT2,r); GEN gamma = gamma_equiv_matrix(vecreverse(w), w); gel(annT2, r) = mkmat22(gen_1,gen_1, gamma,gen_1); } /* form the 3-torsion relations */ annT31 = cgetg(T31->nb+1, t_VEC); vecT31 = hash_to_vec(T31); for (r = 1; r <= T31->nb; r++) { GEN M = path_to_ZM( vecreverse(gel(vecT31,r)) ); GEN gamma = ZM_mul(ZM_mul(M, TAU), SL2_inv(M)); gel(annT31, r) = mkmat2(mkcol3(gen_1,gamma,ZM_sqr(gamma)), mkcol3(gen_1,gen_1,gen_1)); } gel(W,6) = F_index; gel(W,7) = S.E2fromE1; gel(W,8) = annT2; gel(W,9) = annT31; gel(W,10)= singlerel; return W; } static GEN cusp_to_P1Q(GEN c) { return c[2]? gdivgs(stoi(c[1]), c[2]): mkoo(); } static GEN mspathgens_i(GEN W) { GEN R, r, g, section, gen, annT2, annT31; long i, l; checkms(W); W = get_msN(W); section = msN_get_section(W); gen = ms_get_genindex(W); l = lg(gen); g = cgetg(l,t_VEC); for (i = 1; i < l; i++) { GEN p = gel(section,gen[i]); gel(g,i) = mkvec2(cusp_to_P1Q(gel(p,1)), cusp_to_P1Q(gel(p,2))); } annT2 = msN_get_annT2(W); annT31= msN_get_annT31(W); if (ms_get_N(W) == 1) { R = cgetg(3, t_VEC); gel(R,1) = mkvec( mkvec2(gel(annT2,1), gen_1) ); gel(R,2) = mkvec( mkvec2(gel(annT31,1), gen_1) ); } else { GEN singlerel = msN_get_singlerel(W); long j, nbT2 = lg(annT2)-1, nbT31 = lg(annT31)-1, nbE1 = ms_get_nbE1(W); R = cgetg(nbT2+nbT31+2, t_VEC); l = lg(singlerel); r = cgetg(l, t_VEC); for (i = 1; i <= nbE1; i++) gel(r,i) = mkvec2(gel(singlerel, i), utoi(i)); for (; i < l; i++) gel(r,i) = mkvec2(gen_1, utoi(i)); gel(R,1) = r; j = 2; for (i = 1; i <= nbT2; i++,j++) gel(R,j) = mkvec( mkvec2(gel(annT2,i), utoi(i + nbE1)) ); for (i = 1; i <= nbT31; i++,j++) gel(R,j) = mkvec( mkvec2(gel(annT31,i), utoi(i + nbE1 + nbT2)) ); } return mkvec2(g,R); } GEN mspathgens(GEN W) { pari_sp av = avma; return gerepilecopy(av, mspathgens_i(W)); } /* Modular symbols in weight k: Hom_Gamma(Delta, Q[x,y]_{k-2}) */ /* A symbol phi is represented by the {phi(g_i)}, {phi(g'_i)}, {phi(g''_i)} * where the {g_i, g'_i, g''_i} are the Z[\Gamma]-generators of Delta, * g_i corresponds to E1, g'_i to T2, g''_i to T31. */ /* FIXME: export. T^1, ..., T^n */ static GEN RgX_powers(GEN T, long n) { GEN v = cgetg(n+1, t_VEC); long i; gel(v, 1) = T; for (i = 1; i < n; i++) gel(v,i+1) = RgX_mul(gel(v,i), T); return v; } /* g = [a,b;c,d] a mat2. Return (X^{k-2} | g)(X,Y)[X = 1]. */ static GEN voo_act_Gl2Q(GEN g, long k) { GEN mc = stoi(-coeff(g,2,1)), d = stoi(coeff(g,2,2)); return RgX_to_RgC(gpowgs(deg1pol_shallow(mc, d, 0), k-2), k-1); } struct m_act { long dim, k, p; GEN q; GEN(*act)(struct m_act *,GEN); }; /* g = [a,b;c,d]. Return (P | g)(X,Y)[X = 1] = P(dX - cY, -b X + aY)[X = 1], * for P = X^{k-2}, X^{k-3}Y, ..., Y^{k-2} */ GEN RgX_act_Gl2Q(GEN g, long k) { GEN a,b,c,d, V1,V2,V; long i; if (k == 2) return matid(1); a = gcoeff(g,1,1); b = gcoeff(g,1,2); c = gcoeff(g,2,1); d = gcoeff(g,2,2); V1 = RgX_powers(deg1pol_shallow(gneg(c), d, 0), k-2); /* d - c Y */ V2 = RgX_powers(deg1pol_shallow(a, gneg(b), 0), k-2); /*-b + a Y */ V = cgetg(k, t_MAT); gel(V,1) = RgX_to_RgC(gel(V1, k-2), k-1); for (i = 1; i < k-2; i++) { GEN v1 = gel(V1, k-2-i); /* (d-cY)^(k-2-i) */ GEN v2 = gel(V2, i); /* (-b+aY)^i */ gel(V,i+1) = RgX_to_RgC(RgX_mul(v1,v2), k-1); } gel(V,k-1) = RgX_to_RgC(gel(V2, k-2), k-1); return V; /* V[i+1] = X^i | g */ } /* z in Z[Gl2(Q)], return the matrix of z acting on V */ static GEN act_ZGl2Q(GEN z, struct m_act *T, hashtable *H) { GEN S = NULL, G, E; pari_sp av; long l, j; /* paranoia: should not occur */ if (typ(z) == t_INT) return scalarmat_shallow(z, T->dim); G = gel(z,1); l = lg(G); E = gel(z,2); av = avma; for (j = 1; j < l; j++) { GEN M, g = gel(G,j), n = gel(E,j); if (typ(g) == t_INT) /* = 1 */ M = n; /* n*Id_dim */ else { /*search in H succeeds because of preload*/ M = H? (GEN)hash_search(H,g)->val: T->act(T,g); if (is_pm1(n)) { if (signe(n) < 0) M = RgM_neg(M); } else M = RgM_Rg_mul(M, n); } if (!S) { S = M; continue; } S = gadd(S, M); if (gc_needed(av,1)) { if(DEBUGMEM>1) pari_warn(warnmem,"act_ZGl2Q, j = %ld",j); S = gerepileupto(av, S); } } return gerepilecopy(av, S); } static GEN _RgX_act_Gl2Q(struct m_act *S, GEN z) { return RgX_act_Gl2Q(z, S->k); } /* acting on (X^{k-2},...,Y^{k-2}) */ GEN RgX_act_ZGl2Q(GEN z, long k) { struct m_act T; T.k = k; T.dim = k-1; T.act=&_RgX_act_Gl2Q; return act_ZGl2Q(z, &T, NULL); } /* First pass, identify matrices in Sl_2 to convert to operators; * insert operators in hashtable. This allows GC in act_ZGl2Q */ static void hash_preload(GEN M, struct m_act *S, hashtable *H) { if (typ(M) != t_INT) { ulong h = H->hash(M); hashentry *e = hash_search2(H, M, h); if (!e) hash_insert2(H, M, S->act(S,M), h); } } /* z a sparse operator */ static void hash_vecpreload(GEN z, struct m_act *S, hashtable *H) { GEN G = gel(z,1); long i, l = lg(G); for (i = 1; i < l; i++) hash_preload(gel(G,i), S, H); } static void ZGl2QC_preload(struct m_act *S, GEN v, hashtable *H) { GEN val = gel(v,2); long i, l = lg(val); for (i = 1; i < l; i++) hash_vecpreload(gel(val,i), S, H); } /* Given a sparse vector of elements in Z[G], convert it to a (sparse) vector * of operators on V (given by t_MAT) */ static void ZGl2QC_to_act(struct m_act *S, GEN v, hashtable *H) { GEN val = gel(v,2); long i, l = lg(val); for (i = 1; i < l; i++) gel(val,i) = act_ZGl2Q(gel(val,i), S, H); } /* For all V[i] in Z[\Gamma], find the P such that P . V[i]^* = 0; * write P in basis X^{k-2}, ..., Y^{k-2} */ static GEN ZGV_tors(GEN V, long k) { long i, l = lg(V); GEN v = cgetg(l, t_VEC); for (i = 1; i < l; i++) { GEN a = ZSl2_star(gel(V,i)); gel(v,i) = ZM_ker(RgX_act_ZGl2Q(a,k)); } return v; } static long set_from_index(GEN W11, long i) { if (i <= W11[1]) return 1; if (i <= W11[2]) return 2; if (i <= W11[3]) return 3; if (i <= W11[4]) return 4; if (i <= W11[5]) return 5; return 6; } /* det M = 1 */ static void treat_index(GEN W, GEN M, long index, GEN v) { GEN W11 = gel(W,11); long shift = W11[3]; /* #F + #E2 + T32 */ switch(set_from_index(W11, index)) { case 1: /*F*/ { GEN F_index = gel(W,6), ind = gel(F_index, index); long j, l = lg(ind); for (j = 1; j < l; j++) { GEN IND = gel(ind,j), M0 = gel(IND,2); long index = mael(IND,1,1); treat_index(W, ZM_mul(M,M0), index, v); } break; } case 2: /*E2, E2[r] + gamma * E1[s] = 0 */ { long r = index - W11[1]; GEN z = gel(msN_get_E2fromE1(W), r); index = E2fromE1_c(z); M = G_ZG_mul(M, E2fromE1_Zgamma(z)); /* M * (-gamma) */ gel(v, index) = ZG_add(gel(v, index), M); break; } case 3: /*T32, T32[i] = -T31[i] */ { long T3shift = W11[5] - W11[2]; /* #T32 + #E1 + #T2 */ index += T3shift; index -= shift; gel(v, index) = ZG_add(gel(v, index), to_famat_shallow(M,gen_m1)); break; } default: /*E1,T2,T31*/ index -= shift; gel(v, index) = ZG_add(gel(v, index), to_famat_shallow(M,gen_1)); break; } } static void treat_index_trivial(GEN v, GEN W, long index) { GEN W11 = gel(W,11); long shift = W11[3]; /* #F + #E2 + T32 */ switch(set_from_index(W11, index)) { case 1: /*F*/ { GEN F_index = gel(W,6), ind = gel(F_index, index); long j, l = lg(ind); for (j = 1; j < l; j++) { GEN IND = gel(ind,j); treat_index_trivial(v, W, mael(IND,1,1)); } break; } case 2: /*E2, E2[r] + gamma * E1[s] = 0 */ { long r = index - W11[1]; long s = E2fromE1_c(gel(msN_get_E2fromE1(W), r)); v[s]--; break; } case 3: case 5: case 6: /*T32,T2,T31*/ break; case 4: /*E1*/ v[index-shift]++; break; } } static GEN M2_log(GEN W, GEN M) { GEN a = gcoeff(M,1,1), b = gcoeff(M,1,2); GEN c = gcoeff(M,2,1), d = gcoeff(M,2,2); GEN u, v, D, V; long index, s; W = get_msN(W); V = zerovec(ms_get_nbgen(W)); D = subii(mulii(a,d), mulii(b,c)); s = signe(D); if (!s) return V; if (is_pm1(D)) { /* shortcut, no need to apply Manin's trick */ if (s < 0) { b = negi(b); d = negi(d); } M = Gamma0N_decompose(W, mkmat22(a,b, c,d), &index); treat_index(W, M, index, V); } else { GEN U, B, P, Q, PQ, C1,C2; long i, l; (void)bezout(a,c,&u,&v); B = addii(mulii(b,u), mulii(d,v)); /* [u,v;-c,a] [a,b; c,d] = [1,B; 0,D], i.e. M = U [1,B;0,D] */ U = mkmat22(a,negi(v), c,u); /* {1/0 -> B/D} as \sum g_i, g_i unimodular paths */ PQ = ZV_allpnqn( gboundcf(gdiv(B,D), 0) ); P = gel(PQ,1); l = lg(P); Q = gel(PQ,2); C1 = gel(U,1); for (i = 1; i < l; i++, C1 = C2) { GEN M; C2 = ZM_ZC_mul(U, mkcol2(gel(P,i), gel(Q,i))); if (!odd(i)) C1 = ZC_neg(C1); M = Gamma0N_decompose(W, mkmat2(C1,C2), &index); treat_index(W, M, index, V); } } return V; } /* express +oo->q=a/b in terms of the Z[G]-generators, trivial action */ static void Q_log_trivial(GEN v, GEN W, GEN q) { GEN Q, W3 = gel(W,3), p1N = msN_get_p1N(W); ulong c,d, N = p1N_get_N(p1N); long i, lx; Q = Q_log_init(N, q); lx = lg(Q); c = 0; for (i = 1; i < lx; i++, c = d) { long index; d = Q[i]; if (c && !odd(i)) c = N - c; index = W3[ p1_index(c,d,p1N) ]; treat_index_trivial(v, W, index); } } static void M2_log_trivial(GEN V, GEN W, GEN M) { GEN p1N = gel(W,1), W3 = gel(W,3); ulong N = p1N_get_N(p1N); GEN a = gcoeff(M,1,1), b = gcoeff(M,1,2); GEN c = gcoeff(M,2,1), d = gcoeff(M,2,2); GEN u, v, D; long index, s; D = subii(mulii(a,d), mulii(b,c)); s = signe(D); if (!s) return; if (is_pm1(D)) { /* shortcut, not need to apply Manin's trick */ if (s < 0) d = negi(d); index = W3[ p1_index(umodiu(c,N),umodiu(d,N),p1N) ]; treat_index_trivial(V, W, index); } else { GEN U, B, P, Q, PQ; long i, l; if (!signe(c)) { Q_log_trivial(V,W,gdiv(b,d)); return; } (void)bezout(a,c,&u,&v); B = addii(mulii(b,u), mulii(d,v)); /* [u,v;-c,a] [a,b; c,d] = [1,B; 0,D], i.e. M = U [1,B;0,D] */ U = mkvec2(c, u); /* {1/0 -> B/D} as \sum g_i, g_i unimodular paths */ PQ = ZV_allpnqn( gboundcf(gdiv(B,D), 0) ); P = gel(PQ,1); l = lg(P); Q = gel(PQ,2); for (i = 1; i < l; i++, c = d) { d = addii(mulii(gel(U,1),gel(P,i)), mulii(gel(U,2),gel(Q,i))); if (!odd(i)) c = negi(c); index = W3[ p1_index(umodiu(c,N),umodiu(d,N),p1N) ]; treat_index_trivial(V, W, index); } } } static GEN cusp_to_ZC(GEN c) { switch(typ(c)) { case t_INFINITY: return mkcol2(gen_1,gen_0); case t_INT: return mkcol2(c,gen_1); case t_FRAC: return mkcol2(gel(c,1),gel(c,2)); case t_VECSMALL: return mkcol2(stoi(c[1]), stoi(c[2])); default: pari_err_TYPE("mspathlog",c); return NULL; } } static GEN path2_to_M2(GEN p) { return mkmat2(cusp_to_ZC(gel(p,1)), cusp_to_ZC(gel(p,2))); } static GEN path_to_M2(GEN p) { if (lg(p) != 3) pari_err_TYPE("mspathlog",p); switch(typ(p)) { case t_MAT: RgM_check_ZM(p,"mspathlog"); break; case t_VEC: p = path2_to_M2(p); break; default: pari_err_TYPE("mspathlog",p); } return p; } /* Expresses path p as \sum x_i g_i, where the g_i are our distinguished * generators and x_i \in Z[\Gamma]. Returns [x_1,...,x_n] */ GEN mspathlog(GEN W, GEN p) { pari_sp av = avma; checkms(W); return gerepilecopy(av, M2_log(W, path_to_M2(p))); } /** HECKE OPERATORS **/ /* [a,b;c,d] * cusp */ static GEN cusp_mul(long a, long b, long c, long d, GEN cusp) { long x = cusp[1], y = cusp[2]; long A = a*x+b*y, B = c*x+d*y, u = cgcd(A,B); if (u != 1) { A /= u; B /= u; } return mkcol2s(A, B); } /* f in Gl2(Q), act on path (zm), return path_to_M2(f.path) */ static GEN Gl2Q_act_path(GEN f, GEN path) { long a = coeff(f,1,1), b = coeff(f,1,2); long c = coeff(f,2,1), d = coeff(f,2,2); GEN c1 = cusp_mul(a,b,c,d, gel(path,1)); GEN c2 = cusp_mul(a,b,c,d, gel(path,2)); return mkmat2(c1,c2); } static GEN init_act_trivial(GEN W) { return const_vecsmall(ms_get_nbE1(W), 0); } static GEN mspathlog_trivial(GEN W, GEN p) { GEN v; W = get_msN(W); v = init_act_trivial(W); M2_log_trivial(v, W, path_to_M2(p)); return v; } /* map from W1=Hom(Delta_0(N1),Q) -> W2=Hom(Delta_0(N2),Q), weight 2, * trivial action. v a Gl2_Q or a t_VEC of Gl2_Q (\sum v[i] in Z[Gl2(Q)]). * Return the matrix attached to the action of v. */ static GEN getMorphism_trivial(GEN WW1, GEN WW2, GEN v) { GEN T, section, gen, W1 = get_msN(WW1), W2 = get_msN(WW2); long j, lv, d2; if (ms_get_N(W1) == 1) return cgetg(1,t_MAT); if (ms_get_N(W2) == 1) return zeromat(0, ms_get_nbE1(W1)); section = msN_get_section(W2); gen = msN_get_genindex(W2); d2 = ms_get_nbE1(W2); T = cgetg(d2+1, t_MAT); lv = lg(v); for (j = 1; j <= d2; j++) { GEN w = gel(section, gen[j]); GEN t = init_act_trivial(W1); pari_sp av = avma; long l; for (l = 1; l < lv; l++) M2_log_trivial(t, W1, Gl2Q_act_path(gel(v,l), w)); gel(T,j) = t; avma = av; } return shallowtrans(zm_to_ZM(T)); } static GEN RgV_sparse(GEN v, GEN *pind) { long i, l, k; GEN w = cgetg_copy(v,&l), ind = cgetg(l, t_VECSMALL); for (i = k = 1; i < l; i++) { GEN c = gel(v,i); if (typ(c) == t_INT) continue; gel(w,k) = c; ind[k] = i; k++; } setlg(w,k); setlg(ind,k); *pind = ind; return w; } static int mat2_isidentity(GEN M) { GEN A = gel(M,1), B = gel(M,2); return A[1] == 1 && A[2] == 0 && B[1] == 0 && B[2] == 1; } /* path a mat22/mat22s, return log(f.path)^* . f in sparse form */ static GEN M2_logf(GEN Wp, GEN path, GEN f) { pari_sp av = avma; GEN ind, L; long i, l; if (f) path = Gl2Q_act_path(f, path); else if (typ(gel(path,1)) == t_VECSMALL) path = path2_to_M2(path); L = M2_log(Wp, path); L = RgV_sparse(L,&ind); l = lg(L); for (i = 1; i < l; i++) gel(L,i) = ZSl2_star(gel(L,i)); if (f) ZGC_G_mul_inplace(L, mat2_to_ZM(f)); return gerepilecopy(av, mkvec2(ind,L)); } static hashtable * Gl2act_cache(long dim) { return set_init(dim*10); } /* f zm/ZM in Gl_2(Q), acts from the left on Delta, which is generated by * (g_i) as Z[Gamma1]-module, and by (G_i) as Z[Gamma2]-module. * We have f.G_j = \sum_i \lambda_{i,j} g_i, \lambda_{i,j} in Z[Gamma1] * For phi in Hom_Gamma1(D,V), g in D, phi | f is in Hom_Gamma2(D,V) and * (phi | f)(G_j) = phi(f.G_j) | f * = phi( \sum_i \lambda_{i,j} g_i ) | f * = \sum_i phi(g_i) | (\lambda_{i,j}^* f) * = \sum_i phi(g_i) | \mu_{i,j}(f) * More generally * (\sum_k (phi |v_k))(G_j) = \sum_i phi(g_i) | \Mu_{i,j} * with \Mu_{i,j} = \sum_k \mu{i,j}(v_k) * Return the \Mu_{i,j} matrix as vector of sparse columns of operators on V */ static GEN init_dual_act(GEN v, GEN W1, GEN W2, struct m_act *S) { GEN section = ms_get_section(W2), gen = ms_get_genindex(W2); /* HACK: the actions we consider in dimension 1 are trivial and in * characteristic != 2, 3 => torsion generators are 0 * [satisfy e.g. (1+gamma).g = 0 => \phi(g) | 1+gamma = 0 => \phi(g) = 0 */ long j, lv = lg(v), dim = S->dim == 1? ms_get_nbE1(W2): lg(gen)-1; GEN T = cgetg(dim+1, t_VEC); hashtable *H = Gl2act_cache(dim); for (j = 1; j <= dim; j++) { pari_sp av = avma; GEN w = gel(section, gen[j]); /* path_to_zm( E1/T2/T3 element ) */ GEN t = NULL; long k; for (k = 1; k < lv; k++) { GEN tk, f = gel(v,k); if (typ(gel(f,1)) != t_VECSMALL) f = ZM_to_zm(f); if (mat2_isidentity(f)) f = NULL; tk = M2_logf(W1, w, f); /* mu_{.,j}(v[k]) as sparse vector */ t = t? ZGCs_add(t, tk): tk; } gel(T,j) = gerepilecopy(av, t); } for (j = 1; j <= dim; j++) { ZGl2QC_preload(S, gel(T,j), H); ZGl2QC_to_act(S, gel(T,j), H); } return T; } /* modular symbol given by phi[j] = \phi(G_j) * \sum L[i]*phi[i], L a sparse column of operators */ static GEN dense_act_col(GEN col, GEN phi) { GEN s = NULL, colind = gel(col,1), colval = gel(col,2); long i, l = lg(colind), lphi = lg(phi); for (i = 1; i < l; i++) { long a = colind[i]; GEN t; if (a >= lphi) break; /* happens if k=2: torsion generator t omitted */ t = gel(phi, a); /* phi(G_a) */ t = RgM_RgC_mul(gel(colval,i), t); s = s? RgC_add(s, t): t; } return s; } /* modular symbol given by \phi( G[ind[j]] ) = val[j] * \sum L[i]*phi[i], L a sparse column of operators */ static GEN sparse_act_col(GEN col, GEN phi) { GEN s = NULL, colind = gel(col,1), colval = gel(col,2); GEN ind = gel(phi,2), val = gel(phi,3); long a, l = lg(ind); if (lg(gel(phi,1)) == 1) return RgM_RgC_mul(gel(colval,1), gel(val,1)); for (a = 1; a < l; a++) { GEN t = gel(val, a); /* phi(G_i) */ long i = zv_search(colind, ind[a]); if (!i) continue; t = RgM_RgC_mul(gel(colval,i), t); s = s? RgC_add(s, t): t; } return s; } static int phi_sparse(GEN phi) { return typ(gel(phi,1)) == t_VECSMALL; } /* phi in Hom_Gamma1(Delta, V), return the matrix whose colums are the * \sum_i phi(g_i) | \mu_{i,j} = (phi|f)(G_j), * see init_dual_act. */ static GEN dual_act(long dimV, GEN act, GEN phi) { long l = lg(act), j; GEN v = cgetg(l, t_MAT); GEN (*ACT)(GEN,GEN) = phi_sparse(phi)? sparse_act_col: dense_act_col; for (j = 1; j < l; j++) { pari_sp av = avma; GEN s = ACT(gel(act,j), phi); gel(v,j) = s? gerepileupto(av,s): zerocol(dimV); } return v; } /* in level N > 1 */ static void msk_get_st(GEN W, long *s, long *t) { GEN st = gmael(W,3,3); *s = st[1]; *t = st[2]; } static GEN msk_get_link(GEN W) { return gmael(W,3,4); } static GEN msk_get_inv(GEN W) { return gmael(W,3,5); } /* \phi in Hom(Delta, V), \phi(G_k) = phi[k]. Write \phi as * \sum_{i,j} mu_{i,j} phi_{i,j}, mu_{i,j} in Q */ static GEN getMorphism_basis(GEN W, GEN phi) { GEN R, Q, Ls, T0, T1, Ts, link, basis, inv = msk_get_inv(W); long i, j, r, s, t, dim, lvecT; if (ms_get_N(W) == 1) return ZC_apply_dinv(inv, gel(phi,1)); lvecT = lg(phi); basis = msk_get_basis(W); dim = lg(basis)-1; R = zerocol(dim); msk_get_st(W, &s, &t); link = msk_get_link(W); for (r = 2; r < lvecT; r++) { GEN Tr, L; if (r == s) continue; Tr = gel(phi,r); /* Phi(G_r), r != 1,s */ L = gel(link, r); Q = ZC_apply_dinv(gel(inv,r), Tr); /* write Phi(G_r) as sum_{a,b} mu_{a,b} Phi_{a,b}(G_r) */ for (j = 1; j < lg(L); j++) gel(R, L[j]) = gel(Q,j); } Ls = gel(link, s); T1 = gel(phi,1); /* Phi(G_1) */ gel(R, Ls[t]) = gel(T1, 1); T0 = NULL; for (i = 2; i < lg(link); i++) { GEN L; if (i == s) continue; L = gel(link,i); for (j =1 ; j < lg(L); j++) { long n = L[j]; /* phi_{i,j} = basis[n] */ GEN mu_ij = gel(R, n); GEN phi_ij = gel(basis, n), pols = gel(phi_ij,3); GEN z = RgC_Rg_mul(gel(pols, 3), mu_ij); T0 = T0? RgC_add(T0, z): z; /* += mu_{i,j} Phi_{i,j} (G_s) */ } } Ts = gel(phi,s); /* Phi(G_s) */ if (T0) Ts = RgC_sub(Ts, T0); /* solve \sum_{j!=t} mu_{s,j} Phi_{s,j}(G_s) = Ts */ Q = ZC_apply_dinv(gel(inv,s), Ts); for (j = 1; j < t; j++) gel(R, Ls[j]) = gel(Q,j); /* avoid mu_{s,t} */ for (j = t; j < lg(Q); j++) gel(R, Ls[j+1]) = gel(Q,j); return R; } /* a = s(g_i) for some modular symbol s; b in Z[G] * return s(b.g_i) = b^* . s(g_i) */ static GEN ZGl2Q_act_s(GEN b, GEN a, long k) { if (typ(b) == t_INT) { if (!signe(b)) return gen_0; switch(typ(a)) { case t_POL: a = RgX_to_RgC(a, k-1); /*fall through*/ case t_COL: a = RgC_Rg_mul(a,b); break; default: a = scalarcol_shallow(b,k-1); } } else { b = RgX_act_ZGl2Q(ZSl2_star(b), k); switch(typ(a)) { case t_POL: a = RgX_to_RgC(a, k-1); /*fall through*/ case t_COL: a = RgM_RgC_mul(b,a); break; default: a = RgC_Rg_mul(gel(b,1),a); } } return a; } static int checksymbol(GEN W, GEN s) { GEN t, annT2, annT31, singlerel; long i, k, l, nbE1, nbT2, nbT31; k = msk_get_weight(W); W = get_msN(W); nbE1 = ms_get_nbE1(W); singlerel = gel(W,10); l = lg(singlerel); if (k == 2) { for (i = nbE1+1; i < l; i++) if (!gequal0(gel(s,i))) return 0; return 1; } annT2 = msN_get_annT2(W); nbT2 = lg(annT2)-1; annT31 = msN_get_annT31(W);nbT31 = lg(annT31)-1; t = NULL; for (i = 1; i < l; i++) { GEN a = gel(s,i); a = ZGl2Q_act_s(gel(singlerel,i), a, k); t = t? gadd(t, a): a; } if (!gequal0(t)) return 0; for (i = 1; i <= nbT2; i++) { GEN a = gel(s,i + nbE1); a = ZGl2Q_act_s(gel(annT2,i), a, k); if (!gequal0(a)) return 0; } for (i = 1; i <= nbT31; i++) { GEN a = gel(s,i + nbE1 + nbT2); a = ZGl2Q_act_s(gel(annT31,i), a, k); if (!gequal0(a)) return 0; } return 1; } GEN msissymbol(GEN W, GEN s) { long k, nbgen; checkms(W); k = msk_get_weight(W); nbgen = ms_get_nbgen(W); switch(typ(s)) { case t_VEC: /* values s(g_i) */ if (lg(s)-1 != nbgen) return gen_0; break; case t_COL: if (msk_get_sign(W)) { GEN star = gel(msk_get_starproj(W), 1); if (lg(star) == lg(s)) return gen_1; } if (k == 2) /* on the dual basis of (g_i) */ { if (lg(s)-1 != nbgen) return gen_0; } else { GEN basis = msk_get_basis(W); return (lg(s) == lg(basis))? gen_1: gen_0; } break; case t_MAT: { long i, l = lg(s); GEN v = cgetg(l, t_VEC); for (i = 1; i < l; i++) gel(v,i) = msissymbol(W,gel(s,i))? gen_1: gen_0; return v; } default: return gen_0; } return checksymbol(W,s)? gen_1: gen_0; } /* map op: W1 = Hom(Delta_0(N1),V) -> W2 = Hom(Delta_0(N2),V), given by * \sum v[i], v[i] in Gl2(Q) */ static GEN getMorphism(GEN W1, GEN W2, GEN v) { struct m_act S; GEN B1, M, act; long a, l, k = msk_get_weight(W1); if (k == 2) return getMorphism_trivial(W1,W2,v); S.k = k; S.dim = k-1; S.act = &_RgX_act_Gl2Q; act = init_dual_act(v,W1,W2,&S); B1 = msk_get_basis(W1); l = lg(B1); M = cgetg(l, t_MAT); for (a = 1; a < l; a++) { pari_sp av = avma; GEN phi = dual_act(S.dim, act, gel(B1,a)); GEN D = getMorphism_basis(W2, phi); gel(M,a) = gerepilecopy(av, D); } return M; } static GEN msendo(GEN W, GEN v) { return getMorphism(W, W, v); } static GEN endo_project(GEN W, GEN e, GEN H) { if (msk_get_sign(W)) e = Qevproj_apply(e, msk_get_starproj(W)); if (H) e = Qevproj_apply(e, Qevproj_init0(H)); return e; } static GEN mshecke_i(GEN W, ulong p) { GEN v = ms_get_N(W) % p? Tp_matrices(p): Up_matrices(p); return msendo(W,v); } GEN mshecke(GEN W, long p, GEN H) { pari_sp av = avma; GEN T; checkms(W); if (p <= 1) pari_err_PRIME("mshecke",stoi(p)); T = mshecke_i(W,p); T = endo_project(W,T,H); return gerepilecopy(av, T); } static GEN msatkinlehner_i(GEN W, long Q) { long N = ms_get_N(W); GEN v; if (Q == 1) return matid(msk_get_dim(W)); if (Q == N) return msendo(W, mkvec(mat2(0,1,-N,0))); if (N % Q) pari_err_DOMAIN("msatkinlehner","N % Q","!=",gen_0,stoi(Q)); v = WQ_matrix(N, Q); if (!v) pari_err_DOMAIN("msatkinlehner","gcd(Q,N/Q)","!=",gen_1,stoi(Q)); return msendo(W,mkvec(v)); } GEN msatkinlehner(GEN W, long Q, GEN H) { pari_sp av = avma; GEN w; long k; checkms(W); k = msk_get_weight(W); if (Q <= 0) pari_err_DOMAIN("msatkinlehner","Q","<=",gen_0,stoi(Q)); w = msatkinlehner_i(W,Q); w = endo_project(W,w,H); if (k > 2 && Q != 1) w = RgM_Rg_div(w, powuu(Q,(k-2)>>1)); return gerepilecopy(av, w); } static GEN msstar_i(GEN W) { return msendo(W, mkvec(mat2(-1,0,0,1))); } GEN msstar(GEN W, GEN H) { pari_sp av = avma; GEN s; checkms(W); s = msstar_i(W); s = endo_project(W,s,H); return gerepilecopy(av, s); } #if 0 /* is \Gamma_0(N) cusp1 = \Gamma_0(N) cusp2 ? */ static int iscuspeq(ulong N, GEN cusp1, GEN cusp2) { long p1, q1, p2, q2, s1, s2, d; p1 = cusp1[1]; p2 = cusp2[1]; q1 = cusp1[2]; q2 = cusp2[2]; d = Fl_mul(smodss(q1,N),smodss(q2,N), N); d = ugcd(d, N); s1 = q1 > 2? Fl_inv(smodss(p1,q1), q1): 1; s2 = q2 > 2? Fl_inv(smodss(p2,q2), q2): 1; return Fl_mul(s1,q2,d) == Fl_mul(s2,q1,d); } #endif /* return E_c(r) */ static GEN get_Ec_r(GEN c, long k) { long p = c[1], q = c[2], u, v; GEN gr; (void)cbezout(p, q, &u, &v); gr = mat2(p, -v, q, u); /* g . (1:0) = (p:q) */ return voo_act_Gl2Q(sl2_inv(gr), k); } /* N > 1; returns the modular symbol attached to the cusp c := p/q via the rule * E_c(path from a to b in Delta_0) := E_c(b) - E_c(a), where * E_c(r) := 0 if r != c mod Gamma * v_oo | gamma_r^(-1) * where v_oo is stable by T = [1,1;0,1] (i.e x^(k-2)) and * gamma_r . (1:0) = r, for some gamma_r in SL_2(Z) * */ static GEN msfromcusp_trivial(GEN W, GEN c) { GEN section = ms_get_section(W), gen = ms_get_genindex(W); GEN S = ms_get_hashcusps(W); long j, ic = cusp_index(c, S), l = ms_get_nbE1(W)+1; GEN phi = cgetg(l, t_COL); for (j = 1; j < l; j++) { GEN vj, g = gel(section, gen[j]); /* path_to_zm(generator) */ GEN c1 = gel(g,1), c2 = gel(g,2); long i1 = cusp_index(c1, S); long i2 = cusp_index(c2, S); if (i1 == ic) vj = (i2 == ic)? gen_0: gen_1; else vj = (i2 == ic)? gen_m1: gen_0; gel(phi, j) = vj; } return phi; } static GEN msfromcusp_i(GEN W, GEN c) { GEN section, gen, S, phi; long j, ic, l, k = msk_get_weight(W); if (k == 2) { long N = ms_get_N(W); return N == 1? cgetg(1,t_COL): msfromcusp_trivial(W, c); } k = msk_get_weight(W); section = ms_get_section(W); gen = ms_get_genindex(W); S = ms_get_hashcusps(W); ic = cusp_index(c, S); l = lg(gen); phi = cgetg(l, t_COL); for (j = 1; j < l; j++) { GEN vj = NULL, g = gel(section, gen[j]); /* path_to_zm(generator) */ GEN c1 = gel(g,1), c2 = gel(g,2); long i1 = cusp_index(c1, S); long i2 = cusp_index(c2, S); if (i1 == ic) vj = get_Ec_r(c1, k); if (i2 == ic) { GEN s = get_Ec_r(c2, k); vj = vj? gsub(vj, s): gneg(s); } if (!vj) vj = zerocol(k-1); gel(phi, j) = vj; } return getMorphism_basis(W, phi); } GEN msfromcusp(GEN W, GEN c) { pari_sp av = avma; long N; checkms(W); N = ms_get_N(W); switch(typ(c)) { case t_INFINITY: c = mkvecsmall2(1,0); break; case t_INT: c = mkvecsmall2(smodis(c,N), 1); break; case t_FRAC: c = mkvecsmall2(smodis(gel(c,1),N), smodis(gel(c,2),N)); break; default: pari_err_TYPE("msfromcusp",c); } return gerepilecopy(av, msfromcusp_i(W,c)); } static GEN mseisenstein_i(GEN W) { GEN M, S = ms_get_hashcusps(W), cusps = gel(S,3); long i, l = lg(cusps); if (msk_get_weight(W)==2) l--; M = cgetg(l, t_MAT); for (i = 1; i < l; i++) gel(M,i) = msfromcusp_i(W, gel(cusps,i)); return Qevproj_star(W, QM_image(M)); } GEN mseisenstein(GEN W) { pari_sp av = avma; checkms(W); return gerepilecopy(av, Qevproj_init(mseisenstein_i(W))); } /* upper bound for log_2 |charpoly(T_p|S)|, where S is a cuspidal subspace of * dimension d, k is the weight */ #if 0 static long TpS_char_bound(ulong p, long k, long d) { /* |eigenvalue| <= 2 p^(k-1)/2 */ return d * (2 + (log2((double)p)*(k-1))/2); } #endif static long TpE_char_bound(ulong p, long k, long d) { /* |eigenvalue| <= 2 p^(k-1) */ return d * (2 + log2((double)p)*(k-1)); } GEN mscuspidal(GEN W, long flag) { pari_sp av = avma; GEN S, E, M, T, TE, chE; long bit; forprime_t F; ulong p, N; pari_timer ti; E = mseisenstein(W); N = ms_get_N(W); (void)u_forprime_init(&F, 2, ULONG_MAX); while ((p = u_forprime_next(&F))) if (N % p) break; if (DEBUGLEVEL) timer_start(&ti); T = mshecke(W, p, NULL); if (DEBUGLEVEL) timer_printf(&ti,"Tp, p = %ld", p); TE = Qevproj_apply(T, E); /* T_p | E */ if (DEBUGLEVEL) timer_printf(&ti,"Qevproj_init(E)"); bit = TpE_char_bound(p, msk_get_weight(W), lg(TE)-1); chE = QM_charpoly_ZX_bound(TE, bit); chE = ZX_radical(chE); M = RgX_RgM_eval(chE, T); S = Qevproj_init(QM_image(M)); return gerepilecopy(av, flag? mkvec2(S,E): S); } /** INIT ELLSYM STRUCTURE **/ /* V a vector of ZM. If all of them have 0 last row, return NULL. * Otherwise return [m,i,j], where m = V[i][last,j] contains the value * of smallest absolute value */ static GEN RgMV_find_non_zero_last_row(long offset, GEN V) { long i, lasti = 0, lastj = 0, lV = lg(V); GEN m = NULL; for (i = 1; i < lV; i++) { GEN M = gel(V,i); long j, n, l = lg(M); if (l == 1) continue; n = nbrows(M); for (j = 1; j < l; j++) { GEN a = gcoeff(M, n, j); if (!gequal0(a) && (!m || abscmpii(a, m) < 0)) { m = a; lasti = i; lastj = j; if (is_pm1(m)) goto END; } } } END: if (!m) return NULL; return mkvec2(m, mkvecsmall2(lasti+offset, lastj)); } /* invert the d_oo := (\gamma_oo - 1) operator, acting on * [x^(k-2), ..., y^(k-2)] */ static GEN Delta_inv(GEN doo, long k) { GEN M = RgX_act_ZGl2Q(doo, k); M = RgM_minor(M, k-1, 1); /* 1st column and last row are 0 */ return ZM_inv_denom(M); } /* The ZX P = \sum a_i x^i y^{k-2-i} is given by the ZV [a_0, ..., a_k-2]~, * return Q and d such that P = doo Q + d y^k-2, where d in Z and Q */ static GEN doo_decompose(GEN dinv, GEN P, GEN *pd) { long l = lg(P); *pd = gel(P, l-1); P = vecslice(P, 1, l-2); return shallowconcat(gen_0, ZC_apply_dinv(dinv, P)); } static GEN get_phi_ij(long i,long j,long n, long s,long t,GEN P_st,GEN Q_st,GEN d_st, GEN P_ij, GEN lP_ij, GEN dinv) { GEN ind, pols; if (i == s && j == t) { ind = mkvecsmall(1); pols = mkvec(scalarcol_shallow(gen_1, lg(P_st)-1)); /* x^{k-2} */ } else { GEN d_ij, Q_ij = doo_decompose(dinv, lP_ij, &d_ij); GEN a = ZC_Z_mul(P_ij, d_st); GEN b = ZC_Z_mul(P_st, negi(d_ij)); GEN c = RgC_sub(RgC_Rg_mul(Q_ij, d_st), RgC_Rg_mul(Q_st, d_ij)); if (i == s) { /* j != t */ ind = mkvecsmall2(1, s); pols = mkvec2(c, ZC_add(a, b)); } else { ind = mkvecsmall3(1, i, s); pols = mkvec3(c, a, b); /* image of g_1, g_i, g_s */ } pols = Q_primpart(pols); } return mkvec3(mkvecsmall3(i,j,n), ind, pols); } static GEN mskinit_trivial(GEN WN) { long dim = ms_get_nbE1(WN); return mkvec3(WN, gen_0, mkvec2(gen_0,mkvecsmall2(2, dim))); } /* sum of #cols of the matrices contained in V */ static long RgMV_dim(GEN V) { long l = lg(V), d = 0, i; for (i = 1; i < l; i++) d += lg(gel(V,i)) - 1; return d; } static GEN mskinit_nontrivial(GEN WN, long k) { GEN annT2 = gel(WN,8), annT31 = gel(WN,9), singlerel = gel(WN,10); GEN link, basis, monomials, Inv; long nbE1 = ms_get_nbE1(WN); GEN dinv = Delta_inv(ZG_neg( ZSl2_star(gel(singlerel,1)) ), k); GEN p1 = cgetg(nbE1+1, t_VEC), remove; GEN p2 = ZGV_tors(annT2, k); GEN p3 = ZGV_tors(annT31, k); GEN gentor = shallowconcat(p2, p3); GEN P_st, lP_st, Q_st, d_st; long n, i, dim, s, t, u; gel(p1, 1) = cgetg(1,t_MAT); /* dummy */ for (i = 2; i <= nbE1; i++) /* skip 1st element = (\gamma_oo-1)g_oo */ { GEN z = gel(singlerel, i); gel(p1, i) = RgX_act_ZGl2Q(ZSl2_star(z), k); } remove = RgMV_find_non_zero_last_row(nbE1, gentor); if (!remove) remove = RgMV_find_non_zero_last_row(0, p1); if (!remove) pari_err_BUG("msinit [no y^k-2]"); remove = gel(remove,2); /* [s,t] */ s = remove[1]; t = remove[2]; /* +1 because of = x^(k-2), but -1 because of Manin relation */ dim = (k-1)*(nbE1-1) + RgMV_dim(p2) + RgMV_dim(p3); /* Let (g_1,...,g_d) be the Gamma-generators of Delta, g_1 = g_oo. * We describe modular symbols by the collection phi(g_1), ..., phi(g_d) * \in V := Q[x,y]_{k-2}, with right Gamma action. * For each i = 1, .., d, let V_i \subset V be the Q-vector space of * allowed values for phi(g_i): with basis (P^{i,j}) given by the monomials * x^(j-1) y^{k-2-(j-1)}, j = 1 .. k-1 * (g_i in E_1) or the solution of the torsion equations (1 + gamma)P = 0 * (g_i in T2) or (1 + gamma + gamma^2)P = 0 (g_i in T31). All such P * are chosen in Z[x,y] with Q_content 1. * * The Manin relation (singlerel) is of the form \sum_i \lambda_i g_i = 0, * where \lambda_i = 1 if g_i in T2 or T31, and \lambda_i = (1 - \gamma_i) * for g_i in E1. * * If phi \in Hom_Gamma(Delta, V), it is defined by phi(g_i) := P_i in V * with \sum_i P_i . \lambda_i^* = 0, where (\sum n_i g_i)^* := * \sum n_i \gamma_i^(-1). * * We single out gamma_1 / g_1 (g_oo in Pollack-Stevens paper) and * write P_{i,j} \lambda_i^* = Q_{i,j} (\gamma_1 - 1)^* + d_{i,j} y^{k-2} * where d_{i,j} is a scalar and Q_{i,j} in V; we normalize Q_{i,j} to * that the coefficient of x^{k-2} is 0. * * There exist (s,t) such that d_{s,t} != 0. * A Q-basis of the (dual) space of modular symbols is given by the * functions phi_{i,j}, 2 <= i <= d, 1 <= j <= k-1, mapping * g_1 -> d_{s,t} Q_{i,j} - d_{i,j} Q_{s,t} + [(i,j)=(s,t)] x^{k-2} * If i != s * g_i -> d_{s,t} P_{i,j} * g_s -> - d_{i,j} P_{s,t} * If i = s, j != t * g_i -> d_{s,t} P_{i,j} - d_{i,j} P_{s,t} * And everything else to 0. Again we normalize the phi_{i,j} such that * their image has content 1. */ monomials = matid(k-1); /* represent the monomials x^{k-2}, ... , y^{k-2} */ if (s <= nbE1) /* in E1 */ { P_st = gel(monomials, t); lP_st = gmael(p1, s, t); /* P_{s,t} lambda_s^* */ } else /* in T2, T31 */ { P_st = gmael(gentor, s - nbE1, t); lP_st = P_st; } Q_st = doo_decompose(dinv, lP_st, &d_st); basis = cgetg(dim+1, t_VEC); link = cgetg(nbE1 + lg(gentor), t_VEC); gel(link,1) = cgetg(1,t_VECSMALL); /* dummy */ n = 1; for (i = 2; i <= nbE1; i++) { GEN L = cgetg(k, t_VECSMALL); long j; /* link[i][j] = n gives correspondance between phi_{i,j} and basis[n] */ gel(link,i) = L; for (j = 1; j < k; j++) { GEN lP_ij = gmael(p1, i, j); /* P_{i,j} lambda_i^* */ GEN P_ij = gel(monomials,j); L[j] = n; gel(basis, n) = get_phi_ij(i,j,n, s,t, P_st, Q_st, d_st, P_ij, lP_ij, dinv); n++; } } for (u = 1; u < lg(gentor); u++,i++) { GEN V = gel(gentor,u); long j, lV = lg(V); GEN L = cgetg(lV, t_VECSMALL); gel(link,i) = L; for (j = 1; j < lV; j++) { GEN lP_ij = gel(V, j); /* P_{i,j} lambda_i^* = P_{i,j} */ GEN P_ij = lP_ij; L[j] = n; gel(basis, n) = get_phi_ij(i,j,n, s,t, P_st, Q_st, d_st, P_ij, lP_ij, dinv); n++; } } Inv = cgetg(lg(link), t_VEC); gel(Inv,1) = cgetg(1, t_MAT); /* dummy */ for (i = 2; i < lg(link); i++) { GEN M, inv, B = gel(link,i); long j, lB = lg(B); if (i == s) { B = vecsplice(B, t); lB--; } /* remove phi_st */ M = cgetg(lB, t_MAT); for (j = 1; j < lB; j++) { GEN phi_ij = gel(basis, B[j]), pols = gel(phi_ij,3); gel(M, j) = gel(pols, 2); /* phi_ij(g_i) */ } if (i <= nbE1 && i != s) /* maximal rank k-1 */ inv = ZM_inv_denom(M); else /* i = s (rank k-2) or from torsion: rank k/3 or k/2 */ inv = Qevproj_init(M); gel(Inv,i) = inv; } return mkvec3(WN, gen_0, mkvec5(basis, mkvecsmall2(k, dim), mkvecsmall2(s,t), link, Inv)); } static GEN add_star(GEN W, long sign) { GEN s = msstar_i(W); GEN K = sign? QM_ker_r(gsubgs(s, sign)): cgetg(1,t_MAT); gel(W,2) = mkvec3(stoi(sign), s, Qevproj_init(K)); return W; } /* WN = msinit_N(N) */ static GEN mskinit(ulong N, long k, long sign) { GEN W, WN = msinit_N(N); if (N == 1) { GEN basis, M = RgXV_to_RgM(mfperiodpolbasis(k, 0), k-1); GEN T = cgetg(1, t_VECSMALL), ind = mkvecsmall(1); long i, l = lg(M); basis = cgetg(l, t_VEC); for (i = 1; i < l; i++) gel(basis,i) = mkvec3(T, ind, mkvec(gel(M,i))); W = mkvec3(WN, gen_0, mkvec5(basis, mkvecsmall2(k, l-1), mkvecsmall2(0,0), gen_0, Qevproj_init(M))); } else W = k == 2? mskinit_trivial(WN) : mskinit_nontrivial(WN, k); return add_star(W, sign); } GEN msinit(GEN N, GEN K, long sign) { pari_sp av = avma; GEN W; long k; if (typ(N) != t_INT) pari_err_TYPE("msinit", N); if (typ(K) != t_INT) pari_err_TYPE("msinit", K); k = itos(K); if (k < 2) pari_err_DOMAIN("msinit","k", "<", gen_2,K); if (odd(k)) pari_err_IMPL("msinit [odd weight]"); if (signe(N) <= 0) pari_err_DOMAIN("msinit","N", "<=", gen_0,N); W = mskinit(itou(N), k, sign); return gerepilecopy(av, W); } /* W = msinit, xpm integral modular symbol of weight 2, c t_FRAC * Return image of c> */ static GEN Q_xpm(GEN W, GEN xpm, GEN c) { pari_sp av = avma; GEN v; W = get_msN(W); v = init_act_trivial(W); Q_log_trivial(v, W, c); /* oo -> (a:b), c = a/b */ return gerepileuptoint(av, ZV_zc_mul(xpm, v)); } static GEN eval_single(GEN s, long k, GEN B, long v) { long i, l; GEN A = cgetg_copy(s,&l); for (i=1; i 0? 2: 1; } /* E minimal */ static GEN ellperiod(GEN E, long s) { GEN w = ellR_omega(E,DEFAULTPREC); if (s == 1) w = gel(w,1); else if (nb_components(E) == 2) w = gneg(gel(w,2)); else w = mkcomplex(gen_0, gneg(gmul2n(imag_i(gel(w,2)), 1))); return w; } /* Let W = msinit(conductor(E), 2), xpm an integral modular symbol with the same * eigenvalues as L_E. There exist a unique C such that * C*L(E,(D/.),1)_{xpm} = L(E,(D/.),1) / w1(E_D) != 0, for all D fundamental, * sign(D) = s, and such that E_D has rank 0. Return C * ellperiod(E,s) */ static GEN ell_get_Cw(GEN LE, GEN W, GEN xpm, long s) { long f, NE = ms_get_N(W); const long bit = 64; for (f = 1;; f++) { /* look for chi with conductor f coprime to N(E) and parity s * such that L(E,chi,1) != 0 */ pari_sp av = avma; GEN vchi, vx, G; long l, i; if ((f & 3) == 2 || ugcd(NE,f) != 1) continue; vx = allxpm(W, xpm, f); if (!vx) continue; G = znstar0(utoipos(f),1); vchi = chargalois(G,NULL); l = lg(vchi); for (i = 1; i < l; i++) { pari_sp av2 = avma; GEN tau, z, S, L, chi = gel(vchi,i); long o = zncharisodd(G,chi); if ((s > 0 && o) || (s < 0 && !o) || itos(zncharconductor(G, chi)) != f) continue; S = seval(G, chi, vx); if (!S) { avma = av2; continue; } L = lfuntwist(LE, mkvec2(G, zncharconj(G,chi))); z = lfun(L, gen_1, bit); tau = znchargauss(G, chi, gen_1, bit); return gdiv(gmul(z, tau), S); /* C * w */ } avma = av; } } static GEN ell_get_scale(GEN LE, GEN W, long sign, GEN x) { if (sign) return ell_get_Cw(LE, W, gel(x,1), sign); else { GEN Cwp = ell_get_Cw(LE, W, gel(x,1), 1); GEN Cwm = ell_get_Cw(LE, W, gel(x,2),-1); return mkvec2(Cwp, Cwm); } } /* E minimal */ static GEN msfromell_scale(GEN x, GEN Cw, GEN E, long s) { GEN B = int2n(32); if (s) { GEN C = gdiv(Cw, ellperiod(E,s)); return ZC_Q_mul(gel(x,1), bestappr(C,B)); } else { GEN xp = gel(x,1), Cp = gdiv(gel(Cw,1), ellperiod(E, 1)), L; GEN xm = gel(x,2), Cm = gdiv(gel(Cw,2), ellperiod(E,-1)); xp = ZC_Q_mul(xp, bestappr(Cp,B)); xm = ZC_Q_mul(xm, bestappr(Cm,B)); if (signe(ell_get_disc(E)) > 0) L = mkmat2(xp, xm); /* E(R) has 2 connected components */ else L = mkmat2(gsub(xp,xm), gmul2n(xm,1)); return mkvec3(xp, xm, L); } } /* v != 0 */ static GEN Flc_normalize(GEN v, ulong p) { long i, l = lg(v); for (i = 1; i < l; i++) if (v[i]) { if (v[i] != 1) v = Flv_Fl_div(v, v[i], p); return v; } return NULL; } /* K \cap Ker M [F_l vector spaces]. K = NULL means full space */ static GEN msfromell_ker(GEN K, GEN M, ulong l) { GEN B, Ml = ZM_to_Flm(M, l); if (K) Ml = Flm_mul(Ml, K, l); B = Flm_ker(Ml, l); if (!K) K = B; else if (lg(B) < lg(K)) K = Flm_mul(K, B, l); return K; } /* K = \cap_p Ker(T_p - a_p), 2-dimensional. Set *xl to the 1-dimensional * Fl-basis such that star . xl = sign . xl if sign != 0 and * star * xl[1] = xl[1]; star * xl[2] = -xl[2] if sign = 0 */ static void msfromell_l(GEN *pxl, GEN K, GEN star, long sign, ulong l) { GEN s = ZM_to_Flm(star, l); GEN a = gel(K,1), Sa = Flm_Flc_mul(s,a,l); GEN b = gel(K,2); GEN t = Flv_add(a,Sa,l), xp, xm; if (zv_equal0(t)) { xm = a; xp = Flv_add(b,Flm_Flc_mul(s,b,l), l); } else { xp = t; t = Flv_sub(a, Sa, l); xm = zv_equal0(t)? Flv_sub(b, Flm_Flc_mul(s,b,l), l): t; } /* xp = 0 on Im(S - 1), xm = 0 on Im(S + 1) */ if (sign > 0) *pxl = mkmat(Flc_normalize(xp, l)); else if (sign < 0) *pxl = mkmat(Flc_normalize(xm, l)); else *pxl = mkmat2(Flc_normalize(xp, l), Flc_normalize(xm, l)); } /* return a primitive symbol */ static GEN msfromell_ratlift(GEN x, GEN q) { GEN B = sqrti(shifti(q,-1)); GEN r = FpM_ratlift(x, q, B, B, NULL); if (r) r = Q_primpart(r); return r; } static int msfromell_check(GEN x, GEN vT, GEN star, long sign) { long i, l; GEN sx; if (!x) return 0; l = lg(vT); for (i = 1; i < l; i++) { GEN T = gel(vT,i); if (!gequal0(ZM_mul(T, x))) return 0; /* fail */ } sx = ZM_mul(star,x); if (sign) return ZV_equal(gel(sx,1), sign > 0? gel(x,1): ZC_neg(gel(x,1))); else return ZV_equal(gel(sx,1),gel(x,1)) && ZV_equal(gel(sx,2),ZC_neg(gel(x,2))); } GEN msfromell(GEN E0, long sign) { pari_sp av = avma, av2; GEN T, Cw, E, NE, star, q, vT, xl, xr, W, x = NULL, K = NULL; long lE, single; ulong p, l, N; forprime_t S, Sl; if (typ(E0) != t_VEC) pari_err_TYPE("msfromell",E0); lE = lg(E0); if (lE == 1) return cgetg(1,t_VEC); single = (typ(gel(E0,1)) != t_VEC); E = single ? E0: gel(E0,1); NE = ellQ_get_N(E); /* must make it integral for ellap; we have minimal model at hand */ T = obj_check(E, Q_MINIMALMODEL); if (lg(T) != 2) E = gel(T,3); N = itou(NE); av2 = avma; W = gerepilecopy(av2, mskinit(N,2,0)); star = msk_get_star(W); init_modular_small(&Sl); /* loop for p <= count_Manin_symbols(N) / 6 would be enough */ (void)u_forprime_init(&S, 2, ULONG_MAX); vT = cgetg(1, t_VEC); l = u_forprime_next(&Sl); while( (p = u_forprime_next(&S)) ) { GEN M; if (N % p == 0) continue; av2 = avma; M = RgM_Rg_sub_shallow(mshecke_i(W, p), ellap(E, utoipos(p))); M = gerepilecopy(av2, M); vT = shallowconcat(vT, mkvec(M)); /* for certification at the end */ K = msfromell_ker(K, M, l); if (lg(K) == 3) break; } if (!p) pari_err_BUG("msfromell: ran out of primes"); /* mod one l should be enough */ msfromell_l(&xl, K, star, sign, l); x = ZM_init_CRT(xl, l); q = utoipos(l); xr = msfromell_ratlift(x, q); /* paranoia */ while (!msfromell_check(xr, vT, star, sign) && (l = u_forprime_next(&Sl)) ) { GEN K = NULL; long i, lvT = lg(vT); for (i = 1; i < lvT; i++) { K = msfromell_ker(K, gel(vT,i), l); if (lg(K) == 3) break; } if (i >= lvT) { x = NULL; continue; } msfromell_l(&xl, K, star, sign, l); ZM_incremental_CRT(&x, xl, &q, l); xr = msfromell_ratlift(x, q); } /* linear form = 0 on all Im(Tp - ap) and Im(S - sign) if sign != 0 */ Cw = ell_get_scale(lfuncreate(E), W, sign, xr); if (single) x = msfromell_scale(xr, Cw, E, sign); else { /* assume all E0[i] isogenous, given by minimal models */ GEN v = cgetg(lE, t_VEC); long i; for (i=1; ik, D = S->dim; GEN a = gcoeff(f,1,1), b = gcoeff(f,1,2); GEN c = gcoeff(f,2,1), d = gcoeff(f,2,2); GEN u, z, q = S->q, mat = cgetg(D+1, t_MAT); a = modii(a,q); c = modii(c,q); z = FpX_powu(deg1pol(c,a,0), k-2, q); /* (a+cx)^(k-2) */ /* u := (b+dx) / (a+cx) mod (q,x^D) = (b/a +d/a*x) / (1 - (-c/a)*x) */ if (!equali1(a)) { GEN ai = Fp_inv(a,q); b = Fp_mul(b,ai,q); c = Fp_mul(c,ai,q); d = Fp_mul(d,ai,q); } u = deg1pol_shallow(d, b, 0); /* multiply by 1 / (1 - (-c/a)*x) */ if (signe(c)) { GEN C = Fp_neg(c,q), v = cgetg(D+2,t_POL); v[1] = evalsigne(1)|evalvarn(0); gel(v, 2) = gen_1; gel(v, 3) = C; for (j = 4; j < D+2; j++) { GEN t = Fp_mul(gel(v,j-1), C, q); if (!signe(t)) { setlg(v,j); break; } gel(v,j) = t; } u = FpXn_mul(u, v, D, q); } for (j = 1; j <= D; j++) { gel(mat,j) = RgX_to_RgC(z, D); /* (a+cx)^(k-2) * ((b+dx)/(a+cx))^(j-1) */ if (j != D) z = FpXn_mul(z, u, D, q); } return shallowtrans(mat); } static GEN moments_act(struct m_act *S, GEN f) { pari_sp av = avma; return gerepilecopy(av, moments_act_i(S,f)); } static GEN init_moments_act(GEN W, long p, long n, GEN q, GEN v) { struct m_act S; long k = msk_get_weight(W); S.p = p; S.k = k; S.q = q; S.dim = n+k-1; S.act = &moments_act; return init_dual_act(v,W,W,&S); } static void clean_tail(GEN phi, long c, GEN q) { long a, l = lg(phi); for (a = 1; a < l; a++) { GEN P = FpV_red(gel(phi,a), q); /* phi(G_a) = vector of moments */ long j, lP = lg(P); for (j = c; j < lP; j++) gel(P,j) = gen_0; /* reset garbage to 0 */ gel(phi,a) = P; } } /* concat z to all phi[i] */ static GEN concat2(GEN phi, GEN z) { long i, l; GEN v = cgetg_copy(phi,&l); for (i = 1; i < l; i++) gel(v,i) = shallowconcat(gel(phi,i), z); return v; } static GEN red_mod_FilM(GEN phi, ulong p, long k, long flag) { long a, l; GEN den = gen_1, v = cgetg_copy(phi, &l); if (flag) { phi = Q_remove_denom(phi, &den); if (!den) { den = gen_1; flag = 0; } } for (a = 1; a < l; a++) { GEN P = gel(phi,a), q = den; long j; for (j = lg(P)-1; j >= k+1; j--) { q = muliu(q,p); gel(P,j) = modii(gel(P,j),q); } q = muliu(q,p); for ( ; j >= 1; j--) gel(P,j) = modii(gel(P,j),q); gel(v,a) = P; } if (flag) v = gdiv(v, den); return v; } /* denom(C) | p^(2(k-1) - v_p(ap)) */ static GEN oms_dim2(GEN W, GEN phi, GEN C, GEN ap) { long t, i, k = mspadic_get_weight(W); long p = mspadic_get_p(W), n = mspadic_get_n(W); GEN phi1 = gel(phi,1), phi2 = gel(phi,2); GEN v, q = mspadic_get_q(W); GEN act = mspadic_get_actUp(W); t = signe(ap)? Z_lval(ap,p) : k-1; phi1 = concat2(phi1, zerovec(n)); phi2 = concat2(phi2, zerovec(n)); for (i = 1; i <= n; i++) { phi1 = dual_act(k-1, act, phi1); phi1 = dual_act(k-1, act, phi1); clean_tail(phi1, k + i*t, q); phi2 = dual_act(k-1, act, phi2); phi2 = dual_act(k-1, act, phi2); clean_tail(phi2, k + i*t, q); } C = gpowgs(C,n); v = RgM_RgC_mul(C, mkcol2(phi1,phi2)); phi1 = red_mod_FilM(gel(v,1), p, k, 1); phi2 = red_mod_FilM(gel(v,2), p, k, 1); return mkvec2(phi1,phi2); } /* flag = 0 iff alpha is a p-unit */ static GEN oms_dim1(GEN W, GEN phi, GEN alpha, long flag) { long i, k = mspadic_get_weight(W); long p = mspadic_get_p(W), n = mspadic_get_n(W); GEN q = mspadic_get_q(W); GEN act = mspadic_get_actUp(W); phi = concat2(phi, zerovec(n)); for (i = 1; i <= n; i++) { phi = dual_act(k-1, act, phi); clean_tail(phi, k + i, q); } phi = gmul(lift_shallow(gpowgs(alpha,n)), phi); phi = red_mod_FilM(phi, p, k, flag); return mkvec(phi); } /* lift polynomial P in RgX[X,Y]_{k-2} to a distribution \mu such that * \int (Y - X z)^(k-2) d\mu(z) = P(X,Y) * Return the t_VEC of k-1 first moments of \mu: \int z^i d\mu(z), 0<= i < k-1. * \sum_j (-1)^(k-2-j) binomial(k-2,j) Y^j \int z^(k-2-j) d\mu(z) = P(1,Y) * Input is P(1,Y), bin = vecbinomial(k-2): bin[j] = binomial(k-2,j-1) */ static GEN RgX_to_moments(GEN P, GEN bin) { long j, k = lg(bin); GEN Pd, Bd; if (typ(P) != t_POL) P = scalarpol(P,0); P = RgX_to_RgC(P, k-1); /* deg <= k-2 */ settyp(P, t_VEC); Pd = P+1; /* Pd[i] = coeff(P,i) */ Bd = bin+1;/* Bd[i] = binomial(k-2,i) */ for (j = 1; j < k-2; j++) { GEN c = gel(Pd,j); if (odd(j)) c = gneg(c); gel(Pd,j) = gdiv(c, gel(Bd,j)); } return vecreverse(P); } static GEN RgXC_to_moments(GEN v, GEN bin) { long i, l; GEN w = cgetg_copy(v,&l); for (i=1; iq): zerocol(S->dim); } return v; } GEN msomseval(GEN W, GEN phi, GEN path) { struct m_act S; pari_sp av = avma; GEN v, Wp; long n, vden; checkmspadic(W); if (typ(phi) != t_COL || lg(phi) != 4) pari_err_TYPE("msomseval",phi); vden = itos(gel(phi,2)); phi = gel(phi,1); n = mspadic_get_n(W); Wp= mspadic_get_Wp(W); S.k = mspadic_get_weight(W); S.p = mspadic_get_p(W); S.q = powuu(S.p, n+vden); S.dim = n + S.k - 1; S.act = &moments_act; path = path_to_M2(path); v = omseval_int(&S, phi, M2_logf(Wp,path,NULL), NULL); return gerepilecopy(av, v); } /* W = msinit(N,k,...); if flag < 0 or flag >= k-1, allow all symbols; * else commit to v_p(a_p) <= flag (ordinary if flag = 0)*/ GEN mspadicinit(GEN W, long p, long n, long flag) { pari_sp av = avma; long a, N, k; GEN P, C, M, bin, Wp, Tp, q, pn, actUp, teich, pas; checkms(W); N = ms_get_N(W); k = msk_get_weight(W); if (flag < 0) flag = 1; /* worst case */ else if (flag >= k) flag = k-1; bin = vecbinomial(k-2); Tp = mshecke(W, p, NULL); if (N % p == 0) { if ((N/p) % p == 0) pari_err_IMPL("mspadicinit when p^2 | N"); /* a_p != 0 */ Wp = W; M = gen_0; flag = (k-2) / 2; /* exact valuation */ /* will multiply by matrix with denominator p^(k-2)/2 in mspadicint. * Except if p = 2 (multiply by alpha^2) */ if (p == 2) n += k-2; else n += (k-2)/2; pn = powuu(p,n); /* For accuracy mod p^n, oms_dim1 require p^(k/2*n) */ q = powiu(pn, k/2); } else { /* p-stabilize */ long s = msk_get_sign(W); GEN M1, M2; Wp = mskinit(N*p, k, s); M1 = getMorphism(W, Wp, mkvec(mat2(1,0,0,1))); M2 = getMorphism(W, Wp, mkvec(mat2(p,0,0,1))); if (s) { GEN SW = msk_get_starproj(W), SWp = msk_get_starproj(Wp); M1 = Qevproj_apply2(M1, SW, SWp); M2 = Qevproj_apply2(M2, SW, SWp); } M = mkvec2(M1,M2); n += Z_lval(Q_denom(M), p); /*den. introduced by p-stabilization*/ /* in supersingular case: will multiply by matrix with denominator p^k * in mspadicint. Except if p = 2 (multiply by alpha^2) */ if (flag) { if (p == 2) n += 2*k-2; else n += k; } pn = powuu(p,n); /* For accuracy mod p^n, supersingular require p^((2k-1-v_p(a_p))*n) */ if (flag) /* k-1 also takes care of a_p = 0. Worst case v_p(a_p) = flag */ q = powiu(pn, 2*k-1 - flag); else q = pn; } actUp = init_moments_act(Wp, p, n, q, Up_matrices(p)); if (p == 2) C = gen_0; else { pas = matpascal(n); teich = teichmullerinit(p, n+1); P = gpowers(utoipos(p), n); C = cgetg(p, t_VEC); for (a = 1; a < p; a++) { /* powb[j+1] = ((a - w(a)) / p)^j mod p^n */ GEN powb = Fp_powers(diviuexact(subui(a, gel(teich,a)), p), n, pn); GEN Ca = cgetg(n+2, t_VEC); long j, r, ai = Fl_inv(a, p); /* a^(-1) */ gel(C,a) = Ca; for (j = 0; j <= n; j++) { GEN Caj = cgetg(j+2, t_VEC); GEN atij = gel(teich, Fl_powu(ai,j,p));/* w(a)^(-j) = w(a^(-j) mod p) */ gel(Ca,j+1) = Caj; for (r = 0; r <= j; r++) { GEN c = Fp_mul(gcoeff(pas,j+1,r+1), gel(powb, j-r+1), pn); c = Fp_mul(c,atij,pn); /* binomial(j,r)*b^(j-r)*w(a)^(-j) mod p^n */ gel(Caj,r+1) = mulii(c, gel(P,j+1)); /* p^j * c mod p^(n+j) */ } } } } return gerepilecopy(av, mkvecn(8, Wp,Tp, bin, actUp, q, mkvecsmall3(p,n,flag), M, C)); } #if 0 /* assume phi an ordinary OMS */ static GEN omsactgl2(GEN W, GEN phi, GEN M) { GEN q, Wp, act; long p, k, n; checkmspadic(W); Wp = mspadic_get_Wp(W); p = mspadic_get_p(W); k = mspadic_get_weight(W); n = mspadic_get_n(W); q = mspadic_get_q(W); act = init_moments_act(Wp, p, n, q, M); phi = gel(phi,1); return dual_act(k-1, act, gel(phi,1)); } #endif static GEN eigenvalue(GEN T, GEN x) { long i, l = lg(x); for (i = 1; i < l; i++) if (!isintzero(gel(x,i))) break; if (i == l) pari_err_DOMAIN("mstooms", "phi", "=", gen_0, x); return gdiv(RgMrow_RgC_mul(T,x,i), gel(x,i)); } /* p coprime to ap, return unit root of x^2 - ap*x + p^(k-1), accuracy p^n */ GEN mspadic_unit_eigenvalue(GEN ap, long k, GEN p, long n) { GEN sqrtD, D = subii(sqri(ap), shifti(powiu(p,k-1),2)); if (absequaliu(p,2)) { n++; sqrtD = Zp_sqrt(D, p, n); if (mod4(sqrtD) != mod4(ap)) sqrtD = negi(sqrtD); } else sqrtD = Zp_sqrtlift(D, ap, p, n); /* sqrtD = ap (mod p) */ return gmul2n(gadd(ap, cvtop(sqrtD,p,n)), -1); } /* W = msinit(N,k,...); phi = T_p/U_p - eigensymbol */ GEN mstooms(GEN W, GEN phi) { pari_sp av = avma; GEN Wp, bin, Tp, c, alpha, ap, phi0, M; long k, p, vden; checkmspadic(W); if (typ(phi) != t_COL) { if (!is_Qevproj(phi)) pari_err_TYPE("mstooms",phi); phi = gel(phi,1); if (lg(phi) != 2) pari_err_TYPE("mstooms [dim_Q (eigenspace) > 1]",phi); phi = gel(phi,1); } Wp = mspadic_get_Wp(W); Tp = mspadic_get_Tp(W); bin = mspadic_get_bin(W); k = msk_get_weight(Wp); p = mspadic_get_p(W); M = mspadic_get_M(W); phi = Q_remove_denom(phi, &c); ap = eigenvalue(Tp, phi); vden = c? Z_lvalrem(c, p, &c): 0; if (typ(M) == t_INT) { /* p | N */ GEN c1; alpha = ap; alpha = ginv(alpha); phi0 = mseval(Wp, phi, NULL); phi0 = RgXC_to_moments(phi0, bin); phi0 = Q_remove_denom(phi0, &c1); if (c1) { vden += Z_lvalrem(c1, p, &c1); c = mul_denom(c,c1); } if (umodiu(ap,p)) /* p \nmid a_p */ phi = oms_dim1(W, phi0, alpha, 0); else { phi = oms_dim1(W, phi0, alpha, 1); phi = Q_remove_denom(phi, &c1); if (c1) { vden += Z_lvalrem(c1, p, &c1); c = mul_denom(c,c1); } } } else { /* p-stabilize */ GEN M1, M2, phi1, phi2, c1; if (typ(M) != t_VEC || lg(M) != 3) pari_err_TYPE("mstooms",W); M1 = gel(M,1); M2 = gel(M,2); phi1 = RgM_RgC_mul(M1, phi); phi2 = RgM_RgC_mul(M2, phi); phi1 = mseval(Wp, phi1, NULL); phi2 = mseval(Wp, phi2, NULL); phi1 = RgXC_to_moments(phi1, bin); phi2 = RgXC_to_moments(phi2, bin); phi = Q_remove_denom(mkvec2(phi1,phi2), &c1); phi1 = gel(phi,1); phi2 = gel(phi,2); if (c1) { vden += Z_lvalrem(c1, p, &c1); c = mul_denom(c,c1); } /* all polynomials multiplied by c p^vden */ if (umodiu(ap, p)) { alpha = mspadic_unit_eigenvalue(ap, k, utoipos(p), mspadic_get_n(W)); alpha = ginv(alpha); phi0 = gsub(phi1, gmul(lift_shallow(alpha),phi2)); phi = oms_dim1(W, phi0, alpha, 0); } else { /* p | ap, alpha = [a_p, -1; p^(k-1), 0] */ long flag = mspadic_get_flag(W); if (!flag || (signe(ap) && Z_lval(ap,p) < flag)) pari_err_TYPE("mstooms [v_p(ap) > mspadicinit flag]", phi); alpha = mkmat22(ap,gen_m1, powuu(p, k-1),gen_0); alpha = ginv(alpha); phi = oms_dim2(W, mkvec2(phi1,phi2), gsqr(alpha), ap); phi = Q_remove_denom(phi, &c1); if (c1) { vden += Z_lvalrem(c1, p, &c1); c = mul_denom(c,c1); } } } if (vden) c = mul_denom(c, powuu(p,vden)); if (p == 2) alpha = gsqr(alpha); if (c) alpha = gdiv(alpha,c); if (typ(alpha) == t_MAT) { /* express in basis (omega,-p phi(omega)) */ gcoeff(alpha,2,1) = gdivgs(gcoeff(alpha,2,1), -p); gcoeff(alpha,2,2) = gdivgs(gcoeff(alpha,2,2), -p); /* at the end of mspadicint we shall multiply result by [1,0;0,-1/p]*alpha * vden + k is the denominator of this matrix */ } /* phi is integral-valued */ return gerepilecopy(av, mkcol3(phi, stoi(vden), alpha)); } /* HACK: the v[j] have different lengths */ static GEN FpVV_dotproduct(GEN v, GEN w, GEN p) { long j, l = lg(v); GEN T = cgetg(l, t_VEC); for (j = 1; j < l; j++) gel(T,j) = FpV_dotproduct(gel(v,j),w,p); return T; } /* \int (-4z)^j given \int z^j */ static GEN twistmoment_m4(GEN v) { long i, l; GEN w = cgetg_copy(v, &l); for (i = 1; i < l; i++) { GEN c = gel(v,i); if (i > 1) c = gmul2n(c, (i-1)<<1); gel(w,i) = odd(i)? c: gneg(c); } return w; } /* \int (4z)^j given \int z^j */ static GEN twistmoment_4(GEN v) { long i, l; GEN w = cgetg_copy(v, &l); for (i = 1; i < l; i++) { GEN c = gel(v,i); if (i > 1) c = gmul2n(c, (i-1)<<1); gel(w,i) = c; } return w; } /* W an mspadic, phi eigensymbol, p \nmid D. Return C(x) mod FilM */ GEN mspadicmoments(GEN W, GEN PHI, long D) { pari_sp av = avma; long na, ia, b, lphi, aD = labs(D), pp, p, k, n, vden; GEN Wp, Dact, vL, v, C, pn, phi; struct m_act S; hashtable *H; checkmspadic(W); Wp = mspadic_get_Wp(W); p = mspadic_get_p(W); k = mspadic_get_weight(W); n = mspadic_get_n(W); C = mspadic_get_C(W); if (typ(PHI) != t_COL || lg(PHI) != 4 || typ(gel(PHI,1)) != t_VEC) PHI = mstooms(W, PHI); vden = itos( gel(PHI,2) ); phi = gel(PHI,1); lphi = lg(phi); if (p == 2) { na = 2; pp = 4; } else { na = p-1; pp = p; } pn = powuu(p, n + vden); S.p = p; S.k = k; S.q = pn; S.dim = n+k-1; S.act = &moments_act; H = Gl2act_cache(ms_get_nbgen(Wp)); if (D == 1) Dact = NULL; else { GEN gaD = utoi(aD), Dk = Fp_pows(stoi(D), 2-k, pn); if (!sisfundamental(D)) pari_err_TYPE("mspadicmoments", stoi(D)); if (D % p == 0) pari_err_DOMAIN("mspadicmoments","p","|", stoi(D), utoi(p)); Dact = cgetg(aD, t_VEC); for (b = 1; b < aD; b++) { GEN z = NULL; long s = kross(D, b); if (s) { pari_sp av2 = avma; GEN d; z = moments_act_i(&S, mkmat22(gaD,utoipos(b), gen_0,gaD)); d = s > 0? Dk: Fp_neg(Dk, pn); z = equali1(d)? gerepilecopy(av2, z) : gerepileupto(av2, FpM_Fp_mul(z, d, pn)); } gel(Dact,b) = z; } } vL = cgetg(na+1,t_VEC); /* first pass to precompute log(paths), preload matrices and allow GC later */ for (ia = 1; ia <= na; ia++) { GEN path, La; long a = (p == 2 && ia == 2)? -1: ia; if (Dact) { /* twist by D */ La = cgetg(aD, t_VEC); for (b = 1; b < aD; b++) { GEN Actb = gel(Dact,b); if (!Actb) continue; /* oo -> a/pp + b/|D|*/ path = mkmat22(gen_1, addii(mulss(a, aD), muluu(pp, b)), gen_0, muluu(pp, aD)); gel(La,b) = M2_logf(Wp,path,NULL); ZGl2QC_preload(&S, gel(La,b), H); } } else { path = mkmat22(gen_1,stoi(a), gen_0, utoipos(pp)); La = M2_logf(Wp,path,NULL); ZGl2QC_preload(&S, La, H); } gel(vL,ia) = La; } v = cgetg(na+1,t_VEC); /* second pass, with GC */ for (ia = 1; ia <= na; ia++) { pari_sp av2 = avma; GEN vca, Ca = gel(C,ia), La = gel(vL,ia), va = cgetg(lphi, t_VEC); long i; if (!Dact) vca = omseval_int(&S, phi, La, H); else { /* twist by D */ vca = cgetg(lphi,t_VEC); for (b = 1; b < aD; b++) { GEN T, Actb = gel(Dact,b); if (!Actb) continue; T = omseval_int(&S, phi, gel(La,b), H); for (i = 1; i < lphi; i++) { GEN z = FpM_FpC_mul(Actb, gel(T,i), pn); gel(vca,i) = b==1? z: ZC_add(gel(vca,i), z); } } } if (p != 2) { for (i=1; i= n) break; /* S = log(1+x)^j */ gel(vlog,j+1) = mspadicint(oms,teichi,S); S = S? RgXn_mul(S, L, n+1): L; } m = j; u = utoipos(p == 2? 5: 1+p); logu = glog(cvtop(u, utoipos(p), 4*m), 0); X = gdiv(pol_x(0), logu); s = cgetg(m+1, t_VEC); s2 = oms_is_supersingular(oms)? cgetg(m+1, t_VEC): NULL; bin = pol_1(0); for (j = 0; j < m; j++) { /* bin = binomial(x/log(1+p+O(p^(4*n))), j) mod x^m */ GEN a, v = mspadicint_RgXlog(bin, vlog); int done = 1; gel(s,j+1) = a = gel(v,1); if (!gequal0(a) || valp(a) > 0) done = 0; else setlg(s,j+1); if (s2) { gel(s2,j+1) = a = gel(v,2); if (!gequal0(a) || valp(a) > 0) done = 0; else setlg(s2,j+1); } if (done || j == m-1) break; bin = RgXn_mul(bin, gdivgs(gsubgs(X, j), j+1), m); } s = RgV_to_ser(s,0,lg(s)+1); if (s2) { s2 = RgV_to_ser(s2,0,lg(s2)+1); s = mkvec2(s, s2); } if (kross(oms_get_D(oms), p) >= 0) return gerepilecopy(av, s); return gerepileupto(av, gneg(s)); } void mspadic_parse_chi(GEN s, GEN *s1, GEN *s2) { if (!s) *s1 = *s2 = gen_0; else switch(typ(s)) { case t_INT: *s1 = *s2 = s; break; case t_VEC: if (lg(s) == 3) { *s1 = gel(s,1); *s2 = gel(s,2); if (typ(*s1) == t_INT && typ(*s2) == t_INT) break; } default: pari_err_TYPE("mspadicL",s); *s1 = *s2 = NULL; } } /* oms from mspadicmoments * r-th derivative of L(f,chi^s,psi) in direction - s \in Z_p \times \Z/(p-1)\Z, s-> chi^s=<\chi>^s_1 omega^s_2) - Z -> Z_p \times \Z/(p-1)\Z par s-> (s, s mod p-1). */ GEN mspadicL(GEN oms, GEN s, long r) { pari_sp av = avma; GEN s1, s2, z, S; long p, n, teich; checkoms(oms); p = oms_get_p(oms); n = oms_get_n(oms); mspadic_parse_chi(s, &s1,&s2); teich = umodiu(subii(s2,s1), p==2? 2: p-1); S = xlog1x(n, itos(s1), r, &teich); z = mspadicint(oms, teich, S); if (lg(z) == 2) z = gel(z,1); if (kross(oms_get_D(oms), p) < 0) z = gneg(z); return gerepilecopy(av, z); } /****************************************************************************/ struct siegel { GEN V, Ast; long N; /* level */ long oo; /* index of the [oo,0] path */ long k1, k2; /* two distinguished indices */ long n; /* #W, W = initial segment [in siegelstepC] already normalized */ }; static void siegel_init(struct siegel *C, GEN M) { GEN CPI, CP, MM, V, W, Ast; GEN m = gel(M,11), M2 = gel(M,2), S = msN_get_section(M); GEN E2fromE1 = msN_get_E2fromE1(M); long m0 = lg(M2)-1; GEN E2 = vecslice(M2, m[1]+1, m[2]);/* E2 */ GEN E1T = vecslice(M2, m[3]+1, m0); /* E1,T2,T31 */ GEN L = shallowconcat(E1T, E2); long i, l = lg(L), n = lg(E1T)-1, lE = lg(E2); Ast = cgetg(l, t_VECSMALL); for (i = 1; i < lE; ++i) { long j = E2fromE1_c(gel(E2fromE1,i)); Ast[n+i] = j; Ast[j] = n+i; } for (; i<=n; ++i) Ast[i] = i; MM = cgetg (l,t_VEC); for (i = 1; i < l; i++) { GEN c = gel(S, L[i]); long c12, c22, c21 = ucoeff(c,2,1); if (!c21) { gel(MM,i) = gen_0; continue; } c22 = ucoeff(c,2,2); if (!c22) { gel(MM,i) = gen_m1; continue; } c12 = ucoeff(c,1,2); gel(MM,i) = gdivgs(stoi(c12), c22); /* right extremity > 0 */ } CP = indexsort(MM); CPI = cgetg(l, t_VECSMALL); V = cgetg(l, t_VEC); W = cgetg(l, t_VECSMALL); for (i = 1; i < l; ++i) { gel(V,i) = mat2_to_ZM(gel(S, L[CP[i]])); CPI[CP[i]] = i; } for (i = 1; i < l; ++i) W[CPI[i]] = CPI[Ast[i]]; C->V = V; C->Ast = W; C->n = 0; C->oo = 2; C->N = ms_get_N(M); } static double ZMV_size(GEN v) { long i, l = lg(v); GEN z = cgetg(l, t_VECSMALL); for (i = 1; i < l; i++) z[i] = gexpo(gel(v,i)); return ((double)zv_sum(z)) / (4*(l-1)); } /* apply permutation perm to struct S. Don't follow k1,k2 */ static void siegel_perm0(struct siegel *S, GEN perm) { long i, l = lg(S->V); GEN V2 = cgetg(l, t_VEC), Ast2 = cgetg(l, t_VECSMALL); GEN V = S->V, Ast = S->Ast; for (i = 1; i < l; i++) gel(V2,perm[i]) = gel(V,i); for (i = 1; i < l; i++) Ast2[perm[i]] = perm[Ast[i]]; S->oo = perm[S->oo]; S->Ast = Ast2; S->V = V2; } /* apply permutation perm to full struct S */ static void siegel_perm(struct siegel *S, GEN perm) { siegel_perm0(S, perm); S->k1 = perm[S->k1]; S->k2 = perm[S->k2]; } /* cyclic permutation of lg = l-1 moving a -> 1, a+1 -> 2, etc. */ static GEN rotate_perm(long l, long a) { GEN p = cgetg(l, t_VECSMALL); long i, j = 1; for (i = a; i < l; i++) p[i] = j++; for (i = 1; i < a; i++) p[i] = j++; return p; } /* a1 < c1 <= a2 < c2*/ static GEN basic_op_perm(long l, long a1, long a2, long c1, long c2) { GEN p = cgetg(l, t_VECSMALL); long i, j = 1; p[a1] = j++; for (i = c1; i < a2; i++) p[i] = j++; for (i = a1+1; i < c1; i++) p[i] = j++; p[a2] = j++; for (i = c2; i < l; i++) p[i] = j++; for (i = 1; i < a1; i++) p[i] = j++; for (i = a2+1; i < c2; i++) p[i] = j++; return p; } static GEN basic_op_perm_elliptic(long l, long a1) { GEN p = cgetg(l, t_VECSMALL); long i, j = 1; p[a1] = j++; for (i = 1; i < a1; i++) p[i] = j++; for (i = a1+1; i < l; i++) p[i] = j++; return p; } static GEN ZM2_rev(GEN T) { return mkmat2(gel(T,2), ZC_neg(gel(T,1))); } /* In place, V = vector of consecutive paths, between x <= y. * V[x..y-1] <- g*V[x..y-1] */ static void path_vec_mul(GEN V, long x, long y, GEN g) { long j; GEN M; if (x == y) return; M = gel(V,x); gel(V,x) = ZM_mul(g,M); for (j = x+1; j < y; j++) /* V[j] <- g*V[j], optimized */ { GEN Mnext = gel(V,j); /* Mnext[,1] = M[,2] */ GEN gM = gel(V,j-1), u = gel(gM,2); if (!ZV_equal(gel(M,2), gel(Mnext,1))) u = ZC_neg(u); gel(V,j) = mkmat2(u, ZM_ZC_mul(g,gel(Mnext,2))); M = Mnext; } } static long prev(GEN V, long i) { return (i == 1)? lg(V)-1: i-1; } static long next(GEN V, long i) { return (i == lg(V)-1)? 1: i+1; } static GEN ZM_det2(GEN u, GEN v) { GEN a = gel(u,1), c = gel(u,2); GEN b = gel(v,1), d = gel(v,2); return subii(mulii(a,d), mulii(b,c)); } static GEN ZM2_det(GEN T) { return ZM_det2(gel(T,1),gel(T,2)); } static void fill1(GEN V, long a) { long p = prev(V,a), n = next(V,a); GEN u = gmael(V,p,2), v = gmael(V,n,1); if (signe(ZM_det2(u,v)) < 0) v = ZC_neg(v); gel(V,a) = mkmat2(u, v); } /* a1 < a2 */ static void fill2(GEN V, long a1, long a2) { if (a2 != a1+1) { fill1(V,a1); fill1(V,a2); } /* non adjacent, reconnect */ else /* parabolic */ { long p = prev(V,a1), n = next(V,a2); GEN u, v, C = gmael(V,a1,2), mC = ZC_neg(C); /* = \pm V[a2][1] */ u = gmael(V,p,2); v = (signe(ZM_det2(u,C)) < 0)? mC: C; gel(V,a1) = mkmat2(u,v); v = gmael(V,n,1); u = (signe(ZM_det2(C,v)) < 0)? mC: C; gel(V,a2) = mkmat2(u,v); } } /* DU = det(U), return g = T*U^(-1) or NULL if not in Gamma0(N); if N = 0, * only test whether g is integral */ static GEN ZM2_div(GEN T, GEN U, GEN DU, long N) { GEN a=gcoeff(U,1,1), b=gcoeff(U,1,2), c=gcoeff(U,2,1), d=gcoeff(U,2,2); GEN e=gcoeff(T,1,1), f=gcoeff(T,1,2), g=gcoeff(T,2,1), h=gcoeff(T,2,2); GEN A, B, C, D, r; C = dvmdii(subii(mulii(d,g), mulii(c,h)), DU, &r); if (r != gen_0 || (N && smodis(C,N))) return NULL; A = dvmdii(subii(mulii(d,e), mulii(c,f)), DU, &r); if (r != gen_0) return NULL; B = dvmdii(subii(mulii(a,f), mulii(b,e)), DU, &r); if (r != gen_0) return NULL; D = dvmdii(subii(mulii(a,h), mulii(g,b)), DU, &r); if (r != gen_0) return NULL; return mkmat22(A,B,C,D); } static GEN get_g(struct siegel *S, long a1) { long a2 = S->Ast[a1]; GEN a = gel(S->V,a1), ar = ZM2_rev(gel(S->V,a2)), Dar = ZM2_det(ar); GEN g = ZM2_div(a, ar, Dar, S->N); if (!g) { GEN tau = mkmat22(gen_0,gen_m1, gen_1,gen_m1); /*[0,-1;1,-1]*/ g = ZM2_div(ZM_mul(ar, tau), ar, Dar, 0); } return g; } /* input V = (X1 a X2 | X3 a^* X4) + Ast * a1 = index of a * a2 = index of a^*, inferred from a1. We must have a != a^* * c1 = first cut [ index of first path in X3 ] * c2 = second cut [ either in X4 or X1, index of first path ] * Assume a < a^* (cf Paranoia below): c1 or c2 must be in * ]a,a^*], and the other in the "complement" ]a^*,a] */ static void basic_op(struct siegel *S, long a1, long c1, long c2) { long l = lg(S->V), a2 = S->Ast[a1]; GEN g; if (a1 == a2) { /* a = a^* */ g = get_g(S, a1); path_vec_mul(S->V, a1+1, l, g); siegel_perm(S, basic_op_perm_elliptic(l, a1)); /* fill the hole left at a1, reconnect the path */ fill1(S->V, a1); return; } /* Paranoia: (a,a^*) conjugate, call 'a' the first one */ if (a2 < a1) lswap(a1,a2); /* Now a1 < a2 */ if (c1 <= a1 || c1 > a2) lswap(c1,c2); /* ensure a1 < c1 <= a2 */ if (c2 < a1) { /* if cut c2 is in X1 = X11|X12, rotate to obtain (a X2 | X3 a^* X4 X11|X12): then a1 = 1 */ GEN p = rotate_perm(l, a1); siegel_perm(S, p); a1 = 1; /* = p[a1] */ a2 = S->Ast[1]; /* > a1 */ c1 = p[c1]; c2 = p[c2]; } /* Now a1 < c1 <= a2 < c2; a != a^* */ g = get_g(S, a1); if (S->oo >= c1 && S->oo < c2) /* W inside [c1..c2[ */ { /* c2 -> c1 excluding a1 */ GEN gi = SL2_inv(g); /* g a^* = a; gi a = a^* */ path_vec_mul(S->V, 1, a1, gi); path_vec_mul(S->V, a1+1, c1, gi); path_vec_mul(S->V, c2, l, gi); } else { /* c1 -> c2 excluding a2 */ path_vec_mul(S->V, c1, a2, g); path_vec_mul(S->V, a2+1, c2, g); } siegel_perm(S, basic_op_perm(l, a1,a2, c1,c2)); /* fill the holes left at a1,a2, reconnect the path */ fill2(S->V, a1, a2); } /* a = a^* (elliptic case) */ static void basic_op_elliptic(struct siegel *S, long a1) { long l = lg(S->V); GEN g = get_g(S, a1); path_vec_mul(S->V, a1+1, l, g); siegel_perm(S, basic_op_perm_elliptic(l, a1)); /* fill the hole left at a1 (now at 1), reconnect the path */ fill1(S->V, 1); } /* input V = W X a b Y a^* Z b^* T, W already normalized * X = [n+1, k1-1], Y = [k2+1, Ast[k1]-1], * Z = [Ast[k1]+1, Ast[k2]-1], T = [Ast[k2]+1, oo]. * Assume that X doesn't start by c c^* or a b a^* b^*. */ static void siegelstep(struct siegel *S) { if (S->Ast[S->k1] == S->k1) { basic_op_elliptic(S, S->k1); S->n++; } else if (S->Ast[S->k1] == S->k1+1) { basic_op(S, S->k1, S->Ast[S->k1], 1); /* 1: W starts there */ S->n += 2; } else { basic_op(S, S->k2, S->Ast[S->k1], 1); /* 1: W starts there */ basic_op(S, S->k1, S->k2, S->Ast[S->k2]); basic_op(S, S->Ast[S->k2], S->k2, S->Ast[S->k1]); basic_op(S, S->k1, S->Ast[S->k1], S->Ast[S->k2]); S->n += 4; } } /* normalize hyperbolic polygon */ static void mssiegel(struct siegel *S) { pari_sp av = avma; long k, t, nv; #ifdef COUNT long countset[16]; for (k = 0; k < 16; k++) countset[k] = 0; #endif nv = lg(S->V)-1; if (DEBUGLEVEL>1) err_printf("nv = %ld, expo = %.2f\n", nv,ZMV_size(S->V)); t = 0; while (S->n < nv) { if (S->Ast[S->n+1] == S->n+1) { S->n++; continue; } if (S->Ast[S->n+1] == S->n+2) { S->n += 2; continue; } if (S->Ast[S->n+1] == S->n+3 && S->Ast[S->n+2] == S->n+4) { S->n += 4; continue; } k = nv; while (k > S->n) { if (S->Ast[k] == k) { k--; continue; } if (S->Ast[k] == k-1) { k -= 2; continue; } if (S->Ast[k] == k-2 && S->Ast[k-1] == k-3) { k -= 4; continue; } break; } if (k != nv) { siegel_perm0(S, rotate_perm(nv+1, k+1)); S->n += nv-k; } for (k = S->n+1; k <= nv; k++) if (S->Ast[k] <= k) { t = S->Ast[k]; break; } S->k1 = t; S->k2 = t+1; #ifdef COUNT countset[ ((S->k1-1 == S->n) | ((S->k2 == S->Ast[S->k1]-1) << 1) | ((S->Ast[S->k1] == S->Ast[S->k2]-1) << 2) | ((S->Ast[S->k2] == nv) << 3)) ]++; #endif siegelstep(S); if (gc_needed(av,2)) { if(DEBUGMEM>1) pari_warn(warnmem,"mspolygon, n = %ld",S->n); gerepileall(av, 2, &S->V, &S->Ast); } } if (DEBUGLEVEL>1) err_printf("expo = %.2f\n", ZMV_size(S->V)); #ifdef COUNT for (k = 0; k < 16; k++) err_printf("%3ld: %6ld\n", k, countset[k]); #endif } /* return a vector of char */ static GEN Ast2v(GEN Ast) { long j = 0, k, l = lg(Ast); GEN v = const_vec(l-1, NULL); for (k=1; k < l; k++) { char *sj; if (gel(v,k)) continue; j++; sj = stack_sprintf("$%ld$", j); gel(v,k) = (GEN)sj; if (Ast[k] != k) gel(v,Ast[k]) = (GEN)stack_sprintf("$%ld^*$", j); } return v; }; static GEN M2Q(GEN p) { GEN c = gel(p,1); return gdiv(gel(c,1), gel(c,2)); } static void maybe_decorate(pari_str *s, GEN Ast, GEN G, long k, double c) { GEN g; if (Ast[k] != k) return; g = ZM_sqr(gel(G,k)); if (gequal1(g) || gequalm1(g)) /* \pm Id */ str_printf(s, " node [midway] {$\\circ$}\n"); else str_printf(s, ";\n \\draw (start) arc (180:60:%.4f) node {$\\bullet$} arc (120:0:%.4f)\n",2*c/3,2*c/3); } static GEN polygon2tex(GEN V, GEN Ast, GEN G) { pari_sp av = avma; GEN v = Ast2v(Ast); pari_str s; long j, l = lg(V), flag; double c; str_init(&s, 1); flag = (l <= 16); str_puts(&s, "\n\\begin{tikzpicture}[scale=10]\n"); str_puts(&s, "\\draw (0,0.5)--(0,0) node [very near start, right] {$1^*$} node [below] {$0$}"); for (j=4; j < l; j++) { GEN a = M2Q(gel(V,j-1)), b = M2Q(gel(V,j)); c = gtodouble(gsub(b,a)) / 2; str_printf(&s, "node (start) {} arc (180:0:%.4f)\n", c); if (flag) { long sb = itos(numer_i(b)); long sa = itos(denom_i(b)); str_printf(&s, "node [midway, above] {%s} node [below]{$\\frac{%ld}{%ld}$}\n", (char*)gel(v,j-1), sb, sa); } maybe_decorate(&s,Ast,G,j-1,c); } c = (1- gtodouble(M2Q(gel(V,l-1)))) / 2; str_printf(&s, "node (start) {} arc (180:0:%.4f)\n", c); if (flag) str_printf(&s, "node [midway, above] {%s}", (char*)gel(v,l-1)); maybe_decorate(&s,Ast,G,l-1,c); str_printf(&s,"node [below] {$1$} -- (1,0.5) node [very near end, left] {$1$};"); str_printf(&s, "\n\\end{tikzpicture}"); return gerepileuptoleaf(av, strtoGENstr(s.string)); } static GEN circle2tex(GEN Ast, GEN G) { pari_sp av = avma; GEN v = Ast2v(Ast); pari_str s; long u, n = lg(Ast)-1; const double ang = 360./n; if (n > 30) { v = const_vec(n, (GEN)""); gel(v,1) = (GEN)"$(1,\\infty)$"; } str_init(&s, 1); str_puts(&s, "\n\\begingroup\n\ \\def\\geo#1#2{(#2:1) arc (90+#2:270+#1:{tan((#2-#1)/2)})}\n\ \\def\\sgeo#1#2{(#2:1) -- (#1:1)}\n\ \\def\\unarc#1#2#3{({#1 * #3}:1.2) node {#2}}\n\ \\def\\cell#1#2{({#1 * #2}:0.95) circle(0.05)}\n\ \\def\\link#1#2#3#4#5{\\unarc{#1}{#2}{#5}\\geo{#1*#5}{#3*#5}\\unarc{#3}{#4}{#5}}\n\ \\def\\slink#1#2#3#4#5{\\unarc{#1}{#2}{#5}\\sgeo{#1*#5}{#3*#5}\\unarc{#3}{#4}{#5}}"); str_puts(&s, "\n\\begin{tikzpicture}[scale=4]\n"); str_puts(&s, "\\draw (0, 0) circle(1);\n"); for (u=1; u <= n; u++) { if (Ast[u] == u) { str_printf(&s,"\\draw\\unarc{%ld}{%s}{%.4f}; \\draw\\cell{%ld}{%.4f};\n", u, v[u], ang, u, ang); if (ZM_isscalar(gpowgs(gel(G,u),3), NULL)) str_printf(&s,"\\fill \\cell{%ld}{%.4f};\n", u, ang); } else if(Ast[u] > u) str_printf(&s, "\\draw \\%slink {%ld}{%s}{%ld}{%s}{%.4f};\n", Ast[u]-u==n/2? "s": "", u, v[u], Ast[u], v[Ast[u]], ang); } str_printf(&s, "\\end{tikzpicture}\\endgroup"); return gerepileuptoleaf(av, strtoGENstr(s.string)); } GEN mspolygon(GEN M, long flag) { pari_sp av = avma; struct siegel T; long i, l; GEN v, G, msN; if (typ(M) == t_INT) { long N = itos(M); if (N <= 0) pari_err_DOMAIN("msinit","N", "<=", gen_0,M); msN = msinit_N(N); } else { checkms(M); msN = get_msN(M); } if (flag < 0 || flag > 3) pari_err_FLAG("mspolygon"); if (ms_get_N(msN) == 1) { GEN S = mkS(); T.V = mkvec2(matid(2), S); T.Ast = mkvecsmall2(1,2); G = mkvec2(S, mkTAU()); } else { siegel_init(&T, msN); l = lg(T.V); if (flag & 1) { long oo2 = 0; mssiegel(&T); for (i = 1; i < l; i++) { GEN c = gel(T.V, i); GEN c22 = gcoeff(c,2,2); if (!signe(c22)) { oo2 = i; break; } } if (!oo2) pari_err_BUG("mspolygon"); siegel_perm0(&T, rotate_perm(l, oo2)); } G = cgetg(l, t_VEC); for (i = 1; i < l; i++) gel(G,i) = get_g(&T, i); } if (flag & 2) v = mkvec5(T.V, T.Ast, G, polygon2tex(T.V,T.Ast,G), circle2tex(T.Ast,G)); else v = mkvec3(T.V, T.Ast, G); return gerepilecopy(av, v); } #if 0 static int iselliptic(GEN Ast, long i) { return i == Ast[i]; } static int isparabolic(GEN Ast, long i) { long i2 = Ast[i]; return (i2 == i+1 || i2 == i-1); } #endif /* M from msinit, F QM maximal rank */ GEN mslattice(GEN M, GEN F) { pari_sp av = avma; long i, ivB, j, k, l, lF; GEN D, U, G, A, vB, m, d; checkms(M); if (!F) F = gel(mscuspidal(M, 0), 1); else { if (is_Qevproj(F)) F = gel(F,1); if (typ(F) != t_MAT) pari_err_TYPE("mslattice",F); } lF = lg(F); if (lF == 1) return cgetg(1, t_MAT); D = mspolygon(M,0); k = msk_get_weight(M); F = vec_Q_primpart(F); if (typ(F)!=t_MAT || !RgM_is_ZM(F)) pari_err_TYPE("mslattice",F); G = gel(D,3); l = lg(G); A = gel(D,2); vB = cgetg(l, t_COL); d = mkcol2(gen_0,gen_1); m = mkmat2(d, d); for (i = ivB = 1; i < l; i++) { GEN B, vb, g = gel(G,i); if (A[i] < i) continue; gel(m,2) = SL2_inv2(g); vb = mseval(M, F, m); if (k == 2) B = vb; else { long lB; B = RgXV_to_RgM(vb, k-1); /* add coboundaries */ B = shallowconcat(B, RgM_Rg_sub(RgX_act_Gl2Q(g, k), gen_1)); /* beware: the basis for RgX_act_Gl2Q is (X^(k-2),...,Y^(k-2)) */ lB = lg(B); for (j = 1; j < lB; j++) gel(B,j) = vecreverse(gel(B,j)); } gel(vB, ivB++) = B; } setlg(vB, ivB); vB = shallowmatconcat(vB); if (ZM_equal0(vB)) return gerepilecopy(av, F); (void)QM_ImQ_hnfall(vB, &U, 0); if (k > 2) U = rowslice(U, 1, lgcols(U)-k); /* remove coboundary part */ U = Q_remove_denom(U, &d); F = ZM_hnf(ZM_mul(F, U)); if (d) F = RgM_Rg_div(F, d); return gerepileupto(av, F); } /**** Petersson scalar product ****/ /* oo -> g^(-1) oo */ static GEN cocycle(GEN g) { return mkmat22(gen_1, gcoeff(g,2,2), gen_0, negi(gcoeff(g,2,1))); } /* C = vecbinome(k-2) */ static GEN bil(GEN P, GEN Q, GEN C) { long i, n = lg(C)-2; /* k - 2 */ GEN s; if (!n) return gmul(P,Q); if (typ(P) != t_POL) P = scalarpol_shallow(P,0); if (typ(Q) != t_POL) Q = scalarpol_shallow(Q,0); s = gen_0; for (i = 0; i <= n; i++) { GEN t = gdiv(gmul(RgX_coeff(P,i), RgX_coeff(Q, n-i)), gel(C,i+1)); s = odd(i)? gsub(s, t): gadd(s, t); } return s; } static void mspetersson_i(GEN W, GEN F, GEN G, GEN *pvf, GEN *pvg, GEN *pC) { GEN WN = get_msN(W), annT2, annT31, section, c, vf, vg; long i, n1, n2, n3; annT2 = msN_get_annT2(WN); annT31 = msN_get_annT31(WN); section = msN_get_section(WN); if (ms_get_N(WN) == 1) { vf = cgetg(3, t_VEC); vg = cgetg(3, t_VEC); gel(vf,1) = mseval(W, F, gel(section,1)); gel(vf,2) = gneg(gel(vf,1)); n1 = 0; } else { GEN singlerel = msN_get_singlerel(WN); GEN gen = msN_get_genindex(WN); long l = lg(gen); vf = cgetg(l, t_VEC); vg = cgetg(l, t_VEC); /* generators of Delta ordered as E1,T2,T31 */ for (i = 1; i < l; i++) gel(vf, i) = mseval(W, F, gel(section,gen[i])); n1 = ms_get_nbE1(WN); /* E1 */ for (i = 1; i <= n1; i++) { c = cocycle(gcoeff(gel(singlerel,i),2,1)); gel(vg, i) = mseval(W, G, c); } } n2 = lg(annT2)-1; /* T2 */ for (i = 1; i <= n2; i++) { c = cocycle(gcoeff(gel(annT2,i), 2,1)); gel(vg, i+n1) = gmul2n(mseval(W, G, c), -1); } n3 = lg(annT31)-1; /* T31 */ for (i = 1; i <= n3; i++) { GEN f; c = cocycle(gcoeff(gel(annT31,i), 2,1)); f = mseval(W, G, c); c = cocycle(gcoeff(gel(annT31,i), 3,1)); gel(vg, i+n1+n2) = gdivgs(gadd(f, mseval(W, G, c)), 3); } *pC = vecbinome(msk_get_weight(W) - 2); *pvf = vf; *pvg = vg; } /* Petersson product on Hom_G(Delta_0, V_k) */ GEN mspetersson(GEN W, GEN F, GEN G) { pari_sp av = avma; GEN vf, vg, C; long k, l, tG, tF; checkms(W); if (!F) F = matid(msdim(W)); if (!G) G = F; tF = typ(F); tG = typ(G); if (tF == t_MAT && tG != t_MAT) pari_err_TYPE("mspetersson",G); if (tG == t_MAT && tF != t_MAT) pari_err_TYPE("mspetersson",F); mspetersson_i(W, F, G, &vf, &vg, &C); l = lg(vf); if (tF != t_MAT) { /* , two symbols */ GEN s = gen_0; for (k = 1; k < l; k++) s = gadd(s, bil(gel(vf,k), gel(vg,k), C)); return gerepileupto(av, s); } else if (F != G) { /* <(f_1,...,f_m), (g_1,...,g_n)> */ long iF, iG, lF = lg(F), lG = lg(G); GEN M = cgetg(lG, t_MAT); for (iG = 1; iG < lG; iG++) { GEN c = cgetg(lF, t_COL); gel(M,iG) = c; for (iF = 1; iF < lF; iF++) { GEN s = gen_0; for (k = 1; k < l; k++) s = gadd(s, bil(gmael(vf,k,iF), gmael(vg,k,iG), C)); gel(c,iF) = s; /* M[iF,iG] = */ } } return gerepilecopy(av, M); } else { /* <(f_1,...,f_n), (f_1,...,f_n)> */ long iF, iG, n = lg(F)-1; GEN M = zeromatcopy(n,n); for (iG = 1; iG <= n; iG++) for (iF = iG+1; iF <= n; iF++) { GEN s = gen_0; for (k = 1; k < l; k++) s = gadd(s, bil(gmael(vf,k,iF), gmael(vg,k,iG), C)); gcoeff(M,iF,iG) = s; /* */ gcoeff(M,iG,iF) = gneg(s); } return gerepilecopy(av, M); } }