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 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<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 := 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 [<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(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 := [<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 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
  <P,x>, 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 -> 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 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 
    [[<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|
        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,  [<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);
   M`k    := k;
   M`N    := N;
   M`eps  := eps;
   M`sign := sign;
   M`F    := F;
   M`mlist:= mlist;
   M`quot := rec<CQuotient |  
                 Sgens := Sgens, 
                 Squot := Squot, 
                 Scoef := Scoef,
                 Tgens := Tgens, 
                 Tquot := Tquot>;  
   // 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<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 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 <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 := 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, <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(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<M|V>);
   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<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     := 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 <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>);
   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 <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 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<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 M`k-1;
   R<X,Y> := PolynomialRing(M`F,2);
   i := i-1;
   return ConvFromModularSymbol(M,<X^i*Y^(M`k-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 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 <R.1^w*R.2^(M`k-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 := 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 <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 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 := [<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()|M`N>;
//   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 -> 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,<M`mlist`R!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 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<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 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,