1822-Explicit Presentation for Modular Symbols of Level 65.sagews

Computer Methods for Classical Modular Forms

Sage Days 20.25 talk (March 23, 2010)

William Stein, University of Washington

(see http://8tb.us/home/boothby/cover/samples/ for more)

Project: Include code in Sage for computing pictures as above for any modular form.

This talk is mainly about approaches to solving the following:

      PROBLEM:

               Compute the Hecke module $M_k(N,\varepsilon)$ for $k\geq 2$ an
               integer, and $\varepsilon:(\ZZ/N\ZZ)^* \to \CC^*$ a Dirichlet
               character.

(We are also interested in $M_k(\Gamma_H(N))$, and many other related spaces,
but will focus this talk mostly on the above problem.)

Practical Complexity: For this talk, I mostly only care about "practical complexity", so
constants and real-world implementations matter a lot.  Thus on
practical problems, an exponential time algorithm that completes in
100 hours is far more "practical" than a polynomial time algorithm that would
take years to finish (... implementing).

 

Acknowledgement: "Practical" also means open -- you get to look at the code, can change anything, etc.  Everything I'll show you in this talk today is open source and freely available in any recent copy of Sage... thanks to the hard work of:

• Jordi Quer
• David Loeffler
• Craig Citro
• Alex Ghitza
• John Cremona
• Jon Bober
• Robert Bradshaw
• Robert Miller
• Jen Balakrishnan
• me
• ... and many, many others

Notation for Dirichlet characters:  We represent a character by giving the images of "canonical" generators for $(\ZZ/N\ZZ)^*$, which are minimal at each cyclic prime power factor.

G = DirichletGroup(4*5)
G.unit_gens()

[11, 17]

show(list(G))

\newcommand{\Bold}[1]{\mathbf{#1}}\left[\left[1, 1\right], \left[-1, 1\right], \left[1, \zeta_{4}\right], \left[-1, \zeta_{4}\right], \left[1, -1\right], \left[-1, -1\right], \left[1, -\zeta_{4}\right], \left[-1, -\zeta_{4}\right]\right]

The Dimension

1. Dimension formula for any $N,k, \varepsilon$, with $k\geq 2$.   Proofs start with Riemann-Hurwitz and get complicated.
2. General formulas in Cohen-Oesterle, but with no proofs. Jordi Quer finally gave proofs in around 2007 when he visited Seattle to work on Sage.
3. Dimension formulas incredibly useful double checks on all other algorithms...
4. Computing dimension for $k=1$ is comparably much harder.

@interact
def f(N=(13,(1..500)), k=(2..100)):
    for e in DirichletGroup(N).galois_orbits():
        eps = e[0]
        html("dim $M_{%s}(%s,%s) = %s$"%(k,N,latex(eps),
                             dimension_modular_forms(eps,k)))
    html("<br><hr>")                         
    for e in DirichletGroup(N).galois_orbits():
        eps = e[0]
        html("dim $S_{%s}(%s,%s) = %s$"%(k,N,latex(eps),
                             dimension_cusp_forms(eps,k)))

Two Subproblems: 

Modular Forms = Eisenstein Subspace $\bigoplus$ Cuspidal Subspace

Eisenstein Series

• There is an explicit easy-to-compute direct formula for a basis of $q$-expansions for ${\rm Eis}_k(N,\varepsilon)$, which one can deduce by staring at Miyake's book on Modular Forms long enough.

The formulas:

Suppose $\chi$ and $\psi$ are primitive Dirichlet characters with conductors $L$ and $R$, respectively. Let $$E_{k,\chi,\psi}(q) = c_0 + \sum_{m \geq 1} \left( \sum_{n|m} \psi(n) \cdot \chi(m/n) \cdot n^{k-1}\right) q^{m} \in \QQ(\chi, \psi)[[q]],$$ where $$c_0 = \begin{cases} 0 & {\rm if }\,\,\, L>1, \\ - \frac{B_{k,\psi}}{2k} & {\rm if }\,\,\, L=1. \end{cases}$$

Theorem (Miyake, Chapter 7):

Suppose $t$ is a positive integer and $\chi$, $\psi$ are as above and that $k$ is a positive integer such that $\chi(-1)\psi(-1) = (-1)^k$. Except when $k=2$ and $\chi=\psi=1$, the power series $E_{k,\chi,\psi}(q^t)$ defines an element of $M_k(RLt,\chi \psi)$. If $\chi=\psi=1$, $k=2$, $t>1$, and $E_2(q) = E_{k,\chi,\psi}(q)$, then $E_2(q) - t E_2(q^t)$ is a modular form in $M_2(\Gamma_0(t))$.

Moreover, the set of Eisenstein series in $M_k(N, \varepsilon)$ above with $RLt \mid N$ and $\chi \psi = \varepsilon$ form a basis for the Eisenstein subspace ${\rm Eis}_k(N,\varepsilon)$.

@interact
def f(N=(1..30), k=(12,(2..100))):
    for e in DirichletGroup(N).galois_orbits():
        eps = e[0]
        html("Eis$_{%s}(%s,%s)$:"%(k,N,latex(list(eps.values_on_gens()))))
        for E in EisensteinForms(eps, k).eisenstein_series():
            view((E.parameters(), E))

Important Trick: we can compute generalized Bernoulli numbers $B_{k,\varepsilon}$ very quickly in terms of the classical $B_c$ for $c\leq k$ (see page 656 of Cohen's Number Theory and Diophantine Equations, section 9). The $B_c$ can be computed very quickly, due to fast code in PARI and also very different fast parallel code by David Harvey. 

@interact
def f(N=(1..30), k=(12,(2..100))):
    for e in DirichletGroup(N).galois_orbits():
        eps = e[0]
        html('$B_{%s,%s} = %s$'%(k, latex(eps), latex(eps.bernoulli(k))))

Project: The code for computing $B_{k,\varepsilon}$ in Sage could probably easily be made faster with some more work.  See the comments in the source code.

Computing Cusp Forms

There are (at least) four distinct algorithms for computing classical cuspidal modular forms. Each has unique advantages over all of the others, and it is important to master all of them (no software has yet done so...).

1. Explicit Generators
2. The Hijikata (Eichler-Selberg) trace formula
3. Modular symbols
4. Rational quaternion algebras and supersingular elliptic curves

1. Cusp Forms: Explicit Generators

Unique Advantages:

1. Fast for large weight (using fast univariate polynomial multiplication, thanks to FLINT)
2. Easy to understand and implement once generators are known.

 

Prototype: We have an isomorphism of rings:

$$\bigoplus_k M_k({\rm SL}_2(\ZZ)) \cong \CC[E_4, E_6]$$

where $E_4$ and $E_6$ are the weight 4 and 6 Eisenstein series.

show(eisenstein_series_qexp(4,10))

\newcommand{\Bold}[1]{\mathbf{#1}}\frac{1}{240} + q + 9q^{2} + 28q^{3} + 73q^{4} + 126q^{5} + 252q^{6} + 344q^{7} + 585q^{8} + 757q^{9} + O(q^{10})

time E = eisenstein_series_qexp(4,10^4)

Time: CPU 0.10 s, Wall: 0.39 s

time Delta = delta_qexp(10^5)

Time: CPU 0.22 s, Wall: 0.89 s

Project: Clean up some of these $q$-expansion functions.  For example, look at this mess below, where the names and input parameters for $\Delta$ and $\eta$ couldn't be more inconsistent!!  Also, there should be an easy way to find them all, e.g., modforms.[tab].

delta_qexp(10,'q',ZZ)          # by William Stein

q - 24*q^2 + 252*q^3 - 1472*q^4 + 4830*q^5 - 6048*q^6 - 16744*q^7 + 84480*q^8 - 113643*q^9 + O(q^10)

qexp_eta(ZZ[['q']], 10)        # by David Loeffler

1 - q - q^2 + q^5 + q^7 + O(q^10)

Project: One can find ring generators for other levels.  Code something systematic for small levels in Sage.  

2. Cusp Forms: The Hijikata (Eichler-Selberg) Trace Formula

Unique Advantages:

1. Impressive "practical complexity" to compute ${\rm Tr}(T_p)$ with $p$ big or with weight big.
2. Relevant to certain algorithms for computing $p$-adic modular forms.

 

See http://wstein.org/Tables/hijikata.html for the 2-page formula itself, which is an enormous sum over class numbers and solutions to certain quadratic equations modulo various integers.  There is a non-optimized PARI and Magma implementation at that page, which I wrote in 1998 (!).

 

#g = get_remote_file('http://wstein.org/Tables/hijikata.gp')
g = '/home/wstein/sd20.2/hijikata.gp'

len(open(g).read().splitlines())     # short!

135

gp.read(g)

tr(n,k,N) = tr(T_n on S_k(Gamma_0(N))). (If n | N then n must be prime.)

gp.eval('tr(17, 24, 1)')

'254028147597540'

CuspForms(1,24).T(17).trace()

254028147597540

time CuspForms(1,12).T(100000).trace()   # this just uses explicit expansion of Delta.

-2983637890141033828147200000 Time: CPU 0.01 s, Wall: 0.01 s

time gp.eval('tr(100000,12,1)')   # in theory should be asymptotically much better than explicit q-exp

'-2983637890141033828147200000' Time: CPU 0.00 s, Wall: 15.75 s

Project: Port this to Sage and make it fast enough to be useful.  Challenge: Verify that $p | \tau(p)$ for $p=7758337633$.

3. Cusp Forms: Modular Symbols

Unique Advantages:

1. Fairly general (any $N,k,\varepsilon$ with $k\geq 2$).
2. Yields unique information about special values of $L$-functions.
3. Also, unique information about abelian varieties (and motives when $k>2$) that no other method gives.

Basic Idea: 

1. Compute the homology group $H_1(X_0(N),\QQ, \{{\rm cusps}\})$ in terms of an explicit presentation (Manin symbols) got from images of the path $\{0,\infty\}$ from $0$ to $\infty$.     
2. This amounts to sparse linear algebra over $\QQ$.
3. Then compute Hecke operators using an explicit formula (and Heilbronn matrices -- see Merel's SLNM 1585). 
4. Use dense linear algebra to decompose new subspace of $S_k$ as a direct sum of simple modules.
5. Compute system of Hecke eigenvalues associated to each summand.
6. Everything generalizes to higher weight, nontrivial characters, etc.

@interact
def f(N=(30,(1..400))):
    P = P1List(N)
    def pnt(n,m):
        if m == 0: return infinity
        return float(n)/m
    def geodesic(alpha, beta):
        if alpha == infinity or beta == infinity: return line([(0,0),(0,.5)])
        return disk(((alpha+beta)/2,0), abs((alpha-beta)/2), 
                    (0,pi), thickness=1, fill=False)         
    G = Graphics()
    for i in range(len(P)):
        a,b,c,d = P.lift_to_sl2z(i)
        G += geodesic(pnt(b,d), pnt(a,c))
    G.ymin(0); G.set_aspect_ratio(1)
    html("The %s generating Manin symbols $\\gamma\\{0,\\infty\\}$ of level %s"%(len(P), N))
    show(G, figsize=10)

Explicit Presentation for Modular Symbols of Level 65

M = ModularSymbols(65,sign=1); M

Modular Symbols space of dimension 8 for Gamma_0(65) of weight 2 with sign 1 over Rational Field

M.manin_generators()

[(0,1), (1,0), (1,1), (1,2), (1,3), (1,4), (1,5), (1,6), (1,7), (1,8), (1,9), (1,10), (1,11), (1,12), (1,13), (1,14), (1,15), (1,16), (1,17), (1,18), (1,19), (1,20), (1,21), (1,22), (1,23), (1,24), (1,25), (1,26), (1,27), (1,28), (1,29), (1,30), (1,31), (1,32), (1,33), (1,34), (1,35), (1,36), (1,37), (1,38), (1,39), (1,40), (1,41), (1,42), (1,43), (1,44), (1,45), (1,46), (1,47), (1,48), (1,49), (1,50), (1,51), (1,52), (1,53), (1,54), (1,55), (1,56), (1,57), (1,58), (1,59), (1,60), (1,61), (1,62), (1,63), (1,64), (5,1), (5,2), (5,3), (5,4), (5,6), (5,7), (5,8), (5,9), (5,11), (5,12), (5,13), (5,18), (5,23), (13,1), (13,2), (13,3), (13,4), (13,5)]

M.manin_basis()

[1, 67, 68, 73, 75, 79, 80, 83]

M.manin_gens_to_basis()[:5]

[-1 0 0 0 0 0 0 0] [ 1 0 0 0 0 0 0 0] [ 0 0 0 0 0 0 0 0] [ 0 1 -2 1 -1 0 1 -1] [ 0 0 -1 1 -1 0 1 -1]

Hecke Operators at Level 65

M = ModularSymbols(65,sign=1); M

Modular Symbols space of dimension 8 for Gamma_0(65) of weight 2 with sign 1 over Rational Field

T2 = M.hecke_matrix(2); T2       # wrt rational homology

[ 3 -1 2 -1 1 0 -1 1] [ 0 0 1 1 1 0 0 0] [ 0 1/2 -1 1 1/2 1/2 3/2 -2] [ 0 3/2 0 0 1/2 1/2 1/2 -1] [ 0 3/2 2 -1 3/2 1/2 -3/2 1] [ 0 -1/2 0 0 1/2 1/2 5/2 0] [ 0 0 1 0 0 2 0 1] [ 0 0 0 -1 0 1 0 2]

show(factor(T2.charpoly()))

\newcommand{\Bold}[1]{\mathbf{#1}}(x + 1) \cdot (x - 3)^{3} \cdot (x^{2} - 3) \cdot (x^{2} + 2 x - 1)

M.integral_hecke_matrix(2)      # integral homology

[ 3 -2 2 -1 2 1 0 1] [ 0 1 2 0 1 1 0 1] [ 0 1 -1 1 0 0 1 -2] [ 0 3 0 0 -1 -1 -1 -1] [ 0 3 2 -1 0 -1 -3 1] [ 0 -1 0 0 1 1 3 0] [ 0 0 1 0 0 2 0 1] [ 0 0 0 -1 0 1 0 2]

Decomposition of New Subspace at Level 65

M = ModularSymbols(65,sign=1); M

Modular Symbols space of dimension 8 for Gamma_0(65) of weight 2 with sign 1 over Rational Field

N = M.cuspidal_subspace().new_subspace(); N

Modular Symbols subspace of dimension 5 of Modular Symbols space of dimension 8 for Gamma_0(65) of weight 2 with sign 1 over Rational Field

N.decomposition()

[ Modular Symbols subspace of dimension 1 of Modular Symbols space of dimension 8 for Gamma_0(65) of weight 2 with sign 1 over Rational Field, Modular Symbols subspace of dimension 2 of Modular Symbols space of dimension 8 for Gamma_0(65) of weight 2 with sign 1 over Rational Field, Modular Symbols subspace of dimension 2 of Modular Symbols space of dimension 8 for Gamma_0(65) of weight 2 with sign 1 over Rational Field ]

Systems of Hecke Eigenvalues

M = ModularSymbols(65,sign=1)
D = M.cuspidal_subspace().new_subspace().decomposition()
f = D[1]; f

Modular Symbols subspace of dimension 2 of Modular Symbols space of dimension 8 for Gamma_0(65) of weight 2 with sign 1 over Rational Field

There is a clever "compressed representation" of the Fourier coefficients $a_n$ that arises naturally out of the modular symbols algorithm.

A, v = f.compact_system_of_eigenvalues(prime_range(30))
print v
A

(1, 1/2*alpha + 1) [-2 2] [-1 2] [ 1 0] [ 4 -4] [ 3 -2] [-1 0] [ 0 -4] [ 1 2] [ 1 -2] [-4 8]

The product gives the actual eigenvalues:

A*v

(alpha, alpha + 1, 1, -2*alpha, -alpha + 1, -1, -2*alpha - 4, alpha + 3, -alpha - 1, 4*alpha + 4, 3*alpha + 9, 6*alpha + 6, -2*alpha - 8, 5*alpha + 1, 2*alpha, -6*alpha - 12, 3*alpha + 9, -8, -2, -7*alpha - 5, -6*alpha - 6, 6*alpha + 6, -2*alpha - 8, 6, 4*alpha + 2)

Here is a more impressive example.  The rows of $A$ are small, but the actual Hecke eigenvalues expressed in terms of a power basis are huge (see below).

f = ModularSymbols(389,sign=1)[5]
A, v = f.compact_system_of_eigenvalues(prime_range(30))
A

[ -73/8 -3/4 -61/4 73/4 -9/4 -33/4 29/4 -21/4 49/4 7/2 11/4 -1/4 -1 0 33/2 -17 -1/2 -16 16 -1/4] [ 4 0 4 -4 0 0 0 2 -2 0 0 0 0 0 -4 4 0 4 -4 0] [ 6 0 6 -8 0 4 -4 2 -6 -2 0 0 0 0 -8 8 0 8 -8 0] [ 8 0 10 -12 4 8 -8 2 -10 -4 -4 0 0 0 -12 12 0 10 -10 0] [ 12 0 20 -24 4 8 -8 8 -16 -4 -4 4 0 0 -20 20 4 20 -20 0] [ 14 4 20 -28 8 18 -16 2 -18 -8 -8 0 4 -4 -20 20 0 20 -20 2] [ 169/4 -1/2 105/2 -125/2 -3/2 29/2 -13/2 57/2 -77/2 -3 -11/2 5/2 6 4 -57 58 1 52 -60 9/2] [ -17/4 5/2 -21/2 17/2 -1/2 -5/2 5/2 -9/2 13/2 3 3/2 -1/2 -2 -4 13 -6 -1 -12 12 -1/2] [-195/4 -13/2 -117/2 159/2 -23/2 -77/2 55/2 -27/2 99/2 15 29/2 -3/2 -8 10 55 -66 -3 -52 64 -11/2] [ -37/2 1 -25 33 -5 -5 5 -17 17 2 3 -5 0 0 26 -28 -2 -26 26 -1]

A[1]

(4, 0, 4, -4, 0, 0, 0, 2, -2, 0, 0, 0, 0, 0, -4, 4, 0, 4, -4, 0)

versus

(A*v)[1]

-20146763/1097385680*alpha^19 + 20466323/219477136*alpha^18 + 119884773/274346420*alpha^17 - 753611053/274346420*alpha^16 - 381358355/109738568*alpha^15 + 3611475535/109738568*alpha^14 + 6349339639/1097385680*alpha^13 - 56878934241/274346420*alpha^12 + 71555185319/1097385680*alpha^11 + 163330998525/219477136*alpha^10 - 223188336749/548692840*alpha^9 - 169878973265/109738568*alpha^8 + 265944624817/274346420*alpha^7 + 199655892261/109738568*alpha^6 - 1167579836501/1097385680*alpha^5 - 619178000979/548692840*alpha^4 + 261766056911/548692840*alpha^3 + 4410485304/13717321*alpha^2 - 14646077211/274346420*alpha - 1604641167/68586605

(A*v)[2]

252247073/1097385680*alpha^19 - 89195955/109738568*alpha^18 - 6876716517/1097385680*alpha^17 + 13248231811/548692840*alpha^16 + 3639338697/54869284*alpha^15 - 32041245347/109738568*alpha^14 - 367946611589/1097385680*alpha^13 + 2037515640679/1097385680*alpha^12 + 816602908511/1097385680*alpha^11 - 368120881159/54869284*alpha^10 - 19595857657/1097385680*alpha^9 + 763403091515/54869284*alpha^8 - 403904463001/137173210*alpha^7 - 1743458092745/109738568*alpha^6 + 5304122556631/1097385680*alpha^5 + 9907751136883/1097385680*alpha^4 - 704659646193/274346420*alpha^3 - 236814032039/109738568*alpha^2 + 88326575941/274346420*alpha + 37821733103/274346420

Project:  I optimized the hell out of the compact_system_of_eigenvalues command in several special cases, e.g., trivial character.  But the general case is still slow.  Fix that.

4. Cusp Forms: Rational Quaternion Algebras and Supersingular Elliptic Curves

 

The algorithm is the same as what John Voight talked about on Monday, but specialized to the case of $\QQ$.  In short, suppose $p^{2r+1}$ exactly divides the level $N$.  Compute an Eichler order $R$ of level $N$ in the quaternion algebra ramified at $p,\infty$, enumerate the right ideal classes in $R$, and compute Hecke operators acting on the free abelian group on these ideal classes. 

 

Unique Advantages:

1. As a biproduct, provides explicit theta series basis for (a certain subspace of) $S_2(\Gamma_0(N))$, for certain $N$, and sometimes theta series can be efficiently computed to high precision.
2. Provides unique input (the Monodromy Pairing) needed by my algorithm for computing Tamagawa numbers of modular abelian varieties, and also for computing Kolyvagin cohomology classes. 
3. There is a canonical basis and the matrices of Hecke operators are extremely sparse, e.g., $T_{p}$ has (at most) $p+1$ nonzero entries per row, no matter what the level!!!
4. Mestre Method of Graphs variant is very fast and easy to implement (when applicable).

Example: Level $65 = 5\cdot 13$

B = BrandtModule(5, 13); B

Brandt module of dimension 6 of level 5*13 of weight 2 over Rational Field

B.quaternion_algebra()

Quaternion Algebra (-2, -5) with base ring Rational Field

B.order_of_level_N()

Order of Quaternion Algebra (-2, -5) with base ring Rational Field with basis (1/2 + 1/2*j + 7/2*k, 1/4*i + 1/2*j + 41/4*k, j + 7*k, 13*k)

B.right_ideals()

(Fractional ideal (2 + 2*j + 14*k, i + 2*j + 41*k, 4*j + 28*k, 52*k), Fractional ideal (2 + 2*j + 14*k, 2*i + 4*j + 30*k, 8*j + 4*k, 52*k), Fractional ideal (2 + 6*j + 42*k, i + 2*j + 41*k, 8*j + 56*k, 104*k), Fractional ideal (2 + 6*j + 94*k, i + 6*j + 17*k, 8*j + 56*k, 104*k), Fractional ideal (2 + 14*j + 46*k, i + 14*j + 73*k, 16*j + 112*k, 208*k), Fractional ideal (2 + 14*j + 254*k, i + 30*j + 185*k, 32*j + 224*k, 416*k))

T2 = B.hecke_matrix(2); T2

[0 1 1 1 0 0] [3 0 0 0 0 0] [1 0 0 1 1 0] [1 0 1 0 1 0] [0 0 1 1 0 1] [0 0 0 0 3 0]

factor(T2.charpoly())

(x - 3) * (x + 1) * (x^2 - 3) * (x^2 + 2*x - 1)

B.brandt_series(2)

[1/2 + q + O(q^2) 1/2 + O(q^2) 1/2 + O(q^2) 1/2 + O(q^2) 1/2 + O(q^2) 1/2 + O(q^2)] [ 1/6 + O(q^2) 1/6 + q + O(q^2) 1/6 + O(q^2) 1/6 + O(q^2) 1/6 + O(q^2) 1/6 + O(q^2)] [ 1/2 + O(q^2) 1/2 + O(q^2) 1/2 + q + O(q^2) 1/2 + O(q^2) 1/2 + O(q^2) 1/2 + O(q^2)] [ 1/2 + O(q^2) 1/2 + O(q^2) 1/2 + O(q^2) 1/2 + q + O(q^2) 1/2 + O(q^2) 1/2 + O(q^2)] [ 1/2 + O(q^2) 1/2 + O(q^2) 1/2 + O(q^2) 1/2 + O(q^2) 1/2 + q + O(q^2) 1/2 + O(q^2)] [ 1/6 + O(q^2) 1/6 + O(q^2) 1/6 + O(q^2) 1/6 + O(q^2) 1/6 + O(q^2) 1/6 + q + O(q^2)]

The Method of Graphs: For prime level

Compute Hecke action directly on the free abelian group on supersingular $j$-invariants in characteristic $p$.

@interact
def f(p=(131,tuple(prime_range(11,400)))):
    X = SupersingularModule(p)
    T2 = X.hecke_matrix(2)
    G = DiGraph(T2)
    G.plot(graph_border=True, vertex_labels=False, vertex_size=30).show(figsize=6)

Project: Implement computation of Eichler orders when $p^{2r+1}$ divides the level, with $r\geq 1$.

Project: Higher weight?

Closing Remarks

1. I have only scratched the surface.  There are many, many other relevant algorithms, e.g., SLNM 1585 gives an algorithm for computing weight 1 forms in some cases.  There are also algorithms for half-integral weight forms.
2. Edixhoven et al. have a theoretical algorithm for computing $a_p$ at level 1 in time polynomial in $\log(p)$ via a Schoff-like approach.  They posted a book about this to the arxiv, and seek feedback.
3. I wrote a book on computing modular forms that the AMS published.  It is now totally free at my website.