Special case for gen_legendre_P with n == m gives incorrect results for x in (-1,1)

[jcwomack]

I am using the gen_legendre_P function (an instance of Func_assoc_legendre_P from sage/functions/orthogonal_polys.py) to evaluate associated Legendre functions / Ferrers functions in SageMath 8.1.

There appears to be a discrepancy in the results I obtain, depending on whether I use gen_legendre_P.eval_poly() or directly call gen_legendre_P() in some cases. I think this is because the _eval_ method first tries to call the _eval_special_values_ method, before using eval_poly.

With this input

x = SR.var('x')
print(gen_legendre_P.eval_poly(1, 1, x))
print(gen_legendre_P(1, 1, x))
print(gen_legendre_P.eval_poly(1, 1, 0.5))
print(gen_legendre_P(1, 1, 0.5))

I obtain

-sqrt(-x^2 + 1)
sqrt(x^2 - 1)
-0.866025403784439
5.30287619362453e-17 + 0.866025403784439*I

The result from eval_poly agrees with Mathematica, i.e.

LegendreP[1, 1, 0.5]
-0.866025

Based on the above output, it seems to me that gen_legendre_P.eval_poly(1, 1, cos(theta)) will always be real while gen_legendre_P(1, 1, cos(theta)) will be complex (unless |cos(theta)| = 1), since cos(theta) is in the interval [-1,1].

Looking at the code for Func_assoc_legendre_P._eval_special_values_, I suspect the culprit is the n == m case, which returns

factorial(2*m)/2**m/factorial(m) * (x**2-1)**(m/2)

This discrepancy also seems to be present in spherical_harmonic when n == m (an instance of SphericalHarmonic from sage/functions/special.py), which is built using gen_legendre_P.

After ​discussion in the sage-devel mailing list it appears that this is because the n == m case in _eval_special_values_ is based on ​https://dlmf.nist.gov/14.7#E15, but this is not defined in (-infinity,1].

For the spherical harmonics, where the argument x = cos(theta), x will always be in the range [-1, 1 ], where special case used in _eval_special_values_ is not defined.

On the sage-devel mailing list, Howard Cohl suggested that the correct formula for x in [-1, 1] is

P_m^m(x)=(-1)^m (2m)!/(2^m m!) (1-x^2)^(m/2)

See: ​https://groups.google.com/d/msg/sage-devel/IDtiGF6HB28/QWwnAeLJBAAJ

According to Howard Cohl this is formally a Ferrers function (defined on (-1,1) ), rather than an associated Legendre polynomial. However, the existing code for Func_assoc_legendre_P does not seem to make any distinction between Ferrers and associated Legendre functions.

My proposed fix would be to have Func_assoc_legendre_P._eval_special_values_ choose between two n == m special cases, based on whether -1 <= x <= 1 (above expression) or > 1 (current expression).

This raises the question of whether Func_assoc_legendre_P is correctly defined, as at present it would seem to cover both Ferrers functions and associated Legendre functions.

In my experience with the physics/chemistry literature, the spherical harmonics are universally defined in terms of "associated Legendre functions", even though the argument is x = cos(theta). DLMF suggests these are defined in terms of Ferrers functions of the first kind (​https://dlmf.nist.gov/14.30.E1). Wolfram Mathematica does not seem to distinguish. Possibly it is worth flagging in the docstring for Func_assoc_legendre_P that the class seems to cover both functions.

[egourgoulhon]

A (quite severe IMHO) consequence of this bug is

sage: theta, phi = var('theta phi')
sage: spherical_harmonic(1, 1, theta, phi)
-1/4*sqrt(3)*sqrt(2)*sqrt(cos(theta)^2 - 1)*e^(I*phi)/sqrt(pi)

which is plain wrong: the term sqrt(cos(theta)^2-1) (which is imaginary for theta real) should actually be sin(theta).

[egourgoulhon]

For reference, here are some related issues reported by users:

• ​https://groups.google.com/forum/#!msg/sage-support/I_d_meMxRbM/bDhZ_wtoAgAJ
• ​https://ask.sagemath.org/question/49694
• ​https://ask.sagemath.org/question/48776

[egourgoulhon]

Replying to egourgoulhon:

For reference, here are some related issues reported by users:

• ​https://groups.google.com/forum/#!msg/sage-support/I_d_meMxRbM/bDhZ_wtoAgAJ
• ​https://ask.sagemath.org/question/49694
• ​https://ask.sagemath.org/question/48776

One more:

• ​https://ask.sagemath.org/question/49891/

[egourgoulhon]

For reference, here are some related issues reported by users:

• ​https://groups.google.com/forum/#!msg/sage-support/I_d_meMxRbM/bDhZ_wtoAgAJ
• ​https://ask.sagemath.org/question/49694
• ​https://ask.sagemath.org/question/48776

One more:

• ​https://ask.sagemath.org/question/49891/

Yet one more:

• ​https://ask.sagemath.org/question/51994/