\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ \\ \\ The program modsym.gp \\ \\ \\ 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 made \\ the modifications described below: \\ \\ Les modifications sont de type : dans le programme modsym.gp, \\ les relations binaires n'étaient pas toutes prises en compte, \\ par exemple, certaines relations du type 2*x = 0 (voir par \\ exemple le cas N=6, N=106). Nous avons donc dû enlever \\ l'astuce qui permettait de réduire les matrices de taille \\ N à une taille N/2, ce qui ralentit bien sûr le programme. \\ Nous avons d'autre part rajouté un signe +1 ou -1 afin \\ de calculer aussi l'application de l'espace des formes \\ modulaires dans la partie + ou - de l'homologie. \\ Enfin, nous avons rajouté deux procédures ellsym et \\ xpm qui permettent de calculer les composantes x^\pm E \\ du symbole modulaire associée à une courbe elliptique \\ donnée dans GP/Pari par ellinit ou plutôt à la forme \\ modulaire asssociée. Le but est de calculer des valeurs \\ de fonctions L p-adiques associées à une courbe elliptique \\ et à un discriminant de corps quadratiques (positif ou négatif). \\ \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ talk=1; \\ verbosity flag. {mnormal(cond,v)= lift(v*Mod(1,cond)); } {generatemsymbols(cond, num,ret,curn)= \\ num = [\Gamma:\Gamma_0(N)] = N * Prod_{p|cond} (1+p^-1) num = cond; v = factor(cond); for(n=1,length(v~), num = num*(1+1/v[n,1]); ); \\ initialize M-symbol list ret = vector(num,c,0); curn = 1; \\ generate M-symbols in three lists: \\ list 1: (c:1) for 0 <= c < N for(c=0, cond-1, ret[curn] = [c,1]; curn = curn + 1; ); \\ list 2: (1:d) for 0 <= d < N and gcd(d,N)>1 for(d=0, cond-1, if(gcd(d,cond)>1, ret[curn] = [1,d]; curn = curn + 1; ,); ); \\ list 3: (c:d) with c|N, c<>1,N, gcd(c,d)=1, gcd(d,N)>1 v = curn; for(d=2, cond-1, if(gcd(d,cond)>1, fordiv(cond,c, if(c<>1 && c<>cond && gcd(c,d)==1, if(0 == sum(n=v,curn-1, ((ret[n][1]*d-c*ret[n][2])%cond)==0 ), ret[curn] = [c,d]; curn = curn + 1; ,); ,); ); ,); ); if (curn -num - 1, print("pas le bon nombre de symboles")); ret; } {hashmsymbols(cond,msymbols, ret,tmp) = ret = matrix(cond,cond,c1,d1,0); for(n=1,cond-1, if(gcd(n,cond)==1, for(k=1,length(msymbols), tmp = [msymbols[k][1]*n,msymbols[k][2]*n]; tmp = mnormal(cond,tmp); ret[tmp[1]+1,tmp[2]+1] = k; ); ,); ); ret; } {dohash(cond,hash,v, tmp)= v=mnormal(cond,v); hash[v[1]+1,v[2]+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(cond,msymbols,hash,signe, mask1,n1,n2,n3,s1,s2,s3,relat,nsym,freebasis,symtofree, tempr) = nsym = length(msymbols); relat = matid(nsym); mask1 = vector(nsym,n,1); for(n=1,nsym, n1 = msymbols[n]; n2 = doS(msymbols[n]); n1 = dohash(cond,hash,n1); n2 = dohash(cond,hash,n2); r = relat[n1,] + relat[n2,]; minc = vecmax(abs(r)); for(k=1,nsym, if(r[k] && abs(minc)>=abs(r[k]), minc = r[k]; mincloc = k; ,); ); if(minc, mask1[mincloc] = 0; r = r/minc; relat = relat - relat[,mincloc]*r; ,); ); for(n=1,nsym, n1 = msymbols[n]; n2 = doT1(msymbols[n]); n3 = doT2(msymbols[n]); n1 = dohash(cond,hash,n1); n2 = dohash(cond,hash,n2); n3 = dohash(cond,hash,n3); r = relat[n1,] + relat[n2,] + relat[n3,]; minc = vecmax(abs(r)); for(k=1,nsym, if(r[k] && abs(minc)>=abs(r[k]), minc = r[k]; mincloc = k; ,); ); if(minc, mask1[mincloc] = 0; r = r/minc; relat = relat - relat[,mincloc]*r; ,); ); for(n=1,nsym, n1 = msymbols[n]; n2 = doV(msymbols[n]); n1 = dohash(cond,hash,n1); n2 = dohash(cond,hash,n2); r = relat[n1,] -signe*relat[n2,]; minc = vecmax(abs(r)); for(k=1,nsym, if(r[k] && abs(minc)>=abs(r[k]), minc = r[k]; mincloc = k; ,); ); if(minc, mask1[mincloc] = 0; r = r/minc; relat = relat - relat[,mincloc]*r; ,); ); n1 = vector(nsym,n,mask1[n]*n); n2 = eval(Set(simplify(n1))); freebasis = vector(length(n2)-1,n,n2[n+1]); symtofree = matrix(nsym,length(freebasis),n1,n2,relat[n1,freebasis[n2]]); [freebasis,symtofree~]; } cusps; {initcusps() = cusps=[];} {testcuspeq(cond,cusp1,cusp2, p1,q1,p2,q2,s1,s2,n)= p1 = cusp1[1]; q1 = cusp1[2]; p2 = cusp2[1]; q2 = cusp2[2]; s1 = if(q1>2, lift(Mod(1,q1)/Mod(p1,q1)), 1); s2 = if(q2>2, lift(Mod(1,q2)/Mod(p2,q2)), 1); foo=Mod(s1*q2-s2*q1, gcd(q1*q2,cond))==0; foo; } {getcusp(cond,v,signe, ret,negv) = ret = 0; negv = [-v[1],v[2]]; for(n=1,length(cusps), if(testcuspeq(cond,v,cusps[n][1]), ret=[n,cusps[n][2]], if(testcuspeq(cond,negv,cusps[n][1]), ret=[n,signe*cusps[n][2]])); ); if(ret==0, cusps = concat(cusps, [[v, 1-((signe < 0) && testcuspeq(cond,v,negv))]]); ret = length(cusps); ret = [ret,cusps[ret][2]]; ,); ret; } {delta(cond,msymbol, v,c1,c2,ret)= v = bezout(msymbol[1],msymbol[2]); if(v[3]<>1, print("msymbol not coprime:", msymbol);1/0 ,); c1 = [v[2],msymbol[1]]; c2 = [-v[1],msymbol[2]]; ret=[c2/gcd(c2[1],c2[2]), c1/gcd(c1[1],c1[2])]; ret; } {kerdelta(cond,msymbols,hash,freebasis,signe, rels,trels,m,k,kinv,ksz)= initcusps(); rels = vector(length(freebasis),n, delta(cond,msymbols[freebasis[n]]) ); for(n=1,length(rels), rels[n] = [getcusp(cond,rels[n][1],signe),getcusp(cond,rels[n][2],signe)]; ); m = matrix(length(cusps),length(freebasis),n1,n2,0); for(n=1,length(rels), m[rels[n][1][1],n] = m[rels[n][1][1],n] + rels[n][1][2]; m[rels[n][2][1],n] = m[rels[n][2][1],n] - rels[n][2][2]; ); k = matker(m); ksz = matsize(k); kinvr = (matsupplement(k))^(-1); [k, vecextract(kinvr,2^ksz[2]-1,2^ksz[1]-1)]; } {modulartofree(cond,hash,symtofree,v, q,hashval,n1,s1,sym,mods,ret)= if(v[2]==0, v=[1,cond]); 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 = vector(matsize(symtofree)[2],n,0)~; for(n=2,length(q), mods = [(-1)^(n)*q[n],q[n-1],1]; n1 = dohash(cond,hash,mods); sym[n1] = sym[n1]+1; ); ret=symtofree*sym; ret; } {tpmatrix(cond,msymbols,hash,relations,hdata,p, freebasis,symtofree,hbase,hdual, sz,ret,col,sym,cusp,a)= if(!isprime(p) || cond%p == 0, print("Error, this program can only compute T_p for p prime and coprime to N."); return([]); ); freebasis = relations[1]; symtofree = relations[2]; hbase = hdata[1]; hdual = hdata[2]; sz = matsize(hbase); ret = matrix(sz[2],sz[2],n1,n2,0); for(c=1,sz[2], col = vector(sz[1],n,0)~; for(r=1,sz[1], a = hbase[r,c]; if(a, sym = msymbols[freebasis[r]]; cusp = delta(cond,sym); col = col - a*modulartofree(cond,hash,symtofree,cusp[1]*[p,0;0,1]); col = col + a*modulartofree(cond,hash,symtofree,cusp[2]*[p,0;0,1]); for(n=0,p-1, col = col - a*modulartofree(cond,hash,symtofree,cusp[1]*[1,0;n,p]); col = col + a*modulartofree(cond,hash,symtofree,cusp[2]*[1,0;n,p]); ); ,); ); ret[,c] = hdual*col; ); ret; } {gammap(cond,msymbols,hash,relations,hdata,p, freebasis,symtofree,hbase,hdual, sz,ret,col,sym,cusp,a)= freebasis = relations[1]; symtofree = relations[2]; hbase = hdata[1]; hdual = hdata[2]; sz = matsize(hbase); col = vector(sz[1],n,0)~; for(n=1,p-1, col = col + modulartofree(cond,hash,symtofree,[n,p]); ); hdual*col; } defaultmodsymspace=0; {modsym(N,signe, symbols, hash, relations, hdata, t1,t2,t3,t4)= gettime(); if(talk,print1("1. Generating M-symbols \t\t(")); symbols = generatemsymbols(N); if(talk,print(t1=gettime()," ms)"); print1("2. Hashing M-symbols \t\t\t(")); gettime(); hash = hashmsymbols(N,symbols); if(talk,print(t2=gettime()," ms)"); print1("3. Quotienting out by relations \t(")); gettime(); relations = msymbolrelations(N,symbols,hash,signe); if(talk,print(t3=gettime()," ms)"); print1("4. Computing the kernel of delta \t(")); gettime(); hdata = kerdelta(N,symbols,hash,relations[1],signe); if(talk,print(t4=gettime()," ms)")); if(talk,print("Total time.......................\t(", t1+t2+t3+t4," ms)")); defaultmodsymspace=[N,symbols,hash,relations,hdata,signe]; } {T(p, M=0)= if(type(M)=="t_INT",M=defaultmodsymspace; if(M==0,print("T: Error -- no default space."))); if(matsize(M)[2] != 6, print("T: invalid modular symbols object.")); tpmatrix(M[1],M[2],M[3],M[4],M[5],p,M[6]); } {ellsym (ell,signe, modele, cond, M, r, temp, res)= modele = ellglobalred(ell); cond = modele [1]; M = modsym (cond, signe); r = matsize (M[5][1]) [2]; print("r = ", r); res = if (r > 1, temp = [;]; forprime (p = 2, 10000, if (matrank (temp) == r-1, break,); print ("p = ", p); if (cond % p, temp = concat (temp, T (p, M) - ellap(ell, p)*matid(r)),) ); mattranspose (matker (mattranspose (temp))) * M[5][2], M[5][2] ); defaultellsyms = [cond, M[3], M[4][2], res, signe, ell] } {xpm (a, b, parms)= if (parms,, parms = defaultellsyms); (parms [4] * modulartofree (parms [1], parms [2], parms [3], [a,b], parms [5]))[1]; } {f(N,p)= modsym(N,1); print(factor(charpoly(T(p),x))); modsym(N,-1); print(factor(charpoly(T(p),x))); }