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

[jcwomack]

The current implementation of gen_legendre_P is not accurate. No distinction between Ferrers functions and Legendre functions has been made and some special cases are not correctly implemented. This heavily affects spherical harmonics. For a quick fix, to make spherical harmonics work, I propose to temporarily restrict gen_legendre_P to Ferrers functions. In a follow-up #31637, we will make the distinction complete.

---

Old Description

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/

[gh-mjungmath]

Since this features is obviously highly demanded, and within two years not resolved, I took the freedom to shift it's priority to "critical".

[gh-mjungmath]

What about separately adding Func_ferrers to the library with the appropriate convention then?

One can easily refer to ​https://dlmf.nist.gov/14.7#E8 and ​https://dlmf.nist.gov/14.7#E14 respectively and clarify the convention.

Also, I found a nice stackexchange post explaining the origin of this difference: ​https://math.stackexchange.com/a/2986444

[gh-mjungmath]

I already worked on something. Will push tomorrow.

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​568d0fa	Trac #25034: follow mathematica convention on legendre polynomials
​6d6ead1	Trac #25034: doctest for spherical harmonics added
​571632f	Trac #25034: doctest further improved

[git]

Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:

​4fcc574

Trac #25034: formula for n=m fixed

[gh-mjungmath]

I have uploaded a proposal. Please let me know in case I have missed any details.

What have I done?

• The main modification is as proposed:

           if m == 0:
               return legendre_P(n, x)
           if n == m:
  -            return (-1)**m*(1-x**2)**(m/2) * ZZ(2*m-1).multifactorial(2)
  +            return (-1)**m*factorial(2*m)/(2**m*factorial(m)) * (1-x**2)**(m/2)
  

• Furthermore I have added the feature of negative m which wasn't there yet.
• The special case m+1=n has been added.

You are invited to contribute more doctests, especially for spherical harmonics. So far, I have only collected some of the code snippets posted above that have not worked. They are working now.

Follow-Ups

• As for resolving convention conflicts, I propose to wrap this up in a follow-up ticket (#31637).

• Moreover, imho the documentation for orthogonal polynomials is a mess and urgently needs to be cleaned up. I suggest to attack this in another follow-up ticket (#31636).

[egourgoulhon]

Thank you so much for tackling this!

Shall the ticket be set to needs_review? (in order to draw it to the patchbots' attention)

[gh-mjungmath]

Replying to egourgoulhon:

Thank you so much for tackling this!

Shall the ticket be set to needs_review? (in order to get the attention of the patchbots)

Would you mind to add some more doctests for spherical coordinates before? I think, the current doctesting is a little bit sparse considering that the bug has not been caught. At least I would feel more comfortable if some identities or important values are checked that haven't worked before.

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​a4edb24

Trac #25034: abs(m) exceeding n gives zero immediately

[gh-mjungmath]

Replying to gh-mjungmath:

Would you mind to add some more doctests for spherical coordinates before? I think, the current doctesting is a little bit sparse considering that the bug has not been caught. At least I would feel more comfortable if some identities or important values are checked that haven't worked before.

Alright, technically speaking spherical harmonics do not belong to this ticket. Let's finish this one first before messing with spherical harmonics. I created a follow-up: #31639.

[gh-mjungmath]

Patchbot is green at least.

[tscrim]

The order in which you perform tests is important. Here:

if n in ZZ and m in ZZ and (x in ZZ or not SR(x).is_numeric()) and n >= 0:

you should have the n >= 0 before the x tests as it is quicker.

I am also not sure about doing the multiplication first in the quotient here:

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

I feel like

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

should give better performance since we know the first division is exact. (I really would prefer //, but that won't work if m is a symbolic variable unfortunately...)

Also, a while-we-are-at-it:

-return ZZ(0)
+return ZZ.zero()

I probably would also be good to add support for negative n by replacing it with -n-1. This will get rid of the error for gen_legendre_P(-2,-2,x).

I don't understand why you feed this back through the full evaluation instead of this:

-return x*(2*m+1)*gen_legendre_P(m, m, x)
+return x * (2*m+1) * self._eval_special_values_(m, m, x)

A small tweak in the language:

-Print the first associated Legendre polynomials::
+We give the first associated Legendre polynomials::

I would also add a space around the = to make it easier to read the output.

Please put the REFERENCES: at the end of the docstring (and add the S).

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​7c532a1	Trac #25034: minor improvements
​6b365d7	Trac #25034: slight refactoring of _eval_

[gh-mjungmath]

Thank you Travis for the feedback! I hope I took everything in account now. In addition, I refactored _eval_ a bit and completely removed _eval_special_values_. Please check this out.

[gh-mjungmath]

I know that this is not absolutely precise at the moment. But at least this implementation is more or less consistent and makes spherical harmonics work properly. More elaborate thoughts to resolve this are very much appreciated in #31637.

[gh-mjungmath]

Patchbot green again.

[tscrim]

I am -1 on removing _eval_special_values_ as it has a clear purpose and is useful for future maintenance by having a more logical grouping of the code.

[gh-mjungmath]

Now it is grouped into "integer" and "non-integer" case. Unfortunately, this didn't work well with _eval_special_values_ because it would lead to many repeated if-statements. This is the most concise pattern I could think of. Do you have another idea?

Moreover, the x == 0 case seemingly didn't make sense for _evalf_. Please correct me if I am wrong here.

[tscrim]

I thought the previous version was fine. I don't like forcing it all together, as opposed to what many other functions do. I also don't really think the SR(x).is_numeric() test being run so frequently is good.

Additionally, with your current code, if passing symbolic m == n, it no longer gets the special value AFAICS.

[gh-mjungmath]

Replying to tscrim:

Additionally, with your current code, if passing symbolic m == n, it no longer gets the special value AFAICS.

True. However, this particular formula holds indeed if m and n are integers. But I am not sure whether this extends accordingly to arbitrary m and n using Gamma functions...

Playing around with the numeric tool, it turned out that the case abs(m) > n only yields zero for the integer case (try for example gen_legendre_P(2.,3.*I,.5)).

Grouping it the same way as before would therefore lead to many redundant if-statements, as said above. I don't think this is desirable either.

[git]

Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:

​801298c

Trac #25034: allow negative n

[gh-mjungmath]

What about this?

[gh-mjungmath]

At least the patchbot likes it. This version is very close to the version before, but now with negative n. This is still not optimal but I think more preciseness should be elaborated in the follow-up ticket.

What do you think?

[tscrim]

I think this is a lot better.

This is a bit cumbersome way of saying "let m be a positive integer":

+    For `n` being a non-negative integer, negative integer values for `m`
+    with `|m| \leq n` can be obtained by:
+
+    .. MATH::
+
+            P^{-|m|}_n(x) = (-1)^{|m|} \frac{(n-|m|)!}{(n+|m|)!} P_n^{|m|}(x).

I am not sure about removing the _eval_special_values_ from _evalf_. Can we add doctests for floating point input with special values?

I think we should also add doctests to _eval_special_values_ for symbolic m input as well.

[gh-mjungmath]

Replying to tscrim:

I am not sure about removing the _eval_special_values_ from _evalf_. Can we add doctests for floating point input with special values?

I think it is reasonable to assume that mpmath already uses optimized algorithms, isn't it? Moreover, without the above patch we have something like that:

sage: %%time 
....: gen_legendre_P(1.2,1.1,0)
CPU times: user 2.78 ms, sys: 0 ns, total: 2.78 ms
Wall time: 2.8 ms
2.14354692507259*cos(1.15000000000000*pi)*gamma(33/20)/(sqrt(pi)*gamma(749113601384401/713441525128001))

which already looks clunky. Moreover if we want the answer in terms of floating numbers, we get:

sage: %%time 
....: gen_legendre_P(1.2,1.1,0).n()
CPU times: user 166 ms, sys: 7.01 ms, total: 173 ms
Wall time: 204 ms
-0.996322549249593

whereas compared to with patch:

sage: %%time 
....: gen_legendre_P(1.2,1.1,0)
CPU times: user 302 µs, sys: 20 µs, total: 322 µs
Wall time: 326 µs
-0.996322549249593

For the case m=n however, the mpmath version seems to be slightly slower than evaluating the special case (now everything with the new patch):

sage: %%time 
....: gen_legendre_P._evalf_(1,1,.2)
CPU times: user 1.29 ms, sys: 0 ns, total: 1.29 ms
Wall time: 1.29 ms
-0.979795897113271
sage: %%time 
....: gen_legendre_P._eval_special_values_(1,1,.2)
CPU times: user 124 µs, sys: 0 ns, total: 124 µs
Wall time: 127 µs
-0.979795897113271

But treading these cases differently would dismantle _eval_special_values_ again.

I think we should also add doctests to _eval_special_values_ for symbolic m input as well.

I'll do.

[tscrim]

We need to also be more careful if the type of the output is changing as this can break other code in the wild. It might be good to see how other special functions behave for different inputs to make sure things match if we are going to change stuff.

[git]

Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:

​0b7d917

Trac #25034: Turn gen_legendre_P into Ferrers functions only

[gh-mjungmath]

I had a little chat with Howard Cohl. He said that one has to be extremely careful which formula goes under what condition. In particular, one shouldn't use formulas for cases that are not stated in ​https://dlmf.nist.gov/14.

So I tried to be very careful and always referred to the corresponding formula in the code. But please double check that for me.

Other than that, this is ready for review again.

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​d203066

Trac #25034: wrong indent in docstring

[egourgoulhon]

Thanks for this work! It's nice that you made references to the DLMF for formulas in the code.

In the doctests Check whether :trac:25034 yields correct results compared to Maxima in special.py, shouldn't there be markers # abs tol 1e-14 ? Same remark for

        sage: gen_legendre_P(-5/3,3,1.+I)
        0.238163715352606 + 0.0548443903534220*I

in the docstring of Func_assoc_legendre_P.

[gh-mjungmath]

Replying to egourgoulhon:

In the doctests Check whether :trac:25034 yields correct results compared to Maxima in special.py, shouldn't there be markers # abs tol 1e-14 ? Same remark for

        sage: gen_legendre_P(-5/3,3,1.+I)
        0.238163715352606 + 0.0548443903534220*I

in the docstring of Func_assoc_legendre_P.

Done.

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​c746088

Trac #25034: add abs tol

[gh-mjungmath]

I hope this will make it in Sage 9.3.

[egourgoulhon]

Thanks again! It's great that Sage has decent spherical harmonics at last!

Travis, do you agree with the positive review?

[tscrim]

You have renamed the public function eval_poly, so this need a deprecated_function_alias.

As per comment:39, I consider this to be a regression because you removed _eval_special_values_ from _evalf_:

sage: gen_legendre_P(x,x,.2)                                                                              
gen_legendre_P(x, x, 0.200000000000000)

versus before:

(-0.960000000000000)^(1/2*x)*factorial(2*x)/(2^x*factorial(x))

[gh-mjungmath]

Replying to tscrim:

You have renamed the public function eval_poly, so this need a deprecated_function_alias.

Right!

As per comment:39, I consider this to be a regression because you removed _eval_special_values_ from _evalf_:

sage: gen_legendre_P(x,x,.2)                                                                              
gen_legendre_P(x, x, 0.200000000000000)

versus before:

(-0.960000000000000)^(1/2*x)*factorial(2*x)/(2^x*factorial(x))

This has nothing to do with removing _eval_special_values_. In fact:

sage: assume(x, 'integer')                                                      
sage: gen_legendre_P(x, x, .2)                                                  
0.960000000000000^(1/2*x)*(-1)^x*factorial(2*x)/(2^x*factorial(x))

It is the line

+        if m.is_integer() and n.is_integer():
             ...
             if m == n:

that makes sure that this formula is only used when x is an integer.

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​82a447c

Trac #25034: eval_poly deprecation + doc improved

[gh-mjungmath]

Here we go. Moreover, I have added some examples in the visible doc that show some floating number examples like the one above.

[tscrim]

Okay, thanks for the explanation. This will be good. Just a few more minor things:

• Case |m| > |n| for integers: needs a ::
• Add a period at the end of the .. MATH:: in eval_gen_poly.
• Some formatting: Evaluate the Ferrers function `P(n, m, x)` for `m` and `n` being concrete
• I would move the paragraph These expressions ... `where |m| \leq |n|`. to the main part of the documentation so it is more visible.

Only the first is a must-do, but the rest are good while-we-are-at-it things.

Once these are done, it will be a positive review from me.

[gh-DaveWitteMorris]

I think references should go in the master bibliography file (src/doc/en/reference/references/index.rst), instead of the docstring. See #27497.

[gh-mjungmath]

Replying to tscrim:

• I would move the paragraph These expressions ... `where |m| \leq |n|`. to the main part of the documentation so it is more visible.

I agree. However, the documentation of the whole file needs to be refactored imo (see #31636). Thus, I'll leave it there for now.

As for the rest: thank you for taking such a thorough look! :)

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​0cda09b

Trac #25034: formatting + typos

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​0b14d02

Trac #25034: move reference to master bib

[gh-mjungmath]

Replying to gh-DaveWitteMorris:

I think references should go in the master bibliography file (src/doc/en/reference/references/index.rst), instead of the docstring. See #27497.

Is that alright?

Branch pushed to git repo; I updated commit sha1. New commits:

​0b14d02

Trac #25034: move reference to master bib

[gh-mjungmath]

Patchbot is green. :)

[mkoeppe]

Looks like all sufficient conditions for positive review are satisfied.

[gh-mjungmath]

Thank you Matthias! :)