Skip to content

Commit

Permalink
Trac #24553: Make legendre_P() a GinacFunction
Browse files Browse the repository at this point in the history 
Pynac-0.7.16 offers `legendre_P()`. The performance increase:
{{{
                            Python+Pari        Pynac
                            -----------        -----
sage: _=legendre_P(10,x)      510 µs           29 µs
sage: _=legendre_P(100,x)     2.3 ms           290 µs
sage: _=legendre_P(1000,x)    28 ms            3.5 ms
sage: _=legendre_P(10000,x)   1.1 s            110 ms
sage: _=legendre_P(50000,x)   30 s             2.1 s
}}}
overall a speedup of a comfortable order of magnitude.

URL: https://trac.sagemath.org/24553
Reported by: rws
Ticket author(s): Ralf Stephan
Reviewer(s): Travis Scrimshaw
  • Loading branch information
@vbraun
Release Manager authored and vbraun committed Mar 4, 2018
2 parents bef891a + 6c0ef30 commit e188244
Showing 1 changed file with 63 additions and 135 deletions.
198 changes: 63 additions & 135 deletions src/sage/functions/orthogonal_polys.py
View file
Expand Up @@ -1134,7 +1134,69 @@ def _derivative_(self, n, x, diff_param):
chebyshev_U = Func_chebyshev_U()


class Func_legendre_P(BuiltinFunction):
class Func_legendre_P(GinacFunction):
r"""
EXAMPLES::
sage: legendre_P(4, 2.0)
55.3750000000000
sage: legendre_P(1, x)
x
sage: legendre_P(4, x+1)
35/8*(x + 1)^4 - 15/4*(x + 1)^2 + 3/8
sage: legendre_P(1/2, I+1.)
1.05338240025858 + 0.359890322109665*I
sage: legendre_P(0, SR(1)).parent()
Symbolic Ring
sage: legendre_P(0, 0)
1
sage: legendre_P(1, x)
x
sage: legendre_P(4, 2.)
55.3750000000000
sage: legendre_P(5.5,1.00001)
1.00017875754114
sage: legendre_P(1/2, I+1).n()
1.05338240025858 + 0.359890322109665*I
sage: legendre_P(1/2, I+1).n(59)
1.0533824002585801 + 0.35989032210966539*I
sage: legendre_P(42, RR(12345678))
2.66314881466753e309
sage: legendre_P(42, Reals(20)(12345678))
2.6632e309
sage: legendre_P(201/2, 0).n()
0.0561386178630179
sage: legendre_P(201/2, 0).n(100)
0.056138617863017877699963095883
sage: R.<x> = QQ[]
sage: legendre_P(4,x)
35/8*x^4 - 15/4*x^2 + 3/8
sage: legendre_P(10000,x).coefficient(x,1)
0
sage: var('t,x')
(t, x)
sage: legendre_P(-5,t)
35/8*t^4 - 15/4*t^2 + 3/8
sage: legendre_P(4, x+1)
35/8*(x + 1)^4 - 15/4*(x + 1)^2 + 3/8
sage: legendre_P(4, sqrt(2))
83/8
sage: legendre_P(4, I*e)
35/8*e^4 + 15/4*e^2 + 3/8
sage: n = var('n')
sage: derivative(legendre_P(n,x), x)
(n*x*legendre_P(n, x) - n*legendre_P(n - 1, x))/(x^2 - 1)
sage: derivative(legendre_P(3,x), x)
15/2*x^2 - 3/2
sage: derivative(legendre_P(n,x), n)
Traceback (most recent call last):
...
RuntimeError: derivative w.r.t. to the index is not supported yet
"""
def __init__(self):
r"""
Init method for the Legendre polynomials of the first kind.
Expand All @@ -1150,140 +1212,6 @@ def __init__(self):
'maple':'LegendreP',
'giac':'legendre'})

def _eval_(self, n, x, *args, **kwds):
r"""
Return an evaluation of this Legendre P expression.
EXAMPLES::
sage: legendre_P(4, 2.0)
55.3750000000000
sage: legendre_P(1, x)
x
sage: legendre_P(4, x+1)
35/8*(x + 1)^4 - 15/4*(x + 1)^2 + 3/8
sage: legendre_P(1/2, I+1.)
1.05338240025858 + 0.359890322109665*I
sage: legendre_P(0, SR(1)).parent()
Symbolic Ring
"""
ret = self._eval_special_values_(n, x)
if ret is not None:
return ret
if n in ZZ:
ret = self.eval_pari(n, x)
if ret is not None:
return ret

def _eval_special_values_(self, n, x):
"""
Special values known.
EXAMPLES::
sage: legendre_P(0, 0)
1
sage: legendre_P(1, x)
x
"""
if n == 0 or n == -1 or x == 1:
return ZZ(1)
if n == 1 or n == -2:
return x

def _evalf_(self, n, x, parent=None, **kwds):
"""
EXAMPLES::
sage: legendre_P(4, 2.)
55.3750000000000
sage: legendre_P(5.5,1.00001)
1.00017875754114
sage: legendre_P(1/2, I+1).n()
1.05338240025858 + 0.359890322109665*I
sage: legendre_P(1/2, I+1).n(59)
1.0533824002585801 + 0.35989032210966539*I
sage: legendre_P(42, RR(12345678))
2.66314881466753e309
sage: legendre_P(42, Reals(20)(12345678))
2.6632e309
sage: legendre_P(201/2, 0).n()
0.0561386178630179
sage: legendre_P(201/2, 0).n(100)
0.056138617863017877699963095883
"""
ret = self._eval_special_values_(n, x)
if ret is not None:
return ret

import mpmath
from sage.libs.mpmath.all import call as mpcall
return mpcall(mpmath.legenp, n, 0, x, parent=parent)

def eval_pari(self, n, arg, **kwds):
"""
Use Pari to evaluate legendre_P for integer, symbolic, and
polynomial argument.
EXAMPLES::
sage: R.<x> = QQ[]
sage: legendre_P(4,x)
35/8*x^4 - 15/4*x^2 + 3/8
sage: legendre_P(10000,x).coefficient(x,1)
0
sage: var('t,x')
(t, x)
sage: legendre_P(-5,t)
35/8*t^4 - 15/4*t^2 + 3/8
sage: legendre_P(4, x+1)
35/8*(x + 1)^4 - 15/4*(x + 1)^2 + 3/8
sage: legendre_P(4, sqrt(2))
83/8
sage: legendre_P(4, I*e)
35/8*e^4 + 15/4*e^2 + 3/8
"""
if n < 0:
n = - n - 1
P = parent(arg)
if P in (ZZ, QQ, RR, CC, SR):
from sage.libs.pari.all import pari
R = PolynomialRing(QQ, 'x')
pol = R(pari.pollegendre(n))
return sum(b * arg**a for a, b in enumerate(pol))
elif is_PolynomialRing(P):
from sage.libs.pari.all import pari
if arg == P.gen():
return P(pari.pollegendre(n))
else:
R = PolynomialRing(QQ, 'x')
pol = R(pari.pollegendre(n))
pol = pol.subs({pol.parent().gen():arg})
pol = pol.change_ring(P.base_ring())
return pol

def _derivative_(self, n, x, *args,**kwds):
"""
Return the derivative of legendre_P.
EXAMPLES::
sage: n = var('n')
sage: derivative(legendre_P(n,x), x)
(n*x*legendre_P(n, x) - n*legendre_P(n - 1, x))/(x^2 - 1)
sage: derivative(legendre_P(3,x), x)
15/2*x^2 - 3/2
sage: derivative(legendre_P(n,x), n)
Traceback (most recent call last):
...
NotImplementedError: Derivative w.r.t. to the index is not supported.
"""
diff_param = kwds['diff_param']
if diff_param == 0:
raise NotImplementedError("Derivative w.r.t. to the index is not supported.")
else:
return (n*legendre_P(n-1, x) - n*x*legendre_P(n, x))/(1 - x**2)

legendre_P = Func_legendre_P()

class Func_legendre_Q(BuiltinFunction):
Expand Down

0 comments on commit e188244

Please sign in to comment.