summaryrefslogtreecommitdiffstats
diff options
context:
space:
mode:
authorFrédéric Chapoton <chapoton@math.univ-lyon1.fr>2015-03-08 21:26:43 +0100
committerFrédéric Chapoton <chapoton@math.univ-lyon1.fr>2015-03-08 21:26:43 +0100
commit7c4edd96afe2c0cbbdcf43d499206db3dd269b1c (patch)
treedac34df3381d69ba11b83e9096068e69d0fc321c
parentUpdated Sage version to 6.6.beta3 (diff)
parentAdopt to new spacing in doctests. (diff)
Merge branch 'u/mraum/ticket/16448' into 6.6.b3
-rw-r--r--src/sage/modular/jacobi/__init__.py1
-rw-r--r--src/sage/modular/jacobi/all.py30
-rw-r--r--src/sage/modular/jacobi/classical.py221
-rw-r--r--src/sage/modular/jacobi/classical_weak.py772
-rw-r--r--src/sage/modular/jacobi/higherrank.py992
-rw-r--r--src/sage/modular/jacobi/higherrank_dimension.py483
-rw-r--r--src/sage/modular/jacobi/test_classical.py117
-rw-r--r--src/sage/modular/jacobi/test_classical_weak.py386
-rw-r--r--src/sage/modular/jacobi/test_higherrank.py516
-rw-r--r--src/sage/modular/jacobi/test_higherrank_dimension.py212
-rw-r--r--src/sage/modular/jacobi/test_vector_valued.py541
-rw-r--r--src/sage/modular/jacobi/vector_valued.py674
12 files changed, 4945 insertions, 0 deletions
diff --git a/src/sage/modular/jacobi/__init__.py b/src/sage/modular/jacobi/__init__.py
new file mode 100644
index 00000000..c9fecac
--- /dev/null
+++ b/src/sage/modular/jacobi/__init__.py
@@ -0,0 +1 @@
+import all
diff --git a/src/sage/modular/jacobi/all.py b/src/sage/modular/jacobi/all.py
new file mode 100644
index 00000000..89bb9d9
--- /dev/null
+++ b/src/sage/modular/jacobi/all.py
@@ -0,0 +1,30 @@
+from classical_weak import (
+ classical_jacobi_reduce_fe_index,
+ classical_weak_jacobi_fe_indices,
+
+ classical_weak_jacobi_forms
+)
+
+from classical import (
+ classical_jacobi_fe_indices,
+
+ classical_jacobi_forms
+)
+
+from higherrank import (
+ higherrank_jacobi_fe_indices,
+ higherrank_jacobi_reduce_fe_index,
+ higherrank_jacobi_r_classes,
+
+ higherrank_jacobi_forms
+)
+
+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
new file mode 100644
index 00000000..636d7c5
--- /dev/null
+++ b/src/sage/modular/jacobi/classical.py
@@ -0,0 +1,221 @@
+r"""
+Fourier expansions of classical Jacobi forms.
+
+AUTHOR:
+
+- Martin Raum
+
+EXAMPLES:
+
+To compute a basis of weight `k` and index `m \in \ZZ` Jacobi forms,
+we call ``classical_jacobi_forms``. This compute weight `10`, index
+`2` Jacobi forms up to precision `5`.
+
+::
+
+ sage: from sage.modular.jacobi.all import classical_jacobi_forms
+ sage: classical_jacobi_forms(10, 2, 5)
+ [{(0, 0): 823680,
+ (1, 0): -132447744,
+ (1, 1): -42172416,
+ (1, 2): -329472,
+ (2, 0): -47908194048,
+ (2, 1): -27622932480,
+ (2, 2): -4157936640,
+ (3, 0): -1497954991104,
+ (3, 1): -1045412020224,
+ (3, 2): -318012963840,
+ (4, 0): -17346153217536,
+ (4, 1): -13216244760576,
+ (4, 2): -5574898188288},
+ {(1, 0): -22464,
+ (1, 1): 9984,
+ (1, 2): 1248,
+ (2, 0): -19968,
+ (2, 1): 149760,
+ (2, 2): -149760,
+ (3, 0): 913536,
+ (3, 1): -1707264,
+ (3, 2): 1123200,
+ (4, 0): 359424,
+ (4, 1): 4253184,
+ (4, 2): -2715648}]
+
+For description of Fourier expansions, see the documentation of
+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,
+# 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.combinat.dict_addition import dict_linear_combination
+from sage.matrix.all import matrix
+from sage.misc.all import isqrt, cached_function
+from sage.modular.jacobi.classical_weak import classical_weak_jacobi_fe_indices, classical_weak_jacobi_forms
+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
+
+def classical_jacobi_fe_indices(m, prec, reduced=False):
+ r"""
+ Indices `(n,r)` of Fourier expansions of Jacobi forms of index `m`.
+
+ INPUT:
+
+ - `m` -- A positive integer.
+
+ - ``prec`` -- A non-negative integer.
+
+ - ``reduce`` -- A boolean (default: ``False``). If ``True``
+ restrict to `0 \le r \le m`.
+
+ OUTPUT:
+
+ A generator of pairs of integers `(n,r)`.
+
+ EXAMPLES::
+
+ sage: from sage.modular.jacobi.all import *
+ sage: list(classical_jacobi_fe_indices(2, 3, True))
+ [(1, 0), (1, 1), (1, 2), (2, 0), (2, 1), (2, 2), (0, 0)]
+ sage: list(classical_jacobi_fe_indices(2, 3, False))
+ [(1, 0), (1, 1), (1, -1), (1, 2), (1, -2), (2, 0), (2, 1), (2, -1), (2, 2), (2, -2), (2, 3), (2, -3), (0, 0), (2, 4), (2, -4)]
+ """
+ fm = Integer(4*m)
+
+ if reduced :
+ # positive definite forms
+ for n in range(1, prec):
+ for r in range(min(m + 1, isqrt(fm * n - 1) + 1)):
+ yield (n, r)
+
+ # 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 :
+ # positive definite forms
+ for n in range(1, prec):
+ yield(n, 0)
+ for r in range(1, isqrt(fm * n - 1) + 1):
+ yield (n, r)
+ yield (n, -r)
+
+ # indefinite forms
+ if prec > 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.
+
+ 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)
+ [
+ (1, 0, 0),
+ (0, 0, 1)
+ ]
+ sage: _classical_jacobi_forms_as_weak_jacobi_forms(12, 1)
+ [
+ (1, 0, 0),
+ (0, 1, 0)
+ ]
+ """
+ prec = m//4 + 1
+ weak_forms = classical_weak_jacobi_forms(k, m, prec, algorithm)
+ indices = list(classical_weak_jacobi_fe_indices(m, prec, reduced=True))
+ 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"):
+ r"""
+ Compute Fourier expansions of classical Jacobi forms of weight `k`
+ and index `m`.
+
+ INPUT:
+
+ - `k` -- An integer.
+
+ - `m` -- A non-negative integer.
+
+ - ``prec`` -- A non-negative integer that corresponds to a
+ precision of the `q`-expansion.
+
+ - ``algorithm`` -- Default: ''skoruppa''. Only ''skoruppa'' is
+ implemented.
+
+ OUTPUT:
+
+ A list of dictionaries, mapping indices `(n,r)` to rationals.
+
+ EXAMPLES:
+
+ sage: from sage.modular.jacobi.classical import *
+ sage: k = 4; m = 2; prec = 5
+ sage: classical_jacobi_forms(k, m, prec)
+ [{(0, 0): 40320,
+ (1, 0): 3386880,
+ (1, 1): 2580480,
+ (1, 2): 564480,
+ (2, 0): 23143680,
+ (2, 1): 18063360,
+ (2, 2): 11289600,
+ (3, 0): 51932160,
+ (3, 1): 54190080,
+ (3, 2): 33868800,
+ (4, 0): 138862080,
+ (4, 1): 108380160,
+ (4, 2): 95477760}]
+
+ TESTS:
+
+ See ``test_classical.py``.
+ """
+ weak_forms = classical_weak_jacobi_forms(k, m, prec, algorithm)
+ coords = _classical_jacobi_forms_as_weak_jacobi_forms(k, m)
+
+ return [dict_linear_combination(zip(weak_forms, cs)) for cs in coords]
diff --git a/src/sage/modular/jacobi/classical_weak.py b/src/sage/modular/jacobi/classical_weak.py
new file mode 100644
index 00000000..70217f3
--- /dev/null
+++ b/src/sage/modular/jacobi/classical_weak.py
@@ -0,0 +1,772 @@
+r"""
+Fourier expansions of classical weak Jacobi forms.
+
+AUTHOR:
+
+- Martin Raum
+
+REFERENCES:
+
+.. [Sk84] Skoruppa, Uber den Zusammenhang zwischen Jacobiformen und
+ Modulformen halbganzen Gewichts, 1984, University of Bonn.
+
+EXAMPLES:
+
+To compute a basis of weight `k` and index `m \in \ZZ` weak Jacobi
+forms, we call ``classical_weak_jacobi_forms``. This compute weight
+`9`, index `2` Jacobi forms up to precision `5`.
+
+::
+
+ sage: from sage.modular.jacobi.all import *
+ sage: jforms = classical_weak_jacobi_forms(9, 2, 5)
+ sage: jforms
+ [{(0, 1): 660,
+ (1, 1): -172260,
+ (2, 1): -89901900,
+ (3, 1): -3699449160,
+ (4, 1): -56862841860}]
+
+
+Fourier expansions are represented as dictionaries. Fourier
+coefficients of Jacobi forms are index by pairs `(n,r)` of integers.
+We call `(n,r)` reduced if `0 \le r \le m`. Fourier expansions are
+given in terms of Fourier coefficients of reduced index. They can be
+accessed directly; If a reduced pair does not occur as a key and `n`
+does not exceed the prescribed precision, then the corresponding
+Fourier coefficient vanishes.
+
+To access Fourier coefficients of non reduced index, we compute the attached reduction by ``classical_jacobi_reduce_fe_index``::
+
+ sage: classical_jacobi_reduce_fe_index((2, 3), 2)
+ ((1, 1), -1)
+
+As a result, we obtain `(n', r')` and a sign `s`. The Fourier
+coefficient at `(n,r)` equals `s^k` times the coefficient at `(n',
+r')`. In the specific examples (`k` is odd), the `(2,3)`-th Fourier
+coefficient of the first basis element is equal to::
+
+ sage: -jforms[0][(1,1)]
+ 172260
+"""
+
+#===============================================================================
+#
+# 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,
+# 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.combinat.all import number_of_partitions
+from sage.functions.all import binomial, factorial
+from sage.matrix.all import matrix
+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 import big_oh
+import operator
+
+
+def classical_jacobi_reduce_fe_index((n, r), m):
+ r"""
+ Reduce an index `(n, r)` of the Fourier expansion of Jacobi forms
+ of index `m` suchthat `0 \le r \le m`.
+
+ INPUT:
+
+ - `(n, r)` -- A pair of integers.
+
+ - `m` -- A positive integer.
+
+ OUTPUT:
+
+ A pair `((n', r'), s)`, where `(n', r')` is an index that is
+ equivalent to `(n,r)` and `s` is a sign indicating whether `r -
+ r'` or `r + r'` is divisible by `2m`.
+
+ EXAMPLES::
+
+ sage: from sage.modular.jacobi.all import *
+ sage: classical_jacobi_reduce_fe_index((2,2), 1)
+ ((1, 0), 1)
+ """
+ rred = r % (2*m)
+
+ if rred > m:
+ sgn = -1
+ 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):
+ r"""
+ Indices `(n,r)` of Fourier expansions of weak Jacobi forms of index `m`.
+
+ INPUT:
+
+ - `m` -- A positive integer.
+
+ - ``prec`` -- A non-negative integer.
+
+ - ``reduce`` -- A boolean (default: ``False``). If ``True``
+ restrict to `0 \le r \le m`.
+
+ OUTPUT:
+
+ A generator of pairs `(n,r)`.
+
+ EXAMPLES::
+
+ sage: from sage.modular.jacobi.all import *
+ sage: list(classical_weak_jacobi_fe_indices(2, 3, True))
+ [(1, 0), (1, 1), (1, 2), (2, 0), (2, 1), (2, 2), (0, 0), (0, 1), (0, 2)]
+ sage: list(classical_weak_jacobi_fe_indices(2, 3, False))
+ [(1, 0), (1, 1), (1, -1), (1, 2), (1, -2), (2, 0), (2, 1), (2, -1), (2, 2), (2, -2), (2, 3), (2, -3), (0, -1), (1, 3), (2, -4), (0, 0), (2, 4), (1, -3), (0, 1), (0, -2), (0, 2)]
+ """
+ fm = Integer(4*m)
+
+ if reduced :
+ # positive definite forms
+ for n in range(1, prec):
+ 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 ) :
+ yield (n, r)
+ else :
+ # positive definite forms
+ for n in range(1, prec):
+ yield(n, 0)
+ for r in range(1, isqrt(fm * n - 1) + 1):
+ yield (n, r)
+ yield (n, -r)
+
+ # indefinite forms
+ 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)
+
+ raise StopIteration
+
+## TODO: Implement caching by truncation.
+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
+ the q-expansion.
+
+ - ``algorithm`` -- Default: ''skoruppa''. Only ''skoruppa'' is implemented.
+
+ EXAMPLES::
+
+ sage: from sage.modular.jacobi.all import *
+ sage: classical_weak_jacobi_forms(2, 1, 5)
+ [{(0, 0): -4,
+ (0, 1): 2,
+ (1, 0): -984,
+ (1, 1): 496,
+ (2, 0): -14512,
+ (2, 1): 8238,
+ (3, 0): -106016,
+ (3, 1): 67024,
+ (4, 0): -574488,
+ (4, 1): 385026}]
+
+ TESTS:
+
+ See ``test_classical_weak.py``.
+ """
+ if algorithm != "skoruppa":
+ raise NotImplementedError("Algorithm {} is not implemented.".format(algorithm))
+ factory = ClassicalWeakJacobiFormsFactory(m, prec)
+
+ 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) :
+ r"""
+ 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
+ ``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 :
+ nmb_modular_forms = m + 1
+ start_weight = k
+ 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)) :
+ 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
+
+_classical_weak_jacobi_forms_factory_cache = WeakValueDictionary()
+
+def ClassicalWeakJacobiFormsFactory(m, prec):
+ r"""
+ INPUT:
+
+ - `m` -- A non-negative integer.
+
+ - ``prec`` -- An integer.
+
+ EXAMPLES::
+
+ sage: from sage.modular.jacobi.classical_weak import ClassicalWeakJacobiFormsFactory
+ sage: fcty = ClassicalWeakJacobiFormsFactory(2, 5)
+
+ TESTS::
+
+ sage: fcty = ClassicalWeakJacobiFormsFactory(2, 5)
+ sage: assert fcty is ClassicalWeakJacobiFormsFactory(2, 5)
+ """
+ if (m, prec) in _classical_weak_jacobi_forms_factory_cache:
+ return _classical_weak_jacobi_forms_factory_cache[(m, prec)]
+
+ for (mc, pc) in _classical_weak_jacobi_forms_factory_cache:
+ if m == mc and pc > prec:
+ return _classical_weak_jacobi_forms_factory_cache[(m, pc)]._truncate_cache(prec)
+
+ fcty = ClassicalWeakJacobiForms_factory(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) :
+ r"""
+ INPUT:
+
+ - `m` -- A non-negative integer.
+
+ - ``prec`` -- An integer.
+
+ EXAMPLES::
+
+ sage: from sage.modular.jacobi.classical_weak import ClassicalWeakJacobiForms_factory
+ sage: fcty = ClassicalWeakJacobiForms_factory(2, 5)
+ """
+ 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) :
+ r"""
+ The index of the Jacobi forms that are computed by this factory.
+
+ OUTPUT:
+
+ An integer.
+
+ EXAMPLES::
+
+ sage: from sage.modular.jacobi.classical_weak import ClassicalWeakJacobiForms_factory
+ sage: fcty = ClassicalWeakJacobiForms_factory(2, 5)
+ sage: fcty.jacobi_index()
+ 2
+ """
+ return self.__jacobi_index
+
+ def precision(self) :
+ r"""
+ The precision of Fourier expansions that are computed.
+
+ OUTPUT:
+
+ An integer.
+
+ EXAMPLES::
+
+ sage: from sage.modular.jacobi.classical_weak import ClassicalWeakJacobiForms_factory
+ sage: fcty = ClassicalWeakJacobiForms_factory(2, 5)
+ sage: fcty.precision()
+ 5
+ """
+ return self.__prec
+
+ def _power_series_ring(self) :
+ r"""
+ An auxiliary power series ring that is cached in the factory.
+
+ OUTPUT:
+
+ An power series ring.
+
+ EXAMPLES::
+
+ sage: from sage.modular.jacobi.classical_weak import ClassicalWeakJacobiForms_factory
+ sage: fcty = ClassicalWeakJacobiForms_factory(2, 5)
+ sage: fcty._power_series_ring()
+ Power Series Ring in q over Rational Field
+ """
+ return self.__power_series_ring
+
+ def _power_series_ring_ZZ(self) :
+ r"""
+ An auxiliary power series ring over `\ZZ` that is cached in the factory.
+
+ OUTPUT:
+
+ An power series ring.
+
+ EXAMPLES::
+
+ sage: from sage.modular.jacobi.classical_weak import ClassicalWeakJacobiForms_factory
+ sage: fcty = ClassicalWeakJacobiForms_factory(2, 5)
+ sage: fcty._power_series_ring_ZZ()
+ Power Series Ring in q over Integer Ring
+ """
+ return self.__power_series_ring_ZZ
+
+ def _truncate_cache(self, prec):
+ r"""
+ A new factory instance whose cache is truncated to a given precision.
+
+ INPUT:
+
+ - ``prec`` -- An integer.
+
+ OUTPUT:
+
+ A factory for classical weak Jacobi forms.
+
+ EXAMPLES::
+
+ sage: from sage.modular.jacobi.classical_weak import ClassicalWeakJacobiForms_factory
+ sage: fcty = ClassicalWeakJacobiForms_factory(2, 5)
+ sage: fcty_tr = fcty._truncate_cache(3)
+ sage: assert isinstance(fcty_tr, ClassicalWeakJacobiForms_factory)
+ """
+ other = ClassicalWeakJacobiFormsFactory(self.jacobi_index(), prec)
+
+ try:
+ other._wronskian_adjoint_even = \
+ matrix([[e.truncate(prec) for e in r]
+ for r in self._wronskian_adjoint_even.rows()])
+ except AttributeError:
+ pass
+
+ try:
+ other._wronskian_adjoint_odd = \
+ matrix([[e.truncate(prec) for e in r]
+ for r in self._wronskian_adjoint_odd.rows()])
+ except AttributeError:
+ pass
+
+ try:
+ other.__wronskian_invdeterminant_even = \
+ self.__wronskian_invdeterminant_even.truncate(prec)
+ except AttributeError:
+ pass
+
+ try:
+ other.__wronskian_invdeterminant_odd = \
+ self.__wronskian_invdeterminant_odd.truncate(prec)
+ except AttributeError:
+ pass
+
+ return other
+
+ 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 `\Z`.
+
+ - ``weight_parity`` -- An integer (default: `0`).
+
+ TESTS::
+
+ sage: from sage.modular.jacobi.classical_weak import ClassicalWeakJacobiForms_factory
+ sage: fcty = ClassicalWeakJacobiForms_factory(2, 5)
+ sage: R = fcty._power_series_ring_ZZ()
+ sage: adj = [[R(1), R(0)], [R(7), R(-2)]]
+ sage: fcty._set_wronskian_adjoint(adj, 0)
+ sage: assert all(all(e is f for (e,f) in zip(*rr)) for rr in zip(fcty._wronskian_adjoint(0),adj))
+
+ """
+ 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 :
+ self._wronskian_adjoint_even = wronskian_adjoint
+ else :
+ self._wronskian_adjoint_odd = wronskian_adjoint
+
+ 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:
+
+ A matrix over a power series ring.
+
+ EXAMPLES::
+
+ sage: from sage.modular.jacobi.classical_weak import ClassicalWeakJacobiForms_factory
+ sage: fcty = ClassicalWeakJacobiForms_factory(2, 5)
+ sage: fcty._wronskian_adjoint(0)
+ [[48 - 720*q - 8400*q^3 + 5040*q^4 + O(q^5),
+ -60 + 260*q + 2436*q^3 - 5180*q^4 + O(q^5),
+ 12 - 20*q - 84*q^3 + 140*q^4 + O(q^5)],
+ [3072*q^2 + O(q^5),
+ 32 - 960*q^2 + 2592*q^4 + O(q^5),
+ -8 + 48*q^2 - 72*q^4 + O(q^5)],
+ [-960*q^2 - 4032*q^3 + O(q^5),
+ -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 :
+ return self._wronskian_adjoint_even
+ else :
+ return self._wronskian_adjoint_ood
+
+ 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 :
+ 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)) :
+ 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() )
+
+ 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]
+ 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]
+ 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]
+ 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]
+ 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]
+
+ Wadj = matrix(PS, [ [adj00, adj01, adj02],
+ [adj10, adj11, adj12],
+ [adj20, adj21, adj22] ])
+
+ 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]
+ adj03 = -W[0,3]*W[1,2]*W[2,1] + W[0,2]*W[1,3]*W[2,1] + W[0,3]*W[1,1]*W[2,2] - W[0,1]*W[1,3]*W[2,2] - W[0,2]*W[1,1]*W[2,3] + W[0,1]*W[1,2]*W[2,3]
+
+ adj10 = -W[0,2]*W[1,1]*W[3,0] + W[0,1]*W[1,2]*W[3,0] + W[0,2]*W[1,0]*W[3,1] - W[0,0]*W[1,2]*W[3,1] - W[0,1]*W[1,0]*W[3,2] + W[0,0]*W[1,1]*W[3,2]
+ adj11 = -W[0,3]*W[1,1]*W[3,0] + W[0,1]*W[1,3]*W[3,0] + W[0,3]*W[1,0]*W[3,1] - W[0,0]*W[1,3]*W[3,1] - W[0,1]*W[1,0]*W[3,3] + W[0,0]*W[1,1]*W[3,3]
+ adj12 = -W[0,3]*W[1,2]*W[3,0] + W[0,2]*W[1,3]*W[3,0] + W[0,3]*W[1,0]*W[3,2] - W[0,0]*W[1,3]*W[3,2] - W[0,2]*W[1,0]*W[3,3] + W[0,0]*W[1,2]*W[3,3]
+ adj13 = -W[0,3]*W[1,2]*W[3,1] + W[0,2]*W[1,3]*W[3,1] + W[0,3]*W[1,1]*W[3,2] - W[0,1]*W[1,3]*W[3,2] - W[0,2]*W[1,1]*W[3,3] + W[0,1]*W[1,2]*W[3,3]
+
+ adj20 = -W[0,2]*W[2,1]*W[3,0] + W[0,1]*W[2,2]*W[3,0] + W[0,2]*W[2,0]*W[3,1] - W[0,0]*W[2,2]*W[3,1] - W[0,1]*W[2,0]*W[3,2] + W[0,0]*W[2,1]*W[3,2]
+ adj21 = -W[0,3]*W[2,1]*W[3,0] + W[0,1]*W[2,3]*W[3,0] + W[0,3]*W[2,0]*W[3,1] - W[0,0]*W[2,3]*W[3,1] - W[0,1]*W[2,0]*W[3,3] + W[0,0]*W[2,1]*W[3,3]
+ adj22 = -W[0,3]*W[2,2]*W[3,0] + W[0,2]*W[2,3]*W[3,0] + W[0,3]*W[2,0]*W[3,2] - W[0,0]*W[2,3]*W[3,2] - W[0,2]*W[2,0]*W[3,3] + W[0,0]*W[2,2]*W[3,3]
+ adj23 = -W[0,3]*W[2,2]*W[3,1] + W[0,2]*W[2,3]*W[3,1] + W[0,3]*W[2,1]*W[3,2] - W[0,1]*W[2,3]*W[3,2] - W[0,2]*W[2,1]*W[3,3] + W[0,1]*W[2,2]*W[3,3]
+
+ adj30 = -W[1,2]*W[2,1]*W[3,0] + W[1,1]*W[2,2]*W[3,0] + W[1,2]*W[2,0]*W[3,1] - W[1,0]*W[2,2]*W[3,1] - W[1,1]*W[2,0]*W[3,2] + W[1,0]*W[2,1]*W[3,2]
+ adj31 = -W[1,3]*W[2,1]*W[3,0] + W[1,1]*W[2,3]*W[3,0] + W[1,3]*W[2,0]*W[3,1] - W[1,0]*W[2,3]*W[3,1] - W[1,1]*W[2,0]*W[3,3] + W[1,0]*W[2,1]*W[3,3]
+ adj32 = -W[1,3]*W[2,2]*W[3,0] + W[1,2]*W[2,3]*W[3,0] + W[1,3]*W[2,0]*W[3,2] - W[1,0]*W[2,3]*W[3,2] - W[1,2]*W[2,0]*W[3,3] + W[1,0]*W[2,2]*W[3,3]
+ adj33 = -W[1,3]*W[2,2]*W[3,1] + W[1,2]*W[2,3]*W[3,1] + W[1,3]*W[2,1]*W[3,2] - W[1,1]*W[2,3]*W[3,2] - W[1,2]*W[2,1]*W[3,3] + W[1,1]*W[2,2]*W[3,3]
+
+ Wadj = matrix(PS, [ [adj00, adj01, adj02, adj03],
+ [adj10, adj11, adj12, adj13],
+ [adj20, adj21, adj22, adj23],
+ [adj30, adj31, adj32, adj33] ])
+ 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 :
+ self._wronskian_adjoint_even = wronskian_adjoint
+ else :
+ self._wronskian_adjoint_odd = wronskian_adjoint
+
+ return wronskian_adjoint
+
+ def _set_wronskian_invdeterminant(self, wronskian_invdeterminant, weight_parity = 0) :
+ r"""
+ 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::
+
+ sage: from sage.modular.jacobi.classical_weak import ClassicalWeakJacobiForms_factory
+ sage: fcty = ClassicalWeakJacobiForms_factory(2, 5)
+ sage: invdet = fcty._power_series_ring_ZZ()(13)
+ 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 :
+ wronskian_invdeterminant = self.__power_series_ring_ZZ(wronskian_invdeterminant)
+
+ if weight_parity % 2 == 0 :
+ self._wronskian_invdeterminant_even = wronskian_invdeterminant
+ else :
+ self._wronskian_invdeterminant_odd = wronskian_invdeterminant
+
+ 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.
+
+ INPUT:
+
+ - ``weight_parity`` -- An integer (default: `0`).
+
+ OUTPUT:
+
+ A power series.
+
+ EXAMPLES::
+
+ sage: from sage.modular.jacobi.classical_weak import ClassicalWeakJacobiForms_factory
+ sage: fcty = ClassicalWeakJacobiForms_factory(2, 5)
+ 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 :
+ return self._wronskian_invdeterminant_even
+ else :
+ return self._wronskian_invdeterminant_odd
+
+ except AttributeError :
+ m = self.jacobi_index()
+ if weight_parity % 2 == 0 :
+ pw = (m + 1) * (2 * m + 1)
+ 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 :
+ self._wronskian_invdeterminant_even = wronskian_invdeterminant
+ else :
+ self._wronskian_invdeterminant_odd = wronskian_invdeterminant
+
+ return wronskian_invdeterminant
+
+ 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.
+
+ ALGORITHM:
+
+ We combine the theta decomposition and the heat operator as in
+ [Sk84]. This yields a bijections of the space of weak Jacobi
+ 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: 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)
+ sage: exp_gcd = gcd(expansion.values())
+ 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 :
+ PS = self.__power_series_ring_ZZ
+ 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
+ return dict()
+
+ f_divs = dict()
+ 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) :
+ 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 :
+ ## 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.
+ 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))
+ / factorial(i + k + j - 1)
+ * factorial(2*self.jacobi_index() + k - 1) )
+ for j in range(i + 1) ) )
+ 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)
+ * prod(2*l + 1 for l in xrange(i + 1))
+ / factorial(i + k + j)
+ * 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) )
+ series = self._wronskian_invdeterminant(k) * series
+
+ 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
new file mode 100644
index 00000000..4849c4b
--- /dev/null
+++ b/src/sage/modular/jacobi/higherrank.py
@@ -0,0 +1,992 @@
+r"""
+Using restriction to scalar indices, we compute Jacobi forms of arbitrary index.
+
+AUTHOR:
+
+- Martin Raum
+
+REFERENCE:
+
+.. [Ra] Martin Raum, Computation of Jacobi forms degree 1 and higher rank index.
+
+EXAMPLES:
+
+To compute a basis of weight `k` and index `m` Jacobi forms, we call
+``higherrank_jacobi_forms``. The index `m` is a quadratic form over
+`\ZZ`. This compute weight `10` Jacobi forms up to precision `20`.
+
+::
+
+ sage: from sage.modular.jacobi.all import *
+ sage: m = QuadraticForm(ZZ, 2, [1,1,1])
+ sage: jforms = higherrank_jacobi_forms(10, m, 5)
+ sage: jforms
+ [{(0, (0, 0)): 1,
+ (1, (0, 0)): 0,
+ (1, (1, 1)): -45,
+ (2, (0, 0)): -59130,
+ (2, (1, 1)): -12672,
+ (3, (0, 0)): -1416960,
+ (3, (1, 1)): -558045,
+ (4, (0, 0)): -14346720,
+ (4, (1, 1)): -7165440},
+ {(0, (0, 0)): 0,
+ (1, (0, 0)): 1,
+ (1, (1, 1)): -1/6,
+ (2, (0, 0)): -15,
+ (2, (1, 1)): 5/3,
+ (3, (0, 0)): 90,
+ (3, (1, 1)): -4/3,
+ (4, (0, 0)): -248,
+ (4, (1, 1)): -155/3}]
+
+
+Fourier expansions are represented as dictionaries. They are similar to the ones used for classical Jacobi forms. Fourier
+coefficients of Jacobi forms are index by pairs `(n,r)`, where `n` is an integer and `r` is a tuple. The notion of reduced `r` depends on a choice and is dictated by the following output::
+
+ sage: (r_classes, _) = higherrank_jacobi_r_classes(m)
+ sage: r_classes
+ [[(0, 0)],
+ [(1, 1), (-1, -1), (0, 1), (0, -1), (1, 0), (-1, 0)]]
+
+This is a list of lists, the first elements of which are, by convention, reduced. They are guarenteed to have minimal length in their class in `\ZZ^l / m \ZZ^l`, where `l` is the rank of `m` and `m` is identified with its Gram matrix::
+
+ sage: ZZ^m.dim() / m.matrix().row_module()
+ Finitely generated module V/W over Integer Ring with invariants (3)
+
+Fourier expansions are given in terms of Fourier coefficients of
+reduced index. They can be accessed directly; If a reduced pair does
+not occur as a key and `n` does not exceed the prescribed precision,
+then the corresponding Fourier coefficient vanishes.
+
+To access Fourier coefficients of non reduced index, we compute the attached reduction by ``classical_jacobi_reduce_fe_index``::
+
+ sage: m_adj = QuadraticForm(2 * m.matrix().adjoint())
+ sage: m_span = m.matrix().row_module()
+ sage: higherrank_jacobi_reduce_fe_index((4, (2,1)), m, r_classes, m_adj, m_span)
+ ((3, (0, 0)), 1)
+
+As a result, we obtain `(n', r')` and a sign `s`. The Fourier
+coefficient at `(n,r)` equals `s^k` times the coefficient at `(n',
+r')`. In the specific examples (`k` is odd), the `(4, (2,1))`-th Fourier
+coefficient of the first basis element is equal to::
+
+ sage: jforms[0][(3,(0,0))]
+ -1416960
+"""
+
+#===============================================================================
+#
+# 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,
+# 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.matrix.all import matrix, zero_matrix, identity_matrix
+from sage.misc.cython import cython_lambda
+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.rings.all import ZZ, QQ
+from sage.quadratic_forms.all import QuadraticForm
+from sage.sets.all import Set
+
+import operator
+from random import Random
+
+
+def higherrank_jacobi_reduce_fe_index((n, r), m, r_classes, m_adj, m_span):
+ r"""
+ Reduce a Fourier index `(n, r)`.
+
+ INPUT:
+
+ - `(n, r)` -- A pair of an integer and a tuple of integers.
+
+ - `m` -- A quadratic form over `\Z`.
+
+ - ``r_classes`` -- A list of lists of vectors.
+
+ - `m_adj` -- A quadratic form over `\Z`.
+
+ - ``m_span`` -- The row (or column) span `m`.
+
+ OUTPUT:
+
+ - A pair `((n', r'), s)` where `(n', r')` is the reduced index and
+ `s = \pm 1` tells whether r or -r is equivalent to r modulo `m
+ \Z^l`.
+
+ EXAMPLES::
+
+ sage: from sage.modular.jacobi.higherrank import higherrank_jacobi_r_classes
+ sage: from sage.modular.jacobi.higherrank import higherrank_jacobi_reduce_fe_index
+ sage: m = QuadraticForm(matrix([[2]]))
+ sage: m_adj = QuadraticForm(2 * m.matrix().adjoint())
+ sage: m_span = m.matrix().row_module()
+ sage: r_classes = higherrank_jacobi_r_classes(m)[0]
+ sage: higherrank_jacobi_reduce_fe_index((1,(2,)), m, r_classes, m_adj, m_span)
+ ((0, (0,)), 1)
+
+ TESTS:
+
+ See ``test_higherrank.py``.
+ """
+ try:
+ (rred, sgn) = _higherrank_jacobi_reduce_fe_index__r__cache[(m,r)]
+ except KeyError:
+ (rred, sgn) = _higherrank_jacobi_reduce_fe_index__r(r, r_classes, m_span)
+ _higherrank_jacobi_reduce_fe_index__r__cache[(m,r)] = (rred, sgn)
+
+ nred = n - (m_adj(r) - m_adj(rred)) // (2*m.det())
+
+ return ((nred, rred), sgn)
+
+_higherrank_jacobi_reduce_fe_index__r__cache = {}
+def _higherrank_jacobi_reduce_fe_index__r(r, r_classes, m_span):
+ r"""
+ Find a representative in `r_classes` that is equivalent modulo `m
+ \Z^l` and `\pm` to `r`.
+
+ INPUT:
+
+ - `r` -- A tuple of integers.
+
+ - ``r_classes`` -- A list of lists of vectors.
+
+ - ``m_span`` -- The row (or column) span `m`.
+
+ OUTPUT:
+
+ - A pair `(r', s)`, where `r'` is reduced and `s = \pm 1` tells
+ whether r or -r is equivalent to r modulo `m \Z^l`.
+
+ EXAMPLES::
+
+ sage: from sage.modular.jacobi.higherrank import higherrank_jacobi_r_classes
+ sage: from sage.modular.jacobi.higherrank import _higherrank_jacobi_reduce_fe_index__r
+ sage: m = QuadraticForm(matrix([[2]]))
+ sage: m_span = m.matrix().row_module()
+ sage: r_classes = higherrank_jacobi_r_classes(m)[0]
+ sage: _higherrank_jacobi_reduce_fe_index__r((2,), r_classes, m_span)
+ ((0,), 1)
+
+ TESTS:
+
+ See ``test_higherrank.py``.
+ """
+ for r_class in r_classes:
+ rred = r_class[0]
+
+ if vector(r) - vector(rred) in m_span:
+ return (rred, 1)
+ if vector(r) + vector(rred) in m_span:
+ return (rred, -1)
+ else :
+ raise RuntimeError( "Could not find reduced r" )
+
+def higherrank_jacobi_fe_indices(m, prec, r_classes, reduced=False):
+ r"""
+ Indices `(n, r)` of Fourier expansions of Jacobi forms of index
+ `m`, where `n` is an integer and `r` is a tuple.
+
+ INPUT:
+
+ - `m` -- A quadratic form over `\Z`.
+
+ - ``prec`` -- A nonnegative integer.
+
+ - ``r_classes`` -- A list of list of integers.
+
+ - ``reduced`` -- Boolean. Default: ``False``.
+
+ OUTPUT:
+
+ - A generator of pairs `(n, r)`, where `n` is an integer and `r` is a tupel.
+
+ EXAMPLES::
+
+ sage: from sage.modular.jacobi.higherrank import higherrank_jacobi_r_classes
+ sage: from sage.modular.jacobi.higherrank import higherrank_jacobi_fe_indices
+ sage: m = QuadraticForm(matrix([[2]]))
+ sage: r_classes = higherrank_jacobi_r_classes(m)[0]
+ sage: list(higherrank_jacobi_fe_indices(m, 2, r_classes, reduced=True))
+ [(0, (0,)), (1, (0,)), (1, (1,))]
+ sage: list(higherrank_jacobi_fe_indices(m, 2, r_classes, reduced=False))
+ [(0, (0,)), (1, (0,)), (1, (1,)), (1, (-1,)), (1, (2,)), (1, (-2,))]
+
+ TESTS::
+
+ sage: higherrank_jacobi_fe_indices(m, 2, r_classes, reduced=True)
+ <generator object ...
+
+ See also ``test_higherrank.py``.
+ """
+ m_adj = QuadraticForm(2 * m.matrix().adjoint())
+
+ if reduced:
+ for n in range(0,prec):
+ for r_class in r_classes:
+ r = r_class[0]
+ if m_adj(r) <= 2*m.det()*n:
+ yield (n, r)
+ 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] :
+ yield (n, tuple(r))
+
+ raise StopIteration
+
+def higherrank_jacobi_r_classes(m):
+ r"""
+ Let `l` be the dimension of `m`. For each element of `(\Z^l / m
+ \Z^l) \pm m` of minimal norm, we compute all representatives that
+ minimize the norm with respect to the adjoint of `m`.
+
+ INPUT:
+
+ - `m` -- A quadratic form over `\Z`.
+
+ OUTPUT:
+
+ - A list of lists of vectors.
+
+ EXAMPLES::
+
+ sage: from sage.modular.jacobi.higherrank import higherrank_jacobi_r_classes
+ sage: m = QuadraticForm(matrix([[2]]))
+ sage: higherrank_jacobi_r_classes(m)
+ ([[(0,)], [(1,), (-1,)]], [[1], [1, 1]])
+
+ TESTS:
+
+ See ``test_higherrank.py``.
+ """
+ m_mat = m.matrix()
+ 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]
+ max_norm = max(m_adj(r) for r in canonical_reps)
+
+ current_max_length = 5
+ short_vectors = m_adj.short_vector_list_up_to_length(current_max_length)
+
+ r_classes = []
+ r_classes_reduction_signs = []
+ for r_can in canonical_reps:
+ r_class_found = False
+ for r_class in r_classes:
+ if (vector(r_class[0]) - vector(r_can) in m_span
+ or vector(r_class[0]) + vector(r_can) in m_span):
+ r_class_found = True
+ break
+ if r_class_found: continue
+
+ r_classes.append([])
+ r_class = r_classes[-1]
+ r_classes_reduction_signs.append([])
+ r_class_reduction_signs = r_classes_reduction_signs[-1]
+
+ # travers short vectors with respect to m_adj until we find a
+ # vector that is equivalent to \pm r_can.modulo m \Z^l
+ for length in range(0, max_norm + 1):
+ if length >= current_max_length:
+ current_max_length += 5
+ short_vectors = m_adj.short_vector_list_up_to_length(current_max_length)
+
+ for r in short_vectors[length]:
+ if r_can - r in m_span:
+ r_class.append(tuple(r))
+ r_class_reduction_signs.append(1)
+ elif r_can + r in m_span:
+ r_class.append(tuple(r))
+ r_class_reduction_signs.append(-1)
+
+ 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])
+
+
+ return (r_classes, r_classes_reduction_signs)
+
+def higherrank_jacobi_forms(k, m, prec, algorithm="restriction"):
+ r"""
+ Compute the Fourier expansions of a basis of Jacobi forms (over
+ `\QQ`) of weight `k` and index `m` (an quadratic form over `\Z`)
+ up to given precision.
+
+ ALGORITHM:
+
+ See [Ra]. The algorithm in [Ra] is applied for precision
+ ``relation_prec``. After this, the remaining Fourier coefficients
+ are determined using as few restrictions as possible.
+
+ INPUT:
+
+ - `k` -- An integer. Only `k` modulo `2` is used.
+
+ - `m` -- A quadratic form.
+
+ - ``prec`` -- A nonnegative integer.
+
+ - ``algorithm`` -- Only "restriction" is implemented.
+
+ OUTPUT:
+
+ A list of dictionaries, which describes the Fourier expansion of
+ Jaocib forms.
+
+ EXAMPLES::
+
+ sage: from sage.modular.jacobi.higherrank import higherrank_jacobi_forms
+ sage: k = 10
+ sage: m = QuadraticForm(matrix(2, [2,1,1,2]))
+ sage: jforms = higherrank_jacobi_forms(k, m, 3)
+ sage: Sequence(jforms, cr=True)
+ [
+ {(1, (1, 1)): -45, (2, (1, 1)): -12672, (1, (0, 0)): 0, (2, (0, 0)): -59130, (0, (0, 0)): 1},
+ {(1, (1, 1)): -1/6, (2, (1, 1)): 5/3, (1, (0, 0)): 1, (2, (0, 0)): -15, (0, (0, 0)): 0}
+ ]
+
+ We access these the Fourier coefficients by means of the indices
+ `n` and `r`, typical for Jacaobi forms.
+
+ ::
+
+ sage: n = 2; r = (1, 1)
+ sage: jforms[0][(n,r)]
+ -12672
+
+ This works, since `r = (1,1)` is a reduced vector. For general
+ `r`, we have to invoke index reduction, to find a corresponding
+ index of the dictionary.
+
+ ::
+
+ sage: from sage.modular.jacobi.higherrank import higherrank_jacobi_r_classes
+ sage: from sage.modular.jacobi.higherrank import higherrank_jacobi_reduce_fe_index
+ sage: n = 2; r = (2, 2)
+ sage: m_adj = QuadraticForm(2 * m.matrix().adjoint())
+ sage: r_classes = higherrank_jacobi_r_classes(m)[0]
+ sage: m_span = m.matrix().row_module()
+ sage: ((nred, rred), _) = higherrank_jacobi_reduce_fe_index((n,r), m, r_classes, m_adj, m_span)
+ sage: jforms[0][(nred, rred)]
+ -45
+
+ TESTS:
+
+ See ``test_higherrank.py``.
+ """
+ if algorithm != "restriction":
+ raise NotImplementedError("Algorithm {} is not implemented.".format(algorithm))
+ rand = Random()
+
+
+ dim = jacobi_dimension(k, m)
+ 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
+
+ max_relation_rst_index = max_rst_index
+ relation_rst_vectors = []
+
+
+ while True:
+ try:
+ prec__max = max(prec, relation_prec)
+ jforms = _higherrank_jacobi_forms__restriction(
+ k, prec__max, relation_prec, dim,
+ *_restriction_relation_matrices(k, m,
+ prec__max, relation_prec,
+ rst_vectors, relation_rst_vectors,
+ r_classes, m_span)
+ )
+
+ if prec__max < prec:
+ return jforms
+ else:
+ return [dict(((n,r), c) for ((n,r), c) in phi.items()
+ if n < prec)
+ for phi in jforms]
+
+ except ValueError, err:
+ if len(err.args) < 2 or err.args[1] != "INSUFFICIENT RELATIONS":
+ raise
+
+ relation_prec += minimal_prec
+ max_relation_rst_index += 2
+
+
+ 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 )
+ )
+ 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 `\Z^l` find a complete
+ set of restriction vectors.
+
+ INPUT:
+
+ - `m` -- A quadratic form.
+
+ - ``r_classes`` -- A list of lists of tuples of integers.
+
+ - ``r_classes_reduction_signs`` -- A list of lists of `\pm 1`.
+
+ - ``m_span`` -- The row (or column) span `m`.
+
+ OUTPUT:
+
+ - A set of pairs, the first of which is a vector corresponding to
+ an element in `\Z^l`, and the second of which is an integer.
+
+ EXAMPLES::
+
+ sage: from sage.modular.jacobi.higherrank import _complete_set_of_restriction_vectors
+ sage: from sage.modular.jacobi.higherrank import higherrank_jacobi_r_classes
+ sage: m = QuadraticForm(matrix(2, [2,1,1,2]))
+ sage: m_span = m.matrix().row_module()
+ sage: (r_classes, r_classes_reduction_signs) = higherrank_jacobi_r_classes(m)
+ sage: _complete_set_of_restriction_vectors(m, r_classes, r_classes_reduction_signs, m_span)
+ [((1, 0), 0), ((-2, 1), 1)]
+
+ TESTS:
+
+ See ``test_higherrank.py``
+ """
+ r_classes = [map(vector, r_class) for r_class in r_classes]
+
+ length_inc = 5
+ max_length = 5
+ cur_length = 1
+ short_vectors = m.short_vector_list_up_to_length(max_length+1, True)
+
+ pm_fixed_point_indices = [ix for (ix, r_class) in enumerate(r_classes)
+ if 2*vector(r_class[0]) in m_span]
+
+ rst_vectors = []
+ nmb_rst_vectors_even = 0
+ nmb_rst_vectors_odd = 0
+ rst_kernel_even = FreeModule(QQ, len(r_classes))
+ rst_kernel_odd_matrix = identity_matrix(QQ, len(r_classes))
+ for ix in pm_fixed_point_indices:
+ rst_kernel_odd_matrix[ix,ix] = 0
+ rst_kernel_odd = rst_kernel_odd_matrix.row_module()
+
+
+ while (nmb_rst_vectors_even < len(r_classes)) :
+ while len(short_vectors[cur_length]) == 0:
+ cur_length += 1
+ if max_length < cur_length:
+ max_length += length_inc
+ short_vectors = m.short_vector_list_up_to_length(max_length+1, True)
+
+ s = vector( short_vectors[cur_length].pop() )
+
+
+ rst_imgs_even = {}
+ rst_imgs_odd = {}
+ for (cl_ix, (r_class, r_signs)) \
+ in enumerate(zip(r_classes, r_classes_reduction_signs)):
+
+ for (r, r_sign) in zip(r_class, r_signs):
+ r_rst = s.dot_product(r)
+
+ if r_rst not in rst_imgs_even:
+ rst_imgs_even[r_rst] = zero_vector(len(r_classes))
+ rst_imgs_even[r_rst][cl_ix] += 1
+
+ if cl_ix not in pm_fixed_point_indices:
+ if r_rst not in rst_imgs_odd:
+ rst_imgs_odd[r_rst] = zero_vector(len(r_classes))
+ rst_imgs_odd[r_rst][cl_ix] += r_sign
+
+ for r_rst in rst_imgs_even.keys():
+ if (rst_kernel_even.basis_matrix() * rst_imgs_even[r_rst]).is_zero():
+ continue
+ if r_rst in rst_imgs_odd:
+ contributes_to_odd = not (rst_kernel_odd.basis_matrix()
+ * rst_imgs_odd[r_rst]).is_zero()
+
+ if (nmb_rst_vectors_odd + len(pm_fixed_point_indices)
+ <= nmb_rst_vectors_even
+ and
+ (r_rst not in rst_imgs_odd or not contributes_to_odd)):
+
+ continue
+
+ rst_vectors.append((s, r_rst))
+
+ rst_kernel_even = rst_kernel_even.intersection(
+ matrix(rst_imgs_even[r_rst]).right_kernel())
+ nmb_rst_vectors_even += 1
+ if r_rst in rst_imgs_odd and contributes_to_odd:
+ rst_kernel_odd = rst_kernel_odd.intersection(
+ matrix(rst_imgs_odd[r_rst]).right_kernel())
+ nmb_rst_vectors_odd += 1
+
+
+ return rst_vectors
+
+def _restriction_relation_matrices(k, m, prec, relation_prec,
+ rst_vectors, relation_rst_vectors,
+ r_classes, m_span):
+ r"""
+ INPUT:
+
+ - `k` -- An integer. Only `k` modulo `2` is used.
+
+ - `m` -- A quadratic form.
+
+ - ``prec`` -- A nonnegative integer.
+
+ - ``relation_prec`` -- A nonnegative integer. Precision less than
+ or equal to ``prec`` up to which relations
+ of coefficients are computed.
+
+ - ``rst_vectors`` -- A list of vectors.
+
+ - ``relation_rst_vectors`` -- Compute relations for a give set of
+ restriciton vectors.
+
+ - ``r_classes`` -- A list of lists of vectors.
+
+ - ``m_span`` -- The row (or column) span `m`.
+
+ OUTPUT:
+
+ - A quintuple. See `meth:_restriction_matrix` and
+ `meth:_relation_matrix` for a more detailed description.
+
+ EXAMPLES::
+
+ sage: from sage.modular.jacobi.higherrank import _restriction_relation_matrices
+ sage: from sage.modular.jacobi.higherrank import _complete_set_of_restriction_vectors
+ sage: from sage.modular.jacobi.higherrank import higherrank_jacobi_r_classes
+ sage: m = QuadraticForm(matrix(2, [2,1,1,2]))
+ sage: m_span = m.matrix().row_module()
+ sage: prec = 2; relation_prec = 1
+ sage: (r_classes, r_classes_reduction_signs) = higherrank_jacobi_r_classes(m)
+ sage: rst_vectors_with_image = _complete_set_of_restriction_vectors(m, r_classes, r_classes_reduction_signs, m_span)
+ sage: rst_vectors = [vector(s) for s in Set(tuple(s) for (s, _) in rst_vectors_with_image)]
+ sage: relation_rst_vectors = rst_vectors + [vector(ZZ, [2,0])]
+ sage: Sequence(_restriction_relation_matrices(1, m, prec, relation_prec, rst_vectors, relation_rst_vectors, r_classes, m_span), cr=True)
+ [
+ [ 0 1 0]
+ [ 2 0 0]
+ [ 1 0 0]
+ [ 2 1 0]
+ [ 0 0 -2]
+ [ 0 0 1]
+ [ 2 0 0]
+ [ 1 0 0], [((1, 0), 1, 0, 3), ((-2, 1), 3, 3, 5)], {1: {(1, 0): 0, (0, 0): 2, (1, 1): 1}, 3: {(1, 2): 2, (1, 0): 0, (1, 3): 3, (1, 1): 1, (0, 0): 4}},
+ <BLANKLINE>
+ [(0, (0, 0)), (1, (0, 0)), (1, (1, 1))], [1], [(0, (0, 0))]
+ ]
+
+ TESTS:
+
+ Tested implicitely by ``meth:higherrank_jacobi_forms``. See also ``test_higherrank.py``.
+ """
+ (restriction_matrix__big, row_groups, row_labels, column_labels) = \
+ _restriction_matrix(k, m, prec, rst_vectors, False, r_classes, m_span)
+ (relation_matrix, column_labels_relations) = \
+ _relation_matrix(k, m, relation_prec, relation_rst_vectors, r_classes, m_span)
+
+ restriction_matrix__big = restriction_matrix__big.change_ring(QQ)
+ relation_matrix = relation_matrix.change_ring(QQ)
+
+ 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) :
+ 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.
+
+ INPUT:
+
+ - `k` -- An integer. Only `k` modulo `2` is used.
+
+ - `m` -- A quadratic form.
+
+ - ``prec`` -- A nonnegative integer.
+
+ - ``rst_vectors`` -- A list of vectors.
+
+ - ``find_relation`` -- A boolean. If ``True``, then the restrictions to
+ nonreduced indices will also be computed.
+
+ - ``r_classes`` -- A list of lists of vectors.
+
+ - ``m_span`` -- The row (or column) span `m`.
+
+ OUTPUT:
+
+ A quadruple ``(restriction_matrix, row_groups, row_labels, column_labels)``.
+
+ - ``restriction_matrix`` -- A matrix that describes the
+ restriction of Fourier expansion of Jacobi forms of index `m` so
+ classical Jacobi forms. It acts on column vectors.
+ ``row_groups``, ``row_labels``, ``column_labels`` describe which
+ entry corresponds to which Fourier index.
+
+ - ``row_groups`` -- A list of of quadruples ``(s, m_rst, start,
+ length)``. This means that rows ``start`` to ``start + length``
+ contain the image of Fourier expansions to `s z`. The index of
+ this image corresponds to classical Jacobi forms of index `m_rst`.
+
+ - ``row_labels`` -- A dictionary, which assignes to each index `m_rst`
+ of classical Jacobi forms that occur in ``row_groups`` a
+ labelling. The values of this dictionary are dictionaries
+ themselves, which map pairs `(n,r)` of integers to integer
+ indices ``ix``. This means that for every restriction whose
+ image has index `m` the row ``start + ix`` (``start`` was given
+ above) corresponds to the Fourier coefficients of index `(n,r)`.
+
+ - ``column_labels`` -- A dictionary that maps Fourier indices of
+ Jacobi forms of index `m` to integer indices `ix`. This means
+ that that the `ix`-th column of the restrictio matrix
+ corresponds to Fourier coefficients of index `(n,r)`.
+
+ EXAMPLES::
+
+ sage: from sage.modular.jacobi.higherrank import _restriction_matrix
+ sage: from sage.modular.jacobi.higherrank import _complete_set_of_restriction_vectors
+ sage: from sage.modular.jacobi.higherrank import higherrank_jacobi_r_classes
+ sage: m = QuadraticForm(ma