More lerchphi strangeness

[GoogleCodeExporter]

I am not quite a sure about this bug report as about the previous one. 
According to my computations, when passing across the branch cut from 1 to oo 
from above, on the standard branch, the lerch transcendent should pick up a 
term 2*pi*I*log(z)**(s-1)/gamma(s)/z**a [where the standard choice of argument 
is made for z and also for log(z)], at least for re(a) > 0. Indeed:

In [1]: def a(z, s, a, eps=1e-15): return mpmath.lerchphi(z + eps*1j, s, a) - 
mpmath.lerchphi(z - eps*1j, s, a)
In [2]: def b(z, s, a): return (2*pi*I*log(z)**(s-1)/z**a/gamma(s)).n()

In [3]: a(1.7, 5.5, 3.5) - b(1.7, 5.5, 3.5)
Out[3]: 2.16840434497101e-19⋅ⅈ
In [8]: a(2.5, 1.5, 3.9) - b(2.5, 1.5, 3.9)
Out[8]: 2.77555756156289e-17⋅ⅈ
In [9]: a(2.5, 1.5+1j, 3.9) - b(2.5, 1.5+1j, 3.9)
Out[9]: -1.38777878078145e-17 + 5.55111512312578e-17⋅ⅈ

However:

In [10]: a(2.5, 1.5, 3.9+1j) - b(2.5, 1.5+1j, 3.9+1j)
Out[10]: -0.55002171991196 + 0.166746405436317⋅ⅈ

Indeed something weird is going on:

In [42]: a(2.5, 1.5, 3.9 - 1e-10j, 1e-15)
Out[42]: (0.4646322489177 + 0.19040815320768j)
In [43]: a(2.5, 1.5, 3.9 + 1e-10j, 1e-15)
Out[43]: (-0.4646322489177 + 0.19040815320768j)
In [44]: a(2.5, 1.5, 3.9 + 0*1e-10j, 1e-15)
Out[44]: (0.0 + 0.19040815319107j)

[i.e. lerchphi does not seem to be continuous in a]

Increasing mpmath.mp.prec to 500 does not make this go away.

Now I must admit that I am not entirely sure how the lerch transcendent should 
or could be continued analytically, and even less sure about what is going on 
when computing it numerically. However, I think that

1) lerchphi(z, s, a) should be continuous in a, certainly for fixed s and z and 
re(a) > 0.
2) the branching behaviour I mentioned first should hold

In summary, it seems that for non-real a, lerchphi is evaluated on a 
non-principal branch.

Original issue reported on code.google.com by ness...@googlemail.com on 22 Jul 2011 at 11:50

[GoogleCodeExporter]

The implementation of lerchphi is very simple. It uses a recurrence to force 
Re(a) > 0, and then uses an integral representation. This is all described in 
the documentation.

I'm not sure if I see the reason why it should be continuous in a (probably if 
z and/or s are integers, but generally?).

If I remember correctly, I tried to make the function agree with 
HurwitzLerchPhi in Mathematica, but I'm not sure if that's always the case.

There is an older function LerchPhi in Mathematica with a rather dumb branch 
structure, and they introduced HurwitzLerchPhi in 7.0 specifically to address 
that. So I think the Mathematica developers thought carefully about the branch 
cuts, and it would probably be reasonable to copy what they did.

But I'm not sure if the current function accomplishes that. It's unfortunately 
not documented precisely how Mathematica defines the analytic continuation, so 
the only way to find out is by trial and error.

Original comment by fredrik....@gmail.com on 22 Jul 2011 at 5:02

[GoogleCodeExporter]

I'm not claiming that it should be continuous in a everywhere, but certainly 
for Re(a) > 0, Re(s) > 0, and z not in [1, oo), I think it should. Indeed the 
well-known integral representation

http://upload.wikimedia.org/math/6/a/e/6ae504b23019219cc2bad00fe0a9d47a.png
(phi(z, s, a) = 1/gamma(s) int_0^infty t^{s-1}e^{-at}/(1 - ze^{-t}} dt)

which provides an analytic continuation of the series definition into said 
domain seems continuous to me. Also the integral representation you are 
providing in the docs seems to have this property, although I'm not entirely 
sure about the branches of arctan there.

Again I'm not claiming to be an expert at all, I'm just working to add lerchphi 
etc to sympy, and trying to understand its behaviour on the way :-).

Original comment by ness...@googlemail.com on 22 Jul 2011 at 9:40

[GoogleCodeExporter]

Yes, it might be the arctangent in the integrand that's to blame. The integral 
representation used is the one taken from the DLMF, although it's stated only 
to be valid for Re(s) > 0, |z| < 1 there.

The Wikipedia article contains another version of the "Hermite-like integral 
representation" with the branch dependence on the parameter essentially 
factored out. I'm not able to work on this right now; you could try it out if 
you like.

Original comment by fredrik....@gmail.com on 23 Jul 2011 at 1:02