\documentclass{article}
\relpenalty=10000
\binoppenalty=10000
\bibliographystyle{amsplain}
\include{macros}
\newcommand{\be}{\mathbf{e}}

\title{{\large\sc Computing Characteristic Polynomials of Hecke 
Operators When the Weight is Large}\footnote{This paper was written while 
the author was an invited professor at the University of Rennes, France and an
NSF Postdoctoral Fellow at Harvard University.}}
\author{William A. Stein\\{\sf was@math.harvard.edu}\\{\sf http://math.harvard.edu/\~{\mbox{}}was}}
\begin{document}
\maketitle
\begin{abstract}
We describe an algorithm to compute the Hecke operator $T_{p,k}$
acting on weight-$k$ modular forms on $\SL_2(\Z)$. 
This algorithm is optimized for the case when~$p$ is large.
\end{abstract}

\section{Introduction}

We describe an algorithm to compute the Hecke operator $T_{p,k}$ on
weight-$k$ modular forms on $\SL_2(\Z)$.  The algorithm can be
generalized to level~$N$, given a good method to compute a basis of
$q$-expansions to ``low precision'' for the space of modular forms of
weight~$k$ and level~$N$.  By ``low precision'' we mean to order
$O(q^{np+1})$, where $T_1,T_2,\ldots, T_n$ generate the Hecke algebra
$\T$ as a ring.  In contrast to previously published algorithms, this
algorithm is targeted at the case when~$p$ is large.  Each of the key
ideas which are used in the algorithm boil down to fairly elementary
linear algebra. 

\subsection{Basic Facts About Modular Forms}
We briefly recall everything that one needs to know about modular
forms in order to read this paper.  Fix an even integer~$k$.  Let
$M_k$ denote the space of weight-$k$ modular forms for $\SL_2(\Z)$ and
$S_k$ the subspace of cusp forms, as defined in
\cite{serre:arithmetic}, except Serre's $M_{k}$ is my
$M_{2k}$.  Thus $M_k$ is a $\C$-vector space of dimension 
$1 - \frac{k}{2} +\lfloor\frac{k}{4}\rfloor + \lfloor\frac{k}{3}\rfloor$,
and $M_k$ is equipped with a ring $\T=\Z[\ldots T_{p,k} \ldots]
\subset \End(M_k)$ of \defn{Hecke operators}.  Moreover, there is an
injective $q$-expansion map $M_k\hookrightarrow \C[[q]]$.  For
example, when $k \geq 4$ there is an {\em Eisenstein series} $E_k$, which lies in $M_k$.
The first two Eisenstein series are
$$E_4(q) = \frac{1}{240} + \sum_{n\geq 1} \sigma_4(n) q^n
\,\,\text{ and }\,\,
E_6(q) = \frac{1}{504} + \sum_{n\geq 1} \sigma_6(n) q^n,$$
where $q=e^{2\pi{}iz}$, and
the {\em $\Delta$ function} is
$$\Delta(q) = q\prod_{n=1}^{\infty}(1-q^n)^{24}.$$
For every prime number~$p$, the {\em Hecke operator} $T_{p,k}$ acts on $M_k$ by
$$T_{p,k}\left(\sum_{n\geq 0} a_n q^n\right)
   = \sum_{n\geq 0} a_{np}q^n + p^{k-1} a_n q^{np}.$$
\begin{proposition}\label{prop:basis}
The set of modular forms
$E_4^a E_6^b$ is a basis for $M_k$, 
where~$a$ and~$b$ range through nonnegative integers such
that $4a+6b=k$.
Moreover, $S_k$ is the subspace of $M_k$ 
of elements whose $q$-expansion has valuation at least~$1$.
\end{proposition}



\section{An Algorithm to Compute $T_p$ When~$p$ is Large}
\subsection{The Naive Approach}
Let~$k$ be an even positive integer and~$p$ be a prime.  Our goal is
to compute the characteristic polynomial of the Hecke operator
$T_{p,k}$ acting on $S_k$.  Note that when~$k$ and~$p$ are both
reasonably large, e.g., $k=886$ and $p=59$ (see
Section~\ref{sec:examples}), then the coefficients of the
characteristic polynomial are huge---they are on the
order of $p^{k/2-1}$.  One way to compute the characteristic polynomial
of $T_{p,k}$ is to first 
use the definition to compute the matrix $[T_{p,k}]$ of $T_{p,k}$ on the basis constructed
in Proposition~\ref{prop:basis}, where $E_4$
and $E_6$ are computed to precision $p\dim M_k$, and to then
compute the characteristic polynomial of $[T_{p,k}]$ using,
e.g., a modular algorithm.
The difficulty with this approach is that the coefficients
of the $q$-expansions involved quickly become enormous, so both storing
them and computing with them is  costly and the components 
of $[T_{p,k}]$ are also huge.

\subsection{Our Algorithm: The Eigenform Method}
We now describe another solution to the above problem.
A conjecture of Maeda asserts that $S_k$ is spanned by the
$\Gal(\Qbar/\Q)$-conjugates of a single eigenform $f=\sum b_n q^n$.
This conjecture has been checked for all $k<??$\edit{ask Kevin for the
reference!}.  For simplicity of exposition, we assume this conjecture,
though it is easy to get around this assumption. 
We will refer to this eigenform~$f$ as ``Maeda's eigenform''.

Victor Miller shows in his thesis that $S_k$ has a unique basis 
$f_1, \ldots, f_d\in \Z[[q]]$ with $a_i(f_j)=\delta_{ij}$, i.e., 
the left $d\times d$ block of coefficients is the identity matrix.
Again, the requirement that there is such a basis can be avoided,
but for simplicity of exposition we assume there is such a basis.
We refer to the basis $f_1, \ldots, f_d$ as ``Miller's basis''.

Using the definition, compute the matrix $[T_2]$ of $T_2$ with respect to 
Miller's basis.   (See Section~\ref{sec:examples} for a comparison of the
size of the entries of $[T_2]$ with respect to Miller's basis and
other bases.)  Then compute an eigenvector $\e=(e_1,\ldots,e_d)\in K^d$ 
for $[T_2]$, where~$K$ is a number field 
(see Algorithm~\ref{alg:eigenvector} below). 
The components of $T_2$ are not very large,  so the numbers 
involved in computing~$\be$ are also not so large. 
Observe that, up to Galois conjugacy, 
$$b_p = \sum_{i=1}^d  (e_i/e_1) a_p(f_i),$$
where $b_p$ are the cofficients of Maeda's eigenform~$f$ mentioned
above.    For example, since the $f_i$ are Miller's basis, if $p\leq d$
then $b_p = e_i/e_1$.
Once we have computed $b_p$, we can compute the characteristic 
polynomial of $T_p$, because it is 
the minimal polynomial of $b_p$.

Unfortunately, computing $b_p$ directly is extremely costly because
the divisions $a_p(f_j)/a_1(f_j)$ lead to massive coefficient explosion,
and the same remark applies to computing the minimal polynomial of $b_p$.
Instead we compute the reductions $\overline{b}_p$ 
modulo~$\ell$ and the  characterisic polynomial of $\overline{b}_p$ 
modulo~$\ell$ for many primes~$\ell$, then recover $b_p$ using the 
Chinese Remainder Theorem.  Deligne's bound\edit{reference, and say
what it is.} on the 
magnitude of Fourier coefficients tells us how many primes we need 
to work modulo.  

More precisely, the reduction modulo~$\ell$ steps are
proceed as follows.
The field~$K$ can be viewed as 
$\Q[x]/(f(x))$ where $f(x)\in \Z[x]$ is the characteristic
polynomial of $T_2$.  We work only modulo primes such that
\begin{enumerate}
\item $f(x)$ has no repeated roots modulo~$\ell$,
\item $\ell$ does not divide any denominator involved in $\be$, and
\item the image of $e_1$ in $\F_\ell[x]/(f(x))$ is invertible.
\end{enumerate}
For each such prime, we compute the image $\overline{b}_p$ of $b_p$ in
the reduced Artin ring $\F_\ell[x]/(f(x))$.  Then the characterisic
polynomial of $T_p$ modulo~$p$ equals the minimal polynomial of
$\overline{b}_p$.  This modular arithmetic is fast and
requires negligible storage.  Most of the time is
spent doing the Chinese Remainder Theorem computations,
which we do repeatedly each time we do a few computations of the
characteristic polynomial of $T_p$ modulo $\ell$.

\subsection{How to Write Down an Eigenvector over an Extension Field}
The following algorithm, which was suggested to the author
by H.~Lenstra, produces an eigenvector defined over an 
extension of the base field.

\begin{algorithm}\label{alg:eigenvector}
Let~$A$ be an $n\times n$ matrix over an arbitrary field~$k$ and 
suppose that the characteristic polynomial $f(x)=x^n+\cdots+a_1 x + a_0$ 
of~$A$ is irreducible.   Let~$\alpha$ be a root of $f(x)$
in an algebraic closure~$\kbar$ of~$k$. 
Factor $f(x)$ over $k(\alp)$ as
 $f(x) = (x-\alp) g(x)$. 
Then for any element $v\in k^n$ the vector
$g(A)v$ is either $0$ or it is an eigenvector of~$A$ with eigenvalue~$\alp$. 
The vector $g(A)v$ can be computed by finding 
$Av$, $A(Av)$, $A(A(Av))$, and then using that
  $$g(x)=x^{n-1}+c_{n-2} x^{n-2}+\cdots+c_1 x+ c_0,$$
where the coefficients $c_i$ are determined by the recurrence
$$c_0 = - a_0/\alp,\qquad  c_i = (c_{i-1}-a_i)/\alp.$$

We prove below that $g(A)v\neq 0$ for all vectors~$v$ not 
in a proper subspace of $k^n$.  Thus with high probability, a
``randomly chosen''~$v$ will have the property that $g(A)v\neq 0$.  
Alternatively, if $v_1,\ldots v_n$ form a basis for $k^n$, then 
$g(A)v_i$ must be nonzero for some~$i$.
\end{algorithm}
\begin{proof}
By the Cayley-Hamilton theorem \cite[XIV.3]{lang:algebra}
we have that $f(A)=0$.  Consequently, for any $v\in k^n$,
we have $(A-\alp)g(A)v=0$ so that $A g(A)v = \alp v$. 
Since $f$ is irreducible it is the polynomial of least
degree satisfied by $A$ and so $g(A)\neq 0$. 
Therefore $g(A)v\neq 0$ for all $v$ not in the proper
closed subset $\ker(g(A))$. 
\end{proof}


\section{Some Examples}\label{sec:examples}
Write this.\edit{Remarks about how big the numbers are when doing the $k=886$, 
$p=59$ computation, and how long it takes.}
Compare sizes with respect to various bases. 


\section{Magma Implementation of the Algorithm}
For the reader's convenience, we give a complete \magma{}
implementation of the algorithm described in this paper.  
The implementation below only works when $p \leq \dim S_k$.  

\begin{verbatim}
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<q> := PowerSeriesRing(RationalField());
   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;
\end{verbatim}

\bibliography{biblio}
\end{document}