symbolic Legendre / associated Legendre functions / polynomials

[rwst]

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)
  • https://en.wikipedia.org/wiki/Legendre_polynomials
  • http://mathworld.wolfram.com/LegendrePolynomial.html
• PQ(n,m,x)
  • https://en.wikipedia.org/wiki/Associated_Legendre_polynomials
  • http://mathworld.wolfram.com/AssociatedLegendrePolynomial.html
• functions PQ
  • https://en.wikipedia.org/wiki/Legendre_function
  • http://mathworld.wolfram.com/LegendreFunctionoftheFirstKind.html
  • http://mathworld.wolfram.com/LegendreFunctionoftheSecondKind.html

See also http://ask.sagemath.org/question/27230/problem-with-hypergeometric/

Component: symbolics

Keywords: orthogonal

Author: Ralf Stephan, Stefan Reiterer

Branch/Commit: 7238198

Reviewer: Marc Mezzarobba, Travis Scrimshaw

Issue created by migration from https://trac.sagemath.org/ticket/16813

[sagetrac-maldun]

comment:1

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,

[rwst]

comment:2

Fine! See also http://trac.sagemath.org/wiki/symbolics/functions

[rwst]

comment:3

OK I have a question. What is the equivalent recursive algorithm to https://github.com/sagemath/sage-prod/blob/master/src/sage/functions/orthogonal_polys.py#L812-834
for Legendre polynomials?

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

[sagetrac-maldun]

comment:4

Replying to @rwst:

OK I have a question. What is the equivalent recursive algorithm to https://github.com/sagemath/sage-prod/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

[sagetrac-maldun]

comment:5

Maybe this could help you:

I already implemented all Orthopolys one time: https://github.com/sagemath/sage-prod/files/10650394/trac_9706_ortho_polys.patch.gz

they only would need cleanup/restructuring. Maybe you can reuse some of the implemented methods (like recursions and derivatives)

[rwst]

comment:6

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

[rwst]

comment:7

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

[rwst]

Branch: u/rws/symbolic_legendre___associated_legendre_functions___polynomials

[rwst]

Commit: 0f86b77

[rwst]

comment:9

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

[rwst]

comment:10

Replying to @sagetrac-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