Opened 3 years ago

Closed 8 weeks ago

Last modified 7 weeks ago

#25034 closed defect (fixed)

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

Reported by: jcwomack Owned by:
Priority: blocker Milestone: sage-9.3
Component: misc Keywords: legendre, special function, spherical harmonic
Cc: rws, slelievre, fredrik.johansson, tscrim Merged in:
Authors: Michael Jung Reviewers: Eric Gourgoulhon, Travis Scrimshaw
Report Upstream: N/A Work issues:
Branch: 0b14d02 (Commits, GitHub, GitLab) Commit:
Dependencies: Stopgaps:

Status badges

Description (last modified by gh-mjungmath)

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.

Change History (63)

comment:1 Changed 3 years ago by jcwomack

  • Description modified (diff)

comment:2 Changed 3 years ago by slelievre

  • Cc rws slelievre added
  • Description modified (diff)

comment:3 Changed 3 years ago by 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).

comment:7 Changed 12 months ago by slelievre

  • Cc fredrik.johansson added
  • Description modified (diff)
  • Milestone changed from sage-8.2 to sage-9.2

comment:8 Changed 8 months ago by mkoeppe

  • Milestone changed from sage-9.2 to sage-9.3

comment:9 Changed 2 months ago by gh-mjungmath

  • Priority changed from major to critical

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

comment:10 Changed 2 months ago by gh-mjungmath

  • Cc tscrim added

comment:11 Changed 2 months ago by 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

Last edited 2 months ago by gh-mjungmath (previous) (diff)

comment:12 Changed 2 months ago by gh-mjungmath

I already worked on something. Will push tomorrow.

comment:13 Changed 2 months ago by gh-mjungmath

  • Branch set to public/legendre_25034

comment:14 Changed 2 months ago by git

  • Commit set to 571632f25b4cf431062d4a867ee4bbd54ac15b14

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

568d0faTrac #25034: follow mathematica convention on legendre polynomials
6d6ead1Trac #25034: doctest for spherical harmonics added
571632fTrac #25034: doctest further improved

comment:15 Changed 2 months ago by git

  • Commit changed from 571632f25b4cf431062d4a867ee4bbd54ac15b14 to 4fcc57437dcd43de0e6f32626b7522c0080b57a3

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

4fcc574Trac #25034: formula for n=m fixed

comment:16 Changed 2 months ago by 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).
Last edited 2 months ago by gh-mjungmath (previous) (diff)

comment:17 follow-up: Changed 2 months ago by 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)

Version 0, edited 2 months ago by egourgoulhon (next)

comment:18 in reply to: ↑ 17 ; follow-up: Changed 2 months ago by 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.

Last edited 2 months ago by gh-mjungmath (previous) (diff)

comment:19 Changed 2 months ago by gh-mjungmath

  • Authors set to Michael Jung

comment:20 Changed 2 months ago by git

  • Commit changed from 4fcc57437dcd43de0e6f32626b7522c0080b57a3 to a4edb24ed3ffdd3aa667d3200452464963c137e0

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

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

comment:21 in reply to: ↑ 18 Changed 2 months ago by gh-mjungmath

  • Status changed from new to needs_review

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.

comment:22 Changed 2 months ago by gh-mjungmath

Patchbot is green at least.

comment:23 Changed 2 months ago by 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).

comment:24 Changed 2 months ago by git

  • Commit changed from a4edb24ed3ffdd3aa667d3200452464963c137e0 to 6b365d7722273cb8f0c3c0bc66de195c67509bb9

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

7c532a1Trac #25034: minor improvements
6b365d7Trac #25034: slight refactoring of _eval_

comment:25 Changed 2 months ago by 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.

Last edited 2 months ago by gh-mjungmath (previous) (diff)

comment:26 Changed 2 months ago by 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.

comment:27 Changed 2 months ago by gh-mjungmath

Patchbot green again.

comment:28 Changed 2 months ago by 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.

comment:29 Changed 2 months ago by 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.

comment:30 follow-up: Changed 2 months ago by 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.

comment:31 in reply to: ↑ 30 Changed 2 months ago by 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.

Last edited 2 months ago by gh-mjungmath (previous) (diff)

comment:32 Changed 2 months ago by git

  • Commit changed from 6b365d7722273cb8f0c3c0bc66de195c67509bb9 to 801298c0384e772e3e4ffa52e865b7ae12a97ba7

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

801298cTrac #25034: allow negative n

comment:33 Changed 2 months ago by gh-mjungmath

What about this?

comment:34 Changed 2 months ago by 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?

comment:35 follow-up: Changed 2 months ago by 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.

comment:36 in reply to: ↑ 35 Changed 2 months ago by 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.

Last edited 2 months ago by gh-mjungmath (previous) (diff)

comment:37 Changed 2 months ago by gh-mjungmath

  • Description modified (diff)

comment:38 Changed 2 months ago by gh-mjungmath

  • Status changed from needs_review to needs_work

comment:39 Changed 2 months ago by 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.

comment:40 Changed 2 months ago by git

  • Commit changed from 801298c0384e772e3e4ffa52e865b7ae12a97ba7 to 0b7d917305b77064de0ee01282edd732fc8bd02e

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

0b7d917Trac #25034: Turn gen_legendre_P into Ferrers functions only

comment:41 Changed 2 months ago by gh-mjungmath

  • Status changed from needs_work to needs_review

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.

Last edited 2 months ago by gh-mjungmath (previous) (diff)

comment:42 Changed 2 months ago by git

  • Commit changed from 0b7d917305b77064de0ee01282edd732fc8bd02e to d203066273fa50d46c13e6bae2f1809d639ac505

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

d203066Trac #25034: wrong indent in docstring

comment:43 follow-up: Changed 2 months ago by 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.

comment:44 in reply to: ↑ 43 Changed 2 months ago by gh-mjungmath

  • Reviewers set to Eric Gourgoulhon, Travis Scrimshaw

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.

comment:45 Changed 2 months ago by git

  • Commit changed from d203066273fa50d46c13e6bae2f1809d639ac505 to c746088e3d32b359fbacb500bf29fd35a3e33bed

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

c746088Trac #25034: add abs tol

comment:46 Changed 2 months ago by gh-mjungmath

I hope this will make it in Sage 9.3.

comment:47 Changed 2 months ago by egourgoulhon

  • Status changed from needs_review to positive_review

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

Travis, do you agree with the positive review?

comment:48 Changed 2 months ago by mkoeppe

  • Priority changed from critical to blocker

comment:49 follow-up: Changed 2 months ago by tscrim

  • Status changed from positive_review to needs_work

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

comment:50 in reply to: ↑ 49 Changed 2 months ago by 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.

comment:51 Changed 2 months ago by git

  • Commit changed from c746088e3d32b359fbacb500bf29fd35a3e33bed to 82a447c9bea6e082fbc1ff734ae9c8e5f1de5e35

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

82a447cTrac #25034: eval_poly deprecation + doc improved

comment:52 Changed 2 months ago by 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.

comment:53 Changed 2 months ago by gh-mjungmath

  • Status changed from needs_work to needs_review

comment:54 follow-up: Changed 2 months ago by 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.

comment:55 follow-up: Changed 2 months ago by 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.

comment:56 in reply to: ↑ 54 Changed 2 months ago by 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! :)

comment:57 Changed 2 months ago by git

  • Commit changed from 82a447c9bea6e082fbc1ff734ae9c8e5f1de5e35 to 0cda09b4d306be36eda9e634c44a07ca71a38526

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

0cda09bTrac #25034: formatting + typos

comment:58 Changed 2 months ago by git

  • Commit changed from 0cda09b4d306be36eda9e634c44a07ca71a38526 to 0b14d02cbccf9d46149589acba7ee05dc2f12397

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

0b14d02Trac #25034: move reference to master bib

comment:59 in reply to: ↑ 55 Changed 2 months ago by 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:

0b14d02Trac #25034: move reference to master bib

comment:60 Changed 2 months ago by gh-mjungmath

Patchbot is green. :)

comment:61 Changed 2 months ago by mkoeppe

  • Status changed from needs_review to positive_review

Looks like all sufficient conditions for positive review are satisfied.

comment:62 Changed 8 weeks ago by vbraun

  • Branch changed from public/legendre_25034 to 0b14d02cbccf9d46149589acba7ee05dc2f12397
  • Resolution set to fixed
  • Status changed from positive_review to closed

comment:63 Changed 7 weeks ago by gh-mjungmath

  • Commit 0b14d02cbccf9d46149589acba7ee05dc2f12397 deleted

Thank you Matthias! :)

Note: See TracTickets for help on using tickets.