Lerch transcendental function

[albinahlback]

I know how to efficiently calculate Lerch transcendental function for integer s \le 1 (maybe for s = 2 as well) when |x| < 1 and as well as when x 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.

[fredrik-johansson]

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.

[fredrik-johansson]

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).

[albinahlback]

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?

[fredrik-johansson]

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.

[fredrik-johansson]

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.

[albinahlback]

Is that with correct error bounds? If so, can you share the proofs for those error bounds used?

[fredrik-johansson]

Yes. I'm going to write things down somewhere.

[fredrik-johansson]

I put it on my blog: https://fredrikj.net/blog/2022/02/computing-the-lerch-transcendent/

[albinahlback]

Truly great work!!