freeze;
 
/****-*-magma-* EXPORT DATE: 2004-03-08 ************************************
                                                                            
                     MODFORM: Modular Forms in MAGMA
 
                              William A. Stein        
                         
   FILE: operators.m
   
   01/04/04 (WAS) - fixed bug in the require statement in 
          "intrinsic AtkinLehnerOperator(M::ModFrm, q::RngIntElt)"
       
     

   04/06/03 (WAS) - added caching of eigenvector for T2 when computing level 1 charpolys.

   04/06/03 (WAS) - found bug in charpoly routine: In a part of HeckePolynomial it 
                    was finding charpoly of Tn instead of T2 which led to a hang.

   04/05/03 (WAS) - fixed bug in HeckeOperator (it said BaseField instead
                    of BaseRing at a certain point).

   $Header: /home/was/magma/packages/ModFrm/code/RCS/operators.m,v 1.9 2002/10/20 06:30:02 was Exp was $

   $Log: operators.m,v $
   Revision 1.9  2002/10/20 06:30:02  was
   Fixed serious bug in function WqOnNewSpace(M, q, prec) that meant Wq was
   computed wrt wrong basis.

   Revision 1.8  2002/05/04 18:33:24  was
   Made AtkinLehnerOperator much much more efficient by not computing T_p.

   Revision 1.7  2002/04/25 05:34:11  was
   some improvements to HeckePolynomial in the level 1 case.

   Revision 1.6  2002/04/25 00:21:40  was
   *** empty log message ***

   Revision 1.5  2001/07/13 02:33:01  was
   Changed a name.

   Revision 1.4  2001/06/08 06:49:36  was
   Added AtkinLehnerEigenvalue

   Revision 1.3  2001/05/30 18:56:17  was
   Created.

   Revision 1.2  2001/05/16 04:11:40  was
   *** empty log message ***

   Revision 1.1  2001/05/16 03:52:06  was
   Initial revision

      
 ***************************************************************************/

forward 
   ApplyHeckeOperator, 
   ComputeHeckePolynomialUsingHeckeOperator,
   Tn, Tp;

import "eisenstein.m":
   AtkinLehnerOnEisensteinModularFormsSpace;

import "level1.m": 
   Level1CharacteristicPolynomialOfTp,
   EigenvectorOfMatrixWithCharpoly;

import "modular_symbols.m": 
   MF_ModularSymbols;

import "predicates.m": 
   AssociatedSpaceOverZ, 
   SpaceType, 
   SpaceTypeParam;

intrinsic HeckeOperator(M::ModFrm, n::RngIntElt) -> AlgMatElt
{Matrix representing the n-th Hecke operator on M.}
   require IsGamma0(M) : 
     "Hecke operator computation currently only supported for Gamma_0(N).";
   requirege n, 1;
   if not assigned M`hecke_operator then
      M`hecke_operator := [];
   end if;
   if exists(i) { i : i in [1..#M`hecke_operator]
                             | M`hecke_operator[i][1] eq n } then
      return M`hecke_operator[i][2];
   end if;
   if n eq 1 then
      Tn := MatrixAlgebra(BaseRing(M),Dimension(M))!1;
   else
      Tn := MatrixAlgebra(BaseRing(M),Dimension(M))!
          &cat[Eltseq(ApplyHeckeOperator(n,M.i)) : i in [1..Dimension(M)]];
   end if;
   Append(~M`hecke_operator,<n,Tn>);
   return Tn;
end intrinsic;

function ApplyHeckeOperator(n, f)
   // The image of f under the n-th Hecke operator.
   assert Type(n) eq RngIntElt;
   assert Type(f) eq ModFrmElt;
   assert IsGamma0(Parent(f));
   if n eq 1 then
      return f;
   end if;
   prec := (n+1)*PrecisionBound(AmbientSpace(Parent(f)));
   Tnf := Tn(Level(f), DirichletGroup(Level(f),Integers())!1, Weight(f), n,
               PowerSeries(f, prec));
   return Parent(f)!Tnf;
end function;

intrinsic AtkinLehnerEigenvalue(f::ModFrmElt, q::RngIntElt) -> RngElt
{The eigenvalue of the Atkin-Lehner involution W_q on the cuspidal newform f.}
   if assigned f`elliptic_curve and not assigned f`mf_modular_symbols then
      require false : "Not yet programmed for elliptic curve newforms.";
   end if;
   require assigned f`is_newform and f`is_newform and 
             assigned f`mf_modular_symbols : "f must be a cuspidal newform.";
   Wq := AtkinLehner(f`mf_modular_symbols,q);
   assert IsDiagonal(Wq);
   return Wq[1,1];
end intrinsic;


function CoerceUnivariatePolynomialElement(R, f)
   assert Type(R) eq RngUPol;
   assert Type(f) eq RngUPolElt;
   return &+[R.1^i*(R!(Integers()!Coefficient(f,i))) : i in [0..Degree(f)]];
end function;


function GaloisConjugatesOfPolynomial(f)
   assert Type(f) eq RngUPolElt;
   F := BaseRing(Parent(f));
   if Degree(F) eq 1 then
      return [f];
   end if;
   assert Type(F) eq FldCyc;
   X := Parent(f).1;
   ans := [f];
   m := CyclotomicOrder(F);
   for n in [2..EulerPhi(m)] do
      Append(~ans,&+[Conjugate(Coefficient(f,j),n)*X^j : j in [0..Degree(f)]]); 
   end for;
   return ans;
end function;

function ComputeHeckePolynomialOfModularSymbolsSpace(M, modsym, n, Proof)
   assert Type(M) eq ModFrm;
   assert Type(modsym) eq ModSym;
   assert Type(n) eq RngIntElt and n ge 1;
   assert Type(Proof) eq BoolElt;

   if Dimension(M) eq 0 then
      return PolynomialRing(BaseRing(M))!1;
   end if;
   f := HeckePolynomial(modsym, n : Proof := Proof);
   normf := &*GaloisConjugatesOfPolynomial(f);
   return normf;

end function;

function CharpolyOfEigenvalue(a)
   if Type(Parent(a)) eq FldRat then
      return PolynomialRing(RationalField()).1 - a;
   end if;
   return CharacteristicPolynomial(a);
end function;

function ComputeHeckePolynomialUsingEisensteinSeries(M,n)
   assert Type(M) eq ModFrm;
   assert IsEisenstein(M);
   assert GCD(Level(M),n) eq 1;

   if Dimension(M) eq 0 then
      return PolynomialRing(BaseRing(M))!1;
   end if;

// This is incorrect, because the n-th coefficient of the 
// Eisenstein series need not be the Eigenvalue of T_n,
// except in some special cases.   
   print "Computing Hecke polynomial on Eisenstein series not yet implemented. (Returning 1.)";
   R := PolynomialRing(CyclotomicField(EulerPhi(Level(M))));
   return R!1;

   E := EisensteinSeries(M);
   F := &*[R|R.1 - R!Coefficient(PowerSeries(f,n+1),n) : f in E];
   return F;
end function;


function ComputeHeckePolynomialUsingModularSymbols(M,n, Proof)
   assert Type(M) eq ModFrm;
   assert Type(n) eq RngIntElt;
   assert Type(Proof) eq BoolElt;

   if Dimension(M) eq 0 then
      return PolynomialRing(BaseRing(M))!1;
   end if;
   if IsEisenstein(M) then
      if IsNew(M) then
         if IsGamma0(M) then
             return ComputeHeckePolynomialUsingHeckeOperator(M,n,Proof);
         else
             error "Not yet programmed -- I haven't coded an algorithm to compute "*
               "the char poly of T_n on the new Eisenstein series when n divides "*
               "the level (and we're not in the Gamma_0(N) case).";

         end if;
      end if;
      modsym := MF_ModularSymbols(M,0);
      return &*[ComputeHeckePolynomialOfModularSymbolsSpace(M,m,n,Proof) 
                  : m in modsym];
   elif IsCuspidal(M) then
      modsym := MF_ModularSymbols(M,+1);
      return &*[ComputeHeckePolynomialOfModularSymbolsSpace(M,m,n,Proof) 
                  : m in modsym];
   end if;

   return ComputeHeckePolynomialUsingModularSymbols(
                             EisensteinSubspace(M),n,Proof) * 
                ComputeHeckePolynomialUsingModularSymbols(
                             CuspidalSubspace(M),n,Proof) ;
end function;



function ComputeHeckePolynomialUsingHeckeOperator(M, n, Proof)
   assert Type(M) eq ModFrm;
   assert Type(n) eq RngIntElt;
   assert Type(Proof) eq BoolElt;
   Tn := HeckeOperator(M,n);
   if Type(BaseRing(M)) in {RngInt, FldRat} then
      return CharacteristicPolynomial(Tn : Al := "Modular", Proof := Proof);
   end if;
   return CharacteristicPolynomial(Tn);
end function;


intrinsic HeckePolynomial(M::ModFrm, n::RngIntElt :
                  Proof := false) -> AlgMatElt
{The characteristic polynomial of the n-th Hecke operator on M.}
   require n ge 1 : "Argument 2 must be at least 1.";
   if not assigned M`hecke_polynomial then
      M`hecke_polynomial := [];   // sequence of pairs <n,f(x)>
   end if;
 
   if exists(i) { i : i in [1..#M`hecke_polynomial]
                             | M`hecke_polynomial[i][1] eq n } then
      return M`hecke_polynomial[i][2];
   end if;
 
   vprint ModularForms, 2: "Computing Hecke polynomial of index", n;
   if n eq 1 then
      f := (PolynomialRing(BaseRing(M)).1-1)^Dimension(M);
   elif assigned M`made_from_newform then
      f := CharpolyOfEigenvalue(
              Coefficient(PowerSeries(M`made_from_newform,n+1),n));
   elif Characteristic(BaseRing(M)) eq 0 and Level(M) eq 1 and n le Dimension(CuspidalSubspace(M)) and Proof eq false then
      A := AssociatedSpaceOverZ(CuspidalSubspace(M));
      assert Dimension(A) eq Dimension(CuspidalSubspace(M));
      d := Dimension(A);
      T2 := HeckeOperator(A,2);   // yes, T2, because my algorithm uses T2!!
      if exists(i) { i : i in [1..#A`hecke_polynomial]
                             | A`hecke_polynomial[i][1] eq n } then
          f2 := A`hecke_polynomial[i][2];
      else
         vprint ModularForms, 2: 
           "Computing characteristic polynomial of T_2 using standard algorithm.";
          t := Cputime();
          f2  := CharacteristicPolynomial(T2 : Proof := false, Al := "Modular");
          vprintf ModularForms, 2: "Time = %o seconds\n", Cputime(t);
          Append(~A`hecke_polynomial, <2, f2>);
      end if;

      if not assigned A`t2_eigenvector then
         T2Q := MatrixAlgebra(RationalField(),Nrows(T2))!Eltseq(T2);
         f2Q := PolynomialRing(RationalField())!f2;
         A`t2_eigenvector := EigenvectorOfMatrixWithCharpoly(T2Q, f2Q);
      end if;
      
      if n gt 2 then
         cuspidal_part := Level1CharacteristicPolynomialOfTp(Weight(A),n,d,A`t2_eigenvector);
      else
         cuspidal_part := f2;
      end if;        
      eisenstein_part := ComputeHeckePolynomialUsingEisensteinSeries(EisensteinSubspace(M), n);
      f :=  cuspidal_part * eisenstein_part;

   else
      MC := AssociatedSpaceOverZ(M);
      if GCD(Level(M),n) eq 1 then
         f_cuspidal := ComputeHeckePolynomialUsingModularSymbols(
                            CuspidalSubspace(MC), n, Proof);      
         f_eisenstein := ComputeHeckePolynomialUsingEisensteinSeries(
                            EisensteinSubspace(MC), n);
         f := f_cuspidal * f_eisenstein;
      else
         f := ComputeHeckePolynomialUsingModularSymbols(MC, n, Proof);
      end if;
   end if;
   f := CoerceUnivariatePolynomialElement(PolynomialRing(BaseRing(M)), f);
   Append(~M`hecke_polynomial, <n, f>);
   return f;
end intrinsic;

function Tp(N, eps, k, p, r, f)
   assert Type(p) eq RngIntElt;
   assert Type(r) eq RngIntElt;
   assert Type(f) eq RngSerPowElt;

   if r gt 1 then
      return Tp(N,eps,k, p,1,Tp(N,eps,k,p,r-1,f)) 
                     - Evaluate(eps,p)*p^(k-1)*Tp(N,eps,k,p,r-2,f);
   end if;
   if r eq 1 then
      R := Parent(f);
      q := R.1;
      prec := Floor(((AbsolutePrecision(f)-1)/p) + 1);
      h := &+[R|Coefficient(f,n*p)*q^n + 
                     Evaluate(eps,p)*p^(k-1)*Coefficient(f,n)*q^(n*p) :
                    n in [0..prec-1]] + O(q^prec);
      return h; 
   end if;
   if r eq 0 then
      return f;
   end if;
end function;

function Tn(N, eps, k, n, f)
   // The image T_n(f) of f under the Hecke operator T_n on level-N
   // weight-k modular forms with character eps.
   assert Type(N) eq RngIntElt and N ge 1;
   assert Type(n) eq RngIntElt and n ge 1;
   assert Type(k) eq RngIntElt and k ge 2;
   assert Type(eps) eq GrpDrchElt;
   assert Type(f) eq RngSerPowElt;

   for p in Factorization(n) do
      f := Tp(N, eps, k, p[1],p[2], f);
   end for;
   return f;
end function;

/*
intrinsic UOperator(f::ModFrmElt, p::RngIntElt) -> ModFrmElt
{The modular forms f(q^p).}
   error "not written.";
end intrinsic;
*/

// write each element of B2 as a linear combination
// of the elements of B1.
function PowerSeriesCoordinates(B1, B2)
   if #B1 eq 0 then
      return [];
   end if;
   K := FieldOfFractions(BaseRing(Parent(B1[1])));
   prec := AbsolutePrecision(B1[1]);
   V := VectorSpace(K, prec);
   C1 := [V![Coefficient(f,i) : i in [0..prec-1]] : f in B1];
   C2 := [V![Coefficient(f,i) : i in [0..prec-1]] : f in B2];
   W := VectorSpaceWithBasis(C2);
   D1 := [Coordinates(W,v) : v in C1];
   return D1;
end function;

// The matrix A is written wrt the basis B1 of power series.  The function
// B2 return the matrix A written wrt to the basis B2.
function ChangeBasis(A, B1, B2)
   assert Type(A) eq AlgMatElt;
   assert Type(B1) eq SeqEnum;
   assert Type(B2) eq SeqEnum;
   assert #B1 eq #B2;
   assert #B1 eq Degree(Parent(A));

   n := #B1;
   if n eq 0 then
      return A;
   end if;
   K := FieldOfFractions(BaseRing(Parent(A)));
   C := PowerSeriesCoordinates(B1, B2);
   S := MatrixAlgebra(K, n)!(&cat[Eltseq(v) : v in C]);
   return S^(-1)*A*S;
end function;


// The matrix Wq of Atkin Lehner q on Basis(M).
// Here M must be completely new.
function WqOnNewSpace(M, q, prec) 
   assert Type(M) eq ModFrm;
   assert Type(q) eq RngIntElt;
   assert IsPrime(q);
   assert Type(prec) eq RngIntElt;
   assert M subset NewSubspace(M);

   N := Level(M);
   if Dimension(M) eq 0 or Valuation(N, q) lt 1 then
      return MatrixAlgebra(BaseRing(M),Dimension(M))!1;
   end if;
/* Don't do this -- it is very very very slow, though conceptually simple.
   elif Valuation(N, q) eq 1 then   
      // use that Wq = -Tq/q^(k/2-1) on newforms when q is prime.
      return -(1/q^((Weight(M) div 2)-1))* HeckeOperator(M,q);
*/
   D := [];
   B := [];
   for i in [1..NumberOfNewformClasses(M)] do
      f := Newform(M,i);
      eps := AtkinLehnerEigenvalue(f,q);
      D := D cat [eps : i in [1..Dimension(Parent(f))]];
      B := B cat [PowerSeriesRing(BaseRing(M))!PowerSeries(g,prec) : g in Basis(Parent(f))];
   end for;
   assert #D eq Dimension(M);
   Wq := DiagonalMatrix(MatrixAlgebra(Rationals(),Dimension(M)), D);
   // Now change from the B basis that Wq is written wrt to Basis(M),
   // which is what we want.
   return ChangeBasis(Wq, B, [PowerSeries(f,prec) : f in Basis(M)]);   
end function;


function Frob(n, f)
   assert Type(n) eq RngIntElt;
   assert Type(f) eq RngSerPowElt;
   x := Parent(f).1;
   return Evaluate(f,x^n) + O(x^AbsolutePrecision(f));
end function;


// Matrix of Wq on image of M at level Level(M)*t on the basis B, which
// is returned as the second part.
// Here q should be prime and prec is the precision of computations.
function WqOnImageOfNewSpace(M, q, t, prec)
   assert Type(M) eq ModFrm;
   assert Type(q) eq RngIntElt;
   assert IsPrime(q);
   assert Type(prec) eq RngIntElt;

   K := FieldOfFractions(BaseRing(M));
   k := Weight(M);   
   assert Type(k) eq RngIntElt;
   WqOld := WqOnNewSpace(M, q, prec);
   x := PowerSeriesRing(K).1;
   B := [];
   WqB := [];
   qq := q^Valuation(Level(M)*t,q);
   for f in Basis(M) do
      for d in Divisors(t) do
         ps := PowerSeries(f,prec);
         Append(~B, Frob(d,ps));
         alpha := Valuation(Level(M)*t, q);
         beta := alpha - Valuation(Level(M),q);
         gamma := Valuation(d, q);
         dd := Integers()!(q^(beta-2*gamma)*d);
         assert IsEven(k);
         pow := (k div 2)*(beta-2*gamma);
         g := Evaluate(PowerSeries(f*WqOld,prec),x^dd)*q^pow + O(x^prec);
         Append(~WqB, g);
      end for;
   end for;
   C := PowerSeriesCoordinates(B, WqB);
   Wq := MatrixAlgebra(K,#B)!(&cat[Eltseq(v) : v in C]);
   return <Wq, B>;
end function;


function compute_Wq_cuspidal(M, q)
   if Dimension(M) eq 0 then
      return MatrixAlgebra(RationalField(),0)!1;
   end if;
   assert Type(M) eq ModFrm;
   assert Type(q) eq RngIntElt;
   assert IsPrime(q);
   assert IsCuspidal(M);

   K := FieldOfFractions(BaseRing(M));
   if Type(BaseRing(M)) ne RngInt then
      MZ := AssociatedSpaceOverZ(M);      
      Wq := compute_Wq_cuspidal(MZ);
      return MatrixAlgebra(FieldOfFractions(K),Dimension(M))!Wq;
   end if;

   N := Level(M);
   prec := PrecisionBound(M) + 1;
   type := SpaceType(M);    assert type in {"cusp", "cusp_new"};  // the types we know about.
   param := SpaceTypeParam(M);   assert param eq 0 or IsPrime(param);  // all we know about.
   B := [];
   Wq := MatrixAlgebra(K,Dimension(M))!0;
   for d in Divisors(N) do 
      NS := false;
      if type eq "cusp" then
         NS := NewSubspace(CuspidalSubspace(ModularForms(M,d)));
      elif type eq "cusp_new" then
         if param eq 0 then
            if d eq N then
               NS := M;
            end if;
         else 
            if Valuation(d,param) eq Valuation(N,param) then
               NS := NewSubspace(CuspidalSubspace(ModularForms(M,d)));
            end if;
         end if;
      end if;
      if Type(NS) ne BoolElt then
         W, basis := Explode(WqOnImageOfNewSpace(NS, q, N div d, prec));
         for i in [1..Degree(Parent(W))] do 
            for j in [1..Degree(Parent(W))] do 
               Wq[#B + i, #B + j] := W[i,j];
            end for;
         end for;
         B := B cat basis;
      end if;
   end for;
   C := ChangeBasis(Wq, B, [PowerSeries(f,prec) : f in Basis(M)]);   
   return C;
end function;


function compute_Wq(M, q)
   assert Type(M) eq ModFrm;
   assert Type(q) eq RngIntElt;
   
   S := CuspidalSubspace(M);
   E := EisensteinSubspace(M);
   WS := &*[compute_Wq_cuspidal(S,F[1]) : F in Factorization(q)];
   WE := AtkinLehnerOnEisensteinModularFormsSpace(E,q);
   W := DirectSum(WS,WE);
   BS := Basis(S);
   BE := Basis(E);
   BS_BE := [M!b : b in BS] cat [M!b : b in BE];
   // Change of basis from BS cat BE to Basis(M).
   C := MatrixAlgebra(RationalField(),Dimension(M))!
                      &cat[Eltseq(b) : b in BS_BE];
   return C^(-1)*W*C;
end function;

intrinsic AtkinLehnerOperator(M::ModFrm) -> AlgMatElt
{}
   return AtkinLehnerOperator(M,Level(M));
end intrinsic; 

intrinsic AtkinLehnerOperator(M::ModFrm, q::RngIntElt) -> AlgMatElt
{The matrix representing the q-th Atkin-Lehner involution W_q on M
with respect to Basis(M).  At present M must be a cuspidal space of 
modular forms for Gamma_0(N).  Note that the Atkin-Lehner operator
W_q can send a q-expansion with integer coefficients to one with nontrivial
denominators.}

   if Dimension(M) eq 0 or GCD(q,Level(M)) eq 1 then
      return MatrixAlgebra(BaseRing(M),0)!1;
   end if;

   require GCD(q,Level(M) div q) eq 1 : "Argument 2 (=q) must have the property " *
       "that GCD(q, N/q) = 1, where N is the level of argument 1.";
   require IsGamma0(CuspidalSubspace(M)) 
             : "Argument 1 must be a space of modular forms on Gamma_0.";
   require Type(Weight(M)) eq RngIntElt : 
             "Argument 1 must have integral weight.";
   if Dimension(EisensteinSubspace(M)) ne 0 then
      require q eq Level(M) : "Atkin-Lehner W_q for q not equal to the "*
                              "level is not yet implemented on Eisenstein spaces.";
   end if;

   if not assigned M`atkin_operator then
      M`atkin_operator := [];
   end if;
   if exists(i) { i : i in [1..#M`atkin_operator] | 
        M`atkin_operator[i][1] eq q } then
      return M`atkin_operator[i][2];
   end if;
   Wq := compute_Wq(M,q);       
   Append(~M`atkin_operator, <q, Wq>);
   return Wq;
end intrinsic;

