Ticket #4340: trac-4340.patch

sage/modular/modform/all.py

@@ -20 +20 @@
 
 from vm_basis import victor_miller_basis, delta_qexp
 
-from hecke_operator_on_qexp import (hecke_operator_on_qexp,
-                                    hecke_operator_on_basis)                                    
+from hecke_operator_on_qexp import hecke_operator_on_qexp, hecke_operator_on_basis
+
 from numerical import NumericalEigenforms as numerical_eigenforms
 
 from element import is_ModularFormElement, delta_lseries

sage/modular/modform/vm_basis.py

@@ -12 +12 @@
 
 import math
 
-from sage.matrix.all import MatrixSpace
+from sage.matrix.all import MatrixSpace, Matrix
 from sage.modular.dims import dimension_cusp_forms_gamma0
-from sage.rings.all import QQ, ZZ, Integer, binomial
+from sage.rings.all import QQ, ZZ, Integer, binomial, PowerSeriesRing, O as bigO
 from sage.structure.all import Sequence
 from sage.libs.all import ntl
 from sage.misc.all import verbose
 
 from eis_series import eisenstein_series_qexp
 
-
-
 def victor_miller_basis(k, prec=10, cusp_only=False, var='q'):
     r"""
-    Compute and return the Victor-Miller basis for
-    modular forms of weight k and level 1 to precision
-    $O(q^prec)$.  if \code{cusp\_only} is True, return
-    only a basis for the cuspidal subspace.
+    Compute and return the Victor-Miller basis for modular forms of
+    weight k and level 1 to precision $O(q^prec)$.  if
+    \code{cusp\_only} is True, return only a basis for the cuspidal
+    subspace.
 
     INPUT:
         k -- an integer
@@ -37 +35 @@
         var -- string (default: 'q'
 
     OUTPUT:
-        list -- entries a power series in var defined over Q.
+        A sequence whose entries are power series in ZZ[[var]].
 
     EXAMPLES:
         sage: victor_miller_basis(1, 6)
@@ -89 +87 @@
         q + 50220*q^3 + 87866368*q^4 + 18647219790*q^5 + O(q^6),
         q^2 + 432*q^3 + 39960*q^4 - 1418560*q^5 + O(q^6)
         ]
+
+        sage: victor_miller_basis(40,200)[1:] == victor_miller_basis(40,200,cusp_only=True)
+        True
+        sage: victor_miller_basis(200,40)[1:] == victor_miller_basis(200,40,cusp_only=True)
+        True
     """
     k = Integer(k)
+    if k%2 == 1 or k==2:
+        return Sequence([])
+    elif k < 0:
+        raise ValueError, "k must be non-negative"
+    elif k == 0:
+        return Sequence([PowerSeriesRing(ZZ,var)(1).add_bigoh(prec)], cr=True)
+    e = k.mod(12)
+    if e == 2: e += 12
+    n = (k-e) // 12
 
-    R = QQ[[var]]
-    if k == 0:
-        return Sequence([R(1, prec)], cr=True)
-    elif k%2 == 1 or k < 4:
-        return Sequence([])
+    # If prec is less than or equal to the dimension of the space of
+    # cusp forms, which is just n, then we know the answer, and we
+    # simply return it.
+    if prec <= n:
+        q = PowerSeriesRing(ZZ,var).gen(0)
+        err = bigO(q**prec)
+        ls = [0] * (n+1)
+        if not cusp_only:
+            ls[0] = 1 + err
+        for i in range(1,prec):
+            ls[i] = q**i + err
+        for i in range(prec,n+1):
+            ls[i] = err
+        return Sequence(ls, cr=True)
 
-    kk = k % 12
-    if kk == 2:
-        kk += 12
-    b = None
-    for a in range(15):
-        c = kk - 4*a
-        if c % 6 == 0:
-            b = c // 6
-            break
-    assert not (b is None), "bug in VM basis"
+    F6 = ((-504)*eisenstein_series_qexp(6,prec,var=var)).change_ring(ZZ)
+
+    if e == 0:
+        A = F6.parent()(1)
+    elif e == 4:
+        A = (240*eisenstein_series_qexp(4,prec,var=var)).change_ring(ZZ)
+    elif e == 6:
+        A = F6
+    elif e == 8:
+        A = (480*eisenstein_series_qexp(8,prec,var=var)).change_ring(ZZ)
+    elif e == 10:
+        A = (-264*eisenstein_series_qexp(10,prec,var=var)).change_ring(ZZ)
+    else: # e == 14
+        A = (-24*eisenstein_series_qexp(14,prec,var=var)).change_ring(ZZ)
+
+    if n == 0:
+        return Sequence([A], cr=True)
+
+    F6_squared = F6**2
+    D = delta_qexp(prec,var=var)
+    Fprod = F6_squared
+    Dprod = D
+    R = A.parent()
+
+    if cusp_only:
+        ls = [R(0)] + [A] * n
+        start = 1
+    else:
+        ls = [A] * (n+1)
+        start = 0
+        
+    for i in range(1,n+1):
+        ls[n-i] *= Fprod
+        ls[i] *= Dprod
+
+        Fprod *= F6_squared
+        Dprod *= D
+        Dprod = Dprod.add_bigoh(prec)
+
+    if cusp_only:
+        M = Matrix(ZZ, n, prec, [x.list() for x in ls[1:]])
+        for i in range(n):
+            for j in range(i):
+                M.add_multiple_of_row(j,i,-M[j][i+1])
+    else:
+        M = Matrix(ZZ, n+1, prec, [x.list() for x in ls])
+        for i in range(1,n+1):
+            for j in range(i):
+                M.add_multiple_of_row(j,i,-M[j][i])
+
+    return Sequence([ R(x.list()).add_bigoh(prec) for x in M.rows() ], cr=True)
     
-    F4 = 240*eisenstein_series_qexp(4, prec=prec)
-    F6 = -504*eisenstein_series_qexp(6, prec=prec)
-    if var != 'q':
-        F4 = R(F4)
-        F6 = R(F6)
-    Delta = (F4**3 - F6**2)/R(1728,prec)
-    d = dimension_cusp_forms_gamma0(1, k)
-    m = Delta / (F6*F6)
-    g = m * F6**(2*d + b) * F4**a 
-    G = []
-    for j in range(d):
-        G.append(g)
-        if j < d-1:
-            g *= m
-
-    if not cusp_only:
-        G.insert(0, R(eisenstein_series_qexp(k, prec=prec)))
-
-    M = MatrixSpace(QQ, len(G), prec)
-    # we have to slice since precision in products can increase.
-    e = [g.padded_list(prec) for g in G]
-    A = M(sum(e, []))
-    # this is still provably correct -- the guess is still proven right.
-    # it's just that naive guess based on coefficients is way too big.
-    E = A.echelon_form(height_guess=10**(k))
-    return Sequence([R(list(v), prec) for v in E.rows()], cr=True)
-    
-    
-## def delta_qexp(prec=10, var='q'):
-##     """
-##     Return the q-expansion of Delta.
-##     """
-##     F4 = 240*eisenstein_series_qexp(4, prec=prec)
-##     F6 = -504*eisenstein_series_qexp(6, prec=prec)
-##     R = QQ[[var]]
-##     if var != 'q':
-##         F4 = R(F4)
-##         F6 = R(F6)
-##     return (F4**3 - F6**2)/R(1728, prec)
-
-## def eta_qexp(prec=10, var='q'):
-##     """
-##     Return the q-expansion of the Dedekind eta functions as a power
-##     series with coefficients in ZZ.
-
-##     ALGORITHM:
-##         Compute a simple very explicit modular form whose 8th power 
-##         is Delta.   Then compute the 8th power using NTL polynomial
-##         arithmetic, which is VERY fast.   This function
-##         computes a *million* terms of Delta in under a minute.
-
-##     EXAMPLES:
-##         sage: delta_qexp(7)
-##         q - 24*q^2 + 252*q^3 - 1472*q^4 + 4830*q^5 - 6048*q^6 - 16744*q^7 + O(q^7)
-##         sage: delta_qexp(7,'z')
-##         z - 24*z^2 + 252*z^3 - 1472*z^4 + 4830*z^5 - 6048*z^6 - 16744*z^7 + O(z^7)
-##         sage: delta_qexp(-3)
-##         Traceback (most recent call last):
-##         ...
-##         ValueError: prec must be positive
-##     """
-##     if prec <= 0:
-##         raise ValueError, "prec must be positive"
-##     v = [0] * prec
-##     stop = int((-1+math.sqrt(1+8*prec))/2.0)
-##     for n in range(stop+1):
-## ##
-## {{{
-## pari('eta(q+O(q^200))')
-## ///
-## 1 - q - q^2 + q^5 + q^7 - q^12 - q^15 + q^22 + q^26 - q^35 - q^40 + q^51 + q^57 - q^70 - q^77 + q^92 + q^100 - q^117 - q^126 + q^145 + q^155 - q^176 - q^187 + O(q^200)
-## }}}
-## {{{
-## sloane_find([2,5,7,12,15,22,26])
-## ///
-## Searching Sloane's online database...
-## [[1318, 'Generalized pentagonal numbers: n(3n-1)/2, n=0, +- 1, +- 2,....', [0, 1, 2, 5, 7, 12, 15, 22, 26, 35, 40, 51, 57, 70, 77, 92, 100, 117, 126, 145, 155, 176, 187, 210, 222, 247, 260, 287, 301, 330, 345, 376, 392, 425, 442, 477, 495, 532, 551, 590, 610, 651, 672, 715, 737, 782, 805, 852, 876, 925, 950, 1001, 1027, 1080, 1107, 1162, 1190, 1247, 1276, 1335]]]
-## }}}
-
-from sage.rings.all import Integer
-
 def delta_qexp(prec=10, var='q', K=ZZ):
     """
     Return the q-expansion of Delta as a power series with