Ticket #4340: trac-4340.patch

File trac-4340.patch, 8.4 KB (added by craigcitro, 14 years ago)
  • sage/modular/modform/all.py

    # HG changeset patch
    # User Craig Citro <craigcitro@gmail.com>
    # Date 1225402680 25200
    # Node ID bf5e0a63f8fa631da5a13a184c926abc09792614
    # Parent  67e0603ac89003d15b4c69257b0daad0f44db95f
    Fix trac #4340: massively speed up Victor Miller basis code.
    
    diff -r 67e0603ac890 -r bf5e0a63f8fa sage/modular/modform/all.py
    a b  
    2020
    2121from vm_basis import victor_miller_basis, delta_qexp
    2222
    23 from hecke_operator_on_qexp import (hecke_operator_on_qexp,
    24                                     hecke_operator_on_basis)                                   
     23from hecke_operator_on_qexp import hecke_operator_on_qexp, hecke_operator_on_basis
     24
    2525from numerical import NumericalEigenforms as numerical_eigenforms
    2626
    2727from element import is_ModularFormElement, delta_lseries
  • sage/modular/modform/vm_basis.py

    diff -r 67e0603ac890 -r bf5e0a63f8fa sage/modular/modform/vm_basis.py
    a b  
    1212
    1313import math
    1414
    15 from sage.matrix.all import MatrixSpace
     15from sage.matrix.all import MatrixSpace, Matrix
    1616from sage.modular.dims import dimension_cusp_forms_gamma0
    17 from sage.rings.all import QQ, ZZ, Integer, binomial
     17from sage.rings.all import QQ, ZZ, Integer, binomial, PowerSeriesRing, O as bigO
    1818from sage.structure.all import Sequence
    1919from sage.libs.all import ntl
    2020from sage.misc.all import verbose
    2121
    2222from eis_series import eisenstein_series_qexp
    2323
    24 
    25 
    2624def victor_miller_basis(k, prec=10, cusp_only=False, var='q'):
    2725    r"""
    28     Compute and return the Victor-Miller basis for
    29     modular forms of weight k and level 1 to precision
    30     $O(q^prec)$.  if \code{cusp\_only} is True, return
    31     only a basis for the cuspidal subspace.
     26    Compute and return the Victor-Miller basis for modular forms of
     27    weight k and level 1 to precision $O(q^prec)$.  if
     28    \code{cusp\_only} is True, return only a basis for the cuspidal
     29    subspace.
    3230
    3331    INPUT:
    3432        k -- an integer
     
    3735        var -- string (default: 'q'
    3836
    3937    OUTPUT:
    40         list -- entries a power series in var defined over Q.
     38        A sequence whose entries are power series in ZZ[[var]].
    4139
    4240    EXAMPLES:
    4341        sage: victor_miller_basis(1, 6)
     
    8987        q + 50220*q^3 + 87866368*q^4 + 18647219790*q^5 + O(q^6),
    9088        q^2 + 432*q^3 + 39960*q^4 - 1418560*q^5 + O(q^6)
    9189        ]
     90
     91        sage: victor_miller_basis(40,200)[1:] == victor_miller_basis(40,200,cusp_only=True)
     92        True
     93        sage: victor_miller_basis(200,40)[1:] == victor_miller_basis(200,40,cusp_only=True)
     94        True
    9295    """
    9396    k = Integer(k)
     97    if k%2 == 1 or k==2:
     98        return Sequence([])
     99    elif k < 0:
     100        raise ValueError, "k must be non-negative"
     101    elif k == 0:
     102        return Sequence([PowerSeriesRing(ZZ,var)(1).add_bigoh(prec)], cr=True)
     103    e = k.mod(12)
     104    if e == 2: e += 12
     105    n = (k-e) // 12
    94106
    95     R = QQ[[var]]
    96     if k == 0:
    97         return Sequence([R(1, prec)], cr=True)
    98     elif k%2 == 1 or k < 4:
    99         return Sequence([])
     107    # If prec is less than or equal to the dimension of the space of
     108    # cusp forms, which is just n, then we know the answer, and we
     109    # simply return it.
     110    if prec <= n:
     111        q = PowerSeriesRing(ZZ,var).gen(0)
     112        err = bigO(q**prec)
     113        ls = [0] * (n+1)
     114        if not cusp_only:
     115            ls[0] = 1 + err
     116        for i in range(1,prec):
     117            ls[i] = q**i + err
     118        for i in range(prec,n+1):
     119            ls[i] = err
     120        return Sequence(ls, cr=True)
    100121
    101     kk = k % 12
    102     if kk == 2:
    103         kk += 12
    104     b = None
    105     for a in range(15):
    106         c = kk - 4*a
    107         if c % 6 == 0:
    108             b = c // 6
    109             break
    110     assert not (b is None), "bug in VM basis"
     122    F6 = ((-504)*eisenstein_series_qexp(6,prec,var=var)).change_ring(ZZ)
     123
     124    if e == 0:
     125        A = F6.parent()(1)
     126    elif e == 4:
     127        A = (240*eisenstein_series_qexp(4,prec,var=var)).change_ring(ZZ)
     128    elif e == 6:
     129        A = F6
     130    elif e == 8:
     131        A = (480*eisenstein_series_qexp(8,prec,var=var)).change_ring(ZZ)
     132    elif e == 10:
     133        A = (-264*eisenstein_series_qexp(10,prec,var=var)).change_ring(ZZ)
     134    else: # e == 14
     135        A = (-24*eisenstein_series_qexp(14,prec,var=var)).change_ring(ZZ)
     136
     137    if n == 0:
     138        return Sequence([A], cr=True)
     139
     140    F6_squared = F6**2
     141    D = delta_qexp(prec,var=var)
     142    Fprod = F6_squared
     143    Dprod = D
     144    R = A.parent()
     145
     146    if cusp_only:
     147        ls = [R(0)] + [A] * n
     148        start = 1
     149    else:
     150        ls = [A] * (n+1)
     151        start = 0
     152       
     153    for i in range(1,n+1):
     154        ls[n-i] *= Fprod
     155        ls[i] *= Dprod
     156
     157        Fprod *= F6_squared
     158        Dprod *= D
     159        Dprod = Dprod.add_bigoh(prec)
     160
     161    if cusp_only:
     162        M = Matrix(ZZ, n, prec, [x.list() for x in ls[1:]])
     163        for i in range(n):
     164            for j in range(i):
     165                M.add_multiple_of_row(j,i,-M[j][i+1])
     166    else:
     167        M = Matrix(ZZ, n+1, prec, [x.list() for x in ls])
     168        for i in range(1,n+1):
     169            for j in range(i):
     170                M.add_multiple_of_row(j,i,-M[j][i])
     171
     172    return Sequence([ R(x.list()).add_bigoh(prec) for x in M.rows() ], cr=True)
    111173   
    112     F4 = 240*eisenstein_series_qexp(4, prec=prec)
    113     F6 = -504*eisenstein_series_qexp(6, prec=prec)
    114     if var != 'q':
    115         F4 = R(F4)
    116         F6 = R(F6)
    117     Delta = (F4**3 - F6**2)/R(1728,prec)
    118     d = dimension_cusp_forms_gamma0(1, k)
    119     m = Delta / (F6*F6)
    120     g = m * F6**(2*d + b) * F4**a
    121     G = []
    122     for j in range(d):
    123         G.append(g)
    124         if j < d-1:
    125             g *= m
    126 
    127     if not cusp_only:
    128         G.insert(0, R(eisenstein_series_qexp(k, prec=prec)))
    129 
    130     M = MatrixSpace(QQ, len(G), prec)
    131     # we have to slice since precision in products can increase.
    132     e = [g.padded_list(prec) for g in G]
    133     A = M(sum(e, []))
    134     # this is still provably correct -- the guess is still proven right.
    135     # it's just that naive guess based on coefficients is way too big.
    136     E = A.echelon_form(height_guess=10**(k))
    137     return Sequence([R(list(v), prec) for v in E.rows()], cr=True)
    138    
    139    
    140 ## def delta_qexp(prec=10, var='q'):
    141 ##     """
    142 ##     Return the q-expansion of Delta.
    143 ##     """
    144 ##     F4 = 240*eisenstein_series_qexp(4, prec=prec)
    145 ##     F6 = -504*eisenstein_series_qexp(6, prec=prec)
    146 ##     R = QQ[[var]]
    147 ##     if var != 'q':
    148 ##         F4 = R(F4)
    149 ##         F6 = R(F6)
    150 ##     return (F4**3 - F6**2)/R(1728, prec)
    151 
    152 ## def eta_qexp(prec=10, var='q'):
    153 ##     """
    154 ##     Return the q-expansion of the Dedekind eta functions as a power
    155 ##     series with coefficients in ZZ.
    156 
    157 ##     ALGORITHM:
    158 ##         Compute a simple very explicit modular form whose 8th power
    159 ##         is Delta.   Then compute the 8th power using NTL polynomial
    160 ##         arithmetic, which is VERY fast.   This function
    161 ##         computes a *million* terms of Delta in under a minute.
    162 
    163 ##     EXAMPLES:
    164 ##         sage: delta_qexp(7)
    165 ##         q - 24*q^2 + 252*q^3 - 1472*q^4 + 4830*q^5 - 6048*q^6 - 16744*q^7 + O(q^7)
    166 ##         sage: delta_qexp(7,'z')
    167 ##         z - 24*z^2 + 252*z^3 - 1472*z^4 + 4830*z^5 - 6048*z^6 - 16744*z^7 + O(z^7)
    168 ##         sage: delta_qexp(-3)
    169 ##         Traceback (most recent call last):
    170 ##         ...
    171 ##         ValueError: prec must be positive
    172 ##     """
    173 ##     if prec <= 0:
    174 ##         raise ValueError, "prec must be positive"
    175 ##     v = [0] * prec
    176 ##     stop = int((-1+math.sqrt(1+8*prec))/2.0)
    177 ##     for n in range(stop+1):
    178 ## ##
    179 ## {{{
    180 ## pari('eta(q+O(q^200))')
    181 ## ///
    182 ## 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)
    183 ## }}}
    184 ## {{{
    185 ## sloane_find([2,5,7,12,15,22,26])
    186 ## ///
    187 ## Searching Sloane's online database...
    188 ## [[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]]]
    189 ## }}}
    190 
    191 from sage.rings.all import Integer
    192 
    193174def delta_qexp(prec=10, var='q', K=ZZ):
    194175    """
    195176    Return the q-expansion of Delta as a power series with