#===============================================================================
#
# Copyright (C) 2011 Martin Raum
#
# This program is free software; you can redistribute it and/or
# modify it under the terms of the GNU General Public License
# as published by the Free Software Foundation; either version 3
# of the License, or (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
# General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program; if not, see .
#
#===============================================================================
#===============================================================================
#
# This file concerns the case rank T =2 in the proof of Theorem 3 of the paper
# "Kohnen's limit process for real-analytic Siegel modular forms"
# by Kathrin Bringmann, Martin Raum, and Olav K. Richter.
# For detailed definitions see the paper.
#
###############################################################################
## The definition of the Pochhammer symbol and of the power series expansion
## of hypergeometric series around zero.
pochhammer = lambda a, n: prod(a + k for k in range(n))
hyperexpansion = lambda ass, bss, ord: [prod(pochhammer(a, n) for a in ass) / prod(pochhammer(b, n) for b in bss) / factorial(n) for n in range(ord)]
## the weight k, and the two variables u and v, that Maass uses.
var('k v u')
## We will expand all series up to powers of order 3, since there is cancellation.
## We will only use the first two nonvanishing coeffcients.
hyperord = 3
###############################################################################
##
## the indefinite case: We prove that only one linear combination
## of the four solutions can be mapped to zero. This proves the
## indefinite case, since the coefficients for indefinite indices of holomorphic
## Siegel modular forms vanish.
##
###############################################################################
## The parameters for the hypergeometric series that occur, and
## the corresponding power series expansions
h1s = [(0 , [1/2], [(1 + k)/2, 1 + k/2]),
(-k/2 , [(1 - k)/2], [1/2, 1 - k/2]),
((1 - k)/2, [1 - k/2], [3/2, (3 - k)/2]),
(3/2 - k , [1, 2 - k], [5/2 - k, 2 - k/2, (5 - k)/2]),
]
h1exps = \
[v^e * sum(map( operator.mul, hyperexpansion(ass, bss, hyperord),
[(v/4)^n for n in range(hyperord)] ))
for (e, ass, bss) in h1s]
## The two Meijer G-functions that occur if k > 3. The expressions in terms
## of the hypergeometric series is due to Luke
h1exps.append( gamma(k/2 - 1) * gamma((5 - k)/2) * gamma(5/2 - k)**-1 * v**(3/2 - k) * h1exps[3] \
+ gamma(1 - k/2) * gamma(3/2)**-1 * gamma((3 - k)/2)**-1 * h1exps[2] )
h1exps.append( gamma((k - 3)/2) * gamma(2 - k/2)**-1 * gamma(5/2 - k)**-1 * gamma(1 - k)**-1 * h1exps[3] \
+ gamma((3 - k)/2) * gamma((k - 1)/2) * gamma(1/2)**-1 * gamma(1 - k/2)**-1 * gamma((1 + k)/2)**-1 * h1exps[1] )
## This are the asymptotic expansions for v -> 0 of a(Y, I_2^sk) restricted to u = 0
bs = [ e
+ (2 * v - 2 * (k + 1) * (k + 2)) * diff(e, v)
+ 8 * pi^2 * v * (2 * (k + 1) + (4 * pi^2 + 1) * (k + 4)) * diff(e, v, v)
- 64 * pi^4 * v * (4 * pi^2 + 1) * diff(e, v, v, v)
for e in h1exps]
## We need the initial exponents of the asymptotic expansions:
assert bs[0].coefficients(v)[0][1] == 1
assert bs[1].coefficients(v)[0][1] == -1/2*k - 2
assert bs[2].coefficients(v)[0][1] == -1/2*k - 3/2
assert bs[3].coefficients(v)[0][1] == -k - 1/2
assert bs[4].coefficients(v)[0][1] == -2*k + 1
assert bs[5].coefficients(v)[0][1] == -k - 1/2
## Next, we compute under which conditions the initial coefficient vanishes.
## Except for one value of k, these are all values that we are not interested in.
assert bs[0].coefficients(v)[0][0].solve(k) == [k == -1/2*(32*pi^4 + 12*pi^2 + 5)/(4*pi^4 + 3*pi^2 + 1)]
assert bs[1].coefficients(v)[0][0].solve(k) == [k == -4, k == -2, k == 0]
assert bs[2].coefficients(v)[0][0].solve(k) == [k == -3, k == -1, k == 1]
assert bs[3].coefficients(v)[0][0].solve(k) == [k == (-1/2), k == (1/2), k == (3/2)]
assert bs[4].coefficients(v)[0][0].solve(k) == [k == (1/2), k == (3/2), k == 1, gamma(1/2*k - 1) == 0, gamma(-1/2*k + 5/2) == 0]
assert bs[5].coefficients(v)[0][0].solve(k) == [k == (-1/2), k == (1/2), k == (3/2), gamma(1/2*k - 3/2) == 0]
## if k > 3 then the first, fifth and sixth series are relevant. We see that the initial
## coefficients never vanish and that the initial exponents are different. This proves
## the statement as explained in the paper.
## if k < 0 the second to fourth series are relevant. The second and third have
## vanishing initial coefficients in some cases. We check two more coefficients
## to prove that in all cases all initial exponents are pairwise distinct.
assert bs[1].coefficients(v)[1][1] == -1/2*k - 1
assert bs[1].coefficients(v)[1][0].solve(k) == [k == (4*pi^4 - 8*pi^2 - 1)/(4*pi^4 + 2*pi^2 + 1), k == -2, k == 0]
assert bs[1].coefficients(v)[2][1] == -1/2*k
assert bs[1].coefficients(v)[2][0].solve(k) == [k == 3*(8*pi^6 - 30*pi^4 - 12*pi^2 - 3)/(8*pi^6 + 26*pi^4 + 18*pi^2 + 3), k == 0, k == 1]
assert bs[2].coefficients(v)[1][1] == -1/2*k - 1/2
assert bs[2].coefficients(v)[1][0].solve(k) == [k == 2*(16*pi^6 - 44*pi^4 - 18*pi^2 - 3)/(16*pi^6 + 28*pi^4 + 18*pi^2 + 3), k == -1, k == 1]
assert bs[2].coefficients(v)[2][1] == -1/2*k + 1/2
assert bs[2].coefficients(v)[2][0].solve(k) == [k == 1, k == 4*(8*pi^6 - 38*pi^4 - 15*pi^2 - 5)/(8*pi^6 + 42*pi^4 + 30*pi^2 + 5), k == 2]
###############################################################################
##
## the definite case: We prove that at most one a(Y, I_2) is mapped to
## the exponential exp(-tr Y). If k > 3 we also show that no nonzero
## solution is mapped to zero.
##
###############################################################################
## The integral that gives psi if phi is inserted
integ = lambda e: sum(map(lambda v: v[0]/(v[1] + 1) * u^(v[1] + 1), e.coefficients(u)))
## The two fundamental solutions for the differential equations.
## We assume that k is generic.
phi0 = exp(-u).taylor(u,0,hyperord) * u^(k - 1) \
* sum(map( operator.mul, hyperexpansion([2 * k - 2], [2 * k - 2], hyperord),
[(2 * u)^n for n in range(hyperord)]))
phi1 = exp(-u).taylor(u,0,hyperord) * u^(2 - k) \
* sum(map( operator.mul, hyperexpansion([1], [4 - 2 * k], hyperord),
[(2 * u)^n for n in range(hyperord)]))
## We have three fundamental solutions for g0
g0s = [ u^(1 - k) * integ(u^(-1) * phi0),
u^(1 - k) * integ(u^(-1) * phi1),
u^(1 - k) ]
## And these are the resulting asymptotic expansions of a(Y, I_2) restricted to
## v = 0 when u -> 0
bs = [ u^(2 * k - 1) *
( (- 3/4 - (1 - k) / 2 * u^(-1)) * g0
- 5/4 * diff(g0, u, u)
+ (2 - (2 + k) / 2 * u^(-1)) * diff(g0, u) )
for g0 in g0s ]
## The first coefficient of the first solution vanishes, but the next coefficient never does
assert bs[0].coefficients(u)[0][0] == 0
assert bs[0].coefficients(u)[1][0].simplify_rational() == -1/k
## its index is 2k - 2
assert bs[0].coefficients(u)[1][1] == 2 * k - 2
## The second solution begins with u^0 except in two cases, that do not occur
## for Siegel modular forms:
assert bs[1].coefficients(u)[0][1] == 0
assert bs[1].coefficients(u)[0][0].simplify_rational() == 1/2*(8*k^2 - 26*k + 21)/(k - 2)
assert bs[1].coefficients(u)[0][0].simplify_rational().solve(k) == [k == (3/2), k == (7/4)]
## The third solution begins with u^(k - 2). Again, only in cases that do not
## occur in the range k < 0 or k > 3 the first coefficient vanishes.
assert bs[2].coefficients(u)[0][1] == k - 2
assert bs[2].coefficients(u)[0][0].simplify_rational().solve(k) == [k == 1, k == (4/3)]
## Since the initial coefficients are different, only one linear combination
## will be mapped to the exponential exp(-tr Y), and only the trivial combination will
## be mapped to the zero function.