#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:  sage9.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: 
Description (last modified by )
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 followup #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.30287619362453e17 + 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**21)**(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 sagedevel 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 sagedevel 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!) (1x^2)^(m/2)
See: https://groups.google.com/d/msg/sagedevel/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 4 years ago by
 Description modified (diff)
comment:2 Changed 4 years ago by
 Cc rws slelievre added
 Description modified (diff)
comment:3 Changed 4 years ago by
comment:4 followup: ↓ 5 Changed 3 years ago by
For reference, here are some related issues reported by users:
comment:5 in reply to: ↑ 4 ; followup: ↓ 6 Changed 3 years ago by
comment:6 in reply to: ↑ 5 Changed 2 years ago by
For reference, here are some related issues reported by users:
One more:
Yet one more:
comment:7 Changed 2 years ago by
 Cc fredrik.johansson added
 Description modified (diff)
 Milestone changed from sage8.2 to sage9.2
comment:8 Changed 22 months ago by
 Milestone changed from sage9.2 to sage9.3
comment:9 Changed 17 months ago by
 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 17 months ago by
 Cc tscrim added
comment:11 Changed 16 months ago by
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
comment:12 Changed 16 months ago by
I already worked on something. Will push tomorrow.
comment:13 Changed 16 months ago by
 Branch set to public/legendre_25034
comment:14 Changed 16 months ago by
 Commit set to 571632f25b4cf431062d4a867ee4bbd54ac15b14
comment:15 Changed 16 months ago by
 Commit changed from 571632f25b4cf431062d4a867ee4bbd54ac15b14 to 4fcc57437dcd43de0e6f32626b7522c0080b57a3
Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:
4fcc574  Trac #25034: formula for n=m fixed

comment:16 Changed 16 months ago by
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*(1x**2)**(m/2) * ZZ(2*m1).multifactorial(2) + return (1)**m*factorial(2*m)/(2**m*factorial(m)) * (1x**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.
FollowUps
 As for resolving convention conflicts, I propose to wrap this up in a followup 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 followup ticket (#31636).
comment:17 followup: ↓ 18 Changed 16 months ago by
Thank you so much for tackling this!
Shall the ticket be set to needs_review? (in order to draw it to the patchbots' attention)
comment:18 in reply to: ↑ 17 ; followup: ↓ 21 Changed 16 months ago by
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.
comment:19 Changed 16 months ago by
comment:20 Changed 16 months ago by
 Commit changed from 4fcc57437dcd43de0e6f32626b7522c0080b57a3 to a4edb24ed3ffdd3aa667d3200452464963c137e0
Branch pushed to git repo; I updated commit sha1. New commits:
a4edb24  Trac #25034: abs(m) exceeding n gives zero immediately

comment:21 in reply to: ↑ 18 Changed 16 months ago by
 Status changed from new to needs_review
Replying to ghmjungmath:
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 followup: #31639.
comment:22 Changed 16 months ago by
Patchbot is green at least.
comment:23 Changed 16 months ago by
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)) * (1x**2)**(m/2)
I feel like
(1)**m * factorial(2*m)/factorial(m)/2**m * (1x**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 whileweareatit:
return ZZ(0) +return ZZ.zero()
I probably would also be good to add support for negative n
by replacing it with n1
. 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 16 months ago by
 Commit changed from a4edb24ed3ffdd3aa667d3200452464963c137e0 to 6b365d7722273cb8f0c3c0bc66de195c67509bb9
comment:25 Changed 16 months ago by
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.
comment:26 Changed 16 months ago by
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 16 months ago by
Patchbot green again.
comment:28 Changed 16 months ago by
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 16 months ago by
Now it is grouped into "integer" and "noninteger" case. Unfortunately, this didn't work well with _eval_special_values_
because it would lead to many repeated ifstatements. 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 followup: ↓ 31 Changed 16 months ago by
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 16 months ago by
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 ifstatements, as said above. I don't think this is desirable either.
comment:32 Changed 16 months ago by
 Commit changed from 6b365d7722273cb8f0c3c0bc66de195c67509bb9 to 801298c0384e772e3e4ffa52e865b7ae12a97ba7
Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:
801298c  Trac #25034: allow negative n

comment:33 Changed 16 months ago by
What about this?
comment:34 Changed 16 months ago by
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 followup ticket.
What do you think?
comment:35 followup: ↓ 36 Changed 16 months ago by
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 nonnegative integer, negative integer values for `m` + with `m \leq n` can be obtained by: + + .. MATH:: + + P^{m}_n(x) = (1)^{m} \frac{(nm)!}{(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 16 months ago by
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 symbolicm
input as well.
I'll do.
comment:37 Changed 16 months ago by
 Description modified (diff)
comment:38 Changed 16 months ago by
 Status changed from needs_review to needs_work
comment:39 Changed 16 months ago by
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 16 months ago by
 Commit changed from 801298c0384e772e3e4ffa52e865b7ae12a97ba7 to 0b7d917305b77064de0ee01282edd732fc8bd02e
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

comment:41 Changed 16 months ago by
 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.
comment:42 Changed 16 months ago by
 Commit changed from 0b7d917305b77064de0ee01282edd732fc8bd02e to d203066273fa50d46c13e6bae2f1809d639ac505
Branch pushed to git repo; I updated commit sha1. New commits:
d203066  Trac #25034: wrong indent in docstring

comment:43 followup: ↓ 44 Changed 16 months ago by
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 1e14
? 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 16 months ago by
 Reviewers set to Eric Gourgoulhon, Travis Scrimshaw
Replying to egourgoulhon:
In the doctests Check whether :trac:
25034
yields correct results compared to Maxima inspecial.py
, shouldn't there be markers# abs tol 1e14
? Same remark forsage: gen_legendre_P(5/3,3,1.+I) 0.238163715352606 + 0.0548443903534220*Iin the docstring of
Func_assoc_legendre_P
.
Done.
comment:45 Changed 16 months ago by
 Commit changed from d203066273fa50d46c13e6bae2f1809d639ac505 to c746088e3d32b359fbacb500bf29fd35a3e33bed
Branch pushed to git repo; I updated commit sha1. New commits:
c746088  Trac #25034: add abs tol

comment:46 Changed 16 months ago by
I hope this will make it in Sage 9.3.
comment:47 Changed 16 months ago by
 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 16 months ago by
 Priority changed from critical to blocker
comment:49 followup: ↓ 50 Changed 16 months ago by
 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 16 months ago by
Replying to tscrim:
You have renamed the public function
eval_poly
, so this need adeprecated_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 16 months ago by
 Commit changed from c746088e3d32b359fbacb500bf29fd35a3e33bed to 82a447c9bea6e082fbc1ff734ae9c8e5f1de5e35
Branch pushed to git repo; I updated commit sha1. New commits:
82a447c  Trac #25034: eval_poly deprecation + doc improved

comment:52 Changed 16 months ago by
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 16 months ago by
 Status changed from needs_work to needs_review
comment:54 followup: ↓ 56 Changed 16 months ago by
Okay, thanks for the explanation. This will be good. Just a few more minor things:
Case
m > nfor integers:
needs a::
 Add a period at the end of the
.. MATH::
ineval_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 mustdo, but the rest are good whileweareatit things.
Once these are done, it will be a positive review from me.
comment:55 followup: ↓ 59 Changed 16 months ago by
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 16 months ago by
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 16 months ago by
 Commit changed from 82a447c9bea6e082fbc1ff734ae9c8e5f1de5e35 to 0cda09b4d306be36eda9e634c44a07ca71a38526
Branch pushed to git repo; I updated commit sha1. New commits:
0cda09b  Trac #25034: formatting + typos

comment:58 Changed 16 months ago by
 Commit changed from 0cda09b4d306be36eda9e634c44a07ca71a38526 to 0b14d02cbccf9d46149589acba7ee05dc2f12397
Branch pushed to git repo; I updated commit sha1. New commits:
0b14d02  Trac #25034: move reference to master bib

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

comment:60 Changed 16 months ago by
Patchbot is green. :)
comment:61 Changed 16 months ago by
 Status changed from needs_review to positive_review
Looks like all sufficient conditions for positive review are satisfied.
comment:62 Changed 16 months ago by
 Branch changed from public/legendre_25034 to 0b14d02cbccf9d46149589acba7ee05dc2f12397
 Resolution set to fixed
 Status changed from positive_review to closed
comment:63 Changed 16 months ago by
 Commit 0b14d02cbccf9d46149589acba7ee05dc2f12397 deleted
Thank you Matthias! :)
A (quite severe IMHO) consequence of this bug is
which is plain wrong: the term
sqrt(cos(theta)^21)
(which is imaginary fortheta
real) should actually besin(theta)
.