freeze;
 
/****-*-magma-* EXPORT DATE: 2004-11-23 ************************************
                                                                            
                     MODFORM: Modular Forms in MAGMA
 
                              William A. Stein        
                         
   FILE: q-expansions.m

   08/15/03: (WAS) Added intrinsic CyclotomicEmbedding as requested by Kevin Buzzard.
   
   08/15/03: (WAS) Fixed significant bug found by Kevin Buzzard in 
             MoveIntoPowerSeriesRingOverANumberField
             that would have caused incorrect q-expansions
             for modular forms in some cases.
 
   04/05/03: (WAS) fixed bug in qExpansionOfModularFormsBasis:
             When level 1 one and base ring wasn't Z it would
             not coerce the basis correctly.

   ------------- old RCS comments ------------

   Revision 1.13  2002/09/11 17:36:40  was
   Nothing.

   Revision 1.12  2002/05/30 09:42:09  was
   Made default for PrecisionBound Exact := false, since this speeds up lots of other code.

   Revision 1.11  2002/05/22 21:34:56  was
   Added Bound parameter to CoefficientRing.

   Revision 1.10  2001/12/07 01:12:47  was
   fixed little bug in Coefficient.

   Revision 1.9  2001/12/05 22:11:51  was
   Added CoefficientRing and CoefficientField.

   Revision 1.8  2001/11/22 17:54:18  was
   Added a Coefficient intrinsic.

   Revision 1.7  2001/11/10 18:48:19  was
   Changed from SeqEnum to List in an assertion to ApplyMapsTo_qExpansion.

   Revision 1.6  2001/11/10 18:34:11  was
   Fixed serious bug in Reductions.

   Revision 1.5  2001/10/25 02:53:17  was
   ??

   Revision 1.4  2001/07/26 19:23:10  was
   Added some more verbosity.

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

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

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

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

import "level1.m":      Level1Basis;

import "eisenstein.m" : BasisOfEisensteinSeries;

import "misc.m" :       EchelonPowerSeriesSequence;
import "misc.m" :       SaturatePowerSeriesSequence;
import "misc.m" :       CoercePowerSeries;
import "misc.m" :       ToLowerCaseLetter;

import "modular_symbols.m" : MF_ModularSymbols;

import "newforms.m"   : ComputeEisensteinNewforms;

import "predicates.m":  CoercionMap;
import "predicates.m":  AssociatedSpaceOverZ;
import "predicates.m":  SpaceType;

import "qexp_mappings.m" : ModpReductionMaps;
import "qexp_mappings.m" : pAdicEmbeddingMaps;
import "qexp_mappings.m" : ComplexEmbeddingMaps;

forward BasisOfHalfIntegralWeightForms,
        ClearDenominators,
        ClearDenominatorsAndSaturate,
        CoercePowerSeriesElement,
        Compute_ith_Conjugate,
        Compute_qExpansionsOfForm,
        qExpansionOfModularFormsBasis,
        RestrictionOfScalarsToQ; 

function ExactPrecisionBound(M)
   assert Type(M) eq ModFrm;
   b := Dimension(M);
   q := PowerSeriesRing(BaseRing(M)).1;
   f := M.Dimension(M);
   while PowerSeries(f,b) eq O(q^b) do
      b +:= 1;
   end while;
   return b;
end function;

function ApproximatePrecisionBound(M)
   assert Type(M) eq ModFrm;
   N := Level(M);
   k := Weight(M);
   es := Dimension(EisensteinSubspace(M));
   function idxG0(n)
      return &*[Integers()| t[1]^t[2] + t[1]^(t[2]-1) : t in Factorization(n)];
   end function;
   function idxG1(n)
      return EulerPhi(n)*idxG0(n);
   end function;
   if IsGamma0(M) then
      return Ceiling(idxG0(N)*(k/12))+1 + es;
   else
      return Ceiling(idxG1(N)*(k/12))+1 + es;
   end if;
end function;

intrinsic PrecisionBound(M::ModFrm : Exact := false) -> RngIntElt
{Some integer b such that f + O(q^b) determines
 any modular form f in M.  If the optional paramater Exact 
 is set to true then the smallest integer b such that f+O(q^b) determines
 any modular form f in M is returned.}
   if IsRingOfAllModularForms(M) then
      return Infinity();
   end if;
   if not assigned M`precision_bound then
      if Exact then 
         M`precision_bound := ExactPrecisionBound(M);
      else
         return ApproximatePrecisionBound(M);        
      end if;
   end if;
   return M`precision_bound;
end intrinsic;

/*function ExactNewformPrecisionBound(M)
   error "Not written";
end function;

function ApproximateNewformPrecisionBound(M)
   return PrecisionBound(M : Exact := false);  // very bad!!
end function;

intrinsic NewformPrecisionBound(M::ModFrm : Exact := true) -> RngIntElt
{The smallest integer b such that f + O(q^b) determines
 any newform f in M.  If the optional paramater Exact 
is set to false than some integer b such that f+O(q^b) determines
any modular form f in M is returned, but b need not be minimal.}
   if IsRingOfAllModularForms(M) then
      return Infinity();
   end if;
   if not assigned M`newform_precision_bound then
      if Exact then
         M`newform_precision_bound := ExactNewformPrecisionBound(M);
      else
         return ApproximateNewformPrecisionBound(M);
      end if;
   end if;
   return M`newform_precision_bound;
end intrinsic;
*/

function idxG0(n)
   return 
      &*[Integers()| t[1]^t[2] + t[1]^(t[2]-1) : t in Factorization(n)];
end function;

function idxG1(n)
   return EulerPhi(n)*idxG0(n);
end function;

intrinsic CoefficientRing(f::ModFrmElt : Bound := -1) -> Rng
{A ring that contains the coefficients of f.  If IsNewform(f) is true, then
this is the ring generated by the Fourier coefficients.  The optional
paramater Bound can be used to obtain the ring generated by the coefficients
a_n with n <= Bound.}
   if IsNewform(f) then
      if Bound eq -1 then      
         if IsGamma0(Parent(f)) then
            bnd := Ceiling(Weight(f) * idxG0(Level(f)) / 12);
         else
            bnd := Ceiling(Weight(f) * idxG1(Level(f)) / 12);
         end if;
      else 
         bnd := Max(Bound,1);
      end if;
      qexp := PowerSeries(f,bnd+1);
      return Order([Coefficient(qexp,n) : n in [2..bnd]]);
   else
      return Parent(Coefficient(f,1));
   end if;   
end intrinsic;

intrinsic CoefficientField(f::ModFrmElt) -> Fld
{The field in the which the Fourier coefficients lie.}
   return FieldOfFractions(Parent(Coefficient(f,1)));
end intrinsic;

intrinsic Coefficient(f::ModFrmElt, n::RngIntElt) -> RngElt
{The nth Fourier coefficient of f.}
   requirege n, 0;
   return Coefficient(PowerSeries(f,n+1),n);
end intrinsic;

intrinsic qExpansion(f::ModFrmElt) -> RngSerPowElt
{}
   return PowerSeries(f);
end intrinsic;

intrinsic PowerSeries(f::ModFrmElt) -> RngSerPowElt
{}
   return PowerSeries(f,Precision(Parent(f)));
end intrinsic;

intrinsic qExpansion(f::ModFrmElt, prec::RngIntElt) -> RngSerPowElt
{}
   return PowerSeries(f,prec);
end intrinsic;

intrinsic PowerSeries(f::ModFrmElt, prec::RngIntElt) -> RngSerPowElt
{The q-expansion of the modular form f to absolute precision prec.}
   if not assigned f`q_expansion or 
                AbsolutePrecision(f`q_expansion) lt prec then
      f`q_expansion := Compute_qExpansionsOfForm(f,"normal", 0, prec)[1][1];
   end if;
   R := Parent(f`q_expansion);
   if assigned Parent(f)`q_name then
      name := Parent(f)`q_name;
   else
      name := "q";
   end if;
   AssignNames(~R,[name]);
   return f`q_expansion + O(R.1^prec);
end intrinsic;

function Turn_qExpansionsIntoModularForms(f, expansions)
   ans := [* *];
   for orb in expansions do
      mf_orb := [* *];
      for g in orb do
         M := BaseExtend(Parent(f),BaseRing(Parent(g)));
         Append(~mf_orb, M!g);
         mf_orb[#mf_orb]`degree := #orb;
      end for;
      Append(~ans, mf_orb);
   end for;
   return ans;
end function;

intrinsic Reductions(f::ModFrmElt, p::RngIntElt) -> List
{The mod p reductions of the modular form f.  Because of denominators,
the list of reductions can be empty.} 
   require IsPrime(p) and p gt 0 : "Argument 2 must be prime.";
   require Type(BaseRing(Parent(f))) in 
     {FldRat, RngInt, FldCyc, FldNum} : 
    "The base ring of the parent of argument 1 must be a number field or Z.";

   if not assigned f`modp_reductions then
      f`modp_reductions := [];   // sequence of pairs <p,fbar>
   end if;
   if exists(i) { i : i in [1..#f`modp_reductions]
                             | f`modp_reductions[i][1] eq p } then
      return f`modp_reductions[i][2];
   end if;
   prec := PrecisionBound(AmbientSpace(Parent(f)));
   if IsEisensteinSeries(f) then
      c0 := Coefficient(PowerSeries(f,1),0);
      if Type(Parent(c0)) eq FldRat then
         denom := Denominator(c0);
      else
         denom := LCM([Denominator(b) : b in Eltseq(c0)]);
      end if;
      if denom mod p eq 0 then
         return [* *];
      end if;
   end if;
   expansions := Compute_qExpansionsOfForm(f,"modp", p, prec);
   modp_forms := Turn_qExpansionsIntoModularForms(f, expansions);
   Append(~f`modp_reductions, <p, modp_forms>);

   return (f`modp_reductions[#f`modp_reductions])[2];   
end intrinsic;

intrinsic pAdicEmbeddings(f::ModFrmElt, p::RngIntElt) -> List
{The p-adic embeddings of f.}
   require IsPrime(p) and p gt 0 : "Argument 2 must be prime.";
   require Type(BaseRing(Parent(f))) in 
     {FldRat, RngInt, FldCyc, FldNum} : 
    "The base ring of the parent of argument 1 must be a number field or Z.";
/*
   if not assigned f`padic_embeddings then
      f`padic_embeddings := [];   // sequence of pairs <p,fbar>
   end if;
   if exists(i) { i : i in [1..#f`padic_embeddings]
                             | f`padic_embeddings[i][1] eq p } then
      return f`padic_embeddings[i][2];
   end if;
*/

   prec := PrecisionBound(Parent(f));
   expansions := Compute_qExpansionsOfForm(f,"padic", p, prec);
   padic_forms := Turn_qExpansionsIntoModularForms(f, expansions);
   return padic_forms;

/*
   Append(~f`padic_embeddings, <p, padic_forms>);

   return (f`padic_embeddings[#f`padic_embeddings])[2];   
*/
end intrinsic;

intrinsic ComplexEmbeddings(f::ModFrmElt) -> List
{The complex embeddings of f.}
   prec := PrecisionBound(Parent(f));
   expansions := Compute_qExpansionsOfForm(f,"complex", 0, prec);
   complex_forms := Turn_qExpansionsIntoModularForms(f, expansions);
   return complex_forms;
end intrinsic;

/*
   GeneralizedBernoulli: 
   k   : RngIntElt,   nonnegative integer weight
   eps : GrpDrchElt,  a Dirichlet character (not necessarily primitive) 

   Compute the generalized Bernoulli numbers B_{k,eps}, as defined 
   on page 44 of Diamond-Im: 

        sum_{a=1}^{N} eps(a) t*e^(at)/(e^(N*t)-1) 

                 = sum_{k=0}^{\infty} B_{k,eps}/{k!}*t^k. 

   where N is the modulus of eps. 
*/

function GeneralizedBernoulli(k, eps)
   assert Type(k) eq RngIntElt and k ge 0;
   assert Type(eps) eq GrpDrchElt;
   assert Evaluate(eps,-1) eq (-1)^k;
  
   N    := Modulus(eps);
   K    := BaseRing(eps);
   R<t> := LaurentSeriesRing(K);
   prec := k+5;
   F    := &+[Evaluate(eps,a)*t*Exp(a*t+O(t^prec))/(Exp(N*t+O(t^prec)) -1) 
                     : a in [1..N]];
   Bk   := Coefficient(F,k)*Factorial(k);
   return Bk;
end function;


function qExpansion_EisensteinSeries(f,prec)
   assert Type(f) eq ModFrmElt;
   assert Type(prec) eq RngIntElt;
   assert assigned f`eisenstein;
   /*

   chi is a primitive character of conductor L
   psi is a primitive character of conductor M
   MLt divides N

   E_k(chi,psi,t) is 
        c_0 + sum_{m \geq 1}[sum_{n|m}psi(n)n^{k-1}chi(m/n)] q^{mt}
   c_0=0 if L>1 
    and 
   c_0=L(1-k,psi)/2 if L=1 (that second L is an L-function L)

   */
 
   chi, psi, t := Explode(EisensteinData(f));
   L := Conductor(chi);
   M := Conductor(psi);
   mstop := (prec div t) + 1;

   m := LCM(CyclotomicOrder(BaseRing(chi)), CyclotomicOrder(BaseRing(psi)));
   R := PowerSeriesRing(m le 2 select RationalField() else CyclotomicField(m));
   q := R.1;
   k := Weight(f);
   if L eq 1 then
      c_0 := -GeneralizedBernoulli(k,psi)/(2*k);
   else
      c_0 := 0;   
   end if;

   if k eq 2 and IsTrivial(chi) and IsTrivial(psi) then
      assert t gt 1;
      E2 := PrimitiveEisensteinSeries(2,prec);
      q := Parent(E2).1;
      qexp := E2 - t*Evaluate(E2,q^t); 
   else
       qexp := c_0 + &+[ 
       &+[Evaluate(psi,n)*n^(k-1)*Evaluate(chi,m div n) : 
                 n in Divisors(m)]*q^(m*t) :
                 m in [1..mstop]];
   end if;

/*
   a := Coefficient(qexp,0);
   if Type(Parent(a)) eq FldRat then
      denom := Denominator(a);
   else
      denom := LCM([Denominator(b) : b in Eltseq(a)]);
   end if;
*/

   return qexp + O(q^prec);
end function;

function qExpansion_CreatedFrom(f,prec)
   assert Type(f) eq ModFrmElt;
   assert Type(prec) eq RngIntElt;
   assert assigned f`created_from;

   op, g, h := Explode(f`created_from);
   case op:
      when "*":
         return PowerSeries(g,prec)*PowerSeries(h,prec);
      when "+":
         return PowerSeries(g,prec)+PowerSeries(h,prec);
      when "-":
         return PowerSeries(g,prec)-PowerSeries(h,prec);
      when "scalar":
         return g*PowerSeries(h,prec);
      else:
         error "Invalid created_from attribute.";
   end case;
  
end function;

function qExpansion_Theta(f,prec)
   assert Type(f) eq ModFrmElt;
   assert Type(prec) eq RngIntElt;
   assert assigned f`theta;
 
   error "qExpansion_Theta -- not programmed";
   
end function;

function qExpansion_Gadget(f,prec)
   assert Type(f) eq ModFrmElt;
   assert Type(prec) eq RngIntElt;
   assert assigned f`gadget;

   return f`gadget(prec);
end function;

function qExpansion_EllipticCurve(f,prec)
   assert Type(f) eq ModFrmElt;
   assert Type(prec) eq RngIntElt;
   assert assigned f`elliptic_curve;

   return qEigenform(f`elliptic_curve, prec);
end function;


function Compute_ith_Conjugate(F, g, i)
   assert Type(F) in {FldPad, FldRat, FldFin};
   assert Type(g) eq RngSerPowElt;
   assert Type(i) in RngIntElt;
   assert i ge 1;

   error "Compute_ith_Conjugate -- not yet written.";
end function;

function CreateNumberFieldFromCyclotomicField(K)
   assert Type(K) eq FldCyc;
   f := DefiningPolynomial(K);
   F := NumberField(f);
   phi := hom<K -> F | F.1>;
   return F, phi;
end function;

function MoveIntoPowerSeriesRingOverANumberField(f, qexp)
   assert Type(f) eq ModFrmElt;
   assert Type(qexp) eq RngSerPowElt;
   R := BaseRing(Parent(qexp));
   if Type(R) eq FldRat then  // it already is where we want it
      if not assigned f`number_field_coercion_map then
         f`number_field_coercion_map := hom<RationalField() -> RationalField() | x :-> x>;
         f`cyclotomic_embedding_map := f`number_field_coercion_map;
      end if;        
      return qexp;
   end if;

   if not assigned f`number_field_coercion_map then  
      if Type(R) eq RngUPolRes then
         h := Modulus(R);
         K := NumberField(h);
         L := AbsoluteField(K);   
         poly_ring := Parent(h);
         f`number_field_coercion_map := hom<poly_ring -> L | L!K.1>;
         f`cyclotomic_embedding_map := hom<BaseRing(h) -> L | x :-> L!(K!x)>;
      elif Type(R) eq FldCyc then
         _,f`number_field_coercion_map := CreateNumberFieldFromCyclotomicField(R);
         f`cyclotomic_embedding_map := f`number_field_coercion_map;
      else
         error "Bug in q-expansions.m";
      end if;
   end if;
   prec := AbsolutePrecision(qexp);
   S    := PowerSeriesRing(Codomain(f`number_field_coercion_map));
   phi  := f`number_field_coercion_map;
   return &+[phi(Coefficient(qexp,n))*S.1^n : n in [0..prec-1]] + O(S.1^prec);
end function;

intrinsic CyclotomicEmbedding(f::ModFrmElt) -> Map
{The canonical map from the field Q(eps) generated by the
values of eps into the field K_f generated by the Fourier
coefficients of f.  We assume that f is a newform 
with Dirichlet character eps.}
   if not assigned f`cyclotomic_embedding_map then
      dummy := PowerSeries(f,2);
   end if;
   return f`cyclotomic_embedding_map;
end intrinsic;

function qExpansion_ModularSymbols(f,prec)
   assert Type(f) eq ModFrmElt;
   assert Type(prec) eq RngIntElt;
   assert assigned f`mf_modular_symbols;

   vprintf ModularForms, 2: 
     "Computing q-expansion using %o\n", f`mf_modular_symbols;
   eig := qEigenform(f`mf_modular_symbols,prec); 
   g := MoveIntoPowerSeriesRingOverANumberField(f, eig);
  
   R := BaseRing(Parent(g));
   if assigned f`which_conjugate and Type(R) ne FldRat then
      AssignNames(~R,[ToLowerCaseLetter(f`which_conjugate)]);
   end if;

   return g + O(Parent(g).1^prec);
end function;


function qExpansion_Element(f, prec)
   assert Type(f) eq ModFrmElt;
   assert assigned f`element;
   assert Type(prec) eq RngIntElt;
   assert assigned f`element;

   e := Eltseq(f`element);
   M := Parent(f);
   
   B := qExpansionOfModularFormsBasis(M,prec); 
   v := PowerSeriesRing(BaseRing(M))!0;
   for i in [j : j in [1..#e] | e[j] ne 0] do
      v +:= e[i]*B[i];
   end for;
   return v + O(Parent(v).1^prec);
end function;


function ApplyMapTo_qExpansion(f, psi)
   assert Type(f) eq RngSerPowElt;
   assert Type(psi) eq Map;
   R := Parent(f);
   S<q> := PowerSeriesRing(Codomain(psi));
   prec := AbsolutePrecision(f);
   K := Parent(Coefficient(f,0));
   R := Domain(psi);
   assert R eq K;
   return &+[psi(Coefficient(f,n))*S.1^n : n in [0..prec-1]] + O(S.1^prec);
end function;

function ApplyMapsTo_qExpansion(f, phi)
   assert Type(f) eq RngSerPowElt;
   assert Type(phi) eq List;
   assert #phi gt 0;
   ans := [* *];
   for i in [1..#phi] do 
      orbit_of_maps := phi[i];
      orbit := [* *];  
      for psi in orbit_of_maps do
         Append(~orbit, ApplyMapTo_qExpansion(f,psi));
      end for;
      Append(~ans, orbit);
   end for;
   return ans;
end function;


// The q-expansion of the modular form f to absolute precision prec. 
function Compute_qExpansionsOfForm(f, type, data, prec) 
   assert Type(f) eq ModFrmElt;
   assert Type(prec) eq RngIntElt;
   assert Type(data) eq Map or Type(data) eq RngIntElt;
   vprintf ModularForms, 2 : "Computing q-expansion to precision %o.\n", prec;

   case type:
      when "normal":
         if assigned f`created_from then
            vprint ModularForms, 2: "created from";
            g := qExpansion_CreatedFrom(f,prec);
         elif assigned f`eisenstein then
            vprint ModularForms, 2: "eisenstein";
            g := qExpansion_EisensteinSeries(f,prec);
         elif assigned f`theta then
            vprint ModularForms, 2: "theta";
            g := qExpansion_Theta(f,prec);
         elif assigned f`q_expansion_gadget then
            vprint ModularForms, 2: "gadget";
            g := qExpansion_Gadget(f,prec);
         elif assigned f`elliptic_curve then
            vprint ModularForms, 2: "elliptic curve";
            g := qExpansion_EllipticCurve(f,prec);
         elif assigned f`mf_modular_symbols then
            vprint ModularForms, 2: "modular symbols";
            g := qExpansion_ModularSymbols(f,prec);
         elif assigned f`element then 
            vprint ModularForms, 2: "modular element";
            g := qExpansion_Element(f, prec);
         else
            error "Could not figure out how to compute q-expansion of ", f;    
         end if;
         R := Parent(g);
         return [* [g] *];
      when "modp":
         phi := ModpReductionMaps(f, data);
      when "padic":
         phi := pAdicEmbeddingMaps(f, data);
      when "complex":
         phi := ComplexEmbeddingMaps(f);
      when "map":
         phi := data;
      else      
         assert false;
   end case;
   return ApplyMapsTo_qExpansion(PowerSeries(f,prec), phi);
end function;



function BasisOfWeightOneCuspForms(M,prec)
   assert Type(M) eq ModFrm;
   assert Weight(M) eq 1;
   assert Type(prec) eq RngIntElt and prec ge 0;
   assert IsCuspidal(M);

   error "BasisOfWeightOneCuspForms -- not programmed.";   
end function;


function BasisOfWeightOneForms(M,prec)
   assert Type(M) eq ModFrm;
   assert Weight(M) eq 1;
   assert Type(prec) eq RngIntElt and prec ge 0;

   C := CuspidalSubspace(M);
   C_basis := BasisOfWeightOneCuspForms(C,prec);

   E := EisensteinSubspace(M);
   E_basis := qExpansionOfModularFormsBasis(E,prec);

   return  EchelonPowerSeriesSequence(C_basis cat E_basis);
end function;


function CuspformBasisUsingModularSymbols(M,prec)
   assert Type(M) eq ModFrm;
   assert Weight(M) ge 2;
   assert Type(prec) eq RngIntElt and prec ge 0;
   assert IsCuspidal(M);
   assert Type(BaseRing(M)) in {FldRat, RngInt};

   modsym := MF_ModularSymbols(M,+1);

   if Type(modsym) eq SeqEnum then
      Q := &cat [qExpansionBasis(m,prec) : m in modsym];
   else
      Q := qExpansionBasis(modsym,prec);
   end if;
   if #Q eq 0 then
      return Q;
   end if;
   F := BaseRing(Parent(Q[1]));
   if Type(F) eq FldCyc then
      Q := EchelonPowerSeriesSequence(&cat[
                   RestrictionOfScalarsToQ(f) : f in Q]);   
   end if;
   return Q;
end function;


function BasisUsingBrandtModule(M,prec)
   assert Type(M) eq ModFrm;
   assert Weight(M) ge 2;
   assert Type(prec) eq RngIntElt and prec ge 0;

   error "BasisUsingBrandtModule -- Not programmed.";
end function;


function BasisGotByMultiplyingFormsOfLowerWeight(M,prec)
   assert Type(M) eq ModFrm;
   assert Weight(M) ge 2;
   assert Type(prec) eq RngIntElt and prec ge 0;

   error "BasisGotByMultiplyingFormsOfLowerWeight -- Not programmed.";
end function;


function BasisOfIntegralWeightAtLeast2Forms(M, prec)
   assert Type(M) eq ModFrm;
   assert Weight(M) ge 2;
   assert Type(prec) eq RngIntElt and prec ge 0;
   assert not IsEisenstein(M);  // or there might be an infinite recursion.

   C := CuspidalSubspace(M);
   C_basis := CuspformBasisUsingModularSymbols(C,prec);

   E := EisensteinSubspace(M);
   E_basis := qExpansionOfModularFormsBasis(E,prec);
   return  EchelonPowerSeriesSequence(C_basis cat E_basis);
end function;


function MakeRestrictionOfScalars(E)
   assert Type(E) eq SeqEnum;
   assert #E ne 0;
   prec := #E;
   q := PowerSeriesRing(Parent(E[1][1])).1;
   return [&+[E[n+1][i]*q^n : n in [0..prec-1]] + O(q^prec) : i in [1..#E[1]]];
end function;


function RestrictionOfScalarsFromCycloToQ(f)
   assert Type(f) eq RngSerPowElt;

   R := BaseRing(Parent(f));  
   assert Type(R) in {FldCyc, FldRat, RngInt};
   if Type(R) in {FldRat, RngInt} then
      return [f];
   end if;

   /* f is a q-expansion over a cyclotomic field.
      Compute the sequence of restrictions to Q. */
   prec := AbsolutePrecision(f);
   R    := PowerSeriesRing(Rationals());
   q    := R.1;
   if prec eq 0 then
      return [O(q)];
   end if;
   E := [Eltseq(Coefficient(f,n)) : n in [0..prec-1]];
   return MakeRestrictionOfScalars(E);
end function;

function FldEltseq(a, d)
   e := Eltseq(a);
   return e cat [0 : i in [#e+1..d]];
end function;

function RestrictionOfScalarsFromPolResToBase(f)
   assert Type(f) eq RngSerPowElt;
   d := Degree(Modulus(BaseRing(Parent(f))));
   prec := AbsolutePrecision(f);
   E := [FldEltseq(Coefficient(f,n),d) : n in [0..prec-1]];
   return MakeRestrictionOfScalars(E);
/*   q := PowerSeriesRing(BaseRing(BaseRing(Parent(f)))).1;
   return [&+[E[n+1][i]*q^n : n in [0..prec-1]] + O(q^prec) : i in [1..#E[1]]];
*/
end function;

function RestrictionOfScalarsFromFldNumToBase(f)
   assert Type(f) eq RngSerPowElt;
   prec := AbsolutePrecision(f);
   E := [Eltseq(Coefficient(f,n)) : n in [0..prec-1]];
   return MakeRestrictionOfScalars(E);
end function;

function RestrictionOfScalarsToQ(f)
   assert Type(f) eq RngSerPowElt;

   if Type(BaseRing(Parent(f))) eq RngUPolRes then
      F := RestrictionOfScalarsFromPolResToBase(f);
   elif Type(BaseRing(Parent(f))) eq FldNum then
      F := RestrictionOfScalarsFromFldNumToBase(f);
   else
      F := [f];
   end if;

   return &cat[RestrictionOfScalarsFromCycloToQ(g) : g in F];

end function;


function BasisOfEisensteinSpace(M, prec)
   assert Type(M) eq ModFrm;
   assert Type(prec) eq RngIntElt;
   assert SpaceType(M) in { "eis", "eis_new" };
   if SpaceType(M) eq "eis_new" then
      E := [* *];
      for f in ComputeEisensteinNewforms(M) do
         for g in f do
            Append(~E, g);
         end for;
      end for;
   else
      E := EisensteinSeries(M);  
   end if;
   basis := &cat [RestrictionOfScalarsToQ(PowerSeries(f,prec)) : f in E];
   basis := EchelonPowerSeriesSequence(basis);
   return basis;
end function;


function CoercePowerSeriesElement(R, f)
   assert Type(R) eq RngSerPow;
   assert Type(f) eq RngSerPowElt;
   prec := AbsolutePrecision(f);
   return &+[(BaseRing(R)!Coefficient(f,n))*R.1^n : n in [0..prec-1]] + O(R.1^prec);
end function;


function CoerceAndEchelonBasis(basis, F)
   assert Type(basis) eq SeqEnum;
   R := PowerSeriesRing(F);
   B := [R|CoercePowerSeriesElement(R,f) : f in basis];
   return EchelonPowerSeriesSequence(B);
end function;


function qExpansionOfModularFormsBasis(M, prec) 
   assert Type(M) eq ModFrm;
   assert Type(prec) eq RngIntElt;
   prec := Max(prec, ApproximatePrecisionBound(M));

   R := PowerSeriesRing(BaseRing(M));
   q := R.1;
   if Dimension(M) eq 0 then
      return [R|];
   end if;

   if not assigned M`q_expansion_basis or M`q_expansion_basis[1] lt prec then
      if assigned M`made_from_newform then
         f := PowerSeries(M`made_from_newform,prec);
         basis_Z := ClearDenominatorsAndSaturate(
                       EchelonPowerSeriesSequence(
                       RestrictionOfScalarsToQ(f)));
         basis := [R!f : f in basis_Z];
      elif Level(M) eq 1 and not IsEisenstein(M) then
         basis := Level1Basis(Weight(M), prec);
         if IsCuspidal(M) then
            Remove(~basis,1);
         end if;
         basis := [R!f : f in basis];
      elif Type(BaseRing(M)) ne RngInt then
         assert Dimension(M) eq Dimension(AssociatedSpaceOverZ(M));
         basis_Z := qExpansionOfModularFormsBasis(
                           AssociatedSpaceOverZ(M),prec);  // recursive
         basis := [R!f : f in basis_Z];
      else  // over Z
         basis := [R|];
         if Weight(M) le 0 then
            // nothing
         elif IsEisenstein(M) then
            basis := BasisOfEisensteinSpace(M,prec);
         elif Weight(M) eq 1 then
            if IsEven(DirichletCharacters(M)) then
               // nothing
            else
               basis := BasisOfWeightOneForms(M,prec);
            end if;
         elif not (Weight(M) in Integers()) and (2*Weight(M) in Integers()) then
            basis := BasisOfHalfIntegralWeightForms(M,prec);
         else
            basis := BasisOfIntegralWeightAtLeast2Forms(M,prec);
         end if;
         basis := ClearDenominatorsAndSaturate(basis);
      end if;
      M`q_expansion_basis := <prec, basis>;
   end if;
   basis := [R|f + O(q^prec) : f in M`q_expansion_basis[2]];
   basis := basis cat [O(q^prec) : i in [#basis+1..Dimension(M)]];
   return basis;
end function;

function ClearDenominators(f)
   assert Type(f) eq RngSerPowElt;
   if Type(BaseRing(Parent(f))) eq RngInt then
      return f;
   end if;
   assert Type(BaseRing(Parent(f))) eq FldRat;
   denom := LCM([Denominator(Coefficient(f,n)) : n in [0..AbsolutePrecision(f)-1]]);
   return denom*f;
end function;

function ClearDenominatorsAndSaturate(rational_basis)
   assert Type(rational_basis) eq SeqEnum;
   if #rational_basis eq 0 then
      return rational_basis;
   end if;
   assert Type(BaseRing(Parent(rational_basis[1]))) in {RngInt, FldRat};
   no_denominators := [ClearDenominators(f) : f in rational_basis];
   basis           := SaturatePowerSeriesSequence(no_denominators);
   return basis;
end function;
