Computing Characteristic Polynomials of Hecke

Operators When the Weight is Large∗

William A. Stein was@math.harvard.edu http://math.harvard.edu/˜ was

March 30, 2001

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.

1 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}, …, 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.

1.1 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 [2], except Serre’s M_{k} is my M_{2k}. Thus M_{k} is a C-vector space of dimension 1− \frac{k}{2} +⌊\frac{k}{4} ⌋+⌊\frac{k}{3} ⌋, and M_{k} is equipped with a ring T = Z[… T_{p,k} …] ⊂ End(M_{k}) of Hecke operators. Moreover, there is an injective q-expansion map M_{k} ↪ C[[q]]. For example, when k ≥ 4 there is an Eisenstein series E_{k}, which lies in M_{k}. The first two Eisenstein series are

E_{4}(q) = \frac{1}{240} + ∑_{n≥1} σ_{4}(n)q^{n} and E_{6}(q) = \frac{1}{504} + ∑_{n≥1} σ_{6}(n)q^{n},

∗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.

where q = e^{2πiz}, and the Δ function is

Δ(q) = q ∏_{n=1}^{∞}(1 − q^{n})^{24}.

For every prime number p, the Hecke operator T_{p,k} acts on M_{k} by

T_{p,k} (∑_{n≥0} a_{n}q^{n}) = ∑_{n≥0} a_{np}q^{n} + p^{k−1}a_{n}q^{np}.

Proposition 1.1. 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.

2 An Algorithm to Compute T_{p} When p is Large

2.1 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 3), 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 1.1, 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.

2.2 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(\overline{Q}/Q)-conjugates of a single eigenform f = ∑ b_{n}q^{n}. This conjecture has been checked for all k <??^{1} . For simplicity 1 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}, …, f_{d} ∈ Z[[q]] with a_{i}(f_{j}) = δ_{ij}, i.e., the left d × 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}, …, f_{d} as “Miller’s basis”.

1ask Kevin for the reference!

Using the definition, compute the matrix [T_{2}] of T_{2} with respect to Miller’s basis. (See Section 3 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}, …, e_{d}) ∈ K^{d} for [T_{2}], where K is a number field (see Algorithm 2.1 below). The components of T_{2} are not very large, so the numbers involved in computing e are also not so large. Observe that, up to Galois conjugacy,

b_{p} = ∑_{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 ≤ 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 ℓ and the characterisic polynomial of \overline{b}_{p} modulo ℓ for many primes ℓ, then recover b_{p} using the Chinese Remainder Theorem. Deligne’s bound^{2} on the magnitude of Fourier coefficients tells us how many primes we 2 need to work modulo.

More precisely, the reduction modulo ℓ steps are proceed as follows. The field K can be viewed as Q[x]/(f(x)) where f(x) ∈ Z[x] is the characteristic polynomial of T_{2}. We work only modulo primes such that

1. f(x) has no repeated roots modulo ℓ,

2. ℓ does not divide any denominator involved in e, and

3. the image of e_{1} in F_{ℓ}[x]/(f(x)) is invertible.

For each such prime, we compute the image \overline{b}_{p} of b_{p} in the reduced Artin ring F_{ℓ}[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 ℓ.

2.3 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.

Algorithm 2.1. Let A be an n×n matrix over an arbitrary field k and suppose that the characteristic polynomial f(x) = x^{n} + ⋯ + a_{1}x + a_{0} of A is irreducible.

2reference, and say what it is.

Let α be a root of f(x) in an algebraic closure \overline{k} of k. Factor f(x) over k(α) as f(x) = (x − α)g(x). Then for any element v ∈ k^{n} the vector g(A)v is either 0 or it is an eigenvector of A with eigenvalue α. 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} + ⋯ + c_{1}x + c_{0},

where the coefficients c_{i} are determined by the recurrence

c_{0} = −a_{0}/α, c_{i} = (c_{i−1} − a_{i})/α.

We prove below that g(A)v ≠ 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 ≠ 0. Alternatively, if v_{1}, … v_{n} form a basis for k^{n}, then g(A)v_{i} must be nonzero for some i.

Proof. By the Cayley-Hamilton theorem [1, XIV.3] we have that f(A) = 0. Consequently, for any v ∈ k^{n}, we have (A − α)g(A)v = 0 so that Ag(A)v = αv. Since f is irreducible it is the polynomial of least degree satisfied by A and so g(A) ≠ 0. Therefore g(A)v ≠ 0 for all v not in the proper closed subset ker(g(A)).

3 Some Examples

Write this.^{3} Compare sizes with respect to various bases. 3

4 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 ≤ dim S_{k}.

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 [];

3Remarks about how big the numbers are when doing the k = 886, p = 59 computation, and how long it takes.

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;

References

[1] S. Lang, Algebra, third ed., Addison-Wesley Publishing Co., Reading, Mass., 1993.

[2] J-P. Serre, A Course in Arithmetic, Springer-Verlag, New York, 1973, Translated from the French, Graduate Texts in Mathematics, No. 7.