function DimensionSk(k)
   if IsOdd(k) or k lt 4 then
      return 0;
   end if;
   return -k div 2 + Floor(k/4) + Floor(k/3);
end function;

function VictorMillerBasis(k, prec)
   assert Type(k) eq RngIntElt;
   assert Type(prec) eq RngIntElt;
   if IsOdd(k) or k lt 4 then
      return [];
   end if;
   S := PowerSeriesRing(RationalField());
   q := S.1;
   E4 := Eisenstein(4,q+O(q^prec));
   E6 := Eisenstein(6,q+O(q^prec));

   ab := [<Integers()!((k-6*b)/4),b> : b in [0..Floor(k/6)] 
                                     | (k-6*b)/4 in Integers()];
   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;

   gens := [E4pows[x[1]+1]*E6pows[x[2]+1] : x in ab]; 
   V := VectorSpace(RationalField(),prec);
   vecs := [V![Coefficient(g,n) : n in [0..prec-1]] : g in gens];
   W := sub<V|vecs>;    // reduced row-echelon form
   basis := [S!Eltseq(b) + O(q^prec): b in Basis(W)];
   return [basis[i] : i in [2..#basis]];
end function;

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

   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 function;

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

{Matrix representing the pth Hecke operator  
on weight-k level-1 cusp forms, where B is VictorMiller's basis.}
   require IsPrime(p) : "Argument 2 must be prime.";
   d := #B;
   A := MatrixAlgebra(Rationals(),d)!&cat[ [Coefficient(g,i) : i in [1..d]]    
             where g is Tp(p,k,B[j]) : j in [1..d]];
   return A;
end intrinsic;

function CharpolyOfQuotientModp(a, b, p)
   assert Type(a) eq RngUPolResElt;
   assert Type(b) eq RngUPolResElt;
   assert Type(p) eq RngIntElt and IsPrime(p);

   K := GF(p);
   R := PolynomialRing(K);
   S := Parent(a);
   f := R!Modulus(S);
   if (Discriminant(f) eq 0) or (Coefficient(f,0) eq 0) then
      return false;   
   end if;
   L := quo<R|f>;
   t, a := IsCoercible(R,(S!a));
   if not t then
      return false;
   end if;
   amod := L!a;
   if Coefficient(MinimalPolynomial(amod),0) eq 0 then
      return false;
   end if;
   t, b := IsCoercible(R,(S!b));
   if not t then
      return false;
   end if;
   bmod := L!b;

   return MinimalPolynomial(bmod/amod);
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;


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 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
   if Coefficient(f,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;
   v[Random(1,n)] +:= Random(1,10);
   v[Random(1,n)] +:= Random(1,10);
   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;


intrinsic CharacteristicPolynomialOfTp(k::RngIntElt, 
                                       p::RngIntElt) -> RngUPolElt
{The characteristic polynomial of the nth Hecke operator acting
on weight k modular forms for SL_2(Z).}
   requirege k,2;
   d := DimensionSk(k); 
   require p ge 1 and p le d : "Argument 2 must be between 1 and",d;
   require IsPrime(p) : "Argument 2 must be prime.";

   B  := VictorMillerBasis(k,2*d+1);
   T2 := HeckeOperatorMatrix(k,2,B);
   // Warning -- remove ``Proof = false'' to make this rigorous.
   f  := CharacteristicPolynomial(T2 : 
               Proof := false, Al := "Modular");
   eigen := EigenvectorOfMatrixWithCharpoly(T2, f);

   N := 1; ell := 17; t := 1; M := [* *];
   
   while true do
      f := CharpolyOfQuotientModp(eigen[1],eigen[p],ell);
      if Type(f) eq RngUPolElt and Degree(f) eq d then
         Append(~M,Eltseq(f));  
         t +:= 1;
         if t mod 5 eq 0 then
            s0 := Cputime();
            v := CRT_vec(M); 
            if LiftToZ(v) eq LiftToZ(M[1]) then
               break;
            end if;
            M := [* v *];
         end if;
      end if;
      ell := NextPrime(ell);
   end while;

   return PolynomialRing(Rationals())!LiftToZ(v);
end intrinsic;


intrinsic NewtonSlopes(f::RngUPolElt, p::RngIntElt) -> SeqEnum
{Returns the sequences of slopes of the Newton polygon of f with
 at the prime p.  The slopes are the valuations of the p-adic
 roots of f.}
   // I ripped this code off from PARI.
   n := Degree(f);
   if n le 0 then
      return [];
   end if;

   y := [* *];
   vval := [* *];
   for i in [0..n] do 
      Append(~vval,Type(Parent(Coefficient(f,0))) in [RngInt,FldRat] select
                       Valuation(Coefficient(f,i),p)
                   else
   		       Valuation(Coefficient(f,i)));
      Append(~y,0);
   end for;

   a := 1; ind := 2; 
   while a le n do
      if vval[a] ne Infinity() then
         break;
      end if;
      y[ind] := Infinity();   
      ind +:= 1;
      a +:= 1;
   end while;

   b := a+1;
   while b le n+1 do
      while vval[b] eq Infinity() do
         b +:= 1;
      end while;
      u1 := vval[a] - vval[b];
      u2 := b - a;
      c := b+1;
      while c le n+1 do
         if vval[c] eq Infinity() then
            c +:= 1;
            continue;
         end if;
         r1 := vval[a] - vval[c];
         r2 := c - a;
         if u1*r2 le u2*r1 then
            u1 := r1;
            u2 := r2;
             b := c;
         end if;
         c +:= 1;
      end while;

      while ind le b do 
         if u1 eq 0 then
            y[ind] := 0;
         else
            y[ind] := u1 / u2; 
         end if;
         ind +:= 1;
      end while;
      a := b;
      b := a+1;
   end while;
   Remove(~y,1);
   z := [* *];
   for i in [1..#y] do
      Append(~z,y[#y-i+1]);
   end for;
   return z;
end intrinsic;