Opened 7 years ago
Closed 5 years ago
#16813 closed defect (fixed)
symbolic Legendre / associated Legendre functions / polynomials
Reported by:  rws  Owned by:  

Priority:  major  Milestone:  sage7.5 
Component:  symbolics  Keywords:  orthogonal 
Cc:  Merged in:  
Authors:  Ralf Stephan, Stefan Reiterer  Reviewers:  Marc Mezzarobba, Travis Scrimshaw 
Report Upstream:  N/A  Work issues:  
Branch:  7238198 (Commits, GitHub, GitLab)  Commit:  7238198c59a7a042391198039e4c915c24fda86d 
Dependencies:  Stopgaps: 
Description (last modified by )
Defect because there is no Sage binding for a result returned by Maxima:
sage: hypergeometric([1/2,1/2],[2],4).simplify_hypergeometric() 1/2*I*sqrt(3)*assoc_legendre_p(1/2, 1, 5/3) sage: assoc_legendre_p ... NameError: name 'assoc_legendre_p' is not defined
Existing numeric functions are legendre_P
, legendre_Q
, gen_legendre_P
, and gen_legendre_Q
which correspond to P(n,x) / Q(n,x) and associated Legendre P(n,m,x) / Q(n,m,x).
They should be made symbolic. FLINT has fast code for P(n,x).
 PQ(n,x)
 PQ(n,m,x)
 functions PQ
See also http://ask.sagemath.org/question/27230/problemwithhypergeometric/
Change History (97)
comment:1 followup: ↓ 10 Changed 7 years ago by
comment:2 Changed 7 years ago by
Fine! See also http://trac.sagemath.org/wiki/symbolics/functions
comment:3 followup: ↓ 4 Changed 7 years ago by
OK I have a question. What is the equivalent recursive algorithm to https://github.com/sagemath/sage/blob/master/src/sage/functions/orthogonal_polys.py#L812834 for Legendre polynomials?
The link is valid as long #16812 is not merged.
comment:4 in reply to: ↑ 3 Changed 7 years ago by
Replying to rws:
OK I have a question. What is the equivalent recursive algorithm to https://github.com/sagemath/sage/blob/master/src/sage/functions/orthogonal_polys.py#L812834 for Legendre polynomials?
The link is valid as long #16812 is not merged.
Hi!
You won't have luck to find an equivalent recursion algorithm for Legendre Polynomials, since the recursion algorithm for Chebyshev Polynomials uses the fact that cheby polynomials are cosines in disguise, and thus one is able to build Cheby polyis in O(log N) time. For Legendre polynomials you have to use the classic recursion formula given in https://en.wikipedia.org/wiki/Legendre_polynomials#Recursive_definition
comment:5 followup: ↓ 12 Changed 7 years ago by
Maybe this could help you:
I already implemented all Orthopolys one time: http://trac.sagemath.org/attachment/ticket/9706/trac_9706_ortho_polys.patch
they only would need cleanup/restructuring. Maybe you can reuse some of the implemented methods (like recursions and derivatives)
comment:6 Changed 7 years ago by
Some timings for P(n,z):
sage: legendre_P(100,2.5) 6.39483750487443e66 sage: timeit('legendre_P(100,2.5)') 25 loops, best of 3: 21 ms per loop sage: from mpmath import legenp sage: legenp(100,0,2.5) mpf('6.3948375048744286e+66') sage: timeit('legenp(100,0,2.5)') 625 loops, best of 3: 97.2 µs per loop sage: from scipy.special import eval_legendre sage: eval_legendre(int(100),float(2.5)) 6.3948375048744324e+66 sage: timeit('eval_legendre(int(100),float(2.5))') 625 loops, best of 3: 7.62 µs per loop sage: eval_legendre(int(10^5),float(1.00001)) 3.1548483029540554e+192 sage: timeit('eval_legendre(int(10^6),float(2.5))') 25 loops, best of 3: 11.8 ms per loop sage: eval_legendre(int(10^6),float(2.5)) inf
while legenp
will already bail out at 10^{5} because of F convergence issues.
comment:7 Changed 7 years ago by
With P(n,x) symbolics and algebra, Pari is much better than Maxima
sage: P.<t> = QQ[] sage: timeit('legendre_P(1000,t)') 5 loops, best of 3: 2.8 s per loop sage: timeit('pari.pollegendre(1000,t)') 625 loops, best of 3: 366 µs per loop
comment:8 Changed 7 years ago by
 Branch set to u/rws/symbolic_legendre___associated_legendre_functions___polynomials
comment:9 Changed 7 years ago by
 Commit set to 0f86b77a9aa818add6848cacf9a3f371d49c7d3a
This is a proof of concept patch, and one can already use legendre_P
and see from that and the code how the other three functions will look like. So, now would be a good time for fundamental criticism 8)
New commits:
50da8c5  16813: skeleton P(n,x)

e74539b  16813: P(n,x) refined, documentation

0f86b77  16813: fixes for doctest failures

comment:10 in reply to: ↑ 1 ; followup: ↓ 11 Changed 7 years ago by
Replying to maldun:
A good read on Legendre polynomials is also the bible for ortho polys: Abramowitz and Stegun http://people.math.sfu.ca/~cbm/aands/page_331.htm
This appears outdated, it is replaced by http://dlmf.nist.gov/14
comment:11 in reply to: ↑ 10 Changed 7 years ago by
Replying to rws:
Replying to maldun:
A good read on Legendre polynomials is also the bible for ortho polys: Abramowitz and Stegun http://people.math.sfu.ca/~cbm/aands/page_331.htm
This appears outdated, it is replaced by http://dlmf.nist.gov/14
You can't call a source outdated, which still covers information that the newer source doesn't. I checked your link, and some things from A&S are missign e.g. explicit representation of Legendre Polynomials with their polynomial coefficients. And on another note: I don't see much harm in citing a classic work on this topic ...
comment:12 in reply to: ↑ 5 ; followup: ↓ 13 Changed 7 years ago by
Replying to maldun:
I already implemented all Orthopolys one time: http://trac.sagemath.org/attachment/ticket/9706/trac_9706_ortho_polys.patch
they only would need cleanup/restructuring. Maybe you can reuse some of the implemented methods (like recursions and derivatives)
I am not sure about the derivatives. For P(3,2,x).diff(x)
I get 45*x^2 + 15
(Wolfram agrees) while with your formula (lines 23772395 of the patch) I get (after simplification) 45*x^2  15
.
Update: what's your reference there?
comment:13 in reply to: ↑ 12 Changed 7 years ago by
Replying to rws:
Replying to maldun:
I already implemented all Orthopolys one time: http://trac.sagemath.org/attachment/ticket/9706/trac_9706_ortho_polys.patch
they only would need cleanup/restructuring. Maybe you can reuse some of the implemented methods (like recursions and derivatives)
I am not sure about the derivatives. For
P(3,2,x).diff(x)
I get45*x^2 + 15
(Wolfram agrees) while with your formula (lines 23772395 of the patch) I get (after simplification)45*x^2  15
.Update: what's your reference there?
It seems you are right. from GradshteynRyzhik p.1004 formula 8.7311 we have the relation
P(n,m,x).diff(x) = ((n+1m)*P(n+1,m,x)(n+1)*x*P(n,m,x))/(x**21)
The same relation holds for gen_legendre_Q
I suppose that's an copy/paste/rewrite mistake from my side.
comment:14 Changed 7 years ago by
Also, your recursive functions for Q(n,x)
and Q(n,m,x)
appear to be wrong:
sage: legendre_Q.eval_recursive(2,x).subs(x=3) 13/2*I*pi + 13/2*log(2)  9/2 sage: legendre_Q.eval_recursive(2,x).subs(x=3).n() 0.00545667363964419 + 20.4203522483337*I sage: legendre_Q(2,3.) 0.00545667363964451  20.4203522483337*I
The latter result from mpmath is supported by Wolfram.
As to Q(n,m,x)
:
sage: gen_legendre_Q(2,1,x).subs(x=3) 1/8*sqrt(2)*(72*I*pi + 72*log(2)  50) sage: gen_legendre_Q(2,1,x).subs(x=3).n() 39.9859464434253 + 0.0165114736149170*I sage: gen_legendre_Q(2,1,3.) 39.9859464434253 + 0.0165114736149193*I
Again, Wolfram supports the latter value from mpmath (symbolic as (25 i)/(2 sqrt(2))18 i sqrt(2) ((log(4))/2+1/2 (log(2)i pi))
).
comment:15 followup: ↓ 16 Changed 7 years ago by
OK, I resolved it by using conjugate()
on every logarithm in the Q(n,x)
algorithms (on which the Q(n,m,x)
recurrence is based, too).
Update: it however makes symbolic work tedious and differentiation impossible, at the moment.
See also https://groups.google.com/forum/?hl=en#!topic/sagesupport/bEMPMEYeZKU on derivatives of conjugates in Sage.
comment:16 in reply to: ↑ 15 ; followup: ↓ 17 Changed 7 years ago by
Thanks for resolving this issue! I suppose I wasn't careful enough with complex arguments. But in my defense: I hadn't time to test this codes well enough when I wrote them ... but hopefully they give some useful informations.
concerning complex conjugation: I hope my answer on the mailing list give some clues: https://groups.google.com/forum/?hl=en#!topic/sagesupport/bEMPMEYeZKU
Replying to rws:
OK, I resolved it by using
conjugate()
on every logarithm in theQ(n,x)
algorithms (on which theQ(n,m,x)
recurrence is based, too).Update: it however makes symbolic work tedious and differentiation impossible, at the moment.
See also https://groups.google.com/forum/?hl=en#!topic/sagesupport/bEMPMEYeZKU on derivatives of conjugates in Sage.
comment:17 in reply to: ↑ 16 Changed 7 years ago by
Replying to maldun:
Thanks for resolving this issue!
Unfortunately, while it would be easy to resolve this numerically, the instances of conjugate()
introduced in the recurrence will make symbolic results from the recurrence unreadable and, in case of derivatives, impossible to use. A different way of computing the recurrences is needed, one which does away with usage of conjugate()
.
comment:18 Changed 7 years ago by
Hi!
We have several possible ways out of this: 1) avoid recursion for symbolic argument for Legendre_Q and use another library (maxima, flint, sympy ... ) for evaluation. 2) Let it be, but avoid it as default method. 3) Maybe more elegantly: There is a more closed relation for legendre_Q (but it's not really a recursion):
Q(n,z) = ½P(n,z) ln((z+1)/(z1))  W(n1,z)
with
W(n,z) = Σ_{k=1}^n (1/k) P(k1,z) P(nk,z)
(GradshteynRyzhik p 1019f)
Maybe we could find an recursion for W(n,z)
Hope this could be of some use
Edit there should be a recursion since W(n,z) is the convolution of aseries of holonomic functions, and if I remeber correctly there is an theorem saying that convolutions of holonomic functions are also holonomic, thus should have a recursion.
Edit: The above mentioned relation can also be found here: http://people.math.sfu.ca/~cbm/aands/page_334.htm
comment:19 Changed 7 years ago by
sage: from ore_algebra import * sage: def W(n): return sum([(1/k)*legendre_P(k1,t)*legendre_P(nk,t) for k in range(1,n+1)]) ....: sage: R.<t> = QQ['t'] sage: guess([W(n) for n in range(1,10)], OreAlgebra(R['n'], 'Sn')) (n  3)*Sn^2 + (2*t*n + 5*t)*Sn  n  2
comment:20 Changed 7 years ago by
Nice! I didn't know that sage already supports Ore Algebras. It appears that my holonomic function package on Mathematica stopped to work.
comment:21 Changed 7 years ago by
I always need some time to figure out the final form (that offset of n
!) but: n*W_{n} = (2tnt)*W_{n1}  (n1)W_{n2} (W_{0}=0, W_{1}=1).
comment:22 Changed 7 years ago by
However, this will yield the same result unless the log
has conjugate
associated with it. This shows your recurrence is actually correct but numerical results derived from it by substitution may need a warning about the log branch. I know not enough about calculus, maybe I'll ask again, this time on sagedevel, if there should be a switch for log()
to select the branch in case of numerical evaluation. What do you think?
comment:23 Changed 7 years ago by
I think there must be another different formula, because Wolfram has this for Q(2,x)
:
sage: ((3*x^21)/2*(log(x+1)log(1x))/23*x/2).subs(x=3) 13/2*I*pi + 13/2*log(4)  13/2*log(2)  9/2 sage: ((3*x^21)/2*(log(x+1)log(1x))/23*x/2).subs(x=3).n() 0.00545667363964419  20.4203522483337*I
which yields the correct value without use of conjugate
.
The first few Q(n,x)
from Wolfram are:
Q(0,x) = 1/2 log(x+1)1/2 log(1x) Q(1,x) = x (1/2 log(x+1)1/2 log(1x))1 Q(2,x) = 1/2 (3 x^21) (1/2 log(x+1)1/2 log(1x))(3 x)/2 Q(3,x) = (5 x^2)/21/2 (35 x^2) x (1/2 log(x+1)1/2 log(1x))+2/3
which makes clear that instead of log((1+x)/(1x)).conjugate()
we should just use log(1+x)log(1x)
, of course. Oh well.
comment:24 followup: ↓ 26 Changed 7 years ago by
Oh yeah it's again the non uniqueness of the representation of the complex logarithm
sage: log((x+1)/(1x)).subs(x=3) I*pi + log(2) sage: (log(x+1)log(1x)).subs(x=3).simplify_log() I*pi + log(2) sage: log((x+1)/(1x)).subs(x=3).conjugate() I*pi + log(2)
confusing as hell ...
I think Wolfram uses the log(1+x)log(1x) representation simply by the fact that it is independent of the branch in the following sense: Let log(x) = lnx + i*arg(x) + 2kπi and log(y) = lny + i*arg(y) + 2kπi then
log(x)  log(y) = lnx + i*arg(x) + 2kπi lny + i*arg(y) + 2kπi = = lnx/y + i*(arg(x)  arg(y)) + 0
I.e. if we have the same branch on the logarithm the module of 2kπi cancels out.
That means the formula isn't exactly wrong, it uses simply a different branch of the logarithm. But the representation of log as difference saves us indeed a lot of trouble, and as showed above is independent of the branch we use.
Nevertheless I think we should stick to the recursion with W(n,x), because from a computational view it is a lot better since:
1) The computational complexity is the same (solving a two term recursion)
2) we save computation time since we don't have to simplify expressions containing logarithms but only polynomials which are much simpler to handle and expand.
comment:25 Changed 7 years ago by
 Commit changed from 0f86b77a9aa818add6848cacf9a3f371d49c7d3a to 74ca8ea0fd8ef0eab57ceb1405e0887efd971bc9
comment:26 in reply to: ↑ 24 ; followup: ↓ 27 Changed 7 years ago by
 Status changed from new to needs_review
Replying to maldun:
Nevertheless I think we should stick to the recursion with W(n,x), because from a computational view it is a lot better since:
1) The computational complexity is the same (solving a two term recursion)
2) we save computation time since we don't have to simplify expressions containing logarithms but only polynomials which are much simpler to handle and expand.
Well, I have implemented your recurrence using multivariate polynomials where the generator l
stands for the log
term and gets substituted subsequently. This is already twice as fast as Maxima. However, your intuition was right that the W(n,x)
formula is still faster, my guess because univariate polys are faster than multi. Note that the P(n,x)
have to be computed, too, but nevertheless it's about 10x the speed of Maxima (which BTW uses the wrong log branch as well).
I might add some introductory doc cleanup but the functions themselves are now finished. Please review.
comment:27 in reply to: ↑ 26 ; followup: ↓ 28 Changed 7 years ago by
Replying to rws:
Well, I have implemented your recurrence using multivariate polynomials where the generator
l
stands for thelog
term and gets substituted subsequently. This is already twice as fast as Maxima. However, your intuition was right that theW(n,x)
formula is still faster, my guess because univariate polys are faster than multi. Note that theP(n,x)
have to be computed, too, but nevertheless it's about 10x the speed of Maxima (which BTW uses the wrong log branch as well).I might add some introductory doc cleanup but the functions themselves are now finished. Please review.
If Maxima uses also this branch of the logarithm we should make sure that changing the branch of the logarithm does not interfere with existing code. Have you already testet the complete sage library with
sage testall
?
We should also ask on the mailing list if there are some objections with that.
Personally I'm fine with both, as long as it is consistent, since using another branch of the logarithm is not wrong, but maybe not expected. (Maybe I programmed the recursion that way, since I compared it with Maxima that time, so that the output does not change)
comment:28 in reply to: ↑ 27 ; followup: ↓ 29 Changed 7 years ago by
Replying to maldun:
Have you already testet the complete sage library with
sage testall?
Buildbot task is queued.
comment:29 in reply to: ↑ 28 ; followup: ↓ 31 Changed 7 years ago by
Replying to rws:
Replying to maldun:
Have you already testet the complete sage library with
sage testall?
Buildbot task is queued.
I asked on sagedevel if there are other objections on using the log(1+z)  log(1z) representation: https://groups.google.com/forum/?hl=en#!topic/sagedevel/5_4Pr8GypUA Let's see what the other developers think.
comment:30 Changed 7 years ago by
I'm afk for some time, will see when I get back.
comment:31 in reply to: ↑ 29 Changed 7 years ago by
comment:32 Changed 7 years ago by
Very good! Since there are no objections on sagedevel this ticket only needs review now.
comment:33 Changed 7 years ago by
 Commit changed from 74ca8ea0fd8ef0eab57ceb1405e0887efd971bc9 to daf47c5924ead67cea53c409ad37386095c8e49c
comment:34 Changed 6 years ago by
 Commit changed from daf47c5924ead67cea53c409ad37386095c8e49c to a3c7c1ee6d5043f6d741f96e258960a1438fb929
comment:35 Changed 6 years ago by
 Status changed from needs_review to needs_work
I'm getting several doctest failures like:
File "src/sage/functions/orthogonal_polys.py", line 1426, in sage.functions.orthogonal_polys.Func_legendre_Q._maxima_init_evaled_ Failed example: legendre_Q._maxima_init_evaled_(20,x).coeff(x^10) Expected: 29113619535/131072*log((x + 1)/(x  1)) Got: doctest:1: DeprecationWarning: coeff is deprecated. Please use coefficient instead. See http://trac.sagemath.org/17438 for details. 29113619535/131072*log((x + 1)/(x  1))
comment:36 Changed 6 years ago by
 Commit changed from a3c7c1ee6d5043f6d741f96e258960a1438fb929 to 4e99dfa0a5cc82d168b2d3e90a4c2d94f65a21d5
comment:37 Changed 6 years ago by
 Milestone changed from sage6.4 to sage6.6
 Status changed from needs_work to needs_review
comment:38 Changed 6 years ago by
Sorry if that's a stupid question, but why do you need to override __call__
? In any case I guess a comment explaining the reason would be useful.
comment:39 Changed 6 years ago by
sage: legendre_P(0, 1).n() ... AttributeError: 'int' object has no attribute 'n' sage: legendre_P(1, 1).n() 1.00000000000000
comment:40 Changed 6 years ago by
sage: legendre_P(42, 12345678) 1464081544112412716366892468459853695115358209756840823628776385121841706762162962766108364297738809627122469598911070029409071998031560780937314580877929546448067606814039101080853578172130265741673137107647826759931295833598386722228898959000304822089300102623241891719774710215869248045233381461475507593850687226480251/549755813888 sage: legendre_P(42, RR(12345678)) +infinity sage: legendre_P(42, Reals(100)(12345678)) 2.6631488146675341638308827323e309
→ Consistent so far. But:
sage: legendre_P(42, Reals(20)(12345678)) legendre_P(42, 1.2346e7)
comment:41 Changed 6 years ago by
Perhaps related to the previous comment, I'm no fan of the mechanism used to choose real_parent
. Do you have an example where this code would be useful that could not be handled at the level of Function
or perhaps OrthogonalPolynomial
?
comment:42 Changed 6 years ago by
More on the wishlist side of things: The numerical evaluation methods return nonsense in cases where Maple, say, is accurate.
sage: legendre_P(201/2, 0).n() 365146.687569733 sage: legendre_P(201/2, 0).n(100).n() 0.0561386178630179
comment:43 Changed 6 years ago by
 Status changed from needs_review to needs_work
comment:44 Changed 6 years ago by
sage: legendre_P(x, x, x) ... TypeError: Symbolic function legendre_P takes exactly 2 arguments (3 given)
but:
sage: legendre_P(1, x, x) x
comment:45 Changed 6 years ago by
Why does Func_legendre_P.__call__
contain:
elif algorithm == 'recursive': return self.eval_recursive(n, x)
while Func_legendre_P
has no method eval_recursive
?
comment:46 Changed 6 years ago by
 Commit changed from 4e99dfa0a5cc82d168b2d3e90a4c2d94f65a21d5 to 9df78de81e0c307a1181a6ad4453bc56a93870b2
comment:47 Changed 6 years ago by
 Status changed from needs_work to needs_review
This addresses all issues. Indeed the changes were necessary, thanks for your input.
If you set one of the ortho poly tickets to positive, then please wait with the other, so I can remove merge conflicts.
comment:48 Changed 6 years ago by
Thank you for fixing all this.
One more remark: the main module docstring needs updating (at least to remove the TODO about associated Legendre polynomials, perhaps also to clarify what relies on Maxima and what doesn't).
comment:49 followup: ↓ 63 Changed 6 years ago by
Also, like in the case of #17151, I think having:
sage: legendre_P(0, x).parent() Integer Ring sage: legendre_P(0, SR(1)).parent() Integer Ring
is a bug.
Compare for example:
sage: sin(0).parent() # not sure I like this, but somehow makes sense Integer Ring sage: sin(SR(0)).parent() Symbolic Ring
comment:50 Changed 6 years ago by
legendre_P
does not know that P(n, 1) == 1
.
comment:51 followup: ↓ 53 Changed 6 years ago by
I find it strange that Func_legendre_Q
and Func_assoc_legendre_{P,Q}
derive from OrthogonalPolynomial
since these functions are not polynomials, even for fixed integer m and n. And in fact the same argument could apply to Func_legendre_P
itself, since they seem to work for noninteger n, even though only the case of integer n appears to be documented.
comment:52 followups: ↓ 62 ↓ 67 Changed 6 years ago by
Not sure what to think about that:
sage: legendre_P(10, polygen(CC, 'x')) 46189/256*x^10  109395/256*x^8 + 45045/128*x^6  15015/128*x^4 + 3465/256*x^2  63/256 sage: legendre_P(10, polygen(CC, 'x'), algorithm='pari') 180.425781250000*x^10  427.324218750000*x^8 + 351.914062500000*x^6  117.304687500000*x^4 + 13.5351562500000*x^2  0.246093750000000 sage: legendre_P(10, polygen(CC, 'x'), coerce=False) ... TypeError: arguments must be symbolic expressions sage: legendre_P(10, polygen(CC, 'x'), algorithm='pari', coerce=False) 180.425781250000*x^10  427.324218750000*x^8 + 351.914062500000*x^6  117.304687500000*x^4 + 13.5351562500000*x^2  0.246093750000000 sage: legendre_P(10, polygen(RIF, 'x')) 46189/256*(x)^10  109395/256*(x)^8 + 45045/128*(x)^6  15015/128*(x)^4 + 3465/256*(x)^2  63/256
Is the idea that calls to legendre_P(n, x)
with no algorithm
keyword will always take the branch corresponding to x
in SR
, but the user can obtain faster evaluations at integers etc. by specifying algorithm='pari'
is they know what they are doing?
comment:53 in reply to: ↑ 51 ; followup: ↓ 55 Changed 6 years ago by
Replying to mmezzarobba:
I find it strange that
Func_assoc_legendre_{P,Q}
derive fromOrthogonalPolynomial
since these functions are not polynomials, even for fixed integer m and n.
https://en.wikipedia.org/wiki/Associated_Legendre_polynomials "...they satisfy the orthogonality condition..."
comment:54 Changed 6 years ago by
Of course, I'd rather have it all under functions/holonomic/
...
comment:55 in reply to: ↑ 53 Changed 6 years ago by
Replying to rws:
Replying to mmezzarobba:
I find it strange that
Func_assoc_legendre_{P,Q}
derive fromOrthogonalPolynomial
since these functions are not polynomials, even for fixed integer m and n.https://en.wikipedia.org/wiki/Associated_Legendre_polynomials "...they satisfy the orthogonality condition..."
Sure, but they are not polynomials. It is necessarily wrong to use the class OrthogonalPolynomial
for things that are not, strictly speaking, orthogonal polynomials if really necessary, but then there should at least be a clear warning in the documentation of that class.
Replying to rws:
Of course, I'd rather have it all under functions/holonomic/...
I'm not sure being holonomic or not makes a big difference for a symbolic function. What do you have in mind?
(As for me, I would like to have an implementation of general holonomic functions in Sage at some point, but I'm more thinking of a parent entirely separate from SR
and with conversion methods to and from symbolic functions. Symbolic functions that happen to be holonomic could have a method that returns their representation using a system of linear functional equations, and could use that to implement operations for which no better functionspecific code is available.)
comment:56 Changed 6 years ago by
sage: legendre_Q(1., 0.) (+infinity)*sqrt(pi)*sin(0.500000000000000*pi)
comment:57 Changed 6 years ago by
sage: legendre_Q(1, 1, algorithm='pari') ... AttributeError: 'Func_legendre_Q' object has no attribute 'eval_pari'
but
sage: legendre_Q(1, 1, algorithm=pari) ... TypeError: __call__() got an unexpected keyword argument 'algorithm'
comment:58 Changed 6 years ago by
sage: legendre_Q(1, x) ... UnboundLocalError: local variable 'help3' referenced before assignment
comment:59 Changed 6 years ago by
I find it surprising that:
sage: legendre_Q(0, x) 1/2*log(x + 1)  1/2*log(x + 1) sage: legendre_Q(0, pi) 1/2*log(pi + 1)  1/2*log(pi + 1)
but:
sage: legendre_Q(0, 2) legendre_Q(0, 2)
and even:
sage: legendre_Q(0, SR(2)) legendre_Q(0, 2)
comment:60 Changed 6 years ago by
A fun one:
sage: a = legendre_Q(0, 1/2) sage: a.n() 0.549306144334055 sage: maple(a).evalf() .54930614431.570796327*I
I think the issue is that Maple uses nonstandard branch cuts for the Legendre functions, so that the conversion to Maple is not correct.
comment:61 Changed 6 years ago by
 Status changed from needs_review to needs_work
comment:62 in reply to: ↑ 52 Changed 6 years ago by
Replying to mmezzarobba:
Is the idea that calls to
legendre_P(n, x)
with noalgorithm
keyword will always take the branch corresponding tox
inSR
, but the user can obtain faster evaluations at integers etc. by specifyingalgorithm='pari'
is they know what they are doing?
Effectively, one of the __call__
methods (at this point of time OrthogonalPolynomial.__call__
) dispatches to the resp. eval_...
method. If no match is found (and algorithm
is None
) __call__
returns to somewhere in the BuiltinFunction/Function
classes to either return a held object or call evalf
. Jeroen wanted to implement a default list of algorithm/eval_...
pairs checked automatically but, if this keeps coming up, I'll put it high on my list.
I'm not sure being holonomic or not makes a big difference for a symbolic function. What do you have in mind?
It's just a natural category for special functions, corresponding to a Dfinite series expansion.
comment:63 in reply to: ↑ 49 Changed 6 years ago by
 Dependencies set to #17953
Replying to mmezzarobba:
Also, like in the case of #17151, I think having:
sage: legendre_P(0, x).parent() Integer Ring sage: legendre_P(0, SR(1)).parent() Integer Ringis a bug.
This is now #17953.
comment:64 Changed 6 years ago by
 Description modified (diff)
comment:65 Changed 6 years ago by
From the ask.sagemath issue there is also this discrepancy:
sage: 1/2*sqrt(3)*gen_legendre_P(1/2, 1, 5/3) 0.483843755630126 + 0.369716687246133*I vs. Maple: A := z > exp(I*Pi*z)*hypergeom([z,1/2],[2],4): evalf(A(1/2)); 0.3697166872 + 0.4838437556 I
comment:66 Changed 6 years ago by
 Description modified (diff)
 Type changed from enhancement to defect
comment:67 in reply to: ↑ 52 Changed 6 years ago by
Replying to mmezzarobba:
Is the idea that calls to
legendre_P(n, x)
with noalgorithm
keyword will always take the branch corresponding tox
inSR
, but the user can obtain faster evaluations at integers etc. by specifyingalgorithm='pari'
is they know what they are doing?
No, there is no idea, I just haven't figured out why the BuiltinFunction
mechanism won't give me the original x
, so I klugded something in OrthogonalFunction.__call__
triggered by algorithm
. This will have to be fixed but not today.
comment:68 followup: ↓ 82 Changed 6 years ago by
@mezzarobba: The problem with nonnumeric arguments not getting through unconverted to function::eval
cannot be solved with the current BuiltinFunction
. Some numerics aren't even converted but throw an error, see #17790. I thus think that the given examples with polygen(CC)
as argument cannot be implemented at the moment without kludging a taylormade __call__
method, and should be worked around using substitution. A solution will only be possible with #18832.
comment:69 Changed 6 years ago by
 Commit changed from 9df78de81e0c307a1181a6ad4453bc56a93870b2 to 32cc796b408be2a530c8885e3e100525ba042390
Branch pushed to git repo; I updated commit sha1. New commits:
e437940  Merge branch 'develop' into t/16813/symbolic_legendre___associated_legendre_functions___polynomials

c3e0525  16813: brought introduction up to date

96e5e01  16813: rename base classes

32cc796  16813: no more algorithm kwd; small fixes

comment:70 Changed 6 years ago by
 Commit changed from 32cc796b408be2a530c8885e3e100525ba042390 to 21fc2aa4064e8cd9c5837c8c30887603e5d9b315
Branch pushed to git repo; I updated commit sha1. New commits:
21fc2aa  16813: more cases with Q(n,x)

comment:71 Changed 6 years ago by
 Status changed from needs_work to needs_review
comment:72 Changed 6 years ago by
 Milestone changed from sage6.6 to sage6.9
comment:73 Changed 6 years ago by
 Status changed from needs_review to needs_work
needs rebase, does not apply
comment:74 Changed 6 years ago by
 Commit changed from 21fc2aa4064e8cd9c5837c8c30887603e5d9b315 to 87488259454a950d4041cf72a1125e92b5d52cf5
Branch pushed to git repo; I updated commit sha1. New commits:
8748825  Merge branch 'develop' into t/16813/symbolic_legendre___associated_legendre_functions___polynomials

comment:75 Changed 6 years ago by
I have absolutely no idea why there is a merge conflict in orthogonal_polys.py
when no one did a change there in the develop branch. Will see if this continues.
comment:76 Changed 6 years ago by
 Status changed from needs_work to needs_review
comment:77 Changed 6 years ago by
 Branch changed from u/rws/symbolic_legendre___associated_legendre_functions___polynomials to u/rws/168131
comment:78 Changed 6 years ago by
 Commit changed from 87488259454a950d4041cf72a1125e92b5d52cf5 to f511634dda1469513a78e1a42c81c62df857d0f0
 Dependencies #17953 deleted
 Milestone changed from sage6.9 to sage6.10
Squashed it all into one commit.
New commits:
f511634  16813: symbolic Legendre / associated Legendre functions / polynomials

comment:79 Changed 5 years ago by
 Keywords orthogonal added
 Milestone changed from sage6.10 to sage7.2
 Status changed from needs_review to needs_work
some failing doctests, see bot report
comment:80 Changed 5 years ago by
 Commit changed from f511634dda1469513a78e1a42c81c62df857d0f0 to b70ad8d85172b1a3cf13ca0cb8ef6ba61f982077
comment:81 Changed 5 years ago by
 Dependencies set to #19464
The failure in symbolic/expression_conversions.py
is due to a random bug and we now depend on #19464.
comment:82 in reply to: ↑ 68 Changed 5 years ago by
Replying to rws:
@mezzarobba: The problem with nonnumeric arguments not getting through unconverted to
function::eval
cannot be solved with the currentBuiltinFunction
. Some numerics aren't even converted but throw an error, see #17790. I thus think that the given examples withpolygen(CC)
as argument cannot be implemented at the moment without kludging a taylormade__call__
method, and should be worked around using substitution. A solution will only be possible with #18832.
I have opened #20312 on which the above depends.
comment:83 Changed 5 years ago by
 Commit changed from b70ad8d85172b1a3cf13ca0cb8ef6ba61f982077 to 60b86493763a2914763fb1845cbd2fd0a41e116c
Branch pushed to git repo; I updated commit sha1. New commits:
60b8649  Merge branch 'develop' into t/16813/168131

comment:84 Changed 5 years ago by
 Status changed from needs_work to needs_review
comment:85 Changed 5 years ago by
Maybe it's better to factor out the associated function stuff, regardless if Python or C++.
comment:86 Changed 5 years ago by
 Reviewers set to Marc Mezzarobba
 Status changed from needs_review to positive_review
Time to get this merged, I guess. I'm still not 100% happy with how evaluation behaves (ideally, in the cases where the function is a polynomial, one would like to be able to evaluate it on elements of any ℚalgebra), but getting it right in all cases looks difficult, and imperfect code is more useful than no code at all :)
.
comment:87 Changed 5 years ago by
Agree and thanks for the review.
comment:88 Changed 5 years ago by
 Milestone changed from sage7.2 to sage7.5
 Status changed from positive_review to needs_work
Needs rebase.
comment:89 Changed 5 years ago by
 Commit changed from 60b86493763a2914763fb1845cbd2fd0a41e116c to 7238198c59a7a042391198039e4c915c24fda86d
comment:90 Changed 5 years ago by
 Status changed from needs_work to needs_review
Review needed for the last two commits.
comment:91 Changed 5 years ago by
 Reviewers changed from Marc Mezzarobba to Marc Mezzarobba, Travis Scrimshaw
 Status changed from needs_review to positive_review
comment:92 Changed 5 years ago by
 Status changed from positive_review to needs_work
One failing doctest, see patchbot.
comment:93 Changed 5 years ago by
 Status changed from needs_work to needs_review
Does not fail here. Note the patchbot is running 7.5beta1 not 2, and there were gegenbauer
changes with pynac0.7.0. I'll try to trigger the bots.
comment:94 Changed 5 years ago by
 Status changed from needs_review to positive_review
I ended up accidentally trying this on my beta1 as well. All is good.
comment:95 followup: ↓ 96 Changed 5 years ago by
The dependency is invalid/wontfix. What are you depending on?
comment:96 in reply to: ↑ 95 Changed 5 years ago by
 Dependencies #19464 deleted
Replying to vbraun:
The dependency is invalid/wontfix. What are you depending on?
Actually I don't recall and couldn't find out how floor(x,hold=True)
was fixed. That bug was the reason for the random_expr()
doctest fail in this and some other ticket. I'll remove the dependency.
comment:97 Changed 5 years ago by
 Branch changed from u/rws/168131 to 7238198c59a7a042391198039e4c915c24fda86d
 Resolution set to fixed
 Status changed from positive_review to closed
Hi!
Good to see that someone else is working on making the orthogonal polynomials symbolic, since my research interests shifted heavily in the past years.
A good read on Legendre polynomials is also the bible for ortho polys: Abramowitz and Stegun http://people.math.sfu.ca/~cbm/aands/page_331.htm
There you will find also much information on special values and other properties. I currently have some time left and can help with one thing or the other if you like,