freeze;
 
/****-*-magma-* EXPORT DATE: 2004-03-08 ************************************
                                                                            
                     MODFORM: Modular Forms in MAGMA
 
                              William A. Stein        
                         
   FILE: eisenstein.m
   
   01/04/04 (WAS):
        Fixed minor bug AtkinLehnerOnEisensteinModularFormsSpace(M, q)
        that affected
        "operators.m:function compute_Wq(M, q)"---in some cases WE was over Z,
        so DirectSum(WS,WE) would terminate with an error, since WS is over Q.

   07/06/03(WAS): Fixed bug in ExtendBaseMaximal that was introduced by changing
             how Dirichlet characters work.  I had assumed that the parent
             after base extension had properties that it didn't have.  -- WAS
             Also bug in construction of Eisenstein series -- got the wrong 
             ones in some cases!  
   
   Revision 1.6  2002/09/11 18:22:18  was
   Fixed bug in ComputeAllEisensteinSeries that manifest itself
   when some of the characters eps that define M do not satisfy
   eps(-1) eq (-1)^k.

   Revision 1.5  2002/09/11 18:16:04  was
   nothing.

   Revision 1.4  2002/09/11 18:12:43  was
   Rewrote ExtendBaseMaximal to be more efficient and fix a bug.

   Revision 1.3  2002/01/17 03:54:33  was
   added a comment

   Revision 1.2  2001/10/25 02:53:03  was
   Fixed serious bug in function ExtendBaseMaximal(eps) that occured when
   the user "stupidly" created a character like
         G<a>:=DirichletGroup(27, CyclotomicField(4));
   (Note that 4 doesn't divide EulerPhi(27)!)
   Also fixed a bug in ComputeAllEisensteinSeries that might appear when there
   is more than one character.

   Revision 1.1  2001/05/30 18:51:48  was
   Initial revision


   [Theorem 15, page IV-39 in Ogg's Modular Forms book looks very relevant.]

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

forward EisensteinBasisHelper;

import "misc.m":  EchelonPowerSeriesSequence;

import "predicates.m": SpaceType;

import "q-expansions.m": BasisOfAmbientEisensteinSpace;

function QZeta(n)
   assert Type(n) eq RngIntElt; 
   assert n ge 1;
   if n le 2 then
      return RationalField();
   else
      return CyclotomicField(n);
   end if;
end function;

function ExponentOfZModN(N)
   // This could probably be optimized.
   return Exponent(UnitGroup(Integers(N)));
end function;

/* Let eps : (Z/N)^* ----> Q(zeta_n) be a Dirichlet character.
   This function returns an equal Dirichlet character
       chi : (Z/N)^* ----> Q(zeta_m)
   where m is LCM(n, exponent of (Z/N)^*).
*/
function ExtendBaseMaximal(eps)

   assert Type(eps) eq SeqEnum;
   assert #eps gt 0;
   assert Type(eps[1]) eq GrpDrchElt;

   N := Modulus(eps[1]);
   g := ExponentOfZModN(N);
   if g eq 1 then
      g := 2;
   end if;
   _, r := DistinguishedRoot(Parent(eps[1]));

   R := QZeta(LCM(g,r));
   G := DirichletGroup(Modulus(eps[1]),R);
   ans := [G!e : e in eps];
   return ans;
end function;


function MakeEisensteinSeries(M,chi,psi,t)
   assert Type(M) eq ModFrm;
   assert Type(chi) eq GrpDrchElt;
   assert Type(psi) eq GrpDrchElt;
   assert Type(t) eq RngIntElt;
   assert Level(M) mod (Conductor(chi)*Conductor(psi)*t) eq 0;

   f := HackobjCreateRaw(ModFrmElt);
   f`parent := BaseExtend(EisensteinSubspace(M), QZeta(LCM(Order(chi),Order(psi))));
   f`eisenstein := <MinimalBaseRingCharacter(AssociatedPrimitiveCharacter(chi)), 
                    MinimalBaseRingCharacter(AssociatedPrimitiveCharacter(psi)), 
                    t, chi, psi>;
   return f;
end function;

function IsInSequence(x, seq)
   for i in [1..#seq] do
      if x eq seq[i] then
         return true;
      end if;
   end for;
   return false;
end function;

function ComputeAllEisensteinSeries(M)
   assert Type(M) eq ModFrm;
   assert SpaceType(M) in {"full", "eis"};

   ans := [* *]; 
   N   := Level(M);
   k   := Weight(M);

   eps := ExtendBaseMaximal(DirichletCharacters(M));
   K   := BaseRing(eps[1]);
   G   := Elements(Parent(eps[1]));
   // Divide up G by conductor:
   C   := [[] : i in [1..N]];

   for e in G do
      Append(~C[Conductor(e)],e);
   end for;

/* Find all pairs chi, psi such that
     chi * psibar = eps (nebentypus)
   See [Miyake] Lemma 7.1.1. 
*/
   for L in Divisors(N) do
      GL := C[L];
      for R in Divisors(N div L) do
         GR := C[R];
  
         for chi in GL do
            for psi in GR do
               if Evaluate(chi*psi,-1) eq (-1)^k then
                  if IsGamma1(M) 
                       or ((not IsGamma1(M)) and 
                  IsInSequence(chi*(psi^(-1)),eps)) then
                     for t in Divisors(N div (R*L)) do
                        if Weight(M) eq 2 and t eq 1 and 
                              IsTrivial(chi) and IsTrivial(psi) then
                           // do nothing
                        else
                           Append(~ans,MakeEisensteinSeries(M,chi,psi,t));      
                        end if;
                     end for;
                  end if;
               end if;
            end for;
         end for;

      end for;
   end for;
   return ans;
end function;


intrinsic EisensteinSeries(M::ModFrm) -> List
{List of the Eisenstein series associated to M.}
/*   
   chi is a primitive character of conductor L
   psi is a primitive character of conductor R
   RLt divides N
   chi*psi = eps, or if gamma_1, the parity condition.
  
   This function can definitely be optimized much more.
*/
   require Characteristic(BaseRing(M)) eq 0 : 
           "Argument 1 must have characteristic 0.";
   if not assigned M`eisenstein_series then
      case SpaceType(M):
         when "full", "eis": 
            M`eisenstein_series := ComputeAllEisensteinSeries(M);
         when "cusp", "cusp_new":
            M`eisenstein_series := [* *]; 
         when "eis_new":
            M`eisenstein_series := [* *];
            for class in Newforms(M) do
               for f in class do
                  Append(~M`eisenstein_series, f);
               end for;
            end for;
         else: 
            error "Bug in EisensteinSeries.";
      end case; 
   end if;

   return M`eisenstein_series;
end intrinsic;


function DimensionOfEisensteinSpace(M)
   assert Type(M) eq ModFrm;
   assert IsAmbientSpace(M);

   N := Level(M);
   if IsGamma0(M) then
      if IsOdd(Weight(M)) then
         dim := 0;
      else
         dim := &+[EulerPhi(GCD(d,N div d)) : d in Divisors(N)];
         if Weight(M) eq 2 then
            dim -:= 1;
         end if;
      end if;
   elif IsGamma1(M) then
      if N in {1, 2} then
         dim := N;
      elif N eq 4 then 
         if IsOdd(Weight(M)) then
            dim := 2;
         else
            dim := 3;
         end if;
      else
         dim := &+[EulerPhi(GCD(d,N)) : d in [1..N]] div 2;
      end if;
      if Weight(M) eq 2 then
         dim -:= 1;
      end if;
   else
      dim := #EisensteinSeries(M);
   end if;
   return dim;
end function;

function AtkinLehnerOnEisensteinSeries(E)
   assert Type(E) eq ModFrmElt;
   assert IsEisensteinSeries(E);
   k := Weight(E);
   assert IsEven(k);   // This is not strictly necessary!
                       // we do this to avoid thinking about how
                       // to deal with N^(k/2).
   M := AmbientSpace(Parent(E));
   assert IsGamma0(M) or IsGamma1(M) or 
          (&and [e^2 eq 1 : e in DirichletCharacters(M)]);

   _,_,t,chi,psi := Explode(EisensteinData(E));
   L := Conductor(chi);
   R := Conductor(psi);
   N := Level(M);   

   A := N^(k div 2) / (R*t)^k * chi(-1);
   
   t_new   := N div (R*L*t);
   chi_new := psi;
   psi_new := chi;

   ES := EisensteinSeries(M) ;
   for i in [1..#ES] do F := ES[i];
      eis := EisensteinData(F);            
      if eis[3] eq t_new and eis[4] eq chi_new and eis[5] eq psi_new then
         return A*F, A, i;
      end if;
   end for;
   print "OH NO!!  Didn't find Eisenstein series!";
   assert false;
end function;

function AtkinLehnerOnEisensteinSeriesBasis(M)
   assert Type(M) eq ModFrm;
   ES := EisensteinSeries(M);
   if #ES eq 0 then
      return MatrixAlgebra(RationalField(),0)!1;
   end if;
   k := Weight(M);
   assert IsEven(k);   // This is not strictly necessary!
                       // we do this to avoid thinking about how
                       // to deal with N^(k/2).
   assert IsGamma0(M) or IsGamma1(M) or 
          (&and [e^2 eq 1 : e in DirichletCharacters(M)]);

   chi := EisensteinData(ES[1])[4];
   K := BaseRing(Parent(chi));
   W := MatrixAlgebra(K,#ES)!0;

   for i in [1..#ES] do E := ES[i];
      _, c, j := AtkinLehnerOnEisensteinSeries(E);
      W[i,j] := c;
   end for;
   return W;
end function;

function ChangeOfBasis_EisensteinSeries_To_CanonicalBasis(M)
   assert Type(M) eq ModFrm;
   assert IsEisenstein(M);

   B  := Basis(M);   
   ES := EisensteinSeries(M);
   assert #ES eq Dimension(M);

   if #ES eq 0 then
      return MatrixAlgebra(RationalField(),0)!1;
   end if;

   n := PrecisionBound(M);
   chi := EisensteinData(ES[1])[4];
   ES_qexp := [PowerSeries(e,n) : e in ES];      
   M_qexp := [PowerSeries(f,n) : f in Basis(M)];
   K := BaseRing(Parent(chi));
   Eis := RMatrixSpace(K,#ES,n)! (&cat [[Coefficient(e,i) : i in [0..n-1]] : e in ES_qexp]);
   Std := RMatrixSpace(K,Dimension(M),n)! (&cat [[Coefficient(b,i) : i in [0..n-1]] : 
                   b in M_qexp]);
   return Solution(Std,Eis);
end function;

function AtkinLehnerOnEisensteinModularFormsSpace(M, q)
   assert Type(q) eq RngIntElt;

   assert Type(M) eq ModFrm;
   assert IsEisenstein(M);
   if Dimension(M) eq 0 then  
      // Need field of fractions, since Atkin-Lehner doesn't
      // preserve integrality, and "function compute_Wq(M, q)"
      // in ModFrm/operators.m assumes this. 
      return MatrixAlgebra(FieldOfFractions(BaseRing(M)),0)!1;
   end if;

   assert q eq Level(M);
   k := Weight(M);
   assert IsEven(k);   // This is not strictly necessary!
                       // we do this to avoid thinking about how
                       // to deal with N^(k/2).
   assert IsGamma0(M) or IsGamma1(M) or 
          (&and [e^2 eq 1 : e in DirichletCharacters(M)]);

   W := AtkinLehnerOnEisensteinSeriesBasis(M);
   C := ChangeOfBasis_EisensteinSeries_To_CanonicalBasis(M);
   W := C^(-1)*W*C;
   return MatrixAlgebra(RationalField(),Nrows(W))!Eltseq(W);      
end function;


