lpmn wrong results for |z| > 1 (Trac #1877)

[scipy-gitbot]

Original ticket http://projects.scipy.org/scipy/ticket/1877 on 2013-03-26 by @pv, assigned to @pv.

The lpmn function appears to produce wrong factors of -1 and i for |z| > 1:

>>> special.lpmn(3, 3, 1.5)[0].T[2,1]
5.0311529493745271
>>> mpmath.legenp(2, 1, 1.5)
mpc(real='-1.6342848428326313e-24', imag='-5.0311529493745271')

[scipy-gitbot]

@pv wrote on 2013-03-26

Doesn't affect lpmv which is (arbitrarily?) restricted to |z| < 1

[scipy-gitbot]

@pv wrote on 2013-03-27

More info here: #482

[rgommers]

Not a blocker it looks like to me.

[gertingold]

I do not think that there really is a problem. At least from checking the above example with n=2, m=1 I find agreement with the representation of the associated Legendre function of the second kind P in terms of the hypergeometric function given in http://dlmf.nist.gov/14.3#E8. mpmath actually provides two versions of the Legendre function which are obtained by setting type to 2 (default) or 3. Agreement with http://dlmf.nist.gov/14.3 is obtained for type=3 which differs from type=2 by a factor of (-1)**(m/2). The result given by SciPy is always real in agreement with http://dlmf.nist.gov/14.3. It might be useful to add a remark on the phase factor to the documentation of lpmn. In view of http://dlmf.nist.gov/14.3#E10 it might also be useful to specify in the documentation which function Q is actually meant. How about looking into that during the sprint at EuroSciPy?

There is actually another issue with the Legendre functions. Instead of returning infinity for appropriate arguments, they return 10**300 as generated by the Fortran code. These special cases should probably be handled correctly by the Python code. Otherwise relations like http://dlmf.nist.gov/14.2#E5 and http://dlmf.nist.gov/14.2#E11 will yield incorrect results for arguments very close to 1.

[pv]

The reason why the function definition switches at z=1 indeed seems to be the desire to get real-valued results on the real axis. Fixing just the documentation to state this is probably good enough.

I'm not so sure that this is good behavior for the complex-valued function, however, as it messes up analytic properties by introducing a discontinuity essentially arbitrarily at |z|=1.

[gertingold]

I agree, there is an issue with CLPMN. The results of SciPy are in agreement with http://dlmf.nist.gov/14.21 for |z|>1 but the analytical continuation to |z|<1 is not done correctly. If the analytical continuation is done correctly, approaching the real line -1<x<1 from above or below yields an extra phase factor with respect to the real LPMN in that interval (cf. http://dlmf.nist.gov/14.23#E1). This behavior probably needs to be documented. There is a danger that correcting CLPMN breaks existing code but on the other hand CLPMN is probably not used very often.

[pv]

Fixed in 44e7c61

[pv]

The m < 0 transformation done on the Python side doesn't seem to agree with mpmath, so it may be that there's still something to fix.

[gertingold]

In my opinion basic.py does the correct mapping from negative m to positive m. If it is useful, I can submit a PR with comments in basic.py pointing to the respective relations for Legendre functions.

It should be noted however, that the change from type=2 to type=3 in mpmath not only changes phase factors in the interval (1, inf) but also in the interval (-1, 1):

In [1]: from mpmath import legenp

In [2]: legenp(1, 1, 0.5, type=2)
Out[2]: mpf('-0.8660254037844386')

In [3]: legenp(1, 1, 0.5, type=3)
Out[3]: mpc(real='2.8131368796666611e-25', imag='0.8660254037844386')

In [4]: legenp(1, 1, 1.5, type=2)
Out[4]: mpc(real='-3.6317440951836254e-25', imag='-1.1180339887498949')

In [5]: legenp(1, 1, 1.5, type=3)
Out[5]: mpf('1.1180339887498949')

Could it be that this unexpected behavior in the interval (-1, 1) leads to failing of tests?

A side remark: According to http://dlmf.nist.gov/14.21, the cut is not restricted to [-1, 1] but actually is given by (-inf, 1]. If I should submit a PR as proposed above, I could adapt the documentation of clpmn accordingly in the same PR.

[gertingold]

Just in case it might be useful to somebody: An apparent problem with negative m can be caused by accessing the Legendre function with lpmn(m, n, x)[0][m, n]. The correct form is lpmn(m, n, x)[0][abs(m), n] (I just fell into this trap...).

[pv]

The clpmn function indeed seems to work correctly for negative m. The issue seems to be the (-1)**mf factor in the corresponding transform in lpmn --- this is probably OK for type=2 but not for type=3 function, so it boils down to the definition of the function at |x| > 1 (matching to clpmn seems reasonable).

PR fixing lpmn: gh-2808

About the cut: numerically I don't see it at Re z < -1, but this may be because of the integer order.

[gertingold]

According to DLMF, the connection between positive and negative m is different for Ferrer's function on the cut and for the assosciated Legendre function in the complex plane. According to http://dlmf.nist.gov/14.9#E3, a factor (-1)**mf appears for Ferrer's function while according to http://dlmf.nist.gov/14.9#E13 there is no sign appearing for the associated Legendre function in the complex plane. The latter relation holds also in the complex plane, point (iii) in http://dlmf.nist.gov/14.21.

As the discussion goes on, I am coming more and more to the conclusion that the idea of type=2 and type=3 in mpmath is not a good one. For me, it makes much more sense to define an analytic function in the complex plane with a well defined convention (I would prefer following DLMF here, i.e. type=3 in mpmath). This implies a cut at least on the interval (-1, 1). In addition, one can have a second function which lives (only) on this cut. This function should be related to the Legendre polynomials in the way expected by the user as is the case for lpmn and type=2 in mpmath. The fact that lpmn extends the range of acceptable arguments beyond (-1, 1) is a sort of bonus which on the other hand gives rise to questions about phases. Therefore, it might actually make sense to restrict lpmn to the interval (-1, 1), thus avoiding phase ambiguities.Such a restriction would be consistent with DLMF. On the other hand, there exist situations where due to a Wick rotation one needs to evaluate associated Legendre polynomials for imaginary angles. The Wick rotation effectively maps arguments smaller than 1 into arguments larger than 1. In such a situation, a function lpmn working on the real axis outside the cut is useful.

I am not completely sure why the cut in DLMF includes the complete negative real axis. May be it is due to the fact that sqrt(z**2-1)=sqrt(z+1)*sqrt(z-1) gives rise to cuts starting out from 1 towards negative real values and from -1 towards negative real value. Then the cut is given by the interval (-inf, 1]. On the other hand, I agree that the associated Legendre function should not jump on the real axis below -1.

[pv]

This was addressed in the PRs.