Open albinahlback opened 3 years ago
A few relevant papers:
https://carma.newcastle.edu.au/resources/jon/lerch.pdf
https://gnpalencia.org/research/GNP_Lerch_2018.pdf
https://arxiv.org/pdf/1005.4967.pdf
A general method to compute phi(z,s,a) is to use recurrence relations to ensure Re(a) >= 1 and then evaluate the Hermite integral representation by direct numerical integration.
I would maybe start with this and then add various series expansions and Euler-Maclaurin summation where applicable to speed up common cases.
I believe Algorithm 3 in https://carma.newcastle.edu.au/resources/jon/lerch.pdf should be usable as a general solution.
from mpmath import *
mp.dps = 30
mp.pretty = True
def S(z):
return sign(re(z))
def phi(z,s,a):
z = mpmathify(z)
s = mpmathify(s)
a = mpmathify(a)
res = 0
l = +pi
if isint(a):
res -= z**(-a) * l**(s/2) / s * rgamma(s/2)
if z == 1:
res += 1/(s-1) * pi**0.5 * l**((s-1)/2) * rgamma(s/2)
def f(n):
A = n + a
return z**n / (A**2)**(s/2) * (gammainc(s/2, l*A**2)*rgamma(s/2) + gammainc((s+1)/2, l*A**2)*rgamma((s+1)/2) * S(A))
res += 0.5 * nsum(f, [-inf,inf], method="d")
def g(u):
U = u + log(z)/(2*pi*1j)
return exp(-2*pi*1j*a*u)/((U**2)**((1-s)/2)) * (gammainc((1-s)/2, pi**2*U**2/l)*rgamma(s/2) + 1j*gammainc(1-s/2, pi**2*U**2/l)*rgamma((s+1)/2)*S(U))
res += pi**(s-0.5) / (2*z**a) * nsum(g, [-inf,inf], method="d")
return res
z, s, a = 2+3j, 1-2j, 0.75-0.5j
print(phi(z,s,a))
print(lerchphi(z,s,a))
This needs a bit of work to turn into a complete algorithm, though, and it needs some branch cut fix for z in [0, inf).
Seems good, but the error term of it seems daunting to compute (perhaps you have some deeper knowledge in the error term of the gamma function). And I can't seem to find the rate of convergence, do you happen to know how fast it should converge?
Gamma(a,z) converges like exp(-z) when z goes to +infinity, so the convergence is quite strong; we only need O(sqrt(prec)) terms.
For real positive z with z > Re(a), it is possible to use a bound like |Gamma(a,z)| <= max(1, 2^Re(a)) z^(Re(a)-1) exp(-z). We would also need complex z here, and then things get a bit more complicated. The formulas in https://dlmf.nist.gov/13.7#ii are a possible starting point.
I've recently been working on a restricted version of this for Dirichlet L-functions; I should have some time to look at Lerch when I'm done with that.
There is now a first implementation in 1abaf57c00ef48570a756667920c0ad683edb857 using direct summation for |z| << 1 and the Hankel contour integral for the analytic continuation.
Numerical integration isn't the fastest but it should at least work as a fallback and as a comparison baseline for more efficient methods.
Is that with correct error bounds? If so, can you share the proofs for those error bounds used?
Yes. I'm going to write things down somewhere.
I put it on my blog: https://fredrikj.net/blog/2022/02/computing-the-lerch-transcendent/
Truly great work!!
I know how to efficiently calculate Lerch transcendental function for integer
s \le 1
(maybe fors = 2
as well) when|x| < 1
and as well as whenx
is a root of unity. If someone knows how to efficiently calculate the rest of the cases, it might be worthwhile to implement the function.Edit: Seems to be easy to calculate for positive integers
s
.