\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ \\ The program modsym.gp \\ \\ public functions: modsym, xpm, T \\ \\ Timestamp: 2002/07/01 16:02:04 \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ \\ HISTORY: \\ \\ *) Original program written for PARI 1.x by Joseph L. Wetherell \\ \\ *) July 1999: William Stein updated the program for PARI 2.x. \\ \\ *) March 2001: Dominique Bernardi and Bernadette Perrin-Riou: \\ \\ - compute +/- part of homology (+ remove trick to halve matrix \\ size: some relationns were missing [e.g N=6, 106]) \\ \\ - added ellsym + xpm. \\ (compute components x^\pm of the modular symbol associated to E) \\ \\ *) July 2002: Karim Belabas and Bernadette Perrin-Riou: \\ cleanup + general speedup \\ \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ \\ Notation: N is the conductor (N <= 500 is doable in a < 2 minutes, > 600 \\ will require a _lot_ of memory ) global(talk = 1); \\ verbosity flag. global(hash, defaultmodsyms = 0, defaultellsyms = 0); modsym()=0 \\ forward declaration for addhelp xpm()=0 T()=0 { addhelp(modsym, "modsym(N, {sign = 1}): N integer > 1. Sets technical data associated to the homology of X_0(N), computed via modular symbols [+ part if sign = 1, - part otherwise]. This data is stored in global variables, used by xpm() and T(). N is allowed to be an elliptic curve E in ellinit form, in which case set data associated to the quotient corresponding to E in the (plus or minus part of the) homology of X_0(cond(E))"); addhelp(xpm, "xpm(a,b): a,b integers. Needs global data set by modsym(E, sign), E elliptic curve (see modsym). Returns the image of the path (0, a/b) in H_1(E, Q) into H_1(E, Q) / H_1(E, Q)^(-sign), in the basis given by a generator of the image of H_1(E, Z) in the quotient (canonical up to sign)"); addhelp(T, "T(p): p prime number. Needs global data set by modsym(N, sign). Returns matrix giving the action of the Hecke operator T_p (U_p if p | N) on H_1(X_0(N)) on some basis"); } \\ modsym struct x.cond = x[1] x.symb = x[2] x.hash = x[3] x.symquo= x[4][1] x.quosym= x[4][2] x.GEN = x[4][3] x.hbase = x[5][1] \\ H1 as sym x.hdual = x[5][2] \\ ellsym struct \\ .cond, .hash also usable x.hdualell = x[4] x.signe = x[5] x.ell = x[6] \\ 'MACRO' local to generatemsymbolds: uses v,curn,c,d,ret,N isnew() = { for (i=v,curn, \\ i local if ((d*ret[i][1] - c*ret[i][2])%N == 0, return (0)) ); 1 } generatemsymbols(N) = { local(num,ret,curn,divN,v,listp); \\ num = [Gamma:Gamma_0(N)] = N * Prod_{p|N} (1+p^-1) listp = factor(N)[,1]; num = N * prod(n=1,length(listp), 1+1/listp[n]); \\ initialize M-symbol list ret = vector(num); curn = 0; \\ generate M-symbols in three lists: \\ list 1: (c:1) for 0 <= c < N for(c=0, N-1, curn++; ret[curn] = [c,1]; ); \\ list 2: (1:d) for 0 <= d < N and gcd(d,N)>1 for(d=0, N-1, if (gcd(d,N) > 1, curn++; ret[curn] = [1,d]; ); ); \\ list 3: (c:d) with c|N (c != 1,N), gcd(c,d) = 1, gcd(d,N) > 1 v = curn+1; divN = divisors(N); for (d=2, N-1, if (gcd(d,N) > 1, for (i=2, length(divN)-1, \\ omit 1 (first) and N (last) c = divN[i]; if (gcd(c,d)==1 && isnew(), curn++; ret[curn] = [c,d]; ) ) ) ); if (curn != num, error("not the right number of symbols")); ret; } inithashmsymbols(N, msymbols)= { local(s); hash = matrix(N,N); for (n=1, N-1, if (gcd(n,N) == 1, for (k=1, length(msymbols), s = msymbols[k]; hash[(s[1]*n)%N + 1, (s[2]*n)%N + 1] = k; ) ) ) } dohash(N,v)= hash[v[1]%N + 1,v[2]%N + 1]; doV(v) = [-v[1], v[2]] doS(v) = [-v[2], v[1]] doT1(v) = [-v[1]-v[2], v[1]] doT2(v) = [v[2], -v[2]-v[1]] msymbolrelations(N, msymbols, signe) = { local(aux,nsym,n1,n2,n3,s1,relat,gen,quosym,symquo); nsym = length(msymbols); relat = matrix(nsym,3*nsym); for(n1=1,nsym, s1 = msymbols[n1]; n2 = dohash(N, doS(s1)); relat[n1,3*n1-2] = 1; \\ x + sigma(x) = 0 relat[n2,3*n1-2]++; n2 = dohash(N, doT1(s1)); n3 = dohash(N, doT2(s1)); relat[n1,3*n1-1] = 1; \\ x + tau(x) + tau^2(x) = 0 relat[n2,3*n1-1]++; relat[n3,3*n1-1]++; n2 = dohash(N, doV(s1)); relat[n1,3*n1] = 1; \\ x + eps(x) = 0 relat[n2,3*n1] -= signe; ); /* #if 1 */ relat = matrixqz(relat,-2); \\ very costly rank = nsym - length(relat); aux = mattranspose(mathnf(mattranspose(relat),1)[2]); /* #else *//* \\ 1) requires pari-2.2.0 or more [aux^(-1) later, aux non square] \\ 2) would allow much larger dimensions (hence conductors) _assuming_ we \\ install() and use mathnfspec() relat = mathnf(relat); aux = matsnf(concat(relat, matrix(nsym, nsym-length(relat))), 1); D = aux[3]; for (i=1, length(D), if (D[i,i], rank = i-1; break) ); aux = aux[1]; /* #endif */ symquo = vecextract(aux, Str(".."rank), ".."); quosym = vecextract(aux^(-1), Str(".."rank)); gen = listcreate(nsym); for (n=1, nsym, if (quosym[n,], listput(gen, n)) ); gen = Vec(gen); quosym = vecextract(quosym, gen, ".."); [symquo, quosym, gen] } iscuspeq(N,cusp1,cusp2) = { local(p1,q1,p2,q2,s1,s2,n); p1 = cusp1[1]; p2 = cusp2[1]; q1 = cusp1[2]; q2 = cusp2[2]; s1 = if(q1>2, lift(1/Mod(p1,q1)), 1); s2 = if(q2>2, lift(1/Mod(p2,q2)), 1); (s1*q2-s2*q1) % gcd(q1*q2,N) == 0; } \\ delta({c,d}) = [a/c] - [b/d] where ad-bc = 1 delta(msymbol)= { local(v); v = bezout(msymbol[1],msymbol[2]); if(v[3]<>1, error("msymbol not coprime:", msymbol)); [ [-v[1], msymbol[2]], [v[2], msymbol[1]] ] } \\ return [num(w), s], w = v (s=1) or w = -v (s=-1) \\ append w to cusps if needed getcusp(N,v,signe) = { local(n,mv); mv = [-v[1],v[2]]; for(n=1,length(cusps), if (iscuspeq(N,v, cusps[n]), return([n, 1])); if (iscuspeq(N,mv,cusps[n]), return([n, signe])); ); if ((signe > 0) || !iscuspeq(N,v,mv), listput(cusps, v); [length(cusps), 1]; , [1,0] ); } nbcusp(N) = sumdiv(N, d, eulerphi( gcd(d,N/d) )) kerdelta(N,msymbols,symquo,quosym,gen,signe) = { local(c,t,m,k, ngen,ncusp); ngen = length(gen); ncusp = nbcusp(N); cusps = listcreate(ncusp); m = matrix(ncusp,ngen); for (n=1, ngen, t = delta(msymbols[gen[n]]); c1 = getcusp(N,t[1],signe); c2 = getcusp(N,t[2],signe); m[c1[1], n] = c1[2]; m[c2[1], n] -= c2[2]; ); ncusp = length(cusps); if (!ncusp, return ([quosym,symquo])); m = vecextract(m, Str("1.."ncusp), ".."); m = matker(m * quosym, 1); \\ = Ker: Quo --- delta --> Cusps k = quosym*m; \\ = Ker: Symb --delta|Quo--> Cusps if (k, [k, m^(-1) * symquo] , [[;], [;]] ) } modulartosym(N,nsym,v)= { local(q,n1,sym,mods); if (!v[2], v=[1,N]); q = contfrac(v[1]/v[2]); q[1] = 1; for(n=3,length(q), q[n] = q[n]*q[n-1] + q[n-2]; ); sym = vectorv(nsym); for (n=2,length(q), mods = [(-1)^n * q[n], q[n-1]]; n1 = dohash(N,mods); sym[n1]++; ); sym } \\ needs hash TpMatrix(N,msymbols,gen,hbase,hdual,p) = { local(aux,nsym,cusp); if (!isprime(p) || N%p == 0, error("TpMatrix needs p prime and coprime to N."); ); nsym = length(msymbols); aux = matrix(nsym, length(gen)); for (c=1, length(gen), cusp = delta(msymbols[gen[c]]); aux[,c] = modulartosym(N,nsym,cusp[2]*[p,0;0,1]) - modulartosym(N,nsym,cusp[1]*[p,0;0,1]) + sum(n=0, p-1, modulartosym(N,nsym,cusp[2]*[1,0;n,p]) - modulartosym(N,nsym,cusp[1]*[1,0;n,p])); ); hdual * aux * hbase; } \\ needs hash UpMatrix(N,msymbols,gen,hbase,hdual,p) = { local(aux,nsym,cusp); if(!isprime(p) || N%p, error("UpMatrix needs p prime dividing N.") ); nsym = length(msymbols); aux = matrix(nsym, length(gen)); for (c=1, length(gen), cusp = delta(msymbols[gen[c]]); aux[,c] = sum(n=0, p-1, modulartosym(N,nsym,cusp[2]*[1,0;n,p]) - modulartosym(N,nsym,cusp[1]*[1,0;n,p])); ); hdual * aux * hbase; } modsym(N, signe = 1) = { local(symbols, relations, hdata,t1,t2,t3,t4); if (signe != 1 && signe != -1, error("modsym: sign not in {-1,1}")); if (type(N) != "t_INT", return (ellsym(N, signe))); if (N < 2, error("modsym: needs N > 1")); gettime(); if (talk, print1("1. Generating M-symbols (")); symbols = generatemsymbols(N); if (talk, print(t1=gettime()," ms)"); print1("2. Hashing M-symbols (")); gettime(); inithashmsymbols(N,symbols); if (talk, print(t2=gettime()," ms)"); print1("3. Quotienting out by relations (")); gettime(); relations = msymbolrelations(N,symbols,signe); if (talk, print(t3=gettime()," ms)"); print1("4. Computing Ker(delta) (")); gettime(); hdata = kerdelta(N,symbols,relations[1],relations[2],relations[3],signe); if (talk, print(t4=gettime()," ms)")); if (talk, print("Total time....................... (", t1+t2+t3+t4," ms)")); if (talk, print("genus is "length(hdata[1]))); defaultmodsyms = [N,symbols,hash,relations,hdata,signe]; return /* void */ } T(p, M) = { if (type(M)=="t_INT", M = defaultmodsyms; if (!M, error("T: no default space")) ); if (length(M) != 6, error("T: invalid modular symbols object")); hash = M.hash; gettime(); t = if (M.cond % p, TpMatrix(M.cond,M.symb,M.GEN,M.hbase,M.hdual,p); , UpMatrix(M.cond,M.symb,M.GEN,M.hbase,M.hdual,p); ); if (talk, print(gettime()," ms)")); t } ellsym(ell, signe = 1) = { local(cond, M, r, temp, res); cond = ellglobalred(ell)[1]; modsym (cond, signe); r = length(defaultmodsyms.hbase); res = if (r > 1, temp = [;]; forprime (p = 2, 10000, if (matrank (temp) == r-1, break); if (cond % p, temp = concat(temp, T(p) - ellap(ell, p))) ); matker(temp~)~ * defaultmodsyms.hdual , defaultmodsyms.hdual ); hash = defaultmodsyms.hash; defaultellsyms = [cond, defaultmodsyms.hbase, defaultmodsyms.hash, res[1,], signe, ell]; return /* void */ } /* \\ much slower (GP bug!) xpm(a,b) = { defaultellsyms.hdualell * modulartosym(defaultellsyms.cond, length(defaultellsyms.hdualell), [a,b]); } */ xpm(a,b) = { defaultellsyms[4] * modulartosym(defaultellsyms[1], length(defaultellsyms[4]), [a,b]); } f(N,p) = modsym(N); print(factor(charpoly(T(p))));