Coverage for local/lib/python2.7/site-packages/sage/functions/orthogonal_polys.py : 72%

r""" 

Orthogonal Polynomials 

- The Chebyshev polynomial of the first kind arises as a solution 

to the differential equation 

.. MATH:: 

(1-x^2)\,y'' - x\,y' + n^2\,y = 0 

and those of the second kind as a solution to 

.. MATH:: 

(1-x^2)\,y'' - 3x\,y' + n(n+2)\,y = 0. 

The Chebyshev polynomials of the first kind are defined by the 

recurrence relation 

.. MATH:: 

T_0(x) = 1 \, T_1(x) = x \, T_{n+1}(x) = 2xT_n(x) - T_{n-1}(x). \, 

The Chebyshev polynomials of the second kind are defined by the 

recurrence relation 

.. MATH:: 

U_0(x) = 1 \, U_1(x) = 2x \, U_{n+1}(x) = 2xU_n(x) - U_{n-1}(x). \, 

For integers `m,n`, they satisfy the orthogonality 

relations 

.. MATH:: 

\int_{-1}^1 T_n(x)T_m(x)\,\frac{dx}{\sqrt{1-x^2}} =\left\{ \begin{matrix} 0 &: n\ne m~~~~~\\ \pi &: n=m=0\\ \pi/2 &: n=m\ne 0 \end{matrix} \right. 

and 

.. MATH:: 

\int_{-1}^1 U_n(x)U_m(x)\sqrt{1-x^2}\,dx =\frac{\pi}{2}\delta_{m,n}. 

They are named after Pafnuty Chebyshev (alternative 

transliterations: Tchebyshef or Tschebyscheff). 

- The Hermite polynomials are defined either by 

.. MATH:: 

H_n(x)=(-1)^n e^{x^2/2}\frac{d^n}{dx^n}e^{-x^2/2} 

(the "probabilists' Hermite polynomials"), or by 

.. MATH:: 

H_n(x)=(-1)^n e^{x^2}\frac{d^n}{dx^n}e^{-x^2} 

(the "physicists' Hermite polynomials"). Sage (via Maxima) 

implements the latter flavor. These satisfy the orthogonality 

relation 

.. MATH:: 

\int_{-\infty}^\infty H_n(x)H_m(x)\,e^{-x^2}\,dx ={n!2^n}{\sqrt{\pi}}\delta_{nm} 

They are named in honor of Charles Hermite. 

- Each *Legendre polynomial* `P_n(x)` is an `n`-th degree polynomial. 

It may be expressed using Rodrigues' formula: 

.. MATH:: 

P_n(x) = (2^n n!)^{-1} {\frac{d^n}{dx^n} } \left[ (x^2 -1)^n \right]. 

These are solutions to Legendre's differential equation: 

.. MATH:: 

{\frac{d}{dx}} \left[ (1-x^2) {\frac{d}{dx}} P(x) \right] + n(n+1)P(x) = 0. 

and satisfy the orthogonality relation 

.. MATH:: 

\int_{-1}^{1} P_m(x) P_n(x)\,dx = {\frac{2}{2n + 1}} \delta_{mn} 

The *Legendre function of the second kind* `Q_n(x)` is another 

(linearly independent) solution to the Legendre differential equation. 

It is not an "orthogonal polynomial" however. 

The associated Legendre functions of the first kind 

`P_\ell^m(x)` can be given in terms of the "usual" 

Legendre polynomials by 

.. MATH:: 

\begin{array}{ll} P_\ell^m(x) &= (-1)^m(1-x^2)^{m/2}\frac{d^m}{dx^m}P_\ell(x) \\ &= \frac{(-1)^m}{2^\ell \ell!} (1-x^2)^{m/2}\frac{d^{\ell+m}}{dx^{\ell+m}}(x^2-1)^\ell. \end{array} 

Assuming `0 \le m \le \ell`, they satisfy the orthogonality 

relation: 

.. MATH:: 

\int_{-1}^{1} P_k ^{(m)} P_\ell ^{(m)} dx = \frac{2 (\ell+m)!}{(2\ell+1)(\ell-m)!}\ \delta _{k,\ell}, 

where `\delta _{k,\ell}` is the Kronecker delta. 

The associated Legendre functions of the second kind 

`Q_\ell^m(x)` can be given in terms of the "usual" 

Legendre polynomials by 

.. MATH:: 

Q_\ell^m(x) = (-1)^m(1-x^2)^{m/2}\frac{d^m}{dx^m}Q_\ell(x). 

They are named after Adrien-Marie Legendre. 

- Laguerre polynomials may be defined by the Rodrigues formula 

.. MATH:: 

L_n(x)=\frac{e^x}{n!}\frac{d^n}{dx^n}\left(e^{-x} x^n\right). 

They are solutions of Laguerre's equation: 

.. MATH:: 

x\,y'' + (1 - x)\,y' + n\,y = 0\, 

and satisfy the orthogonality relation 

.. MATH:: 

\int_0^\infty L_m(x) L_n(x) e^{-x}\,dx = \delta_{mn}. 

The generalized Laguerre polynomials may be defined by the 

Rodrigues formula: 

.. MATH:: 

L_n^{(\alpha)}(x) = {\frac{x^{-\alpha} e^x}{n!}}{\frac{d^n}{dx^n}} \left(e^{-x} x^{n+\alpha}\right) . 

(These are also sometimes called the associated Laguerre 

polynomials.) The simple Laguerre polynomials are recovered from 

the generalized polynomials by setting `\alpha =0`. 

They are named after Edmond Laguerre. 

- Jacobi polynomials are a class of orthogonal polynomials. They 

are obtained from hypergeometric series in cases where the series 

is in fact finite: 

.. MATH:: 

P_n^{(\alpha,\beta)}(z) =\frac{(\alpha+1)_n}{n!} \,_2F_1\left(-n,1+\alpha+\beta+n;\alpha+1;\frac{1-z}{2}\right) , 

where `()_n` is Pochhammer's symbol (for the rising 

factorial), (Abramowitz and Stegun p561.) and thus have the 

explicit expression 

.. MATH:: 

P_n^{(\alpha,\beta)} (z) = \frac{\Gamma (\alpha+n+1)}{n!\Gamma (\alpha+\beta+n+1)} \sum_{m=0}^n \binom{n}{m} \frac{\Gamma (\alpha + \beta + n + m + 1)}{\Gamma (\alpha + m + 1)} \left(\frac{z-1}{2}\right)^m . 

They are named after Carl Jacobi. 

- Ultraspherical or Gegenbauer polynomials are given in terms of 

the Jacobi polynomials `P_n^{(\alpha,\beta)}(x)` with 

`\alpha=\beta=a-1/2` by 

.. MATH:: 

C_n^{(a)}(x)= \frac{\Gamma(a+1/2)}{\Gamma(2a)}\frac{\Gamma(n+2a)}{\Gamma(n+a+1/2)} P_n^{(a-1/2,a-1/2)}(x). 

They satisfy the orthogonality relation 

.. MATH:: 

\int_{-1}^1(1-x^2)^{a-1/2}C_m^{(a)}(x)C_n^{(a)}(x)\, dx =\delta_{mn}2^{1-2a}\pi \frac{\Gamma(n+2a)}{(n+a)\Gamma^2(a)\Gamma(n+1)} , 

for `a>-1/2`. They are obtained from hypergeometric series 

in cases where the series is in fact finite: 

.. MATH:: 

C_n^{(a)}(z) =\frac{(2a)^{\underline{n}}}{n!} \,_2F_1\left(-n,2a+n;a+\frac{1}{2};\frac{1-z}{2}\right) 

where `\underline{n}` is the falling factorial. (See 

Abramowitz and Stegun p561) 

They are named for Leopold Gegenbauer (1849-1903). 

For completeness, the Pochhammer symbol, introduced by Leo August 

Pochhammer, `(x)_n`, is used in the theory of special 

functions to represent the "rising factorial" or "upper factorial" 

.. MATH:: 

(x)_n=x(x+1)(x+2)\cdots(x+n-1)=\frac{(x+n-1)!}{(x-1)!}. 

On the other hand, the "falling factorial" or "lower factorial" is 

.. MATH:: 

x^{\underline{n}}=\frac{x!}{(x-n)!} , 

in the notation of Ronald L. Graham, Donald E. Knuth and Oren 

Patashnik in their book Concrete Mathematics. 

.. TODO:: 

Implement Zernike polynomials. 

:wikipedia:`Zernike_polynomials` 

REFERENCES: 

- [AS1964]_ 

- :wikipedia:`Chebyshev_polynomials` 

- :wikipedia:`Legendre_polynomials` 

- :wikipedia:`Hermite_polynomials` 

- http://mathworld.wolfram.com/GegenbauerPolynomial.html 

- :wikipedia:`Jacobi_polynomials` 

- :wikipedia:`Laguerre_polynomia` 

- :wikipedia:`Associated_Legendre_polynomials` 

- [Koe1999]_ 

AUTHORS: 

- David Joyner (2006-06) 

- Stefan Reiterer (2010-) 

- Ralf Stephan (2015-) 

The original module wrapped some of the orthogonal/special functions 

in the Maxima package "orthopoly" and was written by Barton 

Willis of the University of Nebraska at Kearney. 

""" 

#***************************************************************************** 

# Copyright (C) 2006 William Stein <wstein@gmail.com> 

# 2006 David Joyner <wdj@usna.edu> 

# 2010 Stefan Reiterer <maldun.finsterschreck@gmail.com> 

# 

# Distributed under the terms of the GNU General Public License (GPL) 

# 

# This code is distributed in the hope that it will be useful, 

# but WITHOUT ANY WARRANTY; without even the implied warranty of 

# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 

# General Public License for more details. 

# 

# The full text of the GPL is available at: 

# 

# http://www.gnu.org/licenses/ 

#***************************************************************************** 

from six.moves import range 

import warnings 

from sage.misc.latex import latex 

from sage.misc.sage_eval import sage_eval 

from sage.rings.all import ZZ, QQ, RR, CC 

from sage.rings.polynomial.polynomial_element import Polynomial 

from sage.rings.polynomial.polynomial_ring import is_PolynomialRing 

from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing 

from sage.rings.real_mpfr import is_RealField 

from sage.rings.complex_field import is_ComplexField 

from sage.symbolic.ring import SR, is_SymbolicVariable 

from sage.symbolic.function import BuiltinFunction, GinacFunction 

from sage.symbolic.expression import Expression 

from sage.functions.other import factorial, binomial 

from sage.structure.all import parent 

class OrthogonalFunction(BuiltinFunction): 

""" 

Base class for orthogonal polynomials. 

This class is an abstract base class for all orthogonal polynomials since 

they share similar properties. The evaluation as a polynomial 

is either done via maxima, or with pynac. 

Convention: The first argument is always the order of the polynomial, 

the others are other values or parameters where the polynomial is 

evaluated. 

""" 

def __init__(self, name, nargs=2, latex_name=None, conversions={}): 

""" 

:class:`OrthogonalFunction` class needs the same input parameter as 

it's parent class. 

EXAMPLES:: 

sage: from sage.functions.orthogonal_polys import OrthogonalFunction 

sage: new = OrthogonalFunction('testo_P') 

sage: new 

testo_P 

""" 

try: 

self._maxima_name = conversions['maxima'] 

except KeyError: 

self._maxima_name = None 

super(OrthogonalFunction,self).__init__(name=name, nargs=nargs, 

latex_name=latex_name, conversions=conversions) 

def eval_formula(self, *args): 

""" 

Evaluate this polynomial using an explicit formula. 

EXAMPLES:: 

sage: from sage.functions.orthogonal_polys import OrthogonalFunction 

sage: P = OrthogonalFunction('testo_P') 

sage: P.eval_formula(1,2.0) 

Traceback (most recent call last): 

... 

NotImplementedError: no explicit calculation of values implemented 

""" 

raise NotImplementedError("no explicit calculation of values implemented") 

def _eval_special_values_(self, *args): 

""" 

Evaluate the polynomial explicitly for special values. 

EXAMPLES:: 

sage: var('n') 

n 

sage: chebyshev_T(n,-1) 

(-1)^n 

""" 

raise ValueError("no special values known") 

def _eval_(self, n, *args): 

""" 

The :meth:`_eval_()` method decides which evaluation suits best 

for the given input, and returns a proper value. 

EXAMPLES:: 

sage: var('n,x') 

(n, x) 

sage: chebyshev_T(5,x) 

16*x^5 - 20*x^3 + 5*x 

""" 

return None 

def __call__(self, *args, **kwds): 

""" 

This overides the call method from SageObject to avoid problems with coercions, 

since the _eval_ method is able to handle more data types than symbolic functions 

would normally allow. 

Thus we have the distinction between algebraic objects (if n is an integer), 

and else as symbolic function. 

EXAMPLES:: 

sage: chebyshev_T(5, x) 

16*x^5 - 20*x^3 + 5*x 

sage: chebyshev_T(5, x, algorithm='pari') 

16*x^5 - 20*x^3 + 5*x 

sage: chebyshev_T(5, x, algorithm='maxima') 

16*x^5 - 20*x^3 + 5*x 

sage: chebyshev_T(5, x, algorithm='recursive') 

16*x^5 - 20*x^3 + 5*x 

""" 

algorithm = kwds.get('algorithm', None) 

if algorithm == 'pari': 

return self.eval_pari(*args, **kwds) 

elif algorithm == 'recursive': 

return self.eval_recursive(*args, **kwds) 

elif algorithm == 'maxima': 

from sage.calculus.calculus import maxima 

kwds['hold'] = True 

return maxima(self._eval_(*args, **kwds))._sage_() 

return super(OrthogonalFunction,self).__call__(*args, **kwds) 

class ChebyshevFunction(OrthogonalFunction): 

""" 

Abstract base class for Chebyshev polynomials of the first and second kind. 

EXAMPLES:: 

sage: chebyshev_T(3,x) 

4*x^3 - 3*x 

""" 

def __call__(self, n, *args, **kwds): 

""" 

This overides the call method from SageObject to avoid problems with coercions, 

since the _eval_ method is able to handle more data types than symbolic functions 

would normally allow. 

Thus we have the distinction between algebraic objects (if n is an integer), 

and else as symbolic function. 

EXAMPLES:: 

sage: K.<a> = NumberField(x^3-x-1) 

sage: chebyshev_T(5, a) 

16*a^2 + a - 4 

sage: chebyshev_T(5,MatrixSpace(ZZ, 2)([1, 2, -4, 7])) 

[-40799 44162] 

[-88324 91687] 

sage: R.<x> = QQ[] 

sage: parent(chebyshev_T(5, x)) 

Univariate Polynomial Ring in x over Rational Field 

sage: chebyshev_T(5, 2, hold=True) 

chebyshev_T(5, 2) 

sage: chebyshev_T(1,2,3) 

Traceback (most recent call last): 

... 

TypeError: Symbolic function chebyshev_T takes exactly 2 arguments (3 given) 

""" 

# If n is an integer: consider the polynomial as an algebraic (not symbolic) object 

if n in ZZ and not kwds.get('hold', False): 

try: 

return self._eval_(n, *args) 

except Exception: 

pass 

return super(ChebyshevFunction,self).__call__(n, *args, **kwds) 

def _eval_(self, n, x): 

""" 

The :meth:`_eval_()` method decides which evaluation suits best 

for the given input, and returns a proper value. 

EXAMPLES:: 

sage: var('n,x') 

(n, x) 

sage: chebyshev_T(5,x) 

16*x^5 - 20*x^3 + 5*x 

sage: chebyshev_T(64, x) 

2*(2*(2*(2*(2*(2*x^2 - 1)^2 - 1)^2 - 1)^2 - 1)^2 - 1)^2 - 1 

sage: chebyshev_T(n,-1) 

(-1)^n 

sage: chebyshev_T(-7,x) 

64*x^7 - 112*x^5 + 56*x^3 - 7*x 

sage: chebyshev_T(3/2,x) 

chebyshev_T(3/2, x) 

sage: R.<t> = QQ[] 

sage: chebyshev_T(2,t) 

2*t^2 - 1 

sage: chebyshev_U(2,t) 

4*t^2 - 1 

sage: parent(chebyshev_T(4, RIF(5))) 

Real Interval Field with 53 bits of precision 

sage: RR2 = RealField(5) 

sage: chebyshev_T(100000,RR2(2)) 

8.9e57180 

sage: chebyshev_T(5,Qp(3)(2)) 

2 + 3^2 + 3^3 + 3^4 + 3^5 + O(3^20) 

sage: chebyshev_T(100001/2, 2) 

doctest:...: RuntimeWarning: mpmath failed, keeping expression unevaluated 

chebyshev_T(100001/2, 2) 

sage: chebyshev_U._eval_(1.5, Mod(8,9)) is None 

True 

""" 

# n is an integer => evaluate algebraically (as polynomial) 

if n in ZZ: 

n = ZZ(n) 

# Expanded symbolic expression only for small values of n 

if isinstance(x, Expression) and n.abs() < 32: 

return self.eval_formula(n, x) 

return self.eval_algebraic(n, x) 

if isinstance(x, Expression) or isinstance(n, Expression): 

# Check for known identities 

try: 

return self._eval_special_values_(n, x) 

except ValueError: 

# Don't evaluate => keep symbolic 

return None 

# n is not an integer and neither n nor x is symbolic. 

# We assume n and x are real/complex and evaluate numerically 

try: 

import sage.libs.mpmath.all as mpmath 

return self._evalf_(n, x) 

except mpmath.NoConvergence: 

warnings.warn("mpmath failed, keeping expression unevaluated", 

RuntimeWarning) 

return None 

except Exception: 

# Numerical evaluation failed => keep symbolic 

return None 

class Func_chebyshev_T(ChebyshevFunction): 

""" 

Chebyshev polynomials of the first kind. 

REFERENCE: 

- [AS1964]_ 22.5.31 page 778 and 6.1.22 page 256. 

EXAMPLES:: 

sage: chebyshev_T(5,x) 

16*x^5 - 20*x^3 + 5*x 

sage: var('k') 

k 

sage: test = chebyshev_T(k,x) 

sage: test 

chebyshev_T(k, x) 

""" 

def __init__(self): 

""" 

Init method for the chebyshev polynomials of the first kind. 

EXAMPLES:: 

sage: var('n, x') 

(n, x) 

sage: from sage.functions.orthogonal_polys import Func_chebyshev_T 

sage: chebyshev_T2 = Func_chebyshev_T() 

sage: chebyshev_T2(1,x) 

x 

sage: chebyshev_T(x, x)._sympy_() 

chebyshevt(x, x) 

sage: maxima(chebyshev_T(1,x, hold=True)) 

_SAGE_VAR_x 

sage: maxima(chebyshev_T(n, chebyshev_T(n, x))) 

chebyshev_t(_SAGE_VAR_n,chebyshev_t(_SAGE_VAR_n,_SAGE_VAR_x)) 

""" 

ChebyshevFunction.__init__(self, 'chebyshev_T', nargs=2, 

conversions=dict(maxima='chebyshev_t', 

mathematica='ChebyshevT', 

sympy='chebyshevt', 

giac='tchebyshev1')) 

def _latex_(self): 

r""" 

TESTS:: 

sage: latex(chebyshev_T) 

T_n 

""" 

return r"T_n" 

def _print_latex_(self, n, z): 

r""" 

TESTS:: 

sage: latex(chebyshev_T(3, x, hold=True)) 

T_{3}\left(x\right) 

""" 

return r"T_{{{}}}\left({}\right)".format(latex(n), latex(z)) 

def _eval_special_values_(self, n, x): 

""" 

Values known for special values of x. 

For details see [AS1964]_ 22.4 (p. 777) 

EXAMPLES: 

sage: var('n') 

n 

sage: chebyshev_T(n,1) 

1 

sage: chebyshev_T(n,0) 

1/2*(-1)^(1/2*n)*((-1)^n + 1) 

sage: chebyshev_T(n,-1) 

(-1)^n 

sage: chebyshev_T._eval_special_values_(3/2,x) 

Traceback (most recent call last): 

... 

ValueError: no special value found 

sage: chebyshev_T._eval_special_values_(n, 0.1) 

Traceback (most recent call last): 

... 

ValueError: no special value found 

""" 

if x == 1: 

return x 

if x == -1: 

return x**n 

if x == 0: 

return (1+(-1)**n)*(-1)**(n/2)/2 

raise ValueError("no special value found") 

def _evalf_(self, n, x, **kwds): 

""" 

Evaluates :class:`chebyshev_T` numerically with mpmath. 

EXAMPLES:: 

sage: chebyshev_T._evalf_(10,3) 

2.26195370000000e7 

sage: chebyshev_T._evalf_(10,3,parent=RealField(75)) 

2.261953700000000000000e7 

sage: chebyshev_T._evalf_(10,I) 

-3363.00000000000 

sage: chebyshev_T._evalf_(5,0.3) 

0.998880000000000 

sage: chebyshev_T(1/2, 0) 

0.707106781186548 

sage: chebyshev_T(1/2, 3/2) 

1.11803398874989 

sage: chebyshev_T._evalf_(1.5, Mod(8,9)) 

Traceback (most recent call last): 

... 

TypeError: cannot evaluate chebyshev_T with parent Ring of integers modulo 9 

This simply evaluates using :class:`RealField` or :class:`ComplexField`:: 

sage: chebyshev_T(1234.5, RDF(2.1)) 

5.48174256255782e735 

sage: chebyshev_T(1234.5, I) 

-1.21629397684152e472 - 1.21629397684152e472*I 

For large values of ``n``, mpmath fails (but the algebraic formula 

still works):: 

sage: chebyshev_T._evalf_(10^6, 0.1) 

Traceback (most recent call last): 

... 

NoConvergence: Hypergeometric series converges too slowly. Try increasing maxterms. 

sage: chebyshev_T(10^6, 0.1) 

0.636384327171504 

""" 

try: 

real_parent = kwds['parent'] 

except KeyError: 

real_parent = parent(x) 

if not is_RealField(real_parent) and not is_ComplexField(real_parent): 

# parent is not a real or complex field: figure out a good parent 

if x in RR: 

x = RR(x) 

real_parent = RR 

elif x in CC: 

x = CC(x) 

real_parent = CC 

if not is_RealField(real_parent) and not is_ComplexField(real_parent): 

raise TypeError("cannot evaluate chebyshev_T with parent {}".format(real_parent)) 

from sage.libs.mpmath.all import call as mpcall 

from sage.libs.mpmath.all import chebyt as mpchebyt 

return mpcall(mpchebyt, n, x, parent=real_parent) 

def eval_formula(self, n, x): 

""" 

Evaluate ``chebyshev_T`` using an explicit formula. 

See [AS1964]_ 227 (p. 782) for details for the recursions. 

See also [Koe1999]_ for fast evaluation techniques. 

INPUT: 

- ``n`` -- an integer 

- ``x`` -- a value to evaluate the polynomial at (this can be 

any ring element) 

EXAMPLES:: 

sage: chebyshev_T.eval_formula(-1,x) 

x 

sage: chebyshev_T.eval_formula(0,x) 

1 

sage: chebyshev_T.eval_formula(1,x) 

x 

sage: chebyshev_T.eval_formula(2,0.1) == chebyshev_T._evalf_(2,0.1) 

True 

sage: chebyshev_T.eval_formula(10,x) 

512*x^10 - 1280*x^8 + 1120*x^6 - 400*x^4 + 50*x^2 - 1 

sage: chebyshev_T.eval_algebraic(10,x).expand() 

512*x^10 - 1280*x^8 + 1120*x^6 - 400*x^4 + 50*x^2 - 1 

""" 

if n < 0: 

return self.eval_formula(-n, x) 

elif n == 0: 

return parent(x).one() 

res = parent(x).zero() 

for j in range(n // 2 + 1): 

f = factorial(n-1-j) / factorial(j) / factorial(n-2*j) 

res += (-1)**j * (2*x)**(n-2*j) * f 

res *= n/2 

return res 

def eval_algebraic(self, n, x): 

""" 

Evaluate :class:`chebyshev_T` as polynomial, using a recursive 

formula. 

INPUT: 

- ``n`` -- an integer 

- ``x`` -- a value to evaluate the polynomial at (this can be 

any ring element) 

EXAMPLES:: 

sage: chebyshev_T.eval_algebraic(5, x) 

2*(2*(2*x^2 - 1)*x - x)*(2*x^2 - 1) - x 

sage: chebyshev_T(-7, x) - chebyshev_T(7,x) 

0 

sage: R.<t> = ZZ[] 

sage: chebyshev_T.eval_algebraic(-1, t) 

t 

sage: chebyshev_T.eval_algebraic(0, t) 

1 

sage: chebyshev_T.eval_algebraic(1, t) 

t 

sage: chebyshev_T(7^100, 1/2) 

1/2 

sage: chebyshev_T(7^100, Mod(2,3)) 

2 

sage: n = 97; x = RIF(pi/2/n) 

sage: chebyshev_T(n, cos(x)).contains_zero() 

True 

sage: R.<t> = Zp(2, 8, 'capped-abs')[] 

sage: chebyshev_T(10^6+1, t) 

(2^7 + O(2^8))*t^5 + (O(2^8))*t^4 + (2^6 + O(2^8))*t^3 + (O(2^8))*t^2 + (1 + 2^6 + O(2^8))*t + (O(2^8)) 

""" 

if n == 0: 

return parent(x).one() 

if n < 0: 

return self._eval_recursive_(-n, x)[0] 

return self._eval_recursive_(n, x)[0] 

def _eval_recursive_(self, n, x, both=False): 

""" 

If ``both=True``, compute ``(T(n,x), T(n-1,x))`` using a 

recursive formula. 

If ``both=False``, return instead a tuple ``(T(n,x), False)``. 

EXAMPLES:: 

sage: chebyshev_T._eval_recursive_(5, x) 

(2*(2*(2*x^2 - 1)*x - x)*(2*x^2 - 1) - x, False) 

sage: chebyshev_T._eval_recursive_(5, x, True) 

(2*(2*(2*x^2 - 1)*x - x)*(2*x^2 - 1) - x, 2*(2*x^2 - 1)^2 - 1) 

""" 

if n == 1: 

return x, parent(x).one() 

assert n >= 2 

a, b = self._eval_recursive_((n+1)//2, x, both or n % 2) 

if n % 2 == 0: 

return 2*a*a - 1, both and 2*a*b - x 

else: 

return 2*a*b - x, both and 2*b*b - 1 

def _eval_numpy_(self, n, x): 

""" 

Evaluate ``self`` using numpy. 

EXAMPLES:: 

sage: import numpy 

sage: z = numpy.array([1,2]) 

sage: z2 = numpy.array([[1,2],[1,2]]) 

sage: z3 = numpy.array([1,2,3.]) 

sage: chebyshev_T(1,z) 

array([ 1., 2.]) 

sage: chebyshev_T(1,z2) 

array([[ 1., 2.], 

[ 1., 2.]]) 

sage: chebyshev_T(1,z3) 

array([ 1., 2., 3.]) 

sage: chebyshev_T(z,0.1) 

array([ 0.1 , -0.98]) 

""" 

from scipy.special import eval_chebyt 

return eval_chebyt(n, x) 

def _derivative_(self, n, x, diff_param): 

""" 

Return the derivative of :class:`chebyshev_T` in form of the Chebyshev 

polynomial of the second kind :class:`chebyshev_U`. 

EXAMPLES:: 

sage: var('k') 

k 

sage: derivative(chebyshev_T(k,x),x) 

k*chebyshev_U(k - 1, x) 

sage: derivative(chebyshev_T(3,x),x) 

12*x^2 - 3 

sage: derivative(chebyshev_T(k,x),k) 

Traceback (most recent call last): 

... 

NotImplementedError: derivative w.r.t. to the index is not supported yet 

""" 

if diff_param == 0: 

raise NotImplementedError("derivative w.r.t. to the index is not supported yet") 

elif diff_param == 1: 

return n*chebyshev_U(n-1, x) 

raise ValueError("illegal differentiation parameter {}".format(diff_param)) 

chebyshev_T = Func_chebyshev_T() 

class Func_chebyshev_U(ChebyshevFunction): 

""" 

Class for the Chebyshev polynomial of the second kind. 

REFERENCE: 

- [AS1964]_ 22.8.3 page 783 and 6.1.22 page 256. 

EXAMPLES::