Opened 8 years ago

Closed 6 years ago

# symbolic Legendre / associated Legendre functions / polynomials

Reported by: Owned by: Ralf Stephan major sage-7.5 symbolics orthogonal Ralf Stephan, Stefan Reiterer Marc Mezzarobba, Travis Scrimshaw N/A 7238198 7238198c59a7a042391198039e4c915c24fda86d

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

### comment:1 follow-up:  10 Changed 8 years ago by Stefan Reiterer

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,

### comment:3 follow-up:  4 Changed 8 years ago by Ralf Stephan

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#L812-834 for Legendre polynomials?

The link is valid as long #16812 is not merged.

### comment:4 in reply to:  3 Changed 8 years ago by Stefan Reiterer

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#L812-834 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 follow-up:  12 Changed 8 years ago by Stefan Reiterer

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)

Last edited 8 years ago by Stefan Reiterer (previous) (diff)

### comment:6 Changed 8 years ago by Ralf Stephan

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 105 because of F convergence issues.

### comment:7 Changed 8 years ago by Ralf Stephan

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 8 years ago by Ralf Stephan

Branch: → u/rws/symbolic_legendre___associated_legendre_functions___polynomials

### comment:9 Changed 8 years ago by Ralf Stephan

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 ; follow-up:  11 Changed 8 years ago by Ralf Stephan

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 8 years ago by Stefan Reiterer

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 ; follow-up:  13 Changed 8 years ago by Ralf Stephan

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 2377-2395 of the patch) I get (after simplification) `-45*x^2 - 15`.

Last edited 8 years ago by Ralf Stephan (previous) (diff)

### comment:13 in reply to:  12 Changed 8 years ago by Stefan Reiterer

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 2377-2395 of the patch) I get (after simplification) `-45*x^2 - 15`.

It seems you are right. from Gradshteyn-Ryzhik p.1004 formula 8.731-1 we have the relation

```P(n,m,x).diff(x) = ((n+1-m)*P(n+1,m,x)-(n+1)*x*P(n,m,x))/(x**2-1)
```

The same relation holds for gen_legendre_Q

I suppose that's an copy/paste/rewrite mistake from my side.

### comment:14 Changed 8 years ago by Ralf Stephan

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 follow-up:  16 Changed 8 years ago by Ralf Stephan

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.

Last edited 8 years ago by Ralf Stephan (previous) (diff)

### comment:16 in reply to:  15 ; follow-up:  17 Changed 8 years ago by Stefan Reiterer

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/sage-support/bEMPMEYeZKU

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.

### comment:17 in reply to:  16 Changed 8 years ago by Ralf Stephan

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 8 years ago by Stefan Reiterer

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)/(z-1)) - W(n-1,z)
```

with

```W(n,z) = Σ_{k=1}^n (1/k) P(k-1,z) P(n-k,z)
```

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

Last edited 8 years ago by Stefan Reiterer (previous) (diff)

### comment:19 Changed 8 years ago by Ralf Stephan

```sage: from ore_algebra import *
sage: def W(n):
return sum([(1/k)*legendre_P(k-1,t)*legendre_P(n-k,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 8 years ago by Stefan Reiterer

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 8 years ago by Ralf Stephan

I always need some time to figure out the final form (that offset of `n`!) but: n*Wn = (2tn-t)*Wn-1 - (n-1)Wn-2 (W0=0, W1=1).

### comment:22 Changed 8 years ago by Ralf Stephan

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 sage-devel, if there should be a switch for `log()` to select the branch in case of numerical evaluation. What do you think?

### comment:23 Changed 8 years ago by Ralf Stephan

I think there must be another different formula, because Wolfram has this for `Q(2,x)`:

```sage: ((3*x^2-1)/2*(log(x+1)-log(1-x))/2-3*x/2).subs(x=3)
-13/2*I*pi + 13/2*log(4) - 13/2*log(2) - 9/2
sage: ((3*x^2-1)/2*(log(x+1)-log(1-x))/2-3*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(1-x)
Q(1,x) = x (1/2 log(x+1)-1/2 log(1-x))-1
Q(2,x) = 1/2 (3 x^2-1) (1/2 log(x+1)-1/2 log(1-x))-(3 x)/2
Q(3,x) = -(5 x^2)/2-1/2 (3-5 x^2) x (1/2 log(x+1)-1/2 log(1-x))+2/3
```

which makes clear that instead of `log((1+x)/(1-x)).conjugate()` we should just use `log(1+x)-log(1-x)`, of course. Oh well.

### comment:24 follow-up:  26 Changed 8 years ago by Stefan Reiterer

Oh yeah it's again the non uniqueness of the representation of the complex logarithm

```sage: log((x+1)/(1-x)).subs(x=3)
I*pi + log(2)
sage: (log(x+1)-log(1-x)).subs(x=3).simplify_log()
-I*pi + log(2)
sage: log((x+1)/(1-x)).subs(x=3).conjugate()
-I*pi + log(2)
```

confusing as hell ...

I think Wolfram uses the log(1+x)-log(1-x) representation simply by the fact that it is independent of the branch in the following sense: Let log(x) = ln|x| + i*arg(x) + 2kπi and log(y) = ln|y| + i*arg(y) + 2kπi then

```log(x) - log(y) = ln|x| + i*arg(x) + 2kπi- ln|y| + i*arg(y) + 2kπi =
= ln|x/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 8 years ago by git

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

 ​0b67667 `16813: Q(n,x) implementation` ​06df7c8 `16813: implement P(n,m,x)` ​7bd4a70 `Merge branch 'develop' into t/16813/symbolic_legendre___associated_legendre_functions___polynomials` ​74ca8ea `16813: implement Q(n,m,x); fixes, doctests`

### comment:26 in reply to:  24 ; follow-up:  27 Changed 8 years ago by Ralf Stephan

Authors: → Ralf Stephan new → needs_review

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 ; follow-up:  28 Changed 8 years ago by Stefan Reiterer

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.

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 ; follow-up:  29 Changed 8 years ago by Ralf Stephan

Authors: Ralf Stephan → Ralf Stephan, Stefan Reiterer

Have you already testet the complete sage library with

```sage -testall
```

?

### comment:29 in reply to:  28 ; follow-up:  31 Changed 8 years ago by Stefan Reiterer

Have you already testet the complete sage library with

```sage -testall
```

?

I asked on sage-devel if there are other objections on using the log(1+z) - log(1-z) representation: https://groups.google.com/forum/?hl=en#!topic/sage-devel/5_4Pr8GypUA Let's see what the other developers think.

### comment:30 Changed 8 years ago by Ralf Stephan

I'm afk for some time, will see when I get back.

### comment:31 in reply to:  29 Changed 8 years ago by Ralf Stephan

Have you already testet the complete sage library with

```sage -testall
```

?

And tested successfully, see the top of the ticket.

### comment:32 Changed 8 years ago by Stefan Reiterer

Very good! Since there are no objections on sage-devel this ticket only needs review now.

### comment:33 Changed 8 years ago by git

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

 ​93915cc `Merge branch 'develop' into t/16813/symbolic_legendre___associated_legendre_functions___polynomials` ​daf47c5 `16813: adaptations due to new BuiltinFunction code`

### comment:34 Changed 8 years ago by git

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

 ​779ea3a `Merge branch 'develop' into t/16813/symbolic_legendre___associated_legendre_functions___polynomials` ​a3c7c1e `16813: cosmetics`

### comment:35 Changed 8 years ago by Marc Mezzarobba

Status: needs_review → 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:
See http://trac.sagemath.org/17438 for details.
-29113619535/131072*log(-(x + 1)/(x - 1))
```

### comment:36 Changed 8 years ago by git

Commit: a3c7c1ee6d5043f6d741f96e258960a1438fb929 → 4e99dfa0a5cc82d168b2d3e90a4c2d94f65a21d5

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

 ​a062018 `Merge branch 'develop' into t/16813/symbolic_legendre___associated_legendre_functions___polynomials` ​4e99dfa `16813: fix doctests`

### comment:37 Changed 8 years ago by Ralf Stephan

Milestone: sage-6.4 → sage-6.6 needs_work → needs_review

### comment:38 Changed 8 years ago by Marc Mezzarobba

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 8 years ago by Marc Mezzarobba

```sage: legendre_P(0, 1).n()
...
AttributeError: 'int' object has no attribute 'n'
sage: legendre_P(1, 1).n()
1.00000000000000
```

### comment:40 Changed 8 years ago by Marc Mezzarobba

```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 8 years ago by Marc Mezzarobba

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 8 years ago by Marc Mezzarobba

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 8 years ago by Ralf Stephan

Status: needs_review → needs_work

### comment:44 Changed 8 years ago by Marc Mezzarobba

```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 8 years ago by Marc Mezzarobba

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 8 years ago by git

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

 ​47f5fa0 `16813: remove __call__ methods; improve doctests` ​07814ad `Merge branch 'develop' into t/16813/symbolic_legendre___associated_legendre_functions___polynomials` ​2cc59c6 `16813: remove real_parent stuff` ​9df78de `16813: more doctests`

### comment:47 Changed 8 years ago by Ralf Stephan

Status: needs_work → 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 8 years ago by Marc Mezzarobba

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 follow-up:  63 Changed 8 years ago by Marc Mezzarobba

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 8 years ago by Marc Mezzarobba

`legendre_P` does not know that `P(n, 1) == 1`.

### comment:51 follow-up:  53 Changed 8 years ago by Marc Mezzarobba

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 non-integer n, even though only the case of integer n appears to be documented.

Last edited 8 years ago by Marc Mezzarobba (previous) (diff)

### comment:52 follow-ups:  62  67 Changed 8 years ago by Marc Mezzarobba

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?

Last edited 8 years ago by Marc Mezzarobba (previous) (diff)

### comment:53 in reply to:  51 ; follow-up:  55 Changed 8 years ago by Ralf Stephan

I find it strange that `Func_assoc_legendre_{P,Q}`derive from `OrthogonalPolynomial` 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 8 years ago by Ralf Stephan

Of course, I'd rather have it all under `functions/holonomic/`...

### comment:55 in reply to:  53 Changed 8 years ago by Marc Mezzarobba

I find it strange that `Func_assoc_legendre_{P,Q}`derive from `OrthogonalPolynomial` 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.

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 function-specific code is available.)

### comment:56 Changed 8 years ago by Marc Mezzarobba

```sage: legendre_Q(-1., 0.)
(+infinity)*sqrt(pi)*sin(0.500000000000000*pi)
```

### comment:57 Changed 8 years ago by Marc Mezzarobba

```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 8 years ago by Marc Mezzarobba

```sage: legendre_Q(-1, x)
...
UnboundLocalError: local variable 'help3' referenced before assignment
```

### comment:59 Changed 8 years ago by Marc Mezzarobba

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 8 years ago by Marc Mezzarobba

A fun one:

```sage: a = legendre_Q(0, 1/2)
sage: a.n()
0.549306144334055
sage: maple(a).evalf()
.5493061443-1.570796327*I
```

I think the issue is that Maple uses non-standard branch cuts for the Legendre functions, so that the conversion to Maple is not correct.

Last edited 8 years ago by Marc Mezzarobba (previous) (diff)

### comment:61 Changed 8 years ago by Marc Mezzarobba

Status: needs_review → needs_work

### comment:62 in reply to:  52 Changed 8 years ago by Ralf Stephan

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?

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 D-finite series expansion.

### comment:63 in reply to:  49 Changed 8 years ago by Ralf Stephan

Dependencies: → #17953

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.

This is now #17953.

### comment:64 Changed 7 years ago by Ralf Stephan

Description: modified (diff)

### comment:65 Changed 7 years ago by Ralf Stephan

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 7 years ago by Ralf Stephan

Description: modified (diff) enhancement → defect

### comment:67 in reply to:  52 Changed 7 years ago by Ralf Stephan

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?

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 follow-up:  82 Changed 7 years ago by Ralf Stephan

@mezzarobba: The problem with non-numeric 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 taylor-made `__call__` method, and should be worked around using substitution. A solution will only be possible with #18832.

### comment:69 Changed 7 years ago by git

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 7 years ago by git

Commit: 32cc796b408be2a530c8885e3e100525ba042390 → 21fc2aa4064e8cd9c5837c8c30887603e5d9b315

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

 ​21fc2aa `16813: more cases with Q(n,x)`

### comment:71 Changed 7 years ago by Ralf Stephan

Status: needs_work → needs_review

### comment:72 Changed 7 years ago by Ralf Stephan

Milestone: sage-6.6 → sage-6.9

### comment:73 Changed 7 years ago by Frédéric Chapoton

Status: needs_review → needs_work

needs rebase, does not apply

### comment:74 Changed 7 years ago by git

Commit: 21fc2aa4064e8cd9c5837c8c30887603e5d9b315 → 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 7 years ago by Ralf Stephan

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 7 years ago by Ralf Stephan

Status: needs_work → needs_review

### comment:77 Changed 7 years ago by Ralf Stephan

Branch: u/rws/symbolic_legendre___associated_legendre_functions___polynomials → u/rws/16813-1

### comment:78 Changed 7 years ago by Ralf Stephan

Commit: 87488259454a950d4041cf72a1125e92b5d52cf5 → f511634dda1469513a78e1a42c81c62df857d0f0 #17953 sage-6.9 → sage-6.10

Squashed it all into one commit.

New commits:

 ​f511634 `16813: symbolic Legendre / associated Legendre functions / polynomials`

### comment:79 Changed 7 years ago by Frédéric Chapoton

Keywords: orthogonal added sage-6.10 → sage-7.2 needs_review → needs_work

some failing doctests, see bot report

### comment:80 Changed 7 years ago by git

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

 ​fdf55be `Merge branch 'u/rws/16813-1' of git://trac.sagemath.org/sage into tmp08` ​b70ad8d `16813: fix doctests`

### comment:81 Changed 7 years ago by Ralf Stephan

Dependencies: → #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 7 years ago by Ralf Stephan

@mezzarobba: The problem with non-numeric 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 taylor-made `__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 6 years ago by git

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

 ​60b8649 `Merge branch 'develop' into t/16813/16813-1`

### comment:84 Changed 6 years ago by Ralf Stephan

Status: needs_work → needs_review

### comment:85 Changed 6 years ago by Ralf Stephan

Maybe it's better to factor out the associated function stuff, regardless if Python or C++.

### comment:86 Changed 6 years ago by Marc Mezzarobba

Reviewers: → Marc Mezzarobba needs_review → 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 6 years ago by Ralf Stephan

Agree and thanks for the review.

### comment:88 Changed 6 years ago by Travis Scrimshaw

Milestone: sage-7.2 → sage-7.5 positive_review → needs_work

Needs rebase.

### comment:89 Changed 6 years ago by git

Commit: 60b86493763a2914763fb1845cbd2fd0a41e116c → 7238198c59a7a042391198039e4c915c24fda86d

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

 ​cfccbb7 `Merge branch 'develop' into t/16813/16813-1` ​7238198 `fix doctest`

### comment:90 Changed 6 years ago by Ralf Stephan

Status: needs_work → needs_review

Review needed for the last two commits.

### comment:91 Changed 6 years ago by Travis Scrimshaw

Reviewers: Marc Mezzarobba → Marc Mezzarobba, Travis Scrimshaw needs_review → positive_review

### comment:92 Changed 6 years ago by Travis Scrimshaw

Status: positive_review → needs_work

One failing doctest, see patchbot.

### comment:93 Changed 6 years ago by Ralf Stephan

Status: needs_work → needs_review

Does not fail here. Note the patchbot is running 7.5beta1 not 2, and there were `gegenbauer` changes with pynac-0.7.0. I'll try to trigger the bots.

### comment:94 Changed 6 years ago by Travis Scrimshaw

Status: needs_review → positive_review

I ended up accidentally trying this on my beta1 as well. All is good.

### comment:95 follow-up:  96 Changed 6 years ago by Volker Braun

The dependency is invalid/wontfix. What are you depending on?

### comment:96 in reply to:  95 Changed 6 years ago by Ralf Stephan

Dependencies: #19464

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.