/****-*-magma-*-************************************************************* * HECKE: Modular Symbols Calculator * * * * Version 0.5M (February 9, 2000) * * * *Send bug reports and suggestions to William Stein was@math.berkeley.edu. * ***************************************************************************/ /************************************************************************* TODO TODO TODO TODO TODO TODO TODO TODO TODO [ ] Intersections with Eisenstein series NOT programmed! [ ] Add a feature for "tensoring with Qp". [ ] Decomposition of S_3(25,[4];Q) does not work. [ ] Need to be able to compute eigenforms over a single number field, so they can be worked with. > m4:=msym(DirichletCharacter(13,[3]),4,+1); > Print(Decomposition(m4)); Modular symbols factors: 13k4A: dim = 1 cuspidal 13k4B: dim = 2 cuspidal 13k4C: dim = 1 eisenstein 13k4D: dim = 1 eisenstein > f4b:=qEigenform(Factor(m4,2),17); > MinimalPolynomial(Coefficient(f4b,2)); x^2 + (-5*zeta - 5)*x + 2*zeta > f4b; q + a*q^2 + (-3*a + (5*zeta + 5))*q^3 + ((5*zeta + 5)*a - 10*zeta)*q^4 + (-5*zeta*a - 20)*q^5 + ((-10*zeta - 10)*a + 6*zeta)*q^6 + ((zeta + 1)*a + 5*zeta)*q^7 + (7*zeta*a + 10)*q^8 + ((15*zeta + 15)*a - 20*zeta)*q^9 + (5*a - 10*zeta - 10)*q^10 + (-15*a + (29*zeta + 29))*q^11 + (-20*zeta*a + 20)*q^12 + ((-7*zeta - 5)*a + (50*zeta + 45))*q^13 + (10*zeta*a + 2)*q^14 + (10*a - 70*zeta - 70)*q^15 + (15*a - 66*zeta - 66)*q^16 + ((-16*zeta - 16)*a + 75*zeta)*q^17 USELESS! [ ] Torsion bound -- should start new bound where left off from old one. [ ] Analytic invariants, e.g., ImaginaryPeriod: don't get force recomputed when the precision of SOME analytic computation goes up and they are re-called on. They should! [ ] Clean interface!! Seperate files? *************************************************************************/ HeckeBugWatch := false; HeckeVerbose := true; HeckeSuperverbose := false; /************************************************************************ * * * DATA STRUCTURES * * * ************************************************************************/ CManSymList := recformat< k, // weight F, // base field R, // F[X,Y]. p1list, // List of elements of P^1(Z/NZ). n // size of list of Manin symbols = #p1list * (k-1). > ; /************************************************************************ * CQuotient: * * First we motivate the format. * * To save memory and space we mod out by the S, T, and possibly I * * relations in two phases. First mod out by the S,I relations * * x + xS = 0, x - xI = 0. * * by identifying equivalent pairs. Next mod out by the T relations * * x + xT + xT^2 = 0. * ************************************************************************/ CQuotient := recformat< Sgens, // * List of the free generators after modding out // by the 2-term S and possibly I relations. Squot, // * The i-th Manin symbol is equivalent to Scoef, // Scoef[i]*(Squot[i]-th free generator). Tgens, // * List of the free generators after modding // out the S-quotient by the T-relations. Tquot // * The i-th Sgen is equal to Tquot[i], which // is a vector on Tgens. >; // The categories: (someday these will be made a part of Magma.) // Presently they are used to implement polymorphism in some cases. ModSym := 1; ModSymFac := 2; function GetCategory(M) return M`theCategory; end function; procedure SetCategory(M, category) M`theCategory := category; end procedure; /////////////////////////////////////////////////////////////////// // ModSym: Modular symbols object. // /////////////////////////////////////////////////////////////////// declare attributes ModTupFld: theCategory, // * indicate that this is a mod sym object k, // * weight N, // * level F, // * the base field eps, // * the Dirichlet character sign, // * sign -- see doc of ModularSymbols() mlist, // * list of Manin symbols, and other related data quot, // * gives representation of any // Manin symbol in terms of free generators disc, // * discriminant of the Hecke algebra decomp, olddecomp, newdecomp, // * ambient decomposition semidecomp, gammaalp, // * action of character via Gamma_0(N). Z, // * integral modular symbols. ModSymBasis, // * Basis of modular symbols. BoundaryMap, // * boundary map SkF, // * cuspidal modular symbols, over F SkZ, // * cuspidal modular symbols, over integer Z. // (only makes sense when char(F)=0). SkZplus,SkZminus, // * for *-involution. old, // * old spaces sknew, sknewdual, // * new cuspidal modular symbols mknew, mknewdual, // * new cuspidal modular symbols newforms, // * basis of eigenforms for new subspace. oldforms, // * old forms, from lower level. winding, // * list of winding submodules T_images, // * images of standard basis vectors under Hecke. X, // * list of character groups for mestre, // each p exactly dividing N. cusplist; // * the cusps. /////////////////////////////////////////////////////////////////// // ModSymFac: Factor of a modular symbols object. // // This is represented by a subspace of the dual of // // a ModSym object. // /////////////////////////////////////////////////////////////////// declare attributes ModTupFld: theCategory, // * indicates that this is a modsym factor. M, // * ambient space of modular symbols. V, // * Subspace of M Z, // * Lattice -- Integral subspace of ambient space. Zdual, // * Lattice -- Integral subspace of ambient space. qexp, // * q-expansion -- a pair, f(x), coefs. eigen, eigenplus, eigenminus, // * Eigenvector in M tensor Fbar. eigenint, // * EigenformInTermsOfIntegralBasis moddeg, // * Degree of the canonical map. PeriodLattice, realvolume,// * Real volume imagvolume,// * Pure imaginary volume cinf, cinfimag, // * Number of real components. Lvals, // * Critical values of the L-function. LRatio, // * Ratios L(M,1)/Omega*(fudge). LRatioOdd, // * odd part of Ratios L(M,1)/Omega*(fudge). cp, // * Orders of component groups. wq, // * Atkin-Lehner qintbasis, // * integral basis of q-expansions. N, // * the level iso, // * the isogeny class isnew, // * whether or not it is new iscusp, // * whether it is cuspidal or not fast, // * fast matrix of Tn on M. signspace, // * A^+ and A^-. // * there is also the X from ModularSymbols... xdata, // * local phix and mx ZxZalp, VolLEven, VolLOdd, // * For LRatio. PeriodGens, PGfast, RatPeriodMap, RatPeriodLat, RatPeriodConj,RatPeriodMapSign, PeriodMap, PeriodMapPrecision; /************************************************************************ * * * Manin Symbols (basic functions) * * * ************************************************************************/ ///////////////////////////////////////////////////////////////////////// // NormalizePair (written by Allan Steel) // // INPUT: [u,v] in Z/N x Z/N. // // OUTPUT: 1) the *index* of a fixed choice [a,b] of representative // // in the P^1(Z/NZ) equivalence class of [u,v]. // // If gcd(u,v,N)>1 then the index returned is 0. // // 2) a scalar s such that // // [a*s,b*s] = [u,v]. // // The character relation is thus // // eps(s) [a,b] == [u,v]. // ///////////////////////////////////////////////////////////////////////// function NormalizePair(x) Z := IntegerRing(); u := x[1]; v := x[2]; R := Parent(u); N := Modulus(R); if u eq 0 then return [R | 0, 1], v; else u := Z ! u; g, s := XGCD(u, N); if g ne 1 then d := N div g; while GCD(s, N) ne 1 do s := (s + d) mod N; end while; end if; // Multiply (u, v) by s; new u = g u := R!g; v := v*s; minv := Z!v; mint := 1; if g ne 1 then Ng := N div g; vNg := v*Ng; t := R!1; for k := 2 to g do v +:= vNg; t +:= Ng; if Z!v lt minv and IsUnit(t) then minv := Z!v; mint := t; end if; end for; end if; if GCD([Z | u, minv, N]) ne 1 then printf "Bad x=%o, N=%o, u=%o, minv=%o, s=%o, mint=%o\n", x, N, u, minv, s, mint; error ""; end if; // return [R|u, minv], R!InverseMod(s*mint,N); return [R|u, minv], 1/(R!(s*mint)); end if; end function; ////////////////////////////////////////////////////////////////////////// // P1Reduce: // // INPUT: [u,v] in Z/N x Z/N. // // OUTPUT: 1) the *index* of a fixed choice [a,b] of representative // // in the P^1(Z/NZ) equivalence class of [u,v]. // // If gcd(u,v,N)>1 then the index returned is 0. // // 2) a scalar s such that // // [a*s,b*s] = [u,v]. // // so the character relation is // // eps(s) [a,b] == [u,v]. // ////////////////////////////////////////////////////////////////////////// function P1Reduce(x, list) N := Characteristic(Parent(x[1])); if N eq 1 then return 1, 1; end if; if Gcd([x[1], x[2], N]) ne 1 then return 0, 1; end if; n, s := NormalizePair(x); return Index(list, n), s; end function; ////////////////////////////////////////////////////////////////////////// // EnumerateP1: // // Construct a list of distinct normalized representatives for // // the set of equivalence classes of P^1(Z/NZ). // // The output of this function is used by P1Reduce. // ////////////////////////////////////////////////////////////////////////// function EnumerateP1(N) R := (N gt 1) select IntegerRing(N) else quo; return {@ NormalizePair([R|c,1]) : c in [0..N-1] @} join {@ NormalizePair([R|1,d]) : d in [0..N-1] | Gcd(d,N) gt 1 @} join {@ NormalizePair([R|c,d]) : c in Divisors(N), d in [0..N-1] | c ne 1 and c ne N and Gcd(c,d) eq 1 and Gcd(d,N) gt 1 @}; end function; ////////////////////////////////////////////////////////////////////////// // ManinSymbolList: // // Construct a list of distinct Manin symbols. These are // // elements of the Cartesion product: // // {0,...,k-2} x P^1(Z/NZ). // // In fact, we only store a list of the elements of P^1(Z/NZ), // // as the full Cartesion product can be stored using // // the following scheme. // // // // There are Manin symbols 1,...,#({0,..,k-2}xP^1) indexed by i. // // Thus i is an integer between 1 and the number of generating Manin // // symbols. Suppose P^1(N) = {x_1, ..., x_n}. The encoding is // // as follows: // // 1=, 2=<0, x_2>, ..., n=<0,x_n>, // // n+1=,n+2=<1, x_2>, ...,2n=<1,x_n>, // // ... // ////////////////////////////////////////////////////////////////////////// function ManinSymbolList(k,N,F) p1list := EnumerateP1(N); n := (k-1)*#p1list; return rec; end function; /////////////////////////////////////////////////////////////////////////// // ManinSymbolApply // // Apply an element g=[a,b, c,d] of SL(2,Z) to the i-th // // standard Manin symbol. The definition of the action is // // g.(X^i*Y^j[u,v]) := // // (aX+bY)^i*(cX+dY)^j [[u,v]*g]. // // The argument "manin" should be as returned by // // the function ManinSymbolsList above. // // The character eps should be an array which // // defines a Dirichlet character Z/NZ ---> K. // // Thus eps[a] in K and makes sense for a in (Z/NZ)^*. // // If eps is trivial, then the character is assumed trivial. // /////////////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////////////// // Unwind -- Given an int i, this function gives back the // // ith generating Manin symbol. // /////////////////////////////////////////////////////////////////////////// function UnwindManinSymbol(i, mlist) // P^1(N) part. p1list := mlist`p1list; n := #p1list; ind := (i mod n); if ind eq 0 then ind := n; end if; uv := p1list[ind]; w := Integers()!((i - ind)/n); return uv, w, ind; end function; // Wind -- Given w in the range 0<=w<=k-2, and an index ind // into the P^1-list, this function gives back // the number of the generating Manin symbol. function WindManinSymbol(w, ind, mlist) return w*#(mlist`p1list) + ind; end function; intrinsic ManinSymbolApply(g::SeqEnum, i::RngIntElt, mlist::Rec, eps::Rec) -> SeqEnum {Apply g to the ith Manin symbol.} k := mlist`k; uv, w, ind := UnwindManinSymbol(i,mlist); p1list := mlist`p1list; uvg := Parent(uv) ! [g[1]*uv[1]+g[3]*uv[2], g[2]*uv[1]+g[4]*uv[2]]; act_uv, scalar := P1Reduce(uvg, p1list); if act_uv eq 0 then return [<0,1>]; end if; if k eq 2 then if not IsTrivial(eps) then return []; /* if Evaluate(eps,scalar) ne 0 then return []; else return [<0,act_uv>]; end if; */ else return [<1,act_uv>]; end if; end if; // Polynomial part. R :=PolynomialRing(mlist`F); // univariate is fine for computation hP := (g[1]*x+g[2])^w*(g[3]*x+g[4])^(k-2-w); if not IsTrivial(eps)then hP *:= Evaluate(eps,scalar); end if; pol := ElementToSequence(hP); // Put it together n := #p1list; ans := [ : w in [0..#pol-1]]; return [x : x in ans | x[1] ne 0]; end intrinsic; ///////////////////////////////////////////////////////////////// // ModularSymbolApply // // Apply an element g=[a,b, c,d] of SL(2,Z) to the given // // modular symbol. The definition of the action is // // g.(X^i*Y^j{u,v}) := // // (dX-bY)^i*(-cX+aY)^j {g(u),g(v)}. // // A modular symbol is represented as a sequence of pairs // // , where P is a polynomial in X and Y (an element // // of the ring M`mlist`R, and x=[[a,b],[c,d]] is a pair // // of elements of P^1(Q), // // where [a,b] <--> a/b and [c,d] <--> c/d. // // After computing the action, no further reduction is done. // ///////////////////////////////////////////////////////////////// intrinsic ModularSymbolApply( M::ModTupFld, g::SeqEnum, s::Tup ) -> SeqEnum {Apply an element g=[a,b, c,d] of SL(2,Z) to the given modular symbol. The definition of the action is g.(X^i*Y^j\{u,v\}) := (dX-bY)^i*(-cX+aY)^j \{g(u),g(v)\}. A modular symbol is represented as a pair , where P is a polynomial in X and Y (an element of the ring M`mlist`R, and x=[[a,b],[c,d]] is a pair of elements of P^1(Q), where [a,b] <--> a/b and [c,d] <--> c/d.} R := M`mlist`R; h := hom R | g[4]*R.1-g[2]*R.2, -g[3]*R.1+g[1]*R.2>; hP := h(s[1]); a,b:= Explode(s[2][1]); c,d:= Explode(s[2][2]); gx := [[g[1]*a+g[2]*b, g[3]*a+g[4]*b], [g[1]*c+g[2]*d, g[3]*c+g[4]*d]]; return ; end intrinsic; intrinsic ModularSymbolApply(M::ModTupFld, g::SeqEnum, s::SeqEnum) -> SeqEnum {Apply an element g=[a,b, c,d] of SL(2,Z) to the given modular symbol. The definition of the action is g.(X^i*Y^j\{u,v\}) := (dX-bY)^i*(-cX+aY)^j \{g(u),g(v)\}. A modular symbol is represented as a sequence of pairs , where P is a polynomial in X and Y (an element of the ring M`mlist`R, and x=[[a,b],[c,d]] is a pair of elements of P^1(Q), where [a,b] <--> a/b and [c,d] <--> c/d.} return [ModularSymbolApply(M, g,term): term in s]; end intrinsic; ///////////////////////////////////////////////////////////////// // Sparse2Quotient: // // Performs Sparse Gauss elimination on a matrix all of // // whose columns have at most 2 nonzero entries. I just // // use the fairly obvious algorithm. It runs fast // // enough. Example: // // rels := [[1,2], [1,3], [2,3], [4,5]]; // // relc := [[3,-1],[1,1], [1,1], [1,-1]]; // // n := 5; // x1,...,x5 // // corresponds to 3*x1-x2=0, x1+x3=0, x2+x3=0, x4-x5=0. // ///////////////////////////////////////////////////////////////// function Sparse2Quotient (rels, relc, n, F) free := [1..n]; coef := [F|1 : i in [1..n]]; related_to_me := [[] : i in [1..n]]; for i in [1..#rels] do t := rels[i]; c1 := relc[i][1]*coef[t[1]]; c2 := relc[i][2]*coef[t[2]]; // Mod out by the relation // c1*x_free[t[1]] + c2*x_free[t[2]] = 0. die := 0; if c1 eq 0 and c2 eq 0 then // do nothing. elif c1 eq 0 and c2 ne 0 then // free[t[2]] --> 0 die := free[t[2]]; elif c2 eq 0 and c1 ne 0 then die := free[t[1]]; elif free[t[1]] eq free[t[2]] then if c1+c2 ne 0 then // all xi equal to free[t[1]] must now equal to zero. die := free[t[1]]; end if; else // x1 = -c2/c1 * x2. x := free[t[1]]; free[x] := free[t[2]]; coef[x] := -c2/c1; for i in related_to_me[x] do free[i] := free[x]; coef[i] *:= coef[x]; Append (~related_to_me[free[t[2]]], i); end for; Append (~related_to_me[free[t[2]]], x); end if; if die gt 0 then for i in related_to_me[die] do free[i] := 1; coef[i] := 0; end for; free[die] := 1 ; coef[die] := 0; end if; end for; // Enumerate the subscripts of free generators that survived. // x_{i_1} <----> y_1 // x_{i_2} <----> y_2, etc. for i in [1..#free] do if coef[i] eq 0 then free[i] := -1; end if; end for; ykey := {@ x : x in Set(free) | x ne -1 @}; for i in [1..#free] do if free[i] eq -1 then free[i] := 1; else free[i] := Index(ykey,free[i]); end if; end for; return ykey, free, coef; end function; /********************************************************* Quotient: The INPUT is a list of (sparse) relations [[, , ... ], ... ] The first listed spare relation is c_1*e_{i_1} + c_2*e_{i_2} + ... c_r*e_{i_r} = 0. The integer n denotes the total number of basis elements e_i. The field K is the field over which the c_i are defined. The OUTPUT is (1) an indexed set of g free generators, and (2) an expression for each e_i in terms of the free generators. These expressions are elements of the g-dimensional vector space over K. generators, quotient *********************************************************/ function Pivots(B) // find pivots of reduced basis. pivots := []; for b in B do j := 1; while b[j] eq 0 do j +:= 1; end while; Append(~pivots, j); end for; return pivots; end function; forward Quotient2; function Quotient(rels, n, K) ringchoice := Characteristic(K) eq 0; if Type(K) eq FldPr then ringchoice := false; end if; return Quotient2(rels,n,K : ring:=ringchoice); end function; function Quotient2(rels, n, K : ring := true) if HeckeSuperverbose then "# relations = ", #rels; end if; if ring then // POLYNOMIAL ALGORITHM if HeckeSuperverbose then "ring reduction"; end if; R := PolynomialRing(K, n); vars := {@ R.i : i in [1..n] @}; Rels := [ &+[t[1]*R.t[2] : t in r] : r in rels]; Rels := Reduce(Rels); W := quo; gens := {@ i : i in [1..n] | not exists { g : g in Rels | LeadingMonomial(g) eq R.i } @}; Z := VectorSpace(K,#gens); quot := [ &+[Z| LeadingCoefficient(t) * Z.(Index(gens,Index(vars,t/LeadingCoefficient(t)))) : t in Terms(R!W.i)] : i in [1..n] ]; else // VECTOR SPACE ALGORITHM if HeckeSuperverbose then "vector space reduction"; end if; V := VectorSpace(K, n); Rels := [ &+[t[1]*V.t[2] : t in r] : r in rels]; S := RowSpace(RMatrixSpace(K, #Rels, n) ! Rels); m := n - Dimension(S); B := Basis(S); W := VectorSpace(K, m); quot := [ W!0 : i in [1..n]]; pivots := Pivots(B); gens := {@ i : i in [1..n] | i notin pivots @}; for i := 1 to #pivots do quot[pivots[i]] := -W![B[i][j] : j in gens]; end for; for i :=1 to #gens do quot[gens[i]] := W.i; end for; end if; return gens, quot; end function; //////////////////////////////////////////////////////////////////// // CONSTRUCT QUOTIENT BY 2 AND 3 TERM RELS // //////////////////////////////////////////////////////////////////// function ManSym2termQuotient (mlist, eps, sign) n := mlist`n; K := mlist`F; S := [0,-1, 1,0]; I := [-1,0, 0,1]; xS := [ManinSymbolApply(S,i,mlist,eps) : i in [1..n]]; S_rels := [ [i, (xS[i][1])[2]] : i in [1..n]| #xS[i] gt 0]; S_relc := [ [K!1, (xS[i][1])[1]] : i in [1..n]| #xS[i] gt 0]; if sign ne 0 then xI := [ManinSymbolApply(I,i,mlist,eps) : i in [1..n]]; I_rels := [ [i, (xI[i][1])[2]] : i in [1..n]]; I_relc := [ [K!1, -sign*(xI[i][1])[1]] : i in [1..n]]; else I_rels := []; I_relc := []; end if; rels := S_rels cat I_rels; relc := S_relc cat I_relc; return Sparse2Quotient(rels,relc,n,K); end function; function ManSym3termQuotient (mlist, eps, Skey, Squot, Scoef) // Construct the subspace of 3-term relations. n := mlist`n; F := mlist`F; k := mlist`k; T := [0,-1, 1,-1]; TT:= [-1,1, -1,0]; if k eq 2 then mask := [false : i in [1..n]]; rels := []; for j in [1..n] do if not mask[j] then t := ManinSymbolApply(T,j,mlist,eps)[1]; tt := ManinSymbolApply(TT,j,mlist,eps)[1]; mask[t[2]] := true; mask[tt[2]] := true; Append(~rels, [, , ]); end if; end for; else rels := [&cat[ [], [ : x in ManinSymbolApply(T,i,mlist,eps)], [ : x in ManinSymbolApply(TT,i,mlist,eps)] ] : i in [1..n]]; end if; return Quotient(rels, #Skey, F); end function; /********************************************************* * * * MODULARSYMBOLS: * * Construct the space of modular symbols of weight k, * * level N and character eps. * * * *********************************************************/ intrinsic ModularSymbols(N::RngIntElt) -> ModTupFld {Returns the space of modular symbols of weight 2 for Gamma_0(N), over the rational numbers.} return ModularSymbols(DirichletCharacter(N),2,0); end intrinsic; intrinsic ModularSymbols(N::RngIntElt, k::RngIntElt) -> ModTupFld {Returns the space of modular symbols of weight k for Gamma_0(N), over the rational numbers.} return ModularSymbols(DirichletCharacter(N),k,0); end intrinsic; intrinsic ModularSymbols(N::RngIntElt, k::RngIntElt, sign::RngIntElt) -> ModTupFld {Returns the space of modular symbols of weight k for Gamma_0(N), over the rational numbers. If sign=+1 then returns the +1 quotient, if sign=-1 returns the -1 quotient, and if sign=0 returns the full space. If sign is a positive prime, then the full space of modular symbols mod sign isreturned.} if sign gt 1 and IsPrime(sign) then return ModularSymbols(N,k,0,sign); end if; return ModularSymbols(DirichletCharacter(N),k,sign); end intrinsic; intrinsic ModularSymbols(N::RngIntElt, k::RngIntElt, sign::RngIntElt, p::RngIntElt) -> ModTupFld {Returns the the space of modular symbols of weight k for Gamma_0(N), over the finite field Fp.} assert IsPrime(p); return ModularSymbols(DirichletCharacter(N,p),k,sign); end intrinsic; intrinsic ModularSymbols(eps::Rec) -> ModTupFld {Returns the space of modular symbols of weight k for Gamma_1(N) with character eps, over the rational numbers.} return ModularSymbols(eps,2,0); end intrinsic; intrinsic ModularSymbols(eps::Rec, k::RngIntElt) -> ModTupFld {Returns the space of modular symbols of weight k for Gamma_1(N) with character eps, over the rational numbers.} return ModularSymbols(eps,k,0); end intrinsic; forward TrivialSpace; intrinsic ModularSymbols(eps::Rec, k::RngIntElt, sign::RngIntElt) -> ModTupFld {Returns the space of modular symbols of character eps and weight k. The level and base field are specified by the DirichletCharacter eps. The third argument "sign" allows for working in certain quotients. The possible values are -1, 0, or +1, which correspond to the -1 quotient, full space, and +1 quotient.} if not (-1 le sign and sign le 1) then error "sign must be -1, 0, or 1"; end if; if k lt 2 then error "Modular symbols are only defined for weight k >= 2."; end if; N := CharacterLevel(eps); F := FieldOfValues(eps); if HeckeVerbose then tt := Cputime(); printf "Creating M_%o(Gamma_1(%o),eps;F_%o)", k, N, Characteristic(F) eq 0 select 0 else #F; if sign eq 1 then printf "^+"; elif sign eq -1 then printf "^-"; end if; if not HeckeSuperverbose then printf "\n"; end if; end if; if Evaluate(eps,-1) ne (-1)^k then return TrivialSpace(k,eps,sign); end if; if HeckeSuperverbose then t := Cputime(); "Step I. Manin symbols list."; end if; mlist := ManinSymbolList(k,N,F); if HeckeSuperverbose then Cputime(t), " seconds."; end if; if HeckeSuperverbose then t := Cputime(); "Step II. 2-term relations."; end if; Sgens, Squot, Scoef := ManSym2termQuotient(mlist, eps, sign); if HeckeSuperverbose then Cputime(t), " seconds."; end if; if #Sgens lt 1 then return TrivialSpace(k,eps,sign); end if; if HeckeSuperverbose then t := Cputime(); "Step III. 3-term relations."; end if; Tgens, Tquot := ManSym3termQuotient(mlist, eps, Sgens, Squot, Scoef); if HeckeSuperverbose then Cputime(t), " seconds."; end if; if #Tgens lt 1 then return TrivialSpace(k,eps,sign); end if; V := VectorSpace(F,#Tgens); // The reason we use "WithBasis" is that for some reason Magma // distinguishes different elements of the category ModTupFld // created via this constructor (it still says they are equal but // their attributes are kept distinct). If we did not use WithBasis, // then any two space of modular symbols of the same dimension // would be equal. Someday this will be replaced with // MY OWN CATEGORY. It is best to use this hack for now, // so in the future, when I get a category, minimal changes // will be necessary. M := VectorSpaceWithBasis(Basis(V)); SetCategory(M,ModSym); M`k := k; M`N := N; M`eps := eps; M`sign := sign; M`F := F; M`mlist:= mlist; M`quot := rec; // M`quot`Tquot lies in V and really M is almost a V, so there // is a possibility of a "circular reference bug" in Magma // leading to a memory leak! M`quot`Tquot := [V!v : v in M`quot`Tquot]; // move them into M. if HeckeVerbose then ", ",Cputime(tt), "seconds."; end if; return M; end intrinsic; function TrivialSpace(k,eps,sign) F := FieldOfValues(eps); V := VectorSpace(F,0); N := CharacterLevel(eps); M := VectorSpaceWithBasis(Basis(V)); M`k := k; M`N := N; M`eps := eps; M`sign := sign; M`F := F; M`quot := rec; SetCategory(M,ModSym); return M; end function; intrinsic Weight(M::ModTupFld) -> RngIntElt {Return the weight of the space M of modular symbols.} if not assigned M`k then return Weight(ModSymParent(M)); end if; return M`k; end intrinsic; intrinsic Level(M::ModTupFld) -> RngIntElt {Return the level of the space M of modular symbols.} return M`N; end intrinsic; intrinsic Conductor(M::ModTupFld) -> RngIntElt {Return the level of the space M of modular symbols.} return M`N; end intrinsic; intrinsic BaseField(M::ModTupFld) -> . {Return the base field of the space M of modular symbols.} if not assigned M`F then return BaseField(ModSymParent(M)); end if; return M`F; end intrinsic; intrinsic Character(M::ModTupFld) -> SeqEnum {Return the Dirichlet character of the space M of modular symbols.} case GetCategory(M): when ModSym: return M`eps; when ModSymFac: return M`M`eps; else: error "Character -- invalid input type."; end case; end intrinsic; intrinsic IsPlusQuotient(M::ModTupFld) -> SeqEnum {Return true iff working in the plus one quotient of the space M of modular symbols.} if not assigned M`sign then return M`M`sign eq 1; end if; return M`sign eq 1; end intrinsic; intrinsic IsMinusQuotient(M::ModTupFld) -> SeqEnum {Return true iff working in the plus one quotient of the space M of modular symbols.} if not assigned M`sign then return M`M`sign eq -1; end if; return M`sign eq -1; end intrinsic; intrinsic Sign(M::ModTupFld) -> RngIntElt {Returns the sign of M.} if not assigned M`sign then return M`M`sign; end if; return M`sign; end intrinsic; intrinsic IsogenyClass(A::ModTupFld) -> RngIntElt {Returns the isogeny class of A.} if not assigned A`iso then return 0; end if; return A`iso; end intrinsic; /********************************************************* * * * COMPUTE HECKE OPERATORS * * * *********************************************************/ intrinsic HeilbronnMerel(n::RngIntElt) -> SeqEnum {Return Heilbronn matrices of determinant n, as given by Merel. Here n can be composite.} H := []; i := 0; for a in [1..n] do for d in [1..n] do; // ad-bc=n so c=0 and ad=n, or b=(ad-n)/c bc := a*d - n; if bc lt 0 then continue; end if; if bc eq 0 then for b in [0..a-1] do Append(~H,[a,b,0,d]); end for; end if; if bc ne 0 then for c in Divisors(Abs(bc)) do if c lt d and Floor(bc/c) lt a then Append(~H,[a,Floor(bc/c),c,d]); end if; end for; elif 0 lt a then for c in [1..d-1] do Append(~H,[a,0,c,d]); end for; end if; end for; end for; return H; end intrinsic; intrinsic HeilbronnCremona(p::RngIntElt) -> SeqEnum {Return the Heilbronn matrices of determinant p, as defined by Cremona.} assert IsPrime(p); if p eq 2 then return [[1,0, 0,2], [2,0, 0,1], [2,1, 0,1], [1,0, 1,2]]; end if; ans := [[1,0, 0,p]]; for r in [Ceiling(-p/2)..Floor(p/2)] do x1:=p; x2:=-r; y1:=0; y2:=1; a:=-p; b:=r; c:=0; x3:=0; y3:=0; q:=0; Append(~ans, [x1,x2, y1,y2]); while b ne 0 do q := Round(a/b); c := a-b*q; a := -b; b := c; x3 := q*x2-x1; x1 := x2; x2 := x3; y3 := q*y2-y1; y1 := y2; y2 := y3; Append(~ans, [x1,x2, y1, y2]); end while; end for; return ans; end intrinsic; intrinsic Heilbronn(n::RngIntElt) -> SeqEnum {If n is prime, return Cremona's Heilbronn matrices; otherwise return Merel's.} requirege n,1; if IsPrime(n) then return HeilbronnCremona(n); else return HeilbronnMerel(n); end if; end intrinsic; // m is a sequence representing sum alp*e_i, // where e_i runs through the initial list of // generating Manin symbols. intrinsic InitialManinSymbolGenListToSquot(m::SeqEnum, M::ModTupFld, Tmat::ModMatFldElt) -> ModTupFldElt {} Tquot := M`quot`Tquot; V:=Parent(Tquot[1]); Sgens := M`quot`Sgens; Scoef := M`quot`Scoef; Squot := M`quot`Squot; // Tmat is a map from W to V, where W is W:=VectorSpace(Field(V),#Sgens); v := W!0; for t in m do v[Squot[t[2]]] +:= t[1]*Scoef[t[2]]; end for; return v; end intrinsic; intrinsic InitialManinSymbolGenListToM(m::SeqEnum, M::ModTupFld) -> ModTupFldElt {} Scoef := M`quot`Scoef; Tquot := M`quot`Tquot; V:=Parent(Tquot[1]); Squot := M`quot`Squot; return &+[V|t[1]*Scoef[t[2]] *Tquot[Squot[t[2]]] : t in m]; end intrinsic; /////////////////////////////////////////////////////////////// // COMPUTATION of HECKE OPERATORS // /////////////////////////////////////////////////////////////// intrinsic Tn(M::ModTupFld, n::RngIntElt) -> AlgMatElt {} return HeckeOperator(M,n); end intrinsic; intrinsic HeckeOperator(M::ModTupFld, n::RngIntElt) -> AlgMatElt {} requirege n,1; return HeckeOperator(M,Heilbronn(n)); end intrinsic; intrinsic HeckeOperator(M::ModTupFld, Heil::SeqEnum) -> AlgMatElt {Compute the n-th Hecke operator Tn on the space M of modular symbols.} if Dimension(M) eq 0 then return Hom(M,M)!0; end if; F := M`F; mlist := M`mlist; eps := M`eps; quot := M`quot; Sgens := quot`Sgens; Squot := quot`Squot; Tgens := quot`Tgens; Tquot := quot`Tquot; Mn := MatrixAlgebra(F,Dimension(M)); V := Parent(Tquot[1]); T := []; W := VectorSpace(F, #Sgens); Tmat := RMatrixSpace(F, Dimension(W), Dimension(V)) ! [Tquot[i] : i in [1..#Sgens]]; for j in [1..#Tgens] do i := Sgens[Tgens[j]]; v := InitialManinSymbolGenListToSquot( &cat[ ManinSymbolApply(h,i, mlist,eps) : h in Heil], M,Tmat); Append(~T,v); end for; S := RMatrixSpace(F, #T, #Sgens) ! T; return Mn!(S*Tmat); end intrinsic; intrinsic TnSparse(M::ModTupFld, n::RngIntElt, sparsevec::SeqEnum) -> AlgMatElt {} requirege n,1; if HeckeSuperverbose then "TnSparse ", n; end if; return TnSparse(M, Heilbronn(n), sparsevec); end intrinsic; intrinsic TnSparse(M::ModTupFld, Heil::SeqEnum, sparsevec::SeqEnum) -> AlgMatElt {Compute the action of the Hecke operator defined by the Heilbronn matrices Heil on the sparse vector v.} if Dimension(M) eq 0 then return M!0; end if; F := M`F; k := M`k; mlist := M`mlist; eps := M`eps; quot := M`quot; Sgens := quot`Sgens; Tgens := quot`Tgens; Tquot := quot`Tquot; Mn := MatrixAlgebra(F,Dimension(M)); v := M!0; for m in sparsevec do i := Sgens[Tgens[m[2]]]; if k eq 2 then // faster, since always just one entry. x := [ ManinSymbolApply(h,i, mlist,eps)[1] : h in Heil]; else x := &cat[ ManinSymbolApply(h,i, mlist,eps) : h in Heil]; end if; v +:= m[1]*InitialManinSymbolGenListToM(x,M); end for; return v; end intrinsic; /******************************************************* CyclicHeckeImage: The goal of this function is to compute the cyclic Hecke module generated by a vector v. The result is returned as a subspace with Z-basis. Let Tn be the matrix representing the right action of the nth Hecke operator on P. We must compute gens for an F-module who (Z-span) equals span {Tn*v : 1<=n<=b} where b is a bound on the number of Tn necessary to generate the Hecke algebra as a Z-module. We do this by: 1) Computing Mp for p prime using Hecke action on easy elements of M_k(Q). 2) Iteratively computing prod Mp^i*c where n=prod p^i varies from 1 to b. *********************************************************/ intrinsic SparseRepresentation(v::ModTupFldElt) -> SeqEnum {Return sparse representation of vector v.} sp := []; v := Eltseq(v); for i in [1..#v] do if v[i] ne 0 then Append(~sp, ); end if; end for; return sp; end intrinsic; intrinsic HeckeBound(M::ModTupFld) -> RngIntElt {Returns a bound b such that T1, ..., Tb generate the Hecke algebra as a Z-module.} b := HeckeBound(M`k, M`N); if not IsTrivial(M`eps) then " Warning: William Stein has not proved that his bound on the number"; " of Hecke operators needed to generate the Hecke algebra is correct"; " when Character(M) is not trivial."; b *:= 3; end if; return b; end intrinsic; forward idxG0; intrinsic HeckeBound(k::RngIntElt, N::RngIntElt) -> RngIntElt {Returns a bound b such that T1, ..., Tb generate the Hecke algebra as a Z-module.} return Ceiling(k * idxG0(N) / 12); end intrinsic; intrinsic VectorSpaceZBasis(B::SeqEnum) -> ModTupFld {Return vector space with basis the lattice defined by the elements of B.} V := VectorSpace(Parent(B[1][1]),#Eltseq(B[1])); Mat := RMatrixSpace(Rationals(),#B, Dimension(V)); X := Mat!B; if X eq Mat!0 then return VectorSpaceWithBasis([V|]); end if; Latbase := Basis(Lattice(X)); return VectorSpaceWithBasis([V!v : v in Latbase]); end intrinsic; intrinsic HeckeSpan(M::ModTupFld, v::ModTupFldElt) -> ModTupFld {This function computes the Hecke module generated by a vector v. The result is returned as a subspace with Z-basis.} s := SparseRepresentation(v); b := HeckeBound(M); B := [v] cat [TnSparse(M, n, s) : n in [2..b]]; return VectorSpaceZBasis(B); end intrinsic; intrinsic HeckeSpan(M::ModTupFld, V::ModTupFld) -> ModTupFld {This function computes the Hecke module generated by a basis of V. The result is returned as a subspace with Z-basis.} S := [SparseRepresentation(v) : v in Basis(V)]; b := HeckeBound(M); H := [Heilbronn(n) : n in [1..b]]; B := &cat[[TnSparse(M, H[n], s) : n in [1..b]] : s in S]; return VectorSpaceZBasis(B); end intrinsic; /****************************************************** * Computation of T_p on dual space. * ******************************************************/ CFastData := recformat< V, e, VEinv>; intrinsic FastTnData(M::ModTupFld, V::SeqEnum) -> Rec {Compute data needed by FastTn.} // Step 1: find a set of very simple vectors e[1],...,e[n] // in M such that det(v[i]*e[j]) =/= 0. // 1. [Compute e] e = set of positions so that elements of the space // spanned by the columns of V are determined by the entries in these // spots. This is accomplished by row reducing, and setting e equal to // the sequence of column numbers for which the columns are pivotal. n := #V; B := Basis(sub); assert #B eq n; // Find pivot columns. e := Pivots(B); // Step 2: Compute the matrix V*E of inner products. VE := RMatrixSpace(M`F,n,n)![V[i][e[j]] : i in [1..n], j in [1..n]]; VEinv := VE^(-1); return rec; end intrinsic; intrinsic FastTnData(M::ModTupFld, V::ModTupFld) -> Rec {} return FastTnData(M, Basis(V)); end intrinsic; intrinsic FastTn(M::ModTupFld, FastData::Rec, n::RngIntElt) -> AlgMatElt {} return FastTn(M,FastData,n,Heilbronn(n)); end intrinsic; intrinsic FastTn(M::ModTupFld, FastData::Rec, n::RngIntElt, H::SeqEnum) -> AlgMatElt {Compute action of Transpose(Tn) on the Hecke-stable subspace V.} F := M`F; V := FastData`V; n := #V; m := Dimension(M); e := FastData`e; VEinv := FastData`VEinv; // The next step is where all of the time is spent. TE := [TnSparse(M, H, [<1,e[i]>]) : i in [1..n]]; Vmat := RMatrixSpace(F, n, m) ! V; TEmat := RMatrixSpace(F, n, m) ! TE; return MatrixAlgebra(F,n)!(Vmat*Transpose(TEmat)*VEinv); end intrinsic; intrinsic FastTn(FastData::Rec, n::RngIntElt, M::ModTupFld) -> AlgMatElt {} return FastTn(M,FastData,n,Heilbronn(n)); end intrinsic; intrinsic FastTn(M::ModTupFld, V::ModTupFld, n::RngIntElt) -> AlgMatElt, Rec {Compute action of Transpose(Tn) on the Hecke-stable subspace V.} FastData := FastTnData(M, V); return FastTn(M, FastData, n); end intrinsic; intrinsic FastTn(M::ModTupFld, V::ModTupFld, n::RngIntElt, H::SeqEnum) -> AlgMatElt, Rec {Compute action of Transpose(Tn) on the Hecke-stable subspace V.} FastData := FastTnData(V,M); return FastTn(M, FastData, n, H); end intrinsic; intrinsic FastTn( M::ModTupFld, V::ModTupFld, nlist::SeqEnum, H::SeqEnum) -> SeqEnum, Rec {Compute action of Transpose(Tn), n in plist, on the Hecke-stable subspace V.} FastData := FastTnData(M, V); return [FastTn(M, FastData, n, H[n]): n in nlist]; end intrinsic; intrinsic FastTn(M::ModTupFld, V::ModTupFld, nlist::SeqEnum) -> SeqEnum, Rec {Compute action of Transpose(Tn), n in nlist, on the Hecke-stable subspace V.} FastData := FastTnData(M,V); return [FastTn(M, FastData, n): n in nlist]; end intrinsic; /****************************************************** Action of T_n on modular symbols basis. ******************************************************/ function ActionOnModularSymbolsBasis(g, M) // 1. Compute basis of modular symbols for M. B := ModularSymbolsBasis(M); // 2. Apply g to each basis element. gB := [ModularSymbolApply(M, g,B[i]) : i in [1..#B]]; // 3. Map the result back to M. gM := [ConvFromModularSymbol(M,gB[i]) : i in [1..#gB]]; return MatrixAlgebra(M`F,Dimension(M))!gM; end function; // Compute the action of the Atkin-Lehner involution on M, // when this makes sense. intrinsic Wq(M::ModTupFld, q::RngIntElt) -> AlgMatElt {Compute the action of the Atkin-Lehner involution Wq on M, when this makes sense (i.e., trivial character, even weight). The Atkin-Lehner map Wq is rescaled so that it is an involution, except when k>2 and char(F) divides q.} if GetCategory(M) eq ModSymFac then return Restrict(Transpose(Wq(ModSymParent(M),q)),M); end if; N := M`N; k := M`k; assert q ne 0; assert N mod q eq 0; assert k mod 2 eq 0; assert IsTrivial(M`eps); repeat d := Gcd(Integers()!(N/q),q); if d eq 1 then break; end if; q *:= d; until false; assert (N mod q) eq 0; // 1. Compute matrix to act by. d, x, y:= XGCD(q,-Integers()!(N/q)); g := [q*x, y, N, q]; A := ActionOnModularSymbolsBasis(g, M); p := Characteristic(M`F); if p eq 0 or q mod p ne 0 then A /:= q^(Integers()!(k/2)-1); end if; return A; end intrinsic; intrinsic StarInvolution(M::ModTupFld) -> AlgMatElt {Compute the conjugation involution * on V. This involution is defined by the 2x2 matrix [-1,0,0,1]; it sends X^i*Y^j\{u,v\} to (-1)^j*X^i*Y^j \{-u,-v\}.} return ActionOnModularSymbolsBasis([-1,0,0,1], M); end intrinsic; intrinsic Tnpoly(n::RngIntElt, M::ModTupFld) -> SeqEnum {Compute the characteristic polynomial of the p-th Hecke operator Tn on the space M of modular symbols. Uses the modular algorithm, without proof.} return ModularCharpoly(HeckeOperator(n,M)); end intrinsic; intrinsic Restrict (A::AlgMatElt, L::Lat) -> AlgMatElt { Suppose A is a linear transformation of V which leaves the lattice L invariant. Returns A restricted to W, with respect to the basis Basis(W).} Z := Basis(L); V := VectorSpace(Parent(A[1,1]), Dimension(L)); B := [V!Z[i] : i in [1..Dimension(V)]]; return Restrict (A, B); end intrinsic; intrinsic Restrict (A::AlgMatElt, W::ModTupFld) -> AlgMatElt { Suppose A is a linear transformation of V which leaves the subspace W invariant. Returns A restricted to W, with respect to the basis Basis(W). } B := Basis(W); return Restrict (A, B); end intrinsic; intrinsic Restrict (A::AlgMatElt, W::ModTupRng) -> AlgMatElt { Suppose A leaves W invariant. Returns A restricted to W, with respect to the basis Basis(W).} B := Basis(W); return Restrict (A, B); end intrinsic; intrinsic Restrict (A::AlgMatElt, B::SeqEnum) -> AlgMatElt { Suppose A is a linear transformation of V which leaves the subspace W spanned by the vectors listed in B invariant. Returns A restricted to W, with respect to the basis B.} S := RSpaceWithBasis(B); return MatrixAlgebra(Parent(A[1,1]),#B) ! &cat[Coordinates(S, S.i*A) : i in [1..#B]]; end intrinsic; intrinsic KernelOn (A::AlgMatElt, B::SeqEnum) -> ModTupFld {Suppose B is a basis for an n-dimensional subspace of some ambient space and A is an nxn matrix. Then A defines a linear transformation of the space spanned by B. This function returns the kernel of that transformation.} n := #B; if n eq 0 then return B; end if; assert Nrows(A) eq Ncols(A) and Nrows(A) eq n; K := Kernel(A); return VectorSpaceWithBasis(LinearCombinations(Basis(K),B)); end intrinsic; intrinsic KernelOn (A::AlgMatElt, W::ModTupFld) -> ModTupFld {} return KernelOn (A, Basis(W)); end intrinsic; intrinsic HeckeOperatorOn(M::ModTupFld, V::ModTupFld, n::RngIntElt) -> AlgMatElt {Compute the action of the Hecke operator Tn with respect to the basis for the subspace V of M.} return Restrict(HeckeOperator(M,n), V); end intrinsic; intrinsic HeckeOperatorMkZ(M::ModTupFld, n::RngIntElt) -> AlgMatElt {Compute the action of Tn on the integral modular symbols.} return Restrict(HeckeOperator(M,n), IntegralModularSymbols(M)); end intrinsic; intrinsic HeckeOperatorSkZ(M::ModTupFld, n::RngIntElt) -> AlgMatElt {Compute the action of T_p on a subspace of integral modular symbols.} return Restrict(HeckeOperator(M,n), SkZ(M)); end intrinsic; /********************************************************** * * * BOUNDARY SYMBOLS: Compute the boundary map delta. * * * **********************************************************/ /////////////////////////////////////////////////////////////////////////// //////////// IMPLEMENTATION 2: CUSPS FOR NONTRIVIAL CHARACTER ///////////// /////////////////////////////////////////////////////////////////////////// function ReduceCusp(a) d := Gcd(a); return [Integers()|x/d : x in a]; end function; function CuspEquiv(N,a,b) u1, v1 := Explode(ReduceCusp(a)); u2, v2 := Explode(ReduceCusp(b)); s1 := 0; s2 := 0; if [u1,v1] ne [0,1] then s1 := (v1 eq 0 or v1 eq 1) select 1 else InverseMod(u1,Abs(v1)); end if; if [u2,v2] ne [0,1] then s2 := (v2 eq 0 or v2 eq 1) select 1 else InverseMod(u2,Abs(v2)); end if; g := Gcd(v1*v2,N); a := s1*v2 - s2*v1; if a mod g ne 0 then return false, 1; end if; // Now we know the cusps are equivalent. Use the proof of Prop2.2.3 // of Cremona to find a matrix in Gamma_0(N) relating them. dum,s2,r2 := Xgcd(u2,-v2); assert dum eq 1; dum,s1,r1 := Xgcd(u1,-v1); assert dum eq 1; a := s1*v2 - s2*v1; assert a mod g eq 0; // solve x*v1*v2 + a = 0 (mod N). d,x0,y0 := Xgcd(v1*v2,N); // x0*v1*v2 + y0*N = d = g. // so x0*v1*v2 - g = 0 (mod N) x := -x0 * Integers()!(a/g); // now x*v1*v2 + a = 0 (mod N) s1p := s1+x*v1; return true, (u2*s1p-r2*v1) mod N; end function; function CuspToFreeHelper(M, sign, a) if not assigned M`cusplist then M`cusplist := []; end if; list := M`cusplist; F := M`F; eps := M`eps; N := M`N; k := M`k; a := ReduceCusp(a); if a[2] lt 0 then a[1] *:= F!-1; a[2] *:= F!-1; end if; // Check if we've already encountered this cusp. for i in [1..#list] do b := list[i]; equiv, alp := CuspEquiv(N, b[1], a); // [gam_alp(b[1])]=?=[a]. if equiv then if b[2] eq 0 then return ; end if; if IsTrivial(Character(M)) then return <1,i>; else return ; end if; end if; if sign ne 0 then equiv, alp := CuspEquiv(N, b[1], [-a[1],a[2]]); if equiv then if b[2] eq 0 then return ; end if; if IsTrivial(Character(M)) then return ; else return ; end if; end if; end if; end for; // Determine if this cusp class is killed by the relations. c := F!1; if not IsTrivial(Character(M)) then u := a[1]; v := a[2]; g := Gcd(N,v); x := Integers()!(N/Gcd(N,v)); for j in [0..g-1] do alp := 1-j*x; if Gcd(alp,N) eq 1 then if (v*(1-alp)) mod N eq 0 and (u*(1-alp)) mod g eq 0 then if Evaluate(eps,alp) ne 1 then c := F!0; break; end if; end if; end if; end for; end if; Append(~list,); M`cusplist := list; i := #list; if sign ne 0 then // Check is sign relation kills this cusp. ans := CuspToFreeHelper(M,0,[-a[1],a[2]]); if ans[2] eq i and ans[1] ne sign then M`cusplist[i][2] := 0; c := F!0; end if; end if; return ; end function; intrinsic CuspToFree(M::ModTupFld, a::SeqEnum) -> Tup {Given a cusp a=[u,v] this function finds the index i and a scalar alpha such that a = alpha*(i-th standard cusp). It then returns alpha, i.} return CuspToFreeHelper(M,Sign(M),a); end intrinsic; forward LiftToCosetRep; intrinsic BoundaryMap(M::ModTupFld) -> ModMatFldElt {Compute the Boundary map.} if not assigned M`BoundaryMap then if Dimension(M) eq 0 then M`BoundaryMap := RMatrixSpace(M`F,0,0)!0; return M`BoundaryMap; end if; Tgens := M`quot`Tgens; Sgens := M`quot`Sgens; F := M`F; n := #Tgens; D := [ [] : i in [1..n]]; for j in [1..n] do i := Sgens[Tgens[j]]; uv, w, ind := UnwindManinSymbol(i,M`mlist); g := LiftToCosetRep(uv, M`N); if w eq M`k-2 then Append(~D[j], CuspToFree(M,[g[1],g[3]])); end if; if w eq 0 then // substract the second term. t := CuspToFree(M,[g[2],g[4]]); t[1] *:= -1; Append(~D[j], t); end if; end for; if &cat D eq [] then M`BoundaryMap := RMatrixSpace(M`F,Dimension(M),0)!0; return M`BoundaryMap; end if; m := Max([x[2] : x in &cat D]); V := VectorSpace(F,m); Drows := [&+[V| x[1]*V.x[2] : x in D[i]] : i in [1..n]]; M`BoundaryMap := RMatrixSpace(M`F,n,m) ! Drows; end if; return M`BoundaryMap; end intrinsic; intrinsic CuspidalSymbols(M::ModTupFld) -> ModTupFld {Return the subspace of S_k(N,Q) of cuspidal modular symbols.} if not assigned M`SkF then M`SkF := Kernel(BoundaryMap(M)); end if; return M`SkF; end intrinsic; intrinsic CuspidalSymbols(M::ModTupFld) -> ModTupFld {Return the subspace of S_k(N,Q) of cuspidal modular symbols.} if not assigned M`SkF then D := BoundaryMap(M); if HeckeSuperverbose then "Computing kernel of boundary map: "; t := Cputime(); end if; M`SkF := Kernel(D); if HeckeSuperverbose then Cputime(t), " seconds."; end if; end if; return M`SkF; end intrinsic; intrinsic Sk(M::ModTupFld) -> ModTupFld {} return CuspidalSymbols(M); end intrinsic; intrinsic IntegralCuspidalSymbols(M::ModTupFld) -> ModTupFld {} if not assigned M`SkZ then Z := IntegralModularSymbols(M); D := BoundaryMap(M); DofZ := [v*D : v in Basis(Z)]; DZ := RMatrixSpace(M`F,Nrows(D),Ncols(D)) ! DofZ; M`SkZ := VectorSpaceWithBasis( LinearCombinations(Basis(IntegerKernel(DZ)), Basis(Z)) ); end if; return M`SkZ; end intrinsic; intrinsic SkZ(M::ModTupFld) -> ModTupFld {} case GetCategory(M): when ModSym: return IntegralCuspidalSymbols(M); when ModSymFac: return DecompZ(M); else error "SkZ only takes a ModSym or ModSymFac as argument."; end case; end intrinsic; intrinsic SkZPlus(M::ModTupFld) -> ModTupFld {} if not assigned M`SkZplus then S := SkZ(M); star := Restrict(StarInvolution(M),S); M`SkZplus := IntegerKernelOn(star-1, S); end if; return M`SkZplus; end intrinsic; intrinsic SkZMinus(M::ModTupFld) -> ModTupFld {} if not assigned M`SkZminus then S := SkZ(M); star := Restrict(StarInvolution(M),S); M`SkZminus := IntegerKernelOn(star+1, S); end if; return M`SkZminus; end intrinsic; /* function EisensteinSymbols(M) end function; */ /******************************************************************* DEGENERACY COSET REPRESENTATIVES: Let N be a positive integer and M a divisor of N. Let t be a divisor of N/M, and let T be the 2x2 matrix T=[0,0, 0,t]. Find coset reps for Gamma_0(N) \ T Gamma_0(M). FACTS: T^(-1)*(a,b,c,d)*T = (a,bt,c/t,d) T^(-1)*Gamma_0(N)*T is contained in Gamma_0(M) Gamma_0(N)*T is contained in T*Gamma_0(M) SOLUTION: (1) Computes representatives for Gamma_0(N/t,t) inside of Gamma_0(M): COSET EQUIVALENCE: Two right cosets represented by (a,b;c,d) & (a',b';c',d') of Gamma_0(N/t,t) in SL_2(Z) are equivalent iff (a,b)=(a',b') as points of P^1(t), i.e., ab'=ba'(mod t), and (c,d)=(c',d') as points of P^1(N/t). ALGORITHM: (A) Compute the number of cosets. (B) Compute a random element x of Gamma_0(M). (C) Check if x is equivalent to anything generated so far; if not, add x to the list. (D) Continue until the list is as long as the bound computed in A. (2) There is a natural bijection Gamma_0(N)\T*Gamma_0(M) <---> Gamma_0(N/t,t)\Gamma_0(M) given by Tr <---------> r Consequently we obtain coset representatives for Gamma_0(N)\T*Gamma_0(M) by left multiplying by T each coset representative of Gamma_0(N/t,t)\Gamma_0(M) found in step 1. ********************************************************************/ intrinsic DegeneracyCosetReps(M, N, d) -> . {} n := idxG0(N)/idxG0(M); // number of coset representatives. Ndivd := N div d; R := []; // coset reps found so far. max := 4*(n+10); halfmax := Round(max/2); while #R lt n do // try to find another coset representative. cc := M*Random(-halfmax, halfmax); dd := Random(-halfmax, halfmax); g, bb, aa := XGCD(-cc,dd); if g eq 0 then continue; end if; cc div:= g; dd div:= g; if cc mod M ne 0 then continue; end if; // Test if we've found a new coset representative. isnew:=true; for r in R do if ((r[2]*aa - r[1]*bb) mod d eq 0) and ((r[4]*cc - r[3]*dd) mod Ndivd eq 0) then isnew := false; break ; end if; end for; // If it is new add it to the list. if isnew then Append(~R,[aa,bb,cc,dd]); end if; end while; // Return the list left multiplied by T. return [[r[1],r[2],r[3]*d,r[4]*d] : r in R]; end intrinsic; /********************************************************** DEGENERACY AND INCLUSION MAPS -- newforms. If possible compute the map M1 --> M2 corresponding to the integer d. WARNING: There might be a problem here when the characteristic of the base field divides d, because d occurs in the denominators of the DegeneracyReps. d:=1;c:=3;M:=11;N:=M*d*c; A:=Mk(M,2); B:=Mk(N,2); bet:=dm(A,B,d); alp:=dm(B,A,d); bet*alp; // must be a scalar. **********************************************************/ intrinsic DegeneracyMap(M1::ModTupFld, M2::ModTupFld) -> AlgMatElt {If possible compute the map M1 --> M2 corresponding to the integer d.} return DegeneracyMap(M1,M2,1); end intrinsic; intrinsic DegeneracyMap(M1::ModTupFld, M2::ModTupFld, d::RngIntElt) -> AlgMatElt {If possible compute the degeneracy map M1 --> M2 corresponding to the integer d. Here M1 and M2 are spaces of modular symbols, of the same weight and d is a divisor of the quotient of the two levels N1 and N2. If the N1 divides N2 then the "raising" map is computed, if N1 is divisible by N2, then the "lowering" map is computed.} assert M1`F eq M2`F; assert Weight(M1) eq Weight(M2); N1 := Level(M1); N2 := Level(M2); if N1 eq N2 then assert d eq 1; assert M1 eq M2; // we could write down the map between // different presentations of M_k(N), but // we won't for now, since there should // never be a need. F := BaseField(M1); return Hom(M1,M1)!MatrixAlgebra(F,Degree(M1))!1; end if; if N1 mod N2 eq 0 then // N2 divides N1 -- lower assert (N1 div N2) mod d eq 0; B := ModularSymbolsBasis(M1); D := [d,0,0,1]; DB := [ModularSymbolApply(M1,D, B[i]) : i in [1..#B] ]; A := [ConvFromModularSymbol(M2,DB[i]) : i in [1..#DB]]; return Hom(M1, M2) ! A; elif N2 mod N1 eq 0 then// N1 divides N2 -- raise level assert (N2 div N1) mod d eq 0; B := ModularSymbolsBasis(M1); R := DegeneracyCosetReps(N1, N2, d); eps := Character(M1); if IsTrivial(eps) then RB := [ &cat[ModularSymbolApply(M1, r, B[i]) : r in R] : i in [1..#B]]; A := [ConvFromModularSymbol(M2,RB[i]) : i in [1..#RB]]; else A := [ &+[Evaluate(eps,r)^(-1)*ConvFromModularSymbol(M2, ModularSymbolApply(M1, r, B[i])) : r in R] : i in [1..#B]]; end if; return Hom(M1, M2) ! A; else error "DegeneracyMap: N1 must be divisible by N2, or vice-versa."; end if; end intrinsic; intrinsic OldModularSymbols(M::ModTupFld, NN::RngIntElt) -> ModTupFld {Returns the space Mk(NN,eps') associated to Mk(N,eps). Here NN must be a divisor of N.} N := Level(M); if NN eq N then return M; end if; divs := Divisors(N); if not assigned M`old then M`old := SequenceToList(divs); end if; pos := Index(divs,NN); if pos eq 0 then error "OldModularSymbols: NN =",NN,"is not a divisor of the level =",N; end if; if Type(M`old[pos]) eq RngIntElt then /**************************************************************** This object does not know the old modular symbols at level NN. Here's what we do: (1) If this object has a parent, ask it for the modular symbols at level NN. If so, done; otherwise continue with the next step. (2) If not, construct the space of symbols at level NN. (3) If this object has a parent, tell it about the symbols at level NN. ****************************************************************/ // Step 1. if not ModSymIsRoot(M) then // Parent takes care of everything recursively. M`old[pos] := OldModularSymbols(ModSymParent(M),NN); else // Step 2: We've got to do the work. NN := Integers()!NN; eps := Character(M); F := BaseField(M); if not CanReduceModM(eps,NN) then M`old[pos] := BaseExtend(TrivialSpace(Weight(M),eps,Sign(M)),F); else epsNN := AssociatedModMCharacter(eps,NN); M`old[pos] := BaseExtend( ModularSymbols(epsNN,Weight(M),Sign(M)),F); end if; (M`old[pos])`M := M; // tag this as *our* child. end if; end if; return M`old[pos]; end intrinsic; /********************************************************** For each prime p dividing the level there are two maps alp_1,alp_p: Mk(N,e) ----> Mk(N/p,e') The intersection Kernel(alp_1) meet Kernel(alp_p) is called the p-new subspace Mkpnew(N) of Mk(N). If e doesn't factor through (Z/(N/p)Z)^* then e' is by definition 0. ********************************************************/ intrinsic pNewSubspace(M::ModTupFld, p::RngIntElt) -> ModTupFld {Returns the p-new subspace of M.} assert IsPrime(p); N := Level(M); k := Weight(M); assert N mod p eq 0; eps := Character(M); NN := Integers()!(N/p); if not CanReduceModM(eps,NN) then return M; end if; old := OldModularSymbols(M,NN); if Dimension(old) eq 0 then return M; end if; h1 := DegeneracyMap(M, old, 1); hp := DegeneracyMap(M, old, p); return Kernel(h1) meet Kernel(hp); end intrinsic; intrinsic pNewDualSubspace(M::ModTupFld, p::RngIntElt) -> ModTupFld {Returns the p-new subspace of M.} assert IsPrime(p); N := Level(M); k := Weight(M); assert N mod p eq 0; eps := Character(M); NN := Integers()!(N/p); if not CanReduceModM(eps,NN) then return M; end if; old := OldModularSymbols(M,NN); h1 := DegeneracyMap(old, M, 1); hp := DegeneracyMap(old, M, p); return Kernel(Transpose(h1)) meet Kernel(Transpose(hp)); end intrinsic; intrinsic Sknew(M::ModTupFld) -> ModTupFld {} return CuspidalNewSubspace(M); end intrinsic; intrinsic CuspidalNewSubspace(M::ModTupFld) -> ModTupFld {} if not assigned M`sknew then M`sknew := NewSubspace(M) meet CuspidalSymbols(M); end if; return M`sknew; end intrinsic; intrinsic NewSubspace(M::ModTupFld) -> ModTupFld {Returns the p-new subspace of M.} if not assigned M`mknew then if Level(M) eq 1 then M`mknew := M; else M`mknew := &meet[pNewSubspace(M,p[1]) : p in Factorization(Level(M))]; end if; end if; return M`mknew; end intrinsic; intrinsic Mknewdual(M::ModTupFld) -> ModTupFld {} return NewDualSubspace(M); end intrinsic; intrinsic Sknewdual(M::ModTupFld) -> ModTupFld {} return CuspidalNewDualSubspace(M); end intrinsic; intrinsic CuspidalNewDualSubspace(M::ModTupFld) -> ModTupFld {} if not assigned M`sknewdual then M`sknewdual := NewDualSubspace(M); // cut out the cuspidal part using Hecke operators. S := Sknew(M); if Dimension(S) eq 0 then M`sknewdual := S; else p := 2; while Dimension(M`sknewdual) gt Dimension(S) do T := HeckeOperator(M,p); f := ModularCharpoly(Restrict(T,S)); Ton := Restrict(Evaluate(f,Transpose(T)),M`sknewdual); M`sknewdual := KernelOn(Ton,M`sknewdual); p := NextPrime(p); end while; end if; end if; return M`sknewdual; end intrinsic; intrinsic NewDualSubspace(M::ModTupFld) -> ModTupFld {Returns the p-new subspace of M dual.} if not assigned M`mknewdual then M`mknewdual := &meet ([pNewDualSubspace(M,p[1]) : p in Factorization(Level(M))] cat [M]); end if; return M`mknewdual; end intrinsic; /***************************************************************** THE THETA OPERATOR: On q-expansions, the theta operator is q*d/dq. It takes the q-expansion of a mod ell modular form of weight k to a mod ell modular form of weight k + ell + 1. On modular symbols the theta operator is simply multiplication by X^{ell}*Y - X*Y^{ell}. Proof: Let Q(X,Y) = X^{ell}*Y - X*Y^{ell}. Then for any g in GL_2(Q), Q(det(g) g^{-1}(X,Y)) = det(g) Q(X,Y) (mod ell). (This congruence, of course, does not hold in characteristic 0.) To prove the congruence simply multiply everything out and use the fact that elements in F_{ell} satisfy a^{ell} = a. This result was suggested by Kevin Buzzard; see the proof of Proposition~4.5.1, page 61, of his Ph.D. thesis. *****************************************************************/ intrinsic ThetaOperator(M1::ModTupFld, M2::ModTupFld, x::ModTupFldElt) -> ModTupFldElt {Computes the image in M2 of x in M1 under the theta operator q*d/dq : M1 --> M2. Thus M1 and M2 should be spaces of modular symbols mod ell of the same level, and Weight(M2) = Weight(M1)+ell+1.} ell := Characteristic(BaseField(M1)); if not IsPrime(ell) then error "The characteristic must be prime."; end if; if ell ne Characteristic(BaseField(M2)) then error "The characteristics must be equal."; end if; if Level(M1) ne Level(M2) then error "The levels must be equal."; end if; if Character(M1)`ims ne Character(M2)`ims then error "The characters must be equal."; end if; sym := ConvToModularSymbol(M1,x); if #sym eq 0 then return M2!0; end if; P := Parent(sym[1][1]); Q := X^ell*Y - Y^ell*X; for i in [1..#sym] do sym[i][1] *:= Q; end for; return ConvFromModularSymbol(M2,sym); end intrinsic; //////////////////// END THETA OPERATOR //////////////////////////// /********************************************************** CONVERSIONS: Modular symbols <----> Manin symbols **********************************************************/ // WINDING ELEMENT intrinsic WindingElement(M::ModTupFld) -> ModTupFldElt {Returns the winding element \{0,oo\}.} return WindingElement(M,1); end intrinsic; intrinsic WindingSubmodule(M::ModTupFld) -> ModTupFldElt {Returns the winding submodule T\{0,oo\}.} return WindingSubmodule(M,1); end intrinsic; intrinsic WindingElement(M::ModTupFld, i::RngIntElt) -> ModTupFldElt {Returns the winding element X^(i-1)*Y^(k-2-(i-1))*\{0,oo\}.} assert i ge 1 and i le M`k-1; R := PolynomialRing(M`F,2); i := i-1; return ConvFromModularSymbol(M,); end intrinsic; intrinsic WindingSubmodule(M::ModTupFld, i::RngIntElt) -> ModTupFld {Returns the submodule Te spanned by the ith winding element.} assert i ge 1 and i le M`k-1; if not assigned M`winding then M`winding := Seqlist([0 : i in [1..M`k-1]]); end if; if Type(M`winding[i]) eq RngIntElt then M`winding[i] := HeckeSpan(M,WindingElement(M,i)) ; end if; return M`winding[i]; end intrinsic; /********************************************************** ConvToModularSymbol returns i-th freely generating Manin symbol written as a sum of modular symbols sum X^i*Y^(k-2-i)*{alp,beta} **********************************************************/ /********************************************************** LiftToCosetRep x = [u,v] is an element of Z/NZ x Z/NZ such that gcd (u,v,N) = 1. This function finds a 2x2 matrix [a,b, c,d] such that c=u, d=v both modulo N, and so that ad-bc=1. **********************************************************/ function LiftToCosetRep(x, N) c:=Integers()!x[1]; d:=Integers()!x[2]; g, z1, z2 := Xgcd(c,d); // We're lucky: z1*c + z2*d = 1. if g eq 1 then return [z2, -z1, c, d]; end if ; // Have to try harder. if c eq 0 then c +:= N; end if; if d eq 0 then d +:= N; end if; m := c; // compute prime-to-d part of m. repeat g := Gcd(m,d); if g eq 1 then break; end if; m div:= g; until false; // compute prime-to-N part of m. repeat g := Gcd(m,N); if g eq 1 then break; end if; m div:= g; until false; d +:= N*m; g, z1, z2 := Xgcd(c,d); if g eq 1 then return [z2, -z1, c, d]; else error "LiftToCosetRep: ERROR!!!"; end if ; end function; intrinsic ConvToManinSymbol(M::ModTupFld, i::RngIntElt) -> Tup {} mlist := M`mlist; uv, w:= UnwindManinSymbol ( M`quot`Sgens[M`quot`Tgens[i]], mlist); R := mlist`R; return ; end intrinsic; intrinsic ConvToManinSymbol(M::ModTupFld, v::ModTupFldElt) -> SeqEnum {Returns an expression in terms of Manin symbols for the element v in the space M of modular symbols.} assert Degree(v) eq Dimension(M); return [ : i in [1..Dimension(M)] | v[i] ne 0 where x is ConvToManinSymbol(M,i) ]; end intrinsic; function ConvToModSym(M, i) mlist := M`mlist; uv, w, ind := UnwindManinSymbol ( M`quot`Sgens[M`quot`Tgens[i]], mlist); g := LiftToCosetRep(uv, M`N); k := M`k; R := mlist`R; if k eq 2 then return ; // {g(0), g(oo)} end if; h := hom R | g[4]*R.1-g[2]*R.2, -g[3]*R.1+g[1]*R.2>; P := R.1^w*R.2^(k-2-w); hP := h(P); return ; end function; // a good call is "ToModularSymbol(M.i,M)". intrinsic ConvToModularSymbol(M::ModTupFld, v::ModTupFldElt) -> SeqEnum {Returns an expression in terms of modular symbols of the element v in the space M of modular symbols. We represent a point in P^1(Q) by a pair [a,b]. Such a pair corresponds to the point a/b in P^1(Q).} assert v in M; nz := [i : i in [1..Dimension(M)] | v[i] ne 0]; return [ : i in nz | true where x is ConvToModSym(M, i)]; end intrinsic; intrinsic ModularSymbolsBasis(M::ModTupFld) -> SeqEnum {Return the basis of M in terms of modular symbols.} if not assigned M`ModSymBasis then M`ModSymBasis := [ConvToModularSymbol(M, M.i) : i in [1..Dimension(M)]]; end if; return M`ModSymBasis; end intrinsic; /********************************************************** FromManinSymbol: Given a Manin symbol [P(X,Y),[u,v]], this function computes the corresponding element of the space M of Modular symbols. **********************************************************/ intrinsic ConvFromManinSymbol (M::ModTupFld, x::Tup) -> ModTupFldElt {} return ConvFromManinSymbol(M,x[1],x[2]); end intrinsic; intrinsic ConvFromManinSymbol (M::ModTupFld, uv::SeqEnum ) -> ModTupFldElt {} R := M`mlist`R; return ConvFromManinSymbol(M, R.2^(M`k-2), uv); end intrinsic; intrinsic ConvFromManinSymbol (M::ModTupFld, P::RngMPolElt, uv::SeqEnum ) -> ModTupFldElt {Given a Manin symbol [P(X,Y),[u,v]], this function computes the corresponding element of the space M of Modular symbols.} R := M`mlist`R; P := R!P; k := M`k; n := #(M`mlist`p1list); ind,s := P1Reduce(uv, M`mlist`p1list); if not IsTrivial(M`eps) then P *:= Evaluate(M`eps, s); end if; P := R!P; ans := [ : w in [0..k-2]]; ans := [x : x in ans | x[1] ne 0]; return InitialManinSymbolGenListToM(ans, M); end intrinsic; /********************************************************** Given a modular symbol P(X,Y)*{alp,beta}, FromModularSymbol returns the corresponding element of M. The input is a pair , where P is a polynomial in X and Y, and x=[[a,b],[c,d]] is a pair of elements of P^1(Q), where [a,b] <--> a/b and [c,d] <--> c/d. **********************************************************/ // This function returns the n-th continued fraction convergent // of the continued fraction defined by the sequence of integers v. // We assume n>=1. If the continued fraction integers are // v = [a1 a2 a3 ... ak] // then "convergent(2)" is the rational number // a_(k-1) + 1 / ak // and "convergent(k)" is the rational number // a1 + 1/(a2+1/(...) ... ). function convergent(v, n) i := n; x := v[i]; i -:= 1; while i ge 1 do x := v[i] + 1/x; i -:= 1; end while ; return x; end function; // Compute P*{0,a/b} function ConvFromModSym(P, a, b, M) // "a/b = ", a/b; ZN:=quo; // ZN:=IntegerRing(M`N); if b eq 0 then // {0,oo} assert a ne 0; return ConvFromManinSymbol(M, P, [ZN|0,1]); end if; v := ContinuedFraction(a/b); p := 0; q := 1; pp := 1; qq := 0; ans := M!0; for j:= 1 to #v+1 do det := q*pp-p*qq; assert Abs(det) eq 1; g := [pp, p, qq, q]; if det lt 1 then g[2] *:= -1; g[4] *:= -1; end if; R := M`mlist`R; h := hom R | g[1]*R.1+g[2]*R.2, g[3]*R.1+g[4]*R.2>; hP := h(P); ans +:= ConvFromManinSymbol(M, hP, [ZN|g[3], g[4]]); if j le #v then p := pp; q := qq; cn := Rationals()!convergent(v,j); pp := Numerator(cn); qq := Denominator(cn); end if; end for; return ans; end function; intrinsic ConvFromModularSymbol(M::ModTupFld, alpha::SeqEnum, beta::SeqEnum) -> ModTupFldElt {Returns the modular symbol \{alpha,beta\} in M.} return ConvFromModularSymbol(M,); end intrinsic; intrinsic ConvFromModularSymbol(M::ModTupFld, Px::Tup) -> ModTupFldElt {Given a modular symbol P(X,Y)*\{alp,beta\}, FromModularSymbol returns the corresponding element of M. The input is a pair , where P is a polynomial in X and Y, and x=[[a,b],[c,d]] is a pair of elements of P^1(Q), where [a,b] <--> a/b and [c,d] <--> c/d.} // P(X,Y) in F[X,Y] is homogeneous of degree k-2. // Using the relation // P*{a/b,c/d}=P*{a/b,0}+P*{0,c/d}=-P*{0,a/b} + P*{0,c/d} // we break the problem into two subgroups. return -ConvFromModSym(Px[1],Px[2][1][1],Px[2][1][2],M) +ConvFromModSym(Px[1],Px[2][2][1],Px[2][2][2],M); end intrinsic; intrinsic ConvFromModularSymbol(M::ModTupFld, Px::SeqEnum) -> ModTupFldElt {Given a sequence of modular symbols P(X,Y)*\{alp,beta\}, FromModularSymbol returns the corresponding element of M. The input is a sequence of pairs , where P is a polynomial in X and Y, and x=[[a,b],[c,d]] is a pair of elements of P^1(Q), where [a,b] <--> a/b and [c,d] <--> c/d.} return &+[ConvFromModularSymbol(M, Px[i]) : i in [1..#Px]]; end intrinsic; intrinsic TestConv(M::ModTupFld) {Tests the modular symbols conversion routines. Must always return the basis for M.} a := [ConvToModularSymbol(M, M.i) : i in [1..Dimension(M)]]; b := [ConvFromModularSymbol(M, a[i]) : i in [1..Dimension(M)]]; if b ne Basis(M) then error b,"TestConv: Modular symbols conversion test failed!"; else printf "TestConv: test passed for M=%o_%o\n",Level(M),Weight(M); end if; end intrinsic; /********************************************************** INTEGRAL STRUCTURE: Algorithm: Compute basis for the *lattice* generated by enough generating Manin symbols [X^iY^(k-2-i),(u,v)]. **********************************************************/ intrinsic IntegralModularSymbols(M::ModTupFld) -> ModTupFld {Return the sub Z-module M_k(N,Z) of integral modular symbols.} assert Dimension(M) ne 0; if not assigned M`Z then // It suffices to mod out by the free // generating symbols X^iY^(k-2-i).(u,v) // after modding out only by the S(and I) relations. Q := M`quot; gens := [ Q`Scoef[i]*Q`Tquot[Q`Squot[i]] : i in [1..#Q`Squot] ]; n := #Q`Tgens; gens := [g[i] : i in [1..n], g in gens]; B := [M!v : v in Basis(Lattice(n,gens))]; M`Z := VectorSpaceWithBasis(B); end if; return M`Z; end intrinsic; intrinsic MkZ(M::ModTupFld) -> ModTupFld {} return IntegralModularSymbols(M); end intrinsic; ALPHABET := {@ "A", "B", "C", "D", "E", "F", "G", "H", "I", "J", "K", "L", "M", "N", "O", "P", "Q", "R", "S", "T", "U", "V", "W", "X", "Y", "Z" @}; alphabet := {@ "a", "b", "c", "d", "e", "f", "g", "h", "i", "j", "k", "l", "m", "n", "o", "p", "q", "r", "s", "t", "u", "v", "w", "x", "y", "z" @}; intrinsic ToIsogenyCode(n::RngIntElt) -> MonStgElt {Returns the n-th isogeny coding. The coding goes A,B,C,...,Z,AA,BB,...} if n eq 0 then return "?"; end if; return &cat[ALPHABET[((n-1)mod 26)+1]: i in [0..((n-1) div 26)]]; end intrinsic; intrinsic IsogenyCodeToInteger(s::MonStgElt) -> MonStgElt {Returns the number n so that s = ToIsogenyCode(n).} return 26*(#s-1)+Index(ALPHABET,s[1])+Index(alphabet,s[1]); end intrinsic; intrinsic SmallestPrimeNondivisor(N::RngIntElt) -> RngIntElt {} return SmallestPrimeNondivisor(N,2); end intrinsic; intrinsic SmallestPrimeNondivisor(N::RngIntElt,p::RngIntElt) -> RngIntElt {Return the smallest prime number ell not dividing N and such that ell >= p.} if not IsPrime(p) then p := NextPrime(p); end if; while N mod p eq 0 do p := NextPrime(p); end while; return p; end intrinsic; PRIMES := {@ p : p in [1..10000] | IsPrime(p) @}; intrinsic NthPrime(n::RngIntElt) -> RngIntElt {Returns the nth prime.} if n lt #PRIMES then return PRIMES[n]; end if; p:=2; for i in [2..n] do p := NextPrime(p); end for; return p; end intrinsic; intrinsic PrimePos(p::RngIntElt) -> RngIntElt {Returns the integer n so that p is the nth prime.} i := Index(PRIMES,p); if i ne 0 then return i; else ell := NextPrime(PRIMES[#PRIMES]); i := #PRIMES+1; end if; while ell lt p do ell := NextPrime(ell); i +:= 1; end while; return i; end intrinsic; intrinsic PrimePos(p::RngIntElt, N::RngIntElt) -> RngIntElt {Returns the position of the divisor p of N, or 0 if p does not divide N.} if N mod p ne 0 then return 0 ; end if; return Index([x[1] : x in Factorization(N)], p); end intrinsic; intrinsic ModularCharpoly(A::AlgMatElt) -> RngUPolElt {Fast, non-proved, charpoly.} assert Nrows(A) gt 0; if Type(A[1,1]) eq FldRatElt then return CharacteristicPolynomial(A : Al := "Modular", Proof:=false); end if; if Type(A[1,1]) eq FldPadElt then // TODO: need better error message here. M:=MatrixAlgebra(IntegerRing(Parent(A[1,1])),Degree(Parent(A))); return CharacteristicPolynomial(M!A); end if; return CharacteristicPolynomial(A); end intrinsic; intrinsic FactorCharpoly(A::AlgMatElt) -> SeqEnum {} return Factorization(ModularCharpoly(A)); end intrinsic; intrinsic IntegerKernel(A::ModMatFldElt) -> ModTupRng {Compute integer kernel of the not-necessarily-square matrix A.} n := Nrows(A); m := Ncols(A); d := Lcm([Integers()|Denominator(a) : a in Eltseq(A)]); B := RMatrixSpace(Integers(),n,m) ! (d*A); b := Basis(Kernel(B)); return VectorSpaceWithBasis([VectorSpace(Rationals(),n)!x: x in b]); end intrinsic; intrinsic IntegerKernel(A::ModMatRngElt) -> ModTupRng {Compute integer kernel of the not-necessarily-square matrix A.} n := Nrows(A); m := Ncols(A); B := RMatrixSpace(Integers(),n,m) ! A; b := Basis(Kernel(B)); return VectorSpaceWithBasis([VectorSpace(Rationals(),n)!x: x in b]); end intrinsic; intrinsic IntegerKernel(A::AlgMatElt) -> ModTupRng {Compute integer kernel of the matrix A.} if Type(A[1,1]) eq RngIntElt then return Kernel(A); end if; n := Nrows(A); d := LCM([Denominator(a) : a in Eltseq(A)]); B := MatrixAlgebra(Integers(),n) ! (d*A); K := Kernel(B); b := Basis(K); if #b eq 0 then return K; // 0 dimensional space. end if; return VectorSpaceWithBasis([VectorSpace(Rationals(),n)!x: x in b]); end intrinsic; intrinsic IntegerKernelZ(A::AlgMatElt) -> ModTupRng {Compute integer kernel of the matrix A, and return as integral space.} if Type(A[1,1]) eq RngIntElt then return Kernel(A); end if; n:= Nrows(A); d:= LCM([Denominator(a) : a in Eltseq(A)]); B:= MatrixAlgebra(Integers(),n) ! (d*A); return Kernel(B); end intrinsic; intrinsic LinearCombinations(Comb::SeqEnum, B::SeqEnum) -> SeqEnum {Compute the linear combinations of the elements of B defined by the elements of Comb.} n := #B; if n eq 0 then return B; end if; return [&+[B[i]*v[i] : i in [1..n]] : v in Comb]; end intrinsic; intrinsic IntegerKernelOn(A::AlgMatElt, B::ModTupFld) -> ModTupFld {} return IntegerKernelOn(A,Basis(B)); end intrinsic; intrinsic IntegerKernelOn(A::AlgMatElt, B::SeqEnum) -> ModTupFld {Suppose B is a basis for an n-dimensional subspace of some ambient space and A is an nxn matrix. Then A defines a linear transformation of the space spanned by B. This function returns the integer kernel of that transformation, as a subspace with basis.} K := IntegerKernel(A); if Dimension(K) eq 0 then V := VectorSpace(Parent(B[1][1]),#B); return sub; end if; return VectorSpaceWithBasis(LinearCombinations(Basis(K),B)); end intrinsic; intrinsic LatticeIndex(V::Lat, W::Lat) -> FldRatElt {Compute the lattice index [V:W].} Vec := VectorSpace(Rationals(),Degree(V)); B:=[Vec!b : b in Basis(V)]; C:=[Vec!c : c in Basis(W)]; Vspace := VectorSpaceWithBasis(B); D:=[Coordinates(Vspace, c) : c in C]; if Dimension(W) lt Dimension(V) then return 0; end if; n := #D; return Abs(Determinant(MatrixAlgebra(Rationals(),n)!(&cat D))); end intrinsic; /************************************************************* DECOMPOSITION: *************************************************************/ intrinsic ModSymIsRoot(A::ModTupFld) -> BoolElt {Returns true if and only if A is a root in the ModSymParent tree, i.e., if and only if A has no parent.} return not assigned A`M; end intrinsic; intrinsic ModSymParent(A::ModTupFld) -> ModTupFld {Returns the parent of the modular symbols factor A.} return A`M; end intrinsic; intrinsic IsNew (V::ModTupFld) -> BoolElt {} return V`isnew; end intrinsic; intrinsic OldFactor (V::ModTupFld) -> ModTupFld {} return V`old; end intrinsic; intrinsic Dimension(D::SeqEnum) -> SeqEnum {List the dimensions of the factors in D.} return [Dimension(D[i]) : i in [1..#D]]; end intrinsic; intrinsic Sha(D::SeqEnum) -> SeqEnum {List the Shas of factors in D.} return [Sha(Factor(D,i)) : i in [1..NumberOfNewforms(D)]]; end intrinsic; intrinsic NumberOfNewforms(D::SeqEnum) -> RngIntElt {Returns the number of newforms in the decomposition D.} i := 0; for d in D do if IsNew(d) and IsCuspidal(d) then i +:= 1; end if; end for; return i; end intrinsic; intrinsic NumberOfCuspidalFactors(D::SeqEnum) -> RngIntElt {Returns the number of cuspidal factors in the decomposition D.} i := 0; for d in D do if IsCuspidal(d) then i +:= 1; end if; end for; return i; end intrinsic; intrinsic Wq(D::SeqEnum) -> SeqEnum {List the Atkin-Lehner signs of the factors in D.} DecompWq(D); return [D[i]`wq : i in [1..#D]]; end intrinsic; intrinsic J0(N::RngIntElt) -> List {Compute J_0(N). (Really this is a decomposition of H1(X_0(N),Z)^+.} return Decomposition(ModularSymbols(N,2,+1)); end intrinsic; ///////////// SEMI-DECOMPOSITION ///////////////////////////////// /*********** A _fast_ and dirty decomposition, which doesn't identify oldform labels, doesn't compute the old spaces, and might not completely refine the decomposition. However this is _very_ useful when you already know what you are looking for at a given level. *************************/ intrinsic SemiDecomposition(M::ModTupFld) -> List {Compute the SemiDecomposition of M.} if GetCategory(M) eq ModSym then if not assigned M`semidecomp then M`semidecomp := SemiDecomposition(M,M,2,13); end if; return M`semidecomp; end if; end intrinsic; intrinsic SemiDecomposition(M::ModTupFld, stop::RngIntElt) -> List {Compute the SemiDecomposition of M.} return SemiDecomposition(M,M,2,stop); end intrinsic; intrinsic SemiDecomposition(M::ModTupFld, Vdual::ModTupFld, stop::RngIntElt) -> List {Compute the SemiDecomposition of the subspace Vdual inside Mdual.} return SemiDecomposition(M, Vdual, 2, stop); end intrinsic; intrinsic SemiDecomposition(M::ModTupFld, Vdual::ModTupFld, p::RngIntElt, stop::RngIntElt) -> List {Compute the SemiDecomposition of the subspace Vdual inside Mdual, starting with Tp.} p := SmallestPrimeNondivisor(M`N, p); T := FastTn(M, Vdual, p); D := [ VectorSpace(BaseField(M), Dimension(M)) ]; Remove(~D,1); for fac in FactorCharpoly(T) do f := fac[1]; a := fac[2]; fT := Evaluate(f,T); W := KernelOn(fT,Vdual); SetCategory(W,ModSymFac); if (Sign(M) ne 0 and a eq 1) or (Sign(M) eq 0 and a le 2) then // occurs w/ multiplicity 1, so new. W`isnew := true; W`N := M`N; W`M := M; Append(~D,W); elif p gt stop then if (IsPrime(M`N) and M`k eq 1) or M`N eq 1 then // must be new anyways. W`N := M`N; W`isnew := true; else // probably isn't new. W`isnew := false; // WARNING: in fact it could be TRUE here!!!! end if; W`M := M; W`N := 0; Append(~D,W); // need to use degen maps!! else for Sub in SemiDecomposition(M, W, NextPrime(p), stop) do Append(~D, Sub); end for; end if; end for; return D; end intrinsic; ///////////// END SEMI-DECOMPOSITION //////////////////////////// intrinsic FullCuspidalFactor(M::ModTupFld) -> ModTupFld {Returns the sum of all cuspidal factors of M.} // WARNING: My feeling is that there must be a much // more efficient algorithm to compute this; however, // I have *not* been able to find one! W := M; p := 2; while Dimension(W) gt Dimension(Sk(M)) do T := Tn(M,p); f := CharacteristicPolynomial(Restrict(T,Sk(M))); W := Kernel(Evaluate(f,Transpose(T))) meet W; p := NextPrime(p); end while; W`N := M`N; W`M := M; W`iscusp := true; return W; end intrinsic; ///////////// RIGOROUS DECOMPOSITION //////////////////////////// function NewDecompositionRecursive(M, Vdual, p) D := [ Parent(VectorSpace(BaseField(M), 0))| ]; if Dimension(Vdual) eq 0 then return D; end if; T := FastTn(M, Vdual, p); for fac in FactorCharpoly(T) do f := fac[1]; a := fac[2]; fT := Evaluate(f,T); W := KernelOn(fT,Vdual); SetCategory(W,ModSymFac); W`isnew := true; W`N := M`N; W`M := M; W`iscusp := W subset Sknewdual(M); appended := false; if (a eq 1 or (Sign(M) eq 0 and a eq 2 and W`iscusp)) then // irred. Append(~D,W); appended := true; else if p ge 13 then if HeckeSuperverbose then "Trying to prove irreducible with a random sum of Hecke operators."; end if; Tsum := &+[Random(1,10)*FastTn(M,W,n) : n in [1..p]]; fac2 := FactorCharpoly(Tsum); if #fac2 eq 1 and (fac2[1][2] eq 1) or (Sign(M) eq 0 and fac2[1][2] eq 2 and W`iscusp) then Append(~D,W); appended := true; end if; end if; if p ge 37 then if HeckeVerbose then "NewDecomposition: WARNING -- space not breaking up."; end if; end if; if not appended then if HeckeSuperverbose then "Newdecomp--recursing with p = ",NextPrime(p), "and dim = ",Dimension(W); end if; for WW in NewDecompositionRecursive(M, W, NextPrime(p)) do Append(~D, WW); end for; end if; end if; end for; return D; end function; intrinsic OldDecomposition(M::ModTupFld) -> SeqEnum {Compute the Hecke decomposition of the old subspace of M.} assert GetCategory(M) eq ModSym; if not assigned M`olddecomp then D := [ Parent(VectorSpace(BaseField(M), 0))| ]; N := Level(M); NN := Reverse([a : a in Divisors(N)| a ne N]); for i in [1..#NN] do Divs := Divisors(Integers()!(N/NN[i])); MM := OldModularSymbols(M,NN[i]); DD := NewDecomposition(MM); // now take all images of DD in M. for A in DD do Append(~D,&+[sub : d in Divs ]); D[#D]`M := M; D[#D]`N := N; D[#D]`isnew := false; D[#D]`old := A; D[#D]`iscusp:= A`iscusp; SetCategory(D[#D],ModSymFac); end for; end for; D := LabelFactors(D); M`olddecomp := D; end if; return M`olddecomp; end intrinsic; intrinsic Decomposition(M::ModTupFld) -> SeqEnum {Compute the Hecke decomposition of the *NEW* subspace of M.} return NewDecomposition(M); end intrinsic; intrinsic NewDecomposition(M::ModTupFld) -> List {Compute the Hecke decomposition of the *NEW* subspace of M.} if Characteristic(BaseField(M)) ne 0 then error "Decomposition not programmed in char. p"; end if; assert GetCategory(M) eq ModSym; if not assigned M`newdecomp then M`newdecomp := NewDecompositionRecursive(M,Mknewdual(M),2); M`newdecomp := LabelFactors(M`newdecomp); end if; return M`newdecomp; end intrinsic; intrinsic FullDecomposition(M::ModTupFld) -> List {Compute the Hecke decomposition of M.} if Characteristic(BaseField(M)) ne 0 then error "Decomposition not programmed in char. p"; end if; assert GetCategory(M) eq ModSym; if not assigned M`decomp then M`decomp := NewDecomposition(M) cat OldDecomposition(M); M`decomp := LabelFactors(M`decomp); end if; return M`decomp; end intrinsic; ////////////// END RIGOROUS DECOMPOSITION ///////////////////// intrinsic GetV(A::ModTupFld) ->ModTupFld {} if not assigned A`V then M := A`M; V := RSpace(M`F,Dimension(M)); p := SmallestPrimeNondivisor(M`N); while Rank(V) gt Dimension(A) do T := HeckeOperator(M,p); f := DecompCharpolyTp(A,p); fT := Evaluate(f,T); V meet:= Kernel(fT); p := SmallestPrimeNondivisor(M`N,p+1); end while; A`V := V; end if; return A`V; end intrinsic; intrinsic IsCuspidal(A::ModTupFld) -> BoolElt {} if not assigned A`iscusp then error "Can't determine if A is cuspidal."; end if; return A`iscusp; end intrinsic; intrinsic IsEisenstein(A::ModTupFld) -> BoolElt {} return not IsCuspidal(A); end intrinsic; intrinsic WqOn(A::ModTupFld, q::RngIntElt) -> AlgMatElt {Compute the action of the Atkin-Lehner involution Wq on A.} assert IsPrime(q) and A`N mod q eq 0; if not assigned A`wq then A`wq := []; for r in [p[1] : p in Factorization(A`N)] do w := Transpose(Wq(A`M,r)); Append(~A`wq,Restrict(w,A)); end for; end if; return A`wq[PrimePos(Factorization(q)[1][1],A`N)]; end intrinsic; intrinsic SignWN(A::ModTupFld) -> RngIntElt {Returns the sign of the Atkin-Lehner involution on A.} return &*[Integers()|WqOn(A,p)[1,1] where p is fac[1] : fac in Factorization(Level(A))]; end intrinsic; intrinsic DecompWq(D::SeqEnum) {} if assigned D[1]`wq then return; end if; M := D[1]`M; for i in [1..#D] do D[i]`wq := []; end for; for q in [p : p in Divisors(M`N)|IsPrime(p)] do w:=Transpose(Wq(M,q)); for i in [1..#D] do if IsNew(D[i]) then Append(~D[i]`wq,Restrict(w,D[i])); end if; end for; end for; end intrinsic; intrinsic DecompFast(A::ModTupFld, maxp::RngIntElt) {} if not assigned A`fast then A`fast := [ ]; end if; for p in [NthPrime(#A`fast+1)..maxp] do if IsPrime(p) then Append(~A`fast, FastTn(A`M,A, p)); end if; end for; end intrinsic; intrinsic DecompFast(D::SeqEnum, maxp::RngIntElt, NewOnly::BoolElt) {} for i in [1..#D] do if not NewOnly or IsNew(D[i]) then DecompFast(D[i], maxp); end if; end for; end intrinsic; intrinsic LabelHelper(prec::RngIntElt, D::SeqEnum) -> SeqEnum {} prec := NextPrime(prec-1); nprimes := PrimePos(prec); M := D[1]`M; DecompFast(D,prec,true); was := []; for i in [1..#D] do if IsNew(D[i]) then /* Definition of ordering (_extends_ Cremona to higher dimension). 'Smallest' forms come *first* in output lists: (1) If they have different levels then the larger level is smallest. (2) If both are old, compare by iso class. (3) Smallest dimension is lesser (to agree with Cremona). (4) order by Wq eigenvalues, starting with *smallest* p|N and with "+" being less than "-" (5) order by abs(trace(a_p)), with + one being smaller in the the event of a tie. The Eisenstein factors are placed at the end, by setting their dimension during the sort to 10^9. */ if IsTrivial(Character(M)) then seq := [Dimension(D[i])] cat [-D[i]`wq[j][1,1] : j in [1..#D[i]`wq]] cat // give absolute value of trace(a_p), &cat[ // then the sign of trace(a_p). Characteristic(M`F) eq 0 select [Abs(Norm(t)),-Sign(t)] else [t,t] where t is Trace(D[i]`fast[j]) : j in [1..nprimes] ] cat [i]; else /* Away from the case of J_0(N), order instead by dimension and then the trace down to Z of the Hecke eigenvalues a_p. */ seq := [Dimension(D[i])] cat [Trace(Trace(DecompFastTp(D[i],p))) : p in [1..nprimes] | IsPrime(p)] cat [i]; end if; if not IsCuspidal(D[i]) then seq[1] := 10^9; if Characteristic(M`F) ne 0 then seq[1] := -1; // at least near the end... end if; end if; Append(~was,seq); end if; end for; return was; end intrinsic; intrinsic LabelFactors(D::SeqEnum) -> SeqEnum {Label the factors in the decomposition D, then return the sorted decomposition.} if #D eq 0 then return D; end if; for i in [1..#D] do delete D[i]`iso; end for; M := D[1]`M; if IsTrivial(Character(M)) then DecompWq(D); end if; prec := 5; if HeckeVerbose then printf "Sorting and labeling factors at level %o.\n",Level(M); end if; finished := false; repeat was := LabelHelper(prec,D); if #was eq 0 then break; end if; // remove last entry (which identifies factor) and check // for duplicates. n := #was[1]; was2 := [Remove(x,n) : x in was]; if #Set(was2) eq #was then // no duplicates was := Sort(was); for i in [1..#was] do D[Integers()!was[i][n]]`iso := i; end for; finished := true; end if; prec +:= 5; until finished; // old factors start := #was+1; for i in [1..#D] do if not assigned D[i]`iso then assert not IsNew(D[i]); D[i]`iso := start; start +:= 1; end if; end for; // Sort the entries of the SeqEnum D, so that D[i]`iso = i. E := D; for i in [1..#D] do E[D[i]`iso] := D[i]; end for; return E; end intrinsic; intrinsic Print(A::ModTupFld : abbrev:=false) {} if GetCategory(A) eq ModSym then if IsTrivial(A`eps) then printf "Full Modular symbols space of dimension %o, level %o and weight %o ", Dimension(A), Level(A), Weight(A); else printf "Full Modular symbols space of dimension %o, level %o, weight %o and " cat "character %o ", Dimension(A), Level(A), Weight(A), DirichletCharacterImages(A`eps); end if; if IsPlusQuotient(A) then printf "(plus quotient)"; elif IsMinusQuotient(A) then printf "(minus quotient)"; end if; printf "\n"; end if; if GetCategory(A) eq ModSymFac then eps := IsTrivial(Character(A)) select "" else "eps"; type := IsCuspidal(A) select "cuspidal" else "eisenstein"; if IsNew(A) then if IsCuspidal(A) then d := DimensionAbelianVariety(A); else d := Dimension(A); end if; if abbrev then printf "%3ok%1o%o%1o: dim =%2o %10o\n", Level(A), Weight(A), eps, ToIsogenyCode(IsogenyClass(A)), Dimension(A), type; else printf "New %o factor: dimension %o (%o), " cat "level %o, weight %o (%o%o)\n", type, Dimension(A), d, Level(A), Weight(A), Level(A), ToIsogenyCode(IsogenyClass(A)); end if; else if abbrev then printf "%3ok%1o%o%1o: dim =%2o %10o (old: %3o%1o)\n", Level(A), Weight(A), eps, ToIsogenyCode(IsogenyClass(A)), Dimension(A), type, Level(OldFactor(A)), ToIsogenyCode(IsogenyClass(OldFactor(A))), Dimension(A); else printf "Old %o factor: dimension %o, level %o, weight %o " cat "(%o%o --> %o%o)\n", type, Dimension(A), Level(A), Weight(A), Level(OldFactor(A)), ToIsogenyCode(IsogenyClass(OldFactor(A))), Level(A), ToIsogenyCode(IsogenyClass(A)); end if; end if; end if; end intrinsic; intrinsic Print(D::SeqEnum) {Print the factors in D.} if #D lt 1 then "Empty collection of modular symbols factor.\n"; return; end if; "Modular symbols factors:"; if not IsTrivial(Character(D[1])) then "(nontrivial character)"; end if; if IsogenyClass(D[1]) eq 0 then "(The isogeny classes have not been labeled.)"; for i in [1..#D] do Print(D[i]); end for; return; end if; for i in [1..#D] do for j in [1..#D] do if IsogenyClass(D[j]) eq i then Print(D[j] : abbrev:=true); break; end if; end for; end for; end intrinsic; intrinsic NumOldFactors(D::SeqEnum) -> RngIntElt {Return the number of new factors of the decomposition D.} return #[i : i in [1..#D] | not D[i]`isnew]; end intrinsic; intrinsic NewFactors(D::SeqEnum) -> List {Return the indexes of new factors in the list.} return [i : i in [1..#D] | D[i]`isnew]; end intrinsic; intrinsic NewformFactors(D::SeqEnum) -> List {Return the indexes of new factors in the list.} return [i : i in [1..#D] | D[i]`isnew and IsCuspidal(D[i])]; end intrinsic; intrinsic NumNewFactors(D::SeqEnum) -> RngIntElt {Return the number of new factors of the decomposition D.} return #NewFactors(D); end intrinsic; intrinsic NumNewformFactors(D::SeqEnum) -> RngIntElt {Return the number of new factors of the decomposition D.} return #NewformFactors(D); end intrinsic; intrinsic newform(s::MonStgElt) -> ModTupFld {} return ModularFactor(s); end intrinsic; intrinsic ModularFactor(s::MonStgElt) -> ModTupFld {Returns the modular factor "Nk[Weight][IsogenyClass]", where N is the level, k is the weight, and [IsogenyClass] is a letter such as "A", "B", etc. An alternative input format is "N[IsogenyClass]" in which case the weight is assumed to be two.} i := 1; sN := ""; while i le #s and s[i] ge "0" and s[i] le "9" do sN := sN cat s[i]; i +:= 1; end while; N := StringToInteger(sN); if i gt #s then return ModularFactor(N,2,1); end if; if s[i] eq "k" then i +:= 1; sk := ""; while i le #s and s[i] ge "0" and s[i] le "9" do sk := sk cat s[i]; i +:= 1; end while; k := StringToInteger(sk); else k := 2; end if; if i gt #s then iso := 1; else sIso := &cat [s[j] : j in [i..#s]]; iso := IsogenyCodeToInteger(sIso); end if; return ModularFactor(N,k,iso); end intrinsic; intrinsic ModularFactor(N::RngIntElt, k::RngIntElt, iso::RngIntElt) -> ModTupFld {} return ModularFactor(ModularSymbols(N,k),iso); end intrinsic; intrinsic ModularFactor(N::RngIntElt, iso::RngIntElt) -> ModTupFld {} return ModularFactor(ModularSymbols(N),iso); end intrinsic; intrinsic ModularFactor(M::ModTupFld, iso::RngIntElt) -> ModTupFld {} D := Decomposition(M); if iso gt #D then D := FullDecomposition(M); end if; return ModularFactor(D, iso); end intrinsic; intrinsic ModularFactor(D::SeqEnum, iso::RngIntElt) -> ModTupFld {Returns the index into the D list of the new factor with the indicated isogeny class. If the iso classes have not yet been computed yet, they are computed. If there is no factor of class iso this function returns 0.} return Factor(D,iso); end intrinsic; intrinsic Factor(M::ModTupFld, s::MonStgElt) -> ModTupFld {} t:=&cat[s[i] : i in [1..#s] | not ("0" le s[i] and s[i] le "9")]; return Factor(M,IsogenyCodeToInteger(t)); end intrinsic; intrinsic Factor(M::ModTupFld, iso::RngIntElt) -> ModTupFld {} return Factor(Decomposition(M),iso); end intrinsic; intrinsic Factor(D::SeqEnum, s::MonStgElt) -> ModTupFld {} t:=&cat[s[i] : i in [1..#s] | not ("0" le s[i] and s[i] le "9")]; return Factor(D,IsogenyCodeToInteger(t)); end intrinsic; intrinsic Factor(D::SeqEnum, iso::RngIntElt) -> ModTupFld {Returns the index into the D list of the new factor with the indicated isogeny class. If the iso classes have not yet been computed yet, they are computed. If there is no factor of class iso this function returns 0.} requirerange iso,1,#D; for i in [1..#D] do if not assigned D[i]`iso then D:=LabelFactors(D); end if; if D[i]`iso eq iso then return D[i]; end if; end for; printf "There is no factor %o.",iso; error ""; end intrinsic; intrinsic Factor(M::ModTupFld, A::ModTupFld) -> ModTupFld {Give the subspace A of the DUAL of M, attach necessary data to make it a factor of M. For this to make sense, A should be Hecke stable.} A`M := M; A`N := M`N; A`theCategory := ModSymFac; return A; end intrinsic; intrinsic DecompFastTp(A::ModTupFld, p::RngIntElt) -> AlgMatElt {} DecompFast(A,p); // make sure Tp is known. return A`fast[PrimePos(p)]; end intrinsic; intrinsic DecompCharpolyTp(A::ModTupFld, p::RngIntElt) -> RngUPolElt {} return ModularCharpoly(DecompFastTp(A,p)); end intrinsic; intrinsic DecompZ(A::ModTupFld) -> ModTupFld {} M := A`M; if not assigned A`Z then Z := RSpace(Integers(),Dimension(Sk(M))); Zdual := RSpace(Integers(),Dimension(Sk(M))); p := 2; while Rank(Z) gt Dimension(A) do T := HeckeOperatorSkZ(M,p); f := DecompCharpolyTp(A,p); fT := Evaluate(f,T); Z meet:= IntegerKernelZ(fT); Zdual meet:= IntegerKernelZ(Transpose(fT)); p := SmallestPrimeNondivisor(M`N,p+1); end while; A`Z := Z; A`Zdual := Zdual; end if; return A`Z; end intrinsic; intrinsic DecompZdual(A::ModTupFld) -> ModTupFld {} if not assigned A`Zdual then DecompZ(A); // compute A`Zdual as well. end if; return A`Zdual; end intrinsic; /************************************************************* CONGRUENCES: *************************************************************/ intrinsic CuspOrder(A::ModTupFld) -> RngIntElt {Compute the order in A of the difference (0)-(oo).} return CuspOrder(A,[0,1],[1,0],0); end intrinsic; intrinsic CuspOrder(A::ModTupFld, alpha::SeqEnum, beta::SeqEnum) -> RngIntElt {Compute the order in A of the difference (x)-(y) of the cusps defined by the 2-tuples alpha and beta.} return CuspOrder(A,alpha,beta,0); end intrinsic; intrinsic TorsionBound(A::ModTupFld) -> RngIntElt, SeqEnum {Computes the following upper bound on the torsion subgroup of A: gcd \{ #A(F_p) : 3 <= p <= 19, p not dividing N \}} return TorsionBound(A,19); end intrinsic; intrinsic TorsionBound(A::ModTupFld, maxp::RngIntElt) -> RngIntElt {Computes the following upper bound on the torsion subgroup of A: gcd \{ #A(F_p) : 3 <= p <= maxp, p not dividing N \}} assert Weight(A) eq 2 and IsTrivial(Character(A)); S := [ p : p in [3..maxp] | IsPrime(p) and Level(A) mod p ne 0]; if GetCategory(A) eq ModSym then seq := [Integers()|Evaluate(CharacteristicPolynomial( Restrict(Tn(A,p),Sk(A))),p+1) : p in S]; else seq := [Integers()|Evaluate(DecompCharpolyTp(A, p),p+1) : p in S]; end if; bound := Gcd(seq); if Sign(A) eq 0 then b := Integers()!Round(Sqrt(bound)); assert b^2 eq bound; bound := b; end if; return bound, seq; end intrinsic; intrinsic ScaledRationalPeriodMap(A::ModTupFld) -> ModMatFldElt {Compute the rational period mapping. The period mapping is scaled so that the integral modular symbols SkZ are taken surjectively onto the lattice Z^d. Note: the choice of rational period mapping is well-defined only once A is created, and can be different if A is created again; i.e., it is with respect to the non-canonical basis used internally to define A.} d := Dimension(A); P := Transpose(RMatrixSpace(A`M`F, d, Degree(A))!Basis(A)); Z := SkZ(A`M); PofZ := [z*P : z in Basis(Z)]; B := Basis(Lattice(RMatrixSpace(Rationals(), #PofZ, d)!PofZ)); B := [VectorSpace(A`M`F, d)|b : b in B]; I := (RMatrixSpace(A`M`F, d, d)!B)^(-1); return P*I; end intrinsic; intrinsic CuspOrder(A::ModTupFld, alpha::SeqEnum, beta::SeqEnum, i::RngIntElt) -> RngIntElt {Compute the order in A of the cusp defined by X^i*Y^\{k-2-i\}*\{alpha,beta\}.} R := PolynomialRing(A`M`F,2); modsym := ; v := ConvFromModularSymbol(A`M, modsym); w := v*ScaledRationalPeriodMap(A); return Lcm([Denominator(w[i]): i in [1..Degree(w)]]); end intrinsic; intrinsic ModularKernel(A::ModTupFld) -> RngIntElt, ModTupRng {} Z := DecompZ(A); Zdual := DecompZdual(A); if Rank(Z) eq 0 then return 0,[]; end if; D := MatrixAlgebra(Integers(),Dimension(A))! [InnerProduct(v, w) : v in Basis(Z), w in Basis(Zdual)]; S := SmithForm(D); return Determinant(S), [S[i,i] : i in [1..Ncols(S)]]; end intrinsic; intrinsic ModularDegree(A::ModTupFld) -> ModTupRng {Compute the modular degree of A. This is the degree of the canonical map from the dual of A to A. When A corresponds to an elliptic curve, this is the SQUARE of what is frequently referred to as the "modular degree" in the literature.} if not assigned A`moddeg then A`moddeg := ModularKernel(A); end if; return A`moddeg; end intrinsic; intrinsic CongruenceR(A::ModTupFld, prec::RngIntElt) -> RngIntElt, ModTupRng {Returns the congruence number r = [S_k(N,Z) : (This) + (Complement)], along with the structure of the quotient; q-expansions are computed to precision prec.} Q1 := qIntegralBasis(A, prec); M := ModSymParent(A); D := FullDecomposition(M); Q2 := SaturatePolySeq(&cat[qIntegralBasis(B,prec) : B in D | IsCuspidal(B) and (B`iso ne A`iso)]); if #Q2 eq 0 then return 1, []; end if; return IndexInSaturation(Q1 cat Q2); end intrinsic; function TorsionQuotient(L) // L is a free abelian group in some ambient Z^n. // This function returns the structure of the // finite abelian group (Z^n/L)_tor. D := RMatrixSpace(Integers(),Dimension(L),Degree(L))!Basis(L); S := SmithForm(D); return [S[i,i] : i in [1..Min(Nrows(S),Ncols(S))] | S[i,i] gt 1]; end function; intrinsic AbelianIntersection(A::ModTupFld, B::ModTupFld) -> ModTupRng {Computes the invariants of the group-theoretic intersection of the abelian varieties corresponding to the modular factors A and B. If the factors are not distinct, then the 0 group is returned.} assert ModSymParent(A) eq ModSymParent(B); if not IsCuspidal(A) or not IsCuspidal(B) then return [Integers()|]; end if; ZA := DecompZ(A); ZB := DecompZ(B); assert Rank(ZA) ne 0; assert Rank(ZB) ne 0; Z := ZA + ZB; return TorsionQuotient(Z); end intrinsic; intrinsic AbelianIntersection(S::SeqEnum) -> ModTupRng {Returns the group-theoretic intersection of a sequence S of modular factors.} n := #S; ZS := [DecompZ(S[i]) : i in [1..n]]; m := Degree(ZS[1]); D := RMatrixSpace(Integers(),&+[Dimension(ZS[i]) : i in [1..n]], (n-1)*m)!0; /////////////////////////////////////////////////////////////////// // Fill up the matrix D with generators for the image of the map // // f: B_1 + ... + B_n --> Z^m + ... + Z^m // // given by // // f(x_1,...,x_n) = (x_1-x_2, x_2-x_3, ..., x_{n-1}-x_n). // /////////////////////////////////////////////////////////////////// B := Basis(ZS[1]); for i in [1..#B] do for c in [1..m] do D[i,c] := B[i][c]; end for; end for; s := #B+1; t := 0; for j in [2..#ZS-1] do B := Basis(ZS[j]); for i in [1..#B] do for c in [1..m] do D[s,t+c] := -B[i][c]; D[s,t+c+m] := B[i][c]; end for; s +:= 1; end for; t +:= m; end for; B := Basis(ZS[#ZS]); for i in [1..#B] do for c in [1..m] do D[s,t+c] := -B[i][c]; end for; s +:= 1; end for; S := SmithForm(D); return [S[i,i] : i in [1..Min(Nrows(S),Ncols(S))] | S[i,i] gt 1]; end intrinsic; /************************************************************* Rational parts of L-functions *************************************************************/ intrinsic PeriodIndex(W::ModTupFld, A::ModTupFld, Plus::BoolElt) -> FldRatElt {Compute "[Phi(SkZ^+):Phi(W)]". If Plus is false, compute instead "[Phi(SkZ^-):Phi(W)]". It should be the case that W tensor Q is contained in SkZ^+ tensor Q (or SkZ^- tensor Q, when Plus is false).} M := A`M; if Plus then S := SkZPlus(M); else S := SkZMinus(M); end if; LW := PeriodImage(W,A); if Dimension(LW) eq 0 then return Rationals()!0; end if; LS := PeriodImage(S,A); assert Dimension(LS) ne 0; return LatticeIndex(LS, LW); end intrinsic; intrinsic PeriodImage(W::ModTupFld, A::ModTupFld) -> . {Returns the lattice spanned by the basis of W under the period map defined by A.} if Dimension(W) eq 0 then return VectorSpaceWithBasis([A|]); end if; n := #Eltseq(Basis(A)[1]); W := Basis(Lattice(RMatrixSpace(Rationals(), Dimension(W), n)!Basis(W))); Phi := Basis(Lattice(RMatrixSpace(Rationals(), Dimension(A), n)!Basis(A))); PhiW := [InnerProduct(x, y) : x in Phi, y in W]; MatW := RMatrixSpace(Rationals(), #W, #Phi)!PhiW; if MatW eq 0 then return RowSpace(MatW); end if; return Lattice(MatW); end intrinsic; intrinsic LRatio(A::ModTupFld) -> FldRatElt {} return LRatio(A,1); end intrinsic; intrinsic LRatio(A::ModTupFld, s::RngIntElt) -> FldRatElt {For an integer s in the critical strip (1<=s<=k-1), return the quotient L_A(s)*(s-1)! / (2pi)^(s-1)*Omega, which is a rational number. Here Omega is the volume of A(R), if s is odd, and the volume of the -1 eigenspace for conjugation.} M := A`M; if not (1 le s and s le Weight(A)-1) then "LRatioOddPart: s =",s," k =",Weight(A); error "However, s must be a critical value, i.e., an integer between 1 and k-1."; end if; if not assigned A`LRatio then A`LRatio := [Rationals()!-1 : i in [1..M`k-1]]; end if; if true or A`LRatio[s] eq -1 then k := Weight(A); A`LRatio[s] := PeriodIndex(WindingSubmodule(M,s), A, // the following condition is true iff we take the image of // the + part instead of the image of the - part of SkZ. (IsOdd(s) and IsEven(k)) or (IsEven(s) and IsOdd(k)) ); end if; return A`LRatio[s]; end intrinsic; intrinsic LRatioOddPart(A::ModTupFld) -> FldRatElt {Returns the odd part of the rational part of L(A,s); potentially faster than LRatio.} return LRatioOddPart(A,1); end intrinsic; intrinsic LRatioOddPart(A::ModTupFld, s::RngIntElt) -> FldRatElt {} /********************************************************************** Let e=e_s = -X^(s-1)Y^(k-2-s+1){0,\infty} Let v be the dual eigenvector with sign (-1)^(s-1). If = 0, the BSD value is 0 and we are done so assume =/= 0. This function computes Norm() ------------------ * [Z[x]:Z[alp]] Vol(L) plus or minus where L is the lattice generated by the row vectors (,,...,) Here v = v_1 + x*v_2 + x^2*v_3 + ... + x^{n-1}*v_n, x a root of f(x), and c varies through enough modular symbols to generate H_1(X_0(N),Z). Notice that Norm() depends only on f and the parity of s, and [Z[x]:Z[alp]] depends only on f. Thus for all but two values of k we need only compute Norm(). ************************************************************************/ if not IsCuspidal(A) then "LRatioOddPart: A must be cuspidal."; end if; if not (1 le s and s le Weight(A)-1) then "LRatioOddPart: s =",s," k =",Weight(A); error "However, s must be a critical value, i.e., an integer between 1 and k-1."; end if; M := A`M; if not assigned A`LRatioOdd then A`LRatioOdd := [M`F|-1 : i in [1..M`k-1]]; end if; if A`LRatioOdd[s] eq -1 then eig := Eigenvector(A, IsEven(Weight(A)) select (-1)^(s-1) else (-1)^s ); w := WindingElement(M,s); eigw := 0; for i in [1..Degree(eig)] do if w[i] ne 0 then eigw +:= w[i]*eig[i] ; end if; end for; if eigw eq 0 then A`LRatioOdd[s] := 0 ; else Nrm := Abs(Norm(eigw)); // this is a number in the base field. if Type(M`F) ne FldRat then A`LRatioOdd[s] := Nrm; if HeckeVerbose then "LRatio: warning L-value only up to =/= 0 scalar."; end if; return Nrm; end if; // Q = R/g, where R=F[x]. Q := Parent(eig[1]); if Type(Q) eq FldRat then A`ZxZalp := 1; n := 1; vi := [eig]; V := VectorSpace(Rationals(), n); else g := Modulus(Q); n := Degree(g); V := VectorSpace(Rationals(), n); R := PreimageRing(Q); if not assigned A`ZxZalp then sturm := HeckeBound(M); f := qEigenform(A,sturm); A`ZxZalp := Volume([ V![Coefficient(h,i): i in [0..n-1]] where h is R!Coefficient(f,j) : j in [1..sturm]]); end if; // the ith entry of vi is the coefficient of x^(i-1). vi := [[Coefficient(R!eig[j],i) : j in [1..Degree(eig)]] : i in [0..n-1]]; end if; if s mod 2 eq 0 then if not assigned A`VolLEven then Z := Basis(SkZ(M)); A`VolLEven := Volume([V| [DotProd(z,vi[i]) : i in [1..n]] : z in Z]); end if; VolL := A`VolLEven; else if not assigned A`VolLOdd then Z := Basis(SkZ(M)); A`VolLOdd := Volume([V| [DotProd(z,vi[i]) : i in [1..n]] : z in Z]); end if; VolL := A`VolLOdd; end if; A`LRatioOdd[s] := Nrm * A`ZxZalp / VolL; if IsOdd(s) then A`LRatioOdd[s] /:= RealTamagawa(A); else A`LRatioOdd[s] /:= ImaginaryTamagawa(A); end if; end if; end if; A`LRatioOdd[s] := OddPart(A`LRatioOdd[s]); return A`LRatioOdd[s]; end intrinsic; //////////////////////////////////////////////////////////////////////////// // VANISHING OF TWISTS. // // (currently only prime level.) // //////////////////////////////////////////////////////////////////////////// intrinsic LEpsilonVanishes(A::ModTupFld, eps::Rec) -> BoolElt {Returns true iff L(A,eps,1) = 0.} p := Conductor(A); if Weight(A) ne 2 or not IsTrivial(Character(A)) or not IsNew(A) or Modulus(eps) ne p or not IsPrime(p) then "Currently only implemented for new factors A of J_0(p),"; error "and for eps a character modulo p."; end if; F := FieldOfValues(eps); M := ModSymParent(A); V := VectorSpace(F,Dimension(M)); if not IsTrivial(eps) then e := &+[F!(Evaluate(eps,a))* V!ConvFromManinSymbol(M,[IntegerRing(p)|-a,1]) : a in [1..p-1]]; else e := V!WindingElement(M); end if; if Type(F!0) eq FldPrElt then for b in Basis(A) do if Abs(InnerProduct(V!b,e)) gt 10^(-5) then return false; end if; end for; else for b in Basis(A) do if InnerProduct(V!b,e) ne 0 then return false; end if; end for; end if; return true; end intrinsic; intrinsic LEpsilonVanishing(A::ModTupFld) -> SeqEnum {Returns the degrees d dividing p-1 such that L(A,eps,1) = 0.} p := Conductor(A); assert IsPrime(p); v := [d : d in Divisors(p-1) | LEpsilonVanishes(A,DirichletCharacter(p,[d]))]; return v; end intrinsic; /************************************************************ * * * COMPONENT GROUPS * * * ************************************************************/ ///////////////////////////////////////////////////////////////// // The "XGroup" functions depend on the availability of // // David Kohel's functions: // // BrandtHecke, BrandtModules, BrandTheta // // (Whereas, the Mestre method code is all here.) // ///////////////////////////////////////////////////////////////// XVerbose := false; XSuperverbose := false; declare attributes ModED: B, TM, // theta series, to precision prec. p, M, prec; intrinsic XGroup(p::RngIntElt) -> ModED {Return the character group X of the toric part of the closed fiber of the Neron model of J_0(p) over F_p.} return XGroup(p,1); end intrinsic; intrinsic XGroup(p::RngIntElt, M::RngIntElt) -> ModED {Return the character group X of the toric part of the closed fiber of the Neron model of J_0(pM) over F_p. [RIGHT NOW RETURNS THE EISENSTEIN PART AS WELL.]} tm := Cputime(); if not IsPrime(p) then error "CharGroup: ", "p = ", p, " is not prime."; end if; if not (M ge 1 and Gcd(p,M) eq 1) then error "CharGroup: ", "(p,M)=(", p,",",M,")=",Gcd(p,M),">1."; end if; if XVerbose then if M gt 1 then printf "Computing character group of torus of J_0(%o*%o)/F_%o.\n", p, M, p; else printf "Computing character group of torus of J_0(%o).\n", p; end if; t:=Cputime(); end if; B := BrandtModules(p,M); X := RModule(Integers(), Ncols(B[1][1])); X`B := B; if XVerbose then Cputime(t), " seconds."; end if; X`p:=p; X`M:=M; return X; end intrinsic; intrinsic XGroup_p(X::ModED) -> RngIntElt {Return the characteristic p of the character group of X_0(pM)/F_p.} return X`p; end intrinsic; intrinsic XGroup_M(X::ModED) -> RngIntElt {Return the level M of the character group of X_0(pM)/F_p.} return X`M; end intrinsic; intrinsic TXn(X::ModED, n::RngIntElt) -> AlgMatElt {Compute the n-th Hecke operator on the toric part (+Eisenstein series) of the closed fiber of the Neron model of J_0(pM) over F_p. The second argument returned is a diagonal matrix which is double the weight vector.} if not assigned X`TM then if XSuperverbose then "TXn: Computing q-expansion of theta series to",Max(n,17), " terms."; end if ; X`TM := BrandtTheta(X`B, Max(n,17)); X`prec := Max(n,17); elif n gt X`prec then m := n+Ceiling(n/4); if XSuperverbose then "TXn: Computing q-expansion of theta series to",m, " terms."; t := Cputime(); end if ; X`TM := BrandtTheta(X`B,m); if XSuperverbose then "finished. Time = ",Cputime(t); end if ; X`prec := Max(n,m); end if; T,W := BrandtHecke(X`TM, n); return MatrixAlgebra(Integers(), Ncols(T))!T; end intrinsic; intrinsic MonodromyWeights(X::ModED) -> SeqEnum {Returns the monodromy weights.} dummy := TXn(X,2); // so that TM is defined. TM := X`TM; T,W := BrandtHecke(TM,1); return [Integers()!(W[i,i]/2) : i in [1..Ncols(W)]]; end intrinsic; // Using quaternion algebras module. intrinsic PhiX_and_mX(X::ModED, V::ModTupRng) -> RngIntElt,RngIntElt {Computes the quantities PhiX and mX associated to V.} dummy := TXn(X,2); // so that TM is defined. TM := X`TM; n := Degree(Parent(TM)); m := Dimension(V); T1, W1 := BrandtHecke(TM,1); L := LatticeWithGram(W1/2); B := [ L.i - L.n : i in [1..n-1] ]; C := [ L!v : v in Basis(V) ]; S := Eltseq( LLL( RMatrixSpace(Integers(),n-1,m)! &cat[ [ InnerProduct(u,v) : v in C ] : u in B ]) ); PhiX := AbsoluteValue( Determinant( MatrixAlgebra(Integers(),m)! [ S[i] : i in [1..m^2] ] ) ); disc := Determinant( MatrixAlgebra(Integers(),m)! [ InnerProduct(u,v) : u in C, v in C ] ); mX := AbsoluteValue(disc/PhiX); return Integers()!PhiX, Integers()!mX; end intrinsic; intrinsic XDimension(X::ModED) -> RngIntElt {Returns the dimension of the character group X.} return Ncols(TXn(X,2)); end intrinsic; intrinsic XGroup(M::ModTupFld, p::RngIntElt) -> ModED {Returns the character group of the toric part of the closed fiber at p of the space M of modular symbols. This only makes sense when p exactly divides the level of M and M has weight two.} if not assigned M`k then return XGroupV(M,p); // this is probably what was really wanted. end if; if M`k ne 2 then error "XGroup: the weight must be 2."; end if; if not IsPrime(p) or Valuation(M`N,p) ne 1 then error "XGroup: invalid value of p."; end if; fac := {@ x[1] : x in Factorization(M`N) @}; if not assigned M`X then M`X := SequenceToList([0 : i in [1..#Factorization(M`N)]]); end if; i := Index(fac, p); if Type(M`X[i]) eq RngIntElt then M`X[i] := XGroup(p,Integers()!(M`N/p)); end if; return M`X[i]; end intrinsic; intrinsic XGroup(M::ModTupFld) -> ModED {Returns the character group corresponding to the largest prime p which exactly divides the level. If no such p exists, an error results.} for x in Reverse(Factorization(M`N)) do if x[2] eq 1 then return XGroup(M, x[1]); end if; end for; error "N has no square free divisors."; end intrinsic; intrinsic XGroupV(A::ModTupFld) -> ModTupFld {Compute XGroup for xmallest prime exactly dividing the level.} for x in Reverse(Factorization(A`N)) do if x[2] eq 1 then return XGroupV(A, x[1]); end if; end for; error "N has no square free divisors."; end intrinsic; intrinsic XGroupV(A::ModTupFld, p::RngIntElt) -> ModTupRng {Return the factor of the character group corresponding to the p-adic rigid analytic optimal and co-optimal quotient associated to A.} if (not IsCuspidal(A)) or (not A`M`k eq 2) or (not IsNew(A)) then error "XGroupV: weight must be 2 and A must be cuspidal and new."; end if; i := PrimePos(p, A`N); if not assigned A`X then A`X := SequenceToList([0 : j in [1..#Factorization(A`M`N)]]); end if; if Type(A`X[i]) eq RngIntElt then Z := XGroup(A`M,p); A`X[i] := RSpace(Integers(),XDimension(Z)); d := DimensionAbelianVariety(A); p := 2; while Dimension(A`X[i]) gt d do T := TXn(Z, p); // compute p-th Hecke operator. f := DecompCharpolyTp(A, p); A`X[i] := A`X[i] meet IntegerKernelZ(Evaluate(f,T)); p := NextPrime(p); if p ge 23 and HeckeVerbose then "XGroupV: Warning -- space taking more than 23 Hecke"; " operators to break up..."; end if; end while; end if; return A`X[i]; end intrinsic; intrinsic PhiX(A::ModTupFld) -> RngIntElt {Compute the p-adic modular degree of A, at the largest prime which exactly divides the level.} for x in Reverse(Factorization(A`N)) do if x[2] eq 1 then return PhiX(A,x[1]); end if; end for; error "N has no square free divisors."; end intrinsic; intrinsic PhiX(A::ModTupFld, p::RngIntElt) -> RngIntElt {Compute the order of the image of the component group of J_0(N) in the component group of A, all at p.} if not IsNew(A) then error "A must be new."; end if; if Valuation(A`N,p) ne 1 then error "p must exactly divide the level."; end if; if not assigned A`xdata then A`xdata := SequenceToList([[0,0] : i in [1..#Factorization(A`M`N)]]); end if; i := PrimePos(p,A`N); if A`xdata[i] eq [0,0] then if IsPrime(Conductor(A)) then // can use the much faster Mestre module. x, y := PhiX_and_mX(MestreGroup(ModSymParent(A)),MestreGroupV(A)); else x, y := PhiX_and_mX(XGroup(A`M,p), XGroupV(A,p)); end if; A`xdata[i] := [x,y]; end if; return A`xdata[i][1]; end intrinsic; intrinsic XModularDegree(A::ModTupFld) -> RngIntElt {Compute the p-adic modular degree of A, at largest prime which exactly divides the level.} for x in Reverse(Factorization(A`N)) do if x[2] eq 1 then return XModularDegree(A,x[1]); end if; end for; error "N has no square free divisors."; end intrinsic; intrinsic XModularDegree(A::ModTupFld, p::RngIntElt) -> RngIntElt {Compute the p-adic modular degree of A.} if not assigned A`xdata or A`xdata[PrimePos(p,A`N)][1] eq 0 then dummy := PhiX(A,p); end if; return A`xdata[PrimePos(p,A`N)][2]; end intrinsic; /////////////// END ACCESS TO QUATERNION ALGEBRAS CODE //////////// /////////////////////////////////////////////////////////////////// // Access to the Mestre's module. // /////////////////////////////////////////////////////////////////// intrinsic PhiX_and_mX(X::ModTupRng, V::ModTupRng) -> RngIntElt,RngIntElt {Computes the quantity PhiX associated to V.} n := Degree(X); m := Dimension(V); weights:= MonodromyWeights(X); W := MatrixAlgebra(Integers(),n)!0; for i in [1..n] do W[i,i] := weights[i]; end for; L := LatticeWithGram(W); B := [ L!v : v in Basis(X)]; // it is very important to move to L! C := [ L!v : v in Basis(V) ]; S := Eltseq( LLL( RMatrixSpace(Integers(),n-1,m)! &cat[ [ InnerProduct(u,v) : v in C ] : u in B ]) ); PhiX := AbsoluteValue( Determinant( MatrixAlgebra(Integers(),m)! [ S[i] : i in [1..m^2] ] ) ); disc := Determinant( MatrixAlgebra(Integers(),m)! [ InnerProduct(u,v) : u in C, v in C ] ); mX := AbsoluteValue(disc/PhiX); return Integers()!PhiX, Integers()!mX; end intrinsic; intrinsic MestreGroup(M::ModTupFld) -> ModTupRng {Returns the Mestre module of M. The level must be prime and the weight must be two.} if Weight(M) ne 2 then error "MestreGroup: weight must be 2."; end if; assert IsPrime(Conductor(M)); if not assigned M`mestre then M`mestre := Mestre(Conductor(M)); end if; return M`mestre; end intrinsic; intrinsic MestreGroupV(A::ModTupFld) -> ModTupRng {Return the factor of the character group corresponding (maybe) to the p-adic rigid analytic optimal quotient associated to A. The representation is via the Mestre construction.} if ( not IsCuspidal(A) ) or ( not Weight(A) eq 2) or ( not IsNew(A) ) then error "MestreGroupV: weight must be 2 and A must be cuspidal and new."; end if; assert IsPrime(Conductor(A)); if not assigned A`mestre then Z := MestreGroup(ModSymParent(A)); A`mestre := RSpace(Integers(),Degree(Z)); d := DimensionAbelianVariety(A); p := 2; while Dimension(A`mestre) gt d do T := TpD(Z, p); // compute p-th Hecke operator. f := PolynomialRing(Integers())!DecompCharpolyTp(A, p); A`mestre := A`mestre meet Kernel(Evaluate(f,T)); p := NextPrime(p); if p gt 7 and HeckeVerbose then "MestreGroupV: Warning -- space not breaking up.";; end if; end while; end if; return A`mestre; end intrinsic; intrinsic MestreEigenvector(A::ModTupFld) -> . {} return MestreEigenvector(MestreGroup(ModSymParent(A)), MestreGroupV(A)); end intrinsic; intrinsic MestreEigenvector(X::ModTupRng, V::ModTupRng) -> . {} p := 2; T := Restrict(TpD(X,p),V); f := ModularCharpoly(T); while not IsIrreducible(f) do p := NextPrime(p); if p gt 7 then error "Couldn't verify that V is irreducible."; end if; T +:= Random(1,10)*Restrict(TpD(X,p),V); f := ModularCharpoly(T); end while; e := Eigenvector(MatrixAlgebra(Rationals(),Ncols(T))!T,f); F := Parent(e[1]); W := VectorSpace(F,Degree(V)); B := [W!b : b in Basis(V)]; return &+[e[i]*B[i] : i in [1..#B]]; end intrinsic; intrinsic MestrePowerSum(X::ModTupRng, e::ModTupFldElt, n::RngIntElt) -> . {This function returns sum w_i (e_i)^n, a quantity which, by Gross-Zagier, is relevant to the computation of the order of vanishing of L(A,s) at 1. Currently Conductor(A) must be prime.} w := MonodromyWeights(X); assert Degree(e) eq #w; return &+[e[i]^n*w[i] : i in [1..#w]]; end intrinsic; intrinsic MestrePowerSum(A::ModTupFld, n::RngIntElt) -> . {Let v be a choice of eigenvector, corresponding to A, on the module of supersingular points. This function returns sum w_i (v_i)^n, a quantity which, by Gross-Zagier, is relevant to the computation of the order of vanishing of L(A,s) at 1. Currently Conductor(A) must be prime.} if not IsPrime(Conductor(A)) then error "Currently Conductor(A) must be prime."; end if; return MestrePowerSum(MestreGroup(ModSymParent(A)), MestreEigenvector(A),n); end intrinsic; /**************************************************************** * FUNCTIONS FOR COMPUTING ORDERS OF COMPONENT GROUPS * * Ref. Stein, "Component groups of optimal quotients of * * Jacobians", Preprint, 1999. * ****************************************************************/ intrinsic ComponentGroup(A::ModTupFld) -> RngIntElt {Compute the order of the group of *geometric* points of the component group at the largest prime that exactly divides the level.} assert IsNew(A); for x in Reverse(Factorization(A`N)) do if x[2] eq 1 then return ComponentGroup(A,x[1]); end if; end for; error "N has no square free divisors."; end intrinsic; intrinsic ComponentGroup(A::ModTupFld, p::RngIntElt) -> RngIntElt {Compute the order of the group of *geometric* points of the component group at p, assuming that p exactly divides the level. If working in the +1 quotient, then only the odd part of the order is returned.} if not IsPrime(p) then error "ComponentGroup: input p =", p, " is not prime!"; end if; if Valuation(Level(A),p) eq 0 then return 1; end if; if Valuation(A`N,p) ne 1 then "ComponentGroup: WARNING, don't know how to compute component group at",p; return 0; // DO NOT KNOW! not know. end if; phi := PhiX(A,p); mH := ModularDegree(A); if Sign(A`M) ne 0 then mH /:= 2^Valuation(mH,2); else mH := Round(Sqrt(mH)); end if; mX := XModularDegree(A,p); ans := phi * mH / mX; if Sign(A`M) ne 0 then ans /:= 2^Valuation(ans,2); end if; return Integers()!ans; end intrinsic; intrinsic TamagawaNumber(A::ModTupFld) -> RngIntElt {Compute the order of the group of Fp rational points of the component group of A at the largest prime which exactly divides the level of A. WARNING: Stein has not yet nailed down the power of 2 when Wp=+1!} for x in Reverse(Factorization(A`N)) do if x[2] eq 1 then return TamagawaNumber(A,x[1]); end if; end for; error "N has no square free divisors."; end intrinsic; intrinsic TamagawaNumber(A::ModTupFld, p::RngIntElt) -> RngIntElt {} if Valuation(Level(A),p) eq 0 then return 1; end if; if Valuation(A`N,p) ne 1 then return 0; // i.e., do not know. // error "TamagawaNumber: p must exactly divide the level."; end if; cp := ComponentGroup(A,p); w := WqOn(A,p); if w[1,1] eq -1 then return cp; else return 1; end if; end intrinsic; intrinsic Sha(A::ModTupFld) -> RngIntElt {Try to compute the odd part of the order of Sha. Assuming BSD we obtain a number which is either zero or divides the true order of Sha, and only misses primes where the representations are reducible. The following caveats apply: (a) We assume the BSD formula. (b) We use the upper bound on torsion coming from Hecke operators, so the result of this function might be too small. (C) At primes whose square divide the level we do not know how to compute cp; however such cp can only be divisible by p and primes dividing the order of the torsin group. We simply make our conjectural value of Sha coprime to such primes when there is a p whose square divides the level.} if not IsCuspidal(A) then error "There is no Sha associated to noncuspidal A."; end if; fac := Factorization(A`N); sha := CuspOrder(A)^2; // lower bound on torsion. sha *:= LRatio(A); // L(1)/Omega. if sha eq 0 then return 0 ; end if; S := [f[1] : f in fac | f[2] eq 1]; sha /:= &*[Integers()|TamagawaNumber(A,p) : p in S]; T := 2; if #fac ne #S then // not semistable at some prime. if fac[1][1] eq 2 and fac[1][2] gt 1 then T := sha; // add. reduction at 2, all bets are off!! end if; T *:= TorsionBound(A) * Conductor(A); end if; for f in Factorization(T) do sha /:= f[1]^Valuation(sha,f[1]); end for; return Numerator(sha); // might have extra denom left from bad // torsion estimate. end intrinsic; /************************************************************* * * * EIGENVALUES: COMPUTATION OF HECKE EIGENVALUES * * * *************************************************************/ intrinsic Eigenvector(T::AlgMatElt) -> ModTupFldElt {} f := ModularCharpoly(T); assert IsIrreducible(f); return Eigenvector(T,f); end intrinsic; intrinsic Eigenvector(T::AlgMatElt, f::RngUPolElt) -> ModTupFldElt {Let T be an nxn matrix over K with irreducible characteristic polynomial f. This function returns an eigenvector for T over the extension field K[x]/(f(x)).} /****************************************************************** I decided to implement this using a quotient of a polynomial ring instead of the Magma built in number field functions because basic arithmetic in polynomial extensions surprisingly seem to be almost twice as fast... ******************************************************************/ n := Degree(f); K := Parent(T[1,1]); Kn := VectorSpace(K, n); if n eq 1 then return Kn![1]; end if; R := PolynomialRing(K); L := quo; b := 1/a; c := [-b*Coefficient(f,0)]; for i in [1..Degree(f)-1] do Append(~c, (c[i] - Coefficient(f,i))*b); end for; Ln := RSpace(L,n); v := Ln!0; T := RMatrixSpace(L,n,n)!T; repeat v[Random(1,n)] +:= 1; w := c[1]*v; vv := v; for i in [2..#c] do vv := vv*T; w +:= c[i]*vv; end for; until w ne 0; return w; end intrinsic; ////////////////////////////////////////////////////////////// intrinsic ZAlgDisc(gens::SeqEnum) -> FldRatElt {The elements of gens should all lie in a field Q[x]/(f(x)). It is assumed that they generate a subring O of Q[x]/(f(x)), as a *Z-module*! This function returns the discriminant of this Z-module.} /* We use that disc(O) = Vol(L)^2 * disc(f). where Vol(L) is the volume of the lattice in Euclidean space cut out by the vectors (b_0,..,b_n) corresponding to generators b_0 + b_1*x + ... + b_n*x^n of O_f listed in gens. */ f := Modulus(Parent(gens[1])); n := Degree(f); V := VectorSpace(Rationals(),n); return Discriminant(f) * Volume([ V![Coefficient(g,i): i in [0..n-1]] : g in gens])^2; end intrinsic; // added Feb 09, 2000 intrinsic Discriminant(A::ModTupFld) -> RngIntElt {If A is newform modular factor, compute the discriminant of the ring generated by the Fourier coefficients of one of the q-expansions of A.} if GetCategory(A) ne ModSymFac or not IsNew(A) or not IsCuspidal(A) then error "A must be a new modular symbols factor."; end if; if not assigned A`disc then if DimensionAbelianVariety(A) eq 1 then A`disc := 1; else f := qEigenform(A,HeckeBound(ModSymParent(A))); A`disc := ZAlgDisc([Coefficient(f,i) : i in [1..Degree(f)]]); end if; end if; return A`disc; end intrinsic; intrinsic DiscriminantKf(A::ModTupFld) -> RngIntElt {If A is newform modular factor, compute the discriminant of the ring generated by the Fourier coefficients of one of the q-expansions of A.} a :=Coefficient(qEigenform(A,3),1); return Discriminant(NumberField(Modulus(Parent(a)))); end intrinsic; intrinsic SignedSubspace(A::ModTupFld, sign::RngIntElt) -> ModTupFld {Compute the motive B made from the *=sign subspace of the modular motive A.} assert -1 le sign and sign le 1; if not assigned A`signspace then S := Transpose(StarInvolution(A`M)); Aplus := Kernel(S-1) meet A; Aminus := Kernel(S+1) meet A; Aplus := Factor(A`M,Aplus); Aminus := Factor(A`M,Aminus); A`signspace := [* Aplus, Aminus *]; end if; if sign eq 1 then if IsPlusQuotient(A`M) then return A; end if; if IsMinusQuotient(A`M) then error "No +1 subspace since we are working in the -1 quotient."; end if; return A`signspace[1]; elif sign eq -1 then if IsPlusQuotient(A`M) then error "No -1 subspace since we are working in the +1 quotient."; end if; if IsMinusQuotient(A`M) then return A; end if; return A`signspace[2]; else error "SignedSubspace: invalid sign."; end if; end intrinsic; intrinsic Eigenvector(A::ModTupFld, sign::RngIntElt) -> ModTupFldElt {Compute eigenvector for sign subspace of A.} if IsPlusQuotient(A`M) then if sign eq -1 then "Eigenvector: when working in +1 subspace there are "; error "no -1 eigenvectors, so the computation can not continue."; end if; return Eigenvector(A); end if ; if IsMinusQuotient(A`M) then if sign eq +1 then "Eigenvector: when working in -1 subspace there are "; error "no +1 eigenvectors, so the computation can not continue."; end if; return Eigenvector(A); end if ; if sign eq +1 then if not assigned A`eigenplus then A`eigenplus := Eigenvector(SignedSubspace(A,sign)); end if; return A`eigenplus; end if; if sign eq -1 then if not assigned A`eigenminus then A`eigenminus := Eigenvector(SignedSubspace(A,sign)); end if; return A`eigenminus; end if; end intrinsic; intrinsic Eigenvector(A::ModTupFld) -> ModTupFldElt {Returns an eigenvector of the Hecke algebra on A over a polynomial extension of the ground field. A must be new.} if not assigned A`eigen then // 1. Find a linear combination of Hecke operators whose // charpoly on A is irreducible. p := SmallestPrimeNondivisor(A`M`N); T := DecompFastTp(A,p); f := ModularCharpoly(T); while not IsIrreducible(f) do p := NextPrime(p); p := SmallestPrimeNondivisor(A`M`N,p); T +:= DecompFastTp(A,p); f := ModularCharpoly(T); end while; e := Eigenvector(T,f); F := Parent(e[1]); V := VectorSpace(F,Degree(A)); B := [V!b : b in Basis(A)]; A`eigen := &+[e[i]*B[i] : i in [1..#B]]; end if; return A`eigen; end intrinsic; intrinsic HeckeImages(i::RngIntElt, n::RngIntElt, M::ModTupFld) -> SeqEnum {Returns the images of the ith standard basis vector under the Hecke operators Tp for p<=n.} assert 1 le i and i le Dimension(M); if not assigned M`T_images then M`T_images := [[M|] : i in [1..Dimension(M)]]; end if; if #M`T_images[i] lt PrimePos(n) then // generate more images... p := NthPrime(#M`T_images[i]+1); s := SparseRepresentation(M.i); while p le n do Append(~((M`T_images)[i]),TnSparse(M, p, s)); p := NextPrime(p); end while; end if; return M`T_images[i]; end intrinsic; intrinsic Dot(v::ModTupFldElt, w::ModTupFldElt) -> . {Compute the dot product of v and w.} assert Degree(v) eq Degree(w); ans := &+[v[i]*w[i] : i in [1..Degree(v)]]; return ans; end intrinsic; intrinsic DimensionAbelianVariety(A::ModTupFld) -> RngIntElt {Returns the dimension of the abelian variety attached to A.} return Sign(A`M) ne 0 select Dimension(A) else Integers()!(Dimension(A)/2); end intrinsic; procedure RecursivelyCompute_an(~f, start, eps, k, F) Q := PowerSeriesRing(F); for n in [start..Degree(f)] do if not IsPrime(n) then fac := Factorization(n); if #fac eq 1 then // a_{p^r} := a_p * a_{p^{r-1}} - eps(p)p^{k-1} a_{p^{r-2}}. p := fac[1][1]; r := fac[1][2]; coef := Coefficient(f,p) * Coefficient(f,p^(r-1)) - Evaluate(eps,p)*p^(k-1) * Coefficient(f,p^(r-2)); else // a_m*a_r := a_{mr} and we know all a_i for i RngSerPowElt {Attempt to compute the q-expansion of one of the Galois conjugate modular forms associated to A. This function uses theta series, so if successful it can be used to compute a very large number of Fourier coefficients quickly. If there is no prime which exactly divides Level(A), Weight(A) gt 2, the character of A is nontrivial, then the function qEigenform(A,prec) is called.} requirege prec, 1; assert IsNew(A) and IsCuspidal(A); if (Weight(A) gt 2) or not IsTrivial(Character(A)) then "Warning: theta series not be used."; return qEigenform(A,prec); end if; fac := Factorization(Level(A)); if not exists(t){x[1] : x in Factorization(Level(A)) | x[2] eq 1} then "Warning: theta series not be used."; return qEigenform(A,prec); end if; prec := NextPrime(prec-1); if not assigned A`qexp or Degree(A`qexp) lt prec then if HeckeSuperverbose then "qEigenformTheta: Computing q-expansion."; end if; XJ := XGroup(A`M); V := XGroupV(A); dummy := TXn(XJ,prec); // precompute prec Hecke operators. // Find an eigenvector in V. p := SmallestPrimeNondivisor(A`M`N); T := Restrict(TXn(XJ,p),V); f := ModularCharpoly(T); while not IsIrreducible(f) do p := NextPrime(p); p := SmallestPrimeNondivisor(A`M`N,p); T +:= Restrict(TXn(XJ,p),V); f := ModularCharpoly(T); end while; Mn := MatrixAlgebra(Rationals(),Ncols(T)); eig := Eigenvector(Mn!T); F := Parent(eig[1]); S := VectorSpace(F, Degree(V)); B := [S!b : b in Basis(V)]; e := &+[B[i]*eig[i] : i in [1..#B]]; i := 1; while e[i] eq 0 do i +:= 1; end while; MnF := MatrixAlgebra(F,Degree(V)); Q := PowerSeriesRing(F); A`qexp := q; for n in [p : p in [2 .. prec] | IsPrime(p)] do A`qexp +:= q^n * (e*(MnF!TXn(XJ,n)))[i]/e[i]; end for; RecursivelyCompute_an(~A`qexp, 4, Character(A), Weight(A), F); end if; return A`qexp; end intrinsic; intrinsic qEigenform(A::ModTupFld) -> RngSerPowElt {Returns the at least 7 terms of the q-expansion of A.} return qEigenform(A,7); end intrinsic; intrinsic qEigenform(A::ModTupFld, prec::RngIntElt) -> RngSerPowElt {Returns the q-expansion of one of the eigenforms associated to A, computed to precision at least prec.} if not IsCuspidal(A) then error "I haven't yet programmed q-expansions of Eisenstein series."; end if; prec := NextPrime(prec-1); if not assigned A`qexp or Degree(A`qexp) lt prec then if HeckeSuperverbose then "qEigenform: Computing eigenvector."; end if; if not IsNew(A) then A`qexp := qEigenform(OldFactor(A),prec); else // compute q-expansion of newform. if IsCuspidal(A) then eig := Eigenvector(A,IsMinusQuotient(A) select -1 else +1); else eig := Eigenvector(A); end if; // find i s.t. =/= 0. i := 1; while eig[i] eq 0 do i +:= 1; end while; if HeckeSuperverbose then "qEigenform: Computing Hecke images."; end if; Tpei := HeckeImages(i, prec, A`M); Q := PowerSeriesRing(Parent(eig[1])); F := A`M`F; // Will change, I hope... if Type(eig[1]) eq RngUPolResElt then // The Magma divide in a relative extension is // SUPER SUPER AMAZINGLY STUPIDLY **SLOW** gcd, c,d := Xgcd(PolynomialRing(F)!eig[i], Modulus(Parent(eig[i]))); assert gcd eq 1; one_over := Parent(eig[i])!c; else one_over := 1/eig[i]; end if; if not assigned A`qexp then // create f the first time. startrecur := 4; // startrecur is only used below for a_n f := 0; // Record the quotients a_p = / p := 2; for j in [1..#Tpei] do coef := 0; for m in [1..Degree(eig)] do if Tpei[j][m] ne 0 then coef +:= Tpei[j][m] * eig[m]; end if; end for; f +:= coef * q^p; if HeckeSuperverbose then printf "a_%o\t----->%10o\n",p,coef; end if; p := NextPrime(p); end for; f *:= one_over ; // divide by ; f +:= q; if IsEisenstein(A) then f +:= 1; end if; A`qexp := f; else // add on more terms. d := Degree(A`qexp); p := d; while not IsPrime(p) do p -:= 1; end while; // remove terms of f of degree at least p. for i in [p..Degree(A`qexp)] do A`qexp -:= Coefficient(A`qexp,i)*q^i; end for; startrecur := Max(p+1,4); // used below for a_n for j in [PrimePos(p)..#Tpei] do coef := 0; for m in [1..Degree(eig)] do if Tpei[j][m] ne 0 then coef +:= Tpei[j][m] * eig[m]; end if; end for; coef *:= one_over; A`qexp +:= coef * q^p; p := NextPrime(p); end for; end if; // Obtain a_n by the recurence. RecursivelyCompute_an(~A`qexp,startrecur, A`M`eps, A`M`k, Parent(eig[1])); end if; end if; return A`qexp; end intrinsic; function Saturate(B) F := Parent(B[1][1]); M := Transpose(RMatrixSpace(F,#B, #B[1])! &cat B); C := [Eltseq(b) : b in Basis(Kernel(M))]; if #C eq 0 then N := RMatrixSpace(Rationals(),#B[1],#B[1])!0; else N := Transpose(RMatrixSpace(F, #C, #B[1])! &cat C); end if; return IntegerKernel(N); end function; intrinsic SaturatePolySeq(P::SeqEnum) -> SeqEnum {Returns generators for the saturation of the Z-module generated by the sequence of q-expansions. Uses the smallest degree as precision.} if #P eq 0 then return P; end if; R := Parent(P[1]); m := Max([Degree(f) : f in P]); X := [[Integers()!Coefficient(f,i) : i in [0..m]] : f in P]; S := Saturate(X); return [&+[Integers()!(v[i])*q^(i-1) : i in [1..m+1]] : v in Basis(S)]; end intrinsic; intrinsic IndexInSaturation(P::SeqEnum) -> RngIntElt, SeqEnum {Returns the index and structure of the Z-module generated by the sequence of q-expansions in its saturation. Uses the smallest degree as precision.} if #P eq 0 then return 1; end if; R := Parent(P[1]); m := Min([Degree(f) : f in P]); X := [[Coefficient(f,i) : i in [0..m]] : f in P]; D := RMatrixSpace(Integers(),#P,(m+1))!(&cat X); S := SmithForm(D); ind := [S[i,i] : i in [1..Min(Nrows(S),Ncols(S))] | S[i,i] gt 1]; if #ind eq 0 then return 1, ind; end if; return &*ind, ind; end intrinsic; intrinsic EigenformInTermsOfIntegralBasis(A::ModTupFld) -> ModTupFld {Returns the linear combination of qIntegralBasis(A) which gives Eigenform(A).} /***************************************************** * 1) Find pivot columns for the integral basis * * 2) Invert * * 3) Obtain eigenform in terms of integral basis * *****************************************************/ if not assigned A`eigenint then B := qIntegralBasis(A); f := qEigenform(A); F := Parent(Coefficient(f,1)); d := Degree(f); V := VectorSpace(F,d); fvec := V![Coefficient(f,n) : n in [1..d]]; Bvec := [V![Coefficient(g,n) : n in [1..d]] : g in B]; W := VectorSpaceWithBasis(Bvec); A`eigenint := Eltseq(Coordinates(W,fvec)); end if; return A`eigenint; end intrinsic; intrinsic qIntegralBasis(D::SeqEnum) -> SeqEnum {Return an integral basis for the sum of the spaces in the decomposition D.} return qIntegralBasis(D,7); end intrinsic; intrinsic qIntegralBasis(D::SeqEnum, prec::RngIntElt) -> SeqEnum {Return an integral basis for the sum of the cuspidal spaces in the decomposition D.} return SaturatePolySeq(&cat[qIntegralBasis(B,prec) : B in D | IsCuspidal(B)]); end intrinsic; intrinsic qIntegralBasis(A::ModTupFld) -> SeqEnum {} return qIntegralBasis(A,7); end intrinsic; intrinsic qIntegralBasis(A::ModTupFld, prec::RngIntElt) -> SeqEnum {Compute an integral basis for the space spanned by the Galois conjugates associated to A. If A is a space of modular symbols, returns an integral basis for the corresponding space of cusp forms. The base field must be the rationals.} if not assigned A`qintbasis or Degree(A`qintbasis[1]) lt prec then case GetCategory(A): when ModSymFac: if IsNew(A) then f := qEigenform(A,prec); Q := Parent(Coefficient(f,1)); if Type(Q) eq FldRat then A`qintbasis := [f]; else g := Modulus(Q); n := Degree(g); R := PreimageRing(Q); B := [[Coefficient(R!Coefficient(f,i),j) : i in [1..Degree(f)]] : j in [0..n-1]]; if Type(Parent(B[1][1])) ne FldRat then "qIntegralBasis: choice of Z[zeta_n]-lattice may be not well-defined."; C := B; F := A`M`F; else C := Basis(Saturate(B)); F := Integers(); end if; R := PowerSeriesRing(F); A`qintbasis := [ &+[(R!g[n])*q^n : n in [1..Degree(f)]] : g in C]; end if; else Q := qIntegralBasis(OldFactor(A),prec); NN := Level(OldFactor(A)); R := Parent(Q[1]); A`qintbasis := SaturatePolySeq( [Evaluate(f,q^d) : d in Divisors(Level(A) div NN), f in Q]); end if; when ModSym: A`qintbasis := qIntegralBasis(FullDecomposition(A),prec); end case; end if; return A`qintbasis; end intrinsic; /************************************************************* * * * ARITHMETIC OF FACTORS * * * *************************************************************/ intrinsic AddFactors(A::ModTupFld, B::ModTupFld) -> ModTupFld {Compute the sum of two modular symbols factors.} F := Factor(ModSymParent(A), A + B); if IsCuspidal(A) and IsCuspidal(B) then F`iscusp := true; else F`iscusp := false; end if; return F; end intrinsic; intrinsic IntersectFactors(A::ModTupFld, B::ModTupFld) -> ModTupFld {Compute the intersection of two modular symbols factors. If the intersection is zero-dimensional use the function "Intersection" to compute the structure of the corresponding intersection of abelian varieties.} assert ModSymParent(A) eq ModSymParent(B); F := Factor(ModSymParent(A), A meet B); if IsCuspidal(A) and IsCuspidal(B) then F`cusp := true; else F`cusp := false; end if; return F; end intrinsic; /**********END FACTOR ARITHMETIC *****************************/ /************************************************************** CUTTING OUT A SUBSPACE OF THE DUAL USING SUCCESSIVE FASTTP **************************************************************/ intrinsic DualSubspace(M::ModTupFld, polyprimes::SeqEnum) -> ModTupFld {Given a list polyprimes of pairs computes the following subspace of the dual of M: intersect ker(f(Transpose(Tp))). } p := polyprimes[1][1]; f := polyprimes[1][2]; t := Cputime(); if HeckeVerbose then "Computing Ker(Transpose(f(T_", p,")))"; end if; T := Evaluate(f, Transpose(HeckeOperator(M,p))); W := Kernel(T); if HeckeVerbose then Cputime(t), " seconds."; end if; return DualSubspace(M, W, Remove(polyprimes,1)); end intrinsic; intrinsic DualSubspace(M::ModTupFld, W::ModTupFld, polyprimes::SeqEnum ) -> ModTupFld { Given a sequence polyprimes [,...] of pairs a prime and a polynomial, this function computes the following subspace of the subspace W of the dual of M: intersect ker(f(Transpose(Tp))). } for i in [1..#polyprimes] do if HeckeVerbose then "Current dimension = ", Dimension(W); end if; p := polyprimes[i][1]; f := polyprimes[i][2]; t := Cputime(); if HeckeVerbose then "Computing Ker(Transpose(f(T_", p,")))"; end if; T := FastTn(M, W,p); fT := Evaluate(f,T); K := KernelOn(fT,W); if HeckeVerbose then Cputime(t), " seconds."; end if; W := W meet K; end for; if HeckeVerbose then "Current dimension = ", Dimension(W); end if; return W; end intrinsic; /************************************************************* * * * COMPUTATION OF THE RATIONAL PERIOD MAPPING * * r_A : M_k(Q) ----> M_k(Q)/Ker(Phi). * * * *************************************************************/ intrinsic RationalPeriodMapping(A::ModTupFld) -> Map {Returns the map r_A : M_k(Q) ----> M_k(Q)/Ker(Phi).} if not assigned A`RatPeriodMap then M := ModSymParent(A); V := VectorSpace(BaseField(M), Dimension(A)); A`RatPeriodMap := hom < M -> V | [V![InnerProduct(m,a) : a in Basis(A)] : m in Basis(M)]>; end if; return A`RatPeriodMap; end intrinsic; intrinsic SignedRationalPeriodMapping(A::ModTupFld, sign::RngIntElt) -> Map {Returns the map r_A : M_k(Q) ----> M_k(Q)/(Ker(Phi)+(sign quotient)).} if not assigned A`RatPeriodMapSign then M := ModSymParent(A); Bplus := SignedSubspace(A,1); Bminus := SignedSubspace(A,-1); V := VectorSpace(BaseField(M), Dimension(Bplus)); A`RatPeriodMapSign := [* hom < M -> V | [V![InnerProduct(m,b) : b in Basis(Bplus)] : m in Basis(M)]>, hom < M -> V | [V![InnerProduct(m,b) : b in Basis(Bminus)] : m in Basis(M)]> *]; end if; if sign eq 1 then return A`RatPeriodMapSign[1]; else return A`RatPeriodMapSign[2]; end if; end intrinsic; intrinsic RationalPeriodLattice(A::ModTupFld) -> SeqEnum {Returns a basis for the image of S_k(Z) in M_k/Ker(Phi).} if not assigned A`RatPeriodLat then pi := RationalPeriodMapping(A); // Generators for rational period lattice. B := [pi(v) : v in Basis(SkZ(ModSymParent(A)))]; A`RatPeriodLat := VectorSpaceZBasis(B); end if; return A`RatPeriodLat; end intrinsic; intrinsic RationalPeriodConjugation(A::ModTupFld) -> AlgMatElt {Returns matrix of "conjugation" with respect to the basis RationalPeriodLatticeBasis(A).} if not assigned A`RatPeriodConj then pi := RationalPeriodMapping(A); Z := RationalPeriodLattice(A); B := Basis(Z); star := StarInvolution(ModSymParent(A)); C := [Coordinates(Z, pi((b@@pi)*star)) : b in B]; C := &cat[Eltseq(x) : x in C]; A`RatPeriodConj := MatrixAlgebra(BaseField(A), #B)!C; end if; return A`RatPeriodConj; end intrinsic; intrinsic ImaginaryTamagawa(A::ModTupFld) -> RngIntElt {Computes the number of real components of the abelian variety associated to A.} if HeckeBugWatch then "(Warning: Imaginary comp group formula not yet proved.)"; end if; if not assigned A`cinfimag then M:=ModSymParent(A); if Sign(M) ne 0 then "WARNING: Working in + or - quotient, so can not compute imaginary components."; return 1; end if; C := RationalPeriodConjugation(A); Cmod2 := MatrixAlgebra(GF(2),Ncols(C)) ! C; A`cinfimag := 2^(Dimension(Kernel(Cmod2+1))-DimensionAbelianVariety(A)); end if; return A`cinfimag; end intrinsic; intrinsic RealTamagawa(A::ModTupFld) -> RngIntElt {Computes the number of real components of the abelian variety associated to A.} if not assigned A`cinf then M:=ModSymParent(A); if Sign(M) ne 0 then "WARNING: Working in + or - quotient, so can not compute real components."; return 1; end if; C := RationalPeriodConjugation(A); Cmod2 := MatrixAlgebra(GF(2),Ncols(C)) ! C; A`cinf := 2^(Dimension(Kernel(Cmod2-1))-DimensionAbelianVariety(A)); end if; return A`cinf; end intrinsic; /***************** END RATIONAL PERIOD MAPPING *********************/ /************************************************************** * * * NUMERICAL COMPUTATION OF PERIODS * * * **************************************************************/ // Constants EULER := EulerGamma(RealField()); PI := Pi(RealField()); function DefaultPrec(A) prec := 50+Integers()!(Round(Conductor(A) div 10)); if assigned A`qexp then prec := Max(prec,Degree(A`qexp)); end if; if HeckeVerbose then printf "(Using at least %o terms of q-expansions.)\n",prec; end if; return prec; end function; ///////////////////////////////////////////////////////////// // Compute Int_{alp}^{oo} z^m*q^n dz, using // // an explicit formula derived via integration by parts. // // (C = complex numbers) // ///////////////////////////////////////////////////////////// function Int1(alp, m, n, C) c := 2*Pi(C)*C.1*n; return -2*Pi(C)*C.1 * Exp(c*alp) * &+[(-1)^s *(alp^(m-s)) / c^(s+1) * (&*[Integers()|i : i in [m-s+1..m]]) : s in [0..m]]; end function; ///////////////////////////////////////////////////////////// // Compute Int_{alp}^{oo} P(z)*f(q) dz, using // // an explicit formula derived via integration by parts. // ///////////////////////////////////////////////////////////// function Int2(f, alp, m) C := ComplexField(); return &+[Int1(alp,m,n,C)*(C!Coefficient(f,n)) : n in [1..Degree(f)]]; end function; intrinsic PeriodIntegral(f::RngSerPowElt, alpha::FldPrElt ) -> FldPrElt {Computes for alpha any point in the upper half plane.} return PeriodIntegral(f,PolynomialRing(Rationals())!1,alpha); end intrinsic; intrinsic PeriodIntegral(f::RngSerPowElt, P::RngUPolElt, alpha::FldPrElt ) -> FldPrElt {Computes for alpha any point in the upper half plane.} C := ComplexField(); return &+[Int2(f,alpha,m)* (C!Coefficient(P,m)) : m in [0..Degree(P)] | Coefficient(P,m) ne 0]; end intrinsic; //////////////////////////////////////////////////////////////// // In the following, "fast" refers to using the Fricke // // involution trick to speed convergence. This only works // // in some cases, so we have also implemented a standard // // slow method which should work in all cases. // //////////////////////////////////////////////////////////////// intrinsic FastPeriodIntegral(A::ModTupFld, f::RngSerPowElt, g::SeqEnum) -> FldPrElt {Compute .} R := PolynomialRing(Rationals(),2); return FastPeriodIntegral(A, f, ); end intrinsic; intrinsic FastPeriodIntegral(A::ModTupFld, f::RngSerPowElt, Pg::Tup) -> FldPrElt {} assert IsNew(A); P, g := Explode(Pg); if g[3] eq 0 then // g(oo)=oo. return 0; end if; if g[4] lt 0 then g[1] *:= -1; g[2] *:= -1; g[3] *:= -1; g[4] *:= -1; end if; R := Parent(P); phi := hom R | g[1]*R.1+g[2]*R.2, g[3]*R.1+g[4]*R.2>; giP := phi(P) / (g[1]*g[4]-g[2]*g[3]); assert giP eq P; C := ComplexField(); N := Level(A); k := Weight(A); e := SignWN(A); rootN:= Sqrt(Level(A)); PI := Pi(ComplexField()); a := g[1]; b := g[2]; c := g[3] div N; d := g[4]; if k eq 2 then // the formula takes a simpler form when k=2. cn := [ (e-1)*Exp(-2*PI*n/rootN) + Exp(-2*PI*n/(d*rootN))*(Exp(2*PI*i*n*b/d)-e*Exp(2*PI*i*n*c/d)) : n in [1..Degree(f)] ]; return &+[(Coefficient(f,n)/n) * cn[n] : n in [1..Degree(f)]]; else /* The following formula is in my (William Stein) Ph.D. thesis. It generalizes Cremona's Prop 2.10.3 to the case of higher weight modular symbols. */ W := hom R | R.2, -N*R.1>; WP := W(P)/N^(Integers()!(k/2-1)); S := PolynomialRing(Rationals()); ev := hom S | z, 1>; P := ev(P); WP := ev(WP); ans := PeriodIntegral(f, P-e*WP, i/rootN) +PeriodIntegral(f, e*WP, c/d+i/(d*rootN)) +PeriodIntegral(f, -P , b/d+i/(d*rootN)); return ans; end if; end intrinsic; intrinsic SlowPeriodIntegral(A::ModTupFld, f::RngSerPowElt, Pg::Tup) -> FldPrElt {Given a homogeneous polynomial P(X,Y) of degree k-2, a 2x2 matrix g=[a,b,c,d], and a q-expansion f of a weight k modular form, this function returns the period = Int_\{oo,g(oo)\} P(z,1) f(z) dz.} // Resort to VERY and SLOWLY CONVERGING METHOD. /* We have, for any alp in the upper half plane, = = = = . The integrals converge best when Im(alp) and Im(g(alp)) are both large. We choose alp=i/c. ************/ if not IsTrivial(Character(A)) then error "NOT YET PROGRAMMED FOR NONTRIVIAL CHARACTER."; end if; P, g := Explode(Pg); if g[3] eq 0 then // g(oo)=oo. return 0; end if; R := Parent(P); phi := hom R | g[1]*R.1+g[2]*R.2, g[3]*R.1+g[4]*R.2>; giP := phi(P) / (g[1]*g[4]-g[2]*g[3]); C := ComplexField(); if g[3] lt 0 then // g(oo) = (-g)(oo), since g acts through PSL. g[1] *:= -1; g[2] *:= -1; g[3] *:= -1; g[4] *:= -1; end if; alp := (-g[4]+i)/g[3]; galp := (g[1]*alp + g[2])/(g[3]*alp+g[4]); S := PolynomialRing(Rationals()); ev := hom S | z, 1>; P := ev(P); giP := ev(giP); return PeriodIntegral(f,giP,alp) - PeriodIntegral(f,P,galp); end intrinsic; function NextElement(N,k, P,g, fast) R := Parent(P); if fast then // IT IS ABSOLUTELY ESSENTIAL THAT g[4] be VERY SMALL!! if g[2] lt 2*N then g[2] +:= 1; else g[2] := 1; repeat g[4] +:= 1; until Gcd(g[4],N) eq 1; end if; while Gcd(g[4],g[2]) ne 1 do g[2] +:= 1; end while; d, g[1], g[3] := Xgcd(g[4],-g[2]*N); g[3] *:= N; P :=( (1-g[1]*g[4])*X^2 - g[2]*(g[4]-g[1])*X*Y + g[2]^2*Y^2)^((k-2) div 2) ; else if P ne Y^(k-2) then P := R!(P/X) * Y; else P := X^(k-2); // only change g after trying all P's. if g[1] lt g[3] then g[1] +:= 1; else g[3] +:= N; end if; while Gcd(g[1],g[3]) ne 1 do g[1] +:= Sign(g[1]); end while; d, g[4], g[2] := Xgcd(g[1],-g[3]); end if; end if; return P, g; end function; intrinsic PeriodGenerators(A::ModTupFld, fast::BoolElt) -> SeqEnum, BoolElt {} /////////////////////////////////////////////////////////////////////////// // The period map Phi:Mk(N,eps;Q) --- > C^d factors through // // a finite dimensional Q-vector space quotient // // Mk/Ker(Phi). This function returns a sequence of elements // // of Mk(N,eps,Q) which (1) their images in Mk/Ker(Phi) form // // a basis for Mk/Ker(P), and (2) each symbols is of the // // form P{oo,g(oo)} for some polynomial P and matrix g in Gamma_0(N) // // with the lower left entry of g positive. // // The second condition is necessary so that we can compute // // the map Phi on these elements. // // Note also: if g(P)=P then the precision of the period // // computation can be greatly enhanced. // /////////////////////////////////////////////////////////////////////////// assert IsCuspidal(A) and IsNew(A); if not assigned A`PeriodGens then M := ModSymParent(A); N := Level(A); k := Weight(A); pi := RationalPeriodMapping(A); // M_k --> Mk/Ker(Phi) V := Image(pi); S := sub; G := []; X := []; RR := PolynomialRing(BaseField(A),2); P := a^(k-2); g := [1,N,0,1]; star := StarInvolution(M); cnt := 0; while Dimension(S) lt Dimension(V) do P,g := NextElement(N,k,P,g,fast); x := ConvFromModularSymbol(M, ); // image of x+*(x) and x-*(x) plus := pi (x + x*star) ; minus := pi (x - x*star) ; SS := sub; if Dimension(SS) eq Dimension(S) + 2 then Append(~X, plus); Append(~X, minus); Append(~G, ); S := SS; end if; cnt +:= 1; if cnt gt 500+Conductor(A)^(7/5) then "WARNING: Not finding enough symbols!"; if fast then "Switching to slow method."; return PeriodGenerators(A, false); else "Could not find enough invariant symbols."; error "Need a better algorithm."; end if; end if; end while; A`PeriodGens := G; A`PGfast := fast; end if; return A`PeriodGens, A`PGfast; end intrinsic; intrinsic PeriodMapping(A::ModTupFld) -> SeqEnum {} return PeriodMapping(A, DefaultPrec(A)); end intrinsic; intrinsic ApplyPeriodMapping(Phi::SeqEnum, v::ModTupFldElt) -> ModTupFldElt {Apply the period mapping Phi to v in Mk.} return &+[v[i]*Phi[i] : i in [1..Degree(v)]]; end intrinsic; intrinsic PeriodMapping(A::ModTupFld, n::RngIntElt) -> SeqEnum {Computes the complex period mapping associated to A using n terms of the q-expansions. The period map is a homomorphism Phi:Mk(N,eps;Q)--> C^d, where d=dim A. Furthermore, A(C) := C^d/Phi(Sk(N,eps;Z)).} if not IsNew(A) then error "Not yet programmed for A not new."; end if; if not assigned A`PeriodMap or A`PeriodMapPrecision lt n then k := Weight(A); /**************************************************************** * If fast = true then we can use tricks arising from the * * Fricke involution to greatly speed of the convergence * * of the series. If fast = false we really on a simplistic * * method which gives series that do not converge as quickly. * ****************************************************************/ fast := (k mod 2 eq 0) and IsTrivial(Character(A)); G, fast := PeriodGenerators(A, fast); Q := qIntegralBasis(A, n); C := ComplexField(); Cd := VectorSpace(C,DimensionAbelianVariety(A)); if fast and HeckeSuperverbose then "PeriodMapping: Using WN trick."; end if; if fast then P := [Cd![FastPeriodIntegral(A,Q[i],x) : i in [1..#Q]] : x in G]; else P := [Cd![SlowPeriodIntegral(A,Q[i],x) : i in [1..#Q]] : x in G]; end if; M := ModSymParent(A); star := StarInvolution(M); gens := &cat[[v + v*star, v-v*star] where v is ConvFromModularSymbol(M, ) : x in G]; // images of the gens periods := &cat[[2*Cd![Real(z):z in Eltseq(x)], 2*i*Cd![Imaginary(z):z in Eltseq(x)]] : x in P]; pi := RationalPeriodMapping(A); PG := [ pi(g) : g in gens]; V := VectorSpaceWithBasis(PG); A`PeriodMap := [ &+[periods[i]*coord[i] : i in [1..#coord]] where coord is Coordinates(V, pi(b)) : b in Basis(M) ]; A`PeriodMapPrecision := n; // also compute period lattice: if Type(BaseField(A)!0) eq FldRatElt then S := [ pi(b) : b in Basis(SkZ(M))]; L := Lattice(RMatrixSpace(Rationals(),#S,Dimension(V))!S); A`PeriodLattice := [ &+[periods[i]*coord[i] : i in [1..#coord]] where coord is Coordinates(V, V!b): b in Basis(L)]; else // FAKE lattice "fake lattice..."; error "NOT PROGRAMMED."; S := [ pi(b) : b in Basis(Sk(M))]; // L := RSpaceWithBasis(BaseField(A),Dimension(V))(S); L := 0; //??? A`PeriodLattice := [ &+[periods[i]*coord[i] : i in [1..#coord]] where coord is Coordinates(V, V!b): b in Basis(L)]; end if; end if; return A`PeriodMap; end intrinsic; intrinsic PeriodLattice(A::ModTupFld) -> SeqEnum {} return PeriodLattice(A,DefaultPrec(A)); end intrinsic; intrinsic PeriodLattice(A::ModTupFld, n::RngIntElt) -> SeqEnum {Computes the complex period lattice associated to A using n terms of the q-expansions.} if not assigned A`PeriodLattice or A`PeriodMapPrecision lt n then dummy := PeriodMapping(A,n); end if; return A`PeriodLattice; end intrinsic; intrinsic RealVolume(A::ModTupFld) -> FldPrElt {} return RealVolume(A,DefaultPrec(A)); end intrinsic; intrinsic RealVolume(A::ModTupFld, n::RngIntElt) -> FldPrElt {Computes the volume of A(R), R the RealField. This function returns the volume of the identity component times the number RealTamagawa(A) of real components.} if true or not assigned A`realvolume or A`PeriodMapPrecision lt n then r := RationalPeriodMapping(A); B := Basis(SkZ(ModSymParent(A))); S := Basis(VectorSpaceZBasis([r(v) : v in B])); Conj := RationalPeriodConjugation(A); // * on Image(r). P := IntegerKernelOn(Conj-1, S); // plus one subspace. Z := [v@@r : v in Basis(P)]; Phi := PeriodMapping(A,n); C := &cat[Eltseq(ApplyPeriodMapping(Phi,z)) : z in Z]; R := MatrixAlgebra(Parent(C[1]), #Z)!C; R; A`realvolume := RealTamagawa(A)*Determinant(R); if Real(A`realvolume) lt 0 then A`realvolume *:= -1; end if; end if; return A`realvolume; end intrinsic; intrinsic ImaginaryVolume(A::ModTupFld) -> FldPrElt {} return ImaginaryVolume(A,DefaultPrec(A)); end intrinsic; intrinsic ImaginaryVolume(A::ModTupFld, n::RngIntElt) -> FldPrElt {Computes the volume of A(iR), the pure imaginary points.} if not assigned A`imagvolume or A`PeriodMapPrecision lt n then r := RationalPeriodMapping(A); B := Basis(SkZ(ModSymParent(A))); S := Basis(VectorSpaceZBasis([r(v) : v in B])); Conj := RationalPeriodConjugation(A); // * on Image(r). P := IntegerKernelOn(Conj+1, S); // plus one subspace. Z := [v@@r : v in Basis(P)]; Phi := PeriodMapping(A,n); C := &cat[Eltseq(ApplyPeriodMapping(Phi,z)) : z in Z]; R := MatrixAlgebra(Parent(C[1]), #Z)!C; A`imagvolume := ImaginaryTamagawa(A)*Determinant(R); end if; return A`imagvolume; end intrinsic; function ellinv_NormalizePair(w1, w2) tau := w1/w2; if Imaginary(tau) lt 0 then w1 := -w1; tau := -tau ; end if; w1 := w1 - w2*Round(Real(tau)); tau := w1/w2; for i in [1..300] do // Just to stop infinite loop due to rounding if Abs(tau) ge 1 then break; end if; w3:=-w1; w1:=w2; w2:=w3; tau:=w1/w2; w1:=w1-w2*Round(Real(tau)); tau:=w1/w2; end for; return tau, w1, w2; end function; function ellinv_c4c6(w1, w2) tau := w1/w2; zero := 0; one := 1; two := 2; pi := Pi(RealField()); x := two*pi*Real(tau); y := two*pi*Imaginary(tau); q := Exp(-y) * (Cos(x) + ComplexField().1*Sin(x)); f := two*pi/w2; f2 := f*f; f4 := f2*f2; term := one; qpower := one; sum4 := zero; sum6 := zero; tiny := 10^(-30)*1.0; for n in [1..2000] do if Abs(term) lt tiny then break ; end if; n2 := n^2; qpower *:= q; term := n*n2*qpower/(one-qpower); sum4 +:= term; term *:= n2; sum6 +:= term; end for; if Abs(term) gt 10^(-4) then "term = ",term; "ellinv_getc4c6: WARNING: precision very bad!"; end if; c4 := (one + 240*sum4)*f4; c6 := (one - 504*sum6)*f4*f2; return c4, c6; end function; intrinsic EllipticInvariants(A::ModTupFld) -> SeqEnum {} return EllipticInvariants(A,DefaultPrec(A)); end intrinsic; intrinsic EllipticInvariants(A::ModTupFld, n::RngIntElt) -> SeqEnum {Computes w1, w2, c4, c6, and j to precision n.} assert DimensionAbelianVariety(A) eq 1; assert IsCuspidal(A); L := PeriodLattice(A,n); tau, w1, w2 := ellinv_NormalizePair(L[1][1], L[2][1]); c4, c6 := ellinv_c4c6(w1,w2); j := 1728*c4*c4*c4/(c4*c4*c4 - c6*c6); return [w1, w2, c4, c6, j]; end intrinsic; intrinsic AbelianVariety(A::ModTupFld) -> SeqEnum {Returns the abelian variety associated to A.} if DimensionAbelianVariety(A) gt 1 then error "Not yet programmed for dimension > 1"; end if; return EllipticCurve(A); end intrinsic; intrinsic EllipticCurve(A::ModTupFld) -> SeqEnum {Computes an elliptic curve associated to the weight-two rational modular symbols factor A (this curve is well-defined up to isogeny). By the Shimura-Taniyama theorem the computed curve is correct (assuming correctness of our implementation of the algorithms).} if IsPlusQuotient(A) then error "EllipticCurve: A must not be a +1 quotient."; end if; if IsMinusQuotient(A) then error "EllipticCurve: A must not be a -1 quotient."; end if; if Weight(A) ne 2 then error "EllipticCurve: A must be of weight 2."; end if; if not IsCuspidal(A) then error "EllipticCurve: A must be cuspidal."; end if; if not IsTrivial(Character(A)) then error "EllipticCurve: A must be on Gamma_0(N), i.e., have trivial character."; end if; if DimensionAbelianVariety(A) ne 1 then error "EllipticCurve: A must have rational fourier coefficients."; end if; // TODO: something more intelligent. sturm := HeckeBound(ModSymParent(A)); n := Max(Integers()!Round(Conductor(A)/3), sturm); f := qEigenform(A,n); n := Degree(f); // it might be known to much higher precision... w1, w2, c4, c6, j := Explode(EllipticInvariants(A,n)); if HeckeVerbose then "c4 =",c4; "c6 =",c6; end if; c4 := Integers()!Round(c4); c6 := Integers()!Round(c6); // Find an elliptic curve, whose c4 and c6 are close and // whose ap for p good and less than n equal those of A. maxtries := 1000; tries := 1; rad := 1; cur := 0; dir := 0; // 0=right, 1=down, 2=left, 3=up. seen := {@ [c4,c6] @}; // index set of positions seen so far. if HeckeVerbose then "Searching..."; end if; while true do if HeckeSuperverbose then printf "[%o,%o] ",c4,c6; if tries mod 6 eq 0 then ""; end if; end if; E := EllipticCurve([-27*c4,-54*c6]); Include(~seen, [c4,c6]); if Conductor(E) eq Conductor(A) then found := true; if HeckeVerbose then printf "Candidate curve [%o,%o]:\n",c4,c6; end if; EE := MinimalModel(E); for p in [1..sturm] do if IsPrime(p) and (Conductor(A) mod p ne 0) then EFp := BaseExtend(EE,GF(p)); ap := TraceOfFrobenius(EFp) ; if ap ne Coefficient(f,p) then if HeckeVerbose then "Rejecting--the ",p,"th L-series coefficients disagree"; end if; found := false; break; end if; end if; end for; if found eq true then if HeckeVerbose then "By Shimura-Taniyama, this is in the isogeny class."; end if; return E; end if; end if; if tries gt maxtries then if HeckeVerbose then "\nRecomputing period integrals to higher precision"; end if; n := Integers()!Round(Degree(f) + 50); w1, w2, c4, c6, j := Explode(EllipticInvariants(A,n)); if HeckeVerbose then "c4 =",c4; "c6 =",c6; end if; c4 := Integers()!Round(c4); c6 := Integers()!Round(c6); f := qEigenform(A,n); tries := 1; rad := 1; cur := 0; dir := 0; end if; if Index(seen,[c4,c6]) ne 0 then // move to next pair [c4,c6]. tries +:= 1; repeat if cur ge rad then // turn dir := (dir + 1) mod 4; cur := 0; if dir eq 3 then cur := 1; c4 -:= 1; rad +:= 2; continue; end if; end if; case dir: // 0=right, 1=down, 2=left, 3=up. when 0: c4 +:= 1; when 1: c6 +:= 1; when 2: c4 -:= 1; when 3: c6 -:= 1; end case; cur +:= 1; until Index(seen,[c4,c6]) eq 0; end if; end while; // never reached. end intrinsic; intrinsic LAnalytic(A::ModTupFld, j::RngIntElt) -> FldPrElt {} return LAnalytic(A,j,DefaultPrec(A)); end intrinsic; intrinsic LAnalytic(A::ModTupFld, j::RngIntElt, n::RngIntElt) -> FldPrElt {Compute L(A,j) for j a "critical integer", i.e., one of the integers 1,...,k-1. Use at least n terms of the q-expansions.} requirege j,1; /*********************************************************** * * * L(f,j) = -(2pi)^(j-1)/(j-1)! * * *(i)^(j-1) < f, X^(j-1)*Y^(k-2-(j-1)) {0,oo}>.* * L(A,j) := prod L(f,j) over the conjugate f's. * * * ***********************************************************/ k := Weight(A); assert j le k-1; C := ComplexField(); R := PolynomialRing(Rationals(),2); x := -ConvFromModularSymbol(ModSymParent(A), ); Phi := PeriodMapping(A,n); z := ApplyPeriodMapping(Phi, x); z := Eltseq((2*Pi(C))^(j-1)/Factorial(j-1)*(i)^(j-1) * z); // The vector z now contains the values L(g,j) as g varies // over our basis of modular forms with integer Fourier expansion. if #z eq 1 then return z[1]; // nothing more to do. end if; if &+[Abs(zz) : zz in z] eq 0 then return 0; end if; e := EigenformInTermsOfIntegralBasis(A); Q := Parent(e[1]); PC := PolynomialRing(ComplexField()); g := Modulus(Q); P := PreimageRing(Q); e := [P!f : f in e]; roots := Roots(PC!g); // find all conjugates of eigen... assert #roots eq Degree(g); // g must be square free! Lfval := &*[ &+[Evaluate(PC!e[i],roots[j][1])*z[i] : i in [1..#z]] : j in [1..#roots]]; return Lfval; end intrinsic; intrinsic ClassicalPeriod(A::ModTupFld, j::RngIntElt) -> FldPrElt {} return ClassicalPeriod(A, j, DefaultPrec(A)); end intrinsic; intrinsic ClassicalPeriod(A::ModTupFld, j::RngIntElt, n::RngIntElt) -> FldPrElt {Returns the classical period r_j(f) = int_\{0\}^\{ioo\} f(z) z^j dz.} C := ComplexField(); return Factorial(j)*i^(j+1)*(2*Pi(C))^(-(j+1))* LAnalytic(A,j+1,n); end intrinsic; ////////////////////////////////////////////////////////////////// // // // Code for computing derivatives of L-series. // // // ////////////////////////////////////////////////////////////////// // The code for DerivGr is ported straight from Cremona C++ package. function is_approx_zero(x) return Abs(x) lt 10^(-10); end function; intrinsic DerivG1(x::FldPrElt) -> FldPrElt {Compute the G_1 function defined in section 2.13 of Cremona's Algorithms book.} if x gt 708 then // it's what John does... return 0; end if; if x lt 2 then ans := -Log(x) - EULER; p := -1; if HeckeSuperverbose then printf "Computing DerivG1 for x = %o using series\n",x; end if; ok := 0; // The following does not work for large x! (The terms get too big // before they get small, and overflow destroys the result.) for n in [1..5000] do p /:= n ; p *:= -x ; term := p/n; ans +:= term; if is_approx_zero(term/ans) then break; end if; end for; return ans; else // else x>2, use continued fraction form from B-G-Z p.478, // or section 2.13 of Cremona. if HeckeSuperverbose then printf "Computing DeriveG1 for x = %o using continued fraction\n",x; end if; a0 := 0.0; b0 := 1.0; b1 := x; ans:= 0.0; a1 := Exp(-x); for k in [2..20000] do if IsOdd(k) then ca := Floor((k-1)/2); a2 := x*a1+ca*a0; a0 := a1; a1 := a2; b2 := x*b1+ca*b0; b0 := b1; b1 := b2; else ca := Floor(k/2); a2 := a1+ca*a0; a0 := a1; a1 := a2; b2 := b1+ca*b0; b0 := b1; b1 := b2; end if; newans := a2/b2; if is_approx_zero(ans-newans) then return newans; end if; ans:=newans; end for; end if; "In function g1, continued fraction method, reached end of loop!"; return 0; end intrinsic; intrinsic DerivG2(x::FldPrElt) -> FldPrElt {Compute the G_2 function defined in section 2.13 of Cremona's Algorithms book.} if (x gt 20) then return 0; end if; ans := -Log(x) - EULER; if HeckeSuperverbose then "Computing myg2 for x = ",x; end if; ans := ans*ans/2 + Pi(RealField())*Pi(RealField())/12; p := 1; for n in [1..500] do p /:= n ; p *:= -x ; term := (p/n)/n; ans +:= term; if is_approx_zero(term/ans) then return ans; end if; end for; "In function g2, reached end of loop!"; return ans; end intrinsic; intrinsic DerivG3(x::FldPrElt) -> FldPrElt {Compute the G_3 function defined in section 2.13 of Cremona's Algorithms book.} if (x gt 20) then return 0; end if; ans := -Log(x) - EULER; if HeckeSuperverbose then "Computing myg3 for x = ",x; end if; ans := (Pi(RealField())^2 + 2*ans^2) * ans/12 - 1.20205690315959428540/3; p := -1; for n in [1..500] do p *:= -x/n; term := p/(n^3); ans +:= term; if is_approx_zero(term/ans) then return ans; end if; end for; "In function g3, reached end of loop!"; return ans; end intrinsic; intrinsic DerivWt2(f::RngSerPowElt, N::RngIntElt, r::RngIntElt) -> FldPrElt {Compute (d!)*Prod L^(dr)(fi,1) over the d Galois conjugates of f. TODO: When both r > 1 and d>1 the normalization is wrong! It's trivial combinatoric to work out but I have not done it.} requirege r,0; requirege N,1; // Essentially we compute the formula in Proposition 2.13.1 of [Cr]. case r: when 1: Gr := DerivG1; when 2: Gr := DerivG2; when 3: Gr := DerivG3; else: assert "r must be 1,2, or 3."; end case; rootN := Sqrt(N); G := [Gr(2*PI*n/rootN) : n in [1..Degree(f)]]; Q := Parent(Coefficient(f,1)); if Type(Q) ne FldRat then PC := PolynomialRing(ComplexField()); g := Modulus(Q); d := Degree(g); P := PreimageRing(Q); roots := Roots(PC!g); // find all conjugates of eigen... assert #roots eq Degree(g); // g must be square free! ans := &*[ &+[Evaluate(PC!Coefficient(f,n)/n,roots[j][1]) * G[n] : n in [1..Degree(f)]] : j in [1..#roots]]; ans *:= Factorial(d)*2^d*Factorial(r)^d; else ans := 2*Factorial(r)* &+[Coefficient(f,n)/n * G[n] : n in [1..Degree(f)]]; end if; return ans; end intrinsic; function ComputeSpecialValueAndRankWt2(A,n) f := qEigenform(A,n); N := Conductor(A); e := SignWN(A); if e eq -1 and LRatio(A) ne 0 then return LAnalytic(A,1,n), 0; end if; if e eq 1 then a := DerivWt2(f,N,1); if Abs(a) gt 10^(-3) then return a, 1; end if; if HeckeVerbose then "L'(A,1) ~ ", a, " hence computing second derivative."; end if; end if; if e eq -1 then a := DerivWt2(f,N,2); if Abs(a) gt 10^(-3) then return a, 2; end if; if HeckeVerbose then "L''(A,1) ~ ", a, " hence computing second derivative."; end if; end if; if e eq 1 then a := DerivWt2(f,N,3); if Abs(a) gt 10^(-3) then return a, 3; end if; if HeckeVerbose then "WARNING: L'''(A,1) ~ ", a, " even looks zero! Giving up."; end if; end if; "WARNING: not able to determine rank."; return a, -1; end function; intrinsic LeadingCoefficient(A::ModTupFld) -> FldPrElt {} return LeadingCoefficient(A,1,DefaultPrec(A)); end intrinsic; intrinsic LeadingCoefficient(A::ModTupFld, n::RngIntElt) -> FldPrElt {} return LeadingCoefficient(A,1,n); end intrinsic; intrinsic LeadingCoefficient(A::ModTupFld, j::RngIntElt, n::RngIntElt) -> FldPrElt, RngIntElt {Return the leading coefficient of Taylor expansion about the critical integer j and the analytic rank of the the L-series associated to A. The precision used is at least n terms of the q-expansion.} if IsTrivial(Character(A)) and Weight(A) eq 2 then Lval, r := ComputeSpecialValueAndRankWt2(A,n); return Lval/Factorial(r*DimensionAbelianVariety(A)), r; end if; error "Not programmed for weight > 2"; end intrinsic; // Test -- the rank 3 curve: /* D:=CremonaDatabase(); E:=EllipticCurve(D,"5077A"); R := PowerSeriesRing(Rationals()); f := R!CuspForm(E,17); // this function is in Kohel's code... DerivWt2(f, 5077, 1); DerivWt2(f, 5077, 3); */ /************** END OF ANALYTIC COMPUTATION CODE ********************/ /***************************************************************** SHIMURA'S MODULAR FORMS FORMULAS (Extended from a PARI program of Bruce Caskel and Kevin Buzzard.) *****************************************************************/ function mu0(n) return &*([n] cat [1+1/t[1]: t in Factorization(n)]); end function; function mu20(n) if n mod 4 eq 0 then return 0; end if; return &*[Integers()| 1+KroneckerSymbol(-4,t[1]): t in Factorization(n)]; end function; function mu30(n) if (n mod 2 eq 0) or (n mod 9 eq 0) then return 0 ; end if; return &*[Integers()| 1+KroneckerSymbol(-3,t[1]): t in Factorization(n)]; end function; function c0(n) return &+[EulerPhi(Gcd(d, Integers()!(n/d))) : d in Divisors(n)]; end function; function g0(n) return Integers()!(1+mu0(n)/12-mu20(n)/4-mu30(n)/3-c0(n)/2); end function; function mu1(n) if n le 2 then return mu0(n); else return Integers()!(EulerPhi(n)*mu0(n)/2); end if; end function; function mu21(n) if n lt 4 then return mu20(n); else return 0; end if; end function; function mu31(n) if n lt 4 then return mu30(n); else return 0; end if; end function; function c1(n) if n le 3 then return c0(n); end if; if n eq 4 then return 3; end if; return Integers()!(&+[EulerPhi(d)*EulerPhi(Integers()!(n/d))/2 : d in Divisors(n)]); end function; function g1(n) return Integers()!(1+mu1(n)/12-mu21(n)/4-mu31(n)/3-c1(n)/2); end function; function ss0(n,p) assert IsPrime(p) and (n mod p ne 0); return g0(n*p) - 2*g0(n) + 1; end function; function muXNp(n,p) return mu1(n)*mu0(p); end function; function mu2XNp(n,p) return 0; end function; function mu3XNp(n,p) return 0; end function; function cXNp(n,p) return 2*c1(n); end function; function gXNp(n,p) if n lt 4 then return g0(n*p); end if; return 1 + muXNp(n,p)/12 - mu2XNp(n,p)/4 - mu3XNp(n,p)/3 - cXNp(n,p)/2 ; end function; function ss1(n,p) assert IsPrime(p) and (n mod p ne 0); return gXNp(n,p) - 2*g1(n) + 1; end function; function eisen(p) assert IsPrime(p); return Numerator((p-1)/12); end function; function S0(n,k) assert n gt 0; if (k le 0) or (k mod 2 ne 0) then return 0; end if; if k eq 2 then return g0(n); end if; return Integers()!((k-1)*(g0(n)-1) + (k/2-1)*c0(n)+mu20(n)*Floor(k/4)+mu30(n)*Floor(k/3)); end function; function S1(n,k) assert n gt 0; if (k le 0) or ((n lt 3) and (k mod 2 ne 0)) then return 0; end if; if k eq 1 then error "S1, k = 1 not programmed."; end if; if k eq 2 then return g1(n); end if; if n lt 3 then return S0(n,k); end if; a := (k-1)*(g1(n)-1)+(k/2-1)*c1(n); if n eq 4 and k mod 2 ne 0 then a +:= 1/2; elif n eq 3 then a +:= Floor(k/3); end if; return Integers()!a; end function; function idxG0(n) return &*[Integers()| t[1]^t[2] + t[1]^(t[2]-1) : t in Factorization(n)]; end function; function idxG1(n) return EulerPhi(n)*idxG0(n); end function; /************** END OF SHIMURA'S FORMULAS ********************/ /***************************************************************** DIRICHLET CHARACTERS. *****************************************************************/ // This function isn't currently used. intrinsic CompleteCyclotomicField(p::RngIntElt, v::RngIntElt) -> .,. {Returns two arguments, a v-adic field K and an element zeta in K that is a primitive pth root of unity. Here v should be a prime number, or 0; if v=0 then K=the complex numbers.} if v eq 0 then F := CyclotomicField(p); K := ComplexField(); zeta := K!zeta; return K, zeta; end if; if not IsPrime(v) then error "CompleteCyclotomicField: the place v must be 0 or prime; however, it is ",v,"."; end if; F := pAdicField(v); R := PolynomialRing(F); f := Factorization(R!CyclotomicPolynomial(p))[1,1]; K := quo; return K, zeta; end intrinsic; CCharacter := recformat< N, // The modulos of epsilon. cond, // The conductor of epsilon. f, // Map (Z/NZ)^* --> Prod of copies of Z/aZ's. ims, // The image of each of the generators of (Z/NZ)^* // under the identification given by f. pn, // Ordered list of prime powers exactly dividing N. istrivial, // true iff this is the trivial character, // in which case less data is needed to store character. F, std // standard vector representation (see thesis) >; intrinsic IsTrivial(eps::Rec) -> BoolElt { Returns true iff the character eps is the trivial character. } return eps`istrivial eq true; end intrinsic; intrinsic Modulus(eps::Rec) -> RngIntElt { Returns the level of the character eps. } return eps`N; end intrinsic; intrinsic CharacterLevel(eps::Rec) -> RngIntElt { Returns the level of the character eps. } return eps`N; end intrinsic; intrinsic IsOdd(eps::Rec) -> BoolElt { Return true iff the Dirichlet character is odd.} return Evaluate(eps,-1) eq -1; end intrinsic; intrinsic FieldOfValues(eps::Rec) -> RngIntElt { Returns the field in which the Dirichlet character eps takes values. } return eps`F; end intrinsic; intrinsic DirichletCharacterImages(eps::Rec) -> SeqEnum {} return eps`ims; end intrinsic; intrinsic DirichletCharacter(N::RngIntElt, std::ModEDElt, p::RngIntElt) -> Rec {Using the standard vector std, this function creates a Dirichlet character modulo N that takes values in a field of characteristic p.} lift := [Integers()!std[i] : i in [1..Dimension(Parent(std))]]; e := Characteristic(Parent(std[1])); d := Gcd(lift cat [e]); r := e div d; "d = ",d; "e = ",e; pows := [lift[i] div d : i in [1..#lift]]; "pows = ",pows; return DirichletCharacter(N, pows, p, r); end intrinsic; intrinsic DirichletCharacter(N::RngIntElt, pows::SeqEnum, p::RngIntElt, r::RngIntElt : allowcomplex := false) -> Rec {Create any Dirichlet character modulo N in characteristic p. The input "pows" is a list of nonnegative integers. There is an isomorphism Z/NZ ---> (Z/p1^(i1)Z) x (Z/(p2)^(i2)Z) x ... with p1 < p2 < ... the primes dividing N. Passing to units we find that the right hand side is a cyclic groups, except possibly when p1=2 in which case the cyclic group is trivial, or a group of order 2 times a cyclic group (possibly trivial). Let zeta be a primitive rth root of unity in the algebraic closure of the prime field of characteristic p. The pows specify the power of zeta that a primitive root (computed using PrimitiveRoot(p^i)) of each Z/p^iZ maps to. The primitive root used is the one computed using the function PrimitiveRoot(p^i). In the special case p1=2 and i1>=3, the user must supply two entries in the pows sequence. (But if i1=1 or 2, then the user must supply exactly one entry.) The group (Z/2^nZ)^* is generated by -1 and 5. The entries in the pows sequence correspond bijectively, and in order, to the prime divisors of N, listed in increasing order. } if N lt 1 then error "The level N must be >= 1."; end if; /******************************************************************* Record the _standard_ vector format description of the character. Let e be the prime-to-p part of the exponent of (Z/NZ)^*. and set z := exp^{2pi*i/e}. We assume that zeta := exp^{2pi*i/r}. Let m be minimal so that r/m divides e; then, the pows must all be divisible by m, since zeta^pows[i] must be a power of z. Let n = e / (r/m). The standard vector is got by multiplying each pows[i] by n/m; each entry in the standard vector must divide e. *******************************************************************/ e := Exponent(UnitGroup(Integers(N))); m := r div Gcd(r,e); if p ne 0 then e := e div p^Valuation(e,p); end if; n := e div Gcd(r,e); // Magma doesn't allow creating Z/1Z using Integers(1)! std := RModule(quo, #pows)![pows[i]*(n div m) : i in [1..#pows]]; /*************************************************************/ // Compute the internal representation of the character. /*************************************************************/ if p eq 0 then if r gt 2 then F:=CyclotomicField(r); if allowcomplex then F := ComplexField(); zeta := F!zeta; end if; elif r eq 2 then F := Rationals(); zeta := F!-1; else F := Rationals(); zeta := F!1; end if; if Min(pows) lt 0 then "COMPLEX: By entering a negative power, you've elected to"; "work with modular symbols over the Magma complexes. This"; "is a very bad idea with no practical value. Good luck."; C:=ComplexField(); zeta:=C!zeta; F:=C; for i in [1..#pows] do if pows[i] lt 0 then pows[i] := -pows[i]; end if; end for; end if; else if Gcd(p,r) ne 1 then error "Can not create a mod ", p, " character of degree ", r,"."; end if; if r eq 1 then zeta := GF(p)!1; else Zr := IntegerRing(r); s := Order(Zr!p); ps := p^s; F := GF(ps); z := PrimitiveElement(F); zeta := z^(Integers()!((ps-1)/r)); end if; end if; ZnZ := IntegerRing(N); U, f := UnitGroup(ZnZ); Ugens := Generators(U); gens := [(U.i)@@f : i in [1..#Ugens]]; // DON'T iterate Ugens (no order!) ims := []; fac := Factorization(N); pn := [fac[i][1]^fac[i][2] : i in [1..#fac]]; offset := 0; pr := p eq 0 and (Type(F!0) eq FldPrElt); for j := 1 to #fac do phi := pn[j] - fac[j][1]^(fac[j][2]-1); if fac[j][1] mod 2 ne 0 or fac[j][2] le 2 then if #pows lt j+offset then error "Invalid character: not enough data at",fac[j][1]; end if; if (not pr and zeta^(phi*pows[j+offset]) ne 1) or (pr and Abs(zeta^(phi*pows[j+offset])-1) gt 10^(-10)) then print "phi := ",phi; print "pow at ? := ",pows[j+offset]; print "zeta is a prim ",r,"th root of 1"; error "Invalid character definition at ", fac[j][1] ; end if; else offset:=1; if #pows lt j+1 then error "Invalid character: not enough data at 2."; end if; if (not pr and zeta^(2*pows[j]) ne 1) or (pr and Abs(zeta^(2*pows[j])-1) gt 10^(-10)) then error "Invalid character definition at 2"; end if; if (not pr and zeta^(Integers()!((pn[j]/4)*pows[j+1])) ne 1) or (pr and Abs(zeta^(Integers()!((pn[j]/4)*pows[j+1]))-1) gt 10^(-10)) then error "Invalid character definition at 2"; end if; end if; end for; // Find the image of each generator in the field containing zeta. for g in gens do // Compute image of generator under each of the maps // (Z/NZ)^* --> Prod (Z/p^aZ)^* --> C^* // then combine together the result. im := 1; offset := 0; for j := 1 to #fac do p := fac[j][1]; nu := fac[j][2]; pnu := pn[j]; G, h := UnitGroup(IntegerRing(pnu)); v:=ElementToSequence(h(IntegerRing(pnu)!g)); if #v eq 0 then elif #v eq 1 then im *:= zeta^(Integers()!(pows[j+offset]*v[1])); else // pn >= 8 // h(-1) = G.1, h(5) = G.2 (unless Magma is weird...) im *:= zeta^(pows[j]*v[1]+pows[j+1]*v[2]); offset:=1; end if; end for; Append(~ims, im); end for; eps := rec; if &and[a eq 1 : a in ims] then eps`istrivial := true; end if; eps`cond := Conductor(eps); return eps; end intrinsic; intrinsic DirichletCharacter(N::RngIntElt, degs::SeqEnum, characteristic::RngIntElt : allowcomplex := false) -> Rec {Create a Dirichlet character modulo N. The input zeta is a root of unity in some field. The second input "degs" is a list of nonnegative integers. We explain the notation we use for degs. There is an isomorphism Z/NZ ---> (Z/p1^(i1)Z) x (Z/(p2)^(i2)Z) x ... with p1 < p2 < ... the primes dividing N. Passing to units we find that each factor in the right hand side is a cyclic group, except possibly when p1=2 in which case the cyclic group is trivial, or a group of order 2 times a cyclic group (possibly trivial). The degs specify the order that the image of a primitive root of each Z/p^iZ should map to. The primitive root used is the one computed using the function PrimitiveRoot(p^i). In the special case p=2 and i>=3, one should supply supply two entries corresponding to 2 in the degs sequence, because the group (Z/2^nZ)^* is generated by -1 and 5. The entries in the gens sequence correspond bijectively, and in order, to the prime divisors of N, listed in increasing order. Note: Not every character can be described in this way. } if not (characteristic eq 0 or IsPrime(characteristic)) then printf "There is no field of characteristic %o.",characteristic; error ""; end if; // Simply translate to the format for the other // character creation function and then call it. l := Lcm(degs); pows := [Integers()!(l/d) : d in degs]; return DirichletCharacter(N, pows, characteristic, l : allowcomplex := allowcomplex); end intrinsic; intrinsic DirichletCharacter(N::RngIntElt, degs::SeqEnum : allowcomplex := false) -> Rec {Create Dirichlet character.} return DirichletCharacter(N, degs, 0 : allowcomplex:=allowcomplex); end intrinsic; intrinsic DirichletCharacter(N::RngIntElt) -> Rec {Create the trivial Dirichlet character.} if N lt 1 then error "The level N must be >= 1."; end if; return rec; end intrinsic; intrinsic DirichletCharacter(N::RngIntElt, p::RngIntElt) -> Rec {Create the trivial Dirichlet character with values in the prime field of characteristic p>=0.} return rec; end intrinsic; intrinsic Evaluate(eps::Rec, gamma::SeqEnum) -> Rec {Evaluate the character eps on the matrix gamma=[a,b, c,d] in Gamma_0(N). The result is eps(a).} return Evaluate(eps, gamma[1]); end intrinsic; intrinsic Evaluate(eps::Rec, n::RngIntElt) -> Rec {Evaluate the character eps at the integer n.} /*if eps`istrivial then return 1; end if;*/ return Evaluate(eps, quo!n); end intrinsic; intrinsic Evaluate(eps::Rec, n::RngIntResElt) -> Rec {Evaluate the character eps at the integer (n mod N).} if Gcd(n,eps`N) ne 1 then return (eps`F)!0; end if; if eps`istrivial then return (eps`F)!1; end if; v:=ElementToSequence(eps`f(n)); return &*[(eps`ims)[i]^v[i] : i in [1..#v]]; end intrinsic; intrinsic EvaluateAtp(eps::Rec, p::RngIntElt, n::RngIntResElt) -> Rec {Evaluate the character eps at the integer n.} return EvaluateAtp(eps, p, Integers()!n); end intrinsic; intrinsic EvaluateAtp(eps::Rec, p::RngIntElt, n::RngIntElt) -> Rec {Let eps_p denote the character obtained by restricting eps to the p-factor. This function returns eps_p(n). When gcd(p,n)>1 this is 0. When gcd(p,n)=1, let m be an element of (Z/NZ)^* which is 1 mod q for all q neq p and is n at p. Then eps_p(n) = eps(m).} if Gcd(n,p) ne 1 then return 0; end if; if eps`istrivial then return 1; end if; pn:=eps`pn; X := [(pn[i] mod p eq 0) select n else 1 : i in [1..#pn]]; return Evaluate(eps, CRT(X,pn)); end intrinsic; intrinsic Conductor(eps::Rec) -> RngIntElt {} N := Modulus(eps); std := Eltseq(eps`std); R,h := AdditiveGroup(Parent(std[1])); f := 1; // conductor if IsEven(N) then v := Valuation(N,2); if v eq 1 then Remove(~std,1); elif v eq 2 then if std[1] ne 0 then f *:= 4; end if; Remove(~std,1); else // 8 | N if std[2] eq 0 then if std[1] ne 0 then f *:= 4; end if; else f *:= 4*Order(h(std[2])); end if; Remove(~std,1); Remove(~std,1); end if; end if; fac := Factorization(OddPart(N)); for i in [1..#fac] do if std[i] ne 0 then R := IntegerRing(Order(h(std[i]))); p := R!fac[i][1]; for d in [1..fac[i][2]] do if p^d eq p^(d-1) then f *:= fac[i][1]^d; break; end if; end for; end if; end for; return f; end intrinsic; //////////////////////////////////////////////////////////////////////// // Computed associated mod M character. // // Algorithm: Suppose a1,...,an generate (Z/NZ)^* independently. // // Let a1bar,...,anbar be the images of the a1,...,an in // // (Z/MZ)^*. Consider the images eps(a1),...,eps(an) of // // a1,...,an in C^*. If for each i the divisibility // // order(eps(ai)) | order(aibar) // // then it is possible to define a character // // eps': (Z/MZ)^* ---> C^* such that eps = eps' o red. // // This character is defined by sending aibar to eps(ai). // // Conversely, if one of the above divisibilities fails // // then it is evidently impossible for such an eps' to exist. // //////////////////////////////////////////////////////////////////////// intrinsic CanReduceModM(eps::Rec, M::RngIntElt) -> BoolElt {Returns true iff M>=1 divides Modulus(eps) and is divisible by Conductor(eps).} if M lt 1 or Modulus(eps) mod M ne 0 then return false; end if; if IsTrivial(eps) then return true; end if; if M eq 1 then return false; end if; return M mod Conductor(eps) eq 0; end intrinsic; intrinsic AssociatedModMCharacter(eps::Rec, M::RngIntElt) -> Rec {Given a Dirichlet character eps, such that M divides Level(eps) and Conductor(eps) divides M, return the associated mod M Dirichlet character.} if M eq Modulus(eps) then return eps; end if; if IsTrivial(eps) then return DirichletCharacter(M,Characteristic(FieldOfValues(eps))); end if; if M le 1 then error "The conductor of eps does not divide M."; end if; if M gt Modulus(eps) then error "Character extension not yet programmed."; end if; f := eps`f; U := Codomain(f); a := [u@@f : u in Generators(U)]; ZMZ := IntegerRing(M); abar := [ZMZ!ai : ai in a]; // check if conductor divides M. for i in [1..#a] do if Evaluate(eps,a[i])^Order(abar[i]) ne 1 then error "The conductor of eps does not divide M."; end if; end for; // now construct the new character: UM, h:= UnitGroup(ZMZ); red := [h(x) : x in abar]; for i in [1..#Generators(UM)] do // Generators(UM) not ordered! if Index(red,UM.i) eq 0 then error "Magma's unit group function is not as canonical as I hoped."; end if; end for; ims := [Evaluate(eps,a[Index(red,UM.i)]) : i in [1..#Generators(UM)]]; epsM := rec; if &and[a eq 1 : a in ims] then epsM`istrivial := true; end if; return epsM; end intrinsic; function AllStandardCharacters(N) e := Exponent(UnitGroup(Integers(N))); f := Factorization(N); phi := [f[i][1]^(f[i][2]-1)*(f[i][1] - 1) : i in [1..#f]]; r := #f; s := IsEven(N) select r+1 else r; R := RModule(Integers(e),s); a := [R!0 : i in [1..s]]; if N mod 4 eq 0 then a[1][1] := e div 2; end if; if N mod 8 eq 0 then a[2][2] := e div (2^(Valuation(N,2)-2)); end if; if N mod 2 eq 0 then i := 3; else i := 1; end if; for j in [1..#f] do fac := f[j]; if fac[1] eq 2 then continue; end if; a[i][i] := e div phi[j]; i +:= 1; end for; // Iterate through the submodule generated by the elements // of the vector a. A := [R|]; cur := R!0; while true do Append(~A,cur); cur +:= a[1]; // carry i := 1; while cur[i] eq 0 and i lt Dimension(R) do i +:= 1; cur +:= a[i]; end while; if cur[i] eq 0 and i eq Dimension(R) then break; end if; end while; return A; end function; intrinsic AllDirichletCharacters(N::RngIntElt) -> SeqEnum {Returns a sequence consisting of all of the mod N Dirichlet characters.} A := AllStandardCharacters(N); return [DirichletCharacter(N,A[i],0) : i in [1..#A]], A; end intrinsic; intrinsic AllDirichletCharacterClasses(N::RngIntElt) -> SeqEnum {Returns a sequence consisting of the Gal(Qbar/Q)-conjugacy classes of mod N Dirichlet characters.} e := Exponent(UnitGroup(Integers(N))); end intrinsic; intrinsic Print(eps::Rec) {Print a Dirichlet character. This is a temporary function, while I wait for classes or "hack objects".} printf "Mod %o Dirichlet character: ", eps`N; if IsTrivial(eps) then print "trivial."; else print eps`std; end if; end intrinsic; /************** END OF DIRICHLET CHARACTERS ********************/ /************************************************************************ * * * MODELS FOR MODULAR CURVES * * * ************************************************************************/ intrinsic LinearRelations(P::SeqEnum,prec::RngIntElt) -> SeqEnum {Let P be a sequence of power series (or polynomials) over the integers. This function returns a basis for the Z-linear relations satisfied by these series.} X := [[Coefficient(f,i) : i in [0..prec]] : f in P]; M := RMatrixSpace(Integers(),#X, #X[1])! (&cat X); return Kernel(M); end intrinsic; function MonomialsOfDegreeD_recurse(P, D, R) if #P eq 0 then return ; end if; if D eq 1 then return ; end if; Q := [P[i] : i in [2..#P]]; QQ := [MonomialsOfDegreeD_recurse(Q,d,sub) : d in [1..D]]; serpart := [P[1]^D] cat &cat[ [P[1]^d*f : f in QQ[D-d][1]] : d in [0..D-1]]; polpart := [R.1^D] cat &cat[ [R.1^d*m : m in QQ[D-d][2]] : d in [0..D-1]]; return ; end function; intrinsic MonomialsOfDegreeD(P::SeqEnum, D::RngIntElt) -> Tup {} return MonomialsOfDegreeD_recurse(P,D,PolynomialRing(Rationals(),#P)); end intrinsic; intrinsic RelationsOfDegreeD(P::SeqEnum, D::RngIntElt, prec::RngIntElt) -> SeqEnum, SeqEnum {} m := MonomialsOfDegreeD(P,D); return LinearRelations(m[1],prec), m[2]; end intrinsic; intrinsic RelationsOfDegreeD(P::SeqEnum, D::RngIntElt) -> SeqEnum {} prec := Min([Degree(f) : f in P]); m := MonomialsOfDegreeD(P,D); r := Basis(LinearRelations(m[1],prec)); R := PolynomialRing(Integers(),#P); return [R!&+[m[2][i] * b[i] : i in [1..#m[2]]] : b in r]; end intrinsic; /***************************************************************** HIJIKATA TRACE FORMULA FOR Tr(1) = dim S_k(Gamma_1(N),eps). *****************************************************************/ function hij_a(s,k,N) case s: when 0: case k mod 4: when 0: return -1/2; when 2: return 1/2; else return 0; end case; when -1: case k mod 3: when 1: return 0; when 2: return 1/2; when 0: return -1/2; end case; when 1: case k mod 6: when 1: return 0; when 2: return 1/2; when 3: return 1/2; when 4: return 0; when 5: return -1/2; when 0: return -1/2; end case; when -2: return (-1)^k/(4*N); when 2: return 1/(4*N); else error "Hijikata: internal bug. (code a)"; end case; end function; function hij_b(s,f) case Abs(s): when 2: return 1; when 1: return 1/3; when 0: return 1/2; else error "Hijikata: internal bug. (code b)"; end case; end function; function hij_csfp(s,f,p,N,eps) rho := Valuation(f,p); nu := Valuation(N,p); case Abs(s): when 0: case p mod 4: when 2: if nu eq 1 and rho eq 0 then return 1; else return 0; end if; when 3: return 0; when 1: if rho gt 0 then return 0; else // x is a fourth root of 1 in (Z/p^{nu}Z)^* x := Sqrt(IntegerRing(p^nu)!(-1)); return EvaluateAtp(eps, p, x) + EvaluateAtp(eps,p,-x); end if; else error "Hijikata: internal bug. (code csfp)"; end case; when 1: if p eq 2 then return 0 ; elif (rho gt 0 or KroneckerSymbol(-3,p) eq -1) or (p eq 3 and nu gt 1) then return 0; else R := IntegerRing(p^nu); z := Sqrt(R!(-3)); w := R!(1/2); x1 := (1+z)*w; x2 := (1-z)*w; if p ne 3 then return EvaluateAtp(eps, p, s*x1) + EvaluateAtp(eps, p, s*x2); else return EvaluateAtp(eps, p, s*x1); end if; end if; when 2: return &+[EvaluateAtp(eps, p, Integers()!(2/s + i*p^(Ceiling(nu/2)+rho))): i in [1..p^Floor(nu/2)]] + &+[EvaluateAtp(eps, p, Integers()!(2/s + j*p^(Ceiling((nu+1)/2)+rho))): j in [1..p^Floor((nu-1)/2)]]; else error "Hijikata: internal bug. (code csfp)"; end case; end function; function hij_c(s,f,N,eps) return &*[Integers()| hij_csfp(s, f, p[1], N,eps) : p in Factorization(N)]; end function; function hij_bc(s,f,N,eps) return hij_b(s,f) * hij_c(s,f,N,eps); end function; function hij_abc(k, N,eps) return hij_a(0, k, N) * hij_bc(0, 1, N,eps) + hij_a(-1, k, N) * hij_bc(-1, 1, N,eps) + hij_a(1, k, N) * hij_bc(1, 1, N,eps) + hij_a(-2, k, N) * (&+[hij_bc(-2, f, N,eps) : f in [1..N]]) + hij_a(2, k, N) * (&+[hij_bc(2, f, N,eps) : f in [1..N]]); end function; intrinsic DimensionSk(eps::Rec, k::RngIntElt) -> RngIntElt {Compute the dimension of the space of cusp forms of character eps and weight k.} if Characteristic(FieldOfValues(eps)) ne 0 then error "DimensionCuspforms: Character must take values in a field of characteristic zero."; end if; N := CharacterLevel(eps); assert k gt 1; if IsTrivial(eps) then return S0(N,k); end if; if Evaluate(eps, N-1) ne (-1)^k then return 0; end if; return Integers()!(idxG0(N) * (k-1)/12 - hij_abc(k,N,eps)); end intrinsic; function CohenOesterle(eps, k) N := Modulus(eps); facN := Factorization(N); f := Conductor(eps); facf := Factorization(f); gamma_k := 0; if k mod 4 eq 2 then gamma_k eq -1/4; elif k mod 4 eq 0 then gamma_k eq 1/4; end if; mu_k := 0; if k mod 3 eq 2 then mu_k := -1/3; elif k mod 3 eq 0 then mu_k := 1/3; end if; function lambda(r,s,p) if 2*s le r then if IsEven(r) then return p^(r/2) + p^(r/2-1); else return p^((r-1)/2); end if; else return 2*p^(r-s); end if; end function; return -(1/2) * &*[lambda(facN[i][2],facf[i][2],facN[i][1]) : i in [1..#facN]] +gamma_k * &+[Evaluate(eps,x) : x in Integers(N) | x^2+1 eq 0] +mu_k * &+[Evaluate(eps,x) : x in Integers(N) | x^2+x+1 eq 0]; end function; intrinsic DimSk(eps::Rec, k::RngIntElt) -> RngIntElt {Compute the dimension of the space of cusp forms of character eps and weight k.} if Characteristic(FieldOfValues(eps)) ne 0 then error "DimensionCuspforms: Character must take values in a field of characteristic zero."; end if; N := CharacterLevel(eps); assert k gt 1; if IsTrivial(eps) then return S0(N,k); end if; if Evaluate(eps, N-1) ne (-1)^k then return 0; end if; return Integers()!(idxG0(N) * (k-1)/12 + CohenOesterle(eps,k)); end intrinsic; intrinsic DimensionSk(N::RngIntElt, k::RngIntElt) -> RngIntElt {Compute the dimension of the space of cusp forms on Gamma_0(N) of weight k.} return DimensionSk(DirichletCharacter(N), k); end intrinsic; intrinsic DimensionSkG1(N::RngIntElt, k::RngIntElt) -> RngIntElt {Returns dim S_k(Gamma_1(N)).} return S1(N,k); end intrinsic; intrinsic DimensionS2(N::RngIntElt) -> RngIntElt {} return DimensionSk(N,2); end intrinsic; function mumu(N) if Type(N) ne RngIntElt or N lt 1 then error "mumu(N): N must be a positive integer."; end if; p := 1; m := Factorization(N); for m in Factorization(N) do if m[2] gt 2 then p := 0 ; elif m[2] eq 1 then p := -2*p; end if; end for; return p; end function; intrinsic DimensionSknew(N::RngIntElt, k::RngIntElt) -> RngIntElt {Compute the dimension of the new subspace of S_k(Gamma_0(N)).} return &+[DimensionSk(M,k)*mumu(Integers()!(N/M)) : M in Divisors(N)]; end intrinsic; intrinsic DimensionSkG1new(N::RngIntElt, k::RngIntElt) -> RngIntElt {Compute the dimension of the new subspace of S_k(Gamma_0(N)).} return &+[DimensionSkG1(M,k)*mumu(Integers()!(N/M)) : M in Divisors(N)]; end intrinsic; /********************************************************************** * * * Mestre's * * METHOD OF GRAPHS * * Mestre, J.-F., La m\'ethode des graphes. Exemples et applications * * Proceedings of the international conference on class numbers * * and fundamental units of algebraic number fields (Katata) * * 217--242, Nagoya Univ., Nagoya, 1986. * * * **********************************************************************/ SSTABLE := [ // if negative, uses Kronecker test, if positive then is p. [-1,1728], [-2,8000], [-3,0], [-7,255^3], [-11,(-32)^3], [-19,-884736], [-43,-884736000], [-67, (-5280)^3], [-163, -262537412640768000], [15073,5797], [18313,10629], [38833,16622], [69337,2840], [71593,33119], [242713,87914], [448177,18081], [708481,79970], [716857, 389681], [882289,589416], [886537, 212553] ]; intrinsic FindSupersingularJ(p::RngIntElt) -> FldFinElt {Returns a supersingular j-invariant in characteristic p < 10^6.} requirerange p,2,1000000; for s in SSTABLE do if (s[1] lt 0 and KroneckerSymbol(s[1],p) eq -1) or s[1] eq p then return GF(p^2)!s[2]; end if; end for; end intrinsic; /////////////////////////////////////////////////////////////////// // MODULAR POLYNOMIALS Phi_{ell} for ell=2,3,5,7. // /////////////////////////////////////////////////////////////////// R := PolynomialRing(Integers(),2); PHI := [0, // Phi2: x^3 + j^3 - x^2*j^2 + 1488 * x * j * (x + j) - 162000 * (j^2 + x^2) + 40773375 * j * x + 8748000000 * (x + j) - 157464000000000, // Phi3: x*(x+2^15*3*5^3)^3+j*(j+2^15*3*5^3)^3-x^3*j^3 +2^3*3^2*31*x^2*j^2*(x+j)-2^2*3^3*9907*x*j*(x^2+j^2) +2*3^4*13*193*6367*x^2*j^2 +2^16*3^5*5^3*17*263*x*j*(x+j)-2^31*5^6*22973*x*j, 0, // Phi5: x^6 + (-j^5 + 3720*j^4 - 4550940*j^3 + 2028551200*j^2 - 246683410950*j + 1963211489280)*x^5 + (3720*j^5 + 1665999364600*j^4 + 107878928185336800*j^3 + 383083609779811215375*j^2 + 128541798906828816384000*j + 1284733132841424456253440)*x^4 + (-4550940*j^5 + 107878928185336800*j^4 - 441206965512914835246100*j^3 + 26898488858380731577417728000*j^2 - 192457934618928299655108231168000*j + 280244777828439527804321565297868800)*x^3 + (2028551200*j^5 + 383083609779811215375*j^4 + 26898488858380731577417728000*j^3 + 5110941777552418083110765199360000*j^2 + 36554736583949629295706472332656640000*j + 6692500042627997708487149415015068467200)*x^2 + (-246683410950*j^5 + 128541798906828816384000*j^4 - 192457934618928299655108231168000*j^3 + 36554736583949629295706472332656640000*j^2 - 264073457076620596259715790247978782949376*j + 53274330803424425450420160273356509151232000)*x + (j^6 + 1963211489280*j^5 + 1284733132841424456253440*j^4 + 280244777828439527804321565297868800*j^3 + 6692500042627997708487149415015068467200*j^2 + 53274330803424425450420160273356509151232000*j + 141359947154721358697753474691071362751004672000), 0, // Phi7: x^8 + (-j^7 + 5208*j^6 - 10246068*j^5 + 9437674400*j^4 - 4079701128594*j^3 + 720168419610864*j^2 - 34993297342013192*j + 104545516658688000)*x^7 + (5208*j^7 + 312598931380281*j^6 + 177089350028475373552*j^5 + 4460942463213898353207432*j^4 + 16125487429368412743622133040*j^3 + 10685207605419433304631062899228*j^2 + 1038063543615451121419229773824000*j + 3643255017844740441130401792000000)*x^6 + (-10246068*j^7 + 177089350028475373552*j^6 - 18300817137706889881369818348*j^5 + 14066810691825882583305340438456800*j^4 - 901645312135695263877115693740562092344*j^3 + 11269804827778129625111322263056523132928000*j^2 - 40689839325168186578698294668599003971584000000*j + 42320664241971721884753245384947305283584000000000)*x^5 + (9437674400*j^7 + 4460942463213898353207432*j^6 + 14066810691825882583305340438456800*j^5 + 88037255060655710247136461896264828390470*j^4 + 17972351380696034759035751584170427941396480000*j^3 + 308718989330868920558541707287296140145328128000000*j^2 + 553293497305121712634517214392820316998991872000000000*j + 41375720005635744770247248526572116368162816000000000000)*x^4 + (-4079701128594*j^7 + 16125487429368412743622133040*j^6 - 901645312135695263877115693740562092344*j^5 + 17972351380696034759035751584170427941396480000*j^4 - 5397554444336630396660447092290576395211374592000000*j^3 + 72269669689202948469186346100000679630099972096000000000*j^2 - 129686683986501811181602978946723823397619367936000000000000*j + 13483958224762213714698012883865296529472356352000000000000000)*x^3 + (720168419610864*j^7 + 10685207605419433304631062899228*j^6 + 11269804827778129625111322263056523132928000*j^5 + 308718989330868920558541707287296140145328128000000*j^4 + 72269669689202948469186346100000679630099972096000000000*j^3 - 46666007311089950798495647194817495401448341504000000000000*j^2 - 838538082798149465723818021032241603179964268544000000000000000*j + 1464765079488386840337633731737402825128271675392000000000000000000)*x^2 + (-34993297342013192*j^7 + 1038063543615451121419229773824000*j^6 - 40689839325168186578698294668599003971584000000*j^5 + 553293497305121712634517214392820316998991872000000000*j^4 - 129686683986501811181602978946723823397619367936000000000000*j^3 - 838538082798149465723818021032241603179964268544000000000000000*j^2 + 1221349308261453750252370983314569119494710493184000000000000000000*j)*x + (j^8 + 104545516658688000*j^7 + 3643255017844740441130401792000000*j^6 + 42320664241971721884753245384947305283584000000000*j^5 + 41375720005635744770247248526572116368162816000000000000*j^4 + 13483958224762213714698012883865296529472356352000000000000000*j^3 + 1464765079488386840337633731737402825128271675392000000000000000000*j^2) ]; intrinsic TpSS(p::RngIntElt, j::FldFinElt) -> SeqEnum {} F := Parent(j); R := PolynomialRing(F,2); S := PolynomialRing(F); h := hom S | S.1, j>; phi := h(R!PHI[p]); fac := Factorization(phi); if Max([Degree(a[1]) : a in fac]) gt 1 then print "Char = ",Characteristic(F); print "p = ",p; print "j = ",j; print "fac = ",fac; end if; assert Max([Degree(a[1]) : a in fac]) eq 1; return &cat[[-Evaluate(a[1],0) : i in [1..a[2]]] : a in fac]; end intrinsic; // Mestre class. declare attributes ModTupRng: p, M, ssbasis, weights, T, decomp, eisen; intrinsic Mestre(p::RngIntElt) -> ModTupRng {} return Mestre(p,1); end intrinsic; intrinsic Mestre(p::RngIntElt, M::RngIntElt) -> ModTupRng {Returns the Mestre method of graphs module.} assert Gcd(p,M) eq 1; assert IsPrime(p); if M ne 1 then error "Not programmed for M ne 1."; end if; D := RSpace(Integers(),DimensionSk(p,2)+1); B := Basis(D); // It is very important (efficiency-wise!) that the following be // given in Echelon-Form. X := sub; X`p := p; X`M := M; X`T := [* 0,0,0,0,0,0,0 *]; return X; end intrinsic; intrinsic TpD(X::ModTupRng, p::RngIntElt) -> AlgMatElt {Computes Tp on the Mestre graphs module. p must be a prime between 2 and 7.} requirerange p,2,7; assert p in {2,3,5,7}; if Type(X`T[p]) eq RngIntElt then if not assigned X`ssbasis then d := Degree(X); X`ssbasis := {@ FindSupersingularJ(X`p) @}; else d := #X`ssbasis; end if; V := RSpace(Integers(),d); T := MatrixAlgebra(Integers(),d)!0; for i in [1..d] do S := TpSS(p, X`ssbasis[i]); for j in S do Include(~X`ssbasis, j); T[i,Index(X`ssbasis,j)] +:= 1; end for; end for; X`T[p] := T; end if; return X`T[p]; end intrinsic; intrinsic TpX(X::ModTupRng, p::RngIntElt) -> AlgMatElt {Computes Tp on the subspace of degree 0 divisors in the Mestre module. p must be a prime between 2 and 7.} return Restrict(TpD(X,p),X); end intrinsic; intrinsic MestreEisensteinFactor(X::ModTupRng) -> ModTupRng {Computes the Eisenstein factor of the Mestre module X.} if not assigned X`eisen then p := 2; T := TpD(X,p)-(p+1); V := Kernel(T); while Dimension(V) gt 1 and p lt 7 do p := NextPrime(p); T := TpD(X,p)-(p+1); V := V meet Kernel(T); end while; if Dimension(V) gt 1 then error "Not enough Hecke operators to find Eisenstein series."; end if; X`eisen := V; end if; return X`eisen; end intrinsic; intrinsic MonodromyWeights(X::ModTupRng) -> SeqEnum {Returns the monodromy weights.} if not assigned X`weights then eis := Eltseq(Basis(MestreEisensteinFactor(X))[1]); n := Lcm([Integers()|e : e in eis]); X`weights := [Integers()|Abs(Numerator(n/e)) : e in eis]; end if; return X`weights; end intrinsic; intrinsic SupersingularBasis(X::ModTupRng) -> SeqEnum {Returns the basis for D of supersingular j-invariants.} if not assigned X`ssbasis then dummy:=TpD(X,2); end if; return X`ssbasis; end intrinsic; intrinsic Decomposition(X::ModTupRng) -> SeqEnum {Returns the decomposition of X into Hecke stable submodules under T2,T3,T5,T7.} if not assigned X`decomp then p := 2; Decomp := [ RSpace(Integers(),Degree(X)) ]; Irred := [ false ]; while p le 7 and not &and Irred do T := TpD(X,p); for V in [ V : V in Decomp | not Irred[Index(Decomp,V)] ] do fact := Factorization(ModularCharpoly(Restrict(T,V))); if #fact gt 1 then Remove(~Irred,Index(Decomp,V)); Exclude(~Decomp,V); for f in fact do Append(~Decomp,V meet Kernel(Evaluate(f[1],T))); if f[2] eq 1 then Append(~Irred,true); else Append(~Irred,false); end if; end for; elif fact[1][2] eq 1 then Irred[Index(Decomp,V)] := true; end if; end for; p := NextPrime(p); end while; X`decomp := Decomp; end if; return X`decomp; end intrinsic; intrinsic TestDisc(p::RngIntElt) -> BoolElt {Returns false if p does NOT divide the discriminant of the subalgeba of Hecke generated by T2, T3, T5, and T7.} assert IsPrime(p); X :=Mestre(p); for q in [2,3,5,7] do T := MatrixAlgebra(GF(p),Degree(X))!TpD(X,q); if Discriminant(CharacteristicPolynomial(T)) ne 0 then return false; end if; end for; return true; end intrinsic; /****************************************************************** * END METHOD OF GRAPHS * ******************************************************************/ /****************************************************************** * Kaneko and Zagier * * "Supersingular j-invariants, hypergeometric series, * * and Atkin's orthogonal polynomials.", Conference in Honor of * * Oliver Atkin. 1998, AMS/IP Studies in Adv. Math, Vol 7. * ******************************************************************/ /////////////////////////////////////////////////////////////////// // Atkin's polynomial. From Theorem 4 of Kaneko and // /////////////////////////////////////////////////////////////////// intrinsic AtkinPolynomial(n::RngIntElt, p::RngIntElt) -> RngUPolElt {Returns the nth Atkin polynomial modulo p > 2n.} assert p gt 2*n; K := GF(p); P := PolynomialRing(Rationals()); R := PolynomialRing(K); Q := PolynomialRing(K); // Compute the binomial poly's "choose m" for m <= n. // a choose i <---> Evaluate(B[i+1],a) B := [R!BinomialPolynomial(a,i) : i in [0..n]]; return &+[ 12^(3*i)* &+[ (-1)^m * Evaluate(B[i-m+1],K!-1/12) *Evaluate(B[i-m+1],K!-5/12) *Evaluate(B[m+1],K!(n+1/12)) *Evaluate(B[m+1],K!(n-7/12)) /Evaluate(B[m+1],K!(2*n-1)) : m in [0..i] ]*j^(n-i) : i in [0..n] // summing over these i ]; end intrinsic; /////////////////////////////////////////////////////////////////// // Continuing with Kaneko and Zagier: // // Theorem 3: Let n_p = genus(X_0(p))+1 be the degree // // of the pth supersingular polynomial ss_p. Then // // ss_p = n_pth Atkin polynomial (mod p). // /////////////////////////////////////////////////////////////////// intrinsic SupersingularPolynomial(p::RngIntElt) -> RngUPolElt {Returns the pth supersingular polynomial. The roots are the supersingular j-invariants in characteristic p.} return AtkinPolynomial(g0(p)+1,p); end intrinsic; intrinsic SupersingularBasis(p::RngIntElt) -> SeqEnum {Returns a list of the supersingular j-invariants in characteristic p. However, it is *much* faster to use the command SupersingularBasis(Mestre(p)).} Fx := PolynomialRing(GF(p^2)); return [-Evaluate(a[1],0) : a in Factorization(Fx!SupersingularPolynomial(p))]; end intrinsic; /************* END Kaneko and Zagier ****************************/ /****************************************************************** * SERRE MOD PQ * * Joint work with Barry Mazur * ******************************************************************/ intrinsic RibetRaise (A::ModTupFld, ell::RngIntElt) -> RngIntElt {Given a new factor A=A_f and a prime ell, returns Norm((ell+1-a_ell(f))*(ell+1+a_ell(f))). This is an integer divisible by the congruence primes relating A to a form at level N(A)*ell.} f := qEigenform(A,ell+1); a := Coefficient(f,ell); return Norm((ell+1-a)*(ell+1+a)); end intrinsic; function SerreModPQHelper( HN, f, HM, g, HNM) N := Level(HN); M := Level(HM); aM := Coefficient(qEigenform(f,M+1),M); bN := Coefficient(qEigenform(g,N+1),N); plist := Factorization(Norm((1+M-aM)*(1+M+aM))); plist := [p[1] : p in plist | Numerator((N-1)/12) mod p[1] ne 0]; qlist := Factorization(Norm((1+N-bN)*(1+N+bN))); qlist := [q[1] : q in qlist | Numerator((M-1)/12) mod q[1] ne 0]; // locate f and g inside HNM: D := FullDecomposition(HNM); for i in [1..#D] do if not assigned fnew and not IsNew(D[i]) and Level(OldFactor(D[i])) eq N and IsogenyClass(OldFactor(D[i])) eq IsogenyClass(f) then fnew := D[i]; end if; if not assigned gnew and not IsNew(D[i]) and Level(OldFactor(D[i])) eq M and IsogenyClass(OldFactor(D[i])) eq IsogenyClass(g) then gnew := D[i]; end if; end for; // find congruences at level NM: n := NumNewFactors(D); fcon := [Conductor(&*([1] cat AbelianIntersection(fnew,Factor(D,i)))) : i in [1..n]]; gcon := [Conductor(&*([1] cat AbelianIntersection(gnew,Factor(D,i)))) : i in [1..n]]; ans := [<1,0,fcon>, <0,1,gcon>]; for p in plist do for q in qlist do l:=[]; for i in [1..n] do if fcon[i] mod p eq 0 and gcon[i] mod q eq 0 then Append(~l,i); end if; end for; Append(~ans,); end for; end for; return ans; end function; intrinsic SerreModPQ(N::RngIntElt, M::RngIntElt) -> SeqEnum {} if (N lt 11 or N eq 13) or (M lt 11 or M eq 13) then return []; end if; // if HeckeVerbose then "SerreModPQ: Creating and decomposing modular symbols spaces."; s := Cputime(); // end if; HN := ModularSymbols(N,2,0); DN := FullDecomposition(HN); if NumNewformFactors(DN) eq 0 then return []; end if; HM := ModularSymbols(M,2,0); DM := FullDecomposition(HM); if NumNewformFactors(DM) eq 0 then return []; end if; HNM := ModularSymbols(N*M,2,0); DNM := FullDecomposition(HNM); "Time: ", Cputime(s); // if HeckeVerbose then "SerreModPQ: Considering pairs (f,g):"; // end if; ans := []>]; for i in [1..NumNewformFactors(DN)] do for j in [1..NumNewformFactors(DM)] do f := Factor(DN,i); g := Factor(DM,j); // if HeckeVerbose then "f = ", qEigenform(f,7); "g = ", qEigenform(g,7); // end if; Append(~ans,); end for; end for; return ans; end intrinsic; intrinsic SerreModPQTable( name::MonStgElt, NMstart::RngIntElt, NMstop::RngIntElt) {} out := Open(name, "w"); fprintf out, "// Serre's conjecture modulo pq.\n"; fprintf out, "// %o <= NM <= %o\n\n", NMstart, NMstop; fprintf out, "SERRE := [*\n"; delete out; s := [p : p in [11..NMstop] | IsPrime(p)]; for i in [1..#s] do for j in [i+1..#s] do if s[i]*s[j] ge NMstart and s[i]*s[j] le NMstop then cmd := "SerreModPQTableNM(\"" cat name cat "\"," cat IntegerToString(s[i]) cat "," cat IntegerToString(s[j]) cat "); quit; quit;"; go := "echo '" cat cmd cat "' | magma"; print "System(",go,")"; System(go); end if; end for; end for; out := Open(name, "a"); fprintf out, "[]\n*];\n"; delete out; end intrinsic; intrinsic SerreModPQTableNM(name::MonStgElt, N::RngIntElt, M::RngIntElt) {} out := Open(name, "a"); fprintf out, "%o,\n", SerreModPQ(N,M); delete out; end intrinsic; /************* END SERRE MOD PQ ***********************************/ /****************************************************************** * MAZUR's CONNECTEDNESS QUESTION * * "Is the Hecke algebra T_0(p,k), for k > 2, connected?" * * Instead of determining connectedness of the Hecke algebra, we * * determine connectedness of the subtori intersection graph. * ******************************************************************/ intrinsic IsIntersectionGraphConnected(N::RngIntElt, k::RngIntElt) -> BoolElt {Returns True if and only if the subtori intersection graph of J_k(N) is connected.} "Analyzing (",N,",",k,"):"; // 1. Construct modular symbols of weight k for Gamma_0(N). M := ModularSymbols(N,k); // 2. Find eigenspaces D := FullDecomposition(M); Print(D); // 3. Determine the connected component containing D[1]. /* Algorithm: a) Add vertex (1) (that is, D[1]) to the connected component. b) Add each vertex (i)=/=(1) that is adjacent to (1) to the connected component. c) Repeat step b) for each new vertex just added. If no new vertex is added, continue with step d). d) If the connected component equals the whole graph, then the graph is connected. */ if NumberOfCuspidalFactors(D) eq 0 then return true; end if; n := #D ; // number of vertices V := [Integers()|]; // vertices in the connected component of (1) W := [2..n]; A := [1]; // vertices just added while A ne [Integers()|] do B := [Integers()|]; // new added vertex list for i in A do for j in W do if IsCuspidal(D[j]) and #AbelianIntersection(D[i],D[j]) gt 0 then i,"<-->",j,": ",AbelianIntersection(D[i],D[j]); Append(~B,j); Remove(~W,Position(W,j)); end if; end for; end for; V := V cat A; A := B; "V = ",V; end while; return #V eq NumberOfCuspidalFactors(D); end intrinsic; /****************************************************************** * ARTIN's CONJECTURE * * Buzzard-Stein, "New examples of modular Icosahedral Galois * * representations." In preparation, 1999. * ******************************************************************/ intrinsic Artin(N::RngIntElt, degrees::SeqEnum, zeroprimes::SeqEnum) -> ModTupFld, ModTupFld { Given an integer N (the conductor), a sequence [d1,...,dr] of degrees which define a mod N character eps, and a sequence [p1,p2,...,pn] of primes p such that a_p=0, this function computes the subspace ker(Tp1') meet ker(Tp2') ... meet ker(Tpn') of the space M_5(Gamma_1(N),eps;F5bar)^+. (Primes denote transpose.) The first argument returned is the subspace W, and the second is the space of modular symbols. The subspace W is a subspace of the dual and can be computed efficiently using FastTp.} "Artin in weight 5"; "* Creating Dirichlet Character"; eps := DirichletCharacter(N, degrees, 5); "* Creating the space of modular symbols"; M := ModularSymbols(eps,5,+1); "Dimension = ", Dimension(M); "* Computing kernels"; R:=PolynomialRing(GF(5)); return DualSubspace(M, [ : p in zeroprimes]), M; end intrinsic; function ToLatex(x) case Type(x): when FldFinElt: if x eq 0 then return "0"; end if; F := Parent(x); if x in PrimeField(F) then s := Sprintf("%o", x); else s := Sprintf("%o^{%o}", a, Log(a,x)); end if; return s; else error "Type not programmed."; end case; end function; intrinsic TableEigenvalueList(filename::MonStgElt, ap::SeqEnum, colheight::RngIntElt) {} out := Open(filename,"w"); fprintf out, "\\\\begin{array}{|rl|}\\\\hline\n"; p := 2; for i in [1..#ap] do fprintf out, "%o&%o\\\\\\\\\n", p, ToLatex(ap[i]); if i mod colheight eq 0 or i eq #ap then fprintf out, "\\\\hline\\\\end{array}\n"; if i lt #ap then fprintf out, "\\\\begin{array}{|rl|}\\\\hline\n"; end if; end if; p := NextPrime(p); end for; delete out; end intrinsic; intrinsic ArtinOrderTwo(h::RngUPolElt, pstop::RngIntElt) -> SeqEnum {} degs := [ < p, Lcm( [Degree(f[1]) : f in Factorization(PolynomialRing(IntegerRing(p))!h)])> : p in [2..pstop] | IsPrime(p)]; return [a[1] : a in degs | a[2] eq 2]; end intrinsic; intrinsic Artin2(N::RngIntElt, degrees::SeqEnum, zeroprimes::SeqEnum) -> ModTupFld, ModTupFld { Given an integer N that is coprime to 5; a sequence [d1,...,dr] of degrees which define a mod N character eps; and a sequence [p1,p2,...,pn] of primes p such that a_p=0; this function computes the subspace ker(Tp1') meet ker(Tp2') ... meet ker(Tpn') of the space M_2(Gamma_1(5*N),eps_5*eps;F5bar)^+. (Primes denote transpose.) The first argument returned is the subspace W, and the second is the space of modular symbols. The subspace W is a subspace of the dual and can be computed efficiently using FastTp.} "Artin in weight 2"; "* Creating Dirichlet Character"; // insert degree 4 at 5 into the character. i := 0; if N mod 2 eq 0 then i +:= 1; end if; if N mod 8 eq 0 then i +:= 1; end if; if N mod 3 eq 0 then i +:= 1; end if; degrees := [degrees[j] : j in [1..i]] cat [4] cat [degrees[j] : j in [i+1..#degrees]]; "degrees = ",degrees; eps := DirichletCharacter(5*N, degrees, 5); "* Creating the space of weight 2 modular symbols"; M := ModularSymbols(eps,2,+1); "Dimension = ", Dimension(M); "* Computing kernels"; R:=PolynomialRing(GF(5)); return DualSubspace(M, [ : p in zeroprimes]), M; end intrinsic; /************* END Conjecture of Artin ********************************/ /******************************************************************** * EMPIRICAL EVIDENCE FOR THE BIRCH AND SWINNERTON-DYER CONJECTURES * * FOR MODULAR JACOBIANS OF GENUS 2 CURVES. * * Flynn, Leprevost, Schaefer, Stein, Stoll, Wetherell. * * * ********************************************************************/ EVIDENCE_CURVES := [ <"23", "23A">, <"29", "29A">, <"31","31A">, <"35","35B">, <"39", "39B">, <"63", "63B">, <"65,A", "65C">, <"65,B", "65B">, <"67", "67B">, <"73", "73B">, <"85", "85B">, <"87", "87A">, <"93", "93A">, <"103", "103A">, <"107", "107A">, <"115", "115B">, <"117,A", "117B">, <"117,B", "117C">, <"125,A", "125A">,<"125,B", "125B">, <"133,A", "133B">, <"133,B","133A">, <"135", "135D">, <"147", "147D">, <"161", "161B">, <"165", "165A">, <"167", "167A">, <"175", "175E">, <"177", "177A">, <"188", "188B">, <"189", "189E">, <"191", "191A">]; intrinsic EvidenceCurve(n::RngIntElt) -> ModTupFld {Returns the modular symbols factor corresponding to the nth curve from the paper "Empirical evidence..."} requirege n,1; assert n le #EVIDENCE_CURVES; return ModularFactor(EVIDENCE_CURVES[n][2]); end intrinsic; intrinsic EvidenceCurve(s::MonStgElt) -> ModTupFld {Returns the modular symbols factor corresponding to the curve from the paper "Empirical evidence..." labeled s. Here s is given by the notation in that paper. E.g., "29", or "65,A".} for C in EVIDENCE_CURVES do if C[1] eq s then return ModularFactor(C[2]); end if; end for; error "Curve not found."; end intrinsic; intrinsic EvidenceLabel(n::RngIntElt) -> ModTupFld {} return EVIDENCE_CURVES[n][1]; end intrinsic; /********************************************************************* * TABLE CONSTRUCTION FUNCTIONS * * * * Note, to avoid problems with garbage collection, some of these * * functions write their output to a file and spawn sub-magma * * processes to do various bits of the computation. * * * * Characteristic Polynomials: TableCharpoly * * Birch and Swinnerton-Dyer conjecture: TableBSD * *********************************************************************/ intrinsic TableEigenforms(name::MonStgElt, Nstart::RngIntElt, Nstop::RngIntElt, k::RngIntElt, prec::RngIntElt) {Create a table of eigneforms on Gamma_0(N) of weight k.} out := Open(name, "w"); fprintf out, "// Table of weight %o eigenforms on Gamma_0(N) for %o<=N<=%o.\n", k, Nstart,Nstop; fprintf out, "// \n",prec; fprintf out, "// Each ap(x) is in terms of a root of h(x).\n\n\n"; fprintf out,"R := PolynomialRing(Rationals());\n"; fprintf out,"Nstart := %o;\n", Nstart; fprintf out,"Nstop := %o;\n", Nstop; fprintf out,"f := [Parent([<1,x,[1],[x]>])|[] : i in [1..%o]];\n", Nstop; delete out; for N in [Nstart..Nstop] do /*********************************************************** We spawn a new magma subprocess to avoid problems with memory garbage collection. Further, if while processing this level the spawned magma process is killed or halts due to a bug then we simply spawn a process to compute the next level, instead of completely halting. ***********************************************************/ "\nSpawning level ", N, " subprocess."; cmd := "TableEigenforms(\"" cat name cat "\"," cat IntegerToString(N) cat "," cat IntegerToString(k) cat "," cat IntegerToString(prec) cat "); quit; quit;"; s := "echo '" cat cmd cat "' | magma"; print "System(",s,")"; System(s); end for; end intrinsic; intrinsic TableEigenforms(name::MonStgElt, N::RngIntElt, k::RngIntElt, prec::RngIntElt) {Append the level N data for the table of eigenforms to the file name.} out := Open(name, "a"); pdiv := [x[1] : x in Factorization(N)]; M := ModularSymbols(N,k,+1); D := NewDecomposition(M); num := NumNewformFactors(D); if num gt 0 then Print(D); fprintf out, "f[%o] := [\n", N; for i in [1..num] do A := Factor(D,i); f := qEigenform(A,prec); "f = ",f; F := Parent(Coefficient(f,1)); R := PolynomialRing(Rationals()); fprintf out, "<%o, %o, %o, [", i, [WqOn(A,p)[1,1] : p in pdiv], Type(F) eq FldRat select "x-1" else Modulus(F); for p in [a : a in [2..prec] | IsPrime(a)] do fprintf out, "%o", R!Coefficient(f,p); if p lt prec then fprintf out, ","; end if; end for; Flush(out); fprintf out, "]>"; if i lt num then fprintf out,","; else fprintf out,"\n];"; end if; fprintf out, "\n"; A := 0; end for; Flush(out); end if; delete out; end intrinsic; intrinsic TableCharpoly(name::MonStgElt, k::RngIntElt, N1::RngIntElt, N2::RngIntElt, p1::RngIntElt, p2::RngIntElt) {Create a table of factored characteristic polynomials of Tp on S_k(Gamma_0(N)) for N between N1 and N2 and p between p1 and p2.} assert IsPrime(p1) and IsPrime(p2); P := [p : p in [p1..p2] | IsPrime(p)]; out := Open(name, "a"); fprintf out, "//Table of characteristic polynomials of Hecke operators Tp on\n" cat"//S_%o(Gamma_0(N)) for N between %o and %o and p between %o " cat"and %o.\n\n", k,N1,N2,p1,p2; fprintf out,"charpoly%o := [*\n", k; for N in [N1..N2] do M := ModularSymbols(N,k,+1); if Dimension(M) eq 0 then continue; end if; S := Sk(M); if Dimension(S) gt 0 then fprintf out, "\n\n// N=%o\n [*", N; for p in P do fprintf out, "%o", Factor(ModularCharpoly(HeckeOperatorOn(M,S,p))); Flush(out); if p ne p2 then fprintf out, ",\n"; else fprintf out, "\n *]"; if N ne N2 then fprintf out, ","; end if; end if; fprintf out, "\t\t//p=%o\n", p; end for; end if; end for; fprintf out, "\n *];"; end intrinsic; intrinsic TableBSD(name::MonStgElt, N1::RngIntElt, N2::RngIntElt) {Create a table of odd parts of BSD data for N between N1 and N2.} out := Open(name, "w"); fprintf out, "// Table of ODD parts of BSD data for N between %o and %o.\n", N1,N2; fprintf out, "// [Level, Isogeny class, Dimension,Wq's, Order(0-oo), Torsion bound,\n"; fprintf out, "// L(1)/Omega, Modular Degree, cpbar...,ShaAnalytic].\n\n\n"; fprintf out," N1 := %o;\n", N1; fprintf out," N2 := %o;\n", N2; fprintf out,"BSD := [*\n"; delete out; for N in [n : n in [N1..N2]|IsPrime(n)] do /*********************************************************** We spawn a new magma subprocess to avoid problems with memory garbage collection. Further, if while processing this level the spawned magma process is killed or halts due to a bug then we simply spawn a process to compute the next level, instead of completely halting. ***********************************************************/ "\nSpawning level ", N, " subprocess."; cmd := "TableBSD(\"" cat name cat "\"," cat IntegerToString(N) cat "); quit; quit;"; s := "echo '" cat cmd cat "' | magma"; print "System(",s,")"; System(s); if N lt N2 then out := Open(name, "a"); fprintf out, ","; delete out; end if; end for; out := Open(name, "a"); fprintf out,"*];\n"; delete out; end intrinsic; intrinsic TableBSD(name::MonStgElt, N::RngIntElt) {Append the level N of the BSD table to the file "name".} out := Open(name, "a"); pdiv := [x[1] : x in Factorization(N)]; M := ModularSymbols(N,2,+1); D := Decomposition(M); num := NumNewformFactors(D); if num gt 0 then Print(D); fprintf out, "[ // J_0(%o)\n", N; for i in [1..num] do A := Factor(D,i); dat := [N, i, Dimension(A)] cat [WqOn(A,p)[1,1] : p in pdiv] cat [OddPart(CuspOrder(A)), OddPart(TorsionBound(A,DefaultPrec(A))), OddPart(LRatio(A)), OddPart(ModularDegree(A))] cat [ComponentGroup(A,p) : p in pdiv] cat [Sha(A)]; dat; fprintf out, "%o", dat; Flush(out); if i lt num then fprintf out, ","; end if; fprintf out, "\n"; A := 0; end for; fprintf out, " ] // end J_0(%o)\n\n", N; Flush(out); else fprintf out, "[ ] // J_0(%o) (no newforms)\n", N; end if; delete out; end intrinsic; intrinsic TableBSDFull(name::MonStgElt, N1::RngIntElt, N2::RngIntElt) {Create a table of odd parts of BSD data for N between N1 and N2.} out := Open(name, "w"); fprintf out, "// Table of BSD data for N between %o and %o.\n", N1,N2; fprintf out, "// [Level, Isogeny class, Dimension,Wq's, Order(0-oo), Torsion bound,\n"; fprintf out, "// L(1)/Omega, Modular Degree, cpbar...,ShaAnalytic].\n\n\n"; fprintf out," N1 := %o;\n", N1; fprintf out," N2 := %o;\n", N2; fprintf out,"BSD := [*\n"; delete out; for N in [n : n in [N1..N2]|IsPrime(n)] do "\nSpawning level ", N, " subprocess."; cmd := "TableBSDFull(\"" cat name cat "\"," cat IntegerToString(N) cat "); quit; quit;"; s := "echo '" cat cmd cat "' | magma"; print "System(",s,")"; System(s); if N lt N2 then out := Open(name, "a"); fprintf out, ","; delete out; end if; end for; out := Open(name, "a"); fprintf out,"*];\n"; delete out; end intrinsic; intrinsic TableBSDFull(name::MonStgElt, N::RngIntElt) {Append the level N of the BSD table to the file "name".} out := Open(name, "a"); pdiv := [x[1] : x in Factorization(N)]; M := ModularSymbols(N,2); D := Decomposition(M); num := NumNewformFactors(D); if num gt 0 then Print(D); fprintf out, "[ // J_0(%o)\n", N; for i in [1..num] do A := Factor(D,i); dat := [N, i, Dimension(A) div 2] cat [WqOn(A,p)[1,1] : p in pdiv] cat [CuspOrder(A), TorsionBound(A,DefaultPrec(A)), LRatio(A), // slower, but works at 2. ModularDegree(A)] cat [ComponentGroup(A,p) : p in pdiv] cat [Sha(A)]; dat; fprintf out, "%o", dat; Flush(out); if i lt num then fprintf out, ","; end if; fprintf out, "\n"; A := 0; end for; fprintf out, " ] // end J_0(%o)\n\n", N; Flush(out); else fprintf out, "[ ] // J_0(%o) (no newforms)\n", N; end if; delete out; end intrinsic; intrinsic TableLRatio(name::MonStgElt, N::RngIntElt, k1::RngIntElt, k2::RngIntElt) {Create a table of rational parts of special values of L-functions.} out := Open(name, "a"); fprintf out, "// Rational parts of special values of L-functions\n"; fprintf out, "// N = %o, %o <= k <= %o\n", N, k1, k2; fprintf out, "// [dimension, LRatio(1), LRatio(2), ... ].\n\n\n"; fprintf out,"LRatio := [*\n"; for k in [k1..k2] do M := ModularSymbols(N,k); if Dimension(M) eq 0 then fprintf out, "[] // S_%o(%o) (no newforms)\n", k, N; continue; end if; D := Decomposition(M,23); num := NumNewformFactors(D); if num gt 0 then Print(D); fprintf out, "[ // S_%o(%o)\n", k, N; for i in [1..num] do A := Factor(D,i); dat := [Dimension(A)] cat [LRatio(A,j) : j in [1..k-1]]; dat; fprintf out, "%o", dat; Flush(out); if i lt num then fprintf out, ","; end if; fprintf out, "\n"; A := 0; end for; fprintf out, "]%o // end S_%o(%o)\n\n", k lt k2 select "," else " ", k, N; Flush(out); else fprintf out, "[]%o // S_%o(%o) (no newforms)\n", k lt k2 select "," else " ", k, N; end if; end for; fprintf out, "\n *];"; end intrinsic; intrinsic TablePrimeSplittings(name::MonStgElt, P1::RngIntElt, P2::RngIntElt) {Create a table of splittings of J_0(p).} out := Open(name, "w"); fprintf out, "// Table of splittings of J_0(p) for p between %o and %o.\n", P1,P2; fprintf out, "// [Level, dimension 1, dimension 2, ...]\n"; fprintf out,"SPLITTING := [*\n"; delete out; for p in [p : p in [P1..P2]|IsPrime(p)] do /*********************************************************** We spawn a new magma subprocess to avoid problems with memory garbage collection. Further, if while processing this level the spawned magma process is killed or halts due to a bug then we simply spawn a process to compute the next level, instead of completely halting. ***********************************************************/ "\nSpawning level ", p, " subprocess."; cmd := "TablePrimeSplittings(\"" cat name cat "\"," cat IntegerToString(p) cat "); quit; quit;"; s := "echo '" cat cmd cat "' | magma"; print "System(",s,")"; System(s); if p lt P2 then out := Open(name, "a"); fprintf out, ","; delete out; end if; end for; out := Open(name, "a"); fprintf out,"*];\n"; delete out; end intrinsic; intrinsic TablePrimeSplittings(name::MonStgElt, p::RngIntElt) {Append the level p of the prime splitting table.} out := Open(name, "a"); X := Mestre(p); D := Decomposition(X); dims:= Sort(Dimension(D)); dims[1] := p; fprintf out, "%o\n", dims; delete out; end intrinsic; intrinsic TableBSD(name::MonStgElt, N1::RngIntElt, N2::RngIntElt) {Create a table of odd parts of BSD data for N between N1 and N2.} out := Open(name, "w"); fprintf out, "// Table of ODD parts of BSD data for N between %o and %o.\n", N1,N2; fprintf out, "// [Level, Isogeny class, Dimension,Wq's, Order(0-oo), Torsion bound,\n"; fprintf out, "// L(1)/Omega, Modular Degree, cpbar...,ShaAnalytic].\n\n\n"; fprintf out," N1 := %o;\n", N1; fprintf out," N2 := %o;\n", N2; fprintf out,"BSD := [*\n"; delete out; for N in [n : n in [N1..N2]|IsPrime(n)] do /*********************************************************** We spawn a new magma subprocess to avoid problems with memory garbage collection. Further, if while processing this level the spawned magma process is killed or halts due to a bug then we simply spawn a process to compute the next level, instead of completely halting. ***********************************************************/ "\nSpawning level ", N, " subprocess."; cmd := "TableBSD(\"" cat name cat "\"," cat IntegerToString(N) cat "); quit; quit;"; s := "echo '" cat cmd cat "' | magma"; print "System(",s,")"; System(s); if N lt N2 then out := Open(name, "a"); fprintf out, ","; delete out; end if; end for; out := Open(name, "a"); fprintf out,"*];\n"; delete out; end intrinsic; intrinsic TableCongruenceR(name::MonStgElt, N1::RngIntElt, N2::RngIntElt) {Create a table of congruence numbers r.} out := Open(name, "w"); fprintf out, "// Table of congruence numbers r.\n"; fprintf out, "// [Level, Isogeny class, r, structure].\n\n\n"; fprintf out,"N1 := %o;\n", N1; fprintf out,"N2 := %o;\n", N2; fprintf out,"CongruenceR := [*\n"; delete out; for N in [n : n in [N1..N2]] do "\nSpawning level ", N, " subprocess."; cmd := "TableCongruenceR(\"" cat name cat "\"," cat IntegerToString(N) cat "); quit; quit;"; s := "echo '" cat cmd cat "' | magma"; print "System(",s,")"; System(s); if N lt N2 then out := Open(name, "a"); fprintf out, ","; delete out; end if; end for; out := Open(name, "a"); fprintf out,"*];\n"; delete out; end intrinsic; intrinsic TableCongruenceR(name::MonStgElt, N::RngIntElt) {Append the level N of the congruence table to the file "name".} out := Open(name, "a"); M := ModularSymbols(N,2,+1); D := Decomposition(M); num := NumNewformFactors(D); if num gt 0 then Print(D); fprintf out, "[ // J_0(%o)\n", N; for i in [1..num] do A := Factor(D,i); r, struc := CongruenceR(A,Integers()!Round(2.5*HeckeBound(M))); dat := [N, i, r]; dat; fprintf out, "%o", dat; Flush(out); if i lt num then fprintf out, ","; end if; fprintf out, "\n"; A := 0; end for; fprintf out, " ] // end J_0(%o)\n\n", N; Flush(out); else fprintf out, "[ ] // J_0(%o) (no newforms)\n", N; end if; delete out; end intrinsic; intrinsic TableConnected(name::MonStgElt, N::RngIntElt, k::RngIntElt) {} out := Open(name, "a"); fprintf out, "C[%o,%o] := %o;\n", N, k, IsIntersectionGraphConnected(N,k) eq true select 1 else 0; delete out; end intrinsic; intrinsic TableConnected(name::MonStgElt, N1::RngIntElt, N2::RngIntElt, k1::RngIntElt, k2::RngIntElt) {Create a table of connectedness of abelian intersection graph in the "box" T_k(Gamma_0(N)), where N1 <= N <= N2 (N prime) and k1 <= k <= k2.} out := Open(name, "w"); fprintf out, "// Table of connected in the box %o <= N <= %o and %o <= k <= %o.\n", N1,N2, k1, k2; fprintf out, "// C[N,k] := either 1 (true) or 0 (false),\n\n"; fprintf out," N1 := %o;\n", N1; fprintf out," N2 := %o;\n", N2; fprintf out," k1 := %o;\n", k1; fprintf out," k2 := %o;\n", k2; fprintf out," C := RMatrixSpace(Integers(),%o,%o);\n\n", N2, k2; delete out; for N in [n : n in [N1..N2] | IsPrime(n)] do for k in [k : k in [k1..k2] | IsEven(k)] do "\nSpawning (", N, ",",k,") subprocess."; cmd := "TableConnected(\"" cat name cat "\"," cat IntegerToString(N) cat "," cat IntegerToString(k) cat "); quit; quit;"; s := "echo '" cat cmd cat "' | magma"; print "System(",s,")"; System(s); end for; end for; end intrinsic; /************* END Tables ****************************************/ /************************************************************* * * * BASIC ARITHMETIC (not specific to modular symbols * * * *************************************************************/ intrinsic Conductor(N::RngIntElt) -> RngIntElt {} return &*[Integers()|p[1] : p in Factorization(N)]; end intrinsic; intrinsic OddPart(n::.) -> . {Returns the odd part of n.} return n eq 0 select 0 else n / 2^Valuation(n,2); end intrinsic; intrinsic DotProd(a::., b::.) -> . {} return &+[a[i]*b[i] : i in [1..#Eltseq(a)]]; end intrinsic; intrinsic Volume(L::SeqEnum) -> RngIntElt {Compute the absolute value of the determinant of the matrix got by finding a Z-lattice basis for the Z-span of the rows in the sequence L.} n := Degree(L[1]); B := Basis(Lattice(RMatrixSpace(Rationals(), #L, n)!L)); VB := [VectorSpace(Rationals(),n)!b : b in B]; r := #B; B := &cat[Eltseq(b) : b in B]; if r ne n then return 0; end if; Mat := MatrixAlgebra(Rationals(), n)!B; return Abs(Determinant(Mat)); end intrinsic; intrinsic Norm(x::RngUPolResElt) -> . {Compute the norm of an element of a quotient of a polynomial ring.} return Resultant(PreimageRing(Parent(x))!x,Modulus(Parent(x))); end intrinsic; intrinsic ToNth(n::RngIntElt) -> MonStgElt {Returns the word for nth.} requirege n,0; case n: when 1: return "1st"; when 2: return "2nd"; when 3: return "3rd"; else return IntegerToString(n) cat "th"; end case; end intrinsic; intrinsic PrimeSeq(p1::RngIntElt, p2::RngIntElt) -> SeqEnum {The sequence of primes [p1 ... p2].} return [p : p in [p1..p2] | IsPrime(p)]; end intrinsic; /***************** END BASIC ARITHMETIC ************************/ /***************** BASE CHANGE *********************************/ intrinsic BaseExtendCharacter(eps::Rec, F::.) -> Rec {Base extend the character eps to the field F.} eps`F := F; if assigned eps`ims then eps`ims:=[F!eps`ims[i] : i in [1..#eps`ims]]; end if; return eps; end intrinsic; intrinsic BaseExtendMlist(mlist::Rec, F::.) -> Rec { Base extend the mlist to the field F.} mlist`F := F; mlist`R := PolynomialRing(F,2); return mlist; end intrinsic; intrinsic BaseExtendQuot(quot::Rec, F::.) -> Rec { Base extend the quot to the field F.} quot`Scoef:=[F!quot`Scoef[i] : i in [1..#quot`Scoef]]; if #quot`Tquot gt 1 then V := VectorSpace(F,Degree(quot`Tquot[1])); quot`Tquot := [V!quot`Tquot[i] : i in [1..#quot`Tquot]]; end if; return quot; end intrinsic; /* // not needed yet. intrinsic BaseExtendVectorSpaceWithBasis(V::ModTupFld, F::.) -> ModTupFld {Base extend the VectorSpaceWithBasis V to the field F.} W := VectorSpace(F,Degree(V)); return VectorSpaceWithBasis([W|b : b in Basis(V)]); end intrinsic; */ intrinsic BaseExtendDecomposition(D::SeqEnum, F::.) -> SeqEnum {Base extend the decomposition D to the field F.} error "BaseExtendDecomposition: not yet programmed, though not hard!"; /* return [BaseExtendVectorSpaceWithBasis(D[i],F) : i in [1..#D]];*/ end intrinsic; intrinsic BaseExtend(M::ModTupFld, F::.) -> ModTupFld {Given a space M of modular symbols and a field F, attempt to compute M tensor F.} if Type(M) eq Type(F) and BaseField(M) eq F then return M; end if; V := VectorSpace(F,Dimension(M)); N := VectorSpaceWithBasis(Basis(V)); // N = M tensor F N`theCategory := M`theCategory; N`k := M`k; N`N := M`N; N`F := F; N`eps := BaseExtendCharacter(M`eps, F); N`sign := M`sign; if assigned M`mlist then N`mlist := BaseExtendMlist(M`mlist, F); end if; if assigned M`quot then N`quot := BaseExtendQuot(M`quot, F); end if; return N; end intrinsic; /*************RIP**************************************************/ // BUGS // BUGS // BUGS // BUGS // BUGS // BUGS // BUGS // BUGS /* [ ] 1/24/00. Error when winding submodule is trivial. > WindingSubmodule(M,11); [0 0 0] [0 0 0] WindingSubmodule( M: Full Vector space of degree 3 over Rational Field, i: 11 ) HeckeSpan( M: Full Vector space of degree 3 over Rational Field, v: (0 0 0) ) VectorSpaceZBasis( B: [ (0 0 0), (0 0 0) ] ) In file "/home/was/magma/hecke/hecke.m", line 1174, column 28: >> Latbase := Basis(Lattice(X)); ^ Runtime error in 'Lattice': Dimension of lattice must be non-zero I had to fix this in two places; in VectorSpaceZBasis I broke things apart more, and in the intrinsic PeriodImage(W,A) I added if Dimension(W) eq 0 then return VectorSpaceWithBasis([A|]); end if; [x] Old forms are ordered incorrectly, e.g., S_2(Gamma_0(74)) ANS: It seems fine; I must have already fixed this bug. [x] Period lattice computation still not good enough: e.g., > EE:=EllipticCurve(ModularFactor("27A")); WARNING: Not finding enough symbols! Switching to slow method. WARNING: Not finding enough symbols! Could not find enough invariant symbols. Need a better algorithm. ANS: I tried to find this bug later (01/01/00), but hecke worked fine. [ ] Bizarre overflow in real volumes when the dimension of the abelian variety is large. -- or it could just be that this is a hard computation... > A:=ModularFactor("41k4B"); Creating M_4(Gamma_1(41),eps;F_0), 0.56 seconds. Sorting and labeling factors... > DimensionAbelianVariety(A); 7 > RealVolume(A); (Using at least 54 terms of q-expansions.) 0 > RealVolume(A,60); 7.9507713280 E143 + 1.7868374241 E141*i [ ] Magma-bug?? : when computing LRatios and the level is VERY large so that the numbers are super-huge, the real volume is computed as zero, which is of course impossible. The algorithm never exhibits this problem at lower level so I suspect a Magma bug. Here is an example: System( echo 'TableBSD("table/bsd2501-2800",2503); quit; quit;' | magma ) Magma V2.5-1 Tue Aug 10 1999 06:33:33 on yomomma [Seed = 1944814577] Type ? for help. Type -D to quit. Loading startup file "/home/was/magma/Startup/init.m" Creating M_2(Gamma_1(2503),eps;F_0)^+, 49.05 seconds. Decomposition of space of modular symbols. Sorting and labeling factors... Modular symbols factor of dimension 94, level 2503, weight 2 (2503A) Modular symbols factor of dimension 114, level 2503, weight 2 (2503B) Modular symbols factor of dimension 1, level 2503, weight 2 (2503C) (Using at least 300 terms of q-expansions.) warning: currently there is a bug -- off by powers of 2. warning: currently there is a bug -- off by powers of 2. [ 2503, 1, 94, 1, 1, 1, 0, 1, 1, 0 ] (Using at least 300 terms of q-expansions.) warning: currently there is a bug -- off by powers of 2. TableBSD( name: table/bsd2501-2800, N: 2503 ) LRatio( A: Vector space of degree 209, dimension 114 over Rational Fiel... ) LRatio( A: Vector space of degree 209, dimension 114 over Rational Fiel..., s: 1 ) In file "/home/was/magma/Hecke/hecke.m", line 3095, column 40: >> A`LRatio[s] := Nrm * A`ZxZalp / VolL; ^ Runtime error in '/': Division by zero Total time: 44387.160 seconds */ /******************* SOLVED BUGS!!! :) :) ******************************** [X] ModularFActor("7k4B") didn't work, as pointed out in Kevin's tutorial. Fix: Changed ModularFactor(ModTupFld, RngIntElt): D := Decomposition(M); if iso gt #D then D := FullDecomposition(M); end if; return ModularFactor(D, iso); [] NormalizePair was returning s^(-1) instead of s causing the machine to compute M_k(N,eps^(-1)). This affected computation of the boundary map. [ ] Fix up degeneracy maps. Make newforms be correctly identified in all cases. [ ] AUG 15: solved! Problem: cusp definition wasn't quite right. I was using eps^(-1) instead of eps because I changed conventions. N=1825, chi=[4,3], mod 5 and 37 the dimension of Sk given by the program is 733, but it should (trace formula) be 734. time M:=ModularSymbols(DirichletCharacter(1825,[4,3],5),5,+1); time S:=Sk(M); print "dim1=",Dimension(S); Creating M_5(Gamma_1(1825),eps;F_25)... Created M in 67.861 seconds. Time: 67.959 Time: 362.449 dim1= 733 [ ] AUG 11: > L:=ssinvs(32,8); Creating M_8(Gamma_1(32),eps;F_0), 4.859 seconds. Sorting and labeling factors... (Warning: Imaginary comp group formula not yet proved.) 32k8A c4= -1010566113763763.74609201240896538820 - 9.9475982971 E-14*i c6= 5052082.65853537789483562182082029292 + 0.000000095035821023233019194*i j= 1728.0000000000000000000000 + 0.E-25*i o= 6.07084028795 E82 + 3.03542014397 E82*i i= 0.000000007450580596 + 0.003783406689*i > A:=L[1]; > RealVolume(A,150); 1.21416805731 E83 + 6.07084028795 E82*i > ImaginaryVolume(A,150); 0.000000007450580596 + 0.003783406689*i The real and imaginary volumes are COMPLETELY wrong. This is undoubtedly the result of an overflow in Magma somewhere. However, the period lattice seems right: > PeriodLattice(A); [ (-0.00173111669423147007725483719953*i), (0.00086555834711573503862740511631 - 0.00259667504134720511588225579930*i) ] This must be a bug because the lattice must be invariant under multiplication by i. This is probably a bug in Magma arising because certain numbers get too large. :( ANSWER: No! it is was my own carelessness, not a bug in Magma! In the function "function ellinv_NormalizePair(w1, w2)" the line if Abs(tau) lt 1 then should be if Abs(tau) ge 1 then *******************END SOLVED BUGS ********************************/