Sharedwww / Tables / hecke.mOpen in CoCalc
/****-*-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 MtheCategory;
end function;

procedure SetCategory(M, category)
MtheCategory := 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<Integers()|1>;
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=<X^0Y^(k-2),x_1>,  2=<0, x_2>, ..., n=<0,x_n>,                   //
// n+1=<X^1Y^(k-3),x_1>,n+2=<1, x_2>, ...,2n=<1,x_n>,                   //
// ...                                                                  //
//////////////////////////////////////////////////////////////////////////
function ManinSymbolList(k,N,F)
p1list := EnumerateP1(N);
n      := (k-1)*#p1list;
return rec<CManSymList |
k      := k,
F      := F,
R      := PolynomialRing(F,2),
p1list := p1list,
n      := n
>;
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 := mlistp1list;
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*#(mlistp1list) + ind;
end function;

intrinsic ManinSymbolApply(g::SeqEnum, i::RngIntElt,
mlist::Rec, eps::Rec) -> SeqEnum
{Apply g to the ith Manin symbol.}
k := mlistk;
uv, w, ind := UnwindManinSymbol(i,mlist);
p1list := mlistp1list;
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 [<Evaluate(eps,scalar),act_uv>];
/*
if Evaluate(eps,scalar) ne 0 then
return [<Evaluate(eps,scalar)^(-1),act_uv>];
else
return [<0,act_uv>];
end if;
*/
else
return [<1,act_uv>];
end if;
end if;

// Polynomial part.
R<x> :=PolynomialRing(mlistF);    // 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 := [<pol[w+1],  w*n + act_uv> : 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     //
//  <P,x>, where P is a polynomial in X and Y (an element      //
//  of the ring MmlistR, 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
<P,x>, where P is a polynomial in X and Y (an element
of the ring MmlistR, 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  := MmlistR;
h  := hom <R -> 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 <hP, gx>;
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
<P,x>, where P is a polynomial in X and Y (an element
of the ring MmlistR, 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
[[<c_1, i_1>, <c_2, i_2>, ... <c_r,i_r>], ... ]
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<R|Rels>;
gens := {@ i : i in [1..n] | not exists
{ g : g in Rels | LeadingMonomial(g) eq R.i } @};
Z    := VectorSpace(K,#gens);
quot := [ &+[Z|
: 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 := mlistn;
K := mlistF;
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 := mlistn;
F := mlistF;
k := mlistk;
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
t  := ManinSymbolApply(T,j,mlist,eps)[1];
tt := ManinSymbolApply(TT,j,mlist,eps)[1];
Append(~rels,  [<Scoef[j],Squot[j]>,
<t[1]*Scoef[t[2]],Squot[t[2]]>,
<tt[1]*Scoef[tt[2]],Squot[tt[2]]>]);
end if;
end for;
else
rels := [&cat[
[<Scoef[i],Squot[i]>],
[<x[1]*Scoef[x[2]],Squot[x[2]]>
: x in ManinSymbolApply(T,i,mlist,eps)],
[<x[1]*Scoef[x[2]],Squot[x[2]]>
: 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);
Mk    := k;
MN    := N;
Meps  := eps;
Msign := sign;
MF    := F;
Mmlist:= mlist;
Mquot := rec<CQuotient |
Sgens := Sgens,
Squot := Squot,
Scoef := Scoef,
Tgens := Tgens,
Tquot := Tquot>;
// MquotTquot 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!
MquotTquot := [V!v : v in MquotTquot];  // 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));
Mk := k;
MN := N;
Meps := eps;
Msign := sign;
MF := F;
Mquot := rec<CQuotient | Sgens:=[], Squot:=[],
Scoef:=[], Tgens := [], Tquot:=[]>;
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 Mk then return Weight(ModSymParent(M)); end if;
return Mk;
end intrinsic;

intrinsic Level(M::ModTupFld) -> RngIntElt
{Return the level of the space M of modular symbols.}
return MN;
end intrinsic;

intrinsic Conductor(M::ModTupFld) -> RngIntElt
{Return the level of the space M of modular symbols.}
return MN;
end intrinsic;

intrinsic BaseField(M::ModTupFld) -> .
{Return the base field of the space M of modular symbols.}
if not assigned MF then
return BaseField(ModSymParent(M));
end if;
return MF;
end intrinsic;

intrinsic Character(M::ModTupFld) -> SeqEnum
{Return the Dirichlet character of the space M of modular symbols.}
case GetCategory(M):
when ModSym:
return Meps;
when ModSymFac:
return MMeps;
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 Msign then return MMsign eq 1; end if;
return Msign 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 Msign then return MMsign eq -1; end if;
return Msign eq -1;
end intrinsic;

intrinsic Sign(M::ModTupFld) -> RngIntElt
{Returns the sign of M.}
if not assigned Msign then return MMsign; end if;
return Msign;
end intrinsic;

intrinsic IsogenyClass(A::ModTupFld) -> RngIntElt
{Returns the isogeny class of A.}
if not assigned Aiso then
return 0;
end if;
return Aiso;
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;
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 <alp, i> 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 := MquotTquot;
V:=Parent(Tquot[1]);
Sgens := MquotSgens;
Scoef := MquotScoef;
Squot := MquotSquot;
// 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 := MquotScoef;
Tquot := MquotTquot;
V:=Parent(Tquot[1]);
Squot := MquotSquot;
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     := MF;
mlist := Mmlist;
eps   := Meps;

quot  := Mquot;
Sgens := quotSgens;
Squot := quotSquot;
Tgens := quotTgens;
Tquot := quotTquot;

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     := MF;
k     := Mk;
mlist := Mmlist;
eps   := Meps;
quot  := Mquot;
Sgens := quotSgens;
Tgens := quotTgens;
Tquot := quotTquot;

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, <v[i],i>);
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(Mk, MN);
if not IsTrivial(Meps) 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<M|V>);
assert #B eq n;
// Find pivot columns.
e := Pivots(B);

// Step 2: Compute the matrix V*E of inner products.
VE    := RMatrixSpace(MF,n,n)![V[i][e[j]] : i in [1..n], j in [1..n]];
VEinv := VE^(-1);

return rec<CFastData| V:=V, e:=e, VEinv:=VEinv>;
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     := MF;
V     := FastDataV;
n     := #V;
m     := Dimension(M);
e     := FastDatae;
VEinv := FastDataVEinv;
// 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(MF,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 := MN;
k := Mk;
assert q ne 0;
assert N mod q eq 0;
assert k mod 2 eq 0;
assert IsTrivial(Meps);
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(MF);
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 Mcusplist then
Mcusplist := [];
end if;

list  := Mcusplist;
F     := MF;
eps   := Meps;
N     := MN;
k     := Mk;

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 <F!0,1>;
end if;
if IsTrivial(Character(M)) then
return <1,i>;
else
return <Evaluate(eps,alp)^(-1),i>;
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 <F!0,1>;
end if;
if IsTrivial(Character(M)) then
return <sign,i>;
else
return <sign*Evaluate(eps,alp),i>;
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,<a,c>);
Mcusplist := 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
Mcusplist[i][2] := 0;
c := F!0;
end if;
end if;

return <c,i>;
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 MBoundaryMap then
if Dimension(M) eq 0 then
MBoundaryMap := RMatrixSpace(MF,0,0)!0;
return MBoundaryMap;
end if;
Tgens := MquotTgens;
Sgens := MquotSgens;
F     := MF;
n     := #Tgens;
D     := [ [] : i in [1..n]];
for j in [1..n] do
i := Sgens[Tgens[j]];
uv, w, ind := UnwindManinSymbol(i,Mmlist);
g := LiftToCosetRep(uv, MN);
if w eq Mk-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
MBoundaryMap := RMatrixSpace(MF,Dimension(M),0)!0;
return MBoundaryMap;
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]];
MBoundaryMap := RMatrixSpace(MF,n,m) ! Drows;
end if;
return MBoundaryMap;
end intrinsic;

intrinsic CuspidalSymbols(M::ModTupFld) -> ModTupFld
{Return the subspace of S_k(N,Q) of cuspidal modular symbols.}
if not assigned MSkF then
MSkF := Kernel(BoundaryMap(M));
end if;
return MSkF;
end intrinsic;

intrinsic CuspidalSymbols(M::ModTupFld) -> ModTupFld
{Return the subspace of S_k(N,Q) of cuspidal modular symbols.}
if not assigned MSkF then
D := BoundaryMap(M);
if HeckeSuperverbose then
"Computing kernel of boundary map: ";
t := Cputime();
end if;
MSkF := Kernel(D);
if HeckeSuperverbose then
Cputime(t), " seconds.";
end if;
end if;
return MSkF;
end intrinsic;

intrinsic Sk(M::ModTupFld) -> ModTupFld
{}
return CuspidalSymbols(M);
end intrinsic;

intrinsic IntegralCuspidalSymbols(M::ModTupFld) -> ModTupFld
{}
if not assigned MSkZ then
Z     := IntegralModularSymbols(M);
D     := BoundaryMap(M);
DofZ  := [v*D : v in Basis(Z)];
DZ    := RMatrixSpace(MF,Nrows(D),Ncols(D)) ! DofZ;
MSkZ := VectorSpaceWithBasis(
LinearCombinations(Basis(IntegerKernel(DZ)), Basis(Z))
);
end if;
return MSkZ;
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 MSkZplus then
S    := SkZ(M);
star := Restrict(StarInvolution(M),S);
MSkZplus := IntegerKernelOn(star-1, S);
end if;
return MSkZplus;
end intrinsic;

intrinsic SkZMinus(M::ModTupFld) -> ModTupFld
{}
if not assigned MSkZminus then
S    := SkZ(M);
star := Restrict(StarInvolution(M),S);
MSkZminus := IntegerKernelOn(star+1, S);
end if;
return MSkZminus;
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 M1F eq M2F;
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 Mold then
Mold := 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(Mold[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.
Mold[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
Mold[pos] := BaseExtend(TrivialSpace(Weight(M),eps,Sign(M)),F);
else
epsNN := AssociatedModMCharacter(eps,NN);
Mold[pos] := BaseExtend(
ModularSymbols(epsNN,Weight(M),Sign(M)),F);
end if;
(Mold[pos])M := M;  // tag this as *our* child.
end if;
end if;
return Mold[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 Msknew then
Msknew := NewSubspace(M) meet CuspidalSymbols(M);
end if;
return Msknew;
end intrinsic;

intrinsic NewSubspace(M::ModTupFld) -> ModTupFld
{Returns the p-new subspace of M.}
if not assigned Mmknew then
if Level(M) eq 1 then
Mmknew := M;
else
Mmknew := &meet[pNewSubspace(M,p[1])
: p in Factorization(Level(M))];
end if;
end if;
return Mmknew;
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 Msknewdual then
Msknewdual := NewDualSubspace(M);
// cut out the cuspidal part using Hecke operators.
S := Sknew(M);
if Dimension(S) eq 0 then
Msknewdual := S;
else
p := 2;
while Dimension(Msknewdual) gt Dimension(S) do
T := HeckeOperator(M,p);
f := ModularCharpoly(Restrict(T,S));
Ton := Restrict(Evaluate(f,Transpose(T)),Msknewdual);
Msknewdual := KernelOn(Ton,Msknewdual);
p := NextPrime(p);
end while;
end if;
end if;
return Msknewdual;
end intrinsic;

intrinsic NewDualSubspace(M::ModTupFld) -> ModTupFld
{Returns the p-new subspace of M dual.}
if not assigned Mmknewdual then
Mmknewdual := &meet ([pNewDualSubspace(M,p[1])
: p in Factorization(Level(M))] cat
[M]);
end if;
return Mmknewdual;
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<X,Y> := 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 Mk-1;
R<X,Y> := PolynomialRing(MF,2);
i := i-1;
return ConvFromModularSymbol(M,<X^i*Y^(Mk-2-i),[[0,1],[1,0]]>);
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 Mk-1;
if not assigned Mwinding then
Mwinding := Seqlist([0 : i in [1..Mk-1]]);
end if;
if Type(Mwinding[i]) eq RngIntElt then
Mwinding[i] := HeckeSpan(M,WindingElement(M,i)) ;
end if;
return Mwinding[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
**********************************************************/
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 := Mmlist;
uv, w:= UnwindManinSymbol (
MquotSgens[MquotTgens[i]], mlist);
R := mlistR;
return <R.1^w*R.2^(Mk-2-w),uv>;
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 [<v[i]*x[1],x[2]> : i in [1..Dimension(M)]
| v[i] ne 0
where x is ConvToManinSymbol(M,i) ];
end intrinsic;

function ConvToModSym(M, i)
mlist := Mmlist;
uv, w, ind   := UnwindManinSymbol (
MquotSgens[MquotTgens[i]], mlist);
g     := LiftToCosetRep(uv, MN);
k := Mk;
R := mlistR;
if k eq 2 then
return <R!1,[[g[2],g[4]], [g[1],g[3]]]>;   // {g(0), g(oo)}
end if;
h  := hom <R -> 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 <hP, [[g[2],g[4]], [g[1],g[3]]]>;
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 [ <v[i]*x[1],x[2]> : 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 MModSymBasis then
MModSymBasis :=
[ConvToModularSymbol(M, M.i) : i in [1..Dimension(M)]];
end if;
return MModSymBasis;
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   := MmlistR;
return ConvFromManinSymbol(M, R.2^(Mk-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   := MmlistR;
P   := R!P;
k   := Mk;
n   := #(Mmlistp1list);
ind,s := P1Reduce(uv, Mmlistp1list);
if not IsTrivial(Meps) then
P *:= Evaluate(Meps, s);
end if;
P := R!P;
ans := [<MonomialCoefficient(P,R.1^w*R.2^(k-2-w)),
w*n + ind> : 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 <P,x>, 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<Integers()|MN>;
//   ZN:=IntegerRing(MN);
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  := MmlistR;
h  := hom <R -> 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,<MmlistR!1,[alpha,beta]>);
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 <P,x>, 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 <P,x>, 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 MZ 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    := Mquot;
gens := [ QScoef[i]*QTquot[QSquot[i]] : i in [1..#QSquot] ];
n    := #QTgens;
gens := [g[i] :  i in [1..n], g in gens];
B    := [M!v : v in Basis(Lattice(n,gens))];
MZ  := VectorSpaceWithBasis(B);
end if;
return MZ;
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;
// 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<V|V!0>;
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 AM;
end intrinsic;

intrinsic ModSymParent(A::ModTupFld) -> ModTupFld
{Returns the parent of the modular symbols factor A.}
return AM;
end intrinsic;

intrinsic IsNew (V::ModTupFld) -> BoolElt
{}
return Visnew;
end intrinsic;

intrinsic OldFactor (V::ModTupFld) -> ModTupFld
{}
return Vold;
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 Msemidecomp then
Msemidecomp := SemiDecomposition(M,M,2,13);
end if;
return Msemidecomp;
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(MN, 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.
Wisnew := true;
WN     := MN;
WM     := M;
Append(~D,W);
elif p gt stop then
if (IsPrime(MN) and Mk eq 1) or MN eq 1 then
// must be new anyways.
WN     := MN;
Wisnew := true;
else  // probably isn't new.
Wisnew := false;    // WARNING: in fact it could be TRUE here!!!!
end if;
WM     := M;
WN     := 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;
WN := MN;
WM := M;
Wiscusp := 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);
Wisnew := true;
WN := MN;
WM := M;
Wiscusp := W subset Sknewdual(M);
appended := false;
if (a eq 1 or
(Sign(M) eq 0 and a eq 2 and Wiscusp)) then  // irred.
Append(~D,