//genspace.m

/*
absprec := 10000;
S<w> := PowerSeriesRing(Rationals());
prE4 := Eisenstein(4,w+O(w^absprec));
prE6 := Eisenstein(6,w+O(w^absprec));
*/

declare verbose Level1, 2;    // since mine isn't in the kernel yet.

intrinsic BasisSk(k::RngIntElt, K::Fld, prec::RngIntElt) -> SeqEnum
{}
   if IsOdd(k) or k lt 4 then
      return [];
   end if;


   vprint Level1 : "Compute the primitive Eisenstein series.";
   // BUG in MAGMA -- if you compute the Eisenstein series
   // over q then it gives a TOTALLY BOGUS ANSWER!
   R<q> := PowerSeriesRing(K);
   S<w> := PowerSeriesRing(Rationals());
   E4 := R!(Eisenstein(4,w+O(w^prec)));
   E6 := R!(Eisenstein(6,w+O(w^prec)));

   vprint Level1 : "Find all pairs of integers a, b such that 4*a + 6*b = k.";
   // a = (k-6*b)/4
   ab := [<Integers()!((k-6*b)/4),b> : b in [0..Floor(k/6)] | (k-6*b)/4 in Integers()];
 
   t := Cputime(); 
   vprint Level1: "Making power list.";
   E4pows := [1, E4]; 
   for i in [1..Max([x[1] : x in ab])] do
      Append(~E4pows, E4pows[#E4pows]*E4);
   end for;
   E6pows := [1, E6]; 
   for i in [1..Max([x[2] : x in ab])] do
      Append(~E6pows, E6pows[#E6pows]*E6);
   end for;
   vprint Level1: "made pow lists: ", Cputime(t);
   vprint Level1: "Compute generators.";

   gens := [E4pows[x[1]+1]*E6pows[x[2]+1] : x in ab]; 
   vprint Level1: "There are ", #gens, " generators.";

   // make basis *canonical*
   vprint Level1: "Coerce into vectors.";
   V := VectorSpace(K,prec);
   vecs := [V![Coefficient(g,n) : n in [0..prec-1]] : g in gens];

   vprint Level1: "Compute corresponding vector space.";
   W := sub<V|vecs>;   
   vprint Level1: "The space has dimension ", Dimension(W);

   // Find basis of power series.
   basis := [R!Eltseq(b) + O(q^prec): b in Basis(W)];

   // strip eisenstein series and return.
   return [basis[i] : i in [2..#basis]];
end intrinsic;


function BasisMatrix(B)
   K := BaseRing(Parent(B[1]));
   d := #B;
   I := MatrixAlgebra(K,d)!0;
   for r in [1..d] do 
      for c in [1..d] do
         I[r,c] := Coefficient(B[r],c);
      end for;
   end for;
   return I;
end function;

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

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


intrinsic HeckeOperator(n::RngIntElt, k::RngIntElt, f::RngSerPowElt)
      -> RngSerPowElt
{The image T_n(f) of f under the Hecke operator T_n on level-1 
weight-k modular forms.}
   require n ge 1 : "Argument 1 must be at least 1.";
   require k ge 2 : "Argument 2 must be at least 2.";

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

end intrinsic;


intrinsic HeckeOperatorMatrix(k::RngIntElt, 
                              n::RngIntElt, 
                              B::SeqEnum ) -> AlgMatElt

{Matrix representing the Hecke operator T_n (acting on the
right, as usual in MAGMA) on weight-k level-1 modular forms
with respect to the basis B.}
   d := #B;
   mat := BasisMatrix(B);
   I := mat^(-1);
   A := Parent(I)!&cat[ [Coefficient(g,i) : i in [1..d]]    
             where g is HeckeOperator(n,k,B[j]) : j in [1..d]];
   return I*A;
end intrinsic;


intrinsic HeckeOperatorOver(k::RngIntElt,
			   n::RngIntElt,
			   K::Fld) -> AlgMatElt
{}
   require IsEven(k) : "Argument 1 must be even.";
   requirege n, 2 ;
   d := DimensionCuspFormsGamma0(1,k);
   prec := d*n+1;

   vprint Level1: "Computing basis.";   
   B := BasisSk(k,K,prec);
   vprint Level1: "Computing Hecke operator.";   
   T := HeckeOperatorMatrix(k,n,B);

   return T;
end intrinsic;


/****************************************
  Yet another idea...
 *****************************************/
 

intrinsic T2VictorMiller(k::RngIntElt, K::Fld) -> AlgMatElt
{The matrix of T2 over K with respect to the Victor Miller's thesis
basis of S_k(1).}
   vprint Level1: "Victor Miller with field: ", K;
   return HeckeOperatorOver(k,2,K);
end intrinsic;


function CRT_vec(list)
   // input is a list of vectors of the same degree modulo various primes
   assert Type(list) eq List;
   assert #list gt 0;
   d := #list[1];
   assert d gt 0;
   M := [Characteristic(Parent(v[1])) : v in list];
   return [Integers(LCM(M))| CRT([Integers()!v[i] : v in list], M) : i in [1..d] ];
end function;

function CRT_mat(M)
   V := [* *];
   for m in M do
      Append(~V, Eltseq(m));
   end for;
   crt := CRT_vec(V);
   N := Characteristic(Parent(crt[1]));
   return MatrixAlgebra(IntegerRing(N),Degree(Parent(M[1])))!crt;
end function;



function LiftToZ(T)
   if Type(T) eq AlgMatElt then
      return MatrixAlgebra(Rationals(),Degree(Parent(T)))!LiftToZ(Eltseq(T));
   end if;
   assert Type(T) eq SeqEnum;

   d := #T;
   N := Characteristic(Parent(T[1]));

   T0 := [];
   for i in [1..d] do
      if (Integers()!T[i]) ge N/2 then     // or something like this.
         T0[i] := (Integers()!T[i])-N;
      else
         T0[i] := Integers()!T[i];
      end if;
   end for;

   return T0;
end function;


intrinsic T2VictorMiller(k::RngIntElt) -> AlgMatElt
{}
   return T2VictorMiller(k, 123456789,10^100000);
end intrinsic;

intrinsic T2VictorMiller(k::RngIntElt, 
                    pstart::RngIntElt, bound::RngIntElt) -> AlgMatElt
{}
   requirege pstart, 13;

   N := 1;
   p := NextPrime(pstart-1);
   t := 1;

   M := [* T2VictorMiller(k, GF(p)) *];
   T := LiftToZ(M[1]);
   p := NextPrime(p);      

   while N lt 2*bound do
      T2 := T2VictorMiller(k, GF(p));
      Append(~M,  T2);

      if t mod 10 eq 0 then
         M := [* CRT_mat(M) *];         
         L := LiftToZ(M[1]);
         if T eq LiftToZ(M[1]) then
            vprint Level1: "Breaking out after", t, "iterations.";
            break;
         end if;
         T := L;          
      end if;
      N *:= p;
      p  := NextPrime(p);      
      t +:= 1;
   end while;

   T2 := CRT_mat(M);
   return LiftToZ(T2);
end intrinsic;



function EigenvectorOfMatrixWithCharpoly(T, f)
/* Let T be an nxn matrix over K with irreducible characteristic
 polynomial f.  This function returns an eigenvector for T
 over the extension field K[x]/(f(x)). */

   // This is implemented using a quotient of a polynomial ring
   // because this works generically for any field.
   n  := Degree(f);
   K  := Parent(T[1,1]);
   if n eq 1 then
      return VectorSpace(K,n)![1];
   end if;
   R<x> := PolynomialRing(K);
   L<a> := quo<R | f>;
   if Coefficient(MinimalPolynomial(a),0) eq 0 then
      return false;
   end if;
   b    := 1/a;
   c    := [-b*Coefficient(f,0)];
   for i in [1..Degree(f)-1] do
      Append(~c, (c[i] - Coefficient(f,i))*b);
   end for;

   Ln := RSpace(L,n);
   v  := Ln!0;
   T  := RMatrixSpace(L,n,n)!T;

   repeat
      v[Random(1,n)] +:= 1;
      w  := c[1]*v;
      vv := v;
      for i in [2..#c] do 
         vv := vv*T;
         w +:= c[i]*vv;
      end for;
   until w ne 0;

   return w;   
end function;



function VMvector_helper(k, K)

   T2 := T2VictorMiller(k);
   vprint Level1 : "// FOUND T2!!";
   vprint Level1 : "T2 := ", T2;
   f2 := CharacteristicPolynomial(T2 : Al := "Modular", Proof := false);

   T2 := MatrixAlgebra(K, Degree(Parent(T2)))!T2;
   f2 := PolynomialRing(K)!f2;

   v2 := EigenvectorOfMatrixWithCharpoly(T2, f2) ;

   return v2;

end function;

intrinsic VMvector(k::RngIntElt) -> ModTupRngElt
{A vector v that is gives an eigenvector for T2 with respect
to the Victor Miller basis for S_k(1). 
 A cool observation is that this vector v has the property 
 that MinimalPolynomial(v[n]) is 
 the characteristic polynomial of T_n acting on modular forms
 of level 1 and weight k.  WARNING: this function assumes two things that
 are not justified.  First, it uses Proof = false in the computation
 of a characteristic polynomial.  Second, it computes
 the matrix of T2 with respect to the "Victor Miller basis" 
 by computing this matrix modulo lots of primes until the answer
 is the same after adding in 10 more primes.}

   v :=  VMvector_helper(k, RationalField());
   return (1/v[1])*v;

end intrinsic;


intrinsic VMvector(k::RngIntElt, K::Fld) -> SeqEnum
{}
   v := VMvector(k);
   R<x> := PolynomialRing(Rationals());
   f := Modulus(Parent(v[1]));
   RF<y> := PolynomialRing(K);
   L<a> := quo<RF|(RF!f)>;
   return [L!(RF!(R!z)) : z in Eltseq(v)];
end intrinsic;


intrinsic NewtonPolygon(k::RngIntElt, p::RngIntElt,
            acc_factor::RngIntElt) 
      -> SeqEnum
{}
   d := DimensionCuspFormsGamma0(1,k);  
   if p gt d then
      print "The algorithm for p > dimension is not implemented.";
   end if;
   K := pAdicField(p, d*acc_factor);  // Can we know this a priori.
   K := pAdicField(p);  // Can we know this a priori.
   v := VMvector(k, K);
   f := MinimalPolynomial(v[p]);
   return NewtonSlopes(f,p);
end intrinsic;

function CharpolyOfTnModp(T, n, p)
   assert Type(T) eq AlgMatElt;
   assert Type(n) eq RngIntElt and n ge 1 and n le Degree(Parent(T));
   assert Type(p) eq RngIntElt and IsPrime(p);

   d := Degree(Parent(T));
   Tmod := MatrixAlgebra(GF(p),d)!T;   
   f := CharacteristicPolynomial(Tmod);
   FAC := Factorization(f);
   if &*[F[2] : F in FAC] gt 1 then  // not square free
      return false;
   end if;
   v := EigenvectorOfMatrixWithCharpoly(Tmod, f);
   if Type(v) eq BoolElt then
      return false;
   end if;
   // this is worrisome -- if a component in the artin ring is nonzero; 
   // but there's probably a good reason why no components are zero!
   if Coefficient(MinimalPolynomial(v[1]),0) eq 0 then
      return false; 
   end if;
   v := (1/v[1])*v;  
   return MinimalPolynomial(v[n]);
end function;

function XXXCharpolyOfTnModp(T, n, p)
   assert Type(T) eq AlgMatElt;
   assert Type(n) eq RngIntElt and n ge 1 and n le Degree(Parent(T));
   assert Type(p) eq RngIntElt and IsPrime(p);

   d := Degree(Parent(T));
   Tmod := MatrixAlgebra(GF(p),d)!T;   
   f := CharacteristicPolynomial(Tmod);
   FAC := Factorization(f);
   if &*[F[2] : F in FAC] gt 1 then  // not square free
      return false;
   end if;

   P0 := PolynomialRing(GF(p));
   an := [ ];
   for F in FAC do
      // create finite field.
      g := F[1];
      K := GF(p^Degree(g));
      P<x> := PolynomialRing(K);
      R := Roots(P ! g);
      A := MatrixAlgebra(K,d)!Tmod;
      V := [Basis(Kernel(A-r[1]))[1] : r in R];
      V := [v/v[1] : v in V];  // normalize
      h := &*[x - v[n] : v in V];
      Append(~an, P0!h);
   end for;
   return &*an;
end function;


intrinsic ModularMethod(k::RngIntElt, n::RngIntElt, bound::RngIntElt) -> RngUPolElt
{}
   requirege k,2;
   d := DimensionCuspFormsGamma0(1,k);
   require n ge 1 and n le d : "Argument 2 must be between 1 and",d;

   T2 := T2VictorMiller(k);
   N := 1;
//   p := NextPrime(123456789);
   p := 17;
   M := [* *];
   s := Cputime();

   if IsVerbose("Level1") then
      b := 1;
      vN := p;
      q := NextPrime(p);
      while vN lt bound do
         vN *:= q;
         q := NextPrime(q);
         b +:= 1;
      end while;
   end if;
   vprint Level1 : "Roughly ", b, " dots (up to p =",q,")";

   while N lt bound do
      vprintf Level1: ".";
      f := CharpolyOfTnModp(T2,n, p);
      if Type(f) eq RngUPolElt then
         Append(~M,Eltseq(f));  
         N *:= p;
      end if;
      p := NextPrime(p);
   end while;
   vprint Level1: "time = ", Cputime(s);
   vprint Level1 : "Doing final CRT";
   t := Cputime();
   v := CRT_vec(M);
   vprint Level1 : "Time = ", Cputime(t);
   return PolynomialRing(Rationals())!LiftToZ(v);
end intrinsic;


intrinsic ModularMethod2(k::RngIntElt, n::RngIntElt) -> RngUPolElt
{}
   requirege k,2;
   d := DimensionCuspFormsGamma0(1,k);
   require n ge 1 and n le d : "Argument 2 must be between 1 and",d;

   T2 := T2VictorMiller(k);
   N := 1;
   p := 17;
   t := 1;
   M := [* *];
   s := Cputime();
   
   while true do
      vprint Level1: t, "p=",p;
      f := CharpolyOfTnModp(T2,n, p);
      if Type(f) eq RngUPolElt and Degree(f) eq d then
         Append(~M,Eltseq(f));  
         t +:= 1;
         if t mod 5 eq 0 then
            vprintf Level1 : "Doing a CRT: ";
            s0 := Cputime();
            v := CRT_vec(M); 
            vprint Level1 : "time = ", Cputime(s0), " (total =", Cputime(s), ")";
            if LiftToZ(v) eq LiftToZ(M[1]) then
               break;
            end if;
            M := [* v *];
         end if;
      end if;
      p := NextPrime(p);
   end while;

   vprint Level1: "Total time = ", Cputime(s);
   return PolynomialRing(Rationals())!LiftToZ(v);
end intrinsic;