CoCalc Public Fileswww / tables / hecke.m
Author: William A. Stein
1/****-*-magma-*-*************************************************************
2 *                 HECKE:  Modular Symbols Calculator                      *
3 *                                                                         *
4 *                      Version 0.5M  (February 9, 2000)                   *
5 *                                                                         *
6 *Send bug reports and suggestions to William Stein was@math.berkeley.edu. *
7 ***************************************************************************/
8
9/*************************************************************************
10  TODO  TODO  TODO  TODO  TODO  TODO  TODO  TODO  TODO
11
12
13   [ ] Intersections with Eisenstein series NOT programmed!
14
15   [ ] Add a feature for "tensoring with Qp".
16
17   [ ] Decomposition  of S_3(25,[4];Q) does not work.
18
19   [ ] Need to be able to compute eigenforms over a single
20       number field, so they can be worked with.
21
22> m4:=msym(DirichletCharacter(13,[3]),4,+1);
23> Print(Decomposition(m4));
24Modular symbols factors:
25 13k4A:  dim = 1    cuspidal
26 13k4B:  dim = 2    cuspidal
27 13k4C:  dim = 1  eisenstein
28 13k4D:  dim = 1  eisenstein
29> f4b:=qEigenform(Factor(m4,2),17);
30> MinimalPolynomial(Coefficient(f4b,2));
31x^2 + (-5*zeta - 5)*x + 2*zeta
32> f4b;
33q + a*q^2 + (-3*a + (5*zeta + 5))*q^3 + ((5*zeta + 5)*a - 10*zeta)*q^4 +
34    (-5*zeta*a - 20)*q^5 + ((-10*zeta - 10)*a + 6*zeta)*q^6 + ((zeta + 1)*a +
35    5*zeta)*q^7 + (7*zeta*a + 10)*q^8 + ((15*zeta + 15)*a - 20*zeta)*q^9 + (5*a
36    - 10*zeta - 10)*q^10 + (-15*a + (29*zeta + 29))*q^11 + (-20*zeta*a +
37    20)*q^12 + ((-7*zeta - 5)*a + (50*zeta + 45))*q^13 + (10*zeta*a + 2)*q^14 +
38    (10*a - 70*zeta - 70)*q^15 + (15*a - 66*zeta - 66)*q^16 + ((-16*zeta - 16)*a
39    + 75*zeta)*q^17
40
41USELESS!
42
43   [ ] Torsion bound -- should start new bound where left off from
44       old one.
45
46   [ ] Analytic invariants, e.g., ImaginaryPeriod: don't get force
47       recomputed when the precision of SOME analytic computation
48       goes up and they are re-called on.  They should!
49
50   [ ] Clean interface!! Seperate files?
51
52  *************************************************************************/
53
54HeckeBugWatch := false;
55HeckeVerbose := true;
56HeckeSuperverbose := false;
57
58/************************************************************************
59 *                                                                      *
60 *          DATA STRUCTURES                                             *
61 *                                                                      *
62 ************************************************************************/
63
64CManSymList := recformat<
65      k,      // weight
66      F,      // base field
67      R,      // F[X,Y].
68      p1list, // List of elements of P^1(Z/NZ).
69      n       // size of list of Manin symbols = #p1list * (k-1).
70> ;
71
72/************************************************************************
73 *  CQuotient:                                                          *
74 *  First we motivate the format.                                       *
75 *  To save memory and space we mod out by the S, T, and possibly I     *
76 *  relations in two phases.  First mod out by the S,I relations        *
77 *          x + xS = 0,  x - xI = 0.                                    *
78 *  by identifying equivalent pairs.  Next mod out by the T relations   *
79 *          x + xT + xT^2 = 0.                                          *
80 ************************************************************************/
81CQuotient := recformat<
82       Sgens,     // * List of the free generators after modding out
83                  //   by the 2-term S and possibly I relations.
84       Squot,     // * The i-th Manin symbol is equivalent to
85       Scoef,     //     Scoef[i]*(Squot[i]-th free generator).
86       Tgens,     // * List of the free generators after modding
87                  //   out the S-quotient by the T-relations.
88       Tquot      // * The i-th Sgen is equal to Tquot[i], which
89                  //   is a vector on Tgens.
90>;
91
92
93// The categories: (someday these will be made a part of Magma.)
94// Presently they are used to implement polymorphism in some cases.
95ModSym    := 1;
96ModSymFac := 2;
97
98function GetCategory(M)
99    return MtheCategory;
100end function;
101
102procedure SetCategory(M, category)
103    MtheCategory := category;
104end procedure;
105
106///////////////////////////////////////////////////////////////////
107// ModSym: Modular symbols object.                               //
108///////////////////////////////////////////////////////////////////
109declare attributes ModTupFld:
110         theCategory,        // * indicate that this is a mod sym object
111         k,                  // * weight
112         N,                  // * level
113         F,                  // * the base field
114         eps,                // * the Dirichlet character
115         sign,               // * sign -- see doc of ModularSymbols()
116         mlist,              // * list of Manin symbols, and other related data
117         quot,               // * gives representation of any
118                             //   Manin symbol in terms of free generators
119         disc,               // * discriminant of the Hecke algebra
120         decomp, olddecomp, newdecomp,  // * ambient decomposition
121         semidecomp,
122         gammaalp,           // * action of character via Gamma_0(N).
123         Z,                  // * integral modular symbols.
124         ModSymBasis,        // * Basis of modular symbols.
125         BoundaryMap,        // * boundary map
126         SkF,                // * cuspidal modular symbols, over F
127         SkZ,                // * cuspidal modular symbols, over integer Z.
128                             //   (only makes sense when char(F)=0).
129         SkZplus,SkZminus,   // * for *-involution.
130         old,                // * old spaces
131         sknew, sknewdual,   // * new cuspidal modular symbols
132         mknew, mknewdual,   // * new cuspidal modular symbols
133         newforms,           // * basis of eigenforms for new subspace.
134         oldforms,           // * old forms, from lower level.
135         winding,            // * list of winding submodules
136         T_images,           // * images of standard basis vectors under Hecke.
137         X,                  // * list of character groups for
138         mestre,             //   each p exactly dividing N.
139         cusplist;           // * the cusps.
140
141
142///////////////////////////////////////////////////////////////////
143// ModSymFac: Factor of a modular symbols object.                //
144// This is represented by a subspace of the dual of              //
145// a ModSym object.                                              //
146///////////////////////////////////////////////////////////////////
147declare attributes ModTupFld:
148         theCategory, // * indicates that this is a modsym factor.
149         M,         // * ambient space of modular symbols.
150         V,         // * Subspace of M
151         Z,         // * Lattice -- Integral subspace of ambient space.
152         Zdual,     // * Lattice -- Integral subspace of ambient space.
153         qexp,      // * q-expansion -- a pair, f(x), coefs.
154         eigen, eigenplus, eigenminus,    // * Eigenvector in M tensor Fbar.
155         eigenint,  // * EigenformInTermsOfIntegralBasis
156         moddeg,    // * Degree of the canonical map.
157         PeriodLattice,
158         realvolume,// * Real volume
159         imagvolume,// * Pure imaginary volume
160         cinf, cinfimag,     // * Number of real components.
161         Lvals,     // * Critical values of the L-function.
162         LRatio,      // * Ratios L(M,1)/Omega*(fudge).
163         LRatioOdd,   // * odd part of Ratios L(M,1)/Omega*(fudge).
164         cp,        // * Orders of component groups.
165         wq,        // * Atkin-Lehner
166         qintbasis, // * integral basis of q-expansions.
167         N,         // * the level
168         iso,       // * the isogeny class
169         isnew,     // * whether or not it is new
170         iscusp,    // * whether it is cuspidal or not
171         fast,      // * fast matrix of Tn on M.
172         signspace, // * A^+ and A^-.
173                    // * there is also the X from ModularSymbols...
174         xdata,     // * local phix and mx
175         ZxZalp, VolLEven, VolLOdd, // * For LRatio.
176         PeriodGens, PGfast,
177         RatPeriodMap, RatPeriodLat,
178         RatPeriodConj,RatPeriodMapSign,
179         PeriodMap, PeriodMapPrecision;
180
181
182
183/************************************************************************
184 *                                                                      *
185 * Manin Symbols (basic functions)                                      *
186 *                                                                      *
187 ************************************************************************/
188
189/////////////////////////////////////////////////////////////////////////
190// NormalizePair  (written by Allan Steel)                             //
191// INPUT: [u,v] in Z/N x Z/N.                                          //
192// OUTPUT:  1) the *index* of a fixed choice [a,b] of representative   //
193//             in the P^1(Z/NZ) equivalence class of [u,v].            //
194//             If gcd(u,v,N)>1 then the index returned is 0.           //
195//          2) a scalar s such that                                    //
196//                    [a*s,b*s] = [u,v].                               //
197//             The character relation is thus                          //
198//                  eps(s) [a,b] == [u,v].                             //
199/////////////////////////////////////////////////////////////////////////
200function NormalizePair(x)
201   Z := IntegerRing();
202   u := x[1];
203   v := x[2];
204   R := Parent(u);
205   N := Modulus(R);
206   if u eq 0 then
207      return [R | 0, 1], v;
208   else
209      u := Z ! u;
210      g, s := XGCD(u, N);
211      if g ne 1 then
212         d := N div g;
213	 while GCD(s, N) ne 1 do
214	    s := (s + d) mod N;
215	 end while;
216      end if;
217      // Multiply (u, v) by s; new u = g
218      u := R!g;
219      v := v*s;
220      minv := Z!v;
221      mint := 1;
222      if g ne 1 then
223	 Ng := N div g;
224	 vNg := v*Ng;
225	 t := R!1;
226	 for k := 2 to g do
227	    v +:= vNg;
228	    t +:= Ng;
229	    if Z!v lt minv and IsUnit(t) then
230	       minv := Z!v;
231	       mint := t;
232	    end if;
233	 end for;
234      end if;
235      if GCD([Z | u, minv, N]) ne 1 then
236	 printf "Bad x=%o, N=%o, u=%o, minv=%o, s=%o, mint=%o\n",
237	     x, N, u, minv, s, mint;
238	 error "";
239      end if;
240//      return [R|u, minv], R!InverseMod(s*mint,N);
241      return [R|u, minv], 1/(R!(s*mint));
242   end if;
243end function;
244
245
246//////////////////////////////////////////////////////////////////////////
247//  P1Reduce:                                                           //
248//  INPUT: [u,v] in Z/N x Z/N.                                          //
249//  OUTPUT:  1) the *index* of a fixed choice [a,b] of representative   //
250//              in the P^1(Z/NZ) equivalence class of [u,v].            //
251//              If gcd(u,v,N)>1 then the index returned is 0.           //
252//           2) a scalar s such that                                    //
253//                     [a*s,b*s] = [u,v].                               //
254//              so the character relation is                            //
255//                   eps(s) [a,b] == [u,v].                             //
256//////////////////////////////////////////////////////////////////////////
257function P1Reduce(x, list)
258   N := Characteristic(Parent(x[1]));
259   if N eq 1 then
260      return 1, 1;
261   end if;
262   if Gcd([x[1], x[2], N]) ne 1 then
263      return 0, 1;
264   end if;
265   n, s := NormalizePair(x);
266   return Index(list, n), s;
267end function;
268
269
270//////////////////////////////////////////////////////////////////////////
271//  EnumerateP1:                                                        //
272//  Construct a list of distinct normalized representatives for         //
273//  the set of equivalence classes of P^1(Z/NZ).                        //
274//  The output of this function is used by P1Reduce.                    //
275//////////////////////////////////////////////////////////////////////////
276function EnumerateP1(N)
277   R := (N gt 1) select IntegerRing(N) else quo<Integers()|1>;
278   return {@ NormalizePair([R|c,1]) : c in [0..N-1] @} join
279          {@ NormalizePair([R|1,d]) : d in [0..N-1] | Gcd(d,N) gt 1 @} join
280	  {@ NormalizePair([R|c,d]) : c in Divisors(N), d in [0..N-1] |
281                       c ne 1 and c ne N and Gcd(c,d) eq 1
282                       and Gcd(d,N) gt 1 @};
283end function;
284
285
286//////////////////////////////////////////////////////////////////////////
287// ManinSymbolList:                                                     //
288// Construct a list of distinct Manin symbols.  These are               //
289// elements of the Cartesion product:                                   //
290//     {0,...,k-2} x P^1(Z/NZ).                                         //
291// In fact, we only store a list of the elements of P^1(Z/NZ),          //
292// as the full Cartesion product can be stored using                    //
293// the following scheme.                                                //
294//                                                                      //
295// There are Manin symbols 1,...,#({0,..,k-2}xP^1) indexed by i.        //
296// Thus i is an integer between 1 and the number of generating Manin    //
297// symbols.  Suppose P^1(N) = {x_1, ..., x_n}.  The encoding is         //
298// as follows:                                                          //
299//   1=<X^0Y^(k-2),x_1>,  2=<0, x_2>, ..., n=<0,x_n>,                   //
300// n+1=<X^1Y^(k-3),x_1>,n+2=<1, x_2>, ...,2n=<1,x_n>,                   //
301// ...                                                                  //
302//////////////////////////////////////////////////////////////////////////
303function ManinSymbolList(k,N,F)
304   p1list := EnumerateP1(N);
305   n      := (k-1)*#p1list;
306   return rec<CManSymList |
307      k      := k,
308      F      := F,
309      R      := PolynomialRing(F,2),
310      p1list := p1list,
311      n      := n
312   >;
313end function;
314
315
316///////////////////////////////////////////////////////////////////////////
317// ManinSymbolApply                                                      //
318// Apply an element g=[a,b, c,d] of SL(2,Z) to the i-th                  //
319// standard Manin symbol.  The definition of the action is               //
320//       g.(X^i*Y^j[u,v]) :=                                             //
321//                     (aX+bY)^i*(cX+dY)^j [[u,v]*g].                    //
322// The argument "manin" should be as returned by                         //
323// the function ManinSymbolsList above.                                  //
324// The character eps should be an array which                            //
325// defines a Dirichlet character Z/NZ ---> K.                            //
326// Thus eps[a] in K and makes sense for a in (Z/NZ)^*.                   //
327// If eps is trivial, then the character is assumed trivial.             //
328///////////////////////////////////////////////////////////////////////////
329
330
331///////////////////////////////////////////////////////////////////////////
332// Unwind -- Given an int i, this function gives back the                //
333// ith generating Manin symbol.                                          //
334///////////////////////////////////////////////////////////////////////////
335function UnwindManinSymbol(i, mlist)
336   // P^1(N) part.
337   p1list := mlistp1list;
338   n := #p1list;
339   ind := (i mod n);
340   if ind eq 0 then
341      ind := n;
342   end if;
343   uv := p1list[ind];
344   w := Integers()!((i - ind)/n);
345   return uv, w, ind;
346end function;
347
348// Wind -- Given w in the range 0<=w<=k-2, and an index ind
349//   into the P^1-list, this function gives back
350//   the number of the generating Manin symbol.
351function WindManinSymbol(w, ind, mlist)
352   return w*#(mlistp1list) + ind;
353end function;
354
355
356intrinsic ManinSymbolApply(g::SeqEnum, i::RngIntElt,
357               mlist::Rec, eps::Rec) -> SeqEnum
358{Apply g to the ith Manin symbol.}
359   k := mlistk;
360   uv, w, ind := UnwindManinSymbol(i,mlist);
361   p1list := mlistp1list;
362   uvg    := Parent(uv) ! [g[1]*uv[1]+g[3]*uv[2], g[2]*uv[1]+g[4]*uv[2]];
363   act_uv, scalar := P1Reduce(uvg, p1list);
364
365   if act_uv eq 0 then
366      return [<0,1>];
367   end if;
368
369   if k eq 2 then
370      if not IsTrivial(eps) then
371         return [<Evaluate(eps,scalar),act_uv>];
372/*
373         if Evaluate(eps,scalar) ne 0 then
374            return [<Evaluate(eps,scalar)^(-1),act_uv>];
375         else
376            return [<0,act_uv>];
377         end if;
378*/
379      else
380         return [<1,act_uv>];
381      end if;
382   end if;
383
384   // Polynomial part.
385   R<x> :=PolynomialRing(mlistF);    // univariate is fine for computation
386   hP := (g[1]*x+g[2])^w*(g[3]*x+g[4])^(k-2-w);
387
388   if not IsTrivial(eps)then
389      hP *:= Evaluate(eps,scalar);
390   end if;
391   pol := ElementToSequence(hP);
392
393   // Put it together
394   n   := #p1list;
395   ans := [<pol[w+1],  w*n + act_uv> : w in [0..#pol-1]];
396   return [x : x in ans | x[1] ne 0];
397end intrinsic;
398
399
400/////////////////////////////////////////////////////////////////
401//  ModularSymbolApply                                         //
402//  Apply an element g=[a,b, c,d] of SL(2,Z) to the given      //
403//  modular symbol.  The definition of the action is           //
404//       g.(X^i*Y^j{u,v}) :=                                   //
405//                     (dX-bY)^i*(-cX+aY)^j {g(u),g(v)}.       //
406//  A modular symbol is represented as a sequence of pairs     //
407//  <P,x>, where P is a polynomial in X and Y (an element      //
408//  of the ring MmlistR, and x=[[a,b],[c,d]] is a pair       //
409//  of elements of P^1(Q),                                     //
410//     where [a,b] <--> a/b and [c,d] <--> c/d.                //
411//  After computing the action, no further reduction is done.  //
412/////////////////////////////////////////////////////////////////
413
414intrinsic ModularSymbolApply( M::ModTupFld, g::SeqEnum, s::Tup
415             ) -> SeqEnum
416{Apply an element g=[a,b, c,d] of SL(2,Z) to the given
417 modular symbol.  The definition of the action is
418       g.(X^i*Y^j\{u,v\}) :=
419                     (dX-bY)^i*(-cX+aY)^j \{g(u),g(v)\}.
420  A modular symbol is represented as a pair
421  <P,x>, where P is a polynomial in X and Y (an element
422  of the ring MmlistR, and x=[[a,b],[c,d]] is a pair
423  of elements of P^1(Q),
424     where [a,b] <--> a/b and [c,d] <--> c/d.}
425
426   R  := MmlistR;
427   h  := hom <R -> R  |  g[4]*R.1-g[2]*R.2, -g[3]*R.1+g[1]*R.2>;
428   hP := h(s[1]);
429   a,b:= Explode(s[2][1]);
430   c,d:= Explode(s[2][2]);
431   gx := [[g[1]*a+g[2]*b, g[3]*a+g[4]*b],
432          [g[1]*c+g[2]*d, g[3]*c+g[4]*d]];
433   return <hP, gx>;
434end intrinsic;
435
436intrinsic ModularSymbolApply(M::ModTupFld, g::SeqEnum, s::SeqEnum) -> SeqEnum
437{Apply an element g=[a,b, c,d] of SL(2,Z) to the given
438 modular symbol.  The definition of the action is
439       g.(X^i*Y^j\{u,v\}) :=
440                     (dX-bY)^i*(-cX+aY)^j \{g(u),g(v)\}.
441  A modular symbol is represented as a sequence of pairs
442  <P,x>, where P is a polynomial in X and Y (an element
443  of the ring MmlistR, and x=[[a,b],[c,d]] is a pair
444  of elements of P^1(Q),
445     where [a,b] <--> a/b and [c,d] <--> c/d.}
446   return [ModularSymbolApply(M, g,term): term in s];
447end intrinsic;
448
449
450/////////////////////////////////////////////////////////////////
451//  Sparse2Quotient:                                           //
452//   Performs Sparse Gauss elimination on a matrix all of      //
453//   whose columns have at most 2 nonzero entries.  I just     //
454//   use the fairly obvious algorithm.   It runs fast          //
455//   enough.  Example:                                         //
456//   rels := [[1,2], [1,3], [2,3], [4,5]];                     //
457//   relc := [[3,-1],[1,1], [1,1], [1,-1]];                    //
458//   n    := 5;                // x1,...,x5                    //
459//   corresponds to 3*x1-x2=0, x1+x3=0, x2+x3=0, x4-x5=0.      //
460/////////////////////////////////////////////////////////////////
461
462function Sparse2Quotient (rels, relc, n, F)
463   free := [1..n];
464   coef := [F|1 : i in [1..n]];
465   related_to_me := [[] : i in [1..n]];
466   for i in [1..#rels] do
467      t := rels[i];
468      c1 := relc[i][1]*coef[t[1]];
469      c2 := relc[i][2]*coef[t[2]];
470      // Mod out by the relation
471      //    c1*x_free[t[1]] + c2*x_free[t[2]] = 0.
472      die := 0;
473      if c1 eq 0 and c2 eq 0 then
474         // do nothing.
475      elif c1 eq 0 and c2 ne 0 then  // free[t[2]] --> 0
476         die := free[t[2]];
477      elif c2 eq 0 and c1 ne 0 then
478         die := free[t[1]];
479      elif free[t[1]] eq free[t[2]] then
480         if c1+c2 ne 0 then
481            // all xi equal to free[t[1]] must now equal to zero.
482            die := free[t[1]];
483         end if;
484      else   // x1 = -c2/c1 * x2.
485         x := free[t[1]];
486         free[x] := free[t[2]];
487         coef[x] := -c2/c1;
488         for i in related_to_me[x] do
489	     free[i] := free[x];
490             coef[i] *:= coef[x];
491	     Append (~related_to_me[free[t[2]]], i);
492         end for;
493	 Append (~related_to_me[free[t[2]]], x);
494      end if;
495
496      if die gt 0 then
497         for i in related_to_me[die] do
498            free[i] := 1;
499            coef[i] := 0;
500         end for;
501	 free[die] := 1 ;
502	 coef[die] := 0;
503      end if;
504   end for;
505
506   // Enumerate the subscripts of free generators that survived.
507   // x_{i_1}  <----> y_1
508   // x_{i_2}  <----> y_2, etc.
509
510   for i in [1..#free] do
511      if coef[i] eq 0 then
512         free[i] := -1;
513      end if;
514   end for;
515   ykey := {@ x : x in Set(free) | x ne -1 @};
516   for i in [1..#free] do
517      if free[i] eq -1 then
518         free[i] := 1;
519      else
520         free[i] := Index(ykey,free[i]);
521      end if;
522   end for;
523   return ykey, free, coef;
524end function;
525
526
527
528/*********************************************************
529  Quotient:
530  The INPUT is a list of (sparse) relations
531    [[<c_1, i_1>, <c_2, i_2>, ... <c_r,i_r>], ... ]
532  The first listed spare relation is
533     c_1*e_{i_1} + c_2*e_{i_2} + ... c_r*e_{i_r} = 0.
534  The integer n denotes the total number of basis
535  elements e_i.
536  The field K is the field over which the c_i are defined.
537  The OUTPUT is (1) an indexed set of g free generators, and
538  (2) an expression for each e_i in terms of the free
539  generators.  These expressions are elements of the
540  g-dimensional vector space over K.
541       generators,   quotient
542 *********************************************************/
543function Pivots(B)   // find pivots of reduced basis.
544   pivots := [];
545   for b in B do
546      j := 1;
547      while b[j] eq 0 do
548         j +:= 1;
549      end while;
550      Append(~pivots, j);
551   end for;
552   return pivots;
553end function;
554
555forward Quotient2;
556function Quotient(rels, n, K)
557   ringchoice := Characteristic(K) eq 0;
558   if Type(K) eq FldPr then
559      ringchoice := false;
560   end if;
561   return Quotient2(rels,n,K : ring:=ringchoice);
562
563end function;
564function Quotient2(rels, n, K : ring := true)
565   if HeckeSuperverbose then
566      "# relations = ", #rels;
567   end if;
568   if ring then    // POLYNOMIAL ALGORITHM
569      if HeckeSuperverbose then
570         "ring reduction";
571      end if;
572      R    := PolynomialRing(K, n);
573      vars := {@ R.i : i in [1..n] @};
574      Rels := [ &+[t[1]*R.t[2] : t in r] : r in rels];
575      Rels := Reduce(Rels);
576      W    := quo<R|Rels>;
577      gens := {@ i : i in [1..n] | not exists
578                 { g : g in Rels | LeadingMonomial(g) eq R.i } @};
579      Z    := VectorSpace(K,#gens);
580      quot := [ &+[Z|
583                   : t in Terms(R!W.i)] : i in [1..n] ];
584   else            // VECTOR SPACE ALGORITHM
585      if HeckeSuperverbose then
586         "vector space reduction";
587      end if;
588      V := VectorSpace(K, n);
589      Rels := [ &+[t[1]*V.t[2] : t in r] : r in rels];
590      S := RowSpace(RMatrixSpace(K, #Rels, n) ! Rels);
591      m := n - Dimension(S);
592      B := Basis(S);
593      W := VectorSpace(K, m);
594      quot := [ W!0 : i in [1..n]];
595
596      pivots := Pivots(B);
597
598      gens := {@ i : i in [1..n] | i notin pivots @};
599      for i := 1 to #pivots do
600         quot[pivots[i]] := -W![B[i][j] : j in gens];
601      end for;
602      for i :=1 to #gens do
603         quot[gens[i]] := W.i;
604      end for;
605   end if;
606   return gens, quot;
607end function;
608
609
610////////////////////////////////////////////////////////////////////
611//  CONSTRUCT QUOTIENT BY 2 AND 3 TERM RELS                       //
612////////////////////////////////////////////////////////////////////
613
614function ManSym2termQuotient (mlist, eps, sign)
615     n := mlistn;
616     K := mlistF;
617     S := [0,-1,  1,0];
618     I := [-1,0,  0,1];
619     xS := [ManinSymbolApply(S,i,mlist,eps) : i in [1..n]];
620     S_rels := [ [i,   (xS[i][1])[2]] : i in [1..n]| #xS[i] gt 0];
621     S_relc := [ [K!1, (xS[i][1])[1]] : i in [1..n]| #xS[i] gt 0];
622     if sign ne 0 then
623        xI := [ManinSymbolApply(I,i,mlist,eps) : i in [1..n]];
624        I_rels := [ [i,    (xI[i][1])[2]] : i in [1..n]];
625        I_relc := [ [K!1, -sign*(xI[i][1])[1]] : i in [1..n]];
626     else
627        I_rels := [];
628        I_relc := [];
629     end if;
630     rels  := S_rels cat I_rels;
631     relc  := S_relc cat I_relc;
632     return Sparse2Quotient(rels,relc,n,K);
633end function;
634
635function ManSym3termQuotient (mlist, eps, Skey, Squot, Scoef)
636   // Construct the subspace of 3-term relations.
637   n := mlistn;
638   F := mlistF;
639   k := mlistk;
640   T := [0,-1,  1,-1];
641   TT:= [-1,1, -1,0];
642
643   if k eq 2 then
644      mask := [false : i in [1..n]];
645      rels := [];
646      for j in [1..n] do
648            t  := ManinSymbolApply(T,j,mlist,eps)[1];
649            tt := ManinSymbolApply(TT,j,mlist,eps)[1];
652            Append(~rels,  [<Scoef[j],Squot[j]>,
653               <t[1]*Scoef[t[2]],Squot[t[2]]>,
654               <tt[1]*Scoef[tt[2]],Squot[tt[2]]>]);
655         end if;
656      end for;
657   else
658      rels := [&cat[
659               [<Scoef[i],Squot[i]>],
660               [<x[1]*Scoef[x[2]],Squot[x[2]]>
661                 : x in ManinSymbolApply(T,i,mlist,eps)],
662               [<x[1]*Scoef[x[2]],Squot[x[2]]>
663                 : x in ManinSymbolApply(TT,i,mlist,eps)]
664              ]
665             : i in [1..n]];
666   end if;
667
668   return Quotient(rels, #Skey, F);
669end function;
670
671
672/*********************************************************
673 *                                                       *
674 * MODULARSYMBOLS:                                       *
675 * Construct the space of modular symbols of weight k,   *
676 * level N and character eps.                            *
677 *                                                       *
678 *********************************************************/
679
680intrinsic ModularSymbols(N::RngIntElt) -> ModTupFld
681{Returns the space of modular symbols
682 of weight 2 for Gamma_0(N), over the rational numbers.}
683   return ModularSymbols(DirichletCharacter(N),2,0);
684end intrinsic;
685
686intrinsic ModularSymbols(N::RngIntElt, k::RngIntElt) -> ModTupFld
687{Returns the space of modular symbols
688 of weight k for Gamma_0(N), over the rational numbers.}
689   return ModularSymbols(DirichletCharacter(N),k,0);
690end intrinsic;
691
692intrinsic ModularSymbols(N::RngIntElt, k::RngIntElt,
693                               sign::RngIntElt) -> ModTupFld
694{Returns the space of modular symbols
695 of weight k for Gamma_0(N), over the rational numbers.
696 If sign=+1 then returns the +1 quotient, if sign=-1
697 returns the -1 quotient, and if sign=0 returns the full space.
698 If sign is a positive prime, then the full space of modular symbols
699 mod sign isreturned.}
700   if sign gt 1 and IsPrime(sign) then
701      return ModularSymbols(N,k,0,sign);
702   end if;
703   return ModularSymbols(DirichletCharacter(N),k,sign);
704end intrinsic;
705
706intrinsic ModularSymbols(N::RngIntElt, k::RngIntElt,
707                         sign::RngIntElt,
708                         p::RngIntElt) -> ModTupFld
709{Returns the the space of modular symbols
710 of weight k for Gamma_0(N), over the finite field Fp.}
711   assert IsPrime(p);
712   return ModularSymbols(DirichletCharacter(N,p),k,sign);
713end intrinsic;
714
715intrinsic ModularSymbols(eps::Rec) -> ModTupFld
716{Returns the space of modular symbols
717 of weight k for Gamma_1(N) with character eps,
718 over the rational numbers.}
719   return ModularSymbols(eps,2,0);
720end intrinsic;
721
722intrinsic ModularSymbols(eps::Rec, k::RngIntElt) -> ModTupFld
723{Returns the space of modular symbols
724 of weight k for Gamma_1(N) with character eps,
725 over the rational numbers.}
726   return ModularSymbols(eps,k,0);
727end intrinsic;
728
729forward TrivialSpace;
730intrinsic ModularSymbols(eps::Rec, k::RngIntElt,
731                           sign::RngIntElt) -> ModTupFld
732{Returns the space of modular symbols of character eps and weight k.
733 The level and base field are specified by the DirichletCharacter
734 eps.  The third argument "sign" allows for working in certain
735 quotients.  The possible values are -1, 0, or +1, which correspond
736 to the -1 quotient, full space, and +1 quotient.}
737   if not (-1 le sign and sign le 1) then
738      error "sign must be -1, 0, or 1";
739   end if;
740   if k lt 2 then
741      error "Modular symbols are only defined for weight k >= 2.";
742   end if;
743
744   N := CharacterLevel(eps);
745   F := FieldOfValues(eps);
746   if HeckeVerbose then
747      tt := Cputime();
748      printf "Creating M_%o(Gamma_1(%o),eps;F_%o)",
749               k, N, Characteristic(F) eq 0 select 0 else #F;
750      if sign eq 1 then
751         printf "^+";
752      elif sign eq -1 then
753         printf "^-";
754      end if;
755      if not HeckeSuperverbose then
756         printf "\n";
757      end if;
758   end if;
759
760   if Evaluate(eps,-1) ne (-1)^k then
761       return TrivialSpace(k,eps,sign);
762   end if;
763
764   if HeckeSuperverbose then
765      t := Cputime(); "Step I.   Manin symbols list.";
766   end if;
767   mlist := ManinSymbolList(k,N,F);
768   if HeckeSuperverbose then
769      Cputime(t), " seconds.";
770   end if;
771
772   if HeckeSuperverbose then
773      t := Cputime(); "Step II.  2-term relations.";
774   end if;
775   Sgens, Squot, Scoef := ManSym2termQuotient(mlist, eps, sign);
776   if HeckeSuperverbose then
777      Cputime(t), " seconds.";
778   end if;
779
780   if #Sgens lt 1 then
781       return TrivialSpace(k,eps,sign);
782   end if;
783
784   if HeckeSuperverbose then
785      t := Cputime(); "Step III. 3-term relations.";
786   end if;
787   Tgens, Tquot := ManSym3termQuotient(mlist, eps, Sgens, Squot, Scoef);
788   if HeckeSuperverbose then
789      Cputime(t), " seconds.";
790   end if;
791   if #Tgens lt 1 then
792       return TrivialSpace(k,eps,sign);
793   end if;
794
795   V := VectorSpace(F,#Tgens);
796
797   // The reason we use "WithBasis" is that for some reason Magma
798   // distinguishes different elements of the category ModTupFld
799   // created via this constructor (it still says they are equal but
800   // their attributes are kept distinct).  If we did not use WithBasis,
801   // then any two space of modular symbols of the same dimension
802   // would be equal.   Someday this will be replaced with
803   // MY OWN CATEGORY.  It is best to use this hack for now,
804   // so in the future, when I get a category, minimal changes
805   // will be necessary.
806   M := VectorSpaceWithBasis(Basis(V));
807   SetCategory(M,ModSym);
808   Mk    := k;
809   MN    := N;
810   Meps  := eps;
811   Msign := sign;
812   MF    := F;
813   Mmlist:= mlist;
814   Mquot := rec<CQuotient |
815                 Sgens := Sgens,
816                 Squot := Squot,
817                 Scoef := Scoef,
818                 Tgens := Tgens,
819                 Tquot := Tquot>;
820   // MquotTquot lies in V and really M is almost a V, so there
821   // is a possibility of a "circular reference bug" in Magma
822   // leading to a memory leak!
823   MquotTquot := [V!v : v in MquotTquot];  // move them into M.
824   if HeckeVerbose then
825      ", ",Cputime(tt), "seconds.";
826   end if;
827   return M;
828end intrinsic;
829
830function TrivialSpace(k,eps,sign)
831   F := FieldOfValues(eps);
832   V := VectorSpace(F,0);
833   N := CharacterLevel(eps);
834   M := VectorSpaceWithBasis(Basis(V));
835   Mk := k;
836   MN := N;
837   Meps := eps;
838   Msign := sign;
839   MF := F;
840   Mquot := rec<CQuotient | Sgens:=[], Squot:=[],
841               Scoef:=[], Tgens := [], Tquot:=[]>;
842   SetCategory(M,ModSym);
843   return M;
844end function;
845
846intrinsic Weight(M::ModTupFld) -> RngIntElt
847{Return the weight of the space M of modular symbols.}
848   if not assigned Mk then return Weight(ModSymParent(M)); end if;
849   return Mk;
850end intrinsic;
851
852intrinsic Level(M::ModTupFld) -> RngIntElt
853{Return the level of the space M of modular symbols.}
854   return MN;
855end intrinsic;
856
857intrinsic Conductor(M::ModTupFld) -> RngIntElt
858{Return the level of the space M of modular symbols.}
859   return MN;
860end intrinsic;
861
862intrinsic BaseField(M::ModTupFld) -> .
863{Return the base field of the space M of modular symbols.}
864   if not assigned MF then
865      return BaseField(ModSymParent(M));
866   end if;
867   return MF;
868end intrinsic;
869
870intrinsic Character(M::ModTupFld) -> SeqEnum
871{Return the Dirichlet character of the space M of modular symbols.}
872   case GetCategory(M):
873      when ModSym:
874         return Meps;
875      when ModSymFac:
876         return MMeps;
877      else:
878         error "Character -- invalid input type.";
879   end case;
880end intrinsic;
881
882intrinsic IsPlusQuotient(M::ModTupFld) -> SeqEnum
883{Return true iff working in the plus one quotient of the
884  space M of modular symbols.}
885   if not assigned Msign then return MMsign eq 1; end if;
886   return Msign eq 1;
887end intrinsic;
888
889intrinsic IsMinusQuotient(M::ModTupFld) -> SeqEnum
890{Return true iff working in the plus one quotient of
891 the space M of modular symbols.}
892   if not assigned Msign then return MMsign eq -1; end if;
893   return Msign eq -1;
894end intrinsic;
895
896intrinsic Sign(M::ModTupFld) -> RngIntElt
897{Returns the sign of M.}
898   if not assigned Msign then return MMsign; end if;
899   return Msign;
900end intrinsic;
901
902intrinsic IsogenyClass(A::ModTupFld) -> RngIntElt
903{Returns the isogeny class of A.}
904   if not assigned Aiso then
905      return 0;
906   end if;
907   return Aiso;
908end intrinsic;
909
910
911
912
913/*********************************************************
914 *                                                       *
915 *        COMPUTE HECKE OPERATORS                        *
916 *                                                       *
917 *********************************************************/
918intrinsic HeilbronnMerel(n::RngIntElt) -> SeqEnum
919{Return Heilbronn matrices of determinant n, as given by Merel.
920   Here n can be composite.}
921   H := [];
922   i := 0;
923   for a in [1..n] do
924      for d in [1..n] do;
926        bc := a*d - n;
927        if bc lt 0 then
928           continue;
929        end if;
930        if bc eq 0 then
931           for b in [0..a-1] do
932              Append(~H,[a,b,0,d]);
933           end for;
934        end if;
935        if bc ne 0 then
936           for c in Divisors(Abs(bc)) do
937              if c lt d and Floor(bc/c) lt a then
938                 Append(~H,[a,Floor(bc/c),c,d]);
939              end if;
940           end for;
941        elif 0 lt a then
942           for c in [1..d-1] do
943              Append(~H,[a,0,c,d]);
944           end for;
945        end if;
946     end for;
947   end for;
948   return H;
949end intrinsic;
950
951
952intrinsic HeilbronnCremona(p::RngIntElt) -> SeqEnum
953{Return the Heilbronn matrices of determinant p, as defined by Cremona.}
954   assert IsPrime(p);
955   if p eq 2 then
956      return [[1,0, 0,2], [2,0, 0,1], [2,1, 0,1], [1,0, 1,2]];
957   end if;
958   ans := [[1,0, 0,p]];
959   for r in [Ceiling(-p/2)..Floor(p/2)] do
960      x1:=p; x2:=-r; y1:=0; y2:=1; a:=-p; b:=r; c:=0; x3:=0; y3:=0; q:=0;
961      Append(~ans, [x1,x2, y1,y2]);
962      while b ne 0 do
963         q := Round(a/b);
964         c := a-b*q;
965         a := -b;
966         b := c;
967         x3 := q*x2-x1;
968         x1 := x2; x2 := x3; y3 := q*y2-y1; y1 := y2; y2 := y3;
969         Append(~ans, [x1,x2, y1, y2]);
970      end while;
971   end for;
972   return ans;
973end intrinsic;
974
975intrinsic Heilbronn(n::RngIntElt) -> SeqEnum
976{If n is prime, return Cremona's Heilbronn matrices; otherwise
977return Merel's.}
978   requirege n,1;
979   if IsPrime(n) then
980      return HeilbronnCremona(n);
981   else
982      return HeilbronnMerel(n);
983   end if;
984end intrinsic;
985
986// m is a sequence <alp, i> representing sum alp*e_i,
987// where e_i runs through the initial list of
988// generating Manin symbols.
989intrinsic InitialManinSymbolGenListToSquot(m::SeqEnum, M::ModTupFld,
990          Tmat::ModMatFldElt)
991       -> ModTupFldElt
992   {}
993   Tquot := MquotTquot;
994   V:=Parent(Tquot[1]);
995   Sgens := MquotSgens;
996   Scoef := MquotScoef;
997   Squot := MquotSquot;
998   // Tmat is a map from W to V, where W is
999   W:=VectorSpace(Field(V),#Sgens);
1000   v := W!0;
1001   for t in m do
1002      v[Squot[t[2]]] +:= t[1]*Scoef[t[2]];
1003   end for;
1004   return v;
1005end intrinsic;
1006
1007intrinsic InitialManinSymbolGenListToM(m::SeqEnum, M::ModTupFld)
1008       -> ModTupFldElt
1009   {}
1010   Scoef := MquotScoef;
1011   Tquot := MquotTquot;
1012   V:=Parent(Tquot[1]);
1013   Squot := MquotSquot;
1014   return &+[V|t[1]*Scoef[t[2]]
1015                 *Tquot[Squot[t[2]]] : t in m];
1016end intrinsic;
1017
1018///////////////////////////////////////////////////////////////
1019//  COMPUTATION of HECKE OPERATORS                           //
1020///////////////////////////////////////////////////////////////
1021
1022intrinsic Tn(M::ModTupFld, n::RngIntElt) -> AlgMatElt
1023{}
1024   return HeckeOperator(M,n);
1025end intrinsic;
1026
1027intrinsic HeckeOperator(M::ModTupFld, n::RngIntElt) -> AlgMatElt
1028{}
1029   requirege n,1;
1030   return HeckeOperator(M,Heilbronn(n));
1031end intrinsic;
1032
1033intrinsic HeckeOperator(M::ModTupFld, Heil::SeqEnum) -> AlgMatElt
1034{Compute the n-th Hecke operator Tn on the space M of modular symbols.}
1035   if Dimension(M) eq 0 then
1036      return Hom(M,M)!0;
1037   end if;
1038
1039   F     := MF;
1040   mlist := Mmlist;
1041   eps   := Meps;
1042
1043   quot  := Mquot;
1044   Sgens := quotSgens;
1045   Squot := quotSquot;
1046   Tgens := quotTgens;
1047   Tquot := quotTquot;
1048
1049
1050   Mn := MatrixAlgebra(F,Dimension(M));
1051   V := Parent(Tquot[1]);
1052   T := [];
1053
1054   W := VectorSpace(F, #Sgens);
1055   Tmat := RMatrixSpace(F, Dimension(W), Dimension(V)) !
1056                [Tquot[i] : i in [1..#Sgens]];
1057
1058   for j in [1..#Tgens] do
1059      i := Sgens[Tgens[j]];
1060      v := InitialManinSymbolGenListToSquot(
1061                   &cat[ ManinSymbolApply(h,i, mlist,eps)
1062                       : h in Heil],
1063                M,Tmat);
1064      Append(~T,v);
1065   end for;
1066   S := RMatrixSpace(F, #T, #Sgens) ! T;
1067   return Mn!(S*Tmat);
1068end intrinsic;
1069
1070intrinsic TnSparse(M::ModTupFld,  n::RngIntElt, sparsevec::SeqEnum)
1071                      -> AlgMatElt
1072{}
1073   requirege n,1;
1074   if HeckeSuperverbose then
1075       "TnSparse ", n;
1076   end if;
1077   return TnSparse(M, Heilbronn(n), sparsevec);
1078end intrinsic;
1079
1080
1081intrinsic TnSparse(M::ModTupFld, Heil::SeqEnum, sparsevec::SeqEnum)
1082                      -> AlgMatElt
1083{Compute the action of the Hecke operator defined by the
1084Heilbronn matrices Heil on the sparse vector v.}
1085
1086   if Dimension(M) eq 0 then
1087      return M!0;
1088   end if;
1089   F     := MF;
1090   k     := Mk;
1091   mlist := Mmlist;
1092   eps   := Meps;
1093   quot  := Mquot;
1094   Sgens := quotSgens;
1095   Tgens := quotTgens;
1096   Tquot := quotTquot;
1097
1098   Mn := MatrixAlgebra(F,Dimension(M));
1099   v := M!0;
1100   for m in sparsevec do
1101     i := Sgens[Tgens[m[2]]];
1102     if k eq 2 then     // faster, since always just one entry.
1103        x := [ ManinSymbolApply(h,i, mlist,eps)[1] : h in Heil];
1104     else
1105        x := &cat[ ManinSymbolApply(h,i, mlist,eps) : h in Heil];
1106     end if;
1107     v +:= m[1]*InitialManinSymbolGenListToM(x,M);
1108   end for;
1109   return v;
1110end intrinsic;
1111
1112
1113
1114/*******************************************************
1115   CyclicHeckeImage:
1116   The goal of this function is to compute
1117   the cyclic Hecke module generated by a vector v.
1118   The result is returned as a subspace with Z-basis.
1119
1120   Let Tn be the matrix representing the right
1121   action of the nth Hecke operator on P.
1122   We must compute gens for an F-module who (Z-span) equals
1123         span {Tn*v : 1<=n<=b}
1124   where b is a bound on the number of Tn necessary
1125   to generate the Hecke algebra as a Z-module.
1126
1127   We do this by:
1128    1) Computing Mp for p prime using Hecke action on
1129       easy elements of M_k(Q).
1130    2) Iteratively computing prod Mp^i*c where n=prod
1131       p^i varies from 1 to b.
1132*********************************************************/
1133
1134intrinsic SparseRepresentation(v::ModTupFldElt) -> SeqEnum
1135{Return sparse representation of vector v.}
1136   sp := [];
1137   v  := Eltseq(v);
1138   for i in [1..#v] do
1139      if v[i] ne 0 then
1140         Append(~sp, <v[i],i>);
1141      end if;
1142   end for;
1143   return sp;
1144end intrinsic;
1145
1146
1147intrinsic HeckeBound(M::ModTupFld) -> RngIntElt
1148{Returns a bound b such that T1, ..., Tb generate the Hecke
1149 algebra as a Z-module.}
1150   b := HeckeBound(Mk, MN);
1151   if not IsTrivial(Meps) then
1152      "  Warning: William Stein has not proved that his bound on the number";
1153      "  of Hecke operators needed to generate the Hecke algebra is correct";
1154      "  when Character(M) is not trivial.";
1155      b *:= 3;
1156   end if;
1157   return b;
1158end intrinsic;
1159
1160forward idxG0;
1161intrinsic HeckeBound(k::RngIntElt, N::RngIntElt) -> RngIntElt
1162{Returns a bound b such that T1, ..., Tb generate the Hecke
1163 algebra as a Z-module.}
1164   return Ceiling(k * idxG0(N) / 12);
1165end intrinsic;
1166
1167intrinsic VectorSpaceZBasis(B::SeqEnum) -> ModTupFld
1168{Return vector space with basis the lattice defined by the elements of B.}
1169   V   := VectorSpace(Parent(B[1][1]),#Eltseq(B[1]));
1170   Mat := RMatrixSpace(Rationals(),#B, Dimension(V));
1171   X   := Mat!B;
1172   if X eq Mat!0 then
1173      return VectorSpaceWithBasis([V|]);
1174   end if;
1175   Latbase := Basis(Lattice(X));
1176   return VectorSpaceWithBasis([V!v : v in Latbase]);
1177end intrinsic;
1178
1179intrinsic HeckeSpan(M::ModTupFld, v::ModTupFldElt) -> ModTupFld
1180{This function computes
1181the Hecke module generated by a vector v.
1182The result is returned as a subspace with Z-basis.}
1183   s := SparseRepresentation(v);
1184   b := HeckeBound(M);
1185   B := [v] cat [TnSparse(M, n, s) : n in [2..b]];
1186   return VectorSpaceZBasis(B);
1187end intrinsic;
1188
1189intrinsic HeckeSpan(M::ModTupFld, V::ModTupFld) -> ModTupFld
1190{This function computes
1191 the Hecke module generated by a basis of V.
1192 The result is returned as a subspace
1193   with Z-basis.}
1194
1195   S := [SparseRepresentation(v) : v in Basis(V)];
1196   b := HeckeBound(M);
1197   H := [Heilbronn(n) : n in [1..b]];
1198   B := &cat[[TnSparse(M,  H[n], s) : n in [1..b]] : s in S];
1199   return VectorSpaceZBasis(B);
1200end intrinsic;
1201
1202
1203/******************************************************
1204 * Computation of T_p on dual space.                  *
1205 ******************************************************/
1206CFastData := recformat< V, e, VEinv>;
1207
1208intrinsic FastTnData(M::ModTupFld, V::SeqEnum) -> Rec
1209{Compute data needed by FastTn.}
1210
1211   // Step 1: find a set of very simple vectors e[1],...,e[n]
1212   // in M such that det(v[i]*e[j]) =/= 0.
1213   // 1. [Compute e] e = set of positions so that elements of the space
1214   // spanned by the columns of V are determined by the entries in these
1215   // spots.  This is accomplished by row reducing, and setting e equal to
1216   // the sequence of column numbers for which the columns are pivotal.
1217
1218   n := #V;
1219   B := Basis(sub<M|V>);
1220   assert #B eq n;
1221   // Find pivot columns.
1222   e := Pivots(B);
1223
1224   // Step 2: Compute the matrix V*E of inner products.
1225   VE    := RMatrixSpace(MF,n,n)![V[i][e[j]] : i in [1..n], j in [1..n]];
1226   VEinv := VE^(-1);
1227
1228   return rec<CFastData| V:=V, e:=e, VEinv:=VEinv>;
1229end intrinsic;
1230
1231intrinsic FastTnData(M::ModTupFld, V::ModTupFld) -> Rec
1232{}
1233   return FastTnData(M, Basis(V));
1234end intrinsic;
1235
1236intrinsic FastTn(M::ModTupFld, FastData::Rec, n::RngIntElt)
1237                   -> AlgMatElt
1238{}
1239   return FastTn(M,FastData,n,Heilbronn(n));
1240end intrinsic;
1241
1242intrinsic FastTn(M::ModTupFld, FastData::Rec, n::RngIntElt, H::SeqEnum)
1243                   -> AlgMatElt
1244{Compute action of Transpose(Tn) on the Hecke-stable
1245 subspace V.}
1246   F     := MF;
1247   V     := FastDataV;
1248   n     := #V;
1249   m     := Dimension(M);
1250   e     := FastDatae;
1251   VEinv := FastDataVEinv;
1252   // The next step is where all of the time is spent.
1253   TE    := [TnSparse(M, H, [<1,e[i]>]) : i in [1..n]];
1254   Vmat  := RMatrixSpace(F, n, m) ! V;
1255   TEmat := RMatrixSpace(F, n, m) ! TE;
1256   return  MatrixAlgebra(F,n)!(Vmat*Transpose(TEmat)*VEinv);
1257end intrinsic;
1258
1259intrinsic FastTn(FastData::Rec, n::RngIntElt, M::ModTupFld)
1260                -> AlgMatElt
1261{}
1262   return FastTn(M,FastData,n,Heilbronn(n));
1263end intrinsic;
1264
1265intrinsic FastTn(M::ModTupFld, V::ModTupFld, n::RngIntElt)
1266                   -> AlgMatElt, Rec
1267{Compute action of Transpose(Tn) on the Hecke-stable
1268 subspace V.}
1269   FastData := FastTnData(M, V);
1270   return FastTn(M, FastData, n);
1271end intrinsic;
1272
1273intrinsic FastTn(M::ModTupFld, V::ModTupFld, n::RngIntElt, H::SeqEnum)
1274                   -> AlgMatElt, Rec
1275{Compute action of Transpose(Tn) on the Hecke-stable
1276    subspace V.}
1277   FastData := FastTnData(V,M);
1278   return FastTn(M, FastData, n, H);
1279end intrinsic;
1280
1281intrinsic FastTn( M::ModTupFld, V::ModTupFld, nlist::SeqEnum,
1282                     H::SeqEnum)
1283                   -> SeqEnum, Rec
1284   {Compute action of Transpose(Tn), n in plist, on the Hecke-stable
1285    subspace V.}
1286   FastData := FastTnData(M, V);
1287   return [FastTn(M, FastData, n, H[n]): n in nlist];
1288end intrinsic;
1289
1290intrinsic FastTn(M::ModTupFld, V::ModTupFld, nlist::SeqEnum)
1291                   -> SeqEnum, Rec
1292{Compute action of Transpose(Tn), n in nlist, on the Hecke-stable
1293    subspace V.}
1294   FastData := FastTnData(M,V);
1295   return [FastTn(M, FastData, n): n in nlist];
1296end intrinsic;
1297
1298
1299
1300/******************************************************
1301 Action of T_n on modular symbols basis.
1302 ******************************************************/
1303
1304function ActionOnModularSymbolsBasis(g, M)
1305   // 1. Compute basis of modular symbols for M.
1306   B  := ModularSymbolsBasis(M);
1307   // 2. Apply g to each basis element.
1308   gB := [ModularSymbolApply(M, g,B[i]) : i in [1..#B]];
1309   // 3. Map the result back to M.
1310   gM := [ConvFromModularSymbol(M,gB[i]) : i in [1..#gB]];
1311   return MatrixAlgebra(MF,Dimension(M))!gM;
1312end function;
1313
1314
1315
1316// Compute the action of the Atkin-Lehner involution on M,
1317// when this makes sense.
1318intrinsic Wq(M::ModTupFld, q::RngIntElt) -> AlgMatElt
1319{Compute the action of the Atkin-Lehner involution Wq on M,
1320    when this makes sense (i.e., trivial character, even weight).
1321    The Atkin-Lehner map Wq is rescaled so that it is an
1322    involution, except when k>2 and char(F) divides q.}
1323   if GetCategory(M) eq ModSymFac then
1324      return Restrict(Transpose(Wq(ModSymParent(M),q)),M);
1325   end if;
1326   N := MN;
1327   k := Mk;
1328   assert q ne 0;
1329   assert N mod q eq 0;
1330   assert k mod 2 eq 0;
1331   assert IsTrivial(Meps);
1332   repeat
1333      d := Gcd(Integers()!(N/q),q);
1334      if d eq 1 then break; end if;
1335      q *:= d;
1336   until false;
1337   assert (N mod q) eq 0;
1338   // 1. Compute matrix to act by.
1339   d, x, y:= XGCD(q,-Integers()!(N/q));
1340   g := [q*x, y, N, q];
1341   A := ActionOnModularSymbolsBasis(g, M);
1342   p := Characteristic(MF);
1343   if p eq 0 or q mod p ne 0 then
1344      A /:= q^(Integers()!(k/2)-1);
1345   end if;
1346   return A;
1347end intrinsic;
1348
1349intrinsic StarInvolution(M::ModTupFld) -> AlgMatElt
1350{Compute the conjugation involution * on V.  This involution
1351is defined by the 2x2 matrix [-1,0,0,1]; it sends
1352X^i*Y^j\{u,v\} to (-1)^j*X^i*Y^j \{-u,-v\}.}
1353   return ActionOnModularSymbolsBasis([-1,0,0,1], M);
1354end intrinsic;
1355
1356intrinsic Tnpoly(n::RngIntElt, M::ModTupFld) -> SeqEnum
1357{Compute the characteristic polynomial of the p-th Hecke
1358    operator Tn on the space M of modular symbols.  Uses the modular
1359    algorithm, without proof.}
1360   return ModularCharpoly(HeckeOperator(n,M));
1361end intrinsic;
1362
1363intrinsic Restrict (A::AlgMatElt, L::Lat) -> AlgMatElt
1364{  Suppose A is a linear transformation of V
1365   which leaves the lattice L invariant.
1366   Returns A restricted to W, with respect to
1367   the basis Basis(W).}
1368   Z := Basis(L);
1369   V := VectorSpace(Parent(A[1,1]), Dimension(L));
1370   B := [V!Z[i] : i in [1..Dimension(V)]];
1371   return Restrict (A, B);
1372end intrinsic;
1373
1374intrinsic Restrict (A::AlgMatElt, W::ModTupFld) -> AlgMatElt
1375{  Suppose A is a linear transformation of V
1376   which leaves the subspace  W invariant.
1377   Returns A restricted to W, with respect to
1378   the basis Basis(W). }
1379   B := Basis(W);
1380   return Restrict (A, B);
1381end intrinsic;
1382
1383intrinsic Restrict (A::AlgMatElt, W::ModTupRng) -> AlgMatElt
1384{  Suppose A leaves W invariant.
1385   Returns A restricted to W, with respect to
1386   the basis Basis(W).}
1387   B := Basis(W);
1388   return Restrict (A, B);
1389end intrinsic;
1390
1391intrinsic Restrict (A::AlgMatElt, B::SeqEnum) -> AlgMatElt
1392{  Suppose A is a linear transformation of V
1393   which leaves the subspace W spanned by the
1394   vectors listed in B invariant.
1395   Returns A restricted to W, with respect to
1396   the basis B.}
1397   S := RSpaceWithBasis(B);
1398   return MatrixAlgebra(Parent(A[1,1]),#B) !
1399           &cat[Coordinates(S, S.i*A) : i in [1..#B]];
1400end intrinsic;
1401
1402
1403intrinsic KernelOn (A::AlgMatElt, B::SeqEnum) -> ModTupFld
1404{Suppose B is a basis for an n-dimensional subspace
1405    of some ambient space and A is an nxn matrix.
1406    Then A defines a linear transformation of the space
1407    spanned by B.  This function returns the
1408    kernel of that transformation.}
1409   n := #B;
1410   if n eq 0 then
1411      return B;
1412   end if;
1413   assert Nrows(A) eq Ncols(A) and Nrows(A) eq n;
1414   K := Kernel(A);
1415   return VectorSpaceWithBasis(LinearCombinations(Basis(K),B));
1416end intrinsic;
1417
1418intrinsic KernelOn (A::AlgMatElt, W::ModTupFld) -> ModTupFld
1419{}
1420   return KernelOn (A, Basis(W));
1421end intrinsic;
1422
1423intrinsic HeckeOperatorOn(M::ModTupFld, V::ModTupFld, n::RngIntElt)
1424                 -> AlgMatElt
1425  {Compute the action of the Hecke operator Tn with
1426   respect to the basis for the subspace V of M.}
1427  return Restrict(HeckeOperator(M,n), V);
1428end intrinsic;
1429
1430intrinsic HeckeOperatorMkZ(M::ModTupFld, n::RngIntElt) -> AlgMatElt
1431  {Compute the action of Tn on the integral modular symbols.}
1432  return Restrict(HeckeOperator(M,n), IntegralModularSymbols(M));
1433end intrinsic;
1434
1435intrinsic HeckeOperatorSkZ(M::ModTupFld, n::RngIntElt) -> AlgMatElt
1436  {Compute the action of T_p on a subspace of integral modular symbols.}
1437  return Restrict(HeckeOperator(M,n), SkZ(M));
1438end intrinsic;
1439
1440
1441
1442/**********************************************************
1443 *                                                        *
1444 * BOUNDARY SYMBOLS:  Compute the boundary map delta.     *
1445 *                                                        *
1446 **********************************************************/
1447
1448
1449///////////////////////////////////////////////////////////////////////////
1450//////////// IMPLEMENTATION 2: CUSPS FOR NONTRIVIAL CHARACTER /////////////
1451///////////////////////////////////////////////////////////////////////////
1452function ReduceCusp(a)
1453   d := Gcd(a);
1454   return [Integers()|x/d : x in a];
1455end function;
1456
1457function CuspEquiv(N,a,b)
1458   u1, v1 := Explode(ReduceCusp(a));
1459   u2, v2 := Explode(ReduceCusp(b));
1460   s1 := 0;
1461   s2 := 0;
1462   if [u1,v1] ne [0,1] then
1463      s1 := (v1 eq 0 or v1 eq 1) select 1 else InverseMod(u1,Abs(v1));
1464   end if;
1465   if [u2,v2] ne [0,1] then
1466      s2 := (v2 eq 0 or v2 eq 1) select 1 else InverseMod(u2,Abs(v2));
1467   end if;
1468   g := Gcd(v1*v2,N);
1469   a := s1*v2 - s2*v1;
1470   if a mod g ne 0 then
1471      return false, 1;
1472   end if;
1473
1474   // Now we know the cusps are equivalent.  Use the proof of Prop2.2.3
1475   // of Cremona to find a matrix in Gamma_0(N) relating them.
1476   dum,s2,r2 := Xgcd(u2,-v2);
1477   assert dum eq 1;
1478   dum,s1,r1 := Xgcd(u1,-v1);
1479   assert dum eq 1;
1480   a := s1*v2 - s2*v1;
1481   assert a mod g eq 0;
1482   // solve x*v1*v2 + a = 0 (mod N).
1483   d,x0,y0 := Xgcd(v1*v2,N);          // x0*v1*v2 + y0*N = d = g.
1484   // so x0*v1*v2 - g = 0 (mod N)
1485   x := -x0 * Integers()!(a/g);
1486   // now  x*v1*v2 + a = 0 (mod N)
1487   s1p := s1+x*v1;
1488   return true, (u2*s1p-r2*v1) mod N;
1489end function;
1490
1491function CuspToFreeHelper(M, sign, a)
1492   if not assigned Mcusplist then
1493      Mcusplist := [];
1494   end if;
1495
1496   list  := Mcusplist;
1497   F     := MF;
1498   eps   := Meps;
1499   N     := MN;
1500   k     := Mk;
1501
1502   a := ReduceCusp(a);
1503   if a[2] lt 0 then
1504      a[1] *:= F!-1;
1505      a[2] *:= F!-1;
1506   end if;
1507
1508   // Check if we've already encountered this cusp.
1509   for i in [1..#list] do
1510      b          := list[i];
1511      equiv, alp := CuspEquiv(N, b[1], a);   // [gam_alp(b[1])]=?=[a].
1512      if equiv then
1513         if b[2] eq 0 then
1514            return <F!0,1>;
1515         end if;
1516         if IsTrivial(Character(M)) then
1517            return <1,i>;
1518         else
1519            return <Evaluate(eps,alp)^(-1),i>;
1520         end if;
1521      end if;
1522      if sign ne 0 then
1523         equiv, alp := CuspEquiv(N, b[1], [-a[1],a[2]]);
1524         if equiv then
1525            if b[2] eq 0 then
1526               return <F!0,1>;
1527            end if;
1528            if IsTrivial(Character(M)) then
1529               return <sign,i>;
1530            else
1531               return <sign*Evaluate(eps,alp),i>;
1532            end if;
1533         end if;
1534      end if;
1535   end for;
1536
1537   // Determine if this cusp class is killed by the relations.
1538   c := F!1;
1539   if not IsTrivial(Character(M)) then
1540      u := a[1]; v := a[2];
1541      g := Gcd(N,v);
1542      x := Integers()!(N/Gcd(N,v));
1543      for j in [0..g-1] do
1544         alp := 1-j*x;
1545         if Gcd(alp,N) eq 1 then
1546            if (v*(1-alp)) mod N eq 0 and (u*(1-alp)) mod g eq 0 then
1547               if Evaluate(eps,alp) ne 1 then
1548                  c := F!0;
1549                  break;
1550               end if;
1551            end if;
1552         end if;
1553      end for;
1554   end if;
1555
1556   Append(~list,<a,c>);
1557   Mcusplist := list;
1558   i := #list;
1559
1560   if sign ne 0 then
1561      // Check is sign relation kills this cusp.
1562      ans := CuspToFreeHelper(M,0,[-a[1],a[2]]);
1563      if ans[2] eq i and ans[1] ne sign then
1564         Mcusplist[i][2] := 0;
1565         c := F!0;
1566      end if;
1567   end if;
1568
1569   return <c,i>;
1570end function;
1571
1572intrinsic CuspToFree(M::ModTupFld, a::SeqEnum) -> Tup
1573{Given a cusp a=[u,v] this function finds the index i and a
1574 scalar alpha such that a = alpha*(i-th standard cusp).
1575 It then returns alpha, i.}
1576   return CuspToFreeHelper(M,Sign(M),a);
1577end intrinsic;
1578
1579forward LiftToCosetRep;
1580intrinsic BoundaryMap(M::ModTupFld) -> ModMatFldElt
1581{Compute the Boundary map.}
1582   if not assigned MBoundaryMap then
1583      if Dimension(M) eq 0 then
1584         MBoundaryMap := RMatrixSpace(MF,0,0)!0;
1585         return MBoundaryMap;
1586      end if;
1587      Tgens := MquotTgens;
1588      Sgens := MquotSgens;
1589      F     := MF;
1590      n     := #Tgens;
1591      D     := [ [] : i in [1..n]];
1592      for j in [1..n] do
1593         i := Sgens[Tgens[j]];
1594         uv, w, ind := UnwindManinSymbol(i,Mmlist);
1595         g := LiftToCosetRep(uv, MN);
1596         if w eq Mk-2 then
1597            Append(~D[j], CuspToFree(M,[g[1],g[3]]));
1598         end if;
1599         if w eq 0 then // substract the second term.
1600            t := CuspToFree(M,[g[2],g[4]]);
1601            t[1] *:= -1;
1602            Append(~D[j], t);
1603         end if;
1604      end for;
1605      if &cat D eq [] then
1606         MBoundaryMap := RMatrixSpace(MF,Dimension(M),0)!0;
1607         return MBoundaryMap;
1608      end if;
1609      m := Max([x[2] : x in &cat D]);
1610      V := VectorSpace(F,m);
1611      Drows := [&+[V| x[1]*V.x[2] : x in D[i]] : i in [1..n]];
1612      MBoundaryMap := RMatrixSpace(MF,n,m) ! Drows;
1613   end if;
1614   return MBoundaryMap;
1615end intrinsic;
1616
1617
1618
1619
1620
1621
1622
1623intrinsic CuspidalSymbols(M::ModTupFld) -> ModTupFld
1624{Return the subspace of S_k(N,Q) of cuspidal modular symbols.}
1625   if not assigned MSkF then
1626      MSkF := Kernel(BoundaryMap(M));
1627   end if;
1628   return MSkF;
1629end intrinsic;
1630
1631intrinsic CuspidalSymbols(M::ModTupFld) -> ModTupFld
1632{Return the subspace of S_k(N,Q) of cuspidal modular symbols.}
1633   if not assigned MSkF then
1634      D := BoundaryMap(M);
1635      if HeckeSuperverbose then
1636         "Computing kernel of boundary map: ";
1637         t := Cputime();
1638      end if;
1639      MSkF := Kernel(D);
1640      if HeckeSuperverbose then
1641         Cputime(t), " seconds.";
1642      end if;
1643   end if;
1644   return MSkF;
1645end intrinsic;
1646
1647intrinsic Sk(M::ModTupFld) -> ModTupFld
1648{}
1649   return CuspidalSymbols(M);
1650end intrinsic;
1651
1652intrinsic IntegralCuspidalSymbols(M::ModTupFld) -> ModTupFld
1653{}
1654   if not assigned MSkZ then
1655      Z     := IntegralModularSymbols(M);
1656      D     := BoundaryMap(M);
1657      DofZ  := [v*D : v in Basis(Z)];
1658      DZ    := RMatrixSpace(MF,Nrows(D),Ncols(D)) ! DofZ;
1659      MSkZ := VectorSpaceWithBasis(
1660                  LinearCombinations(Basis(IntegerKernel(DZ)), Basis(Z))
1661               );
1662   end if;
1663   return MSkZ;
1664end intrinsic;
1665
1666
1667intrinsic SkZ(M::ModTupFld) -> ModTupFld
1668{}
1669   case GetCategory(M):
1670      when ModSym: return IntegralCuspidalSymbols(M);
1671      when ModSymFac: return DecompZ(M);
1672      else error "SkZ only takes a ModSym or ModSymFac as argument.";
1673   end case;
1674end intrinsic;
1675
1676intrinsic SkZPlus(M::ModTupFld) -> ModTupFld
1677{}
1678   if not assigned MSkZplus then
1679      S    := SkZ(M);
1680      star := Restrict(StarInvolution(M),S);
1681      MSkZplus := IntegerKernelOn(star-1, S);
1682   end if;
1683   return MSkZplus;
1684end intrinsic;
1685
1686intrinsic SkZMinus(M::ModTupFld) -> ModTupFld
1687{}
1688   if not assigned MSkZminus then
1689      S    := SkZ(M);
1690      star := Restrict(StarInvolution(M),S);
1691      MSkZminus := IntegerKernelOn(star+1, S);
1692   end if;
1693   return MSkZminus;
1694end intrinsic;
1695
1696/*
1697function EisensteinSymbols(M)
1698end function;
1699*/
1700
1701
1702/*******************************************************************
1703  DEGENERACY COSET REPRESENTATIVES:
1704    Let N be a positive integer and M a divisor of N.
1705    Let t be a divisor of N/M, and let T be the 2x2 matrix T=[0,0,  0,t].
1706    Find coset reps for Gamma_0(N) \ T Gamma_0(M).
1707
1708  FACTS: T^(-1)*(a,b,c,d)*T = (a,bt,c/t,d)
1709         T^(-1)*Gamma_0(N)*T is contained in Gamma_0(M)
1710         Gamma_0(N)*T is contained in T*Gamma_0(M)
1711
1712
1713  SOLUTION:
1714    (1) Computes representatives for Gamma_0(N/t,t) inside
1715        of Gamma_0(M):
1716        COSET EQUIVALENCE:
1717           Two right cosets represented by (a,b;c,d) &
1718           (a',b';c',d') of Gamma_0(N/t,t) in SL_2(Z) are equivalent iff
1719           (a,b)=(a',b') as points of P^1(t), i.e., ab'=ba'(mod t),
1720           and (c,d)=(c',d') as points of P^1(N/t).
1721        ALGORITHM:
1722           (A) Compute the number of cosets.
1723           (B) Compute a random element x of Gamma_0(M).
1724           (C) Check if x is equivalent to anything generated so
1725               far; if not, add x to the list.
1726           (D) Continue until the list is as long as the bound
1727               computed in A.
1728
1729    (2) There is a natural bijection
1730         Gamma_0(N)\T*Gamma_0(M) <---> Gamma_0(N/t,t)\Gamma_0(M)
1731        given by
1732            Tr  <---------> r
1733        Consequently we obtain coset representatives for
1734        Gamma_0(N)\T*Gamma_0(M) by left multiplying by T each
1735        coset representative of Gamma_0(N/t,t)\Gamma_0(M) found
1736        in step 1.
1737 ********************************************************************/
1738intrinsic DegeneracyCosetReps(M, N, d) -> .
1739{}
1740   n       := idxG0(N)/idxG0(M);      // number of coset representatives.
1741   Ndivd   := N div d;
1742   R       := [];                     // coset reps found so far.
1743   max     := 4*(n+10);
1744   halfmax := Round(max/2);
1745   while #R lt n do
1746      // try to find another coset representative.
1747      cc := M*Random(-halfmax, halfmax);
1748      dd :=   Random(-halfmax, halfmax);
1749      g, bb, aa := XGCD(-cc,dd);
1750      if g eq 0 then continue; end if;
1751      cc div:= g;
1752      dd div:= g;
1753      if cc mod M ne 0 then continue; end if;
1754      // Test if we've found a new coset representative.
1755      isnew:=true;
1756      for r in R do
1757           if ((r[2]*aa - r[1]*bb) mod d eq 0) and
1758               ((r[4]*cc - r[3]*dd) mod Ndivd eq 0) then
1759            isnew := false;
1760            break ;
1761         end if;
1762      end for;
1763      // If it is new add it to the list.
1764      if isnew then
1765         Append(~R,[aa,bb,cc,dd]);
1766      end if;
1767   end while;
1768   // Return the list left multiplied by T.
1769   return [[r[1],r[2],r[3]*d,r[4]*d] : r in R];
1770end intrinsic;
1771
1772
1773/**********************************************************
1774  DEGENERACY AND INCLUSION MAPS -- newforms.
1775  If possible compute the map M1 --> M2 corresponding
1776  to the integer d.
1777
1778  WARNING: There might be a problem here when
1779  the characteristic of the base field divides d,
1780  because d occurs in the denominators of the
1781  DegeneracyReps.
1782
1783 d:=1;c:=3;M:=11;N:=M*d*c;
1784 A:=Mk(M,2); B:=Mk(N,2);
1785 bet:=dm(A,B,d); alp:=dm(B,A,d);
1786 bet*alp;      // must be a scalar.
1787
1788 **********************************************************/
1789intrinsic DegeneracyMap(M1::ModTupFld, M2::ModTupFld)
1790                         -> AlgMatElt
1791{If possible compute the map M1 --> M2 corresponding
1792 to the integer d.}
1793   return DegeneracyMap(M1,M2,1);
1794end intrinsic;
1795
1796intrinsic DegeneracyMap(M1::ModTupFld, M2::ModTupFld, d::RngIntElt)
1797                         -> AlgMatElt
1798{If possible compute the degeneracy map M1 --> M2 corresponding
1799to the integer d.  Here M1 and M2 are spaces of modular symbols,
1800of the same weight and d is a divisor of the quotient of
1801the two levels N1 and N2.  If the N1 divides N2 then the
1802"raising" map is computed, if N1 is divisible by N2,
1803then the "lowering" map is computed.}
1804   assert M1F eq M2F;
1805   assert Weight(M1) eq Weight(M2);
1806   N1 := Level(M1);
1807   N2 := Level(M2);
1808   if N1 eq N2 then
1809      assert d eq 1;
1810      assert M1 eq M2; // we could write down the map between
1811                       // different presentations of M_k(N), but
1812                       // we won't for now, since there should
1813                       // never be a need.
1814      F := BaseField(M1);
1815      return Hom(M1,M1)!MatrixAlgebra(F,Degree(M1))!1;
1816   end if;
1817
1818   if N1 mod N2 eq 0 then  // N2 divides N1 -- lower
1819      assert (N1 div N2) mod d eq 0;
1820      B  := ModularSymbolsBasis(M1);
1821      D  := [d,0,0,1];
1822      DB := [ModularSymbolApply(M1,D, B[i])  : i in [1..#B] ];
1823      A  := [ConvFromModularSymbol(M2,DB[i]) : i in [1..#DB]];
1824      return Hom(M1, M2) ! A;
1825
1826   elif N2 mod N1 eq 0 then// N1 divides N2 -- raise level
1827      assert (N2 div N1) mod d eq 0;
1828      B   := ModularSymbolsBasis(M1);
1829      R   := DegeneracyCosetReps(N1, N2, d);
1830      eps := Character(M1);
1831      if IsTrivial(eps) then
1832         RB := [ &cat[ModularSymbolApply(M1, r, B[i]) : r in R]
1833                                                  : i in [1..#B]];
1834         A := [ConvFromModularSymbol(M2,RB[i]) : i in [1..#RB]];
1835      else
1836         A   := [ &+[Evaluate(eps,r)^(-1)*ConvFromModularSymbol(M2,
1837                  ModularSymbolApply(M1, r, B[i])) : r in R]
1838               : i in [1..#B]];
1839      end if;
1840      return Hom(M1, M2) ! A;
1841   else
1842      error "DegeneracyMap: N1 must be divisible by N2, or vice-versa.";
1843   end if;
1844end intrinsic;
1845
1846intrinsic OldModularSymbols(M::ModTupFld, NN::RngIntElt) -> ModTupFld
1847{Returns the space Mk(NN,eps') associated to Mk(N,eps).
1848Here NN must be a divisor of N.}
1849   N := Level(M);
1850   if NN eq N then return M; end if;
1851   divs := Divisors(N);
1852   if not assigned Mold then
1853      Mold := SequenceToList(divs);
1854   end if;
1855   pos := Index(divs,NN);
1856   if pos eq 0 then
1857      error "OldModularSymbols: NN =",NN,"is not a divisor of the level =",N;
1858   end if;
1859   if Type(Mold[pos]) eq RngIntElt then
1860
1861      /****************************************************************
1862         This object does not know the old modular symbols
1863         at level NN.  Here's what we do:
1864            (1) If this object has a parent, ask it for
1865                the modular symbols at level NN. If so, done;
1866                otherwise continue with the next step.
1867            (2) If not, construct the space of symbols at level NN.
1868            (3) If this object has a parent, tell it about the symbols
1869                at level NN.
1870       ****************************************************************/
1871
1872      // Step 1.
1873      if not ModSymIsRoot(M) then
1874         // Parent takes care of everything recursively.
1875         Mold[pos] := OldModularSymbols(ModSymParent(M),NN);
1876      else
1877         // Step 2: We've got to do the work.
1878         NN    := Integers()!NN;
1879         eps   := Character(M);
1880         F     := BaseField(M);
1881         if not CanReduceModM(eps,NN) then
1882            Mold[pos] := BaseExtend(TrivialSpace(Weight(M),eps,Sign(M)),F);
1883         else
1884            epsNN := AssociatedModMCharacter(eps,NN);
1885            Mold[pos] := BaseExtend(
1886                     ModularSymbols(epsNN,Weight(M),Sign(M)),F);
1887         end if;
1888         (Mold[pos])M := M;  // tag this as *our* child.
1889      end if;
1890   end if;
1891   return Mold[pos];
1892end intrinsic;
1893
1894/**********************************************************
1895   For each prime p dividing the level there are two maps
1896    alp_1,alp_p: Mk(N,e) ----> Mk(N/p,e')
1897   The intersection
1898         Kernel(alp_1) meet Kernel(alp_p)
1899   is called the p-new subspace Mkpnew(N) of Mk(N).
1900   If e doesn't factor through (Z/(N/p)Z)^* then
1901   e' is by definition 0.
1902 ********************************************************/
1903intrinsic pNewSubspace(M::ModTupFld, p::RngIntElt) -> ModTupFld
1904{Returns the p-new subspace of M.}
1905   assert IsPrime(p);
1906   N     := Level(M);
1907   k     := Weight(M);
1908   assert N mod p eq 0;
1909
1910   eps   := Character(M);
1911   NN    := Integers()!(N/p);
1912
1913   if not CanReduceModM(eps,NN) then
1914      return M;
1915   end if;
1916
1917   old  := OldModularSymbols(M,NN);
1918
1919   if Dimension(old) eq 0 then
1920      return M;
1921   end if;
1922   h1  := DegeneracyMap(M, old, 1);
1923   hp  := DegeneracyMap(M, old, p);
1924   return Kernel(h1) meet Kernel(hp);
1925end intrinsic;
1926
1927
1928intrinsic pNewDualSubspace(M::ModTupFld, p::RngIntElt) -> ModTupFld
1929{Returns the p-new subspace of M.}
1930   assert IsPrime(p);
1931   N     := Level(M);
1932   k     := Weight(M);
1933   assert N mod p eq 0;
1934
1935   eps   := Character(M);
1936   NN    := Integers()!(N/p);
1937
1938   if not CanReduceModM(eps,NN) then
1939      return M;
1940   end if;
1941
1942   old  := OldModularSymbols(M,NN);
1943
1944   h1  := DegeneracyMap(old, M, 1);
1945   hp  := DegeneracyMap(old, M, p);
1946   return Kernel(Transpose(h1)) meet Kernel(Transpose(hp));
1947end intrinsic;
1948
1949
1950
1951intrinsic Sknew(M::ModTupFld) -> ModTupFld
1952{}
1953   return CuspidalNewSubspace(M);
1954end intrinsic;
1955
1956intrinsic CuspidalNewSubspace(M::ModTupFld) -> ModTupFld
1957{}
1958   if not assigned Msknew then
1959      Msknew := NewSubspace(M) meet CuspidalSymbols(M);
1960   end if;
1961   return Msknew;
1962end intrinsic;
1963
1964intrinsic NewSubspace(M::ModTupFld) -> ModTupFld
1965{Returns the p-new subspace of M.}
1966   if not assigned Mmknew then
1967      if Level(M) eq 1 then
1968         Mmknew := M;
1969      else
1970         Mmknew := &meet[pNewSubspace(M,p[1])
1971                   : p in Factorization(Level(M))];
1972      end if;
1973   end if;
1974   return Mmknew;
1975end intrinsic;
1976
1977
1978
1979intrinsic Mknewdual(M::ModTupFld) -> ModTupFld
1980{}
1981   return NewDualSubspace(M);
1982end intrinsic;
1983
1984intrinsic Sknewdual(M::ModTupFld) -> ModTupFld
1985{}
1986   return CuspidalNewDualSubspace(M);
1987end intrinsic;
1988
1989intrinsic CuspidalNewDualSubspace(M::ModTupFld) -> ModTupFld
1990{}
1991   if not assigned Msknewdual then
1992      Msknewdual := NewDualSubspace(M);
1993      // cut out the cuspidal part using Hecke operators.
1994      S := Sknew(M);
1995      if Dimension(S) eq 0 then
1996         Msknewdual := S;
1997      else
1998         p := 2;
1999         while Dimension(Msknewdual) gt Dimension(S) do
2000            T := HeckeOperator(M,p);
2001            f := ModularCharpoly(Restrict(T,S));
2002            Ton := Restrict(Evaluate(f,Transpose(T)),Msknewdual);
2003            Msknewdual := KernelOn(Ton,Msknewdual);
2004            p := NextPrime(p);
2005         end while;
2006      end if;
2007   end if;
2008   return Msknewdual;
2009end intrinsic;
2010
2011intrinsic NewDualSubspace(M::ModTupFld) -> ModTupFld
2012{Returns the p-new subspace of