Ticket #9706: 9706_review.patch

File 9706_review.patch, 46.9 KB (added by jdemeyer, 8 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  
    1212
    1313-  The Chebyshev polynomial of the first kind arises as a solution
    1414   to the differential equation
    15    
     15
    1616   .. math::
    1717
    18          (1-x^2)\,y'' - x\,y' + n^2\,y = 0 
     18         (1-x^2)\,y'' - x\,y' + n^2\,y = 0
    1919
    2020
    2121   and those of the second kind as a solution to
    22    
     22
    2323   .. math::
    2424
    25          (1-x^2)\,y'' - 3x\,y' + n(n+2)\,y = 0. 
     25         (1-x^2)\,y'' - 3x\,y' + n(n+2)\,y = 0.
    2626
    2727
    2828   The Chebyshev polynomials of the first kind are defined by the
    2929   recurrence relation
    30    
     30
    3131   .. math::
    3232
    33      T_0(x) = 1 \, T_1(x) = x \, T_{n+1}(x) = 2xT_n(x) - T_{n-1}(x). \, 
     33     T_0(x) = 1 \, T_1(x) = x \, T_{n+1}(x) = 2xT_n(x) - T_{n-1}(x). \,
    3434
    3535
    3636   The Chebyshev polynomials of the second kind are defined by the
    3737   recurrence relation
    38    
     38
    3939   .. math::
    4040
    41      U_0(x) = 1 \, U_1(x) = 2x \, U_{n+1}(x) = 2xU_n(x) - U_{n-1}(x). \, 
     41     U_0(x) = 1 \, U_1(x) = 2x \, U_{n+1}(x) = 2xU_n(x) - U_{n-1}(x). \,
    4242
    4343
    4444
    4545   For integers `m,n`, they satisfy the orthogonality
    4646   relations
    47    
     47
    4848   .. math::
    4949
    50      \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. 
     50     \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.
    5151
    5252
    5353   and
    5454
    55    
     55
    5656   .. math::
    5757
    58      \int_{-1}^1 U_n(x)U_m(x)\sqrt{1-x^2}\,dx =\frac{\pi}{2}\delta_{m,n}. 
     58     \int_{-1}^1 U_n(x)U_m(x)\sqrt{1-x^2}\,dx =\frac{\pi}{2}\delta_{m,n}.
    5959
    6060
    6161
     
    6363   transliterations: Tchebyshef or Tschebyscheff).
    6464
    6565-  The Hermite polynomials are defined either by
    66    
     66
    6767   .. math::
    6868
    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}
    7070
    7171
    7272   (the "probabilists' Hermite polynomials"), or by
    7373
    74    
     74
    7575   .. math::
    7676
    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}
    7878
    7979
    8080   (the "physicists' Hermite polynomials"). Sage (via Maxima)
    8181   implements the latter flavor. These satisfy the orthogonality
    8282   relation
    83    
     83
    8484   .. math::
    8585
    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}
    8787
    8888
    8989
     
    105105   and satisfy the orthogonality relation
    106106
    107107   .. math::
    108    
     108
    109109      \int_{-1}^{1} P_m(x) P_n(x)\,dx = {\frac{2}{2n + 1}} \delta_{mn}
    110110
    111111   The *Legendre function of the second kind* `Q_n(x)` is another
     
    118118
    119119   .. math::
    120120
    121      \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} 
     121     \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}
    122122
    123123
    124124   Assuming `0 \le m \le \ell`, they satisfy the orthogonality
    125125   relation:
    126    
     126
    127127   .. math::
    128128
    129       \int_{-1}^{1} P_k ^{(m)} P_\ell ^{(m)} dx  = \frac{2 (\ell+m)!}{(2\ell+1)(\ell-m)!}\ \delta _{k,\ell}, 
     129      \int_{-1}^{1} P_k ^{(m)} P_\ell ^{(m)} dx  = \frac{2 (\ell+m)!}{(2\ell+1)(\ell-m)!}\ \delta _{k,\ell},
    130130
    131131
    132132   where `\delta _{k,\ell}` is the Kronecker delta.
     
    135135   `Q_\ell^m(x)` can be given in terms of the "usual"
    136136   Legendre polynomials by
    137137
    138    
     138
    139139   .. math::
    140140
    141      Q_\ell^m(x)   =  (-1)^m(1-x^2)^{m/2}\frac{d^m}{dx^m}Q_\ell(x). 
     141     Q_\ell^m(x)   =  (-1)^m(1-x^2)^{m/2}\frac{d^m}{dx^m}Q_\ell(x).
    142142
    143143
    144144
    145145   They are named after Adrien-Marie Legendre.
    146146
    147147-  Laguerre polynomials may be defined by the Rodrigues formula
    148    
     148
    149149   .. math::
    150150
    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).
    152152
    153153
    154154   They are solutions of Laguerre's equation:
    155155
    156    
     156
    157157   .. math::
    158158
    159       x\,y'' + (1 - x)\,y' + n\,y = 0\, 
     159      x\,y'' + (1 - x)\,y' + n\,y = 0\,
    160160
    161161   and satisfy the orthogonality relation
    162162
    163    
     163
    164164   .. math::
    165165
    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}.
    167167
    168168
    169169
    170170   The generalized Laguerre polynomials may be defined by the
    171171   Rodrigues formula:
    172172
    173    
     173
    174174   .. math::
    175175
    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) .
    177177
    178178
    179179   (These are also sometimes called the associated Laguerre
     
    185185-  Jacobi polynomials are a class of orthogonal polynomials. They
    186186   are obtained from hypergeometric series in cases where the series
    187187   is in fact finite:
    188    
     188
    189189   .. math::
    190190
    191      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) , 
     191     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) ,
    192192
    193193
    194194   where `()_n` is Pochhammer's symbol (for the rising
    195195   factorial), (Abramowitz and Stegun p561.) and thus have the
    196196   explicit expression
    197197
    198    
     198
    199199   .. math::
    200200
    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{z-1}{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{z-1}{2}\right)^m .
    202202
    203203
    204204
     
    208208   the Jacobi polynomials `P_n^{(\alpha,\beta)}(x)` with
    209209   `\alpha=\beta=a-1/2` by
    210210
    211    
     211
    212212   .. math::
    213213
    214      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). 
     214     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).
    215215
    216216
    217217   They satisfy the orthogonality relation
    218    
     218
    219219   .. math::
    220220
    221      \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)} , 
     221     \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)} ,
    222222
    223223
    224224   for `a>-1/2`. They are obtained from hypergeometric series
    225225   in cases where the series is in fact finite:
    226226
    227    
     227
    228228   .. math::
    229229
    230      C_n^{(a)}(z) =\frac{(2a)^{\underline{n}}}{n!} \,_2F_1\left(-n,2a+n;a+\frac{1}{2};\frac{1-z}{2}\right) 
     230     C_n^{(a)}(z) =\frac{(2a)^{\underline{n}}}{n!} \,_2F_1\left(-n,2a+n;a+\frac{1}{2};\frac{1-z}{2}\right)
    231231
    232232
    233233   where `\underline{n}` is the falling factorial. (See
     
    242242
    243243.. math::
    244244
    245          (x)_n=x(x+1)(x+2)\cdots(x+n-1)=\frac{(x+n-1)!}{(x-1)!}. 
     245         (x)_n=x(x+1)(x+2)\cdots(x+n-1)=\frac{(x+n-1)!}{(x-1)!}.
    246246
    247247
    248248On the other hand, the "falling factorial" or "lower factorial" is
    249249
    250250.. math::
    251251
    252      x^{\underline{n}}=\frac{x!}{(x-n)!} , 
     252     x^{\underline{n}}=\frac{x!}{(x-n)!} ,
    253253
    254254
    255255in the notation of Ronald L. Graham, Donald E. Knuth and Oren
     
    287287
    288288.. :wikipedia:`Associated_Legendre_polynomials`
    289289
    290 .. [EffCheby] Wolfram Koepf: Effcient Computation of Chebyshev Polynomials 
     290.. [EffCheby] Wolfram Koepf: Effcient Computation of Chebyshev Polynomials
    291291    in Computer Algebra
    292     Computer Algebra Systems: A Practical Guide. 
     292    Computer Algebra Systems: A Practical Guide.
    293293    John Wiley, Chichester (1999): 79-99.
    294294
    295295AUTHORS:
     
    318318import warnings
    319319
    320320from sage.misc.sage_eval import sage_eval
    321 from sage.rings.all import ZZ, RR, CC, RIF, CIF
     321from sage.rings.all import ZZ, RR, CC
     322from sage.rings.real_mpfr import is_RealField
     323from sage.rings.complex_field import is_ComplexField
    322324from sage.calculus.calculus import maxima
    323325
    324326
    325 from sage.symbolic.function import BuiltinFunction, GinacFunction, is_inexact
     327from sage.symbolic.function import BuiltinFunction
    326328from 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
     329from sage.functions.other import factorial, binomial
    335330from sage.structure.coerce import parent
    336331
    337332_done = False
     
    340335    Internal function which checks if Maxima has loaded the
    341336    "orthopoly" package.  All functions using this in this
    342337    file should call this function first.
    343    
     338
    344339    TEST:
    345340
    346     The global starts ``False``:: 
     341    The global starts ``False``::
    347342
    348343        sage: sage.functions.orthogonal_polys._done
    349344        False
    350345
    351346    Then after using one of these functions, it changes::
    352347
    353         sage: from sage.functions.orthogonal_polys import legendre_P 
     348        sage: from sage.functions.orthogonal_polys import legendre_P
    354349        sage: legendre_P(2,x)
    355350        3/2*(x - 1)^2 + 3*x - 2
    356351        sage: sage.functions.orthogonal_polys._done
     
    372367    _done = True
    373368
    374369
    375 
    376370class OrthogonalPolynomial(BuiltinFunction):
    377371    """
    378372    Base class for orthogonal polynomials.
     
    381375    they share similar properties. The evaluation as a polynomial
    382376    is either done via maxima, or with pynac.
    383377
    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.
    386381    """
    387382    def __init__(self, name, nargs=2, latex_name=None, conversions={}):
    388383        """
    389384        :class:`OrthogonalPolynomial` class needs the same input parameter as
    390385        it's parent class.
    391        
     386
    392387        EXAMPLES::
    393388
    394389            sage: from sage.functions.orthogonal_polys import OrthogonalPolynomial
     
    397392            testo_P
    398393        """
    399394        try:
    400             self._maxima_name = conversions['maxima'] 
     395            self._maxima_name = conversions['maxima']
    401396        except KeyError:
    402397            self._maxima_name = None
    403398
    404         super(OrthogonalPolynomial,self).__init__(name=name, nargs=nargs, 
     399        super(OrthogonalPolynomial,self).__init__(name=name, nargs=nargs,
    405400                                 latex_name=latex_name, conversions=conversions)
    406401
    407402    def _maxima_init_evaled_(self, *args):
    408403        r"""
    409404        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.
    430406
    431407        EXAMPLES::
    432408
    433409            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)
    436425            Traceback (most recent call last):
    437426            ...
    438             NotImplementedError: no recursive calculation of values implemented
     427            NotImplementedError: no explicit calculation of values implemented
    439428        """
    440         raise NotImplementedError("no recursive calculation of values implemented")
     429        raise NotImplementedError("no explicit calculation of values implemented")
    441430
    442     def _eval_special_values_(self,*args):
     431    def _eval_special_values_(self, *args):
    443432        """
    444433        Evaluate the polynomial explicitly for special values.
    445434
    446435        EXAMPLES::
    447            
     436
    448437            sage: var('n')
    449438            n
    450439            sage: chebyshev_T(n,-1)
     
    452441        """
    453442        raise ValueError("no special values known")
    454443
    455     def _eval_(self, *args):
     444    def _eval_(self, n, *args):
    456445        """
    457446        The _eval_ method decides which evaluation suits best
    458447        for the given input, and returns a proper value.
    459        
     448
    460449        EXAMPLES::
    461450
     451            sage: var('n,x')
     452            (n, x)
    462453            sage: chebyshev_T(5,x)
    463             16*x^5 - 20*x^3 + 5*x 
    464             sage: var('n')
    465             n
     454            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
    466457            sage: chebyshev_T(n,-1)
    467458            (-1)^n
    468459            sage: chebyshev_T(-7,x)
    469460            64*x^7 - 112*x^5 + 56*x^3 - 7*x
    470461            sage: chebyshev_T(3/2,x)
    471462            chebyshev_T(3/2, x)
    472             sage: x = PolynomialRing(QQ, 'x').gen()
    473             sage: chebyshev_T(2,x)
    474             2*x^2 - 1
    475             sage: chebyshev_U(2,x)
    476             4*x^2 - 1
     463            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
    477468            sage: parent(chebyshev_T(4, RIF(5)))
    478469            Real Interval Field with 53 bits of precision
    479470            sage: RR2 = RealField(5)
     
    482473            sage: chebyshev_T(5,Qp(3)(2))
    483474            2 + 3^2 + 3^3 + 3^4 + 3^5 + O(3^20)
    484475            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
    486477            chebyshev_T(100001/2, 2)
    487478            sage: chebyshev_U._eval_(1.5, Mod(8,9)) is None
    488479            True
    489480        """
    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)
    503482
    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)
    508491
    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
    519498                return None
    520499
    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
    528511            return None
    529512
    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        """
    537515        This overides the call method from SageObject to avoid problems with coercions,
    538516        since the _eval_ method is able to handle more data types than symbolic functions
    539517        would normally allow.
    540518        Thus we have the distinction between algebraic objects (if n is an integer),
    541519        and else as symbolic function.
    542        
    543         EXAMPLES:: 
    544        
    545             sage: K.<a> = NumberField(x^3-x-1) 
    546             sage: chebyshev_T(5, a) 
    547             16*a^2 + a - 4 
     520
     521        EXAMPLES::
     522
     523            sage: K.<a> = NumberField(x^3-x-1)
     524            sage: chebyshev_T(5, a)
     525            16*a^2 + a - 4
    548526            sage: chebyshev_T(5,MatrixSpace(ZZ, 2)([1, 2, -4, 7]))
    549527            [-40799  44162]
    550528            [-88324  91687]
    551529            sage: R.<x> = QQ[]
    552530            sage: parent(chebyshev_T(5, x))
    553531            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
    564545
    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)
    570547
    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):
     548class Func_chebyshev_T(OrthogonalPolynomial):
    579549    """
    580550    Chebyshev polynomials of the first kind.
    581551
     
    609579        super(Func_chebyshev_T,self).__init__("chebyshev_T", nargs=2,
    610580                                      conversions=dict(maxima='chebyshev_t',
    611581                                                       mathematica='ChebyshevT'))
    612    
    613     def _eval_special_values_(self,*args):
     582
     583    def _eval_special_values_(self, n, x):
    614584        """
    615585        Values known for special values of x.
    616586        For details see [ASHandbook]_ 22.4 (p. 777)
     
    621591            n
    622592            sage: chebyshev_T(n,1)
    623593            1
     594            sage: chebyshev_T(n,0)
     595            1/2*(-1)^(1/2*n)*((-1)^n + 1)
    624596            sage: chebyshev_T(n,-1)
    625597            (-1)^n
    626             sage: chebyshev_T(-7, x) - chebyshev_T(7,x)
    627             0
    628598            sage: chebyshev_T._eval_special_values_(3/2,x)
    629599            Traceback (most recent call last):
    630600            ...
    631             ValueError: No special values for non integral indices!
     601            ValueError: no special value found
    632602            sage: chebyshev_T._eval_special_values_(n, 0.1)
    633603            Traceback (most recent call last):
    634604            ...
    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
    640606        """
    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
    649609
    650         if (args[-1] == 0 and args[-1] in CC):
    651             return (1+(-1)**args[0])*(-1)**(args[0]/2)/2
     610        if x == -1:
     611            return x**n
    652612
    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
    655615
    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):
    659619        """
    660620        Evaluates :class:`chebyshev_T` numerically with mpmath.
    661         If the index is an integer we use the recursive formula since
    662         it is faster.
    663621
    664622        EXAMPLES::
    665623
    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))
    667627            2.261953700000000000000e7
    668             sage: chebyshev_T(10,I).n()
     628            sage: chebyshev_T._evalf_(10,I)
    669629            -3363.00000000000
    670             sage: chebyshev_T(5,0.3).n()
     630            sage: chebyshev_T._evalf_(5,0.3)
    671631            0.998880000000000
    672632            sage: chebyshev_T(1/2, 0)
    673633            0.707106781186548
     634            sage: chebyshev_T(1/2, 3/2)
     635            1.11803398874989
    674636            sage: chebyshev_T._evalf_(1.5, Mod(8,9))
    675637            Traceback (most recent call last):
    676638            ...
    677             ValueError: No compatible type!
     639            TypeError: cannot evaluate chebyshev_T with parent Ring of integers modulo 9
    678640
     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
    679657        """
    680         if args[0] in ZZ and args[0] >= 0:
    681             return self._cheb_recur_(*args)[0]
    682        
    683658        try:
    684659            real_parent = kwds['parent']
    685660        except KeyError:
    686             real_parent = parent(args[-1])
     661            real_parent = parent(x)
    687662
    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
    702671
    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
    706675        from sage.libs.mpmath.all import call as mpcall
    707676        from sage.libs.mpmath.all import chebyt as mpchebyt
    708677
    709         return mpcall(mpchebyt,args[0],x,parent=step_parent)
     678        return mpcall(mpchebyt, n, x, parent=real_parent)
    710679
    711     def _maxima_init_evaled_(self, *args):
     680    def _maxima_init_evaled_(self, n, x):
    712681        """
    713682        Evaluate the Chebyshev polynomial ``self`` with maxima.
    714683
    715684        EXAMPLES::
    716685
     686            sage: var('n, x')
     687            (n, x)
    717688            sage: chebyshev_T._maxima_init_evaled_(1,x)
    718689            '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))
    723692
    724693        """
    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_()))
    728695
    729        
    730     def _apply_formula_(self,*args):
     696
     697    def eval_formula(self, n, x):
    731698        """
    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.
    734700        See [ASHandbook]_ 227 (p. 782) for details for the recurions.
    735701        See also [EffCheby]_ for fast evaluation techniques.
    736702
     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
    737710        EXAMPLES::
    738711
    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)
    740719            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) - x
    743             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()
    744723            512*x^10 - 1280*x^8 + 1120*x^6 - 400*x^4 + 50*x^2 - 1
    745724
    746725        """
    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()
    749730
    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(k-j-1) / factorial(j) / factorial(k-2*j)
    762                 help3 = help3 + (-1)**j * (2*x)**(k-2*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(n-1-j) / factorial(j) / factorial(n-2*j)
     734            res += (-1)**j * (2*x)**(n-2*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, 'capped-abs')[]
     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]
    765778        else:
    766             return self._cheb_recur_(k,x)[0]
     779            return self._eval_recursive_(-n, x)[0]
    767780
    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(n-1,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(n-1,x))` using a
     784        recursive formula.
     785        If ``both=False``, return instead a tuple `(T(n,x), False)`.
    788786
     787        EXAMPLES::
    789788
    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):
    791806        """
    792807        Evaluate ``self`` using numpy.
    793808
     
    798813            sage: z2 = numpy.array([[1,2],[1,2]])
    799814            sage: z3 = numpy.array([1,2,3.])
    800815            sage: chebyshev_T(1,z)
    801             array([1, 2])
     816            array([ 1.,  2.])
    802817            sage: chebyshev_T(1,z2)
    803             array([[1, 2],
    804                    [1, 2]])
     818            array([[ 1.,  2.],
     819                   [ 1.,  2.]])
    805820            sage: chebyshev_T(1,z3)
    806821            array([ 1.,  2.,  3.])
    807822            sage: chebyshev_T(z,0.1)
    808823            array([ 0.1 , -0.98])
    809824        """
    810825        from scipy.special import eval_chebyt
    811         return eval_chebyt(args[0],args[-1])
     826        return eval_chebyt(n, x)
    812827
    813     def _derivative_(self, *args, **kwds):
     828    def _derivative_(self, n, x, diff_param):
    814829        """
    815830        Return the derivative of :class:`chebyshev_T` in form of the Chebyshev
    816831        polynomial of the second kind :class:`chebyshev_U`.
     
    828843            ...
    829844            NotImplementedError: derivative w.r.t. to the index is not supported yet
    830845        """
    831         diff_param = kwds['diff_param']
    832         if diff_param == 0:
     846        if diff_param == 0:
    833847            raise NotImplementedError("derivative w.r.t. to the index is not supported yet")
     848        elif diff_param == 1:
     849            return n*chebyshev_U(n-1, x)
     850        else:
     851            raise ValueError("illegal differentiation parameter %s"%diff_param)
    834852
    835         return args[0]*chebyshev_U(args[0]-1,args[1])
    836    
    837853chebyshev_T = Func_chebyshev_T()
    838854
    839855class Func_chebyshev_U(OrthogonalPolynomial):
    840856    """
    841857    Class for the Chebyshev polynomial of the second kind.
    842    
     858
    843859    REFERENCE:
    844860
    845861    - [ASHandbook]_ 22.8.3 page 783 and 6.1.22 page 256.
    846    
     862
    847863    EXAMPLES::
    848    
    849         sage: x = PolynomialRing(QQ, 'x').gen()
    850         sage: chebyshev_U(2,x)
    851         4*x^2 - 1
    852         sage: chebyshev_U(3,x)
    853         8*x^3 - 4*x
     864
     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
    854870    """
    855871    def __init__(self):
    856872        """
    857873        Init method for the chebyshev polynomials of the second kind.
    858874
    859875        EXAMPLES::
    860        
     876
    861877            sage: from sage.functions.orthogonal_polys import Func_chebyshev_U
    862878            sage: chebyshev_U2 = Func_chebyshev_U()
    863879            sage: chebyshev_U2(1,x)
     
    866882        OrthogonalPolynomial.__init__(self, "chebyshev_U", nargs=2,
    867883                                      conversions=dict(maxima='chebyshev_u',
    868884                                                       mathematica='ChebyshevU'))
    869        
    870     def _apply_formula_(self,*args):
     885
     886    def eval_formula(self, n, x):
    871887        """
    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.
    874889        See [ASHandbook]_ 227 (p. 782) for details on the recurions.
    875890        See also [EffCheby]_ for the recursion formulas.
    876891
     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
    877899        EXAMPLES::
    878900
    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)
    880912            True
    881913        """
    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(-n-2, x)
    889916
    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(k-j) / factorial(j) / factorial(k-2*j) # Change to a binomial?
    897                 help3 = help3 + (-1)**j * (2*x)**(k-2*j) * f
    898             return help3
    899            
     917        res = parent(x).zero()
     918        for j in xrange(0, n//2+1):
     919            f = binomial(n-j, j)
     920            res += (-1)**j * (2*x)**(n-2*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(n-1, cos(x)).contains_zero()
     960            True
     961            sage: R.<t> = Zp(2, 6, 'capped-abs')[]
     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]
    900969        else:
    901             return self._cheb_recur_(k,x)[0]
     970            return -self._eval_recursive_(-n-2, x)[0]
    902971
    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(n-1,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)) < 1e-10
    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**2-1
    924            
    925         a, b = self._cheb_recur_((n-1)//2, x, True)
    926         if n % 2 == 0:
    927             return (b+a)*(b-a), both and 2*b*(x*b-a)
    928         else:
    929             return 2*a*(b-x*a), both and (b+a)*(b-a)
     972    def _eval_recursive_(self, n, x, both=False):
     973        """
     974        If ``both=True``, compute `(U(n,x), U(n-1,x))` using a
     975        recursive formula.
     976        If ``both=False``, return instead a tuple `(U(n,x), False)`.
    930977
    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_((n-1)//2, x, True)
     991        if n % 2 == 0:
     992            return (b+a)*(b-a), both and 2*b*(x*b-a)
     993        else:
     994            return 2*a*(b-x*a), both and (b+a)*(b-a)
     995
     996    def _maxima_init_evaled_(self, n, x):
    932997        """
    933998        Uses maxima to evaluate ``self``.
    934999
    9351000        EXAMPLES::
    9361001
     1002            sage: var('n, x')
     1003            (n, x)
    9371004            sage: maxima(chebyshev_U(5,x))
    9381005            32*x^5-32*x^3+6*x
    939             sage: var('n')
    940             n
    9411006            sage: maxima(chebyshev_U(n,x))
    9421007            chebyshev_u(n,x)
    9431008            sage: maxima(chebyshev_U(2,x))
    9441009            4*x^2-1
    9451010        """
    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_()))
    9491012
    950     def _evalf_(self, *args,**kwds):
     1013    def _evalf_(self, n, x, **kwds):
    9511014        """
    9521015        Evaluate :class:`chebyshev_U` numerically with mpmath.
    953         If index is an integer use recursive formula since it is faster,
    954         for chebyshev polynomials.
    9551016
    9561017        EXAMPLES::
    9571018
     
    9621023            sage: chebyshev_U._evalf_(1.5, Mod(8,9))
    9631024            Traceback (most recent call last):
    9641025            ...
    965             ValueError: No compatible type!
     1026            TypeError: cannot evaluate chebyshev_U with parent Ring of integers modulo 9
    9661027        """
    967         if args[0] in ZZ and args[0] >= 0:
    968             return self._cheb_recur_(*args)[0]
    9691028        try:
    9701029            real_parent = kwds['parent']
    9711030        except KeyError:
    972             real_parent = parent(args[-1])
     1031            real_parent = parent(x)
    9731032
    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
    9881041
    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)
    9911044
    9921045        from sage.libs.mpmath.all import call as mpcall
    9931046        from sage.libs.mpmath.all import chebyu as mpchebyu
    9941047
    995         return mpcall(mpchebyu,args[0],args[-1],parent = step_parent)
     1048        return mpcall(mpchebyu, n, x, parent=real_parent)
    9961049
    997     def _eval_special_values_(self,*args):
     1050    def _eval_special_values_(self, n, x):
    9981051        """
    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).
    10001054
    10011055        EXAMPLES::
    1002        
     1056
    10031057            sage: var('n')
    10041058            n
    10051059            sage: chebyshev_U(n,1)
    10061060            n + 1
     1061            sage: chebyshev_U(n,0)
     1062            1/2*(-1)^(1/2*n)*((-1)^n + 1)
    10071063            sage: chebyshev_U(n,-1)
    10081064            (-1)^n*(n + 1)
    1009             sage: chebyshev_U._eval_special_values_(26, Mod(0,9))
     1065            sage: chebyshev_U._eval_special_values_(n, 2)
    10101066            Traceback (most recent call last):
    10111067            ...
    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
    10291069        """
    1030         if (not is_Expression(args[0])) and (not args[0] in ZZ):
    1031             raise ValueError("No special values for non integral indices!")
     1070        if x == 1:
     1071            return x*(n+1)
    10321072
    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)
    10411075
    1042         if (args[-1] == 0 and args[-1] in CC):
    1043             return (1+(-1)**args[0])*(-1)**(args[0]/2)/2
     1076        if x == 0:
     1077            return (1+(-1)**n)*(-1)**(n/2)/2
    10441078
    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")
    10471080
    1048         raise ValueError("Value not found!")
    1049 
    1050     def _eval_numpy_(self, *args):
     1081    def _eval_numpy_(self, n, x):
    10511082        """
    10521083        Evaluate ``self`` using numpy.
    10531084
     
    10581089            sage: z2 = numpy.array([[1,2],[1,2]])
    10591090            sage: z3 = numpy.array([1,2,3.])
    10601091            sage: chebyshev_U(1,z)
    1061             array([2, 4])
     1092            array([ 2.,  4.])
    10621093            sage: chebyshev_U(1,z2)
    1063             array([[2, 4],
    1064                    [2, 4]])
     1094            array([[ 2.,  4.],
     1095                   [ 2.,  4.]])
    10651096            sage: chebyshev_U(1,z3)
    10661097            array([ 2.,  4.,  6.])
    10671098            sage: chebyshev_U(z,0.1)
    10681099            array([ 0.2 , -0.96])
    10691100        """
    10701101        from scipy.special import eval_chebyu
    1071         return eval_chebyu(args[0],args[1])
     1102        return eval_chebyu(n, x)
    10721103
    10731104
    1074     def _derivative_(self, *args, **kwds):
     1105    def _derivative_(self, n, x, diff_param):
    10751106        """
    10761107        Return the derivative of :class:`chebyshev_U` in form of the Chebyshev
    10771108        polynomials of the first and second kind.
    1078        
     1109
    10791110        EXAMPLES::
    10801111
    10811112            sage: var('k')
     
    10891120            ...
    10901121            NotImplementedError: derivative w.r.t. to the index is not supported yet
    10911122        """
    1092         diff_param = kwds['diff_param']
    1093         if diff_param == 0:
     1123        if diff_param == 0:
    10941124            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]**2-1)
     1125        elif diff_param == 1:
     1126            return ((n+1)*chebyshev_T(n+1, x) - x*chebyshev_U(n,x))/(x*x-1)
     1127        else:
     1128            raise ValueError("illegal differentiation parameter %s"%diff_param)
    10981129
    10991130chebyshev_U = Func_chebyshev_U()
    11001131
     
    11221153        sage: gen_laguerre(3,0,x)
    11231154        -1/6*x^3 + 3/2*x^2 - 3*x + 1
    11241155    """
    1125     from sage.functions.all import sqrt
    11261156    _init()
    11271157    return sage_eval(maxima.eval('gen_laguerre(%s,%s,x)'%(ZZ(n),a)), locals={'x':x})
    11281158