Ticket #9706: trac_9706_chebyshev_collection.patch

File trac_9706_chebyshev_collection.patch, 30.1 KB (added by maldun, 7 years ago)

collection of all patches

  • sage/functions/orthogonal_polys.py

    # HG changeset patch
    # User Stefan Reiterer <domors@gmx.net>
    # Date 1386543882 -3600
    # Node ID f66a8803fcf7324f1a4475614bf453709fb8b9d7
    # Parent  f0ee3538887fe739601babb54e177ec5e1133b7a
    trac 9706: Collective patch. Bugfixes, extensions, optimizations, documentation, doctests for chebyshev_T, chebyshev_U and base class for ortho polys
    
    diff --git a/sage/functions/orthogonal_polys.py b/sage/functions/orthogonal_polys.py
    a b  
    261261   (it loads "specfun", but I'm not sure if that is the real reason).
    262262   The next call is usually faster but not always.
    263263
    264 TODO: Implement associated Legendre polynomials and Zernike
    265 polynomials. (Neither is in Maxima.)
    266 http://en.wikipedia.org/wiki/Associated_Legendre_polynomials
    267 http://en.wikipedia.org/wiki/Zernike_polynomials
     264.. TODO::
     265
     266    Implement associated Legendre polynomials and Zernike
     267    polynomials. (Neither is in Maxima.)
     268    :wikipedia:`Associated_Legendre_polynomials`
     269    :wikipedia:`Zernike_polynomials`
    268270
    269271REFERENCES:
    270272
    271 - Abramowitz and Stegun: Handbook of Mathematical Functions,
    272    http://www.math.sfu.ca/ cbm/aands/
     273.. [ASHandbook] Abramowitz and Stegun: Handbook of Mathematical Functions,
     274    http://www.math.sfu.ca/ cbm/aands/
    273275
    274 -  http://en.wikipedia.org/wiki/Chebyshev_polynomials
     276.. :wikipedia:`Chebyshev_polynomials`
    275277
    276 -  http://en.wikipedia.org/wiki/Legendre_polynomials
     278.. :wikipedia:`Legendre_polynomials`
    277279
    278 -  http://en.wikipedia.org/wiki/Hermite_polynomials
     280.. :wikipedia:`Hermite_polynomials`
    279281
    280 - http://mathworld.wolfram.com/GegenbauerPolynomial.html
     282.. http://mathworld.wolfram.com/GegenbauerPolynomial.html
    281283
    282 -  http://en.wikipedia.org/wiki/Jacobi_polynomials
     284.. :wikipedia:`Jacobi_polynomials`
    283285
    284 -  http://en.wikipedia.org/wiki/Laguerre_polynomia
     286.. :wikipedia:`Laguerre_polynomia`
    285287
    286 -  http://en.wikipedia.org/wiki/Associated_Legendre_polynomials
     288.. :wikipedia:`Associated_Legendre_polynomials`
    287289
     290.. [EffCheby] Wolfram Koepf: Effcient Computation of Chebyshev Polynomials
     291    in Computer Algebra
     292    Computer Algebra Systems: A Practical Guide.
     293    John Wiley, Chichester (1999): 79-99.
    288294
    289295AUTHORS:
    290296
    291297- David Joyner (2006-06)
     298- Stefan Reiterer (2010-)
    292299"""
    293300
    294301#*****************************************************************************
    295302#       Copyright (C) 2006 William Stein <wstein@gmail.com>
    296303#                     2006 David Joyner <wdj@usna.edu>
     304#                     2010 Stefan Reiterer <maldun.finsterschreck@gmail.com>
    297305#
    298306#  Distributed under the terms of the GNU General Public License (GPL)
    299307#
     
    307315#                  http://www.gnu.org/licenses/
    308316#*****************************************************************************
    309317
     318import warnings
     319
    310320from sage.misc.sage_eval import sage_eval
    311 from sage.rings.all import ZZ
     321from sage.rings.all import ZZ, RR, CC, RIF, CIF
    312322from sage.calculus.calculus import maxima
    313323
     324
     325from sage.symbolic.function import BuiltinFunction, GinacFunction, is_inexact
     326from sage.symbolic.expression import is_Expression
     327import sage.functions.special
     328from sage.functions.special import MaximaFunction, meval
     329from sage.functions.other import floor, gamma, factorial, abs, binomial
     330from sage.functions.other import sqrt, conjugate
     331from sage.functions.trig import sin, cos
     332from sage.functions.log import ln
     333import sage.symbolic.expression as expression
     334from sage.structure.parent import Parent
     335from sage.structure.coerce import parent
     336
    314337_done = False
    315338def _init():
    316339    """
     
    327350
    328351    Then after using one of these functions, it changes::
    329352
    330         sage: from sage.functions.orthogonal_polys import chebyshev_T
    331         sage: chebyshev_T(2,x)
    332         2*(x - 1)^2 + 4*x - 3
     353        sage: from sage.functions.orthogonal_polys import legendre_P
     354        sage: legendre_P(2,x)
     355        3/2*(x - 1)^2 + 3*x - 2
    333356        sage: sage.functions.orthogonal_polys._done
    334357        True
    335358
     359
    336360    Note that because here we use a Pynac variable ``x``,
    337361    the representation of the function is different from
    338362    its actual doctest, where a polynomial indeterminate
     
    348372    _done = True
    349373
    350374
    351 def chebyshev_T(n,x):
     375
     376class OrthogonalPolynomial(BuiltinFunction):
    352377    """
    353     Returns the Chebyshev function of the first kind for integers
    354     `n>-1`.
     378    Base class for orthogonal polynomials.
     379
     380    This class is an abstract base class for all orthogonal polynomials since
     381    they share similar properties. The evaluation as a polynomial
     382    is either done via maxima, or with pynac.
     383
     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.
     386    """
     387    def __init__(self, name, nargs=2, latex_name=None, conversions={}):
     388        """
     389        :class:`OrthogonalPolynomial` class needs the same input parameter as
     390        it's parent class.
     391       
     392        EXAMPLES::
     393
     394            sage: from sage.functions.orthogonal_polys import OrthogonalPolynomial
     395            sage: new = OrthogonalPolynomial('testo_P')
     396            sage: new
     397            testo_P
     398        """
     399        try:
     400            self._maxima_name = conversions['maxima']
     401        except KeyError:
     402            self._maxima_name = None
     403
     404        super(OrthogonalPolynomial,self).__init__(name=name, nargs=nargs,
     405                                 latex_name=latex_name, conversions=conversions)
     406
     407    def _maxima_init_evaled_(self, *args):
     408        r"""
     409        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.
     430
     431        EXAMPLES::
     432
     433            sage: from sage.functions.orthogonal_polys import OrthogonalPolynomial
     434            sage: new = OrthogonalPolynomial('testo_P')
     435            sage: new._apply_formula_(1,2.0)
     436            Traceback (most recent call last):
     437            ...
     438            NotImplementedError: no recursive calculation of values implemented
     439        """
     440        raise NotImplementedError("no recursive calculation of values implemented")
     441
     442    def _eval_special_values_(self,*args):
     443        """
     444        Evaluate the polynomial explicitly for special values.
     445
     446        EXAMPLES::
     447           
     448            sage: var('n')
     449            n
     450            sage: chebyshev_T(n,-1)
     451            (-1)^n
     452        """
     453        raise ValueError("no special values known")
     454
     455    def _eval_(self, *args):
     456        """
     457        The _eval_ method decides which evaluation suits best
     458        for the given input, and returns a proper value.
     459       
     460        EXAMPLES::
     461
     462            sage: chebyshev_T(5,x)
     463            16*x^5 - 20*x^3 + 5*x 
     464            sage: var('n')
     465            n
     466            sage: chebyshev_T(n,-1)
     467            (-1)^n
     468            sage: chebyshev_T(-7,x)
     469            64*x^7 - 112*x^5 + 56*x^3 - 7*x
     470            sage: chebyshev_T(3/2,x)
     471            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
     477            sage: parent(chebyshev_T(4, RIF(5)))
     478            Real Interval Field with 53 bits of precision
     479            sage: RR2 = RealField(5)
     480            sage: chebyshev_T(100000,RR2(2))
     481            8.9e57180
     482            sage: chebyshev_T(5,Qp(3)(2))
     483            2 + 3^2 + 3^3 + 3^4 + 3^5 + O(3^20)
     484            sage: chebyshev_T(100001/2, 2)
     485            doctest:500: RuntimeWarning: Warning: mpmath returns NoConvergence exception! Use other method instead.
     486            chebyshev_T(100001/2, 2)
     487            sage: chebyshev_U._eval_(1.5, Mod(8,9)) is None
     488            True
     489        """
     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
     503
     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
     508
     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:
     519                return None
     520
     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:
     528            return None
     529
     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        """
     537        This overides the call method from SageObject to avoid problems with coercions,
     538        since the _eval_ method is able to handle more data types than symbolic functions
     539        would normally allow.
     540        Thus we have the distinction between algebraic objects (if n is an integer),
     541        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
     548            sage: chebyshev_T(5,MatrixSpace(ZZ, 2)([1, 2, -4, 7]))
     549            [-40799  44162]
     550            [-88324  91687]
     551            sage: R.<x> = QQ[]
     552            sage: parent(chebyshev_T(5, x))
     553            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)
     564
     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.
     570
     571        EXAMPLES::
     572       
     573            sage: chebyshev_T._old_maxima_(-7,x) is None
     574            True
     575        """
     576        None
     577
     578class Func_chebyshev_T(OrthogonalPolynomial):
     579    """
     580    Chebyshev polynomials of the first kind.
     581
     582    REFERENCE:
     583
     584    - [ASHandbook]_ 22.5.31 page 778 and 6.1.22 page 256.
     585
     586    EXAMPLES::
     587
     588       sage: chebyshev_T(3,x)
     589       4*x^3 - 3*x
     590       sage: chebyshev_T(5,x)
     591       16*x^5 - 20*x^3 + 5*x
     592       sage: var('k')
     593       k
     594       sage: test = chebyshev_T(k,x)
     595       sage: test
     596       chebyshev_T(k, x)
     597    """
     598    def __init__(self):
     599        """
     600        Init method for the chebyshev polynomials of the first kind.
     601
     602        EXAMPLES::
     603
     604            sage: from sage.functions.orthogonal_polys import Func_chebyshev_T
     605            sage: chebyshev_T2 = Func_chebyshev_T()
     606            sage: chebyshev_T2(1,x)
     607            x
     608        """
     609        super(Func_chebyshev_T,self).__init__("chebyshev_T", nargs=2,
     610                                      conversions=dict(maxima='chebyshev_t',
     611                                                       mathematica='ChebyshevT'))
     612   
     613    def _eval_special_values_(self,*args):
     614        """
     615        Values known for special values of x.
     616        For details see [ASHandbook]_ 22.4 (p. 777)
     617
     618        EXAMPLES:
     619
     620            sage: var('n')
     621            n
     622            sage: chebyshev_T(n,1)
     623            1
     624            sage: chebyshev_T(n,-1)
     625            (-1)^n
     626            sage: chebyshev_T(-7, x) - chebyshev_T(7,x)
     627            0
     628            sage: chebyshev_T._eval_special_values_(3/2,x)
     629            Traceback (most recent call last):
     630            ...
     631            ValueError: No special values for non integral indices!
     632            sage: chebyshev_T._eval_special_values_(n, 0.1)
     633            Traceback (most recent call last):
     634            ...
     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!
     640        """
     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]
     649
     650        if (args[-1] == 0 and args[-1] in CC):
     651            return (1+(-1)**args[0])*(-1)**(args[0]/2)/2
     652
     653        if args[0] < 0 and args[0] in ZZ:
     654            return self._eval_(-args[0],args[-1])
     655
     656        raise ValueError("Value not found!")
     657   
     658    def _evalf_(self, *args,**kwds):
     659        """
     660        Evaluates :class:`chebyshev_T` numerically with mpmath.
     661        If the index is an integer we use the recursive formula since
     662        it is faster.
     663
     664        EXAMPLES::
     665
     666            sage: chebyshev_T(10,3).n(75)
     667            2.261953700000000000000e7
     668            sage: chebyshev_T(10,I).n()
     669            -3363.00000000000
     670            sage: chebyshev_T(5,0.3).n()
     671            0.998880000000000
     672            sage: chebyshev_T(1/2, 0)
     673            0.707106781186548
     674            sage: chebyshev_T._evalf_(1.5, Mod(8,9))
     675            Traceback (most recent call last):
     676            ...
     677            ValueError: No compatible type!
     678
     679        """
     680        if args[0] in ZZ and args[0] >= 0:
     681            return self._cheb_recur_(*args)[0]
     682       
     683        try:
     684            real_parent = kwds['parent']
     685        except KeyError:
     686            real_parent = parent(args[-1])
     687
     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
     702
     703        if not x_set:
     704            raise ValueError("No compatible type!")
     705               
     706        from sage.libs.mpmath.all import call as mpcall
     707        from sage.libs.mpmath.all import chebyt as mpchebyt
     708
     709        return mpcall(mpchebyt,args[0],x,parent=step_parent)
     710
     711    def _maxima_init_evaled_(self, *args):
     712        """
     713        Evaluate the Chebyshev polynomial ``self`` with maxima.
     714
     715        EXAMPLES::
     716
     717            sage: chebyshev_T._maxima_init_evaled_(1,x)
     718            'x'
     719            sage: var('n')
     720            n
     721            sage: maxima(chebyshev_T(n,x))
     722            chebyshev_t(n,x)
     723
     724        """
     725        n = args[0]
     726        x = args[1]
     727        return maxima.eval('chebyshev_t({0},{1})'.format(n,x))
     728
     729       
     730    def _apply_formula_(self,*args):
     731        """
     732        Applies explicit formulas for :class:`chebyshev_T`.
     733        This is much faster for numerical evaluation than maxima!
     734        See [ASHandbook]_ 227 (p. 782) for details for the recurions.
     735        See also [EffCheby]_ for fast evaluation techniques.
     736
     737        EXAMPLES::
     738
     739            sage: chebyshev_T._apply_formula_(2,0.1) == chebyshev_T._evalf_(2,0.1)
     740            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)
     744            512*x^10 - 1280*x^8 + 1120*x^6 - 400*x^4 + 50*x^2 - 1
     745
     746        """
     747        k = args[0]
     748        x = args[1]
     749
     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
     765        else:
     766            return self._cheb_recur_(k,x)[0]
     767
     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
     788
     789
     790    def _eval_numpy_(self, *args):
     791        """
     792        Evaluate ``self`` using numpy.
     793
     794        EXAMPLES::
     795
     796            sage: import numpy
     797            sage: z = numpy.array([1,2])
     798            sage: z2 = numpy.array([[1,2],[1,2]])
     799            sage: z3 = numpy.array([1,2,3.])
     800            sage: chebyshev_T(1,z)
     801            array([1, 2])
     802            sage: chebyshev_T(1,z2)
     803            array([[1, 2],
     804                   [1, 2]])
     805            sage: chebyshev_T(1,z3)
     806            array([ 1.,  2.,  3.])
     807            sage: chebyshev_T(z,0.1)
     808            array([ 0.1 , -0.98])
     809        """
     810        from scipy.special import eval_chebyt
     811        return eval_chebyt(args[0],args[-1])
     812
     813    def _derivative_(self, *args, **kwds):
     814        """
     815        Return the derivative of :class:`chebyshev_T` in form of the Chebyshev
     816        polynomial of the second kind :class:`chebyshev_U`.
     817
     818        EXAMPLES::
     819
     820            sage: var('k')
     821            k
     822            sage: derivative(chebyshev_T(k,x),x)
     823            k*chebyshev_U(k - 1, x)
     824            sage: derivative(chebyshev_T(3,x),x)
     825            12*x^2 - 3
     826            sage: derivative(chebyshev_T(k,x),k)
     827            Traceback (most recent call last):
     828            ...
     829            NotImplementedError: derivative w.r.t. to the index is not supported yet
     830        """
     831        diff_param = kwds['diff_param']
     832        if diff_param == 0:
     833            raise NotImplementedError("derivative w.r.t. to the index is not supported yet")
     834
     835        return args[0]*chebyshev_U(args[0]-1,args[1])
     836   
     837chebyshev_T = Func_chebyshev_T()
     838
     839class Func_chebyshev_U(OrthogonalPolynomial):
     840    """
     841    Class for the Chebyshev polynomial of the second kind.
    355842   
    356843    REFERENCE:
    357844
    358     - AS 22.5.31 page 778 and AS 6.1.22 page 256.
    359    
    360     EXAMPLES::
    361    
    362         sage: x = PolynomialRing(QQ, 'x').gen()
    363         sage: chebyshev_T(2,x)
    364         2*x^2 - 1
    365     """
    366     _init()
    367     return sage_eval(maxima.eval('chebyshev_t(%s,x)'%ZZ(n)), locals={'x':x})
    368 
    369 def chebyshev_U(n,x):
    370     """
    371     Returns the Chebyshev function of the second kind for integers `n>-1`.
    372    
    373     REFERENCE:
    374 
    375     - AS, 22.8.3 page 783 and AS 6.1.22 page 256.
     845    - [ASHandbook]_ 22.8.3 page 783 and 6.1.22 page 256.
    376846   
    377847    EXAMPLES::
    378848   
    379849        sage: x = PolynomialRing(QQ, 'x').gen()
    380850        sage: chebyshev_U(2,x)
    381851        4*x^2 - 1
     852        sage: chebyshev_U(3,x)
     853        8*x^3 - 4*x
    382854    """
    383     _init()
    384     return sage_eval(maxima.eval('chebyshev_u(%s,x)'%ZZ(n)), locals={'x':x})
     855    def __init__(self):
     856        """
     857        Init method for the chebyshev polynomials of the second kind.
     858
     859        EXAMPLES::
     860       
     861            sage: from sage.functions.orthogonal_polys import Func_chebyshev_U
     862            sage: chebyshev_U2 = Func_chebyshev_U()
     863            sage: chebyshev_U2(1,x)
     864            2*x
     865        """
     866        OrthogonalPolynomial.__init__(self, "chebyshev_U", nargs=2,
     867                                      conversions=dict(maxima='chebyshev_u',
     868                                                       mathematica='ChebyshevU'))
     869       
     870    def _apply_formula_(self,*args):
     871        """
     872        Applies explicit formulas for :class:`chebyshev_U`.
     873        This is much faster for numerical evaluation than maxima.
     874        See [ASHandbook]_ 227 (p. 782) for details on the recurions.
     875        See also [EffCheby]_ for the recursion formulas.
     876
     877        EXAMPLES::
     878
     879            sage: chebyshev_U._apply_formula_(2,0.1) == chebyshev_U._evalf_(2,0.1)
     880            True
     881        """
     882        k = args[0]
     883        x = args[1]
     884           
     885        if k == 0:
     886            return 1
     887        if k == 1:
     888            return 2*x
     889
     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           
     900        else:
     901            return self._cheb_recur_(k,x)[0]
     902
     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)
     930
     931    def _maxima_init_evaled_(self, *args):
     932        """
     933        Uses maxima to evaluate ``self``.
     934
     935        EXAMPLES::
     936
     937            sage: maxima(chebyshev_U(5,x))
     938            32*x^5-32*x^3+6*x
     939            sage: var('n')
     940            n
     941            sage: maxima(chebyshev_U(n,x))
     942            chebyshev_u(n,x)
     943            sage: maxima(chebyshev_U(2,x))
     944            4*x^2-1
     945        """
     946        n = args[0]
     947        x = args[1]
     948        return maxima.eval('chebyshev_u({0},{1})'.format(n,x))
     949
     950    def _evalf_(self, *args,**kwds):
     951        """
     952        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
     956        EXAMPLES::
     957
     958            sage: chebyshev_U(5,-4+3.*I)
     959            98280.0000000000 - 11310.0000000000*I
     960            sage: chebyshev_U(10,3).n(75)
     961            4.661117900000000000000e7
     962            sage: chebyshev_U._evalf_(1.5, Mod(8,9))
     963            Traceback (most recent call last):
     964            ...
     965            ValueError: No compatible type!
     966        """
     967        if args[0] in ZZ and args[0] >= 0:
     968            return self._cheb_recur_(*args)[0]
     969        try:
     970            real_parent = kwds['parent']
     971        except KeyError:
     972            real_parent = parent(args[-1])
     973
     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
     988
     989        if not x_set:
     990            raise ValueError("No compatible type!")
     991
     992        from sage.libs.mpmath.all import call as mpcall
     993        from sage.libs.mpmath.all import chebyu as mpchebyu
     994
     995        return mpcall(mpchebyu,args[0],args[-1],parent = step_parent)
     996
     997    def _eval_special_values_(self,*args):
     998        """
     999        Special values that known. [ASHandbook]_ 22.4 (p.777).
     1000
     1001        EXAMPLES::
     1002       
     1003            sage: var('n')
     1004            n
     1005            sage: chebyshev_U(n,1)
     1006            n + 1
     1007            sage: chebyshev_U(n,-1)
     1008            (-1)^n*(n + 1)
     1009            sage: chebyshev_U._eval_special_values_(26, Mod(0,9))
     1010            Traceback (most recent call last):
     1011            ...
     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
     1029        """
     1030        if (not is_Expression(args[0])) and (not args[0] in ZZ):
     1031            raise ValueError("No special values for non integral indices!")
     1032
     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)
     1041
     1042        if (args[-1] == 0 and args[-1] in CC):
     1043            return (1+(-1)**args[0])*(-1)**(args[0]/2)/2
     1044
     1045        if args[0] < 0 and args[0] in ZZ:
     1046            return -self._eval_(-args[0]-2,args[-1])
     1047
     1048        raise ValueError("Value not found!")
     1049
     1050    def _eval_numpy_(self, *args):
     1051        """
     1052        Evaluate ``self`` using numpy.
     1053
     1054        EXAMPLES::
     1055
     1056            sage: import numpy
     1057            sage: z = numpy.array([1,2])
     1058            sage: z2 = numpy.array([[1,2],[1,2]])
     1059            sage: z3 = numpy.array([1,2,3.])
     1060            sage: chebyshev_U(1,z)
     1061            array([2, 4])
     1062            sage: chebyshev_U(1,z2)
     1063            array([[2, 4],
     1064                   [2, 4]])
     1065            sage: chebyshev_U(1,z3)
     1066            array([ 2.,  4.,  6.])
     1067            sage: chebyshev_U(z,0.1)
     1068            array([ 0.2 , -0.96])
     1069        """
     1070        from scipy.special import eval_chebyu
     1071        return eval_chebyu(args[0],args[1])
     1072
     1073
     1074    def _derivative_(self, *args, **kwds):
     1075        """
     1076        Return the derivative of :class:`chebyshev_U` in form of the Chebyshev
     1077        polynomials of the first and second kind.
     1078       
     1079        EXAMPLES::
     1080
     1081            sage: var('k')
     1082            k
     1083            sage: derivative(chebyshev_U(k,x),x)
     1084            ((k + 1)*chebyshev_T(k + 1, x) - x*chebyshev_U(k, x))/(x^2 - 1)
     1085            sage: derivative(chebyshev_U(3,x),x)
     1086            24*x^2 - 4
     1087            sage: derivative(chebyshev_U(k,x),k)
     1088            Traceback (most recent call last):
     1089            ...
     1090            NotImplementedError: derivative w.r.t. to the index is not supported yet
     1091        """
     1092        diff_param = kwds['diff_param']
     1093        if diff_param == 0:
     1094            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)
     1098
     1099chebyshev_U = Func_chebyshev_U()
     1100
    3851101
    3861102def gen_laguerre(n,a,x):
    3871103    """
    3881104    Returns the generalized Laguerre polynomial for integers `n > -1`.
    389     Typically, a = 1/2 or a = -1/2.
     1105    Typically, `a = 1/2` or `a = -1/2`.
    3901106   
    391     REFERENCE:
     1107    REFERENCES:
    3921108
    393     - table on page 789 in AS.
     1109    - Table on page 789 in [ASHandbook]_.
    3941110   
    3951111    EXAMPLES::
    3961112   
     
    4931209   
    4941210    REFERENCE:
    4951211
    496     - AS 22.5.40 and 22.5.41, page 779.
     1212    - [ASHandbook]_ 22.5.40 and 22.5.41, page 779.
    4971213   
    4981214    EXAMPLES::
    4991215   
     
    5281244   
    5291245    REFERENCE:
    5301246
    531     - table on page 789 in AS.
     1247    - Table on page 789 in [ASHandbook]_.
    5321248   
    5331249    EXAMPLES::
    5341250   
     
    5431259
    5441260def laguerre(n,x):
    5451261    """
    546     Returns the Laguerre polynomial for integers `n > -1`.
     1262    Return the Laguerre polynomial for integers `n > -1`.
    5471263   
    5481264    REFERENCE:
    5491265
    550     - AS 22.5.16, page 778 and AS page 789.
     1266    - [ASHandbook]_ 22.5.16, page 778 and page 789.
    5511267   
    5521268    EXAMPLES::
    5531269   
     
    5691285   
    5701286    REFERENCE:
    5711287
    572     - AS 22.5.35 page 779.
     1288    - [ASHandbook]_ 22.5.35 page 779.
    5731289   
    5741290    EXAMPLES::
    5751291   
     
    5921308def legendre_Q(n,x):
    5931309    """
    5941310    Returns the Legendre function of the second kind for integers
    595     `n>-1`.
     1311    `n > -1`.
    5961312   
    5971313    Computed using Maxima.
    5981314   
     
    6201336   
    6211337    REFERENCE:
    6221338
    623     - AS 22.5.27
     1339    - [ASHandbook]_ 22.5.27
    6241340   
    6251341    EXAMPLES::
    6261342