summaryrefslogtreecommitdiffstats
diff options
context:
space:
mode:
authorFrédéric Chapoton <chapoton@math.univ-lyon1.fr>2015-06-27 23:15:07 +0200
committerFrédéric Chapoton <chapoton@math.univ-lyon1.fr>2015-06-27 23:15:07 +0200
commit23d2ce8e1466ed8c5612d13564c21a20d147b035 (patch)
treebb8399c71cf93d1c51461d61996fafa8ba185385
parentMerge branch 'public/ticket/16448' into 6.8.b6 (diff)
trac #16448 pyflakes and pep8 work
-rw-r--r--src/sage/modular/jacobi/all.py7
-rw-r--r--src/sage/modular/jacobi/classical.py42
-rw-r--r--src/sage/modular/jacobi/classical_weak.py387
-rw-r--r--src/sage/modular/jacobi/higherrank.py120
-rw-r--r--src/sage/modular/jacobi/higherrank_dimension.py158
-rw-r--r--src/sage/modular/jacobi/test_classical.py10
-rw-r--r--src/sage/modular/jacobi/test_classical_weak.py94
-rw-r--r--src/sage/modular/jacobi/test_higherrank.py39
-rw-r--r--src/sage/modular/jacobi/test_higherrank_dimension.py44
-rw-r--r--src/sage/modular/jacobi/test_vector_valued.py79
-rw-r--r--src/sage/modular/jacobi/vector_valued.py136
11 files changed, 569 insertions, 547 deletions
diff --git a/src/sage/modular/jacobi/all.py b/src/sage/modular/jacobi/all.py
index 89bb9d9..744ff8d 100644
--- a/src/sage/modular/jacobi/all.py
+++ b/src/sage/modular/jacobi/all.py
@@ -1,13 +1,11 @@
from classical_weak import (
classical_jacobi_reduce_fe_index,
classical_weak_jacobi_fe_indices,
-
- classical_weak_jacobi_forms
+ classical_weak_jacobi_forms
)
from classical import (
classical_jacobi_fe_indices,
-
classical_jacobi_forms
)
@@ -15,7 +13,6 @@ from higherrank import (
higherrank_jacobi_fe_indices,
higherrank_jacobi_reduce_fe_index,
higherrank_jacobi_r_classes,
-
higherrank_jacobi_forms
)
@@ -23,8 +20,6 @@ from vector_valued import (
vector_valued_modular_forms,
vector_valued_modular_forms_weakly_holomorphic,
vector_valued_modular_forms_weakly_holomorphic_with_principal_part,
-
theta_decomposition,
-
stably_equivalent_positive_definite_quadratic_form
)
diff --git a/src/sage/modular/jacobi/classical.py b/src/sage/modular/jacobi/classical.py
index 401491a..e9834ba 100644
--- a/src/sage/modular/jacobi/classical.py
+++ b/src/sage/modular/jacobi/classical.py
@@ -46,19 +46,19 @@ classical weak Jacobi forms.
"""
#===============================================================================
-#
+#
# Copyright (C) 2010-2014 Martin Raum
-#
+#
# This program is free software; you can redistribute it and/or
# modify it under the terms of the GNU General Public License
# as published by the Free Software Foundation; either version 3
# of the License, or (at your option) any later version.
-#
-# This program is distributed in the hope that it will be useful,
+#
+# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
# General Public License for more details.
-#
+#
# You should have received a copy of the GNU General Public License
# along with this program; if not, see <http://www.gnu.org/licenses/>.
#
@@ -72,7 +72,8 @@ from sage.rings.all import Integer, ZZ
# We do not implement this separately, because this is the same
# reduction as in the case of weak Jacobi forms.
-from sage.modular.jacobi.classical_weak import classical_jacobi_reduce_fe_index
+# from sage.modular.jacobi.classical_weak import classical_jacobi_reduce_fe_index
+
def classical_jacobi_fe_indices(m, prec, reduced=False):
r"""
@@ -102,7 +103,7 @@ def classical_jacobi_fe_indices(m, prec, reduced=False):
"""
fm = Integer(4 * m)
- if reduced :
+ if reduced:
# positive definite forms
for n in range(1, prec):
for r in range(min(m + 1, isqrt(fm * n - 1) + 1)):
@@ -110,9 +111,9 @@ def classical_jacobi_fe_indices(m, prec, reduced=False):
# indefinite forms
for r in xrange(0, min(m+1, isqrt((prec-1) * fm) + 1) ):
- if fm.divides(r**2) :
- yield (r**2 // fm, r)
- else :
+ if fm.divides(r ** 2):
+ yield (r ** 2 // fm, r)
+ else:
# positive definite forms
for n in range(1, prec):
yield(n, 0)
@@ -122,25 +123,26 @@ def classical_jacobi_fe_indices(m, prec, reduced=False):
# indefinite forms
if prec > 0:
- yield (0,0)
+ yield (0, 0)
for n in xrange(1, prec):
if (fm * n).is_square():
rt_fmm = isqrt(fm * n)
yield(n, rt_fmm)
yield(n, -rt_fmm)
-
+
raise StopIteration
+
@cached_function
def _classical_jacobi_forms_as_weak_jacobi_forms(k, m, algorithm="skoruppa"):
r"""
The coordinates of Jacobi forms with respect to a basis of
classical weak Jacobi forms computed using the given algorithm.
-
+
INPUT:
-
+
- `k` -- An integer.
-
+
- `m` -- A non-negative integer.
- ``algorithm`` -- Default: ''skoruppa''. Only ''skoruppa'' is implemented.
@@ -148,9 +150,9 @@ def _classical_jacobi_forms_as_weak_jacobi_forms(k, m, algorithm="skoruppa"):
OUTPUT:
A list of vectors.
-
+
EXAMPLES::
-
+
sage: from sage.modular.jacobi.classical import _classical_jacobi_forms_as_weak_jacobi_forms
sage: _classical_jacobi_forms_as_weak_jacobi_forms(10, 1)
[
@@ -169,7 +171,7 @@ def _classical_jacobi_forms_as_weak_jacobi_forms(k, m, algorithm="skoruppa"):
weak_index_matrix = \
matrix(ZZ, [ [ (f[(n,r)] if (n,r) in f else 0) for (n,r) in indices
if 4*m*n - r**2 < 0 ] for f in weak_forms] )
-
+
return weak_index_matrix.left_kernel().echelonized_basis()
def classical_jacobi_forms(k, m, prec, algorithm="skoruppa"):
@@ -178,9 +180,9 @@ def classical_jacobi_forms(k, m, prec, algorithm="skoruppa"):
and index `m`.
INPUT:
-
+
- `k` -- An integer.
-
+
- `m` -- A non-negative integer.
- ``prec`` -- A non-negative integer that corresponds to a
diff --git a/src/sage/modular/jacobi/classical_weak.py b/src/sage/modular/jacobi/classical_weak.py
index b1a7aca..f44f3c7 100644
--- a/src/sage/modular/jacobi/classical_weak.py
+++ b/src/sage/modular/jacobi/classical_weak.py
@@ -51,19 +51,19 @@ coefficient of the first basis element is equal to::
"""
#===============================================================================
-#
+#
# Copyright (C) 2010-2014 Martin Raum
-#
+#
# This program is free software; you can redistribute it and/or
# modify it under the terms of the GNU General Public License
# as published by the Free Software Foundation; either version 3
# of the License, or (at your option) any later version.
-#
-# This program is distributed in the hope that it will be useful,
+#
+# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
# General Public License for more details.
-#
+#
# You should have received a copy of the GNU General Public License
# along with this program; if not, see <http://www.gnu.org/licenses/>.
#
@@ -76,7 +76,7 @@ from sage.misc.all import cached_function
from sage.misc.all import prod, isqrt
from sage.misc.weak_dict import WeakValueDictionary
from sage.modular.all import ModularForms
-from sage.rings.all import (PowerSeriesRing, Integer, ZZ, QQ, GF)
+from sage.rings.all import (PowerSeriesRing, Integer, ZZ, QQ)
from sage.rings import big_oh
import operator
@@ -111,9 +111,9 @@ def classical_jacobi_reduce_fe_index((n, r), m):
rred = 2*m - rred
else:
sgn = 1
-
+
nred = n - (r**2 - rred**2)//(4*m)
-
+
return ((nred, rred), sgn)
def classical_weak_jacobi_fe_indices(m, prec, reduced=False):
@@ -143,17 +143,17 @@ def classical_weak_jacobi_fe_indices(m, prec, reduced=False):
"""
fm = Integer(4*m)
- if reduced :
+ if reduced:
# positive definite forms
for n in range(1, prec):
- for r in range(min(m + 1, isqrt(fm * n - 1) + 1)) :
+ for r in range(min(m + 1, isqrt(fm * n - 1) + 1)):
yield (n, r)
# indefinite forms
- for n in xrange(0, prec) :
- for r in xrange( isqrt(fm * n - 1) + 1 if n != 0 else 0, m + 1 ) :
+ for n in xrange(0, prec):
+ for r in xrange( isqrt(fm * n - 1) + 1 if n != 0 else 0, m + 1 ):
yield (n, r)
- else :
+ else:
# positive definite forms
for n in range(1, prec):
yield(n, 0)
@@ -162,26 +162,27 @@ def classical_weak_jacobi_fe_indices(m, prec, reduced=False):
yield (n, -r)
# indefinite forms
- for n in xrange(0, min(m // 4 + 1, prec)) :
- if n == 0 :
+ for n in xrange(0, min(m // 4 + 1, prec)):
+ if n == 0:
r_iteration = range(-m + 1, m + 1)
- else :
- r_iteration = range( -m + 1, -isqrt(fm * n - 1) ) \
- + range( isqrt(fm * n - 1) + 1, m + 1 )
- for r in r_iteration :
- for l in range( (- r - isqrt(r**2 - 4 * m * (n - (prec - 1))) - 1) // (2 * m) + 1,
- (- r + isqrt(r**2 - 4 * m * (n - (prec - 1)))) // (2 * m) + 1 ) :
- yield (n + l * r + m * l**2, r + 2 * m * l)
+ else:
+ r_iteration = range(-m + 1, -isqrt(fm * n - 1))
+ r_iteration += range(isqrt(fm * n - 1) + 1, m + 1)
+ for r in r_iteration:
+ for l in range((- r - isqrt(r**2 - 4 * m * (n - (prec - 1))) - 1) // (2 * m) + 1,
+ (- r + isqrt(r**2 - 4 * m * (n - (prec - 1)))) // (2 * m) + 1 ):
+ yield (n + l * r + m * l ** 2, r + 2 * m * l)
raise StopIteration
+
## TODO: Implement caching by truncation.
-def classical_weak_jacobi_forms(k, m, prec, algorithm="skoruppa") :
+def classical_weak_jacobi_forms(k, m, prec, algorithm="skoruppa"):
r"""
INPUT:
-
+
- `k` -- An integer.
-
+
- `m` -- A non-negative integer.
- ``prec`` -- A non-negative integer that corresponds to a precision of
@@ -215,47 +216,49 @@ def classical_weak_jacobi_forms(k, m, prec, algorithm="skoruppa") :
return [ factory.from_taylor_expansion(fs, k, is_integral=True)
for fs in _classical_weak_jacobi_taylor_coefficients(k, m)]
+
@cached_function
-def _classical_weak_jacobi_taylor_coefficients(k, m) :
+def _classical_weak_jacobi_taylor_coefficients(k, m):
r"""
- A product basis of the echelon bases of
-
+ A product basis of the echelon bases of
+
- `M_k, M_{k + 2}, ..., M_{k + 2 m}` etc. if ``weight`` is even,
-
+
- `M_{k + 1}, ..., M_{k + 2 m - 3}` if ``weight`` is odd.
-
+
INPUT:
-
+
- `k` -- An integer.
-
+
- `m` -- A non-negative integer.
OUTPUT:
- A list of lists of functions. When called with with an arguement
+ A list of lists of functions. When called with with an argument
``prec``, each of them returns a `q`-expansion of that precision.
-
+
EXAMPLES::
-
+
sage: from sage.modular.jacobi.classical_weak import _classical_weak_jacobi_taylor_coefficients
sage: _classical_weak_jacobi_taylor_coefficients(12, 1)
[[<bound method ModularFormElement.qexp of 1 + 196560*q^2 + 16773120*q^3 + 398034000*q^4 + 4629381120*q^5 + O(q^6)>, <function <lambda> at ...>], [<bound method ModularFormElement.qexp of q - 24*q^2 + 252*q^3 - 1472*q^4 + 4830*q^5 + O(q^6)>, <function <lambda> at ...>], [<function <lambda> at ...>, <bound method ModularFormElement.qexp of 1 - 24*q - 196632*q^2 - 38263776*q^3 - 1610809368*q^4 - 29296875024*q^5 + O(q^6)>]]
"""
- R = PowerSeriesRing(ZZ, 'q'); q = R.gen()
-
- if k % 2 == 0 :
+ R = PowerSeriesRing(ZZ, 'q')
+ q = R.gen()
+
+ if k % 2 == 0:
nmb_modular_forms = m + 1
start_weight = k
- else :
+ else:
nmb_modular_forms = m - 1
start_weight = k + 1
-
+
modular_forms = list()
- for (i,k) in enumerate(range(start_weight, start_weight + 2 * nmb_modular_forms, 2)) :
+ for (i, k) in enumerate(range(start_weight, start_weight + 2 * nmb_modular_forms, 2)):
modular_forms += [ [lambda p: big_oh.O(q**p) for _ in range(i)] + [b.qexp] + [lambda p: big_oh.O(q**p) for _ in range(nmb_modular_forms - 1 - i)]
for b in ModularForms(1, k).echelon_basis() ]
-
- return modular_forms
+
+ return modular_forms
_classical_weak_jacobi_forms_factory_cache = WeakValueDictionary()
@@ -288,15 +291,15 @@ def ClassicalWeakJacobiFormsFactory(m, prec):
_classical_weak_jacobi_forms_factory_cache[(m, prec)] = fcty
return fcty
+
class ClassicalWeakJacobiForms_factory:
r"""
A factory for Jacobi form of degree 1 and scalar index `m`.
"""
-
- def __init__(self, m, prec) :
+ def __init__(self, m, prec):
r"""
INPUT:
-
+
- `m` -- A non-negative integer.
- ``prec`` -- An integer.
@@ -308,11 +311,11 @@ class ClassicalWeakJacobiForms_factory:
"""
self.__prec = prec
self.__jacobi_index = m
-
+
self.__power_series_ring_ZZ = PowerSeriesRing(ZZ, 'q')
self.__power_series_ring = PowerSeriesRing(QQ, 'q')
- def jacobi_index(self) :
+ def jacobi_index(self):
r"""
The index of the Jacobi forms that are computed by this factory.
@@ -329,7 +332,7 @@ class ClassicalWeakJacobiForms_factory:
"""
return self.__jacobi_index
- def precision(self) :
+ def precision(self):
r"""
The precision of Fourier expansions that are computed.
@@ -346,7 +349,7 @@ class ClassicalWeakJacobiForms_factory:
"""
return self.__prec
- def _power_series_ring(self) :
+ def _power_series_ring(self):
r"""
An auxiliary power series ring that is cached in the factory.
@@ -362,8 +365,8 @@ class ClassicalWeakJacobiForms_factory:
Power Series Ring in q over Rational Field
"""
return self.__power_series_ring
-
- def _power_series_ring_ZZ(self) :
+
+ def _power_series_ring_ZZ(self):
r"""
An auxiliary power series ring over `\ZZ` that is cached in the factory.
@@ -409,8 +412,7 @@ class ClassicalWeakJacobiForms_factory:
pass
try:
- other._wronskian_adjoint_odd = \
- matrix([[e.truncate(prec) for e in r]
+ other._wronskian_adjoint_odd = matrix([[e.truncate(prec) for e in r]
for r in self._wronskian_adjoint_odd.rows()])
except AttributeError:
pass
@@ -429,14 +431,14 @@ class ClassicalWeakJacobiForms_factory:
return other
- def _set_wronskian_adjoint(self, wronskian_adjoint, weight_parity = 0) :
+ def _set_wronskian_adjoint(self, wronskian_adjoint, weight_parity=0):
r"""
Set the cache for the adjoint of the wronskian. See _wronskian_adjoint.
-
+
INPUT:
-
+
- ``wronskian_adjoint`` -- A list of lists of power series over `\ZZ`.
-
+
- ``weight_parity`` -- An integer (default: `0`).
TESTS::
@@ -452,24 +454,24 @@ class ClassicalWeakJacobiForms_factory:
wronskian_adjoint = [ [ e if e in self.__power_series_ring_ZZ else self.__power_series_ring_ZZ(e)
for e in row ]
for row in wronskian_adjoint ]
-
- if weight_parity % 2 == 0 :
+
+ if weight_parity % 2 == 0:
self._wronskian_adjoint_even = wronskian_adjoint
- else :
+ else:
self._wronskian_adjoint_odd = wronskian_adjoint
-
- def _wronskian_adjoint(self, weight_parity = 0) :
+
+ def _wronskian_adjoint(self, weight_parity=0):
r"""
The matrix `W^\# \pmod{p}`, mentioned on page 142 of [Sk84]_.
This matrix is represented by a list of lists of q-expansions.
-
+
The q-expansion is shifted by `-(m + 1) (2*m + 1) / 24` in the
case of even weights, and it is shifted by `-(m - 1) (2*m - 3)
/ 24` otherwise. This is compensated by the missing q-powers
returned by _wronskian_invdeterminant.
-
+
INPUT:
-
+
- ``weight_parity`` -- An integer (default: `0`).
OUTPUT:
@@ -491,80 +493,79 @@ class ClassicalWeakJacobiForms_factory:
-2 - 162*q + 1020*q^2 - 550*q^3 + O(q^5),
2 + 18*q - 60*q^2 + 22*q^3 + O(q^5)]]
"""
- try :
- if weight_parity % 2 == 0 :
+ try:
+ if weight_parity % 2 == 0:
return self._wronskian_adjoint_even
- else :
+ else:
return self._wronskian_adjoint_ood
- except AttributeError :
+ except AttributeError:
PS = self.__power_series_ring_ZZ
m = self.jacobi_index()
-
+
twom = 2 * m
- frmsq = twom ** 2
-
+
thetas = dict( ((i, j), dict())
for i in xrange(m + 1) for j in xrange(m + 1) )
## We want to calculate \hat \theta_{j,l} = sum_r (2 m r + j)**2l q**(m r**2 + j r)
## in the case of even weight, and \hat \theta_{j,l} = sum_r (2 m r + j)**(2l + 1) q**(m r**2 + j r),
- ## otherwise.
- for r in xrange(isqrt((self.__prec - 1 + m)//m) + 2) :
- for j in (xrange(m + 1) if weight_parity % 2 == 0 else range(1, m)) :
- fact_p = (twom*r + j)**2
- fact_m = (twom*r - j)**2
- if weight_parity % 2 == 0 :
+ ## otherwise.
+ for r in xrange(isqrt((self.__prec - 1 + m)//m) + 2):
+ for j in (xrange(m + 1) if weight_parity % 2 == 0 else range(1, m)):
+ fact_p = (twom * r + j) ** 2
+ fact_m = (twom * r - j) ** 2
+ if weight_parity % 2 == 0:
coeff_p = 2
coeff_m = 2
- else :
- coeff_p = 2 * (twom*r + j)
- coeff_m = - 2 * (twom*r - j)
-
- for l in (xrange(m + 1) if weight_parity % 2 == 0 else range(1, m)) :
+ else:
+ coeff_p = 2 * (twom * r + j)
+ coeff_m = - 2 * (twom * r - j)
+
+ for l in (xrange(m + 1) if weight_parity % 2 == 0 else range(1, m)):
thetas[(j,l)][m*r**2 + r*j] = coeff_p
thetas[(j,l)][m*r**2 - r*j] = coeff_m
coeff_p = coeff_p * fact_p
coeff_m = coeff_m * fact_m
- if weight_parity % 2 == 0 :
- thetas[(0,0)][0] = 1
-
- thetas = dict( ( k, PS(th).add_bigoh(self.__prec) )
- for (k,th) in thetas.iteritems() )
-
+ if weight_parity % 2 == 0:
+ thetas[(0, 0)][0] = 1
+
+ thetas = dict((k, PS(th).add_bigoh(self.__prec))
+ for (k,th) in thetas.iteritems())
+
W = matrix( PS, m + 1 if weight_parity % 2 == 0 else (m - 1),
[ thetas[(j, l)]
for j in (xrange(m + 1) if weight_parity % 2 == 0 else range(1, m))
for l in (xrange(m + 1) if weight_parity % 2 == 0 else range(1, m)) ] )
-
-
+
+
## Since the adjoint of matrices with entries in a general ring
## is extremely slow for matrices of small size, we hard code the
## the cases `m = 2` and `m = 3`. The expressions are obtained by
## computing the adjoint of a matrix with entries `w_{i,j}` in a
## polynomial algebra.
- if W.nrows() == 1 :
- Wadj = matrix(PS, [[ 1 ]])
- elif W.nrows() == 2 :
- Wadj = matrix(PS, [ [ W[1,1], -W[0,1]],
- [-W[1,0], W[0,0]] ])
-
- elif W.nrows() == 3 and self.__prec > 10**5 :
- adj00 = W[1,1] * W[2,2] - W[2,1] * W[1,2]
+ if W.nrows() == 1:
+ Wadj = matrix(PS, [[1]])
+ elif W.nrows() == 2:
+ Wadj = matrix(PS, [[W[1, 1], -W[0, 1]],
+ [-W[1, 0], W[0, 0]]])
+
+ elif W.nrows() == 3 and self.__prec > 10 ** 5:
+ adj00 = W[1, 1] * W[2, 2] - W[2, 1] * W[1, 2]
adj01 = - W[1,0] * W[2,2] + W[2,0] * W[1,2]
- adj02 = W[1,0] * W[2,1] - W[2,0] * W[1,1]
+ adj02 = W[1,0] * W[2,1] - W[2,0] * W[1,1]
adj10 = - W[0,1] * W[2,2] + W[2,1] * W[0,2]
- adj11 = W[0,0] * W[2,2] - W[2,0] * W[0,2]
+ adj11 = W[0,0] * W[2,2] - W[2,0] * W[0,2]
adj12 = - W[0,0] * W[2,1] + W[2,0] * W[0,1]
- adj20 = W[0,1] * W[1,2] - W[1,1] * W[0,2]
+ adj20 = W[0,1] * W[1,2] - W[1,1] * W[0,2]
adj21 = - W[0,0] * W[1,2] + W[1,0] * W[0,2]
- adj22 = W[0,0] * W[1,1] - W[1,0] * W[0,1]
+ adj22 = W[0,0] * W[1,1] - W[1,0] * W[0,1]
+
+ Wadj = matrix(PS, [[adj00, adj01, adj02],
+ [adj10, adj11, adj12],
+ [adj20, adj21, adj22]])
- Wadj = matrix(PS, [ [adj00, adj01, adj02],
- [adj10, adj11, adj12],
- [adj20, adj21, adj22] ])
-
- elif W.nrows() == 4 and self.__prec > 10**5 :
+ elif W.nrows() == 4 and self.__prec > 10 ** 5:
adj00 = -W[0,2]*W[1,1]*W[2,0] + W[0,1]*W[1,2]*W[2,0] + W[0,2]*W[1,0]*W[2,1] - W[0,0]*W[1,2]*W[2,1] - W[0,1]*W[1,0]*W[2,2] + W[0,0]*W[1,1]*W[2,2]
adj01 = -W[0,3]*W[1,1]*W[2,0] + W[0,1]*W[1,3]*W[2,0] + W[0,3]*W[1,0]*W[2,1] - W[0,0]*W[1,3]*W[2,1] - W[0,1]*W[1,0]*W[2,3] + W[0,0]*W[1,1]*W[2,3]
adj02 = -W[0,3]*W[1,2]*W[2,0] + W[0,2]*W[1,3]*W[2,0] + W[0,3]*W[1,0]*W[2,2] - W[0,0]*W[1,3]*W[2,2] - W[0,2]*W[1,0]*W[2,3] + W[0,0]*W[1,2]*W[2,3]
@@ -589,31 +590,34 @@ class ClassicalWeakJacobiForms_factory:
[adj10, adj11, adj12, adj13],
[adj20, adj21, adj22, adj23],
[adj30, adj31, adj32, adj33] ])
- else :
+ else:
Wadj = W.adjoint()
-
- if weight_parity % 2 == 0 :
- wronskian_adjoint = [ [ Wadj[i,r] for i in xrange(m + 1) ]
- for r in xrange(m + 1) ]
- else :
- wronskian_adjoint = [ [ Wadj[i,r] for i in xrange(m - 1) ]
- for r in xrange(m - 1) ]
-
- if weight_parity % 2 == 0 :
+
+ if weight_parity % 2 == 0:
+ wronskian_adjoint = [[Wadj[j, r] for j in xrange(m + 1)]
+ for r in xrange(m + 1)]
+ else:
+ wronskian_adjoint = [[Wadj[j, r] for j in xrange(m - 1)]
+ for r in xrange(m - 1)]
+
+ if weight_parity % 2 == 0:
self._wronskian_adjoint_even = wronskian_adjoint
- else :
+ else:
self._wronskian_adjoint_odd = wronskian_adjoint
-
+
return wronskian_adjoint
-
- def _set_wronskian_invdeterminant(self, wronskian_invdeterminant, weight_parity = 0) :
+
+ def _set_wronskian_invdeterminant(self, wronskian_invdeterminant,
+ weight_parity=0):
r"""
- Set the cache for the inverse determinant of the Wronskian. See _wronskian_adjoint.
-
+ Set the cache for the inverse determinant of the Wronskian.
+
+ See _wronskian_adjoint.
+
INPUT:
-
+
- ``wronskian_invdeterminant`` -- A power series over `\ZZ`.
-
+
- ``weight_parity`` -- An integer (default: `0`).
TESTS::
@@ -624,21 +628,23 @@ class ClassicalWeakJacobiForms_factory:
sage: fcty._set_wronskian_invdeterminant(invdet, 0)
sage: assert fcty._wronskian_invdeterminant(0) is invdet
"""
- if not wronskian_invdeterminant in self.__power_series_ring_ZZ :
+ if not wronskian_invdeterminant in self.__power_series_ring_ZZ:
wronskian_invdeterminant = self.__power_series_ring_ZZ(wronskian_invdeterminant)
-
- if weight_parity % 2 == 0 :
+
+ if weight_parity % 2 == 0:
self._wronskian_invdeterminant_even = wronskian_invdeterminant
- else :
+ else:
self._wronskian_invdeterminant_odd = wronskian_invdeterminant
-
- def _wronskian_invdeterminant(self, weight_parity = 0) :
+
+ def _wronskian_invdeterminant(self, weight_parity=0):
r"""
The inverse determinant of `W`, which in the considered cases
- is always a negative power of the eta function. See [Sk84]_ for details.
-
+ is always a negative power of the eta function.
+
+ See [Sk84]_ for details.
+
INPUT:
-
+
- ``weight_parity`` -- An integer (default: `0`).
OUTPUT:
@@ -652,31 +658,32 @@ class ClassicalWeakJacobiForms_factory:
sage: fcty._wronskian_invdeterminant(0)
1 + 15*q + 135*q^2 + 920*q^3 + 5220*q^4 + O(q^5)
"""
- try :
- if weight_parity % 2 == 0 :
+ try:
+ if weight_parity % 2 == 0:
return self._wronskian_invdeterminant_even
- else :
+ else:
return self._wronskian_invdeterminant_odd
- except AttributeError :
+ except AttributeError:
m = self.jacobi_index()
- if weight_parity % 2 == 0 :
+ if weight_parity % 2 == 0:
pw = (m + 1) * (2 * m + 1)
- else :
+ else:
pw = (m - 1) * (2 * m - 1)
-
- wronskian_invdeterminant = self.__power_series_ring_ZZ \
- ( [ number_of_partitions(n) for n in xrange(self.__prec) ] ) \
- .add_bigoh(self.__prec) ** pw
-
- if weight_parity % 2 == 0 :
+
+ eta = self.__power_series_ring_ZZ([number_of_partitions(n)
+ for n in xrange(self.__prec)]).add_bigoh(self.__prec)
+
+ wronskian_invdeterminant = eta ** pw
+
+ if weight_parity % 2 == 0:
self._wronskian_invdeterminant_even = wronskian_invdeterminant
- else :
+ else:
self._wronskian_invdeterminant_odd = wronskian_invdeterminant
return wronskian_invdeterminant
-
- def from_taylor_expansion(self, fs, k, is_integral=False) :
+
+ def from_taylor_expansion(self, fs, k, is_integral=False):
r"""
The Fourier expansion of a weak Jacobi form computed from its
corrected Taylor coefficients.
@@ -688,22 +695,22 @@ class ClassicalWeakJacobiForms_factory:
forms of weight `k` and index `m` with the product of spaces
of elliptic modular forms `M_k \times S_{k+2} \times .. \times
S_{k+2m}`.
-
+
INPUT:
-
+
- ``fs`` -- A list of functions that given an integer `p` return the
q-expansion of a modular form with rational coefficients
up to precision `p`. These modular forms correspond to
the components of the above product.
-
+
- `k` -- An integer. The weight of the weak Jacobi form to be computed.
-
+
- ``is_integral`` -- A boolean. If ``True``, the ``fs`` have integral
coefficients.
TESTS::
-
- sage: from sage.modular.jacobi.classical_weak import *
+
+ sage: from sage.modular.jacobi.classical_weak import *
sage: factory = ClassicalWeakJacobiFormsFactory(3, 10)
sage: R.<q> = ZZ[[]]
sage: expansion = factory.from_taylor_expansion([lambda p: 0 + O(q^p), lambda p: CuspForms(1, 12).gen(0).qexp(p)], 9, True)
@@ -711,62 +718,64 @@ class ClassicalWeakJacobiForms_factory:
sage: sorted([ (12 * n - r**2, c/exp_gcd) for ((n, r), c) in expansion.items()])
[(-4, 0), (-1, 0), (8, 1), (11, -2), (20, -14), (23, 32), (32, 72), (35, -210), (44, -112), (47, 672), (56, -378), (59, -728), (68, 1736), (71, -1856), (80, -1008), (83, 6902), (92, -6400), (95, -5792), (104, 10738), (107, -6564)]
"""
- if is_integral :
+ if is_integral:
PS = self.__power_series_ring_ZZ
- else :
+ else:
PS = self.__power_series_ring
-
+
if (k % 2 == 0 and not len(fs) == self.jacobi_index() + 1) \
- or (k % 2 != 0 and not len(fs) == self.jacobi_index() - 1) :
- raise ValueError( "fs (which has length {0}) must be a list of {1} Fourier expansions" \
- .format(len(fs), self.jacobi_index() + 1 if k % 2 == 0 else self.jacobi_index() - 1) )
-
- if self.__prec <= 0: # there are no Fourier indices below the precision
+ or (k % 2 != 0 and not len(fs) == self.jacobi_index() - 1):
+ msg = "fs (which has length {0}) must be a list of {1} Fourier expansions"
+ raise ValueError(msg.format(len(fs),
+ self.jacobi_index() + 1 if k % 2 == 0
+ else self.jacobi_index() - 1))
+
+ if self.__prec <= 0: # there are no Fourier indices below the precision
return dict()
-
+
f_divs = dict()
- for (i, f) in enumerate(fs) :
+ for (i, f) in enumerate(fs):
f_divs[(i, 0)] = PS(f(self.__prec), self.__prec + 10)
m = self.jacobi_index()
-
- for i in (xrange(m + 1) if k % 2 == 0 else xrange(m - 1)) :
- for j in xrange(1, m - i + 1) :
+
+ for i in (xrange(m + 1) if k % 2 == 0 else xrange(m - 1)):
+ for j in xrange(1, m - i + 1):
f_divs[(i,j)] = f_divs[(i, j - 1)].derivative().shift(1)
-
+
phi_divs = list()
- for i in (xrange(m + 1) if k % 2 == 0 else xrange(m - 1)) :
- if k % 2 == 0 :
+ for i in (xrange(m + 1) if k % 2 == 0 else xrange(m - 1)):
+ if k % 2 == 0:
## This is (13) on page 131 of Skoruppa (change of variables n -> i, r -> j).
## The additional factor (2m + k - 1)! is a renormalization to make coefficients integral.
## The additional factor (4m)^i stems from the fact that we have used d / d\tau instead of
## d^2 / dz^2 when treating the theta series. Since these are annihilated by the heat operator
- ## the additional factor compensates for this.
+ ## the additional factor compensates for this.
phi_divs.append( sum( f_divs[(j, i - j)]
- * ( (4 * self.jacobi_index())**i
- * binomial(i,j) / 2**i # 2**(self.__jacobi_index - i + 1)
- * prod(2*l + 1 for l in xrange(i))
+ * ( (4 * self.jacobi_index()) ** i
+ * binomial(i, j) / 2 ** i # 2**(self.__jacobi_index - i + 1)
+ * prod(2 * l + 1 for l in xrange(i))
/ factorial(i + k + j - 1)
- * factorial(2*self.jacobi_index() + k - 1) )
+ * factorial(2*self.jacobi_index() + k - 1) )
for j in range(i + 1) ) )
- else :
+ else:
phi_divs.append( sum( f_divs[(j, i - j)]
- * ( (4 * self.jacobi_index())**i
- * binomial(i,j) / 2**(i + 1) # 2**(self.__jacobi_index - i + 1)
+ * ( (4 * self.jacobi_index()) ** i
+ * binomial(i, j) / 2 ** (i + 1) # 2**(self.__jacobi_index - i + 1)
* prod(2*l + 1 for l in xrange(i + 1))
/ factorial(i + k + j)
- * factorial(2*self.jacobi_index() + k - 1) )
+ * factorial(2*self.jacobi_index() + k - 1) )
for j in range(i + 1) ) )
-
+
phi_coeffs = dict()
- for r in (xrange(m + 1) if k % 2 == 0 else xrange(1, m)) :
- if k % 2 == 0 :
- series = sum( map(operator.mul, self._wronskian_adjoint(k)[r], phi_divs) )
- else :
- series = sum( map(operator.mul, self._wronskian_adjoint(k)[r - 1], phi_divs) )
+ for r in (xrange(m + 1) if k % 2 == 0 else xrange(1, m)):
+ if k % 2 == 0:
+ series = sum(map(operator.mul, self._wronskian_adjoint(k)[r], phi_divs))
+ else:
+ series = sum(map(operator.mul, self._wronskian_adjoint(k)[r - 1], phi_divs))
series = self._wronskian_invdeterminant(k) * series
- for n in xrange(self.__prec) :
+ for n in xrange(self.__prec):
phi_coeffs[(n, r)] = series[n]
return phi_coeffs
diff --git a/src/sage/modular/jacobi/higherrank.py b/src/sage/modular/jacobi/higherrank.py
index 316df62..3617924 100644
--- a/src/sage/modular/jacobi/higherrank.py
+++ b/src/sage/modular/jacobi/higherrank.py
@@ -100,7 +100,7 @@ from sage.misc.flatten import flatten
from sage.modular.jacobi.classical import (classical_jacobi_forms,
classical_jacobi_fe_indices, classical_jacobi_reduce_fe_index)
from sage.modular.jacobi.higherrank_dimension import jacobi_dimension
-from sage.modules.all import FreeModule, vector, span, zero_vector
+from sage.modules.all import FreeModule, vector, zero_vector
from sage.rings.all import ZZ, QQ
from sage.quadratic_forms.all import QuadraticForm
from sage.sets.all import Set
@@ -196,7 +196,7 @@ def _higherrank_jacobi_reduce_fe_index__r(r, r_classes, m_span):
return (rred, 1)
if vector(r) + vector(rred) in m_span:
return (rred, -1)
- else :
+ else:
raise RuntimeError( "Could not find reduced r" )
def higherrank_jacobi_fe_indices(m, prec, r_classes, reduced=False):
@@ -244,11 +244,11 @@ def higherrank_jacobi_fe_indices(m, prec, r_classes, reduced=False):
r = r_class[0]
if m_adj(r) <= 2*m.det()*n:
yield (n, r)
- else :
+ else:
short_vectors = m_adj.short_vector_list_up_to_length(2*m.det()*(prec - 1) + 1)
for n in range(0, prec):
- for length in range(0, 2*m.det()*n + 1):
- for r in short_vectors[length] :
+ for length in range(2 * m.det() * n + 1):
+ for r in short_vectors[length]:
yield (n, tuple(r))
raise StopIteration
@@ -282,8 +282,7 @@ def higherrank_jacobi_r_classes(m):
m_span = m_mat.row_module()
m_adj = QuadraticForm(2 * m_mat.adjoint())
-
- canonical_reps = [r.lift() for r in m_span.ambient_module() / m_span]
+ canonical_reps = [r.lift() for r in m_span.ambient_module() / m_span]
max_norm = max(m_adj(r) for r in canonical_reps)
current_max_length = 5
@@ -298,7 +297,8 @@ def higherrank_jacobi_r_classes(m):
or vector(r_class[0]) + vector(r_can) in m_span):
r_class_found = True
break
- if r_class_found: continue
+ if r_class_found:
+ continue
r_classes.append([])
r_class = r_classes[-1]
@@ -320,11 +320,13 @@ def higherrank_jacobi_r_classes(m):
r_class.append(tuple(r))
r_class_reduction_signs.append(-1)
- if len(r_class) != 0: break
+ if len(r_class) != 0:
+ break
for cl_ix in range(len(r_classes_reduction_signs)):
if r_classes_reduction_signs[cl_ix][0] == -1:
- r_classes_reduction_signs[cl_ix] = map(operator.neg, r_classes_reduction_signs[cl_ix])
+ r_classes_reduction_signs[cl_ix] = map(operator.neg,
+ r_classes_reduction_signs[cl_ix])
return (r_classes, r_classes_reduction_signs)
@@ -401,20 +403,17 @@ def higherrank_jacobi_forms(k, m, prec, algorithm="restriction"):
raise NotImplementedError("Algorithm {} is not implemented.".format(algorithm))
rand = Random()
-
dim = jacobi_dimension(k, m)
- if dim == 0: return []
-
+ if dim == 0:
+ return []
(r_classes, r_classes_reduction_signs) = higherrank_jacobi_r_classes(m)
m_span = m.matrix().row_module()
-
rst_vectors_with_image = _complete_set_of_restriction_vectors(m, r_classes, r_classes_reduction_signs, m_span)
rst_vectors = [vector(s)
for s in Set(tuple(s) for (s, _) in rst_vectors_with_image)]
-
max_rst_index = max([m(s) for s in rst_vectors])
minimal_prec = 2 + (k + max_rst_index) // 12 + 5
relation_prec = minimal_prec
@@ -448,18 +447,19 @@ def higherrank_jacobi_forms(k, m, prec, algorithm="restriction"):
relation_prec += minimal_prec
max_relation_rst_index += 2
-
- if len(relation_rst_vectors): relation_rst_vectors = rst_vectors
+ if len(relation_rst_vectors):
+ relation_rst_vectors = rst_vectors
short_vectors = (
flatten( m.short_vector_list_up_to_length(max_relation_rst_index+1, True)[max_relation_rst_index-2:],
- max_level = 1 )
+ max_level=1)
)
try:
relation_rst_vectors += rand.sample(short_vectors, m.det())
except ValueError:
relation_rst_vectors += short_vectors
+
def _complete_set_of_restriction_vectors(m, r_classes, r_classes_reduction_signs, m_span):
r"""
Given classes ``r_classes`` of elements in `\ZZ^l` find a complete
@@ -514,7 +514,7 @@ def _complete_set_of_restriction_vectors(m, r_classes, r_classes_reduction_signs
rst_kernel_odd = rst_kernel_odd_matrix.row_module()
- while (nmb_rst_vectors_even < len(r_classes)) :
+ while (nmb_rst_vectors_even < len(r_classes)):
while len(short_vectors[cur_length]) == 0:
cur_length += 1
if max_length < cur_length:
@@ -639,7 +639,8 @@ def _restriction_relation_matrices(k, m, prec, relation_prec,
return ( restriction_matrix__big, row_groups, row_labels, column_labels,
relation_matrix, column_labels_relations )
-def _restriction_matrix(k, m, prec, rst_vectors, find_relations, r_classes, m_span) :
+
+def _restriction_matrix(k, m, prec, rst_vectors, find_relations, r_classes, m_span):
r"""
A matrix that maps the Fourier expansion of a Jacobi form of given precision
to their restrictions with respect to the elements of S.
@@ -726,21 +727,19 @@ def _restriction_matrix(k, m, prec, rst_vectors, find_relations, r_classes, m_sp
if len(rst_vectors) == 0:
return (zero_matrix(ZZ, 0, len(column_labels)), [], {}, column_labels)
+ rst_jacobi_indices = [m(s) for s in rst_vectors]
+ rst_indices = dict((m_rst,
+ list(classical_jacobi_fe_indices(
+ m_rst, prec, reduced=not find_relations)))
+ for m_rst in Set(rst_jacobi_indices))
- rst_jacobi_indices = [ m(s) for s in rst_vectors ]
- rst_indices = dict( (m_rst,
- list(classical_jacobi_fe_indices(
- m_rst, prec, reduced = not find_relations)) )
- for m_rst in Set(rst_jacobi_indices) )
-
- row_groups = [ len(rst_indices[m_rst]) for m_rst in rst_jacobi_indices ]
- row_groups = [ (s, m_rst, sum(row_groups[:i]), row_groups[i])
- for ((i, s), m_rst) in zip(enumerate(rst_vectors), rst_jacobi_indices) ]
- row_labels = dict( (m_rst, dict( (nr, i) for (i, nr) in enumerate(rst_indices[m_rst]) ))
- for m_rst in Set(rst_jacobi_indices) )
-
+ row_groups = [len(rst_indices[mm_rst]) for mm_rst in rst_jacobi_indices]
+ row_groups = [(s, mm_rst, sum(row_groups[:i]), row_groups[i])
+ for ((i, s), mm_rst) in zip(enumerate(rst_vectors), rst_jacobi_indices) ]
+ row_labels = dict((m_rst, dict( (nr, i) for (i, nr) in enumerate(rst_indices[m_rst])))
+ for m_rst in Set(rst_jacobi_indices))
- reductions = dict( (nr,[]) for nr in column_labels )
+ reductions = dict((nr, []) for nr in column_labels)
for nr in higherrank_jacobi_fe_indices(m, prec, r_classes, reduced=False):
(nrred, sgn) = higherrank_jacobi_reduce_fe_index(nr, m, r_classes, m_adj, m_span)
reductions[nrred].append((nr, sgn))
@@ -754,18 +753,19 @@ def _restriction_matrix(k, m, prec, rst_vectors, find_relations, r_classes, m_sp
else:
cython_dot_products = False
-
restriction_matrix = zero_matrix(ZZ, row_groups[-1][2] + row_groups[-1][3],
len(column_labels))
for (col, nrred) in enumerate(column_labels):
- for ((n,r), sgn) in reductions[nrred]:
+ for ((n, r), sgn) in reductions[nrred]:
for (ix, (s, m_rst, start, length)) in enumerate(row_groups):
row_labels_dict = row_labels[m_rst]
- if cython_dot_products: rst_r = dot_products[ix](*r)
- else: rst_r = s.dot_product(vector(r))
+ if cython_dot_products:
+ rst_r = dot_products[ix](*r)
+ else:
+ rst_r = s.dot_product(vector(r))
- try :
+ try:
restriction_matrix[start + row_labels_dict[(n, rst_r)], col] \
+= 1 if k == 0 else sgn
except KeyError:
@@ -773,7 +773,8 @@ def _restriction_matrix(k, m, prec, rst_vectors, find_relations, r_classes, m_sp
return (restriction_matrix, row_groups, row_labels, column_labels)
-def _relation_matrix(k, m, prec, rst_vectors, r_classes, m_span) :
+
+def _relation_matrix(k, m, prec, rst_vectors, r_classes, m_span):
r"""
Deduce relations of the coefficients of a Jacobi form from their
specialization to Jacobi form of scalar index.
@@ -827,8 +828,7 @@ def _relation_matrix(k, m, prec, rst_vectors, r_classes, m_span) :
Tested implicitely by ``meth:higherrank_jacobi_forms``. See also ``test_higherrank.py``.
"""
k = k % 2
- m_adj = QuadraticForm(2 * m.matrix().adjoint())
-
+ # m_adj = QuadraticForm(2 * m.matrix().adjoint())
## relations computed from restrictions
(mat, row_groups, row_labels, column_labels) = \
@@ -839,24 +839,25 @@ def _relation_matrix(k, m, prec, rst_vectors, r_classes, m_span) :
row_labels_dict = row_labels[m_rst]
for (nr, ix) in row_labels_dict.items():
(nrred, sgn) = classical_jacobi_reduce_fe_index(nr, m_rst)
- if nrred == nr: continue
+ if nrred == nr:
+ continue
rel = (mat.row(start + row_labels_dict[nrred])
- (1 if k == 0 else sgn) * mat.row(start + ix))
- if rel != 0: relations.append(rel)
-
+ if rel != 0:
+ relations.append(rel)
## Forces zeros in case of odd k
- if k == 1:
- for (ix, (n,r)) in enumerate(column_labels):
- if 2*vector(r) in m_span:
+ if k == 1:
+ for (ix, (n, r)) in enumerate(column_labels):
+ if 2 * vector(r) in m_span:
rel = zero_vector(len(column_labels))
rel[ix] = 1
relations.append(rel)
-
return (matrix(len(relations), len(column_labels), relations), column_labels)
+
def _higherrank_jacobi_forms__restriction(
k, prec, relation_prec, dim,
restriction_matrix__big, row_groups, row_labels, column_labels,
@@ -869,7 +870,6 @@ def _higherrank_jacobi_forms__restriction(
- `k` -- An integer. Only `k` modulo `2` is used.
-
- ``prec`` -- A nonnegative integer.
- ``relation_prec`` -- A nonnegative integer.
@@ -919,11 +919,11 @@ def _higherrank_jacobi_forms__restriction(
row_labels__small = dict()
for (m_rst, row_labels_dict) in row_labels.items():
- row_labels_dict__small = dict()
+ row_labels_dict__small = {}
label_nmb = 0
for (nr, i) in row_labels_dict.items():
- if nr[0] < relation_prec :
+ if nr[0] < relation_prec:
row_labels_dict__small[nr] = (label_nmb, i)
label_nmb += 1
@@ -933,7 +933,7 @@ def _higherrank_jacobi_forms__restriction(
for ((s, m_rst, start, _), (start_small, length_small)) in zip(row_groups, row_groups__small):
row_labels_dict = row_labels__small[m_rst]
row_indices__sub = length_small * [None]
- for (_,(i, i_pre)) in row_labels_dict.items():
+ for (_, (i, i_pre)) in row_labels_dict.items():
row_indices__sub[i] = start + i_pre
row_indices__small += row_indices__sub
@@ -943,9 +943,9 @@ def _higherrank_jacobi_forms__restriction(
[column_labels.index(nr) for nr in column_labels_relations] )
- rst_jacobi_indices = [ m_rst for (_, m_rst, _, _) in row_groups ]
- rst_indices = dict( (m_rst, list(classical_jacobi_fe_indices(m_rst, prec, reduced=True)))
- for m_rst in Set(rst_jacobi_indices) )
+ rst_jacobi_indices = [m_rst for (_, m_rst, _, _) in row_groups]
+ # rst_indices = dict( (m_rst, list(classical_jacobi_fe_indices(m_rst, prec, reduced=True)))
+ # for m_rst in Set(rst_jacobi_indices) )
rst_jacobi_forms = dict( (m_rst, classical_jacobi_forms(k, m_rst, prec))
for m_rst in Set(rst_jacobi_indices) )
@@ -957,11 +957,10 @@ def _higherrank_jacobi_forms__restriction(
v = vector(ZZ, len(row_labels_dict))
for (nr, i) in row_labels_dict.items():
(nrred, sgn) = classical_jacobi_reduce_fe_index(nr, m_rst)
- v[i] = ((1 if k%2 == 0 else sgn) * phi[nr]) if nr in phi else 0
+ v[i] = ((1 if k % 2 == 0 else sgn) * phi[nr]) if nr in phi else 0
rst_jacobi_vectors.append(vector(
- start*[0] + v.list() + (nmb_rst_coords - start - length)*[0] ))
-
+ start * [0] + v.list() + (nmb_rst_coords - start - length) * [0] ))
rst_jacobi_matrix__big = \
matrix(len(rst_jacobi_vectors), nmb_rst_coords, rst_jacobi_vectors).transpose()
@@ -978,7 +977,7 @@ def _higherrank_jacobi_forms__restriction(
if jacobi_expansions_space.dimension() < dim:
raise RuntimeError( "There is a bug in the implementation of the restriction method. Dimensions: {} < {}!".format(jacobi_expansions_space.dimension(), dim) )
- if jacobi_expansions_space.dimension() > dim :
+ if jacobi_expansions_space.dimension() > dim:
raise ValueError( "Could not construct enough restrictions to determine Fourier expansions uniquely", "INSUFFICIENT RELATIONS" )
@@ -987,6 +986,5 @@ def _higherrank_jacobi_forms__restriction(
restriction_matrix * jacobi_expansions_space.basis_matrix().transpose() )
jacobi_expansions__big = restriction_matrix__big.solve_right( rst_jacobi_matrix__big * restriction_coordinates ).transpose()
-
- return [dict((nr,c) for (nr,c) in zip(column_labels, phi))
+ return [dict((nr, c) for (nr, c) in zip(column_labels, phi))
for phi in jacobi_expansions__big.rows()]
diff --git a/src/sage/modular/jacobi/higherrank_dimension.py b/src/sage/modular/jacobi/higherrank_dimension.py
index 18d870f..35a0073 100644
--- a/src/sage/modular/jacobi/higherrank_dimension.py
+++ b/src/sage/modular/jacobi/higherrank_dimension.py
@@ -1,4 +1,4 @@
-r"""
+r"""
A dimension formula for vector-valued modular forms, and functions
that apply it to the case of Jacobi forms.
@@ -16,31 +16,30 @@ AUTHOR:
"""
#===============================================================================
-#
+#
# Copyright (C) 2012-2014 Martin Raum
-#
+#
# This program is free software; you can redistribute it and/or
# modify it under the terms of the GNU General Public License
# as published by the Free Software Foundation; either version 3
# of the License, or (at your option) any later version.
-#
-# This program is distributed in the hope that it will be useful,
+#
+# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
# General Public License for more details.
-#
+#
# You should have received a copy of the GNU General Public License
# along with this program; if not, see <http://www.gnu.org/licenses/>.
#
#===============================================================================
-from sage.functions.all import exp, sqrt, sign
+from sage.functions.all import exp, sqrt
from sage.matrix.all import diagonal_matrix, identity_matrix, matrix
from sage.misc.all import sum, mrange, prod, cython_lambda
-from sage.modules.all import vector
+from sage.modules.all import vector
from sage.rings.all import ComplexIntervalField, ZZ, QQ, lcm
-from sage.rings.all import moebius, gcd, QuadraticField, fundamental_discriminant, kronecker_symbol
-from sage.quadratic_forms.all import QuadraticForm, BinaryQF_reduced_representatives
+from sage.quadratic_forms.all import QuadraticForm
from sage.symbolic.all import I, pi
from copy import copy
import operator
@@ -49,13 +48,13 @@ import operator
def jacobi_dimension(k, m):
r"""
INPUT:
-
+
- `k` -- An integer.
-
+
- `m` -- A quadratic form or an even symmetric matrix (over `\ZZ`).
-
+
TESTS::
-
+
sage: from sage.modular.jacobi.classical import _classical_jacobi_forms_as_weak_jacobi_forms
sage: from sage.modular.jacobi.higherrank_dimension import jacobi_dimension
sage: all(len(_classical_jacobi_forms_as_weak_jacobi_forms(k, m)) == jacobi_dimension(k, matrix([[2 * m]])) for k in range(8, 16) for m in range(1, 10) ) # long time
@@ -106,14 +105,14 @@ def vector_valued_dimension(k, L):
Compute the dimension of the space of weight `k` vector valued
modular forms for the Weil representation attached to the lattice
`L`.
-
+
See [Borcherds, Borcherds - Reflection groups of Lorentzian
lattices] for a proof of the formula that we use here.
-
+
INPUT:
-
+
- `k` -- A half-integer.
-
+
- ``L`` -- An quadratic form.
OUTPUT:
@@ -130,26 +129,26 @@ def vector_valued_dimension(k, L):
sage: vector_valued_dimension(3, QuadraticForm(-matrix(2, [2, 0, 0, 4])))
1
"""
- if 2 * k not in ZZ :
- raise ValueError( "Weight must be half-integral" )
- if k <= 0 :
+ if 2 * k not in ZZ:
+ raise ValueError( "Weight must be half-integral" )
+ if k <= 0:
return 0
- if k < 2 :
+ if k < 2:
raise NotImplementedError( "Weight <2 is not implemented." )
- if L.matrix().rank() != L.matrix().nrows() :
+ if L.matrix().rank() != L.matrix().nrows():
raise ValueError( "The lattice (={0}) must be non-degenerate.".format(L) )
- if L.dim() % 2 != ZZ(2 * k) % 2 :
+ if L.dim() % 2 != ZZ(2 * k) % 2:
return 0
(discriminant_form_exponents, disc_quadratic, disc_bilinear) = _discriminant_form(L)
(singls, pairs) = _discriminant_form_pmone(L, discriminant_form_exponents)
plus_basis = ZZ(L.dim() + 2*k) % 4 == 0
- if plus_basis :
+ if plus_basis:
subspace_dimension = len(singls + pairs)
- else :
+ else:
subspace_dimension = len(pairs)
CC_prec = 50 + subspace_dimension * 2
@@ -162,7 +161,7 @@ def vector_valued_dimension(k, L):
## This function overestimates the number of eigenvalues, if it is not correct
- def eigenvalue_multiplicity(mat, ev) :
+ def eigenvalue_multiplicity(mat, ev):
mat = matrix(CC, mat - ev * identity_matrix(subspace_dimension))
return len(filter( lambda row: all( e.contains_zero() for e in row), _qr(mat).rows() ))
@@ -219,32 +218,31 @@ def _discriminant_form(L):
See test_higherrank_dimension:test__discrimant_form.
"""
- if L.matrix().rank() != L.matrix().nrows() :
- raise ValueError( "The lattice (={0}) must be non-degenerate.".format(L) )
-
+ if L.matrix().rank() != L.matrix().nrows():
+ raise ValueError("The lattice (={0}) must be non-degenerate.".format(L))
## The bilinear and the quadratic form attached to L
- quadratic = lambda x: L(x) // 2
- bilinear = lambda x,y: L(x + y) - L(x) - L(y)
+ # quadratic = lambda x: L(x) // 2
+ # bilinear = lambda x,y: L(x + y) - L(x) - L(y)
## A dual basis for L
(elementary_divisors, dual_basis_pre, _) = L.matrix().smith_form()
elementary_divisors = elementary_divisors.diagonal()
dual_basis = map(operator.div, list(dual_basis_pre), elementary_divisors)
-
+
L_level = ZZ(lcm([ b.denominator() for b in dual_basis ]))
-
+
(elementary_divisors, _, discriminant_basis_pre) = (L_level * matrix(dual_basis)).change_ring(ZZ).smith_form()
elementary_divisors = filter( lambda d: d not in ZZ, (elementary_divisors / L_level).diagonal() )
- elementary_divisors_inv = map(ZZ, [ed**-1 for ed in elementary_divisors])
+ elementary_divisors_inv = map(ZZ, [ed ** -1 for ed in elementary_divisors])
discriminant_basis = matrix(map( operator.mul,
discriminant_basis_pre.inverse().rows()[:len(elementary_divisors)],
elementary_divisors )).transpose()
## This is a form over QQ, so that we cannot use an instance of QuadraticForm
discriminant_form = discriminant_basis.transpose() * L.matrix() * discriminant_basis
- if prod(elementary_divisors_inv) > 100 :
+ if prod(elementary_divisors_inv) > 100:
disc_den = discriminant_form.denominator()
disc_bilinear_pre = \
cython_lambda( ', '.join( ['int a{0}'.format(i) for i in range(discriminant_form.nrows())]
@@ -253,7 +251,7 @@ def _discriminant_form(L):
for i in range(discriminant_form.nrows())
for j in range(discriminant_form.nrows())) )
disc_bilinear = lambda *a: disc_bilinear_pre(*a) / disc_den
- else :
+ else:
disc_bilinear = lambda *xy: vector(ZZ, xy[:discriminant_form.nrows()]) * discriminant_form * vector(ZZ, xy[discriminant_form.nrows():])
disc_quadratic = lambda *a: disc_bilinear(*(2 * a)) / 2
@@ -293,26 +291,26 @@ def _discriminant_form_pmone(L, discriminant_form_exponents):
See test_higherrank_dimension:test__discrimant_form_pmone.
"""
## red gives a normal form for elements in the discriminant group
- red = lambda x : map(operator.mod, x, discriminant_form_exponents)
+ red = lambda x: map(operator.mod, x, discriminant_form_exponents)
## singls and pairs are elements of the discriminant group that are, respectively,
## fixed and not fixed by negation.
singls = list()
pairs = list()
- for x in mrange(discriminant_form_exponents) :
+ for x in mrange(discriminant_form_exponents):
y = red(map(operator.neg, x))
- for (e, f) in zip(x, y) :
- if e < f :
+ for (e, f) in zip(x, y):
+ if e < f:
si = -1
break
- elif e > f :
+ elif e > f:
si = 1
break
else:
singls.append(x)
continue
- if si == 1 :
+ if si == 1:
pairs.append(x)
return (singls, pairs)
@@ -323,7 +321,7 @@ def _weil_representation(CC, L, singls, pairs, plus_basis,
disc_bilinear):
r"""
Construct the Weil representation with values in a complex field
- (or intervall field).
+ (or interval field).
INPUT:
@@ -373,35 +371,36 @@ def _weil_representation(CC, L, singls, pairs, plus_basis,
See test_higherrank_dimension:test__weil_representation.
"""
- zeta_order = ZZ(lcm([8, 12] + map(lambda ex: 2*ex, discriminant_form_exponents)))
+ zeta_order = ZZ(lcm([8, 12] + map(lambda ex: 2 * ex,
+ discriminant_form_exponents)))
zeta = CC(exp(2 * pi * I / zeta_order))
- sqrt2 = CC(sqrt(2))
- drt = CC(sqrt(abs(L.det())))
+ sqrt2 = CC(sqrt(2))
+ drt = CC(sqrt(abs(L.det())))
- Tmat = diagonal_matrix(CC, [zeta**(zeta_order*disc_quadratic(*a))
- for a in (singls + pairs if plus_basis else pairs)])
+ Tmat = diagonal_matrix(CC, [zeta ** (zeta_order * disc_quadratic(*a))
+ for a in (singls + pairs if plus_basis
+ else pairs)])
neg = lambda v: map(operator.neg, v)
- if plus_basis :
- Smat = zeta**(zeta_order / 8 * L.dim()) / drt \
- * matrix( CC, [ [zeta**(-zeta_order * disc_bilinear(*(gamma + delta))) for delta in singls]
+ if plus_basis:
+ Smat = zeta ** (zeta_order / 8 * L.dim()) / drt * matrix(CC, [[zeta**(-zeta_order * disc_bilinear(*(gamma + delta))) for delta in singls]
+ [sqrt2 * zeta**(-zeta_order * disc_bilinear(*(gamma + delta))) for delta in pairs]
- for gamma in singls] \
+ for gamma in singls]
+ [ [sqrt2 * zeta**(-zeta_order * disc_bilinear(*(gamma + delta))) for delta in singls]
+ [zeta**(-zeta_order * disc_bilinear(*(gamma + delta))) + zeta**(-zeta_order * disc_bilinear(*(gamma + neg(delta)))) for delta in pairs]
for gamma in pairs] )
- else :
- Smat = zeta**(zeta_order / 8 * L.dim()) / drt \
- * matrix( CC, [ [ zeta**(-zeta_order * disc_bilinear(*(gamma + delta)))
- - zeta**(-zeta_order * disc_bilinear(*(gamma + neg(delta))))
+ else:
+ Smat = zeta ** (zeta_order / 8 * L.dim()) / drt * matrix( CC, [ [ zeta**(-zeta_order * disc_bilinear(*(gamma + delta)))
+ - zeta ** (-zeta_order * disc_bilinear(*(gamma + neg(delta))))
for delta in pairs]
- for gamma in pairs ] )
+ for gamma in pairs ])
return (Smat, Tmat)
-def _qr(mat) :
+
+def _qr(mat):
r"""
Compute the R matrix in QR decomposition using Housholder reflections.
@@ -460,31 +459,32 @@ def _qr(mat) :
n = mat.ncols()
cur_row = 0
- for j in range(0, n) :
- if all( mat[i,j].contains_zero() for i in xrange(cur_row + 1, m) ) :
- if not mat[cur_row,j].contains_zero() :
+ for j in range(0, n):
+ if all(mat[i, j].contains_zero() for i in xrange(cur_row + 1, m)):
+ if not mat[cur_row, j].contains_zero():
cur_row += 1
continue
- s = sum( (abs(mat[i,j]))**2 for i in xrange(cur_row, m) )
- if s.contains_zero() :
- raise ArithmeticError( "Cannot handle sums of elements that are too imprecise" )
-
- p = sqrt(s)
- if (s - p * mat[cur_row,j]).contains_zero() :
- raise ArithmeticError( "Cannot handle sums of elements that are too imprecise" )
- kappa = 1 / (s - p * mat[cur_row,j])
+ s = sum((abs(mat[i, j])) ** 2 for i in xrange(cur_row, m))
+ if s.contains_zero():
+ raise ArithmeticError("Cannot handle sums of elements that are too imprecise")
- mat[cur_row,j] -= p
- for k in range(j + 1, n) :
- y = sum(mat[i,j].conjugate() * mat[i,k] for i in xrange(cur_row, m)) * kappa
+ p = sqrt(s)
+ if (s - p * mat[cur_row, j]).contains_zero():
+ raise ArithmeticError("Cannot handle sums of elements that are too imprecise")
+ kappa = 1 / (s - p * mat[cur_row, j])
+
+ mat[cur_row, j] -= p
+ for k in range(j + 1, n):
+ y = sum(mat[i, j].conjugate() * mat[i, k]
+ for i in xrange(cur_row, m)) * kappa
for i in range(cur_row, m):
- mat[i,k] -= mat[i,j] * y
+ mat[i, k] -= mat[i, j] * y
+
+ mat[cur_row, j] = p
+ for i in range(cur_row + 1, m):
+ mat[i, j] = CC(0)
- mat[cur_row,j] = p
- for i in range(cur_row + 1, m) :
- mat[i,j] = CC(0)
-
cur_row += 1
-
+
return mat
diff --git a/src/sage/modular/jacobi/test_classical.py b/src/sage/modular/jacobi/test_classical.py
index c57cc23..a10e766 100644
--- a/src/sage/modular/jacobi/test_classical.py
+++ b/src/sage/modular/jacobi/test_classical.py
@@ -7,19 +7,19 @@ AUTHOR:
"""
#===============================================================================
-#
+#
# Copyright (C) 2014 Martin Raum
-#
+#
# This program is free software; you can redistribute it and/or
# modify it under the terms of the GNU General Public License
# as published by the Free Software Foundation; either version 3
# of the License, or (at your option) any later version.
-#
-# This program is distributed in the hope that it will be useful,
+#
+# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
# General Public License for more details.
-#
+#
# You should have received a copy of the GNU General Public License
# along with this program; if not, see <http://www.gnu.org/licenses/>.
#
diff --git a/src/sage/modular/jacobi/test_classical_weak.py b/src/sage/modular/jacobi/test_classical_weak.py
index 4f7a728..bcfe5f2 100644
--- a/src/sage/modular/jacobi/test_classical_weak.py
+++ b/src/sage/modular/jacobi/test_classical_weak.py
@@ -8,19 +8,19 @@ AUTHOR:
"""
#===============================================================================
-#
+#
# Copyright (C) 2010-2014 Martin Raum
-#
+#
# This program is free software; you can redistribute it and/or
# modify it under the terms of the GNU General Public License
# as published by the Free Software Foundation; either version 3
# of the License, or (at your option) any later version.
-#
-# This program is distributed in the hope that it will be useful,
+#
+# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
# General Public License for more details.
-#
+#
# You should have received a copy of the GNU General Public License
# along with this program; if not, see <http://www.gnu.org/licenses/>.
#
@@ -29,7 +29,7 @@ AUTHOR:
from sage.all import (PolynomialRing, LaurentPolynomialRing,
PowerSeriesRing,
ZZ, QQ, gcd,
- factorial, gegenbauer,
+ factorial,
vector, span, FreeModule,
matrix,
ModularForms, Gamma1
@@ -47,9 +47,9 @@ def test_classical_weak_jacobi_forms():
Test classical weak Jacobi forms for correctness. See individual tests
for more details.
- .. NOTE:
+ .. NOTE::
- This is a test generator to be used by nosetest.
+ This is a test generator to be used by nosetest.
TESTS::
@@ -60,13 +60,12 @@ def test_classical_weak_jacobi_forms():
prec = 20
for k in [10, 11, 15, 18]:
- for m in range(1, 4):
-
+ for m in range(1, 4):
yield (_test_classical_weak_jacobi_forms__taylor_coefficients,
k, m, prec)
- if k%2 == 0:
- nu_bound = 20
+ if k % 2 == 0:
+ # nu_bound = 20
yield (_test_classical_weak_jacobi_forms__taylor_coefficient_modularity,
10, k, m, prec)
@@ -74,7 +73,7 @@ def test_classical_weak_jacobi_forms():
yield (_test_classical_weak_jacobi_forms__multiplication,
k, m, k_mod, prec)
- for torsion_point in [0, 1/2]:
+ for torsion_point in [0, QQ.one()/2]:
yield (_test_classical_weak_jacobi_forms__torsion_point,
torsion_point, k, m, prec)
@@ -84,9 +83,9 @@ def _test_classical_weak_jacobi_forms__taylor_coefficients(k, m, prec):
coefficients.
INPUT:
-
+
- `k -- An integer.
-
+
- `m` -- A non-negative integer.
- ``prec`` -- A non-negative integer that corresponds to a precision of
@@ -105,7 +104,8 @@ def _test_classical_weak_jacobi_forms__taylor_coefficients(k, m, prec):
_taylor_coefficients(phi, k, m, prec),
_predicted_taylor_coefficients(tcs, prec) ) )
-def _taylor_coefficients(expansion, k, m, prec) :
+
+def _taylor_coefficients(expansion, k, m, prec):
r"""
Normalized Taylor coefficients of a Jacobi form.
@@ -131,16 +131,16 @@ def _taylor_coefficients(expansion, k, m, prec) :
sage: _taylor_coefficients(phi, 4, 1, 2)
[1 + 240*q + O(q^2), q + O(q^2)]
"""
- R = PowerSeriesRing(ZZ, 'q'); q = R.gen(0)
+ R = PowerSeriesRing(ZZ, 'q')
projs = list()
- for pw in (range(0, 2*m + 1, 2) if k % 2 == 0 else range(1, 2*m - 1, 2)):
- proj = dict( (n, 0) for n in range(prec) )
+ for pw in (range(0, 2 * m + 1, 2) if k % 2 == 0 else range(1, 2*m - 1, 2)):
+ proj = {n: 0 for n in range(prec)}
for (n, r) in classical_weak_jacobi_fe_indices(m, prec):
- ((nred, rred), sign) = classical_jacobi_reduce_fe_index((n,r), m)
- try :
- proj[n] += (sign * r)**pw * expansion[(nred, rred)]
- except (KeyError, ValueError) :
+ ((nred, rred), sign) = classical_jacobi_reduce_fe_index((n, r), m)
+ try:
+ proj[n] += (sign * r) ** pw * expansion[(nred, rred)]
+ except (KeyError, ValueError):
pass
projs.append(proj)
@@ -148,12 +148,13 @@ def _taylor_coefficients(expansion, k, m, prec) :
gcd_projs = [gcd(proj.values()) for proj in projs]
gcd_projs = [g if g != 0 else 1 for g in gcd_projs]
projs = [sorted(proj.iteritems()) for proj in projs]
- projs = [ R([c for (_, c) in proj]).add_bigoh(prec) / gcd_proj
- for (proj, gcd_proj) in zip(projs, gcd_projs) ]
+ projs = [R([c for (_, c) in proj]).add_bigoh(prec) / gcd_proj
+ for (proj, gcd_proj) in zip(projs, gcd_projs)]
return projs
-def _predicted_taylor_coefficients(fs, prec) :
+
+def _predicted_taylor_coefficients(fs, prec):
r"""
Given a list of power series, which are the corrected Taylor coefficients
of a Jacobi form, return the normalized, uncorrected ones, assuming that
@@ -176,14 +177,14 @@ def _predicted_taylor_coefficients(fs, prec) :
sage: _predicted_taylor_coefficients([e4,lambda prec: e4(prec)-e4(prec)], 2)
[1 + 240*q + O(q^2), q + O(q^2)]
"""
- R = PowerSeriesRing(ZZ, 'q'); q = R.gen(0)
+ R = PowerSeriesRing(ZZ, 'q')
diff = lambda f: f.derivative().shift(1)
normalize = lambda f: f / gcd(f.list()) if f != 0 else f
taylor_coefficients = list()
- allf = R(0)
- for f in fs :