#===============================================================================
#
# 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 Equation (15) of 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.
#
###############################################################################
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)]
var('v k')
hyperord = 15
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 differential operator attached to the equation for h1
op = lambda f: - 4 * k**3 * diff(f, v) - 20 * k**2 * v * diff(f, v, v) \
- 32 * k * v**2 * diff(f, v, v, v) - 16 * v**3 * diff(f, v, v, v, v) \
- 10 * k**2 * diff(f, v) + 4 * k * v * diff(f, v) \
- 60 * k * v * diff(f, v, v) + 4 * v**2 * diff(f, v, v) \
- 64 * v**2 * diff(f, v, v, v) + 2 * k * f \
- 2 * k * diff(f, v) + 4 * v * diff(f, v) \
- 28 * v * diff(f, v, v) - f \
+ 4 * diff(f, v)
h1ims = map(lambda e: e.coefficients(v), map(op, h1exps))
## We first check that Sage has not automatically eliminated powers of v
## The first terms that we have to check are
## [0, -k/2 - 4, (-7 - k)/2, -5/2 - k]
assert [h1im[0][1] for h1im in h1ims] == [0, -1/2*k - 1, -1/2*k - 1/2, -k + 1/2]
assert all( [e.simplify_rational() for (e,_) in h1im[:11]] == 11*[0] for h1im in h1ims )