Ticket #9706: 9706_review.patch
File 9706_review.patch, 46.9 KB (added by , 7 years ago) 


sage/functions/orthogonal_polys.py
# HG changeset patch # User Jeroen Demeyer <jdemeyer@cage.ugent.be> # Date 1386614457 3600 # Node ID 99ad4c29d89ab7d9004251ccdeb40de0068c272e # Parent 5ff57dc4416a2f1cfb84842f66711d8b59ea4da6 Symbolic Chebyshev polynomials: reviewer patch diff git a/sage/functions/orthogonal_polys.py b/sage/functions/orthogonal_polys.py
a b 12 12 13 13  The Chebyshev polynomial of the first kind arises as a solution 14 14 to the differential equation 15 15 16 16 .. math:: 17 17 18 (1x^2)\,y''  x\,y' + n^2\,y = 0 18 (1x^2)\,y''  x\,y' + n^2\,y = 0 19 19 20 20 21 21 and those of the second kind as a solution to 22 22 23 23 .. math:: 24 24 25 (1x^2)\,y''  3x\,y' + n(n+2)\,y = 0. 25 (1x^2)\,y''  3x\,y' + n(n+2)\,y = 0. 26 26 27 27 28 28 The Chebyshev polynomials of the first kind are defined by the 29 29 recurrence relation 30 30 31 31 .. math:: 32 32 33 T_0(x) = 1 \, T_1(x) = x \, T_{n+1}(x) = 2xT_n(x)  T_{n1}(x). \, 33 T_0(x) = 1 \, T_1(x) = x \, T_{n+1}(x) = 2xT_n(x)  T_{n1}(x). \, 34 34 35 35 36 36 The Chebyshev polynomials of the second kind are defined by the 37 37 recurrence relation 38 38 39 39 .. math:: 40 40 41 U_0(x) = 1 \, U_1(x) = 2x \, U_{n+1}(x) = 2xU_n(x)  U_{n1}(x). \, 41 U_0(x) = 1 \, U_1(x) = 2x \, U_{n+1}(x) = 2xU_n(x)  U_{n1}(x). \, 42 42 43 43 44 44 45 45 For integers `m,n`, they satisfy the orthogonality 46 46 relations 47 47 48 48 .. math:: 49 49 50 \int_{1}^1 T_n(x)T_m(x)\,\frac{dx}{\sqrt{1x^2}} =\left\{ \begin{matrix} 0 &: n\ne m~~~~~\\ \pi &: n=m=0\\ \pi/2 &: n=m\ne 0 \end{matrix} \right. 50 \int_{1}^1 T_n(x)T_m(x)\,\frac{dx}{\sqrt{1x^2}} =\left\{ \begin{matrix} 0 &: n\ne m~~~~~\\ \pi &: n=m=0\\ \pi/2 &: n=m\ne 0 \end{matrix} \right. 51 51 52 52 53 53 and 54 54 55 55 56 56 .. math:: 57 57 58 \int_{1}^1 U_n(x)U_m(x)\sqrt{1x^2}\,dx =\frac{\pi}{2}\delta_{m,n}. 58 \int_{1}^1 U_n(x)U_m(x)\sqrt{1x^2}\,dx =\frac{\pi}{2}\delta_{m,n}. 59 59 60 60 61 61 … … 63 63 transliterations: Tchebyshef or Tschebyscheff). 64 64 65 65  The Hermite polynomials are defined either by 66 66 67 67 .. math:: 68 68 69 H_n(x)=(1)^n e^{x^2/2}\frac{d^n}{dx^n}e^{x^2/2} 69 H_n(x)=(1)^n e^{x^2/2}\frac{d^n}{dx^n}e^{x^2/2} 70 70 71 71 72 72 (the "probabilists' Hermite polynomials"), or by 73 73 74 74 75 75 .. math:: 76 76 77 H_n(x)=(1)^n e^{x^2}\frac{d^n}{dx^n}e^{x^2} 77 H_n(x)=(1)^n e^{x^2}\frac{d^n}{dx^n}e^{x^2} 78 78 79 79 80 80 (the "physicists' Hermite polynomials"). Sage (via Maxima) 81 81 implements the latter flavor. These satisfy the orthogonality 82 82 relation 83 83 84 84 .. math:: 85 85 86 \int_{\infty}^\infty H_n(x)H_m(x)\,e^{x^2}\,dx ={n!2^n}{\sqrt{\pi}}\delta_{nm} 86 \int_{\infty}^\infty H_n(x)H_m(x)\,e^{x^2}\,dx ={n!2^n}{\sqrt{\pi}}\delta_{nm} 87 87 88 88 89 89 … … 105 105 and satisfy the orthogonality relation 106 106 107 107 .. math:: 108 108 109 109 \int_{1}^{1} P_m(x) P_n(x)\,dx = {\frac{2}{2n + 1}} \delta_{mn} 110 110 111 111 The *Legendre function of the second kind* `Q_n(x)` is another … … 118 118 119 119 .. math:: 120 120 121 \begin{array}{ll} P_\ell^m(x) &= (1)^m(1x^2)^{m/2}\frac{d^m}{dx^m}P_\ell(x) \\ &= \frac{(1)^m}{2^\ell \ell!} (1x^2)^{m/2}\frac{d^{\ell+m}}{dx^{\ell+m}}(x^21)^\ell. \end{array} 121 \begin{array}{ll} P_\ell^m(x) &= (1)^m(1x^2)^{m/2}\frac{d^m}{dx^m}P_\ell(x) \\ &= \frac{(1)^m}{2^\ell \ell!} (1x^2)^{m/2}\frac{d^{\ell+m}}{dx^{\ell+m}}(x^21)^\ell. \end{array} 122 122 123 123 124 124 Assuming `0 \le m \le \ell`, they satisfy the orthogonality 125 125 relation: 126 126 127 127 .. math:: 128 128 129 \int_{1}^{1} P_k ^{(m)} P_\ell ^{(m)} dx = \frac{2 (\ell+m)!}{(2\ell+1)(\ellm)!}\ \delta _{k,\ell}, 129 \int_{1}^{1} P_k ^{(m)} P_\ell ^{(m)} dx = \frac{2 (\ell+m)!}{(2\ell+1)(\ellm)!}\ \delta _{k,\ell}, 130 130 131 131 132 132 where `\delta _{k,\ell}` is the Kronecker delta. … … 135 135 `Q_\ell^m(x)` can be given in terms of the "usual" 136 136 Legendre polynomials by 137 137 138 138 139 139 .. math:: 140 140 141 Q_\ell^m(x) = (1)^m(1x^2)^{m/2}\frac{d^m}{dx^m}Q_\ell(x). 141 Q_\ell^m(x) = (1)^m(1x^2)^{m/2}\frac{d^m}{dx^m}Q_\ell(x). 142 142 143 143 144 144 145 145 They are named after AdrienMarie Legendre. 146 146 147 147  Laguerre polynomials may be defined by the Rodrigues formula 148 148 149 149 .. math:: 150 150 151 L_n(x)=\frac{e^x}{n!}\frac{d^n}{dx^n}\left(e^{x} x^n\right). 151 L_n(x)=\frac{e^x}{n!}\frac{d^n}{dx^n}\left(e^{x} x^n\right). 152 152 153 153 154 154 They are solutions of Laguerre's equation: 155 155 156 156 157 157 .. math:: 158 158 159 x\,y'' + (1  x)\,y' + n\,y = 0\, 159 x\,y'' + (1  x)\,y' + n\,y = 0\, 160 160 161 161 and satisfy the orthogonality relation 162 162 163 163 164 164 .. math:: 165 165 166 \int_0^\infty L_m(x) L_n(x) e^{x}\,dx = \delta_{mn}. 166 \int_0^\infty L_m(x) L_n(x) e^{x}\,dx = \delta_{mn}. 167 167 168 168 169 169 170 170 The generalized Laguerre polynomials may be defined by the 171 171 Rodrigues formula: 172 172 173 173 174 174 .. math:: 175 175 176 L_n^{(\alpha)}(x) = {\frac{x^{\alpha} e^x}{n!}}{\frac{d^n}{dx^n}} \left(e^{x} x^{n+\alpha}\right) . 176 L_n^{(\alpha)}(x) = {\frac{x^{\alpha} e^x}{n!}}{\frac{d^n}{dx^n}} \left(e^{x} x^{n+\alpha}\right) . 177 177 178 178 179 179 (These are also sometimes called the associated Laguerre … … 185 185  Jacobi polynomials are a class of orthogonal polynomials. They 186 186 are obtained from hypergeometric series in cases where the series 187 187 is in fact finite: 188 188 189 189 .. math:: 190 190 191 P_n^{(\alpha,\beta)}(z) =\frac{(\alpha+1)_n}{n!} \,_2F_1\left(n,1+\alpha+\beta+n;\alpha+1;\frac{1z}{2}\right) , 191 P_n^{(\alpha,\beta)}(z) =\frac{(\alpha+1)_n}{n!} \,_2F_1\left(n,1+\alpha+\beta+n;\alpha+1;\frac{1z}{2}\right) , 192 192 193 193 194 194 where `()_n` is Pochhammer's symbol (for the rising 195 195 factorial), (Abramowitz and Stegun p561.) and thus have the 196 196 explicit expression 197 197 198 198 199 199 .. math:: 200 200 201 P_n^{(\alpha,\beta)} (z) = \frac{\Gamma (\alpha+n+1)}{n!\Gamma (\alpha+\beta+n+1)} \sum_{m=0}^n {n\choose m} \frac{\Gamma (\alpha + \beta + n + m + 1)}{\Gamma (\alpha + m + 1)} \left(\frac{z1}{2}\right)^m . 201 P_n^{(\alpha,\beta)} (z) = \frac{\Gamma (\alpha+n+1)}{n!\Gamma (\alpha+\beta+n+1)} \sum_{m=0}^n {n\choose m} \frac{\Gamma (\alpha + \beta + n + m + 1)}{\Gamma (\alpha + m + 1)} \left(\frac{z1}{2}\right)^m . 202 202 203 203 204 204 … … 208 208 the Jacobi polynomials `P_n^{(\alpha,\beta)}(x)` with 209 209 `\alpha=\beta=a1/2` by 210 210 211 211 212 212 .. math:: 213 213 214 C_n^{(a)}(x)= \frac{\Gamma(a+1/2)}{\Gamma(2a)}\frac{\Gamma(n+2a)}{\Gamma(n+a+1/2)} P_n^{(a1/2,a1/2)}(x). 214 C_n^{(a)}(x)= \frac{\Gamma(a+1/2)}{\Gamma(2a)}\frac{\Gamma(n+2a)}{\Gamma(n+a+1/2)} P_n^{(a1/2,a1/2)}(x). 215 215 216 216 217 217 They satisfy the orthogonality relation 218 218 219 219 .. math:: 220 220 221 \int_{1}^1(1x^2)^{a1/2}C_m^{(a)}(x)C_n^{(a)}(x)\, dx =\delta_{mn}2^{12a}\pi \frac{\Gamma(n+2a)}{(n+a)\Gamma^2(a)\Gamma(n+1)} , 221 \int_{1}^1(1x^2)^{a1/2}C_m^{(a)}(x)C_n^{(a)}(x)\, dx =\delta_{mn}2^{12a}\pi \frac{\Gamma(n+2a)}{(n+a)\Gamma^2(a)\Gamma(n+1)} , 222 222 223 223 224 224 for `a>1/2`. They are obtained from hypergeometric series 225 225 in cases where the series is in fact finite: 226 226 227 227 228 228 .. math:: 229 229 230 C_n^{(a)}(z) =\frac{(2a)^{\underline{n}}}{n!} \,_2F_1\left(n,2a+n;a+\frac{1}{2};\frac{1z}{2}\right) 230 C_n^{(a)}(z) =\frac{(2a)^{\underline{n}}}{n!} \,_2F_1\left(n,2a+n;a+\frac{1}{2};\frac{1z}{2}\right) 231 231 232 232 233 233 where `\underline{n}` is the falling factorial. (See … … 242 242 243 243 .. math:: 244 244 245 (x)_n=x(x+1)(x+2)\cdots(x+n1)=\frac{(x+n1)!}{(x1)!}. 245 (x)_n=x(x+1)(x+2)\cdots(x+n1)=\frac{(x+n1)!}{(x1)!}. 246 246 247 247 248 248 On the other hand, the "falling factorial" or "lower factorial" is 249 249 250 250 .. math:: 251 251 252 x^{\underline{n}}=\frac{x!}{(xn)!} , 252 x^{\underline{n}}=\frac{x!}{(xn)!} , 253 253 254 254 255 255 in the notation of Ronald L. Graham, Donald E. Knuth and Oren … … 287 287 288 288 .. :wikipedia:`Associated_Legendre_polynomials` 289 289 290 .. [EffCheby] Wolfram Koepf: Effcient Computation of Chebyshev Polynomials 290 .. [EffCheby] Wolfram Koepf: Effcient Computation of Chebyshev Polynomials 291 291 in Computer Algebra 292 Computer Algebra Systems: A Practical Guide. 292 Computer Algebra Systems: A Practical Guide. 293 293 John Wiley, Chichester (1999): 7999. 294 294 295 295 AUTHORS: … … 318 318 import warnings 319 319 320 320 from sage.misc.sage_eval import sage_eval 321 from sage.rings.all import ZZ, RR, CC, RIF, CIF 321 from sage.rings.all import ZZ, RR, CC 322 from sage.rings.real_mpfr import is_RealField 323 from sage.rings.complex_field import is_ComplexField 322 324 from sage.calculus.calculus import maxima 323 325 324 326 325 from sage.symbolic.function import BuiltinFunction , GinacFunction, is_inexact327 from sage.symbolic.function import BuiltinFunction 326 328 from sage.symbolic.expression import is_Expression 327 import sage.functions.special 328 from sage.functions.special import MaximaFunction, meval 329 from sage.functions.other import floor, gamma, factorial, abs, binomial 330 from sage.functions.other import sqrt, conjugate 331 from sage.functions.trig import sin, cos 332 from sage.functions.log import ln 333 import sage.symbolic.expression as expression 334 from sage.structure.parent import Parent 329 from sage.functions.other import factorial, binomial 335 330 from sage.structure.coerce import parent 336 331 337 332 _done = False … … 340 335 Internal function which checks if Maxima has loaded the 341 336 "orthopoly" package. All functions using this in this 342 337 file should call this function first. 343 338 344 339 TEST: 345 340 346 The global starts ``False``:: 341 The global starts ``False``:: 347 342 348 343 sage: sage.functions.orthogonal_polys._done 349 344 False 350 345 351 346 Then after using one of these functions, it changes:: 352 347 353 sage: from sage.functions.orthogonal_polys import legendre_P 348 sage: from sage.functions.orthogonal_polys import legendre_P 354 349 sage: legendre_P(2,x) 355 350 3/2*(x  1)^2 + 3*x  2 356 351 sage: sage.functions.orthogonal_polys._done … … 372 367 _done = True 373 368 374 369 375 376 370 class OrthogonalPolynomial(BuiltinFunction): 377 371 """ 378 372 Base class for orthogonal polynomials. … … 381 375 they share similar properties. The evaluation as a polynomial 382 376 is either done via maxima, or with pynac. 383 377 384 Convention: The first argument is always the order of the polynomial, 385 the last one is always the value `x` where the polynomial is evaluated. 378 Convention: The first argument is always the order of the polynomial, 379 the others are other values or parameters where the polynomial is 380 evaluated. 386 381 """ 387 382 def __init__(self, name, nargs=2, latex_name=None, conversions={}): 388 383 """ 389 384 :class:`OrthogonalPolynomial` class needs the same input parameter as 390 385 it's parent class. 391 386 392 387 EXAMPLES:: 393 388 394 389 sage: from sage.functions.orthogonal_polys import OrthogonalPolynomial … … 397 392 testo_P 398 393 """ 399 394 try: 400 self._maxima_name = conversions['maxima'] 395 self._maxima_name = conversions['maxima'] 401 396 except KeyError: 402 397 self._maxima_name = None 403 398 404 super(OrthogonalPolynomial,self).__init__(name=name, nargs=nargs, 399 super(OrthogonalPolynomial,self).__init__(name=name, nargs=nargs, 405 400 latex_name=latex_name, conversions=conversions) 406 401 407 402 def _maxima_init_evaled_(self, *args): 408 403 r""" 409 404 Return a string which represents this function evaluated at 410 ``*args`` in Maxima. 411 412 In fact these are thought to be the old wrappers for the orthogonal 413 polynomials. They are used when the other evaluation methods fail, 414 or are not fast enough. It appears that direct computation 415 with pynac is in most cases faster than maxima. Maxima comes into 416 play when all other methods fail. 417 418 EXAMPLES:: 419 420 sage: chebyshev_T(3,x) 421 4*x^3  3*x 422 """ 423 return None 424 425 def _apply_formula_(self, *args): 426 """ 427 Method which uses the three term recursion of the polynomial, 428 or explicit formulas instead of maxima to evaluate the polynomial 429 efficiently, if the `x` argument is not a symbolic expression. 405 ``n, x`` in Maxima. 430 406 431 407 EXAMPLES:: 432 408 433 409 sage: from sage.functions.orthogonal_polys import OrthogonalPolynomial 434 sage: new = OrthogonalPolynomial('testo_P') 435 sage: new._apply_formula_(1,2.0) 410 sage: P = OrthogonalPolynomial('testo_P') 411 sage: P._maxima_init_evaled_(2, 5) is None 412 True 413 """ 414 return None 415 416 def eval_formula(self, *args): 417 """ 418 Evaluate this polynomial using an explicit formula. 419 420 EXAMPLES:: 421 422 sage: from sage.functions.orthogonal_polys import OrthogonalPolynomial 423 sage: P = OrthogonalPolynomial('testo_P') 424 sage: P.eval_formula(1,2.0) 436 425 Traceback (most recent call last): 437 426 ... 438 NotImplementedError: no recursivecalculation of values implemented427 NotImplementedError: no explicit calculation of values implemented 439 428 """ 440 raise NotImplementedError("no recursivecalculation of values implemented")429 raise NotImplementedError("no explicit calculation of values implemented") 441 430 442 def _eval_special_values_(self, *args):431 def _eval_special_values_(self, *args): 443 432 """ 444 433 Evaluate the polynomial explicitly for special values. 445 434 446 435 EXAMPLES:: 447 436 448 437 sage: var('n') 449 438 n 450 439 sage: chebyshev_T(n,1) … … 452 441 """ 453 442 raise ValueError("no special values known") 454 443 455 def _eval_(self, *args):444 def _eval_(self, n, *args): 456 445 """ 457 446 The _eval_ method decides which evaluation suits best 458 447 for the given input, and returns a proper value. 459 448 460 449 EXAMPLES:: 461 450 451 sage: var('n,x') 452 (n, x) 462 453 sage: chebyshev_T(5,x) 463 16*x^5  20*x^3 + 5*x 464 sage: var('n')465 n454 16*x^5  20*x^3 + 5*x 455 sage: chebyshev_T(64, x) 456 2*(2*(2*(2*(2*(2*x^2  1)^2  1)^2  1)^2  1)^2  1)^2  1 466 457 sage: chebyshev_T(n,1) 467 458 (1)^n 468 459 sage: chebyshev_T(7,x) 469 460 64*x^7  112*x^5 + 56*x^3  7*x 470 461 sage: chebyshev_T(3/2,x) 471 462 chebyshev_T(3/2, x) 472 sage: x = PolynomialRing(QQ, 'x').gen()473 sage: chebyshev_T(2, x)474 2* x^2  1475 sage: chebyshev_U(2, x)476 4* x^2  1463 sage: R.<t> = QQ[] 464 sage: chebyshev_T(2,t) 465 2*t^2  1 466 sage: chebyshev_U(2,t) 467 4*t^2  1 477 468 sage: parent(chebyshev_T(4, RIF(5))) 478 469 Real Interval Field with 53 bits of precision 479 470 sage: RR2 = RealField(5) … … 482 473 sage: chebyshev_T(5,Qp(3)(2)) 483 474 2 + 3^2 + 3^3 + 3^4 + 3^5 + O(3^20) 484 475 sage: chebyshev_T(100001/2, 2) 485 doctest: 500: RuntimeWarning: Warning: mpmath returns NoConvergence exception! Use other method instead.476 doctest:...: RuntimeWarning: mpmath failed, keeping expression unevaluated 486 477 chebyshev_T(100001/2, 2) 487 478 sage: chebyshev_U._eval_(1.5, Mod(8,9)) is None 488 479 True 489 480 """ 490 if not is_Expression(args[0]): 491 # If x is no expression and is inexact or n is not an integer > make numerical evaluation 492 if (not is_Expression(args[1])) and (is_inexact(args[1]) or not args[0] in ZZ): 493 try: 494 import sage.libs.mpmath.all as mpmath 495 return self._evalf_(*args) 496 except AttributeError: 497 pass 498 except mpmath.NoConvergence: 499 warnings.warn("Warning: mpmath returns NoConvergence exception! Use other method instead.", 500 RuntimeWarning) 501 except ValueError: 502 pass 481 args_is_symbolic = any(is_Expression(x) for x in args) 503 482 504 # n is not an integer and x is an expression > return symbolic expression. 505 if not args[0] in ZZ: 506 if is_Expression(args[1]): 507 return None 483 # n is an integer => evaluate algebraically (as polynomial) 484 if n in ZZ: 485 n = ZZ(n) 486 # Expanded symbolic expression only for small values of n 487 if args_is_symbolic and n.abs() < 32: 488 return self.eval_formula(n, *args) 489 else: 490 return self.eval_algebraic(n, *args) 508 491 509 # Check for known identities 510 try: 511 return self._eval_special_values_(*args) 512 except ValueError: 513 pass 514 515 #if negative indices are not specified 516 #in _eval_special_values only return symbolic 517 #value 518 if args[0] < 0 and args[0] in ZZ: 492 if args_is_symbolic or is_Expression(n): 493 # Check for known identities 494 try: 495 return self._eval_special_values_(n, *args) 496 except ValueError: 497 # Don't evaluate => keep symbolic 519 498 return None 520 499 521 if args[0] in ZZ: 522 try: 523 return self._apply_formula_(*args) 524 except NotImplementedError: 525 pass 526 527 if self._maxima_name is None: 500 # n is not an integer and neither n nor x is symbolic. 501 # We assume n and x are real/complex and evaluate numerically 502 try: 503 import sage.libs.mpmath.all as mpmath 504 return self._evalf_(n, *args) 505 except mpmath.NoConvergence: 506 warnings.warn("mpmath failed, keeping expression unevaluated", 507 RuntimeWarning) 508 return None 509 except StandardError: 510 # Numerical evaluation failed => keep symbolic 528 511 return None 529 512 530 if args[0] in ZZ: # use maxima as last resort 531 return self._old_maxima_(*args) 532 else: 533 return None 534 535 def __call__(self,*args,**kwds): 536 """ 513 def __call__(self, n, *args, **kwds): 514 """ 537 515 This overides the call method from SageObject to avoid problems with coercions, 538 516 since the _eval_ method is able to handle more data types than symbolic functions 539 517 would normally allow. 540 518 Thus we have the distinction between algebraic objects (if n is an integer), 541 519 and else as symbolic function. 542 543 EXAMPLES:: 544 545 sage: K.<a> = NumberField(x^3x1) 546 sage: chebyshev_T(5, a) 547 16*a^2 + a  4 520 521 EXAMPLES:: 522 523 sage: K.<a> = NumberField(x^3x1) 524 sage: chebyshev_T(5, a) 525 16*a^2 + a  4 548 526 sage: chebyshev_T(5,MatrixSpace(ZZ, 2)([1, 2, 4, 7])) 549 527 [40799 44162] 550 528 [88324 91687] 551 529 sage: R.<x> = QQ[] 552 530 sage: parent(chebyshev_T(5, x)) 553 531 Univariate Polynomial Ring in x over Rational Field 554 """ 555 if 'hold' not in kwds: 556 kwds['hold'] = False 557 if 'coerce' not in kwds: 558 kwds['coerce']=True 559 560 if args[0] in ZZ and kwds['hold'] is False: #check if n is in ZZ>consider polynomial as algebraic structure 561 return self._eval_(*args) # Let eval methode decide which is best 562 else: # Consider OrthogonalPolynomial as symbol 563 return super(OrthogonalPolynomial,self).__call__(*args,**kwds) 532 sage: chebyshev_T(5, 2, hold=True) 533 chebyshev_T(5, 2) 534 sage: chebyshev_T(1,2,3) 535 Traceback (most recent call last): 536 ... 537 TypeError: Symbolic function chebyshev_T takes exactly 2 arguments (3 given) 538 """ 539 # If n is an integer: consider the polynomial as an algebraic (not symbolic) object 540 if n in ZZ and not kwds.get('hold', False): 541 try: 542 return self._eval_(n, *args) 543 except StandardError: 544 pass 564 545 565 def _old_maxima_(self,*args): 566 """ 567 Method which holds the old maxima wrappers as last alternative. 568 It returns None per default, and it only needs to be implemented, 569 if it is necessary. 546 return super(OrthogonalPolynomial,self).__call__(n, *args, **kwds) 570 547 571 EXAMPLES:: 572 573 sage: chebyshev_T._old_maxima_(7,x) is None 574 True 575 """ 576 None 577 578 class Func_chebyshev_T(OrthogonalPolynomial): 548 class Func_chebyshev_T(OrthogonalPolynomial): 579 549 """ 580 550 Chebyshev polynomials of the first kind. 581 551 … … 609 579 super(Func_chebyshev_T,self).__init__("chebyshev_T", nargs=2, 610 580 conversions=dict(maxima='chebyshev_t', 611 581 mathematica='ChebyshevT')) 612 613 def _eval_special_values_(self, *args):582 583 def _eval_special_values_(self, n, x): 614 584 """ 615 585 Values known for special values of x. 616 586 For details see [ASHandbook]_ 22.4 (p. 777) … … 621 591 n 622 592 sage: chebyshev_T(n,1) 623 593 1 594 sage: chebyshev_T(n,0) 595 1/2*(1)^(1/2*n)*((1)^n + 1) 624 596 sage: chebyshev_T(n,1) 625 597 (1)^n 626 sage: chebyshev_T(7, x)  chebyshev_T(7,x)627 0628 598 sage: chebyshev_T._eval_special_values_(3/2,x) 629 599 Traceback (most recent call last): 630 600 ... 631 ValueError: No special values for non integral indices!601 ValueError: no special value found 632 602 sage: chebyshev_T._eval_special_values_(n, 0.1) 633 603 Traceback (most recent call last): 634 604 ... 635 ValueError: Value not found! 636 sage: chebyshev_T._eval_special_values_(26, Mod(9,9)) 637 Traceback (most recent call last): 638 ... 639 ValueError: Value not found! 605 ValueError: no special value found 640 606 """ 641 if (not is_Expression(args[0])) and (not args[0] in ZZ): 642 raise ValueError("No special values for non integral indices!") 643 644 if args[1] == 1: 645 return args[1] 646 647 if args[1] == 1: 648 return args[1]**args[0] 607 if x == 1: 608 return x 649 609 650 if (args[1] == 0 and args[1] in CC):651 return (1+(1)**args[0])*(1)**(args[0]/2)/2610 if x == 1: 611 return x**n 652 612 653 if args[0] < 0 and args[0] in ZZ:654 return self._eval_(args[0],args[1])613 if x == 0: 614 return (1+(1)**n)*(1)**(n/2)/2 655 615 656 raise ValueError(" Value not found!")657 658 def _evalf_(self, *args,**kwds):616 raise ValueError("no special value found") 617 618 def _evalf_(self, n, x, **kwds): 659 619 """ 660 620 Evaluates :class:`chebyshev_T` numerically with mpmath. 661 If the index is an integer we use the recursive formula since662 it is faster.663 621 664 622 EXAMPLES:: 665 623 666 sage: chebyshev_T(10,3).n(75) 624 sage: chebyshev_T._evalf_(10,3) 625 2.26195370000000e7 626 sage: chebyshev_T._evalf_(10,3,parent=RealField(75)) 667 627 2.261953700000000000000e7 668 sage: chebyshev_T (10,I).n()628 sage: chebyshev_T._evalf_(10,I) 669 629 3363.00000000000 670 sage: chebyshev_T (5,0.3).n()630 sage: chebyshev_T._evalf_(5,0.3) 671 631 0.998880000000000 672 632 sage: chebyshev_T(1/2, 0) 673 633 0.707106781186548 634 sage: chebyshev_T(1/2, 3/2) 635 1.11803398874989 674 636 sage: chebyshev_T._evalf_(1.5, Mod(8,9)) 675 637 Traceback (most recent call last): 676 638 ... 677 ValueError: No compatible type!639 TypeError: cannot evaluate chebyshev_T with parent Ring of integers modulo 9 678 640 641 This simply evaluates using :class:`RealField` or :class:`ComplexField`:: 642 643 sage: chebyshev_T(1234.5, RDF(2.1)) 644 5.48174256255782e735 645 sage: chebyshev_T(1234.5, I) 646 1.21629397684152e472  1.21629397684152e472*I 647 648 For large values of ``n``, mpmath fails (but the algebraic formula 649 still works):: 650 651 sage: chebyshev_T._evalf_(10^6, 0.1) 652 Traceback (most recent call last): 653 ... 654 NoConvergence: Hypergeometric series converges too slowly. Try increasing maxterms. 655 sage: chebyshev_T(10^6, 0.1) 656 0.636384327171504 679 657 """ 680 if args[0] in ZZ and args[0] >= 0:681 return self._cheb_recur_(*args)[0]682 683 658 try: 684 659 real_parent = kwds['parent'] 685 660 except KeyError: 686 real_parent = parent( args[1])661 real_parent = parent(x) 687 662 688 x_set = False 689 if hasattr(real_parent,"precision"): # Check if we have a data type with precision 690 x = args[1] 691 step_parent = real_parent 692 x_set = True 693 else: 694 if args[1] in RR: 695 x = RR(args[1]) 696 step_parent = RR 697 x_set = True 698 elif args[1] in CC: 699 x = CC(args[1]) 700 step_parent = CC 701 x_set = True 663 if not is_RealField(real_parent) and not is_ComplexField(real_parent): 664 # parent is not a real or complex field: figure out a good parent 665 if x in RR: 666 x = RR(x) 667 real_parent = RR 668 elif x in CC: 669 x = CC(x) 670 real_parent = CC 702 671 703 if not x_set:704 raise ValueError("No compatible type!")705 672 if not is_RealField(real_parent) and not is_ComplexField(real_parent): 673 raise TypeError("cannot evaluate chebyshev_T with parent %s"%real_parent) 674 706 675 from sage.libs.mpmath.all import call as mpcall 707 676 from sage.libs.mpmath.all import chebyt as mpchebyt 708 677 709 return mpcall(mpchebyt, args[0],x,parent=step_parent)678 return mpcall(mpchebyt, n, x, parent=real_parent) 710 679 711 def _maxima_init_evaled_(self, *args):680 def _maxima_init_evaled_(self, n, x): 712 681 """ 713 682 Evaluate the Chebyshev polynomial ``self`` with maxima. 714 683 715 684 EXAMPLES:: 716 685 686 sage: var('n, x') 687 (n, x) 717 688 sage: chebyshev_T._maxima_init_evaled_(1,x) 718 689 'x' 719 sage: var('n') 720 n 721 sage: maxima(chebyshev_T(n,x)) 722 chebyshev_t(n,x) 690 sage: maxima(chebyshev_T(n, chebyshev_T(n, x))) 691 chebyshev_t(n,chebyshev_t(n,x)) 723 692 724 693 """ 725 n = args[0] 726 x = args[1] 727 return maxima.eval('chebyshev_t({0},{1})'.format(n,x)) 694 return maxima.eval('chebyshev_t({0},{1})'.format(n._maxima_init_(), x._maxima_init_())) 728 695 729 730 def _apply_formula_(self,*args):696 697 def eval_formula(self, n, x): 731 698 """ 732 Applies explicit formulas for :class:`chebyshev_T`. 733 This is much faster for numerical evaluation than maxima! 699 Evaluate ``chebyshev_T`` using an explicit formula. 734 700 See [ASHandbook]_ 227 (p. 782) for details for the recurions. 735 701 See also [EffCheby]_ for fast evaluation techniques. 736 702 703 INPUT: 704 705  ``n``  an integer 706 707  ``x``  a value to evaluate the polynomial at (this can be 708 any ring element) 709 737 710 EXAMPLES:: 738 711 739 sage: chebyshev_T._apply_formula_(2,0.1) == chebyshev_T._evalf_(2,0.1) 712 sage: chebyshev_T.eval_formula(1,x) 713 x 714 sage: chebyshev_T.eval_formula(0,x) 715 1 716 sage: chebyshev_T.eval_formula(1,x) 717 x 718 sage: chebyshev_T.eval_formula(2,0.1) == chebyshev_T._evalf_(2,0.1) 740 719 True 741 sage: chebyshev_T (51,x)742 2*(2*(2*(2*(2*(2*x^2  1)^2  1)*(2*(2*x^2  1)*x  x)  x)*(2*(2*(2*x^2  1)*x  x)^2  1)  x)^2  1)*(2*(2*(2*(2*(2*x^2  1)^2  1)*(2*(2*x^2  1)*x  x)  x)*(2*(2*(2*x^2  1)*x  x)^2  1)  x)*(2*(2*(2*(2*x^2  1)*x  x)^2  1)^2  1)  x)  x743 sage: chebyshev_T. _apply_formula_(10,x)720 sage: chebyshev_T.eval_formula(10,x) 721 512*x^10  1280*x^8 + 1120*x^6  400*x^4 + 50*x^2  1 722 sage: chebyshev_T.eval_algebraic(10,x).expand() 744 723 512*x^10  1280*x^8 + 1120*x^6  400*x^4 + 50*x^2  1 745 724 746 725 """ 747 k = args[0] 748 x = args[1] 726 if n < 0: 727 return self.eval_formula(n, x) 728 elif n == 0: 729 return parent(x).one() 749 730 750 if k == 0: 751 return 1 752 if k == 1: 753 return x 754 755 help1 = 1 756 help2 = x 757 if is_Expression(x) and k <= 25: 758 # Recursion gives more compact representations for large k 759 help3 = 0 760 for j in xrange(0,floor(k/2)+1): 761 f = factorial(kj1) / factorial(j) / factorial(k2*j) 762 help3 = help3 + (1)**j * (2*x)**(k2*j) * f 763 help3 = help3 * k / 2 764 return help3 731 res = parent(x).zero() 732 for j in xrange(0, n//2+1): 733 f = factorial(n1j) / factorial(j) / factorial(n2*j) 734 res += (1)**j * (2*x)**(n2*j) * f 735 res *= n/2 736 return res 737 738 def eval_algebraic(self, n, x): 739 """ 740 Evaluate :class:`chebyshev_T` as polynomial, using a recursive 741 formula. 742 743 INPUT: 744 745  ``n``  an integer 746 747  ``x``  a value to evaluate the polynomial at (this can be 748 any ring element) 749 750 EXAMPLES:: 751 752 sage: chebyshev_T.eval_algebraic(5, x) 753 2*(2*(2*x^2  1)*x  x)*(2*x^2  1)  x 754 sage: chebyshev_T(7, x)  chebyshev_T(7,x) 755 0 756 sage: R.<t> = ZZ[] 757 sage: chebyshev_T.eval_algebraic(1, t) 758 t 759 sage: chebyshev_T.eval_algebraic(0, t) 760 1 761 sage: chebyshev_T.eval_algebraic(1, t) 762 t 763 sage: chebyshev_T(7^100, 1/2) 764 1/2 765 sage: chebyshev_T(7^100, Mod(2,3)) 766 2 767 sage: n = 97; x = RIF(pi/2/n) 768 sage: chebyshev_T(n, cos(x)).contains_zero() 769 True 770 sage: R.<t> = Zp(2, 8, 'cappedabs')[] 771 sage: chebyshev_T(10^6+1, t) 772 (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)) 773 """ 774 if n == 0: 775 return parent(x).one() 776 if n > 0: 777 return self._eval_recursive_(n, x)[0] 765 778 else: 766 return self._ cheb_recur_(k,x)[0]779 return self._eval_recursive_(n, x)[0] 767 780 768 def _cheb_recur_(self,n, x, both=False): 769 """ 770 Generalized recursion formula for Chebyshev polynomials. 771 Implementation suggested by Frederik Johansson. 772 returns (T(n,x), T(n1,x)), or (T(n,x), _) if both=False 773 774 EXAMPLES:: 775 776 sage: chebyshev_T._cheb_recur_(5,x) 777 (2*(2*(2*x^2  1)*x  x)*(2*x^2  1)  x, False) 778 """ 779 if n == 0: 780 return 1, x 781 if n == 1: 782 return x, 1 783 a, b = self._cheb_recur_((n+1)//2, x, both or n % 2) 784 if n % 2 == 0: 785 return 2*a**2  1, both and 2*a*b  x 786 else: 787 return 2*a*b  x, both and 2*b**2  1 781 def _eval_recursive_(self, n, x, both=False): 782 """ 783 If ``both=True``, compute `(T(n,x), T(n1,x))` using a 784 recursive formula. 785 If ``both=False``, return instead a tuple `(T(n,x), False)`. 788 786 787 EXAMPLES:: 789 788 790 def _eval_numpy_(self, *args): 789 sage: chebyshev_T._eval_recursive_(5, x) 790 (2*(2*(2*x^2  1)*x  x)*(2*x^2  1)  x, False) 791 sage: chebyshev_T._eval_recursive_(5, x, True) 792 (2*(2*(2*x^2  1)*x  x)*(2*x^2  1)  x, 2*(2*x^2  1)^2  1) 793 """ 794 if n == 1: 795 return x, parent(x).one() 796 797 assert n >= 2 798 a, b = self._eval_recursive_((n+1)//2, x, both or n % 2) 799 if n % 2 == 0: 800 return 2*a*a  1, both and 2*a*b  x 801 else: 802 return 2*a*b  x, both and 2*b*b  1 803 804 805 def _eval_numpy_(self, n, x): 791 806 """ 792 807 Evaluate ``self`` using numpy. 793 808 … … 798 813 sage: z2 = numpy.array([[1,2],[1,2]]) 799 814 sage: z3 = numpy.array([1,2,3.]) 800 815 sage: chebyshev_T(1,z) 801 array([ 1, 2])816 array([ 1., 2.]) 802 817 sage: chebyshev_T(1,z2) 803 array([[ 1, 2],804 [ 1, 2]])818 array([[ 1., 2.], 819 [ 1., 2.]]) 805 820 sage: chebyshev_T(1,z3) 806 821 array([ 1., 2., 3.]) 807 822 sage: chebyshev_T(z,0.1) 808 823 array([ 0.1 , 0.98]) 809 824 """ 810 825 from scipy.special import eval_chebyt 811 return eval_chebyt( args[0],args[1])826 return eval_chebyt(n, x) 812 827 813 def _derivative_(self, *args, **kwds):828 def _derivative_(self, n, x, diff_param): 814 829 """ 815 830 Return the derivative of :class:`chebyshev_T` in form of the Chebyshev 816 831 polynomial of the second kind :class:`chebyshev_U`. … … 828 843 ... 829 844 NotImplementedError: derivative w.r.t. to the index is not supported yet 830 845 """ 831 diff_param = kwds['diff_param'] 832 if diff_param == 0: 846 if diff_param == 0: 833 847 raise NotImplementedError("derivative w.r.t. to the index is not supported yet") 848 elif diff_param == 1: 849 return n*chebyshev_U(n1, x) 850 else: 851 raise ValueError("illegal differentiation parameter %s"%diff_param) 834 852 835 return args[0]*chebyshev_U(args[0]1,args[1])836 837 853 chebyshev_T = Func_chebyshev_T() 838 854 839 855 class Func_chebyshev_U(OrthogonalPolynomial): 840 856 """ 841 857 Class for the Chebyshev polynomial of the second kind. 842 858 843 859 REFERENCE: 844 860 845 861  [ASHandbook]_ 22.8.3 page 783 and 6.1.22 page 256. 846 862 847 863 EXAMPLES:: 848 849 sage: x = PolynomialRing(QQ, 'x').gen()850 sage: chebyshev_U(2, x)851 4* x^2  1852 sage: chebyshev_U(3, x)853 8* x^3  4*x864 865 sage: R.<t> = QQ[] 866 sage: chebyshev_U(2,t) 867 4*t^2  1 868 sage: chebyshev_U(3,t) 869 8*t^3  4*t 854 870 """ 855 871 def __init__(self): 856 872 """ 857 873 Init method for the chebyshev polynomials of the second kind. 858 874 859 875 EXAMPLES:: 860 876 861 877 sage: from sage.functions.orthogonal_polys import Func_chebyshev_U 862 878 sage: chebyshev_U2 = Func_chebyshev_U() 863 879 sage: chebyshev_U2(1,x) … … 866 882 OrthogonalPolynomial.__init__(self, "chebyshev_U", nargs=2, 867 883 conversions=dict(maxima='chebyshev_u', 868 884 mathematica='ChebyshevU')) 869 870 def _apply_formula_(self,*args):885 886 def eval_formula(self, n, x): 871 887 """ 872 Applies explicit formulas for :class:`chebyshev_U`. 873 This is much faster for numerical evaluation than maxima. 888 Evaluate ``chebyshev_U`` using an explicit formula. 874 889 See [ASHandbook]_ 227 (p. 782) for details on the recurions. 875 890 See also [EffCheby]_ for the recursion formulas. 876 891 892 INPUT: 893 894  ``n``  an integer 895 896  ``x``  a value to evaluate the polynomial at (this can be 897 any ring element) 898 877 899 EXAMPLES:: 878 900 879 sage: chebyshev_U._apply_formula_(2,0.1) == chebyshev_U._evalf_(2,0.1) 901 sage: chebyshev_U.eval_formula(10, x) 902 1024*x^10  2304*x^8 + 1792*x^6  560*x^4 + 60*x^2  1 903 sage: chebyshev_U.eval_formula(2, x) 904 1 905 sage: chebyshev_U.eval_formula(1, x) 906 0 907 sage: chebyshev_U.eval_formula(0, x) 908 1 909 sage: chebyshev_U.eval_formula(1, x) 910 2*x 911 sage: chebyshev_U.eval_formula(2,0.1) == chebyshev_U._evalf_(2,0.1) 880 912 True 881 913 """ 882 k = args[0] 883 x = args[1] 884 885 if k == 0: 886 return 1 887 if k == 1: 888 return 2*x 914 if n < 1: 915 return self.eval_formula(n2, x) 889 916 890 help1 = 1 891 help2 = 2*x 892 if is_Expression(x) and k <= 25: 893 # Recursion gives more compact representations for large k 894 help3 = 0 895 for j in xrange(0,floor(k/2)+1): 896 f = factorial(kj) / factorial(j) / factorial(k2*j) # Change to a binomial? 897 help3 = help3 + (1)**j * (2*x)**(k2*j) * f 898 return help3 899 917 res = parent(x).zero() 918 for j in xrange(0, n//2+1): 919 f = binomial(nj, j) 920 res += (1)**j * (2*x)**(n2*j) * f 921 return res 922 923 def eval_algebraic(self, n, x): 924 """ 925 Evaluate :class:`chebyshev_U` as polynomial, using a recursive 926 formula. 927 928 INPUT: 929 930  ``n``  an integer 931 932  ``x``  a value to evaluate the polynomial at (this can be 933 any ring element) 934 935 EXAMPLES:: 936 937 sage: chebyshev_U.eval_algebraic(5,x) 938 2*((2*x + 1)*(2*x  1)*x  4*(2*x^2  1)*x)*(2*x + 1)*(2*x  1) 939 sage: parent(chebyshev_U(3, Mod(8,9))) 940 Ring of integers modulo 9 941 sage: parent(chebyshev_U(3, Mod(1,9))) 942 Ring of integers modulo 9 943 sage: chebyshev_U(3,x) + chebyshev_U(1,x) 944 0 945 sage: chebyshev_U(1,Mod(5,8)) 946 0 947 sage: parent(chebyshev_U(1,Mod(5,8))) 948 Ring of integers modulo 8 949 sage: R.<t> = ZZ[] 950 sage: chebyshev_U.eval_algebraic(2, t) 951 1 952 sage: chebyshev_U.eval_algebraic(1, t) 953 0 954 sage: chebyshev_U.eval_algebraic(0, t) 955 1 956 sage: chebyshev_U.eval_algebraic(1, t) 957 2*t 958 sage: n = 97; x = RIF(pi/n) 959 sage: chebyshev_U(n1, cos(x)).contains_zero() 960 True 961 sage: R.<t> = Zp(2, 6, 'cappedabs')[] 962 sage: chebyshev_U(10^6+1, t) 963 (2 + O(2^6))*t + (O(2^6)) 964 """ 965 if n == 1: 966 return parent(x).zero() 967 if n >= 0: 968 return self._eval_recursive_(n, x)[0] 900 969 else: 901 return self._cheb_recur_(k,x)[0]970 return self._eval_recursive_(n2, x)[0] 902 971 903 904 def _cheb_recur_(self,n, x, both=False): 905 """ 906 Generalized recursion formula for Chebyshev polynomials. 907 Implementation suggested by Frederik Johansson. 908 returns (U(n,x), U(n1,x)), or (U(n,x), _) if both=False 909 910 EXAMPLES:: 911 912 sage: chebyshev_U._cheb_recur_(3,x) 913 (4*(2*x^2  1)*x, False) 914 sage: chebyshev_U._cheb_recur_(5,x)[0] 915 2*((2*x + 1)*(2*x  1)*x  4*(2*x^2  1)*x)*(2*x + 1)*(2*x  1) 916 sage: abs(pari('polchebyshev(5, 2, 0.1)')  chebyshev_U(5,0.1)) < 1e10 917 True 918 """ 919 920 if n == 0: 921 return 1, both and 2*x 922 if n == 1: 923 return 2*x, both and 4*x**21 924 925 a, b = self._cheb_recur_((n1)//2, x, True) 926 if n % 2 == 0: 927 return (b+a)*(ba), both and 2*b*(x*ba) 928 else: 929 return 2*a*(bx*a), both and (b+a)*(ba) 972 def _eval_recursive_(self, n, x, both=False): 973 """ 974 If ``both=True``, compute `(U(n,x), U(n1,x))` using a 975 recursive formula. 976 If ``both=False``, return instead a tuple `(U(n,x), False)`. 930 977 931 def _maxima_init_evaled_(self, *args): 978 EXAMPLES:: 979 980 sage: chebyshev_U._eval_recursive_(3, x) 981 (4*((2*x + 1)*(2*x  1)  2*x^2)*x, False) 982 sage: chebyshev_U._eval_recursive_(3, x, True) 983 (4*((2*x + 1)*(2*x  1)  2*x^2)*x, ((2*x + 1)*(2*x  1) + 2*x)*((2*x + 1)*(2*x  1)  2*x)) 984 985 """ 986 if n == 0: 987 return parent(x).one(), 2*x 988 989 assert n >= 1 990 a, b = self._eval_recursive_((n1)//2, x, True) 991 if n % 2 == 0: 992 return (b+a)*(ba), both and 2*b*(x*ba) 993 else: 994 return 2*a*(bx*a), both and (b+a)*(ba) 995 996 def _maxima_init_evaled_(self, n, x): 932 997 """ 933 998 Uses maxima to evaluate ``self``. 934 999 935 1000 EXAMPLES:: 936 1001 1002 sage: var('n, x') 1003 (n, x) 937 1004 sage: maxima(chebyshev_U(5,x)) 938 1005 32*x^532*x^3+6*x 939 sage: var('n')940 n941 1006 sage: maxima(chebyshev_U(n,x)) 942 1007 chebyshev_u(n,x) 943 1008 sage: maxima(chebyshev_U(2,x)) 944 1009 4*x^21 945 1010 """ 946 n = args[0] 947 x = args[1] 948 return maxima.eval('chebyshev_u({0},{1})'.format(n,x)) 1011 return maxima.eval('chebyshev_u({0},{1})'.format(n._maxima_init_(), x._maxima_init_())) 949 1012 950 def _evalf_(self, *args,**kwds):1013 def _evalf_(self, n, x, **kwds): 951 1014 """ 952 1015 Evaluate :class:`chebyshev_U` numerically with mpmath. 953 If index is an integer use recursive formula since it is faster,954 for chebyshev polynomials.955 1016 956 1017 EXAMPLES:: 957 1018 … … 962 1023 sage: chebyshev_U._evalf_(1.5, Mod(8,9)) 963 1024 Traceback (most recent call last): 964 1025 ... 965 ValueError: No compatible type!1026 TypeError: cannot evaluate chebyshev_U with parent Ring of integers modulo 9 966 1027 """ 967 if args[0] in ZZ and args[0] >= 0:968 return self._cheb_recur_(*args)[0]969 1028 try: 970 1029 real_parent = kwds['parent'] 971 1030 except KeyError: 972 real_parent = parent( args[1])1031 real_parent = parent(x) 973 1032 974 x_set = False 975 if hasattr(real_parent,"precision"): # Check if we have a data type with precision 976 x = args[1] 977 step_parent = real_parent 978 x_set = True 979 else: 980 if args[1] in RR: 981 x = RR(args[1]) 982 step_parent = RR 983 x_set = True 984 elif args[1] in CC: 985 x = CC(args[1]) 986 step_parent = CC 987 x_set = True 1033 if not is_RealField(real_parent) and not is_ComplexField(real_parent): 1034 # parent is not a real or complex field: figure out a good parent 1035 if x in RR: 1036 x = RR(x) 1037 real_parent = RR 1038 elif x in CC: 1039 x = CC(x) 1040 real_parent = CC 988 1041 989 if not x_set:990 raise ValueError("No compatible type!")1042 if not is_RealField(real_parent) and not is_ComplexField(real_parent): 1043 raise TypeError("cannot evaluate chebyshev_U with parent %s"%real_parent) 991 1044 992 1045 from sage.libs.mpmath.all import call as mpcall 993 1046 from sage.libs.mpmath.all import chebyu as mpchebyu 994 1047 995 return mpcall(mpchebyu, args[0],args[1],parent = step_parent)1048 return mpcall(mpchebyu, n, x, parent=real_parent) 996 1049 997 def _eval_special_values_(self, *args):1050 def _eval_special_values_(self, n, x): 998 1051 """ 999 Special values that known. [ASHandbook]_ 22.4 (p.777). 1052 Values known for special values of x. 1053 See [ASHandbook]_ 22.4 (p.777). 1000 1054 1001 1055 EXAMPLES:: 1002 1056 1003 1057 sage: var('n') 1004 1058 n 1005 1059 sage: chebyshev_U(n,1) 1006 1060 n + 1 1061 sage: chebyshev_U(n,0) 1062 1/2*(1)^(1/2*n)*((1)^n + 1) 1007 1063 sage: chebyshev_U(n,1) 1008 1064 (1)^n*(n + 1) 1009 sage: chebyshev_U._eval_special_values_( 26, Mod(0,9))1065 sage: chebyshev_U._eval_special_values_(n, 2) 1010 1066 Traceback (most recent call last): 1011 1067 ... 1012 ValueError: Value not found! 1013 sage: parent(chebyshev_U(3, Mod(8,9))) 1014 Ring of integers modulo 9 1015 sage: parent(chebyshev_U(3, Mod(1,9))) 1016 Ring of integers modulo 9 1017 sage: chebyshev_U(n, 0) 1018 1/2*(1)^(1/2*n)*((1)^n + 1) 1019 sage: chebyshev_U(3,x) + chebyshev_U(1,x) 1020 0 1021 sage: chebyshev_U._eval_special_values_(1.5, Mod(8,9)) 1022 Traceback (most recent call last): 1023 ... 1024 ValueError: No special values for non integral indices! 1025 sage: chebyshev_U(1,Mod(5,8)) 1026 0 1027 sage: parent(chebyshev_U(1,Mod(5,8))) 1028 Ring of integers modulo 8 1068 ValueError: no special value found 1029 1069 """ 1030 if (not is_Expression(args[0])) and (not args[0] in ZZ):1031 r aise ValueError("No special values for non integral indices!")1070 if x == 1: 1071 return x*(n+1) 1032 1072 1033 if args[0] == 1: 1034 return args[1]*0 1035 1036 if args[1] == 1: 1037 return args[1]*(args[0]+1) 1038 1039 if args[1] == 1: 1040 return args[1]**args[0]*(args[0]+1) 1073 if x == 1: 1074 return x**n*(n+1) 1041 1075 1042 if (args[1] == 0 and args[1] in CC):1043 return (1+(1)** args[0])*(1)**(args[0]/2)/21076 if x == 0: 1077 return (1+(1)**n)*(1)**(n/2)/2 1044 1078 1045 if args[0] < 0 and args[0] in ZZ: 1046 return self._eval_(args[0]2,args[1]) 1079 raise ValueError("no special value found") 1047 1080 1048 raise ValueError("Value not found!") 1049 1050 def _eval_numpy_(self, *args): 1081 def _eval_numpy_(self, n, x): 1051 1082 """ 1052 1083 Evaluate ``self`` using numpy. 1053 1084 … … 1058 1089 sage: z2 = numpy.array([[1,2],[1,2]]) 1059 1090 sage: z3 = numpy.array([1,2,3.]) 1060 1091 sage: chebyshev_U(1,z) 1061 array([ 2, 4])1092 array([ 2., 4.]) 1062 1093 sage: chebyshev_U(1,z2) 1063 array([[ 2, 4],1064 [ 2, 4]])1094 array([[ 2., 4.], 1095 [ 2., 4.]]) 1065 1096 sage: chebyshev_U(1,z3) 1066 1097 array([ 2., 4., 6.]) 1067 1098 sage: chebyshev_U(z,0.1) 1068 1099 array([ 0.2 , 0.96]) 1069 1100 """ 1070 1101 from scipy.special import eval_chebyu 1071 return eval_chebyu( args[0],args[1])1102 return eval_chebyu(n, x) 1072 1103 1073 1104 1074 def _derivative_(self, *args, **kwds):1105 def _derivative_(self, n, x, diff_param): 1075 1106 """ 1076 1107 Return the derivative of :class:`chebyshev_U` in form of the Chebyshev 1077 1108 polynomials of the first and second kind. 1078 1109 1079 1110 EXAMPLES:: 1080 1111 1081 1112 sage: var('k') … … 1089 1120 ... 1090 1121 NotImplementedError: derivative w.r.t. to the index is not supported yet 1091 1122 """ 1092 diff_param = kwds['diff_param'] 1093 if diff_param == 0: 1123 if diff_param == 0: 1094 1124 raise NotImplementedError("derivative w.r.t. to the index is not supported yet") 1095 1096 return ((args[0]+1)*chebyshev_T(args[0]+1,args[1])args[1]* 1097 chebyshev_U(args[0],args[1]))/(args[1]**21) 1125 elif diff_param == 1: 1126 return ((n+1)*chebyshev_T(n+1, x)  x*chebyshev_U(n,x))/(x*x1) 1127 else: 1128 raise ValueError("illegal differentiation parameter %s"%diff_param) 1098 1129 1099 1130 chebyshev_U = Func_chebyshev_U() 1100 1131 … … 1122 1153 sage: gen_laguerre(3,0,x) 1123 1154 1/6*x^3 + 3/2*x^2  3*x + 1 1124 1155 """ 1125 from sage.functions.all import sqrt1126 1156 _init() 1127 1157 return sage_eval(maxima.eval('gen_laguerre(%s,%s,x)'%(ZZ(n),a)), locals={'x':x}) 1128 1158