diff --git a/src/doc/en/reference/modfrm/index.rst b/src/doc/en/reference/modfrm/index.rst index e1f8a483a3f..f8380910357 100644 --- a/src/doc/en/reference/modfrm/index.rst +++ b/src/doc/en/reference/modfrm/index.rst @@ -37,4 +37,17 @@ Design Notes sage/modular/modform/notes +Siegel Modular Forms +==================== + +.. toctree:: + :maxdepth: 2 + + sage/modular/siegel/siegel_modular_form + sage/modular/siegel/siegel_modular_form_prec + sage/modular/siegel/siegel_modular_forms_algebra + sage/modular/siegel/siegel_modular_group + sage/modular/siegel/theta_constant + sage/modular/siegel/fastmult + .. include:: ../footer.txt diff --git a/src/module_list.py b/src/module_list.py index 3bb1999eaa9..aeb84acc894 100644 --- a/src/module_list.py +++ b/src/module_list.py @@ -981,6 +981,9 @@ def uname_specific(name, value, alternative): libraries = ["gmp", "zn_poly"], extra_compile_args=['-std=c99', '-D_XPG6']), + Extension('sage.modular.siegel.fastmult', + sources = ['sage/modular/siegel/fastmult.pyx']), + ################################ ## ## sage.modules diff --git a/src/sage/categories/pushout.py b/src/sage/categories/pushout.py index 0fca565cb00..3ca8a597795 100644 --- a/src/sage/categories/pushout.py +++ b/src/sage/categories/pushout.py @@ -1499,6 +1499,116 @@ def expand(self): return [InfinitePolynomialFunctor((x,), self._order, self._imple) for x in reversed(self._gens)] +class SiegelModularFormsAlgebraFunctor(ConstructionFunctor): + """ + this is needed to make the coercion system happy + + DO NOT REMOVE! + """ + rank = 9 + + def __init__(self, group, weights, degree, default_prec): + r""" + Initialize the functor. + + EXAMPLES:: + + sage: from sage.categories.pushout import SiegelModularFormsAlgebraFunctor + sage: F = SiegelModularFormsAlgebraFunctor('Sp(4,Z)', 'even', 2, 101) + sage: F.rank + 9 + """ + from sage.categories.all import Rings + Functor.__init__(self, Rings(), Rings()) + self._group = group + self._weights = weights + self._degree = degree + self._default_prec = default_prec + + def __call__(self, R): + r""" + Apply the functor ``self`` to the ring ``R``. + + EXAMPLES:: + + sage: from sage.categories.pushout import SiegelModularFormsAlgebraFunctor + sage: F = SiegelModularFormsAlgebraFunctor('Sp(4,Z)', 'all', 2, 41) + sage: F(QQ) + Algebra of Siegel modular forms of degree 2 and all weights on Sp(4,Z) over Rational Field + + TESTS:: + + sage: F(3) + Traceback (most recent call last): + ... + TypeError: The coefficient ring must be a ring + """ + from sage.modular.siegel.all import SiegelModularFormsAlgebra + return SiegelModularFormsAlgebra(coeff_ring=R, group=self._group, weights=self._weights, degree=self._degree, default_prec=self._default_prec) + + def __cmp__(self, other): + r""" + Compare the functor ``self`` and the object ``other``. + + EXAMPLES:: + + sage: from sage.categories.pushout import SiegelModularFormsAlgebraFunctor + sage: F = SiegelModularFormsAlgebraFunctor('Sp(4,Z)', 'even', 2, 101) + sage: G = SiegelModularFormsAlgebraFunctor('Sp(4,Z)', 'even', 2, 101) + sage: F == G + True + sage: F == SiegelModularFormsAlgebraFunctor('Sp(4,Z)', 'all', 2, 41) + False + sage: F == SiegelModularFormsAlgebraFunctor('Sp(4,Z)', 'even', 2, 61) + True + """ + c = cmp(type(self), type(other)) + if c == 0: + # default precision is an implementation detail, not a + # mathematical property, so it is ignored in comparison + c = cmp(self._group, other._group) or cmp(self._weights, other._weights) or cmp(self._degree, other._degree) + return c + + def merge(self, other): + r""" + Merge the two functors ``self`` and ``other``. + + EXAMPLES:: + + sage: from sage.categories.pushout import SiegelModularFormsAlgebraFunctor + sage: F = SiegelModularFormsAlgebraFunctor('Sp(4,Z)', 'even', 2, 101) + sage: G = SiegelModularFormsAlgebraFunctor('Gamma0(2)', 'all', 2, 61) + sage: F == G + False + sage: F.merge(G) is None + True + sage: G = SiegelModularFormsAlgebraFunctor('Sp(4,Z)', 'even', 2, 61) + sage: F.merge(G) is G + True + """ + if not isinstance(other, SiegelModularFormsAlgebraFunctor): + return None + if (self._group == other._group) and (self._weights == other._weights) and (self._degree == other._degree): + if self._default_prec <= other._default_prec: + return self + else: # self._default_prec > other._default_prec + return other + else: + return None + + def __str__(self): + r""" + Return a string describing the functor ``self``. + + EXAMPLES:: + + sage: from sage.categories.pushout import SiegelModularFormsAlgebraFunctor + sage: F = SiegelModularFormsAlgebraFunctor('Gamma0(4)', 'even', 2, 41) + sage: F + SMFAlg{"Gamma0(4)", "even", 2, 41} + """ + return 'SMFAlg{"%s", "%s", %s, %s}' %(self._group, self._weights, self._degree, self._default_prec) + class MatrixFunctor(ConstructionFunctor): """ diff --git a/src/sage/modular/all.py b/src/sage/modular/all.py index d3cc256fa89..982ddfd16ed 100644 --- a/src/sage/modular/all.py +++ b/src/sage/modular/all.py @@ -38,3 +38,5 @@ from .btquotients.all import * from .pollack_stevens.all import * + +from .siegel.all import * diff --git a/src/sage/modular/siegel/__init__.py b/src/sage/modular/siegel/__init__.py new file mode 100644 index 00000000000..2ae28399f5f --- /dev/null +++ b/src/sage/modular/siegel/__init__.py @@ -0,0 +1 @@ +pass diff --git a/src/sage/modular/siegel/all.py b/src/sage/modular/siegel/all.py new file mode 100644 index 00000000000..b77b9da840f --- /dev/null +++ b/src/sage/modular/siegel/all.py @@ -0,0 +1,3 @@ +from siegel_modular_form import SiegelModularForm +from siegel_modular_forms_algebra import SiegelModularFormsAlgebra +from siegel_modular_group import Sp4Z, Sp4Z_Gamma0_constructor as Sp4Z_Gamma0 diff --git a/src/sage/modular/siegel/fastmult.pxd b/src/sage/modular/siegel/fastmult.pxd new file mode 100644 index 00000000000..099b9da21c4 --- /dev/null +++ b/src/sage/modular/siegel/fastmult.pxd @@ -0,0 +1,11 @@ +include 'sage/ext/interrupt.pxi' +include 'sage/ext/stdsage.pxi' +include 'sage/ext/cdefs.pxi' + +from sage.rings.ring cimport Ring +from sage.rings.polynomial.polynomial_ring import is_PolynomialRing +from sage.algebras.group_algebra import GroupAlgebra + +cpdef mult_coeff_int(int a, int b, int c, coeffs_dict1, coeffs_dict2) +cpdef mult_coeff_int(int a, int b, int c, coeffs_dict1, coeffs_dict2) +cpdef mult_coeff_generic(int a, int b, int c, coeffs_dict1, coeffs_dict2, Ring R) diff --git a/src/sage/modular/siegel/fastmult.pyx b/src/sage/modular/siegel/fastmult.pyx new file mode 100644 index 00000000000..90a413d5560 --- /dev/null +++ b/src/sage/modular/siegel/fastmult.pyx @@ -0,0 +1,691 @@ +""" +Low level functions for coefficients of Siegel modular forms +""" + + +# coeffs_dict* are dictionaries of the form (a,b,c) -> Integer, where +# (a,b,c) is GL(2,Z)-reduced, i.e. |b| <= a <= c, and b^2-4ac <= 0, +# and where the keys consist either of all reduced triples in a box of the form +# (0,A)x(0,B)x(0,C) or else of all reduced triples (a,b,c) with 4ac-b^2 below +# a given bound (and c <= if 4ac-b^2=0). + +include "sage/ext/stdsage.pxi" +include "sage/ext/interrupt.pxi" + +from sage.libs.gmp.mpz cimport * +from sage.libs.gmp.mpq cimport * + +from sage.structure.element cimport Element +from sage.rings.integer_ring import ZZ +from sage.rings.integer cimport Integer + +cdef struct int_triple: + int a + int b + int c + +cdef struct long_triple: + long a + long b + long c + +cdef struct int_septuple: + int a + int b + int c + int O + int o + int U + int u + + +cdef struct int_quadruple: + int a + int b + int c + int s + + +cpdef mult_coeff_int(int a, int b, int c, coeffs_dict1, coeffs_dict2): + """ + Returns the value at (a,b,c) of the coefficient dict of the product + of the two forms with dict coeffs_dict1 and coeffs_dict2. + + It is assumed that (a,b,c) is a key in coeffs_dict1 and coeffs_dict2. + + EXAMPLES:: + + sage: from sage.modular.siegel.fastmult import mult_coeff_int + sage: A, B, C, D = SiegelModularFormsAlgebra().gens() + sage: A[(0,0,0)] + 1 + sage: mult_coeff_int(0, 0, 0, A.coeffs(), A.coeffs()) + 1 + sage: mult_coeff_int(2, 2, 2, C.coeffs(), D.coeffs()) + 1 + sage: mult_coeff_int(2, 1, 2, C.coeffs(), D.coeffs()) + 8 + """ + cdef int a1, a2 + cdef int b1, b2 + cdef int c1, c2 + cdef int B1, B2 + + cdef Integer result + + cdef mpz_t tmp + cdef mpz_t left, right + cdef mpz_t mult_coeff + + cdef mpz_t mpz_zero + + cdef long_triple tmp_triple + + result = PY_NEW(Integer) + + mpz_init(tmp) + mpz_init(left) + mpz_init(right) + mpz_init(mult_coeff) + mpz_init(mpz_zero) + + ## We assume that (a,b,c) is semipositive definite! (Havock otherwise.) + ## Reduce the form first so that our search space is smaller + tmp_triple = _reduce_GL(a, b, c) + + sig_on() + + mpz_set_si( mpz_zero, 0) + a = tmp_triple.a + b = tmp_triple.b + c = tmp_triple.c + for a1 from 0 <= a1 < a+1: + a2 = a - a1 + for c1 from 0 <= c1 < c+1: + c2 = c - c1 + mpz_set_si(tmp, 4*a1*c1) + mpz_sqrt(tmp, tmp) + B1 = mpz_get_si(tmp) + + mpz_set_si(tmp, 4*a2*c2) + mpz_sqrt(tmp, tmp) + B2 = mpz_get_si(tmp) + + for b1 from max(-B1, b - B2) <= b1 < min(B1 + 1, b + B2 + 1) : + ## Guarantees that both (a1,b1,c1) and (a2,b2,c2) are + ## positive semidefinite + b2 = b - b1 + + set_coeff(left, a1, b1, c1, coeffs_dict1) + if 0 == mpz_cmp( left, mpz_zero): continue + + set_coeff(right, a2, b2, c2, coeffs_dict2) + if 0 == mpz_cmp( right, mpz_zero): continue + + mpz_mul(tmp, left, right) + mpz_add(mult_coeff, mult_coeff, tmp) + + sig_off() + + mpz_set(result.value, mult_coeff) + mpz_clear(tmp) + mpz_clear(left) + mpz_clear(right) + mpz_clear(mult_coeff) + mpz_clear(mpz_zero) + + return result + + +cdef inline void set_coeff(mpz_t dest, int a, int b, int c, coeffs_dict): + """ + Return the value of coeffs_dict at the triple obtained from + reducing (a,b,c). + + It is assumed that the latter is a valid key + and that (a,b,c) is positive semidefinite. + + """ + cdef long_triple tmp_triple + + mpz_set_si(dest, 0) + + tmp_triple = _reduce_GL(a, b, c) + triple = (tmp_triple.a, tmp_triple.b, tmp_triple.c) + try: + mpz_set(dest, ((coeffs_dict[triple])).value) + except KeyError: + pass + return + + +cpdef mult_coeff_generic(int a, int b, int c, coeffs_dict1, coeffs_dict2, Ring R): + """ + Returns the value at (a,b,c) of the coefficient dict of the product + of the two forms with dict coeffs_dict1 and coeffs_dict2. + + It is assumed that (a,b,c) is a key in coeffs_dict1 and coeffs_dict2. + + EXAMPLES:: + + sage: from sage.modular.siegel.fastmult import mult_coeff_generic + sage: A, B, C, D = SiegelModularFormsAlgebra().gens() + sage: AB = A.satoh_bracket(B) + sage: CD = C.satoh_bracket(D) + sage: R = CD.base_ring() + sage: mult_coeff_generic(3, 2, 7, AB.coeffs(), CD.coeffs(), R) + -7993474656/5*x^4 - 18273219840*x^3*y - 134223332976/5*x^2*y^2 - 120763645536/5*x*y^3 - 17770278912/5*y^4 + """ + cdef int a1, a2 + cdef int b1, b2 + cdef int c1, c2 + cdef int B1, B2 + cdef Integer nt = PY_NEW(Integer) + nmc = R(0) + cdef long_triple tmp_triple + + cdef mpz_t tmp + + mpz_init(tmp) + + ## We assume that (a,b,c) is semipositive definite! (Havock otherwise.) + ## Reduce the form first so that our search space is smaller + tmp_triple = _reduce_GL(a, b, c) + + sig_on() + + a = tmp_triple.a + b = tmp_triple.b + c = tmp_triple.c + for a1 from 0 <= a1 < a+1: + a2 = a - a1 + for c1 from 0 <= c1 < c+1: + c2 = c - c1 + mpz_set_si(tmp, 4*a1*c1) + mpz_sqrt(tmp, tmp) + B1 = mpz_get_si(tmp) + + mpz_set_si(tmp, 4*a2*c2) + mpz_sqrt(tmp, tmp) + B2 = mpz_get_si(tmp) + + for b1 from max(-B1, b - B2) <= b1 < min(B1 + 1, b + B2 + 1) : + ## Guarantees that both (a1,b1,c1) and (a2,b2,c2) are + ## positive semidefinite + b2 = b - b1 + + nl = get_coeff(a1, b1, c1, coeffs_dict1) + if nl is None or nl.is_zero() : continue + + nr = get_coeff(a2, b2, c2, coeffs_dict2) + if nr is None or nr.is_zero() : continue + + nmc += nl*nr + + sig_off() + + mpz_clear(tmp) + + return nmc + + +cdef get_coeff(int a, int b, int c, coeffs_dict): + """ + Return the value of coeffs_dict at the triple obtained from + reducing (a,b,c). + + It is assumed that the latter is a valid key + and that (a,b,c) is positive semidefinite. + + + """ + cdef long_triple tmp_triple + + tmp_triple = _reduce_GL(a, b, c) + triple = (tmp_triple.a, tmp_triple.b, tmp_triple.c) + try: + return coeffs_dict[triple] + except KeyError: + return None + + +cpdef mult_coeff_generic_with_action(int a, int b, int c, coeffs_dict1, coeffs_dict2, Ring R): + """ + Returns the value at (a,b,c) of the coefficient dict of the product + of the two forms with dict coeffs_dict1 and coeffs_dict2. + It is assumed that (a,b,c) is a key in coeffs_dict1 and coeffs_dict2. + + .. TODO:: + + I'm not sure what this does and it's not working + """ + cdef int a1, a2 + cdef int b1, b2 + cdef int c1, c2 + cdef int B1, B2 + cdef Integer nt = PY_NEW(Integer) + nmc = R(0) + cdef long_triple tmp_triple + + cdef mpz_t tmp + + mpz_init(tmp) + + ## We assume that (a,b,c) is semipositive definite! (Havoc otherwise.) + ## Reduce the form first so that our search space is smaller + tmp_triple = _reduce_GL(a, b, c) + + sig_on() + + a = tmp_triple.a + b = tmp_triple.b + c = tmp_triple.c + for a1 from 0 <= a1 < a+1: + a2 = a - a1 + for c1 from 0 <= c1 < c+1: + c2 = c - c1 + mpz_set_si(tmp, 4*a1*c1) + mpz_sqrt(tmp, tmp) + B1 = mpz_get_si(tmp) + + mpz_set_si(tmp, 4*a2*c2) + mpz_sqrt(tmp, tmp) + B2 = mpz_get_si(tmp) + + for b1 from max(-B1, b - B2) <= b1 < min(B1 + 1, b + B2 + 1) : + ## Guarantees that both (a1,b1,c1) and (a2,b2,c2) are + ## positive semidefinite + b2 = b - b1 + + nl = get_coeff_with_action(a1, b1, c1, coeffs_dict1, R) + if nl is None or nl.is_zero() : continue + + nr = get_coeff_with_action(a2, b2, c2, coeffs_dict2, R) + if nr is None or nr.is_zero() : continue + + nmc += nl*nr + + sig_off() + + mpz_clear(tmp) + + return nmc + + +cdef get_coeff_with_action(int a0, int b0, int c0, coeffs_dict, R): + """ + Return the value of coeffs_dict at the triple obtained from + reducing (a,b,c). + + It is assumed that the latter is a valid key + and that (a,b,c) is positive semidefinite. + """ + cdef int_septuple tmp_septuple + cdef int a,b,c,d + + tmp_septuple = _xreduce_GL(a0, b0, c0) + triple = (tmp_septuple.a, tmp_septuple.b, tmp_septuple.c) + a = tmp_septuple.O + b = tmp_septuple.o + c = tmp_septuple.U + d = tmp_septuple.u + + try: + v = coeffs_dict[triple] + except KeyError: + return None + + from sage.algebras.all import GroupAlgebra + if isinstance( R, GroupAlgebra): + B = R.base_ring() + else: + B = R + + from sage.rings.polynomial.polynomial_ring import is_PolynomialRing + if is_PolynomialRing( B): + x = B.gen(0) + y = B.gen(1) + xp = x*a + y*c + yp = x*b + y*d + + if is_PolynomialRing( R): + v = v( xp, yp) + return v + + res = R(0) + for p, chi in v._fs: + if is_PolynomialRing( B): p = p(xp,yp) + G = B.group() + det = G.gen(0) + sig = G.gen(1) + if chi == G.identity(): e = 1 + elif chi == det: e = a*d-b*c + elif chi == sig: e = _sig(a,b,c,d) + else: e = (a*d-b*c) * _sig(a,b,c,d) + res += p*e*R(chi) + return res + + +cdef _sig( int a, int b, int c, int d): + a = a%2 + b = b%2 + c = c%2 + d = d%2 + if 0 == b and 0 == c: return 1 + if 1 == a+d: return 1 + return -1 + + +def reduce_GL(a, b, c): + """ + Return the GL2(ZZ)-reduced form equivalent to (positive semidefinite) + quadratic form ax^2+bxy+cy^2. + + EXAMPLES:: + + sage: from sage.modular.siegel.fastmult import reduce_GL + sage: reduce_GL(2,2,4) + (2, 2, 4) + sage: reduce_GL(3,5,3) + (1, 1, 3) + sage: reduce_GL(1,2,1) + (0, 0, 1) + """ + cdef long_triple res + + # We want to check that (a,b,c) is semipositive definite + # since otherwise we might end up in an infinite loop. + # TODO: the discriminant can become to big + if b*b-4*a*c > 0 or a < 0 or c < 0: + raise NotImplementedError("only implemented for nonpositive discriminant") + + res = _reduce_GL(a, b, c) + return (res.a, res.b, res.c) + + +cdef long_triple _reduce_GL(long a, long b, long c): + """ + Return the GL2(ZZ)-reduced form equivalent to (positive semidefinite) + quadratic form ax^2+bxy+cy^2. + (Following algorithm 5.4.2 found in Cohen's Computational Algebraic Number Theory.) + """ + cdef long_triple res + cdef long twoa, q, r + cdef long tmp + + if b < 0: + b = -b + + while not (b<=a and a<=c): ## while not GL-reduced + twoa = 2 * a + #if b not in range(-a+1,a+1): + if b > a: + ## apply Euclidean step (translation) + q = b / twoa #(2*a) + r = b % twoa #(2*a) + if r > a: + ## want (quotient and) remainder such that -a c: + #(a,c) = (c,a) + tmp = c + c = a + a = tmp + + #b=abs(b) + if b < 0: + b = -b + ## iterate + ## We're finished! Return the GL2(ZZ)-reduced form (a,b,c): + + res.a = a + res.b = b + res.c = c + + return res + + +def xreduce_GL(a, b, c): + """ + Return the GL2(ZZ)-reduced form equivalent to (positive semidefinite) + quadratic form ax^2+bxy+cy^2. + + EXAMPLES:: + + sage: from sage.modular.siegel.fastmult import xreduce_GL + sage: xreduce_GL(3,1,1) + ((1, 1, 3), (0, 1, 1, 0)) + sage: xreduce_GL(1,2,1) + ((0, 0, 1), (0, 1, 1, 1)) + sage: xreduce_GL(3,8,2) + Traceback (most recent call last): + NotImplementedError: only implemented for nonpositive discriminant + sage: xreduce_GL(1,1,1) + ((1, 1, 1), (1, 0, 0, 1)) + sage: xreduce_GL(3,2,9) + ((3, 2, 9), (1, 0, 0, 1)) + sage: xreduce_GL(3,7,9) + ((3, 1, 5), (1, 0, 1, 1)) + """ + cdef int_septuple res + + # We want to check that (a,b,c) is semipositive definite + # since otherwise we might end up in an infinite loop. + if b*b-4*a*c > 0 or a < 0 or c < 0: + raise NotImplementedError("only implemented for nonpositive discriminant") + + res = _xreduce_GL(a, b, c) + return ((res.a, res.b, res.c), (res.O, res.o, res.U, res.u)) + + +cdef int_septuple _xreduce_GL(int a, int b, int c): + """ + Return the GL2(ZZ)-reduced form f and a matrix M such that + ax^2+bxy+cy^2 = f((x,y)*M). + + .. TODO:: + + _xreduce_GL is a factor 20 slower than xreduce_GL. + How can we improve this factor ? + + """ + cdef int_septuple res + cdef int twoa, q, r + cdef int tmp + cdef int O,o,U,u + + # If [A,B,C] is the form to be reduced, then after each reduction + # step we have [A,B,C] =[a,b,c]( (x,y)*matrix(2,2,[O,o, U,u]) ) + O = u = 1 + o = U = 0 + + if b < 0: + b = -b + O = -1 + + while not (b<=a and a<=c): ## while not GL-reduced + twoa = 2 * a + #if b not in range(-a+1,a+1): + if b > a: + ## apply Euclidean step (translation) + q = b / twoa #(2*a) + r = b % twoa #(2*a) + if r > a: + ## want (quotient and) remainder such that -a c: + #(a,c) = (c,a) + tmp = c; c = a; a = tmp + tmp = O; O = o; o = tmp + tmp = U; U = u; u = tmp + + #b=abs(b) + if b < 0: + b = -b + O = -O + U = -U + + ## iterate + ## We're finished! Return the GL2(ZZ)-reduced form (a,b,c) and the O,o,U,u-matrix: + + res.a = a + res.b = b + res.c = c + res.O = O + res.o = o + res.U = U + res.u = u + + return res + + + + + +cdef inline int poly(int e, int f, int g, int h, + int i, int j, int k, int l, + int m, int n, int o, int p) : +## reverse the linebreaks; this won't work + return ( 4*(f*(k*p - l*o) - g*(j*p - l*n) + h*(j*o - k*n)) + - 6*(e*(k*p - l*o) - g*(i*p - l*m) + h*(i*o - k*m)) + + 10*(e*(j*p - l*n) - f*(i*p - l*m) + h*(i*n - j*m)) + - 12*(e*(j*o - k*n) - f*(i*o - k*m) + g*(i*n - j*m)) ) + + +def chi35(disc, A,B,C,D) : + r""" + Returns the odd weight generator the ring of Siegel modular forms of degree 2 + as in Igusa. Uses a formula from Ibukiyama [I]. + + EXAMPLES:: + + sage: A, B, C, D = SiegelModularFormsAlgebra().gens() + sage: from sage.modular.siegel.fastmult import chi35 + sage: chi = chi35(51,A,B,C,D) + sage: chi[(0,0,0)] + Traceback (most recent call last): + KeyError: (0, 0, 0) + sage: chi[(2,1,3)] + 1 + sage: chi[(3,1,4)] + -129421 + + """ + cdef dict Acoeffs, Bcoeffs, Ccoeffs, Dcoeffs + cdef int a1, a2, a3, a4 + cdef int c1, c2, c3, c4 + cdef int b1, b2, b3, b4 + cdef int a, b, c + cdef int B1, B1r, B2, B2r, B3, B4 + + cdef mpz_t coeff, acc, tmp + + cdef mpz_t mpz_zero + + mpz_init(coeff) + mpz_init(acc) + mpz_init(tmp) + mpz_init(mpz_zero) + + + sig_on() + + mpz_set_si( mpz_zero, 0) + coeffs = PY_NEW( dict ) + + prec = min([ disc, A.prec(), B.prec(), C.prec(), D.prec()]) + + if prec != disc: + raise ValueError("Insufficient precision of Igusa forms") + from siegel_modular_form_prec import SiegelModularFormPrecision + prec_class = SiegelModularFormPrecision( prec) + + Acoeffs = A.coeffs() + Bcoeffs = B.coeffs() + Ccoeffs = C.coeffs() + Dcoeffs = D.coeffs() + + for trip in prec_class.positive_forms(): + (a, b, c) = trip + mpz_set_si(coeff, 0) + for c1 from 0 <= c1 < c-1 : + for c2 from 0 <= c2 < c-1-c1 : + for c3 from 1 <= c3 < c-c1-c2 : + c1r = c - c1 + c2r = c1r - c2 + c4 = c2r - c3 + for a1 from 0 <= a1 < a-1 : + for a2 from 0<= a2 < a-1-a1 : + for a3 from 1 <= a3 < a-a1-a2 : + a1r = a - a1 + a2r = a1r - a2 + a4 = a2r - a3 + + mpz_set_si(tmp, 4*a1*c1) + mpz_sqrt(tmp, tmp) + B1 = mpz_get_si(tmp) + + mpz_set_si(tmp, 4*a1r*c1r) + mpz_sqrt(tmp, tmp) + B1r = mpz_get_si(tmp) + + mpz_set_si(tmp, 4*a2*c2) + mpz_sqrt(tmp, tmp) + B2 = mpz_get_si(tmp) + + mpz_set_si(tmp, 4*a2r*c2r) + mpz_sqrt(tmp, tmp) + B2r = mpz_get_si(tmp) + + mpz_set_si(tmp, 4*a3*c3) + mpz_sqrt(tmp, tmp) + B3 = mpz_get_si(tmp) + + mpz_set_si(tmp, 4*a4*c4) + mpz_sqrt(tmp, tmp) + B4 = mpz_get_si(tmp) + + for b1 from max(-B1, b - B1r) <= b1 < min(B1 + 1, b + B1r + 1) : + for b2 from max(-B2, b - b1 - B2r) <= b2 < min(B2 + 1, b - b1 + B2r + 1) : + for b3 from max(-B3, b - b1 - b2 - B4) <= b3 < min(B3 + 1, b - b1 - b2 + B4 + 1) : + b4 = b-(b1+b2+b3) + + set_coeff(acc, a1, b1, c1, Acoeffs) + set_coeff(tmp, a2, b2, c2, Bcoeffs) + mpz_mul(acc, acc, tmp) + set_coeff(tmp, a3, b3, c3, Ccoeffs) + mpz_mul(acc, acc, tmp) + set_coeff(tmp, a4, b4, c4, Dcoeffs) + mpz_mul(acc, acc, tmp) + + mpz_set_si(tmp, poly(a1,a2,a3,a4,b1,b2,b3,b4,c1,c2,c3,c4)) + mpz_mul(acc, acc, tmp) + + mpz_add(coeff, coeff, acc) + if mpz_cmp( coeff, mpz_zero) != 0: + coeffs[trip] = PY_NEW( Integer) + mpz_set((coeffs[trip]).value, coeff) + + for t in coeffs: + coeffs[t] = Integer(coeffs[t]/41472) + + sig_off() + + mpz_clear(coeff) + mpz_clear(acc) + mpz_clear(tmp) + mpz_clear(mpz_zero) + + return coeffs diff --git a/src/sage/modular/siegel/siegel_modular_form.py b/src/sage/modular/siegel/siegel_modular_form.py new file mode 100644 index 00000000000..cc35e941944 --- /dev/null +++ b/src/sage/modular/siegel/siegel_modular_form.py @@ -0,0 +1,1630 @@ +""" +Siegel modular forms + +Implementation of a class describing (scalar and vector valued) Siegel +modular forms of degree 2 and arbitrary group and weight. + +AUTHORS: + +- CC -- Craig Citro + +- MR -- Martin Raum + +- NR -- Nathan Ryan + +- NS -- Nils-Peter Skoruppa + +HISTORY: + +- NS -- Class SiegelModularForm_class. + +- NS -- Factory function SiegelModularForm(). + Code parts concerning the Maass lift are partly based on code + written by David Gruenewald in August 2006 which in turn is + based on the PARI/GP code by Nils-Peter Skoruppa from 2003. + +- CC, NR -- _mult_, _reduce_GL using Cython. + +- NR -- Hecke operators. + +- CC -- Rewriting of fastmult.spyx for arbitrary base rings. + +- NS -- Class Morp. + +- MR, NS -- SiegelModularFormsAlgebra_class, SiegelModularFormsAlgebra + +REFERENCES: + +.. [Sko] Nils-Peter Skoruppa, ... + +.. [I-H] Tomoyoshi Ibukiyama and Shuichi Hayashida, ... + + +USAGE: + +This class allows users to produce Siegel modular forms of degree 2. In particular, one can compute generators of rings for which the generators and known. Currently the rings are of levels 1 to 4. Furthermore, one can compute specific kinds of Siegel modular forms of degree 2: theta series attached to binary quadratic forms, theta constants, Eisenstein series, and Saito-Kurokawa and Maass lifts. Finally one can also compute certain vector-valued Siegel modular forms; in particular, those of weight (k,2). + +For each such Siegel modular form object one can compute the action of Hecke operators on it. One can also add two Siegel modular forms (provided they are both either scalar valued or vector valued and of the same weight) and multiply two Siegel modular forms (in this implementation one can, indeed, multiply scalar-valued and vector-valued Siegel modular forms). + +Siegel modular forms of degree 2 have Fourier expansions supported on positive semi-definite binary quadratic forms. The precision of a Siegel modular form `F` is measured in one of four different ways: box precision `(a,b,c)` means that the support of `F` includes all binary quadratic forms `(x,y,z)` less than `x0` means that the support of `F` contains all positive definite quadratic forms of discrimant `d` such that `-d +# +# Distributed under the terms of the GNU General Public License (GPL) +# http://www.gnu.org/licenses/ +#***************************************************************************** + + +import cPickle +import urllib +from siegel_modular_form_prec import SiegelModularFormPrecision +from sage.rings.all import ZZ +from sage.algebras.all import AlgebraElement +from siegel_modular_group import Sp4Z +from siegel_modular_group import Sp4Z_Gamma0_constructor as Sp4Z_Gamma0 +from sage.modular.modform.constructor import ModularForms + +SMF_DEFAULT_PREC = 101 + + +class SiegelModularForm_class(AlgebraElement): + r""" + Describes a Siegel modular form of degree 2. + """ + def __init__(self, parent, weight, coeffs, prec=SMF_DEFAULT_PREC, name=None): + r""" + Create a Siegel modular form of degree 2. + + INPUT: + + - ``parent`` -- SiegelModularFormsAlgebra (ambient space). + + - ``weight`` -- the weight of the form. + + - ``coeffs`` -- a dictionary whose keys are triples `(a,b,c)` + of integers representing `GL(2,\ZZ)`-reduced quadratic forms + (i.e. forms satisfying `0 \le b \le a \le c`). + + - ``name`` -- an optional string giving a name to the form, + e.g. 'Igusa_4'. This modifies the text and latex + representations of the form. + + OUTPUT: + + - a Siegel modular form + + EXAMPLES:: + + sage: E10 = ModularForms(1, 10).0 + sage: Delta = CuspForms(1, 12).0 + sage: F = SiegelModularForm(E10, Delta); F + Siegel modular form on Siegel Modular Group Sp(4,Z) of weight 10 + + TESTS:: + + sage: TestSuite(F).run() + """ + self.__group = parent.group() + self.__weight = weight + self.__coeffs = coeffs + # TODO: check whether we have enough coeffs + self.__coeff_ring = parent.coeff_ring() + self.__precision = SiegelModularFormPrecision(prec) + self.__prec = _normalized_prec(prec) + self.__name = name + self.__rep_lists = dict() + AlgebraElement.__init__(self, parent) + + def base_ring(self): + """ + Return the ring of coefficients of the Siegel modular form ``self``. + + EXAMPLES:: + + sage: S = SiegelModularFormsAlgebra() + sage: A = S.gen(0) + sage: A.base_ring() + Integer Ring + sage: C = S.gen(2) + sage: onehalfC = C * (-1/2) + sage: onehalfC.base_ring() + Rational Field + sage: B = S.gen(1) + sage: AB = A.satoh_bracket(B) + sage: AB.base_ring() + Multivariate Polynomial Ring in x, y over Rational Field + """ + return self.__coeff_ring + + def _repr_(self): + r""" + Return the plain text representation of ``self``. + + EXAMPLES:: + + sage: C = SiegelModularFormsAlgebra().gen(2) + sage: C._repr_() + 'Igusa_10' + sage: (C^2)._repr_() + 'Siegel modular form on Siegel Modular Group Sp(4,Z) of weight 20' + """ + if self.name() is None: + return 'Siegel modular form on %s of weight %d' % (self.group(), + self.weight()) + else: + return self.name() + + def _latex_(self): + r""" + Return the latex representation of ``self``. + + EXAMPLES:: + + sage: A = SiegelModularFormsAlgebra().gen(0) + sage: A._latex_() + 'Igusa_4' + sage: (A^2)._latex_() + '\\texttt{Siegel modular form on Siegel Modular Group Sp(4,Z) of weight }8' + """ + if self.name() is None: + return r'\texttt{Siegel modular form on %s of weight }%d' % (self.group(), self.weight()) + else: + return self.name() + + def group(self): + r""" + Return the modular group of ``self`` as a string. + + EXAMPLES:: + + sage: C = SiegelModularFormsAlgebra().gen(2) + sage: C.group() + Siegel Modular Group Sp(4,Z) + """ + return self.__group + + def weight(self): + r""" + Return the weight of ``self``. + + EXAMPLES:: + + sage: A, B, C, D = SiegelModularFormsAlgebra().gens() + sage: A.weight() + 4 + sage: B.weight() + 6 + sage: C.weight() + 10 + sage: D.weight() + 12 + sage: (A * B).weight() + 10 + sage: (A*B + C).weight() + 10 + """ + return self.__weight + + def precision(self): + r""" + Return the Precision class object that describes the precision of ``self``. + + EXAMPLES:: + + sage: A = SiegelModularFormsAlgebra().gen(0) + sage: A.precision() + Discriminant precision for Siegel modular form with bound 101 + """ + return self.__precision + + def __getitem__(self, key): + r""" + Return the coefficient indexed by ``key`` in the Fourier expansion + of ``self``. + + INPUT: + + - ``key`` -- a triple of integers `(a, b, c)` defining a quadratic form + + OUTPUT: + + - the coefficient of `q^{(a,b,c)}` as an element of the coefficient + ring of ``self`` + + EXAMPLES:: + + sage: A, B, C, D = SiegelModularFormsAlgebra().gens() + sage: A[(0, 0, 10)] + 272160 + sage: A[(0, 0, 25)] + 3780240 + sage: A[(5, 0, 5)] + 97554240 + sage: A[(0, 0, 26)] + Traceback (most recent call last): + ... + ValueError: precision 101 is not high enough to extract coefficient at (0, 0, 26) + sage: A[(100, 0, 1)] + Traceback (most recent call last): + ... + ValueError: precision 101 is not high enough to extract coefficient at (100, 0, 1) + sage: C[(2, 1, 3)] + 2736 + sage: C[(3, 1, 2)] + 2736 + + Any coefficient indexed by a quadratic form that is not semi-positive + definite is by definition zero:: + + sage: A[(-1, 1, 1)] + 0 + sage: A[(1, 10, 1)] + 0 + + TESTS:: + + sage: A[(1/2, 0, 6)] + Traceback (most recent call last): + ... + TypeError: the key (1/2, 0, 6) must be an integral quadratic form + """ + a, b, c = key + try: + a = ZZ(a) + b = ZZ(b) + c = ZZ(c) + except TypeError: + raise TypeError("the key %s must be an integral quadratic form" % str(key)) + + if b**2 - 4*a*c > 0 or a < 0 or c < 0: + return self.base_ring()(0) + ## otherwise we have a semi-positive definite form + if (a, b, c) in self.__coeffs: + return self.__coeffs[(a, b, c)] + ## otherwise GL2(ZZ)-reduce (a,b,c) and try again + from fastmult import reduce_GL + (a0, b0, c0) = reduce_GL(a, b, c) + if self.precision().is_in_bound((a0, b0, c0)): +## return get_coeff_with_action(a,b,c,self.coeffs(),self.base_ring()) + return self.__coeffs.get((a0, b0, c0), self.base_ring()(0)) + ## otherwise error - precision is too low + raise ValueError('precision %s is not high enough to extract coefficient at %s' % (self.prec(), str(key))) + + def coeffs(self, disc=None): + r""" + Return the dictionary of coefficients of ``self``. + + If ``disc`` is + specified, return the dictionary of coefficients indexed by + semi-positive definite quadratic forms of discriminant ``disc``. + + INPUT: + + - ``disc`` -- optional (default: ``None``); a negative integer giving + the discriminant of a semi-positive definite quadratic form + + EXAMPLES:: + + sage: A, B, C, D = SiegelModularFormsAlgebra().gens(prec=20) + sage: (3, 0, 1) in C.coeffs() + False + sage: (1, 0, 3) in C.coeffs() + True + sage: C.coeffs()[(1, 0, 3)] + -272 + sage: len(D.coeffs(disc=-16).keys()) + 2 + sage: D.coeffs(disc=-16)[(2, 0, 2)] + 17600 + + TESTS:: + + sage: A.coeffs(disc=-20) + Traceback (most recent call last): + ... + ValueError: precision is not high enough to extract coefficients of discriminant -20 + """ + if disc is None: + return self.__coeffs + elif -disc < self.prec(): + if disc < 0 and (disc % 4) in [0, 1]: + from sage.quadratic_forms.binary_qf import BinaryQF_reduced_representatives + forms = [f[:] for f in BinaryQF_reduced_representatives(disc)] + return dict([(Q, self[Q]) for Q in forms]) + else: + return {} + else: + raise ValueError('precision is not high enough to extract coefficients of discriminant %s' % disc) + + def support(self): + r""" + Return the support of ``self``, i.e. the list of reduced quadratic + forms indexing non-zero Fourier coefficients of ``self``. + + EXAMPLES:: + + sage: A, B, C, D = SiegelModularFormsAlgebra().gens(prec=15) + sage: A.support() + [(0, 0, 0), (0, 0, 1), (0, 0, 2), (0, 0, 3), (1, 0, 1), (1, 0, 2), (1, 0, 3), (1, 1, 1), (1, 1, 2), (1, 1, 3), (2, 2, 2)] + sage: (0, 0, 2) in C._keys() + True + sage: C[(0, 0, 2)] + 0 + sage: (0, 0, 2) in C.support() + False + """ + return [k for k in self._keys() if self[k] != 0] + + def max_disc(self): + r""" + Return the largest discriminant corresponding to a non-zero Fourier + coefficient of ``self``. + + Note that since these discriminants are all non-positive, this + corresponds in some sense to the "smallest" non-zero Fourier + coefficient. It is analogous to the valuation of a power series + (smallest degree of a non-zero term). + + EXAMPLES:: + + sage: gens = SiegelModularFormsAlgebra().gens() + sage: [f.max_disc() for f in gens] + [0, 0, -3, -3] + """ + return max([b**2 - 4*a*c for (a, b, c) in self.support()]) + + def _keys(self): + r""" + Return the keys of the dictionary of coefficients of ``self``. + + EXAMPLES:: + + sage: A, B, C, D = SiegelModularFormsAlgebra().gens(prec=15) + sage: A._keys() + [(0, 0, 0), (0, 0, 1), (0, 0, 2), (0, 0, 3), (1, 0, 1), (1, 0, 2), (1, 0, 3), (1, 1, 1), (1, 1, 2), (1, 1, 3), (2, 2, 2)] + """ + return sorted(self.coeffs().keys()) + + def prec(self): + r""" + Return the precision of ``self``. + + EXAMPLES:: + + sage: B = SiegelModularFormsAlgebra().gen(1) + sage: B.prec() + 101 + """ + return self.precision().prec() + + def name(self): + r""" + Return the name of the Siegel modular form ``self``, or None if + ``self`` does not have a custom name. + + EXAMPLES:: + + sage: A, B = SiegelModularFormsAlgebra().gens()[:2] + sage: A.name() + 'Igusa_4' + sage: (A*B).name() is None + True + """ + return self.__name + + def truncate(self, prec): + r""" + Return a Siegel modular form with the same parent as ``self``, but + with Fourier expansion truncated at the new precision ``prec``. + + EXAMPLES:: + + sage: A = SiegelModularFormsAlgebra().gen(0, prec=20); A + Igusa_4 + sage: A.prec() + 20 + sage: len(A._keys()) + 18 + sage: A[(1, 1, 5)] + 1330560 + sage: At = A.truncate(16); At + Igusa_4 + sage: At.prec() + 16 + sage: len(At._keys()) + 14 + sage: At[(1, 1, 5)] + Traceback (most recent call last): + ... + ValueError: precision 16 is not high enough to extract coefficient at (1, 1, 5) + sage: At[(1, 1, 4)] == A[(1, 1, 4)] + True + + TESTS:: + + sage: A.truncate(30) + Traceback (most recent call last): + ... + ValueError: truncated precision 30 cannot be larger than the current precision 20 + """ + if prec > self.prec(): + raise ValueError('truncated precision %s cannot be larger than the current precision %s' % (prec, self.prec())) + else: + return _SiegelModularForm_from_dict(group=self.group(), + weight=self.weight(), + coeffs=self.coeffs(), + prec=prec, parent=None, + name=self.name()) + + ################ + ## Arithmetic ## + ################ + + def _add_(left, right): + r""" + Return the sum of ``left`` and ``right``. + + There are some compatibility conditions that the forms ``left`` + and ``right`` must satisfy in order for addition to work: + + - they must have the same weight + - they must be defined on the same group, or one of them must + be defined on the whole group 'Sp(4,Z)' + + The precision of the sum of ``left`` and ``right`` is the minimum + of the precisions of ``left`` and ``right``. + + EXAMPLES:: + + sage: A, B, C, D = SiegelModularFormsAlgebra().gens() + sage: twoA = A + A + sage: A2 = A * 2 + sage: A2 == twoA + True + sage: A*B + C + Siegel modular form on Siegel Modular Group Sp(4,Z) of weight 10 + sage: E = C._add_(A*B) + sage: E + Siegel modular form on Siegel Modular Group Sp(4,Z) of weight 10 + sage: C[(1, 1, 1)] + 1 + sage: (A*B)[(1, 1, 1)] + 57792 + sage: E[(1, 1, 1)] + 57793 + + sage: F = SiegelModularForm(Sp4Z_Gamma0(2), 4, A.coeffs()) + sage: F + Siegel modular form on Siegel Modular Group Gamma0(2) of weight 4 + sage: A + F + Siegel modular form on Siegel Modular Group Gamma0(2) of weight 4 + sage: F + A == F + A + True + sage: (A + F)[(1, 1, 1)] == A[(1, 1, 1)] + F[(1, 1, 1)] + True + sage: A1 = SiegelModularFormsAlgebra().gen(0, prec=20) + sage: (A1 + F).prec() + 20 + sage: (A1 + F)[(1, 1, 1)] == (A + F)[(1, 1, 1)] + True + + TESTS:: + + sage: A + B + Traceback (most recent call last): + ... + ValueError: cannot add Siegel modular forms of different weights + + sage: F = SiegelModularForm(Sp4Z_Gamma0(2), 4, A.coeffs()) + sage: G = SiegelModularForm(Sp4Z_Gamma0(5), 4, A.coeffs()) + sage: F + G + Traceback (most recent call last): + ... + TypeError: unsupported operand parent(s) for '+': 'Algebra of Siegel modular forms of degree 2 and even weights on Siegel Modular Group Gamma0(2) over Integer Ring' and 'Algebra of Siegel modular forms of degree 2 and even weights on Siegel Modular Group Gamma0(5) over Integer Ring' + sage: Z = A + (-A) + sage: Z + Siegel modular form on Siegel Modular Group Sp(4,Z) of weight 4 + sage: Z.prec() + 101 + sage: Z.coeffs() + {} + """ + # take the minimum of the precisions and truncate the forms if necessary + lp = left.prec() + rp = right.prec() + if lp < rp: + prec = lp + right = right.truncate(prec) + elif rp < lp: + prec = rp + left = left.truncate(prec) + else: + prec = lp + + # figure out the group of the sum = intersection of the groups + # of the two forms + # right now, we only allow two identical groups, or different + # groups where one of them is all of 'Sp(4,Z)' + lg = left.group() + rg = right.group() + if lg is None or lg == Sp4Z: + group = rg + elif rg is None or rg == Sp4Z: + group = lg + elif lg != rg: + raise NotImplementedError("addition of forms on different groups " + "not yet implemented") + else: + group = lg + + # the sum of two Siegel modular forms is another Siegel modular + # form only if the weights are equal + if left.weight() != right.weight(): + raise ValueError('cannot add Siegel modular forms of different ' + 'weights') + else: + wt = left.weight() + + # compute the coefficients of the sum + d = dict() + for f in set(left._keys() + right._keys()): + v = left[f] + right[f] + if not v.is_zero(): + d[f] = v + par = left.parent() + from sage.modular.siegel.siegel_modular_forms_algebra import SiegelModularFormsAlgebra + par2 = SiegelModularFormsAlgebra(coeff_ring=par.coeff_ring(), group=group, weights=par.weights(), degree=par.degree()) + return par2.element_class(parent=par, weight=wt, coeffs=d, prec=prec) + + def _mul_(left, right): + r""" + Return the product of the Siegel modular form ``left`` by the + element ``right``, which can be a scalar or another Siegel modular + form. + + The multiplication of two forms on arbitrary groups is not yet + implemented. At the moment, the forms must either be on the same + group, or one of the groups must be the full `Sp(4,Z)`. + + EXAMPLES:: + + sage: A, B, C, D = SiegelModularFormsAlgebra().gens() + sage: A3 = A * 3 + sage: B2 = 2 * B + sage: Conehalf = C / 2 + sage: Donehalf = (1/2) * D + sage: A3[(1, 1, 1)] / A[(1, 1, 1)] + 3 + sage: B2[(2, 1, 3)] / B[(2, 1, 3)] + 2 + sage: Conehalf[(2, 0, 1)] / C[(2, 0, 1)] + 1/2 + sage: Donehalf[(3, 3, 3)] / D[(3, 3, 3)] + 1/2 + sage: E = C._mul_(D) + sage: E[(0, 0, 5)] + 0 + sage: E[(2, 3, 3)] + 8 + + TESTS:: + + sage: A = SiegelModularFormsAlgebra(default_prec=61).0 + sage: A.prec() + 61 + sage: C = SiegelModularFormsAlgebra(default_prec=81).2 + sage: C.prec() + 81 + sage: A*C == C*A + True + sage: (A*C).prec() + 64 + """ + from sage.modular.siegel.siegel_modular_forms_algebra import SiegelModularFormsAlgebra + par = left.parent() + from sage.rings.all import infinity + if left.prec() is infinity: + tmp = right + right = left + left = tmp + + wt = left.weight() + right.weight() + if wt % 2: + weights = 'all' + else: + weights = 'even' + if right.prec() is infinity: + prec = left.prec() + c = right[(0, 0, 0)] + d = dict() + for f in left.coeffs(): + v = left[f] * c + if not v.is_zero(): + d[f] = v + par = SiegelModularFormsAlgebra(coeff_ring=par.coeff_ring(), + group=par.group(), + weights=weights, + degree=par.degree()) + return par.element_class(parent=par, weight=left.weight(), + coeffs=d, prec=prec) + + lp = left.prec() + rp = right.prec() + if lp < rp: + right = right.truncate(lp) + elif rp < lp: + left = left.truncate(rp) + prec = min(lp - right.max_disc(), rp - left.max_disc()) + + if left.group() is None: + group = right.group() + elif right.group() is None: + group = left.group() + elif left.group() != right.group(): + raise NotImplementedError("multiplication for differing groups " + "not yet implemented") + else: + group = left.group() + + s1, s2, R = left.coeffs(), right.coeffs(), left.base_ring() + d = dict() + _prec = SiegelModularFormPrecision(prec) + if ZZ == R: + from fastmult import mult_coeff_int + for x in _prec: + v = mult_coeff_int(x[0], x[1], x[2], s1, s2) + if not v.is_zero(): + d[x] = v + elif left.parent().base_ring() == left.parent().coeff_ring(): + from fastmult import mult_coeff_generic + for x in _prec: + v = mult_coeff_generic(x[0], x[1], x[2], s1, s2, R) + if not v.is_zero(): + d[x] = v + else: + from fastmult import mult_coeff_generic_with_action + for x in _prec: + v = mult_coeff_generic_with_action(x[0], x[1], x[2], s1, s2, R) + if not v.is_zero(): + d[x] = v + + par = SiegelModularFormsAlgebra(coeff_ring=par.coeff_ring(), group=group, weights=weights, degree=par.degree()) + return par.element_class(parent=par, weight=wt, coeffs=d, prec=prec) + + def _rmul_(self, c): + r""" + Right multiplication -- only by scalars. + + EXAMPLES:: + + sage: A, B, C, D = SiegelModularFormsAlgebra().gens() + sage: onehalfC = C._rmul_(-1/2) + sage: C[(2, 1, 3)] + 2736 + sage: onehalfC[(2, 1, 3)] + -1368 + + TESTS:: + + sage: Z = A._rmul_(0) + sage: Z + Siegel modular form on Siegel Modular Group Sp(4,Z) of weight 4 + sage: Z.coeffs() + {} + """ + d = dict() + for f in self.coeffs(): + v = self[f] * c + if not v.is_zero(): + d[f] = v + par = self.parent() + return par.element_class(parent=par, weight=self.weight(), coeffs=d, prec=self.prec()) + + def _lmul_(self, c): + r""" + Left multiplication -- only by scalars. + + EXAMPLES:: + + sage: A = SiegelModularFormsAlgebra(default_prec=51).0 + sage: A3 = A._lmul_(3) + sage: A[(1, 1, 1)] + 13440 + sage: A3[(1, 1, 1)] + 40320 + """ + return self._rmul_(c) + + def __eq__(left, right): + r""" + Return True if the Siegel modular forms ``left`` and ``right`` + have the same group and weight, and the same coefficients up to + the smaller of the two precisions. + + EXAMPLES:: + + sage: A, B, C, D = SiegelModularFormsAlgebra(default_prec=20).gens() + sage: A2, B2, C2, D2 = SiegelModularFormsAlgebra().gens() + sage: A2 == A + True + sage: A == A + True + sage: A == B + False + + TESTS:: + + sage: E = SiegelModularForm(Sp4Z_Gamma0(2), 6, B.coeffs()) + sage: E == B + False + sage: S = SiegelModularFormsAlgebra(default_prec=21) + sage: S(0) == S(0) + True + sage: S(0) == S(1) + False + """ + if not isinstance(right, SiegelModularForm_class): + return False + if left.group() != right.group(): + return False + if left.weight() != right.weight(): + return False + + # we compare up to the minimum of the two precisions + new_prec = min(left.prec(), right.prec()) + # if new_prec is infinity, then left and right are both + # constant Siegel modular forms + from sage.rings.all import infinity + if new_prec is infinity: + return left[0, 0, 0] == right[0, 0, 0] + + left = left.truncate(new_prec) + right = right.truncate(new_prec) + new_smf_prec = SiegelModularFormPrecision(new_prec) + + lc = left.coeffs() + rc = right.coeffs() + # make coefficients that are implicitly zero be actually zero + for k in new_smf_prec: + if not k in lc: + lc[k] = 0 + if not k in rc: + rc[k] = 0 + + return lc == rc + + def pickle(self, file, format=1, name=None): + r""" + Dump to a file in a portable format. + + EXAMPLES:: + + sage: A = SiegelModularFormsAlgebra(default_prec=20).0 + sage: FILE = open('A.sobj', 'w') + sage: A.pickle(FILE) + + .. TODO:: + + I don't really think I have the right syntax since I could not + load it after doing this. + """ + from sage.rings.all import QQ + if self.base_ring().fraction_field() == QQ: + pol = [0, 1] + coeffs = self.coeffs() + else: + pol = self.base_ring().polynomial().list() + coeffs = dict() + for k in self.coeffs(): + coeffs[k] = self.coeffs()[k].list() + if None == name: + _name = self._repr_() + f1 = 'format %d: [this string, wt, name, pol., prec., dict. of coefficients]' % 1 + data = [f1, self.weight(), _name, pol, self.prec(), coeffs] + cPickle.dump(data, file) + return + + def fourier_jacobi_coeff(self, N): + r""" + Return the `N`-th Fourier-Jacobi coefficient of the Siegel modular + form ``self`` as a dictionary indexed by pairs `(disc,r)` with + `disc<0` and `r^2 \equiv disc \bmod 4N`. + + EXAMPLES:: + + sage: A = SiegelModularFormsAlgebra(default_prec=51).0 + sage: fj = A.fourier_jacobi_coeff(1) + sage: fj[(-8, 0)] + 181440 + """ + prec = self.prec() + NN = 4*N + keys = [(disc, r) for disc in range(prec) + for r in range(NN) if (r**2 - disc) % NN == 0] + coeff = dict() + for disc, r in keys: + (a, b, c) = (-(r**2 - disc)/NN, r, N) + coeff[(-disc, r)] = self[(a, b, c)] + return coeff + + def satoh_bracket(left, right, names=['x', 'y']): + r""" + Return the Satoh bracket [``left``, ``right``], where ``left`` + and ``right`` are scalar-valued Siegel modular forms. + + The result is a vector-valued Siegel modular form whose coefficients + are polynomials in the two variables specified by ``names``. + + EXAMPLES:: + + sage: A, B = SiegelModularFormsAlgebra().gens()[:2] + sage: AB = A.satoh_bracket(B) + sage: AB[(1, 0, 1)] + -20160*x^2 - 40320*y^2 + """ + from sage.rings.all import PolynomialRing + R = PolynomialRing(ZZ, names) + x, y = R.gens() + d_left_coeffs = dict((f, (f[0]*x*x + f[1]*x*y + f[2]*y*y)*left[f]) + for f in left.coeffs()) + d_left = SiegelModularForm(left.group(), left.weight(), + d_left_coeffs, left.prec()) + d_right_coeffs = dict((f, (f[0]*x*x + f[1]*x*y + f[2]*y*y)*right[f]) + for f in right.coeffs()) + d_right = SiegelModularForm(right.group(), right.weight(), + d_right_coeffs, right.prec()) + return d_left*right/left.weight() - d_right*left/right.weight() + + ################# + # Hecke operators + ################# + + def hecke_image(self, ell): + r""" + Return the Siegel modular form which is the image of ``self`` under + the Hecke operator T(``ell``). + + EXAMPLES:: + + sage: A, B, C, D = SiegelModularFormsAlgebra().gens() + sage: Ups20 = -1/1785600*A*B*C - 1/1785600*A^2*D + C^2 + sage: TUps20 = Ups20.hecke_image(2) + sage: TUps20[(1, 1, 1)] / Ups20[(1, 1, 1)] + -840960 + sage: TUps20 + Siegel modular form on Siegel Modular Group Sp(4,Z) of weight 20 + + sage: A.hecke_image(5) + T(5)(Igusa_4) + """ + # TODO: compute precision for T(ell)F + from sage.functions.all import ceil + from sage.rings.arith import gcd + if gcd(ell, self.group().level()) > 1: + raise NotImplementedError("Implementation of Hecke operators T_ell " + "for ell coprime to the level") + try: + # TODO: I am not sure whether this sets the right prec + a, b, c = self.prec() + prec = (ceil(a/ell), ceil(b/ell), ceil(c/ell)) + except TypeError: + prec = ceil(self.prec()/ell/ell) + prec = _normalized_prec(prec) + d = dict((f, self.hecke_coefficient(ell, f)) for f in self._keys() if _is_bounded(f, prec)) + if self.name(): + name = 'T(' + str(ell) + ')(' + self.name() + ')' + else: + name = None + par = self.parent() + from sage.modular.siegel.siegel_modular_forms_algebra import SiegelModularFormsAlgebra + par = SiegelModularFormsAlgebra(coeff_ring=par.coeff_ring(), + group=par.group(), + weights=par.weights(), + degree=par.degree()) + return par.element_class(parent=par, weight=self.weight(), coeffs=d, + prec=prec, name=name) + + def hecke_coefficient(self, ell, (a, b, c)): + r""" + Return the `(a, b, c)`-indexed coefficient of the image of ``self`` + under the Hecke operator T(``ell``). + + EXAMPLES:: + + sage: A, B, C, D = SiegelModularFormsAlgebra().gens() + sage: Ups20 = -1/1785600*A*B*C - 1/1785600*A^2*D + C^2 + sage: Ups20.hecke_coefficient(5, (1, 1, 1)) + 5813608045/992 + sage: Ups20.hecke_coefficient(5, (1, 0, 1)) + 5813608045/248 + sage: Ups20.hecke_coefficient(5, (1, 0, 1)) / Ups20[(1, 0, 1)] + -5232247240500 + """ + from sage.rings.arith import gcd + if gcd(ell, self.group().level()) > 1: + raise NotImplementedError("Implementation of Hecke operators " + "T_ell for ell coprime to the level") + k = self.weight() + coeff = 0 + from sage.rings.all import divisors + from sage.quadratic_forms.binary_qf import BinaryQF + qf = BinaryQF([a, b, c]) + for t1 in divisors(ell): + for t2 in divisors(t1): + cosets = self._get_representatives(ell, t1/t2) + for V in cosets: + aprime, bprime, cprime = qf.matrix_action_right(V)[:] + if aprime % t1 == 0 and bprime % t2 == 0 and cprime % t2 == 0: + try: + coeff = coeff + t1**(k-2)*t2**(k-1)*self[(ell*aprime/t1**2, ell*bprime/t1/t2, ell*cprime/t2**2)] + except KeyError, msg: + raise ValueError('%s' % (self, msg)) + return coeff + + def _get_representatives(self, ell, t): + r""" + A helper function used in hecke_coefficient that computes the right + coset representatives of `\Gamma^0(t)\SL(2,Z)` where `\Gamma^0(t)` + is the subgroup of `SL(2,Z)` where the upper left hand corner is + divisible by `t`. + + EXAMPLES:: + + sage: A = SiegelModularFormsAlgebra().0 + sage: A._get_representatives(5, 3) + [ + [ 0 1] [1 0] [1 1] [1 2] + [-1 0], [0 1], [0 1], [0 1] + ] + + .. NOTE:: + + We use the bijection $\Gamma^0(t)\SL(2,Z) \rightarrow P^1(\Z/t\Z)$ + given by $A \mapsto [1:0]A$. + """ + try: + return self.__rep_lists[(ell, t)] + except KeyError: + from sage.matrix.all import MatrixSpace + M = MatrixSpace(ZZ, 2, 2) + if t == 1: + return [M([1, 0, 0, 1])] + from sage.modular.all import P1List + from sage.rings.all import xgcd + rep_list = [] + for x, y in P1List(t): + e, d, c = xgcd(x, y) + rep = M([x, y, -c, d]) + rep_list.append(rep) + self.__rep_lists[(ell, t)] = rep_list + return rep_list + + +def SiegelModularForm(arg0, arg1=None, arg2=None, prec=None, name=None, + hint=None): + r""" + Construct a Siegel modular form. + + INPUT: + + Supported formats + + 1. SiegelModularForm(f, g): + creates the Siegel modular form VI(f,g) from the elliptic + modular forms f, g as in [Sko]_. + Shortforms: + + - SiegelModularForm(f) for SiegelModularForm(f, 0-form), + - SiegelModularForm(f, 0) for SiegelModularForm(f, 0-form), + - SiegelModularForm(0, f) for SiegelModularForm(0-form, f). + + 2. SiegelModularForm(x): + if x is an element of a commutative ring creates the constant + Siegel modular form with value x. + + 3. SiegelModularForm(string): + creates the Siegel modular form pickled at the location + described by string (url or path_to_file). + + 4. SiegelModularForm(qf): + constructs the degree 2 theta series associated to the given + quadratic form. + + 5. SiegelModularForm(f): + --- TODO: Borcherds lift + + 6. SiegelModularForm(f, g): + --- TODO: Yoshida lift + + 7. SiegelModularForm(repr): + --- TODO: Lift from Weil representation + + 8. SiegelModularForm(char): + --- TODO: Singular weight + + 9. SiegelModularForm(group, weight, dict): + if group is a string, weight an integer and dict a dictionary, + creates a Siegel modular form whose coefficients are contained + in dict. + + EXAMPLES:: + + sage: M4 = ModularForms(1, 4) + sage: M6 = ModularForms(1, 6) + sage: E4 = M4.basis()[0] + sage: F = SiegelModularForm(E4, M6(0), 16); F + Siegel modular form on Siegel Modular Group Sp(4,Z) of weight 4 + + sage: url = 'http://sage.math.washington.edu/home/nils/Siegel-Modular-Forms/data/upsilon-forms_20-32_1000/Upsilon_20.p' # optional -- internet + sage: H = SiegelModularForm(url); H # optional - internet + Upsilon_20 + sage: H[(4, 4, 4)] # optional - internet + 248256200704 + + sage: url = 'http://sage.math.washington.edu/home/nils/Siegel-Modular-Forms/data/upsilon-forms_20-32_1000/Upsilon_28.p' # optional -- internet + sage: L = SiegelModularForm(url); L # optional - internet + Upsilon_28 + sage: L.prec() # optional - internet + 1001 + sage: L[(3, 3, 3)] # optional - internet + -27352334316369546*a^2 + 3164034718941090318*a + 2949217207771097198880 + sage: L[(3, 3, 3)].parent() # optional - internet + Number Field in a with defining polynomial x^3 - x^2 - 294086*x - 59412960 + + sage: l = [(0, 0, 0, 0)]*2 + [(1, 0, 0, 0)] + [(0, 1, 0, 0)] + [(0, 0, 0, 1)] + [(0, 0, 1, 1)] + sage: char_dict = {tuple(l): 1/4} + sage: F = SiegelModularForm(char_dict, prec=100, name="van Geemen F_7"); F + van Geemen F_7 + sage: F.coeffs()[(1, 0, 9)] + -7 + + sage: Q37 = QuadraticForm(ZZ, 4, [1,0,1,1, 2,2,3, 10,2, 6]) + sage: Q37.level() + 37 + sage: F37 = SiegelModularForm(Q37, prec=100); F37 + Siegel modular form on Siegel Modular Group Gamma0(37) of weight 2 + + sage: F37[(2, 1, 3)] + 0 + sage: F37[(2, 1, 7)] + 4 + """ + try: + from sage.structure.coerce import py_scalar_to_element + arg0 = py_scalar_to_element(arg0) + except TypeError: + pass + + from sage.modular.modform.element import is_ModularFormElement + from sage.rings.all import RingElement + from sage.modular.siegel.siegel_modular_group import Sp4Z + if is_ModularFormElement(arg0) and is_ModularFormElement(arg1): + return _SiegelModularForm_as_Maass_spezial_form(arg0, arg1, prec, name) + if is_ModularFormElement(arg0) and (0 == arg1 or arg1 is None): + M = ModularForms(1, arg0.weight() + 2) + return _SiegelModularForm_as_Maass_spezial_form(arg0, M(0), prec, name) + if 0 == arg0 and is_ModularFormElement(arg1): + M = ModularForms(1, arg1.weight() - 2) + return _SiegelModularForm_as_Maass_spezial_form(M(0), arg1, prec, name) + from sage.quadratic_forms.all import QuadraticForm + if isinstance(arg0, QuadraticForm): + return _SiegelModularForm_from_QuadraticForm(arg0, prec, name) + if isinstance(arg0, RingElement) and arg1 is None: + from sage.rings.all import infinity + return _SiegelModularForm_from_dict(group=Sp4Z, weight=0, coeffs={(0, 0, 0): arg0}, prec=infinity) + if isinstance(arg0, str) and arg1 is None: + return _SiegelModularForm_from_file(arg0) + if isinstance(arg0, dict) and arg1 is None: + return _SiegelModularForm_from_theta_characteristics(arg0, prec=prec, name=name, hint=hint) + from sage.rings.all import Integer + from sage.modular.siegel.siegel_modular_group import GSp4_Arithmetic_Subgroups + if isinstance(arg0, GSp4_Arithmetic_Subgroups) and isinstance(arg1, (int, Integer)) and isinstance(arg2, dict): + return _SiegelModularForm_from_dict(group=arg0, weight=arg1, coeffs=arg2, prec=prec, name=name) + raise TypeError("wrong arguments") + + +def _SiegelModularForm_as_Maass_spezial_form(f, g, prec=SMF_DEFAULT_PREC, name=None): + """ + Return the Siegel modular form I(f,g) (Notation as in [Sko]_). + + EXAMPLES:: + + sage: M14 = ModularForms(group=1, weight=14) + sage: E14 = M14.eisenstein_subspace() + sage: f = M14.basis()[0] + sage: S16 = ModularForms(group=1, weight=16).cuspidal_subspace() + sage: g = S16.basis()[0] + sage: from sage.modular.siegel.siegel_modular_form import _SiegelModularForm_as_Maass_spezial_form + sage: IFG = _SiegelModularForm_as_Maass_spezial_form(f, g, prec=100, name=None) + sage: IFG[(2, 1, 3)] + -1080946527072 + + INPUT: + + - f -- modular form of level 1 + - g -- cusp form of level 1 amd wt = wt of f + 2 + - prec -- either a triple (amax,bmac,cmax) or an integer Dmax + """ + k = f.weight() + assert(k+2 == g.weight()) | (f == 0) | (g == 0), "incorrect weights!" + assert(g.q_expansion(1) == 0), "second argument is not a cusp form" + + from sage.modular.siegel.siegel_modular_group import Sp4Z + + if isinstance(prec, tuple): + (amax, bmax, cmax) = prec + amax = min(amax, cmax) + bmax = min(bmax, amax) + clean_prec = (amax, bmax, cmax) + if bmax <= 0: + # no reduced forms below prec + return _SiegelModularForm_from_dict(group=Sp4Z, weight=k, coeffs=dict(), prec=0) + if 1 == amax: + # here prec = (0,0,>=0) + Dtop = 0 + else: + Dtop = 4*(amax-1)*(cmax-1) + else: + clean_prec = max(0, prec) + if 0 == clean_prec: + # no reduced forms below prec + return _SiegelModularForm_from_dict(group=Sp4Z, weight=k, coeffs=dict(), prec=0) + while 0 != clean_prec % 4 and 1 != clean_prec % 4: + clean_prec -= 1 + Dtop = clean_prec - 1 + precision = (Dtop+1)//4 + 1 + # TODO: examine error when called with 1 == prec + if 1 == precision: + precision = 2 + + """ + Create the Jacobi form I(f,g) as in [Sko]. + + It suffices to construct for all Jacobi forms phi only the part + sum_{r=0,1;n} c_phi(r^2-4n) q^n zeta^r. + When, in this code part, we speak of Jacobi form we only mean this part. + + We need to compute Ifg = \sum_{r=0,1; n} c(r^2-4n) q^n zeta^r up to + 4n-r^2 <= Dtop, i.e. n < precision + """ + from sage.rings.all import PowerSeriesRing, QQ + PS = PowerSeriesRing(QQ, name='q') + q = PS.gens()[0] + + ## Create the quasi Dedekind eta^-6 power series: + pari_prec = max(1, precision - 1) + # next line yields error if 0 == pari_prec + from sage.libs.pari.all import pari + from sage.rings.all import O + pari.set_series_precision(pari_prec) + eta_quasi = PS(pari('Vec(eta(q))')) + O(q**precision) + etapow = eta_quasi**-6 + + ## Create the Jacobi forms A=a*etapow and B=b*etapow in stages. + ## Recall a = sum_{s != r mod 2} s^2*(-1)^r*q^((s^2+r^2-1)/4)*zeta^r + ## b = sum_{s != r mod 2} (-1)^r*q^((s^2+r^2-1)/4)*zeta^r + ## r, s run over ZZ but with opposite parities. + ## For r=0, we need s odd, (s^2-1)/4 < precision, with s=2t+1 hence t^2+t < precision. + ## For r=1, we need s even, s^2/4 < precision, with s=2t hence t^2 < precision. + + from sage.misc.all import xsrange + + a1 = -2*sum((2*t)**2 * q**(t**2) for t in xsrange(1, precision) if t*t < precision) + b1 = -2*sum(q**(t**2) for t in xsrange(1, precision) if t*t < precision) + a0 = 2*sum((2*t+1)**2 * q**(t**2+t) for t in xsrange(precision) if t*t + t < precision) + b0 = 2*sum(q**(t**2+t) for t in xsrange(precision) if t*t + t < precision) + b1 = b1 - 1 + + ## Form A and B - the Jacobi forms used in [Sko]'s I map. + + (A0, A1, B0, B1) = (a0*etapow, a1*etapow, b0*etapow, b1*etapow) + + ## Calculate the image of the pair of modular forms (f,g) under + ## [Sko]'s isomorphism I : M_{k} \oplus S_{k+2} -> J_{k,1}. + + fderiv = PS(q * f.qexp(precision).derivative()) + (f, g) = (PS(f.qexp(precision)), PS(g.qexp(precision))) + + ## Finally: I(f,g) is given by the formula below: + Ifg0 = k/2*f*A0 - fderiv*B0 + g*B0 + O(q**precision) + Ifg1 = k/2*f*A1 - fderiv*B1 + g*B1 + O(q**precision) + + ## For applying the Maass' lifting to genus 2 modular forms. + ## we put the coefficients og Ifg into a dictionary Chi + ## so that we can access the coefficient corresponding to + ## discriminant D by going Chi[D]. + + ## Note: Ifg.list()[i].list()[j] gives the coefficient of q^i*zeta^j + + Cphi = {0: 0} # initialise dictionary. Value changed in the loop if we have a 'noncusp form' + qcoeffs0 = Ifg0.list() + qcoeffs1 = Ifg1.list() + for i in xsrange(len(qcoeffs0)): + Cphi[-4*i] = qcoeffs0[i] + Cphi[1-4*i] = qcoeffs1[i] + + ## the most negative discriminant occurs when i is largest + ## and j is zero. That is, discriminant 0^2-4*i + ## Note that i < precision. + maxD = -4*i + + """ + Create the Maass lift F := VI(f,g) as in [Sko]. + """ + + ## The constant term is given by -Cphi[0]*B_{2k}/(4*k) + ## (note in [Sko] this coeff has typos). + ## For nonconstant terms, + ## The Siegel coefficient of q^n * zeta^r * qdash^m is given + ## by the formula \sum_{ a | gcd(n,r,m) } Cphi[D/a^2] where + ## D = r^2-4*n*m is the discriminant. + ## Hence in either case the coefficient + ## is fully deterimined by the pair (D,gcd(n,r,m)). + ## Put (D,t) -> \sum_{ a | t } Cphi[D/a^2] + ## in a dictionary (hash table) maassc. + maassc = dict() + ## First calculate maass coefficients corresponding to strictly + ## positive definite matrices: + from sage.rings.all import is_fundamental_discriminant + from sage.misc.all import isqrt + for disc in [d for d in xsrange(maxD, 0) if is_fundamental_discriminant(d)]: + for s in xsrange(1, isqrt(maxD//disc)+1): + ## add (disc*s^2,t) as a hash key, for each t that divides s + for t in s.divisors(): + maassc[(disc*s**2, t)] = sum([a**(k-1)*Cphi[disc*s**2/a**2] + for a in t.divisors()]) + + ## Compute the coefficients of the Siegel form $F$: + from sage.rings.all import gcd + siegelq = dict() + if isinstance(prec, tuple): + ## Note: m>=n>=r, n>=1 implies m>=n>r^2/4n + for r in xsrange(0, bmax): + for n in xsrange(max(r, 1), amax): + for m in xsrange(n, cmax): + D = r**2 - 4*m*n + g = gcd([n, r, m]) + siegelq[(n, r, m)] = maassc[(D, g)] + bound = cmax + else: + bound = 0 + for n in xsrange(1, isqrt(Dtop//3)+1): + for r in xsrange(n + 1): + bound = max(bound, (Dtop + r*r)//(4*n) + 1) + for m in xsrange(n, (Dtop + r*r)//(4*n) + 1): + D = r**2 - 4*m*n + g = gcd([n, r, m]) + siegelq[(n, r, m)] = maassc[(D, g)] + + ## Secondly, deal with the singular part. + ## Include the coeff corresponding to (0,0,0): + ## maassc = {(0,0): -bernoulli(k)/(2*k)*Cphi[0]} + from sage.rings.all import bernoulli + siegelq[(0, 0, 0)] = -bernoulli(k)/(2*k)*Cphi[0] + + ## Calculate the other discriminant-zero maass coefficients: + from sage.rings.all import sigma + for i in xsrange(1, bound): + ## maassc[(0,i)] = sigma(i, k-1) * Cphi[0] + siegelq[(0, 0, i)] = sigma(i, k-1) * Cphi[0] + + return _SiegelModularForm_from_dict(group=Sp4Z, weight=k, coeffs=siegelq, + prec=clean_prec, name=name) + + +def _SiegelModularForm_from_file(loc): + r""" + Initialize an instance of SiegelModularForm_class from a file. + + EXAMPLE:: + + sage: from sage.modular.siegel.siegel_modular_form import _SiegelModularForm_from_file + sage: url = 'http://sage.math.washington.edu/home/nils/Siegel-Modular-Forms/data/upsilon-forms_20-32_1000/Upsilon_20.p' # optional -- internet + sage: H = _SiegelModularForm_from_file(url); H # optional - internet + Upsilon_20 + sage: H[(4, 4, 4)] # optional - internet + 248256200704 + """ + from sage.modular.siegel.siegel_modular_group import Sp4Z + + if 'http://' == loc[:7]: + f = urllib.urlopen(loc) + else: + f = open(loc, 'r') + data = cPickle.load(f) + f.close() + fmt, wt, name, pol, prec, coeffs = data + if len(pol) > 2: + from sage.rings.all import PolynomialRing, NumberField + R = PolynomialRing(ZZ, 'x') + K = NumberField(R(pol), name='a') + for k in coeffs.iterkeys(): + coeffs[k] = K(coeffs[k]) + return _SiegelModularForm_from_dict(group=Sp4Z, weight=wt, coeffs=coeffs, + prec=prec, name=name) + + +def _SiegelModularForm_from_theta_characteristics(char_dict, prec=SMF_DEFAULT_PREC, name=None, hint=None): + r""" + Return the Siegel modular form + `\sum_{l \in char_dict} \alpha[l] * \prod_i \theta_l[i](8Z)`, + where `\theta_l[i]` denote the theta constant with characteristic `l`. + + INPUT: + + - char_dict -- a dictionary whose keys are *tuples* of theta characteristics + and whose values are in some ring. + - prec -- a precision for Siegel modular forms + - name -- a string describing this modular form + + EXAMPLES:: + + sage: theta_constants = {((1, 1, 0, 0), (0, 0, 1, 1), (1, 1, 0, 0), (0, 0, 1, 1)): 1} + sage: from sage.modular.siegel.siegel_modular_form import _SiegelModularForm_from_theta_characteristics + sage: _SiegelModularForm_from_theta_characteristics(theta_constants, prec=100) + Siegel modular form on None of weight 2 + sage: theta_constants = {((1, 1, 0, 0), (0, 0, 1, 1), (1, 1, 0, 0), (0, 0, 1, 1)): 1} + sage: S = _SiegelModularForm_from_theta_characteristics(theta_constants, prec=100) + sage: S.coeffs()[(2, 0, 10)] + 32 + + .. TODO:: + + Implement a parameter hint = "cusp_form" to prevent computing + singular parts + """ + coeffs = dict() + smf_prec = SiegelModularFormPrecision(prec) + for f in smf_prec: + if hint is 'cusp_form': + a, b, c = f + if 0 == a: + continue + from theta_constant import _compute_theta_char_poly + val = _compute_theta_char_poly(char_dict, f) + if val != 0: + coeffs[f] = val + wt = 0 + for l in char_dict: + wt = max(wt, len(l)/2) + return _SiegelModularForm_from_dict(group=None, weight=wt, + coeffs=coeffs, prec=prec, name=name) + + +def _SiegelModularForm_from_QuadraticForm(Q, prec, name): + """ + Return the theta series of degree 2 for the quadratic form Q. + + INPUT: + + - ``Q`` - a quadratic form. + + - ``prec`` - an integer. + + EXAMPLE:: + + sage: from sage.modular.siegel.siegel_modular_form import _SiegelModularForm_from_QuadraticForm + sage: Q11 = QuadraticForm(ZZ, 4, [3,0,11,0, 3,0,11, 11,0, 11]) + sage: _SiegelModularForm_from_QuadraticForm(Q11, 100, None) + Siegel modular form on Siegel Modular Group Gamma0(11) of weight 2 + """ + N = Q.level() + k = Q.dim()/2 + coeffs = Q.theta_series_degree_2(prec) + return _SiegelModularForm_from_dict(group=Sp4Z_Gamma0(N), weight=k, + coeffs=coeffs, prec=prec, name=name) + + +def _SiegelModularForm_borcherds_lift(f, prec=SMF_DEFAULT_PREC, name=None): + r""" + Returns Borcherds lift. + + EXAMPLES:: + + sage: from sage.modular.siegel.siegel_modular_form import _SiegelModularForm_borcherds_lift + sage: _SiegelModularForm_borcherds_lift(0) + Traceback (most recent call last): + ... + NotImplementedError: Borcherds lift not yet implemented + """ + raise NotImplementedError("Borcherds lift not yet implemented") + + +def _SiegelModularForm_yoshida_lift(f, g, prec=SMF_DEFAULT_PREC, name=None): + r""" + Returns the Yoshida lift. + + EXAMPLES:: + + sage: from sage.modular.siegel.siegel_modular_form import _SiegelModularForm_yoshida_lift + sage: _SiegelModularForm_yoshida_lift(0, 0) + Traceback (most recent call last): + ... + NotImplementedError: Yoshida lift not yet implemented + """ + raise NotImplementedError('Yoshida lift not yet implemented') + + +def _SiegelModularForm_from_weil_representation(gram, prec=SMF_DEFAULT_PREC, name=None): + r""" + Returns a form associated to a Weil representation. + + EXAMPLES:: + + sage: from sage.modular.siegel.siegel_modular_form import _SiegelModularForm_from_weil_representation + sage: _SiegelModularForm_from_weil_representation(matrix([[2, 1], [1, 1]])) + Traceback (most recent call last): + ... + NotImplementedError: Weil representation argument not yet implemented + """ + raise NotImplementedError('Weil representation argument not yet implemented') + + +def _SiegelModularForm_singular_weight(gram, prec=SMF_DEFAULT_PREC, name=None): + r""" + Returns a singular modular form from gram. + + EXAMPLES:: + + sage: from sage.modular.siegel.siegel_modular_form import _SiegelModularForm_singular_weight + sage: _SiegelModularForm_singular_weight(matrix([[2, 1], [1, 1]])) + Traceback (most recent call last): + ... + NotImplementedError: singular form argument not yet implemented + """ + raise NotImplementedError('singular form argument not yet implemented') + + +def _SiegelModularForm_from_dict(group, weight, coeffs, prec, degree=2, parent=None, name=None): + r""" + Create an instance of SiegelModularForm_class(), where the parent + is computed from the coeffs. + + EXAMPLES:: + + sage: d = {(1, 1, 1): 1, (1, 0, 2): 2} + sage: from sage.modular.siegel.siegel_modular_form import _SiegelModularForm_from_dict + sage: S = _SiegelModularForm_from_dict(Sp4Z, 2, d, 1); S + Siegel modular form on Siegel Modular Group Sp(4,Z) of weight 2 + sage: S.coeffs() + {} + sage: S = _SiegelModularForm_from_dict(Sp4Z, 2, d, 4); S + Siegel modular form on Siegel Modular Group Sp(4,Z) of weight 2 + sage: S.coeffs() + {(1, 1, 1): 1} + sage: S = _SiegelModularForm_from_dict(Sp4Z, 2, d, 9); S + Siegel modular form on Siegel Modular Group Sp(4,Z) of weight 2 + sage: S.coeffs() + {(1, 1, 1): 1, (1, 0, 2): 2} + """ + if prec is None: + prec = -min([b**2 - 4*a*c for (a, b, c) in coeffs]) + 1 + + if parent is None: + from sage.structure.all import Sequence + from siegel_modular_forms_algebra import SiegelModularFormsAlgebra + s = Sequence(coeffs.values()) + if weight % 2: + weights = 'all' + else: + weights = 'even' + if len(s) == 0: + coeff_ring = ZZ + else: + coeff_ring = s.universe() + parent = SiegelModularFormsAlgebra(coeff_ring=coeff_ring, group=group, weights=weights, degree=degree) + + smf_prec = SiegelModularFormPrecision(prec) + coeffs_up_to_prec = {} + for k in coeffs: + if smf_prec.is_in_bound(k): + coeffs_up_to_prec[k] = coeff_ring(coeffs[k]) + + return parent.element_class(parent=parent, weight=weight, coeffs=coeffs_up_to_prec, prec=prec, name=name) + + +def _normalized_prec(prec): + r""" + Return a normalized prec for instance of class SiegelModularForm_class. + + EXAMPLES:: + + sage: from sage.modular.siegel.siegel_modular_form import _normalized_prec + sage: _normalized_prec((5, 4, 5)) + (5, 4, 5) + sage: _normalized_prec(101) + 101 + sage: _normalized_prec(infinity) + +Infinity + + .. NOTE:: + + $(a,b,c)$ is within the precison defined + by prec, if $a,b,c < floor(a),floor(b),floor(c)% + respectively $4ac-b^2 < floor(prec)$. + """ + from sage.rings.all import infinity + if prec is infinity: + return prec + from sage.functions.all import floor + if isinstance(prec, tuple): + a, b, c = prec + a = min(floor(a), floor(c)) + b = min(floor(b), a) + if b <= 0: + return (0, 0, 0) + return a, b, c + D = max(0, floor(prec)) + while 0 != D % 4 and 1 != D % 4: + D -= 1 + return D + + +def _is_bounded((a, b, c), prec): + r""" + Return True or False accordingly as (a, b, c) is in the + region defined by prec (see. SiegelModularForm_class for what + this means). + + EXAMPLES:: + + sage: from sage.modular.siegel.siegel_modular_form import _is_bounded + sage: _is_bounded((2, 0, 2), 15) + False + sage: _is_bounded((2, 0, 2), 16) + False + sage: _is_bounded((2, 0, 2), 17) + True + + .. NOTE:: + + It is assumed that prec is normalized accordingly to the output + of _normalized_prec(). + """ + if isinstance(prec, tuple): + return (a, b, c) < prec + else: + D = 4*a*c - b*b + return D < prec if D != 0 else c < (prec+1)/4 + + +def _prec_min(prec1, prec2): + r""" + Returns the min of two precs. + + EXAMPLES:: + + sage: from sage.modular.siegel.siegel_modular_form import _prec_min + sage: _prec_min(100, 200) + 100 + sage: _prec_min((1, 0, 1), (2, 0, 1)) + (1, 0, 1) + + .. NOTE:: + + It is assumed that the arguments are normalized according to the output + of _normalized_prec(). + """ + if type(prec1) != type(prec2): + raise NotImplementedError("addition for differing precs not " + "yet implemented") + return prec1 if prec1 <= prec2 else prec2 diff --git a/src/sage/modular/siegel/siegel_modular_form_prec.py b/src/sage/modular/siegel/siegel_modular_form_prec.py new file mode 100644 index 00000000000..d2a5ae7bb77 --- /dev/null +++ b/src/sage/modular/siegel/siegel_modular_form_prec.py @@ -0,0 +1,502 @@ +""" +Precision for Fourier expansions of Siegel modular forms +""" +from copy import deepcopy + +from sage.rings.integer import Integer +from sage.functions.other import ceil +from sage.misc.functional import isqrt +from sage.structure.sage_object import SageObject +from sage.rings.all import infinity +from sage.misc.latex import latex + +#import operator + +#from sage.functions.other import (floor, ceil) +#from fastmult import reduce_GL + +#load 'fastmult.spyx' +#from fastmult import reduce_GL + + +class SiegelModularFormPrecision (SageObject): + r""" + Stores information on the precision of the Fourier expansion of a Siegel + modular form and provides checking. + """ + def __init__(self, prec): + r""" + INPUT: + + - prec -- an integer or infinity or a triple (aprec, bprec, cprec) of integers. + or an instance of class SiegelModularFormPrecision + + .. NOTE:: + + We distinguish two types of precision bounds: + 'disc' or 'box'. prec indicates that all terms q^(a, b, c) + with GL(2, Z)-reduced (a, b, c) are available such that + either 4ac-b^2 < prec ('disc'-precision) or else + one has componentwise + (a, b, c) < ( aprec, bprec, cprec) + ('box'-precision). + + EXAMPLES:: + + sage: from sage.modular.siegel.siegel_modular_form_prec import SiegelModularFormPrecision + sage: prec = SiegelModularFormPrecision(101) + sage: prec2 = SiegelModularFormPrecision(prec) + sage: prec + Discriminant precision for Siegel modular form with bound 101 + sage: prec2 + Discriminant precision for Siegel modular form with bound 101 + sage: prec == prec2 + True + sage: prec3 = SiegelModularFormPrecision((5, 5, 5)) + sage: prec3 + Box precision for Siegel modular form with bound (5, 5, 5) + sage: prec4 = SiegelModularFormPrecision(infinity) + sage: prec4 + Infinity precision for Siegel modular form with bound +Infinity + sage: prec5 = SiegelModularFormPrecision((5, 0, 5)) + sage: prec5.prec() + (0, 0, 0) + + TESTS:: + + sage: prec = SiegelModularFormPrecision(101) + sage: TestSuite(prec).run() + """ + ## copy constructor + if isinstance(prec, SiegelModularFormPrecision): + self.__type = deepcopy(prec.type()) + self.__prec = deepcopy(prec.prec()) + + elif prec is infinity: + self.__type = 'infinity' + self.__prec = infinity + + elif isinstance(prec, (int, Integer)): + self.__type = 'disc' + prec = max(0, prec) + Dmod = prec % 4 + if Dmod >= 2: + self.__prec = Integer(prec - Dmod + 1) + else: + self.__prec = Integer(prec) + + elif isinstance(prec, tuple) and len(prec) == 3 and all([isinstance(e, (int, Integer)) for e in list(prec)]): + self.__type = 'box' + a, b, c = prec + a = min(a, c) + b = min(b, a) + if b <= 0: + self.__prec = (0, 0, 0) + else: + self.__prec = (Integer(a), Integer(b), Integer(c)) + + else: + raise TypeError("incorrect type %s for prec" % type(prec)) + + def _repr_(self): + r""" + Return the repr of ``self``. + + EXAMPLES:: + + sage: from sage.modular.siegel.siegel_modular_form_prec import SiegelModularFormPrecision + sage: prec = SiegelModularFormPrecision(101) + sage: prec._repr_() + 'Discriminant precision for Siegel modular form with bound 101' + """ + if self.type() == 'infinity': + n = "Infinity" + elif self.type() == 'disc': + n = "Discriminant" + elif self.type() == 'box': + n = "Box" + else: + raise RuntimeError("Unexpected value of self.__type") + + return "%s precision for Siegel modular form with bound %s" % \ + (n, repr(self.__prec)) + + def _latex_(self): + r""" + Return the latex of ``self``. + + EXAMPLES:: + + sage: from sage.modular.siegel.siegel_modular_form_prec import SiegelModularFormPrecision + sage: prec = SiegelModularFormPrecision(101) + sage: prec._latex_() + 'Discriminant precision for Siegel modular form with bound $101$' + + """ + if self.__type == 'infinity': + n = "Infinity" + elif self.__type == 'disc': + n = "Discriminant" + elif self.__type == 'box': + n = "Box" + else: + raise ValueError("Unexpected value of self.__type") + + return "%s precision for Siegel modular form with bound $%s$" % (str(n), str(latex(self.__prec))) + + def type(self): + r""" + Return the type of self: "box", "disc" or "infinity". + + EXAMPLES:: + + sage: from sage.modular.siegel.siegel_modular_form_prec import SiegelModularFormPrecision + sage: prec = SiegelModularFormPrecision((2, 2, 2)) + sage: prec = SiegelModularFormPrecision(infinity) + sage: prec.type() + 'infinity' + sage: prec = SiegelModularFormPrecision(101) + sage: prec.type() + 'disc' + """ + return self.__type + + def prec(self): + r""" + Return the tuple of the maximal box, the integer for the + maximal discriminant and infinity otherwise. + + If ``self`` is of + type "disc" then it will return the largest fundamental + discriminant just below. + + EXAMPLES:: + + sage: from sage.modular.siegel.siegel_modular_form_prec import SiegelModularFormPrecision + sage: prec = SiegelModularFormPrecision(11) + sage: prec.prec() + 9 + sage: prec = SiegelModularFormPrecision(8) + sage: prec.prec() + 8 + """ + return self.__prec + + def is_in_bound(self, t): + r""" + Return ``True`` or ``False`` accordingly as the GL(2, Z)-reduction of + the tuple (a, b, c)=t is within the bounds of this precision. + + EXAMPLES:: + + sage: from sage.modular.siegel.siegel_modular_form_prec import SiegelModularFormPrecision + sage: prec = SiegelModularFormPrecision(101) + sage: prec.is_in_bound((1, 0, 1)) + True + sage: prec.is_in_bound((0, 0, 25)) + True + sage: prec.is_in_bound((0, 0, 26)) + False + sage: prec.is_in_bound((6, 0, 6)) + False + sage: prec = SiegelModularFormPrecision((2, 2, 2)) + sage: prec.is_in_bound((1, 0, 1)) + True + sage: prec.is_in_bound((2, 2, 2)) + False + + .. TODO:: + + fix this + + .. NOTE:: + + It is assumed that (a, b, c) is semi positive definite + """ + (a, b, c) = t + + from fastmult import reduce_GL + if self.__type == 'infinity': + return True + elif self.__type == 'disc': + (a, b, c) = reduce_GL(a, b, c) + if a == 0 and b == 0: + return c < self.get_contents_bound_for_semi_definite_forms() + D = 4 * a * c - b ** 2 + if D < 0: + return False + if D > 0: + return D < self.__prec + elif self.__type == 'box': + (a, b, c) = reduce_GL(a, b, c) + return a < self.__prec[0] and b < self.__prec[1] and c < self.__prec[2] + else: + raise RuntimeError("Unexpected value of self.__type") + + def get_contents_bound_for_semi_definite_forms(self): + r""" + Return the bound for the contents of a semi positive definite + form falling within this precision. + + EXAMPLES:: + + sage: from sage.modular.siegel.siegel_modular_form_prec import SiegelModularFormPrecision + sage: prec = SiegelModularFormPrecision(7) + sage: prec.get_contents_bound_for_semi_definite_forms() + 2 + sage: prec = SiegelModularFormPrecision(101) + sage: prec.get_contents_bound_for_semi_definite_forms() + 26 + + .. NOTE:: + + If (a, b, c) is semi positive definite, then it is + GL(2, Z)-equivalent to a the form (0, 0, g) where g is the + contents (gcd) of a, b, c. + + I don't know what this does. It seems to only do it for singular forms-- NR + """ + if self.__type == 'infinity': + return infinity + elif self.__type == 'disc': + return ceil((self.__prec+1)/4) + elif self.__type == 'box': + return self.__prec[2] + else: + raise RuntimeError("Unexpected value of self.__type") + + def __lt__(self, other): + r""" + Return ``True`` if the set of GL(2, Z)-reduced forms within the + precision ``self`` is contained in but not equal to the + corresponding set of forms within the precision ``other``. + + EXAMPLES:: + + sage: from sage.modular.siegel.siegel_modular_form_prec import SiegelModularFormPrecision + sage: prec1 = SiegelModularFormPrecision(101) + sage: prec2 = SiegelModularFormPrecision(100) + sage: prec3 = SiegelModularFormPrecision(201) + sage: prec1 + Discriminant precision for Siegel modular form with bound 101 + sage: prec2 + Discriminant precision for Siegel modular form with bound 100 + sage: prec2 < prec1 + True + sage: prec3 < prec1 + False + sage: prec1.__lt__(prec3) + True + """ + if not isinstance(other, SiegelModularFormPrecision): + raise NotImplementedError("can only compare with elements of " + "the same class") + + ## TODO: implement comparison of precisions of different types + if self.__type != other.__type: + return False + + if self.__type == 'box': + a, b, c = self.__prec + ao, bo, co = other.__prec + return (a <= ao and b <= bo and c <= co + and any([a < ao, b < bo, c < co])) + + elif self.__type == 'disc': + return self.__prec < other.__prec + + else: + raise RuntimeError("Unexpected value of self.__type") + + def __le__(self, other): + r""" + Return whether ``self <= other``. + + EXAMPLES:: + + sage: from sage.modular.siegel.siegel_modular_form_prec import SiegelModularFormPrecision + sage: prec1 = SiegelModularFormPrecision(101) + sage: prec2 = SiegelModularFormPrecision(100) + sage: prec3 = SiegelModularFormPrecision(101) + sage: prec1 <= prec3 + True + sage: prec1 <= prec2 + False + sage: prec1 <= prec1 + True + sage: prec4 = SiegelModularFormPrecision(infinity) + sage: prec5 = SiegelModularFormPrecision((5, 5, 5)) + sage: prec5.__le__(prec4) + False + + .. TODO:: + + It seems everything should be <= infinity + """ + return self.__lt__(other) or self.__eq__(other) + + def __eq__(self, other): + r""" + Return self == other as defined by having the same prec and type. + + EXAMPLES:: + + sage: from sage.modular.siegel.siegel_modular_form_prec import SiegelModularFormPrecision + sage: prec1 = SiegelModularFormPrecision(101) + sage: prec2 = SiegelModularFormPrecision(100) + sage: prec3 = SiegelModularFormPrecision(101) + sage: prec4 = SiegelModularFormPrecision(infinity) + sage: prec5 = SiegelModularFormPrecision((5, 5, 5)) + sage: prec1 == prec2 + False + sage: prec1 == prec3 + True + sage: prec4 == prec2 + False + sage: prec4 == prec4 + True + sage: prec5.__eq__(prec5) + True + """ + if not isinstance(other, SiegelModularFormPrecision): + raise NotImplementedError("can only compare with elements of the same class") + + return (self.__type == other.type() and self.__prec == other.prec()) + + def __ne__(self, other): + r""" + Return self != other. + + EXAMPLES:: + + sage: from sage.modular.siegel.siegel_modular_form_prec import SiegelModularFormPrecision + sage: prec1 = SiegelModularFormPrecision(101) + sage: prec2 = SiegelModularFormPrecision(100) + sage: prec3 = SiegelModularFormPrecision(101) + sage: prec1.__ne__(prec1) + False + sage: prec1 != prec2 + True + sage: prec1 != prec3 + False + """ + return not self.__eq__(other) + + def __gt__(self, other): + r""" + Return self > other. + + EXAMPLES:: + + sage: from sage.modular.siegel.siegel_modular_form_prec import SiegelModularFormPrecision + sage: prec = SiegelModularFormPrecision(101) + sage: prec2 = SiegelModularFormPrecision(100) + sage: prec > prec2 + True + sage: prec2.__gt__(prec) + False + """ + return other.__lt__(self) + + def __ge__(self, other): + r""" + Return self >= other. + + EXAMPLES:: + + sage: from sage.modular.siegel.siegel_modular_form_prec import SiegelModularFormPrecision + sage: prec = SiegelModularFormPrecision(101) + sage: prec2 = SiegelModularFormPrecision(100) + sage: prec >= prec2 + True + sage: prec.__ge__(prec) + True + """ + return self.__eq__(other) or other.__lt__(self) + + def __iter__(self): + r""" + Iterate over all GL(2, Z)-reduced semi positive forms + which are within the bounds of this precision. + + EXAMPLES:: + + sage: from sage.modular.siegel.siegel_modular_form_prec import SiegelModularFormPrecision + sage: prec = SiegelModularFormPrecision(11) + sage: for k in prec.__iter__(): print(k) + (0, 0, 0) + (0, 0, 1) + (0, 0, 2) + (1, 0, 1) + (1, 0, 2) + (1, 1, 1) + (1, 1, 2) + + .. NOTE:: + + The forms are enumerated in lexicographic order. + """ + if self.__type == 'disc': + bound = self.get_contents_bound_for_semi_definite_forms() + for c in xrange(0, bound): + yield (0, 0, c) + + atop = isqrt(self.__prec // 3) + if 3*atop*atop == self.__prec: + atop -= 1 + for a in xrange(1, atop + 1): + for b in xrange(a+1): + for c in xrange(a, ceil((b**2 + self.__prec)/(4*a))): + yield (a, b, c) + + elif 'box' == self.__type: + (am, bm, cm) = self.__prec + for a in xrange(am): + for b in xrange(min(bm, a+1)): + for c in xrange(a, cm): + yield (a, b, c) + + else: + raise RuntimeError("Unexpected value of self.__type") + + raise StopIteration + + def positive_forms(self): + r""" + Iterate over all GL(2, Z)-reduced strictly positive forms + which are within the bounds of this precision. + + EXAMPLES:: + + sage: from sage.modular.siegel.siegel_modular_form_prec import SiegelModularFormPrecision + sage: prec = SiegelModularFormPrecision(11) + sage: for k in prec.positive_forms(): print(k) + (1, 0, 1) + (1, 0, 2) + (1, 1, 1) + (1, 1, 2) + + .. NOTE:: + + The forms are enumerated in lexicographic order. + """ + if self.__type == 'disc': + atop = isqrt(self.__prec // 3) + if 3*atop*atop == self.__prec: + atop -= 1 + for a in xrange(1, atop + 1): + for b in xrange(a+1): + for c in xrange(a, ceil((b**2 + self.__prec)/(4*a))): + yield (a, b, c) + + elif 'box' == self.__type: + (am, bm, cm) = self.__prec + for a in xrange(am): + for b in xrange(min(bm, a+1)): + for c in xrange(a, cm): + yield (a, b, c) + + else: + raise RuntimeError("Unexpected value of self.__type") + + raise StopIteration diff --git a/src/sage/modular/siegel/siegel_modular_forms_algebra.py b/src/sage/modular/siegel/siegel_modular_forms_algebra.py new file mode 100644 index 00000000000..717085c16b7 --- /dev/null +++ b/src/sage/modular/siegel/siegel_modular_forms_algebra.py @@ -0,0 +1,517 @@ +from siegel_modular_form import (SiegelModularForm, SiegelModularForm_class, + SMF_DEFAULT_PREC) +from siegel_modular_group import Sp4Z +from sage.algebras.algebra import Algebra + +from sage.misc.all import cached_method, cached_function +from sage.rings.all import ZZ, QQ +from sage.structure.factory import UniqueFactory +import sage.groups.group as grp + + +class SiegelModularFormsAlgebraFactory(UniqueFactory): + """ + A factory for creating algebras of Siegel modular forms. It + handles making sure that they are unique as well as handling + pickling. For more details, see + :class:`~sage.structure.factory.UniqueFactory` and + :mod:`~sage.modular.siegel.siegel_modular_forms_algebra`. + + EXAMPLES:: + + sage: S = SiegelModularFormsAlgebra() + sage: S.construction() + (SMFAlg{"Siegel Modular Group Sp(4,Z)", "even", 2, 101}, Integer Ring) + sage: R. = QQ[] + sage: S = SiegelModularFormsAlgebra(coeff_ring=R) + sage: S.construction() + (SMFAlg{"Siegel Modular Group Sp(4,Z)", "even", 2, 101}, Multivariate Polynomial Ring in a, b over Rational Field) + sage: S is loads(dumps(S)) + True + sage: TestSuite(S).run() + """ + def create_key(self, coeff_ring=ZZ, group=Sp4Z, weights='even', degree=2, default_prec=SMF_DEFAULT_PREC): + """ + Create a key which uniquely defines this algebra of Siegel modular + forms. + + TESTS:: + + sage: SiegelModularFormsAlgebra.create_key(coeff_ring=QQ, group=Sp4Z_Gamma0(5), weights='all', degree=2, default_prec=41) + (SMFAlg{"Siegel Modular Group Gamma0(5)", "all", 2, 41}(FractionField(...)), Integer Ring) + """ + from sage.rings.ring import is_Ring + if not is_Ring(coeff_ring): + raise TypeError('The coefficient ring must be a ring') + if not (isinstance(group, grp.Group) or group is None): + raise TypeError('The group must be given by an GSp4_Arithmetic_Subgroup, or None') + if not weights in ('even', 'all'): + raise ValueError("The argument weights must be 'even' or 'all'") + try: + degree = int(degree) + except TypeError: + raise TypeError('The degree must be a positive integer') + if degree < 1: + raise ValueError('The degree must be a positive integer') + if degree == 1: + raise ValueError('Use ModularForms if you want to work with Siegel modular forms of degree 1') + if degree > 2: + raise NotImplementedError('Siegel modular forms of degree > 2 are not yet implemented') + try: + default_prec = int(default_prec) + except TypeError: + raise TypeError('The default precision must be a positive integer') + + from sage.categories.pushout import SiegelModularFormsAlgebraFunctor + F = SiegelModularFormsAlgebraFunctor(group=group, weights=weights, degree=degree, default_prec=default_prec) + + while hasattr(coeff_ring, 'construction'): + C = coeff_ring.construction() + if C is None: + break + F = F * C[0] + coeff_ring = C[1] + return (F, coeff_ring) + + def create_object(self, version, key): + """ + Return the algebra of Siegel modular forms corresponding to the + key ``key``. + + TESTS:: + + sage: key = SiegelModularFormsAlgebra.create_key(coeff_ring=QQ, group=Sp4Z_Gamma0(5), weights='all', degree=2, default_prec=41) + sage: SiegelModularFormsAlgebra.create_object('1.0', key) + Algebra of Siegel modular forms of degree 2 and all weights on Siegel Modular Group Gamma0(5) over Rational Field + + """ + C, R = key + from sage.categories.pushout import CompositeConstructionFunctor, SiegelModularFormsAlgebraFunctor + if isinstance(C, CompositeConstructionFunctor): + F = C.all[-1] + if len(C.all) > 1: + R = CompositeConstructionFunctor(*C.all[:-1])(R) + else: + F = C + if not isinstance(F, SiegelModularFormsAlgebraFunctor): + raise TypeError("We expected a SiegelModularFormsAlgebraFunctor, not %s" % type(F)) + return SiegelModularFormsAlgebra_class(coeff_ring=R, group=F._group, weights=F._weights, degree=F._degree, default_prec=F._default_prec) + +SiegelModularFormsAlgebra = SiegelModularFormsAlgebraFactory('SiegelModularFormsAlgebra') + + +class SiegelModularFormsAlgebra_class(Algebra): + def __init__(self, coeff_ring=ZZ, group=Sp4Z, weights='even', degree=2, default_prec=SMF_DEFAULT_PREC): + r""" + Initialize an algebra of Siegel modular forms of degree ``degree`` + with coefficients in ``coeff_ring``, on the group ``group``. + If ``weights`` is 'even', then only forms of even weights are + considered; if ``weights`` is 'all', then all forms are + considered. + + The parameter ``default_prec`` is the precision to which every + Siegel modular form in `self` is known. + + EXAMPLES:: + + sage: A = SiegelModularFormsAlgebra(QQ) + sage: B = SiegelModularFormsAlgebra(ZZ) + sage: A._coerce_map_from_(B) + True + sage: B._coerce_map_from_(A) + False + sage: A._coerce_map_from_(ZZ) + True + """ + self.__coeff_ring = coeff_ring + self.__group = group + self.__weights = weights + self.__degree = degree + self.__default_prec = default_prec + R = coeff_ring + from sage.algebras.all import GroupAlgebra + if isinstance(R, GroupAlgebra): + R = R.base_ring() + from sage.rings.polynomial.polynomial_ring import is_PolynomialRing + if is_PolynomialRing(R): + self.__base_ring = R.base_ring() + else: + self.__base_ring = R + Algebra.__init__(self, self.__base_ring) + + @cached_method + def gens(self, prec=None): + """ + Return the ring generators of this algebra. + + EXAMPLES:: + + sage: S = SiegelModularFormsAlgebra() + sage: S.gens() + [Igusa_4, Igusa_6, Igusa_10, Igusa_12] + + TESTS:: + + sage: S = SiegelModularFormsAlgebra(group=Sp4Z_Gamma0(2)) + sage: S.gens() + Traceback (most recent call last): + ... + NotImplementedError: Not yet implemented + """ + if prec is None: + prec = self.default_prec() + return _siegel_modular_forms_generators(self, prec=prec) + + def is_commutative(self): + r""" + Return ``True`` since all our algebras are commutative. + + EXAMPLES:: + + sage: S = SiegelModularFormsAlgebra(coeff_ring=QQ) + sage: S.is_commutative() + True + """ + return True + + @cached_method + def ngens(self): + r""" + Return the number of generators of this algebra of Siegel modular + forms. + + EXAMPLES:: + + sage: S = SiegelModularFormsAlgebra(coeff_ring=QQ) + sage: S.ngens() + 4 + """ + return len(self.gens()) + + @cached_method + def gen(self, i, prec=None): + r""" + Return the `i`-th generator of this algebra. + + EXAMPLES:: + + sage: S = SiegelModularFormsAlgebra(coeff_ring=QQ) + sage: S.ngens() + 4 + sage: S.gen(2) + Igusa_10 + """ + return self.gens(prec=prec)[i] + + def coeff_ring(self): + r""" + Return the ring of coefficients of this algebra. + + EXAMPLES:: + + sage: S = SiegelModularFormsAlgebra(coeff_ring=QQ) + sage: S.coeff_ring() + Rational Field + """ + return self.__coeff_ring + + def group(self): + """ + Return the modular group of this algebra. + + EXAMPLES:: + + sage: S = SiegelModularFormsAlgebra(group=Sp4Z_Gamma0(7)) + sage: S.group() + Siegel Modular Group Gamma0(7) + """ + return self.__group + + def weights(self): + """ + Return ``'even'`` or ``'all'`` depending on whether this algebra + contains only even weight forms, or forms of all weights. + + EXAMPLES:: + + sage: S = SiegelModularFormsAlgebra(weights='even') + sage: S.weights() + 'even' + sage: S = SiegelModularFormsAlgebra(weights='all') + sage: S.weights() + 'all' + """ + return self.__weights + + def degree(self): + """ + Return the degree of the modular forms in this algebra. + + EXAMPLES:: + + sage: S = SiegelModularFormsAlgebra() + sage: S.degree() + 2 + """ + return self.__degree + + def default_prec(self): + """ + Return the default precision for the Fourier expansions of + modular forms in this algebra. + + EXAMPLES:: + + sage: S = SiegelModularFormsAlgebra(default_prec=41) + sage: S.default_prec() + 41 + """ + return self.__default_prec + + def _repr_(self): + r""" + Return the plain text representation of this algebra. + + EXAMPLES:: + + sage: S = SiegelModularFormsAlgebra(coeff_ring=QQ) + sage: S._repr_() + 'Algebra of Siegel modular forms of degree 2 and even weights on Siegel Modular Group Sp(4,Z) over Rational Field' + """ + return 'Algebra of Siegel modular forms of degree %s and %s weights on %s over %s' % (self.__degree, self.__weights, self.__group, self.__coeff_ring) + + def _latex_(self): + r""" + Return the latex representation of this algebra. + + EXAMPLES:: + + sage: S = SiegelModularFormsAlgebra(coeff_ring=QQ) + sage: S._latex_() + '\\texttt{Algebra of Siegel modular forms of degree }2\\texttt{ and even weights on Siegel Modular Group Sp(4,Z) over }\\Bold{Q}' + + """ + from sage.misc.all import latex + return r'\texttt{Algebra of Siegel modular forms of degree }%s\texttt{ and %s weights on %s over }%s' % (latex(self.__degree), self.__weights, self.__group, latex(self.__coeff_ring)) + + def _coerce_map_from_(self, other): + r""" + Return ``True`` if it is possible to coerce from ``other`` into the + algebra of Siegel modular forms ``self``. + + EXAMPLES:: + + sage: S = SiegelModularFormsAlgebra(coeff_ring=QQ) + sage: R = SiegelModularFormsAlgebra(coeff_ring=ZZ) + sage: S._coerce_map_from_(R) + True + sage: R._coerce_map_from_(S) + False + sage: S._coerce_map_from_(ZZ) + True + """ + if self.__coeff_ring is other or self.base_ring().has_coerce_map_from(other): + return True + if isinstance(other, SiegelModularFormsAlgebra_class): + if (self.group() == other.group()) or (other.group() == Sp4Z): + if self.coeff_ring().has_coerce_map_from(other.coeff_ring()): + if self.degree() == other.degree(): + return True + return False + + def _element_constructor_(self, x): + r""" + Return the element of the algebra ``self`` corresponding to ``x``. + + EXAMPLES:: + + sage: S = SiegelModularFormsAlgebra(coeff_ring=QQ) + sage: B = SiegelModularFormsAlgebra(coeff_ring=ZZ).1 + sage: S(B) + Igusa_6 + sage: S(1/5) + Siegel modular form on Siegel Modular Group Sp(4,Z) of weight 0 + sage: S(1/5).parent() is S + True + sage: S._element_constructor_(2.67) + Traceback (most recent call last): + ... + TypeError: Unable to construct an element of Algebra of Siegel modular forms of degree 2 and even weights on Siegel Modular Group Sp(4,Z) over Rational Field corresponding to 2.67000000000000 + + sage: S.base_extend(RR)._element_constructor_(2.67) + 2.67000000000000 + """ + if isinstance(x, (int, long)): + x = ZZ(x) + if isinstance(x, float): + from sage.rings.all import RR + x = RR(x) + if isinstance(x, complex): + from sage.rings.all import CC + x = CC(x) + + if isinstance(x.parent(), SiegelModularFormsAlgebra_class): + d = dict((f, self.coeff_ring()(x[f])) for f in x.coeffs()) + return self.element_class(parent=self, weight=x.weight(), coeffs=d, prec=x.prec(), name=x.name()) + + R = self.base_ring() + if R.has_coerce_map_from(x.parent()): + d = {(0, 0, 0): R(x)} + from sage.rings.all import infinity + return self.element_class(parent=self, weight=0, coeffs=d, prec=infinity, name=str(x)) + elif x.parent() is self.coeff_ring(): + from sage.rings.all import infinity + return self.element_class(parent=self, weight=0, + coeffs={(0, 0, 0): x}, prec=infinity) + else: + raise TypeError("Unable to construct an element of %s corresponding to %s" % (self, x)) + + Element = SiegelModularForm_class + + def _an_element_(self): + r""" + Return an element of the algebra ``self``. + + EXAMPLES:: + + sage: S = SiegelModularFormsAlgebra(coeff_ring=QQ) + sage: z = S._an_element_() + sage: z in S + True + """ + return self(0) + + def base_extend(self, R): + r""" + Extends the base ring of the algebra ``self`` to ``R``. + + EXAMPLES:: + + sage: S = SiegelModularFormsAlgebra(coeff_ring=QQ) + sage: S.base_extend(RR) + Algebra of Siegel modular forms of degree 2 and even weights on Siegel Modular Group Sp(4,Z) over Real Field with 53 bits of precision + + """ + #B = self.base_ring() + S = self.coeff_ring() + from sage.rings.polynomial.polynomial_ring import is_PolynomialRing + if is_PolynomialRing(S): + xS = S.base_extend(R) + elif R.has_coerce_map_from(S): + xS = R + else: + raise TypeError("cannot extend to %s" % R) + return SiegelModularFormsAlgebra(coeff_ring=xS, group=self.group(), weights=self.weights(), degree=self.degree(), default_prec=self.default_prec()) + + def construction(self): + r""" + Return the construction of the algebra ``self``. + + EXAMPLES:: + + sage: S = SiegelModularFormsAlgebra(coeff_ring=QQ) + sage: S.construction() + (SMFAlg{"Siegel Modular Group Sp(4,Z)", "even", 2, 101}, Rational Field) + """ + from sage.categories.pushout import SiegelModularFormsAlgebraFunctor + return SiegelModularFormsAlgebraFunctor(self.group(), self.weights(), self.degree(), self.default_prec()), self.coeff_ring() + + +def _siegel_modular_forms_generators(parent, prec=None, coefficient_degree=0): + r""" + Compute the four Igusa generators of the ring of Siegel modular forms + of degree 2 and level 1 and even weight (this happens if weights='even'). + If weights = 'all' you get the Siegel modular forms of degree 2, level 1 + and even and odd weight. + + The parameter ``coefficient_degree`` is 0 if the Siegel modular form + is scalar-valued and otherwise is the degree of the coefficients + of the vector-valued Siegel modular form when coefficients are + considered as being homogeneous polynomials. + + EXAMPLES:: + + sage: A, B, C, D = SiegelModularFormsAlgebra().gens() + sage: C[(2, 0, 9)] + -390420 + sage: D[(4, 2, 5)] + 17689760 + sage: A2, B2, C2, D2 = SiegelModularFormsAlgebra().gens(prec=50) + sage: A[(1, 0, 1)] == A2[(1, 0, 1)] + True + sage: A3, B3, C3, D3 = SiegelModularFormsAlgebra().gens(prec=500) + sage: B2[(2, 1, 3)] == B3[(2, 1, 3)] + True + + TESTS:: + + sage: from sage.modular.siegel.siegel_modular_forms_algebra import _siegel_modular_forms_generators + sage: S = SiegelModularFormsAlgebra() + sage: S.gens() == _siegel_modular_forms_generators(S) + True + """ + group = parent.group() + weights = parent.weights() + if prec is None: + prec = parent.default_prec() + if group == Sp4Z and 0 == coefficient_degree: + from sage.modular.all import ModularForms + E4 = ModularForms(1, 4).gen(0) + M6 = ModularForms(1, 6) + E6 = ModularForms(1, 6).gen(0) + M8 = ModularForms(1, 8) + M10 = ModularForms(1, 10) + Delta = ModularForms(1, 12).cuspidal_subspace().gen(0) + M14 = ModularForms(1, 14) + A = SiegelModularForm(60 * E4, M6(0), prec=prec, name='Igusa_4') + B = SiegelModularForm(-84 * E6, M8(0), prec=prec, name='Igusa_6') + C = SiegelModularForm(M10(0), -Delta, prec=prec, name='Igusa_10') + D = SiegelModularForm(Delta, M14(0), prec=prec, name='Igusa_12') + # TODO: the base_ring of A, B, ... should be ZZ + # Here a hack for now: + a = [A, B, C, D] + b = [] + from sage.rings.all import ZZ + for F in a: + c = F.coeffs() + for f in c.iterkeys(): + c[f] = ZZ(c[f]) + F = parent.element_class(parent=parent, weight=F.weight(), coeffs=c, prec=prec, name=F.name()) + b.append(F) + if weights == 'even': + return b + if weights == 'all': + from fastmult import chi35 + coeffs35 = chi35(prec, b[0], b[1], b[2], b[3]) + R = KleinFourGroupAlgebra(QQ) + det = R(R.group().gen(0)) + E = parent.element_class(parent=parent, weight=35, coeffs=coeffs35, prec=prec, name='Delta_35') + E = E * det + b.append(E) + return b + raise ValueError("weights = '%s': should be 'all' or 'even'" % weights) + if group == Sp4Z and weights == 'even' and 2 == coefficient_degree: + b = _siegel_modular_forms_generators(parent=parent) + c = [] + for F in b: + i = b.index(F) + for G in b[(i + 1):]: + c.append(F.satoh_bracket(G)) + return c + raise NotImplementedError("Not yet implemented") + + +@cached_function +def KleinFourGroupAlgebra(base_ring): + """ + Return the group algebra of the Klein four group over the given ring. + + EXAMPLES:: + + sage: from sage.modular.siegel.siegel_modular_forms_algebra import KleinFourGroupAlgebra + sage: KleinFourGroupAlgebra(ZZ) + Group algebra of group "The Klein 4 group of order 4, as a permutation group" over base ring Integer Ring + """ + from sage.groups.all import KleinFourGroup + from sage.algebras.all import GroupAlgebra + return GroupAlgebra(KleinFourGroup(), base_ring) diff --git a/src/sage/modular/siegel/siegel_modular_group.py b/src/sage/modular/siegel/siegel_modular_group.py new file mode 100644 index 00000000000..d9fd9ddab61 --- /dev/null +++ b/src/sage/modular/siegel/siegel_modular_group.py @@ -0,0 +1,118 @@ +""" +Arithmetic subgroups of {GSp}(4,QQ) +""" +import sage.groups.group as group + + +class GSp4_Arithmetic_Subgroups(group.Group): + """ + Arithmetic subgroups (finite index subgroups of {GSp}(4,QQ)). + + EXAMPLES:: + + sage: sage.modular.siegel.siegel_modular_group.GSp4_Arithmetic_Subgroups() + Generic arithmetic subgroup of GSp(4,Q) + """ + def _repr_(self): + """ + EXAMPLES:: + + sage: sage.modular.siegel.siegel_modular_group.GSp4_Arithmetic_Subgroups() + Generic arithmetic subgroup of GSp(4,Q) + """ + return "Generic arithmetic subgroup of GSp(4,Q)" + + def __reduce__(self): + """ + Helper for pickle + """ + raise NotImplementedError("all subclasses must define a __reduce__ method") + + +class Sp4Z_class(GSp4_Arithmetic_Subgroups): + """ + The full Siegel modular group {Sp}(4,ZZ), regarded as a congruence subgroup of itself. + + EXAMPLES:: + + sage: G = Sp4Z; G + Siegel Modular Group Sp(4,Z) + + sage: G.gens() + Traceback (most recent call last): + ... + NotImplementedError: Number of generators not known. + """ + def level(self): + return 1 + + def _repr_(self): + """ + Return the string representation of self. + """ + return "Siegel Modular Group Sp(4,Z)" + + def __reduce__(self): + return _Sp4Z_ref, () + + +def _Sp4Z_ref(): + """ + Return Sp4Z, (Used for pickling Sp4Z) + + EXAMPLES:: + + sage: sage.modular.siegel.siegel_modular_group._Sp4Z_ref() is Sp4Z + True + """ + return Sp4Z + + +Sp4Z = Sp4Z_class() + +_sp4z_gamma0_cache = {} + + +def Sp4Z_Gamma0_constructor(N): + """ + Return the congruence subgroup Gamma0(N) of the Siegel modular group. + + EXAMPLES:: + + sage: G = Sp4Z_Gamma0(51) ; G # indirect doctest + Siegel Modular Group Gamma0(51) + sage: G == Sp4Z_Gamma0(51) + True + """ + if N == 1: + return Sp4Z + try: + return _sp4z_gamma0_cache[N] + except KeyError: + _sp4z_gamma0_cache[N] = Sp4Z_Gamma0_class(N) + return _sp4z_gamma0_cache[N] + + +class Sp4Z_Gamma0_class(GSp4_Arithmetic_Subgroups): + """ + The congruence subgroup Gamma_0(N) of the Siegel modular group + + EXAMPLES:: + + sage: G = Sp4Z_Gamma0(11); G + Siegel Modular Group Gamma0(11) + sage: loads(G.dumps()) == G + True + """ + + def __init__(self, level): + self.__level = level + + def __repr__(self): + return "Siegel Modular Group Gamma0(%d)" % self.level() + + def level(self): + return self.__level + + def __reduce__(self): + return Sp4Z_Gamma0_constructor, (self.level(),) diff --git a/src/sage/modular/siegel/theta_constant.py b/src/sage/modular/siegel/theta_constant.py new file mode 100644 index 00000000000..4cbe79af840 --- /dev/null +++ b/src/sage/modular/siegel/theta_constant.py @@ -0,0 +1,85 @@ +""" +Theta constants +""" + + +def _compute_theta_char_poly(char_dict, f): + r""" + Return the coefficient at ``f`` of the Siegel modular form + + .. math: + + \sum_{l \in\text{char_dict}} \alpha[l] * \prod_i \theta_l[i](8Z). + + INPUT: + + - ``char_dict`` -- a dictionary whose keys are *tuples* of theta + characteristics and whose values are in some ring. + - ``f`` -- a triple `(a,b,c)` of rational numbers such that the + quadratic form `[a,b,c]` is semi positive definite + (i.e. `a,c \geq 0` and `b^2-4ac \leq 0`). + + EXAMPLES:: + + sage: from sage.modular.siegel.theta_constant import _multiply_theta_char + sage: from sage.modular.siegel.theta_constant import _compute_theta_char_poly + sage: theta_constants = {((1, 1, 0, 0), (0, 0, 1, 1), (1, 1, 0, 0), (0, 0, 1, 1)): 1} + sage: _compute_theta_char_poly(theta_constants, [2, 0, 10]) + 32 + sage: _compute_theta_char_poly(theta_constants, [2, 0, 2]) + 8 + sage: _compute_theta_char_poly(theta_constants, [0, 0, 0]) + 0 + """ + return sum(_multiply_theta_char(list(l), f) * char_dict[l] + for l in char_dict) + + +def _multiply_theta_char(l, f): + r""" + Return the coefficient at ``f`` of the theta series `\prod_t \theta_t` + where `t` runs through the list ``l`` of theta characteristics. + + INPUT: + + - ``l`` -- a list of quadruples `t` in `{0,1}^4` + - ``f`` -- a triple `(a,b,c)` of rational numbers such that the + quadratic form `[a,b,c]` is semi positive definite + (i.e. `a,c \geq 0` and `b^2-4ac \leq 0`). + + EXAMPLES:: + + sage: from sage.modular.siegel.theta_constant import _multiply_theta_char + sage: from sage.modular.siegel.theta_constant import _compute_theta_char_poly + sage: tc = [(1, 1, 0, 0), (0, 0, 1, 1), (1, 1, 0, 0), (0, 0, 1, 1)] + sage: _multiply_theta_char(tc, [2, 0, 6]) + -32 + sage: _multiply_theta_char(tc, [2, 0, 10]) + 32 + sage: _multiply_theta_char(tc, [0, 0, 0]) + 0 + """ + if 0 == len(l): + return (1 if (0, 0, 0) == f else 0) + a, b, c = f + m1, m2, m3, m4 = l[0] + # if the characteristic is not even: + if 1 == (m1 * m3 + m2 * m4) % 2: + return 0 + coeff = 0 + from sage.misc.all import isqrt, xsrange + for u in xsrange(m1, isqrt(a) + 1, 2): + for v in xsrange(m2, isqrt(c) + 1, 2): + if 0 == u and 0 == v: + coeff += _multiply_theta_char(l[1:], (a, b, c)) + continue + ap, bp, cp = (a - u * u, b - 2 * u * v, c - v * v) + if bp * bp - 4 * ap * cp <= 0: + val = (2 if 0 == (u * m3 + v * m4) % 4 else -2) + coeff += val * _multiply_theta_char(l[1:], (ap, bp, cp)) + if u != 0 and v != 0: + ap, bp, cp = (a - u * u, b + 2 * u * v, c - v * v) + if bp * bp - 4 * ap * cp <= 0: + val = (2 if 0 == (u * m3 - v * m4) % 4 else -2) + coeff += val * _multiply_theta_char(l[1:], (ap, bp, cp)) + return coeff