Ticket #9706: trac_9706_orthogonal_polys.patch

File trac_9706_orthogonal_polys.patch, 78.3 KB (added by maldun, 11 years ago)

Latest version of orthogonal polys with scipy support, and changed doctests. Tested in sage-4.6.alpha3

  • sage/functions/orthogonal_polys.py

    # HG changeset patch
    # User maldun <domors@gmx.net>
    # Date 1287233203 -7200
    # Node ID 5fac168ed1fcfdbe0ec7712fd3bbf1309de1eef3
    # Parent  7ae7b5175304fea5b08c9ee313c372a3936baa1e
    Trac 9706: Added new version of orthogonal polynomials. The new versions
               are now BuiltInFunctions instead wrappers arround maxima.
               The numerical evaluation is done with mpmath.
               Symbolic evaluation is done with pynac, instead of maxima for the most cases (except some
               of the legendre functions)
               Many doctest had to be changed, since new functions are now integrated in pynac.
               Some of the old doctests had to be removed due the fact some things don't work in pynac
               like coercon from RR to the symbolic ring.
    
    diff -r 7ae7b5175304 -r 5fac168ed1fc sage/functions/orthogonal_polys.py
    a b  
     1# -*- coding: utf-8 -*-
    12r"""
    23Orthogonal Polynomials
    34
     
    89Maxima-related bug reports and comments on this module to
    910willisb@unk.edu. In your report, please include Maxima and specfun
    1011version information.
     12-----------------------------------------------------------------
     13Update (2010): The orthogonal polynomials are now selfstanding classes.
     14Altough some of the use the old maxima wrappers, most of the evaluation
     15is done with the help of direct methods. In most cases direct computation
     16with pynac works a lot better. Also the orthogonal polys
     17can be handled as symbolic expressions, if the input is not known,
     18what was not possible in maxima.
    1119
    1220
    1321-  The Chebyshev polynomial of the first kind arises as a solution
     
    289297AUTHORS:
    290298
    291299- David Joyner (2006-06)
     300- Stefan Reiterer (2010-)
    292301"""
    293302
    294303#*****************************************************************************
    295304#       Copyright (C) 2006 William Stein <wstein@gmail.com>
    296305#                     2006 David Joyner <wdj@usna.edu>
     306#                     2010 Stefan Reiterer <stefan.reiterer@uni-graz.at>
    297307#
    298308#  Distributed under the terms of the GNU General Public License (GPL)
    299309#
     
    314324from sage.rings.rational_field import RationalField
    315325from sage.rings.real_mpfr import RealField
    316326from sage.misc.sage_eval import sage_eval
    317 from sage.rings.all import QQ, ZZ, CDF, RDF
     327from sage.rings.all import QQ, ZZ, RR, CDF, RDF
     328from sage.symbolic.ring import SR
    318329import sage.rings.commutative_ring as commutative_ring
    319330import sage.rings.ring as ring
    320331from sage.calculus.calculus import maxima
     332from sage.symbolic.function import BuiltinFunction, GinacFunction, is_inexact
     333from sage.symbolic.expression import is_Expression
     334import sage.functions.special
     335from sage.functions.special import MaximaFunction, meval
     336from sage.functions.other import floor, gamma, factorial, abs, binomial
     337from sage.functions.other import sqrt, conjugate
     338from sage.functions.trig import sin, cos
     339from sage.functions.log import ln
     340import sage.symbolic.expression as expression
     341from sage.structure.parent import Parent
     342from sage.structure.coerce import parent
     343
     344from numpy import array as nparray
     345
     346from scipy.special import eval_chebyt, eval_legendre, eval_chebyu
     347from scipy.special import eval_gegenbauer, eval_genlaguerre, eval_hermite
     348from scipy.special import eval_jacobi, eval_laguerre, lpmv
     349
     350#Those imports are commented out, and are found in the _evalf_
     351#functions of the ortho polys. The reason is a problem with
     352#mpmath. gen_legendre_P, and _gen_legendre_Q start to throw
     353#exceptions if mpmath is included globally. Wehn it is done
     354#locally no problems occour so far.
     355#See ticket Nr. 9706 for more details. They should be exchanged,
     356#when this problem with mpmath is solved, because of performance
     357#reasons.
     358 
     359#from sage.libs.mpmath.all import call as mpcall
     360#from sage.libs.mpmath.all import legendre as mplegendre
     361#from sage.libs.mpmath.all import legenp as mplegenp
     362#from sage.libs.mpmath.all import chebyt as mpchebyt
     363#from sage.libs.mpmath.all import chebyu as mpchebyu
     364#from sage.libs.mpmath.all import gegenbauer as mpgegenbauer
     365#from sage.libs.mpmath.all import laguerre as mplaguerre
     366#from sage.libs.mpmath.all import hermite as mphermite
     367#from sage.libs.mpmath.all import jacobi as mpjacobi
     368#from sage.libs.mpmath.all import legenq as mplegenq
     369
    321370
    322371_done = False
    323372def _init():
     
    333382        sage: sage.functions.orthogonal_polys._done
    334383        False
    335384
    336     Then after using one of these functions, it changes::
     385    Then after using one of these functions, it changes:: (The value is now
     386    False for chebyshev_T because chebyshev_T uses clenshaw method instead...)
    337387
    338388        sage: from sage.functions.orthogonal_polys import chebyshev_T
    339389        sage: chebyshev_T(2,x)
    340         2*(x - 1)^2 + 4*x - 3
     390        2*x^2 - 1
     391        sage: sage.functions.orthogonal_polys._done
     392        False
     393        sage: legendre_Q(1,x)
     394        1/2*x*log(-(x + 1)/(x - 1)) - 1
    341395        sage: sage.functions.orthogonal_polys._done
    342396        True
    343 
     397   
    344398    Note that because here we use a Pynac variable ``x``,
    345399    the representation of the function is different from
    346400    its actual doctest, where a polynomial indeterminate
     
    355409    maxima.eval('orthopoly_returns_intervals:false;')
    356410    _done = True
    357411
     412#load /home/maldun/sage/sage-4.5.2/devel/sage-ortho/sage/functions/orthogonal_polys.py
    358413
    359 def chebyshev_T(n,x):
     414class OrthogonalPolynomial(BuiltinFunction):
    360415    """
    361     Returns the Chebyshev function of the first kind for integers
    362     `n>-1`.
     416    Base Class for Orthogonal Polynomials. The evaluation as a polynomial
     417    is either done via maxima, or with pynac due to performance reasons. 
     418    Convention: The first argument is always the order of the polynomial,
     419    he last one is always the value x where the polynomial is evaluated.
     420
     421    """
     422    def __init__(self, name, nargs = 2, latex_name = None, conversions = {}):
     423        try:
     424            self._maxima_name = conversions['maxima']
     425        except KeyError:
     426            self._maxima_name = None
     427
     428        BuiltinFunction.__init__(self, name = name,
     429        nargs = nargs, latex_name = latex_name, conversions = conversions)
     430
     431    def _maxima_init_evaled_(self, *args):
     432        """
     433        Returns a string which represents this function evaluated at
     434        *args* in Maxima.
     435        In fact these are thought to be the old wrappers for the orthogonal
     436        polynomials. These are used when the other evaluation methods fail,
     437        or are not fast enough. It appears that direct computation
     438        with pynac is in most cases faster than maxima. Maxima comes into
     439        play when all other methods fail.
     440        A little switch does the trick...
     441
     442        EXAMPLES::
     443
     444            sage: chebyshev_T(3,x)
     445            4*x^3 - 3*x
     446        """
     447        return None
     448
     449    def _clenshaw_method_(self,*args):
     450        """
     451        The Clenshaw method uses the three term recursion of the polynomial,
     452        or explicit formulas instead of maxima to evaluate the polynomial
     453        efficiently, if the x argument is not a symbolic expression.
     454        The name comes from the Clenshaw algorithm for fast evaluation of
     455        polynomials in chebyshev form.
     456       
     457        comparison: Maxima and Clenshaw algorithm for non-symbolic evaluation:
     458         #sage: time chebyshev_T(50,10) #clenshaw
     459         #CPU times: user 0.00 s, sys: 0.00 s, total: 0.00 s
     460         #Wall time: 0.00 s
     461         #49656746733678312490954442369580252421769338391329426325400124999
     462         #sage: time sage.functions.orthogonal_polys.chebyshev_T(50,10)
     463         #maxima
     464         #CPU times: user 0.00 s, sys: 0.00 s, total: 0.00 s
     465         #Wall time: 0.05 s
     466         #49656746733678312490954442369580252421769338391329426325400124999
     467
     468         #sage: time chebyshev_T(500,10); #clenshaw
     469         #CPU times: user 0.00 s, sys: 0.00 s, total: 0.00 s
     470         #Wall time: 0.00 s
     471         #sage: time sage.functions.orthogonal_polys.chebyshev_T(500,10);
     472         #maxima
     473         #CPU times: user 0.11 s, sys: 0.00 s, total: 0.11 s
     474         #Wall time: 0.77 s
     475        """
     476        raise NotImplementedError(
     477            "No recursive calculation of values implemented (yet)!")
     478
     479    def _eval_special_values_(self,*args):
     480        """
     481        Evals the polynomial explicitly for special values.
     482        EXAMPLES:
     483           
     484            sage: var('n')
     485            n
     486            sage: chebyshev_T(n,-1)
     487            (-1)^n
     488        """
     489        raise ValueError("No special values known!")
     490           
     491
     492    def _eval_(self, *args):
     493        """
     494       
     495        The symbolic evaluation is either done with maxima, or with direct
     496        computaion in pynac, because it seems the evaluation done in maxima
     497        doesn't seem very clever..
     498        For the fast numerical evaluation an other method should be used...
     499        Therefore I suggest Clenshaw's algorithm, which uses the recursion!
     500        The function also checks for special values, and if
     501        the order is an integer and in range!
     502       
     503        performance:
     504        #sage: time chebyshev_T(5,x) #maxima
     505        #CPU times: user 0.00 s, sys: 0.00 s, total: 0.00 s
     506        #Wall time: 0.16 s
     507        #16*(x - 1)^5 + 80*(x - 1)^4 + 140*(x - 1)^3 + 100*(x - 1)^2 + 25*x - 24
     508
     509        #sage: time chebyshev_T(5,x) #clenshaw
     510        #CPU times: user 0.01 s, sys: 0.00 s, total: 0.01 s
     511        #Wall time: 0.01 s
     512        #16*x^5 - 20*x^3 + 5*x
     513
     514        #time chebyshev_T(50,x)
     515        #CPU times: user 0.03 s, sys: 0.00 s, total: 0.03 s
     516        #Wall time: 0.04 s
     517        #562949953421312*(x - 1)^50 + 28147497671065600*(x - 1)^49 +....
     518
     519        #time chebyshev_T(100,x);
     520        #CPU times: user 0.08 s, sys: 0.00 s, total: 0.08 s
     521        #Wall time: 0.08 s
     522
     523        EXAMPLES::
     524            sage: chebyshev_T(5,x)
     525            16*x^5 - 20*x^3 + 5*x 
     526            sage: var('n')
     527            n
     528            sage: chebyshev_T(n,-1)
     529            (-1)^n
     530            sage: chebyshev_T(-7,x)
     531            chebyshev_T(-7, x)
     532            sage: chebyshev_T(3/2,x)
     533            chebyshev_T(3/2, x)
     534   
     535        """
     536       
     537        if not is_Expression(args[0]):
     538           
     539            if not is_Expression(args[-1]) and is_inexact(args[-1]):
     540                try:
     541                    import sage.libs.mpmath.all as mpmath
     542                    return self._evalf_(*args)
     543                except AttributeError:
     544                    pass
     545                except mpmath.NoConvergence:
     546                    print "Warning: mpmath returns NoConvergence!"
     547                    print "Switching to clenshaw_method, but it \
     548                    may not be stable!"
     549                except ValueError:
     550                    pass
     551
     552            #A faster check would be nice...
     553            if args[0] != floor(args[0]):
     554                if not is_Expression(args[-1]):
     555                    try:
     556                        return self._evalf_(*args)
     557                    except AttributeError:
     558                        pass
     559                else:
     560                    return None
     561       
     562            if args[0] < 0:
     563                return None
     564           
     565           
     566        try:
     567            return self._eval_special_values_(*args)
     568        except ValueError:
     569            pass
     570
     571   
     572        if not is_Expression(args[0]):
     573           
     574            try:
     575                return self._clenshaw_method_(*args)
     576            except NotImplementedError:
     577                pass
     578       
     579        if self._maxima_name is None:
     580            return None
     581        else:
     582            _init()
     583            try:
     584                #s = maxima(self._maxima_init_evaled_(*args))
     585                #This above is very inefficient! The older
     586                #methods were much faster...
     587                return self._maxima_init_evaled_(*args)
     588            except TypeError:
     589                return None
     590            if self._maxima_name in repr(s):
     591                return None
     592            else:
     593                return s.sage()
     594
     595
     596class Func_chebyshev_T(OrthogonalPolynomial):
     597    """
     598    Class for the Chebyshev polynomial of the first kind.
    363599   
    364600    REFERENCE:
    365601
     
    367603   
    368604    EXAMPLES::
    369605   
    370         sage: x = PolynomialRing(QQ, 'x').gen()
    371         sage: chebyshev_T(2,x)
    372         2*x^2 - 1
     606       sage: chebyshev_T(3,x)
     607       4*x^3 - 3*x
     608       sage: chebyshev_T(5,x)
     609       16*x^5 - 20*x^3 + 5*x
     610       sage: var('k')
     611       k
     612       sage: test = chebyshev_T(k,x)
     613       sage: test
     614       chebyshev_T(k, x)
    373615    """
    374     _init()
    375     return sage_eval(maxima.eval('chebyshev_t(%s,x)'%ZZ(n)), locals={'x':x})
     616   
     617    def __init__(self):
     618        OrthogonalPolynomial.__init__(self,"chebyshev_T",nargs = 2,
     619conversions =dict(maxima='chebyshev_t',mathematica='ChebyshevT'))
     620   
     621    def _eval_special_values_(self,*args):
     622        """
     623        Values known for special values of x.
     624        For details see A.S. 22.4 (p. 777)
     625       
     626        EXAMPLES:
     627           
     628            sage: var('n')
     629            n
     630            sage: chebyshev_T(n,1)
     631            1
     632            sage: chebyshev_T(n,-1)
     633            (-1)^n
     634        """
     635        if args[-1] == 1:
     636            return 1
     637       
     638        if args[-1] == -1:
     639            return (-1)**args[0]
    376640
    377 def chebyshev_U(n,x):
     641        if (args[-1] == 0):
     642            return (1+(-1)**args[0])*(-1)**(args[0]/2)/2
     643
     644        raise ValueError("Value not found!")
     645   
     646    def _evalf_(self, *args,**kwds):
     647        """
     648        Evals chebyshev_T
     649        numerically with mpmath.
     650        EXAMPLES::
     651            sage: chebyshev_T(10,3).n(75)
     652            2.261953700000000000000e7
     653            sage: chebyshev_T(10,I).n()
     654            -3363.00000000000
     655            sage: chebyshev_T(5,0.3).n()
     656            0.998880000000000
     657        """
     658        try:
     659            step_parent = kwds['parent']
     660        except KeyError:
     661            step_parent = parent(args[-1])
     662
     663        try:
     664            precision = step_parent.prec()
     665        except AttributeError:
     666            precision = RR.prec()
     667
     668        from sage.libs.mpmath.all import call as mpcall
     669        from sage.libs.mpmath.all import chebyt as mpchebyt
     670
     671        return mpcall(mpchebyt,args[0],args[-1],prec = precision)
     672
     673    def _maxima_init_evaled_(self, *args):
     674        n = args[0]
     675        x = args[1]
     676        return sage_eval(maxima.eval('chebyshev_t(%s,x)'%ZZ(n)), locals={'x':x})
     677
     678    def _clenshaw_method_(self,*args):
     679            """
     680            Clenshaw method for chebyshev_T (means use recursions in this case)
     681            This is much faster for numerical evaluation than maxima!
     682            See A.S. 227 (p. 782) for details for the recurions
     683            """
     684
     685            k = args[0]
     686            x = args[1]
     687
     688            if k == 0:
     689                return 1
     690            elif k == 1:
     691                return x
     692            else:
     693                #TODO: When evaluation of Symbolic Expressions works better
     694                #use these explicit formulas instead!
     695                #if -1 <= x <= 1:
     696                #    return cos(k*acos(x))
     697                #elif 1 < x:
     698                #    return cosh(k*acosh(x))
     699                #else: # x < -1
     700                #    return (-1)**(k%2)*cosh(k*acosh(-x))
     701               
     702                help1 = 1
     703                help2 = x
     704                if is_Expression(x):
     705                    #raise NotImplementedError
     706                    help3 = 0
     707                    for j in xrange(0,floor(k/2)+1):
     708                        help3 = \
     709help3 +(-1)**j*(2*x)**(k-2*j)*factorial(k-j-1)/factorial(j)/factorial(k-2*j)
     710                    help3 = help3*k/2
     711                else:
     712                     for j in xrange(0,k-1):
     713                        help3 = 2*x*help2 - help1
     714                        help1 = help2
     715                        help2 = help3
     716               
     717                return help3
     718
     719    def _eval_numpy_(self, *args):
     720        """
     721        EXAMPLES::
     722        sage: import numpy
     723        sage: z = numpy.array([1,2])
     724        sage: z2 = numpy.array([[1,2],[1,2]])
     725        sage: z3 = numpy.array([1,2,3.])
     726        sage: chebyshev_T(1,z)
     727        array([ 1.,  2.])
     728        sage: chebyshev_T(1,z2)
     729        array([[ 1.,  2.],
     730               [ 1.,  2.]])
     731        sage: chebyshev_T(1,z3)
     732        array([ 1.,  2.,  3.])
     733        sage: chebyshev_T(z,0.1)
     734        array([ 0.1 , -0.98])
     735
     736        """
     737
     738        return eval_chebyt(args[0],args[-1])
     739
     740    def _derivative_(self, *args, **kwds):
     741        """
     742        Returns the derivative of chebyshev_T in form of the chebyshev Polynomial
     743        of the second kind chebyshev_U
     744        EXAMPLES::
     745            sage: var('k')
     746            k
     747            sage: derivative(chebyshev_T(k,x),x)
     748            k*chebyshev_U(k - 1, x)
     749            sage: derivative(chebyshev_T(3,x),x)
     750            12*x^2 - 3
     751            sage: derivative(chebyshev_T(k,x),k)
     752            Traceback (most recent call last):
     753            ...
     754            NotImplementedError: Derivative w.r.t. to the index is not supported, yet, and perhaps never will be...
     755 
     756        """
     757        diff_param = kwds['diff_param']
     758        if diff_param == 0:
     759            raise NotImplementedError(
     760"Derivative w.r.t. to the index is not supported, yet, \
     761and perhaps never will be...")
     762        else:
     763            return args[0]*chebyshev_U(args[0]-1,args[1])
     764
     765   
     766chebyshev_T = Func_chebyshev_T()
     767
     768class Func_chebyshev_U(OrthogonalPolynomial):
     769   
    378770    """
    379     Returns the Chebyshev function of the second kind for integers `n>-1`.
     771    Class for the Chebyshev polynomial of the second kind.
    380772   
    381773    REFERENCE:
    382774
     
    387779        sage: x = PolynomialRing(QQ, 'x').gen()
    388780        sage: chebyshev_U(2,x)
    389781        4*x^2 - 1
     782        sage: chebyshev_U(3,x)
     783        8*x^3 - 4*x
    390784    """
    391     _init()
    392     return sage_eval(maxima.eval('chebyshev_u(%s,x)'%ZZ(n)), locals={'x':x})
    393785
    394 def gen_laguerre(n,a,x):
     786    def __init__(self):
     787        OrthogonalPolynomial.__init__(self,"chebyshev_U",
     788nargs = 2,conversions =dict(maxima='chebyshev_u',mathematica='ChebyshevU'))
     789       
     790    def _clenshaw_method_(self,*args):
     791        """
     792        Clenshaw method for chebyshev_U (means use the recursion...)
     793        This is much faster for numerical evaluation than maxima!
     794        See A.S. 227 (p. 782) for details for the recurions
     795        """
     796        k = args[0]
     797        x = args[1]
     798           
     799        if k == 0:
     800            return 1
     801        elif k == 1:
     802            return 2*x
     803        else:
     804            help1 = 1
     805            help2 = 2*x
     806            if is_Expression(x):
     807                #raise NotImplementedError("Maxima is faster here!")
     808                help3 = 0
     809                for j in xrange(0,floor(k/2)+1):
     810                    help3 = \
     811help3 + (-1)**j*(2*x)**(k-2*j)*factorial(k-j)/factorial(j)/factorial(k-2*j)
     812               
     813            else:
     814                for j in xrange(0,k-1):
     815                    help3 = 2*x*help2 - help1
     816                    help1 = help2
     817                    help2 = help3
     818                   
     819            return help3   
     820   
     821    def _maxima_init_evaled_(self, *args):
     822        """
     823        Uses 
     824        EXAMPLES::
     825        sage: x = PolynomialRing(QQ, 'x').gen()
     826        sage: chebyshev_T(2,x)
     827        2*x^2 - 1
     828        """
     829        n = args[0]
     830        x = args[1]
     831        return sage_eval(maxima.eval('chebyshev_u(%s,x)'%ZZ(n)), locals={'x':x})
     832
     833       
     834    def _evalf_(self, *args,**kwds):
     835        """
     836        Evals chebyshev_U
     837        numerically with mpmath.
     838        EXAMPLES::
     839            sage: chebyshev_U(5,-4+3.*I)
     840            98280.0000000000 - 11310.0000000000*I
     841            sage: chebyshev_U(10,3).n(75)
     842            4.661117900000000000000e7
     843        """
     844        try:
     845            step_parent = kwds['parent']
     846        except KeyError:
     847            step_parent = parent(args[-1])
     848
     849        try:
     850            precision = step_parent.prec()
     851        except AttributeError:
     852            precision = RR.prec()
     853
     854        from sage.libs.mpmath.all import call as mpcall
     855        from sage.libs.mpmath.all import chebyu as mpchebyu
     856
     857        return mpcall(mpchebyu,args[0],args[-1],prec = precision)
     858
     859    def _eval_special_values_(self,*args):
     860        """
     861        Special values known. A.S. 22.4 (p.777).
     862        EXAMPLES:
     863       
     864            sage: var('n')
     865            n
     866            sage: chebyshev_U(n,1)
     867            n + 1
     868            sage: chebyshev_U(n,-1)
     869            (n + 1)*(-1)^n
     870        """
     871        if args[-1] == 1:
     872            return (args[0]+1)
     873       
     874        if args[-1] == -1:
     875            return (-1)**args[0]*(args[0]+1)
     876
     877        if (args[-1] == 0):
     878            return (1+(-1)**args[0])*(-1)**(args[0]/2)/2
     879
     880        raise ValueError("Value not found")
     881
     882    def _eval_numpy_(self, *args):
     883        """
     884        EXAMPLES::
     885        sage: import numpy
     886        sage: z = numpy.array([1,2])
     887        sage: z2 = numpy.array([[1,2],[1,2]])
     888        sage: z3 = numpy.array([1,2,3.])
     889        sage: chebyshev_U(1,z)
     890        array([ 2.,  4.])
     891        sage: chebyshev_U(1,z2)
     892        array([[ 2.,  4.],
     893               [ 2.,  4.]])
     894        sage: chebyshev_U(1,z3)
     895        array([ 2.,  4.,  6.])
     896        sage: chebyshev_U(z,0.1)
     897        array([ 0.2 , -0.96])
     898
     899        """
     900
     901        return eval_chebyu(args[0],args[1])
     902
     903
     904    def _derivative_(self, *args, **kwds):
     905        """
     906        Returns the derivative of chebyshev_U in form of the chebyshev
     907        Polynomials of the first and second kind
     908       
     909        EXAMPLES::
     910            sage: var('k')
     911            k
     912            sage: derivative(chebyshev_U(k,x),x)
     913            ((k + 1)*chebyshev_T(k + 1, x) - x*chebyshev_U(k, x))/(x^2 - 1)
     914            sage: derivative(chebyshev_U(3,x),x)
     915            24*x^2 - 4
     916            sage: derivative(chebyshev_U(k,x),k)
     917            Traceback (most recent call last):
     918            ...
     919            NotImplementedError: Derivative w.r.t. to the index is not supported, yet, and perhaps never will be...
     920 
     921        """
     922        diff_param = kwds['diff_param']
     923        if diff_param == 0:
     924            raise NotImplementedError(
     925"Derivative w.r.t. to the index is not supported, \
     926yet, and perhaps never will be...")
     927        else:
     928            return ((args[0]+1)*chebyshev_T(args[0]+1,args[1])-args[1]*
     929                    chebyshev_U(args[0],args[1]))/(args[1]**2-1)
     930       
     931chebyshev_U = Func_chebyshev_U()
     932
     933class Func_gegenbauer(OrthogonalPolynomial):
     934    """
     935    Returns the ultraspherical (or Gegenbauer) polynomial.
     936   
     937    REFERENCE:
     938
     939    - AS 22.5.27
     940   
     941    EXAMPLES::
     942   
     943        sage: x = PolynomialRing(QQ, 'x').gen()
     944        sage: ultraspherical(2,3/2,x)
     945        15/2*x^2 - 3/2
     946        sage: ultraspherical(2,1/2,x)
     947        3/2*x^2 - 1/2
     948        sage: ultraspherical(1,1,x)
     949        2*x     
     950        sage: t = PolynomialRing(RationalField(),"t").gen()
     951        sage: gegenbauer(3,2,t)
     952        32*t^3 - 12*t
     953    """   
     954    def __init__(self):
     955        OrthogonalPolynomial.__init__(self,"gegenbauer",nargs = 3,
     956        conversions =dict(maxima='ultraspherical',mathematica='GegenbauerC'))
     957
     958    def _clenshaw_method_(self,*args):
     959        """
     960        Clenshaw method for gegenbauer poly (means use the recursion...)
     961        This is much faster for numerical evaluation than maxima!
     962        See A.S. 227 (p. 782) for details for the recurions
     963        """
     964        k = args[0]
     965        x = args[-1]
     966        alpha = args[1]
     967
     968        #if is_Expression(alpha) or abs(alpha) > numpy.finfo(float).eps:
     969        #    alpha_zero = False
     970        #else:
     971        #    alpha_zero = True
     972        #This should be discussed. The case 0 is the chebyshev_T.
     973        #(see A.S. 22.4 (p.777)); But since all systems handle
     974        #it this way we perhaps should leave it.
     975
     976        if k == 0:
     977            return 1
     978        elif k == 1:
     979                return 2*x*alpha
     980        else:
     981            help1 = 1
     982            help2 = 2*x*alpha
     983            gamma_alpha = gamma(alpha)
     984
     985            help3 = 0
     986            if is_Expression(x):
     987                for j in xrange(0,floor(k/2)+1):
     988                    help3 = \
     989help3 + (-1)**j*gamma(alpha+k-j)/factorial(j)/factorial(k-2*j)*(2*x)**(k-2*j)/gamma_alpha
     990               
     991            else:
     992                for j in xrange(1,k):
     993                    help3 = 2*(j+alpha)*x*help2 - (j+2*alpha-1)*help1
     994                    help3 = help3/(j+1)
     995                    help1 = help2
     996                    help2 = help3
     997                   
     998            return help3   
     999
     1000    def _maxima_init_evaled_(self, *args):
     1001        """
     1002        Uses 
     1003        EXAMPLES::
     1004        sage: x = PolynomialRing(QQ, 'x').gen()
     1005        sage: ultraspherical(2,1/2,x)
     1006        3/2*x^2 - 1/2
     1007        """
     1008        n = args[0]
     1009        a = args[1]
     1010        x = args[2]
     1011        return sage_eval(maxima.eval('ultraspherical(%s,%s,x)'%(ZZ(n),a)),\
     1012                         locals={'x':x})
     1013
     1014
     1015    def _evalf_(self, *args,**kwds):
     1016        """
     1017        Evals gegenbauer
     1018        numerically with mpmath.
     1019        EXAMPLES::
     1020            sage: gegenbauer(10,2,3.).n(54)
     1021            5.25360702000000e8
     1022        """
     1023           
     1024        try:
     1025            step_parent = kwds['parent']
     1026        except KeyError:
     1027            step_parent = parent(args[-1])
     1028
     1029        try:
     1030            precision = step_parent.prec()
     1031        except AttributeError:
     1032            precision = RR.prec()
     1033
     1034        from sage.libs.mpmath.all import call as mpcall
     1035        from sage.libs.mpmath.all import gegenbauer as mpgegenbauer
     1036
     1037        return mpcall(
     1038            mpgegenbauer,args[0],args[1],args[-1],prec = precision)
     1039
     1040    def _eval_special_values_(self,*args):
     1041        """
     1042        Special values known. A.S. 22.4 (p.777)
     1043        EXAMPLES:
     1044       
     1045            sage: var('n a')
     1046            (n, a)
     1047            sage: gegenbauer(n,1/2,x)
     1048            legendre_P(n, x)
     1049            sage: gegenbauer(n,1,x)
     1050            chebyshev_U(n, x)
     1051            sage: gegenbauer(n,a,1)
     1052            binomial(2*a + n - 1, n)
     1053            sage: gegenbauer(n,a,-1)
     1054            (-1)^n*binomial(2*a + n - 1, n)
     1055            sage: gegenbauer(n,a,0)
     1056            1/2*((-1)^n + 1)*(-1)^(1/2*n)*gamma(a + 1/2*n)/(gamma(a)*gamma(1/2*n))
     1057        """
     1058        #if args[1] == 0 and args[0] != 0:
     1059        #    return args[0]*chebyshev_T(args[0],args[-1])/2
     1060        #This should be discussed
     1061
     1062        if args[1] == 0.5:
     1063            return legendre_P(args[0],args[-1])
     1064
     1065        if args[1] == 1:
     1066            return chebyshev_U(args[0],args[-1])
     1067
     1068        if args[-1] == 1:
     1069            return binomial(args[0] + 2*args[1] - 1,args[0])
     1070       
     1071        if args[-1] == -1:
     1072            return (-1)**args[0]*binomial(args[0] + 2*args[1] - 1,args[0])
     1073
     1074        if (args[-1] == 0):
     1075            return (1+(-1)**args[0])*(-1)**(args[0]/2)/2*\
     1076            gamma(args[1]+args[0]/2)/gamma(args[1])/gamma(args[0]/2)
     1077
     1078        raise ValueError("Value not found")
     1079
     1080    def _eval_numpy_(self, *args):
     1081        """
     1082        EXAMPLES::
     1083        sage: import numpy
     1084        sage: z = numpy.array([1,2])
     1085        sage: z2 = numpy.array([[1,2],[1,2]])
     1086        sage: z3 = numpy.array([1,2,3.])
     1087        sage: gegenbauer(2,0,z)
     1088        array([  2.78134232e-309,   1.94693963e-308])
     1089        sage: gegenbauer(2,1,z)
     1090        array([  3.,  15.])
     1091        sage: gegenbauer(2,1,z2)
     1092        array([[  3.,  15.],
     1093               [  3.,  15.]])
     1094        sage: gegenbauer(2,1,z3)
     1095        array([  3.,  15.,  35.])
     1096        """
     1097
     1098        return eval_gegenbauer(args[0],args[1],args[-1])
     1099
     1100    def _derivative_(self, *args, **kwds):
     1101        """
     1102        Returns the derivative of chebyshev_U in form of the chebyshev
     1103        Polynomials of the first and second kind
     1104       
     1105        EXAMPLES::
     1106            sage: var('k a')
     1107            (k, a)
     1108            sage: derivative(gegenbauer(k,a,x),x)
     1109            (k*x*gegenbauer(k, a, x) - (2*a + k - 1)*gegenbauer(k - 1, a, x))/(x^2 - 1)
     1110            sage: derivative(gegenbauer(4,3,x),x)
     1111            960*x^3 - 240*x
     1112            sage: derivative(gegenbauer(k,a,x),a)
     1113            Traceback (most recent call last):
     1114            ...
     1115            NotImplementedError: Derivative w.r.t. to the index is not supported, yet,and perhaps never will be...
     1116        """
     1117        diff_param = kwds['diff_param']
     1118        if diff_param in [0,1]:
     1119            raise NotImplementedError(
     1120"Derivative w.r.t. to the index is not supported, yet,\
     1121and perhaps never will be...")
     1122        else:
     1123            return (-args[0]*args[-1]*gegenbauer(args[0],args[1],args[2])+\
     1124(args[0] + 2*args[1]-1)*gegenbauer(args[0]-1,args[1],args[2]))/(1-args[-1]**2)
     1125
     1126gegenbauer = Func_gegenbauer()
     1127ultraspherical = Func_gegenbauer()
     1128
     1129class Func_gen_laguerre(OrthogonalPolynomial):
    3951130    """
    3961131    Returns the generalized Laguerre polynomial for integers `n > -1`.
    3971132    Typically, a = 1/2 or a = -1/2.
     
    4141149        sage: gen_laguerre(3,0,x)
    4151150        -1/6*x^3 + 3/2*x^2 - 3*x + 1
    4161151    """
    417     from sage.functions.all import sqrt
    418     _init()
    419     return sage_eval(maxima.eval('gen_laguerre(%s,%s,x)'%(ZZ(n),a)), locals={'x':x})
     1152    def __init__(self):
     1153        OrthogonalPolynomial.__init__(self,"gen_laguerre",nargs = 3,
     1154conversions =dict(maxima='gen_laguerre',mathematica='LaguerreL'))
     1155       
     1156    def _clenshaw_method_(self,*args):
     1157        """
     1158        Clenshaw method for gen_laguerre polynomial (means use the recursion...)
     1159        This is much faster for numerical evaluation than maxima!
     1160        See A.S. 227 (p. 782) for details for the recurions.
     1161        Warning: The clanshaw method for the laguerre Polynomials
     1162        should only used for exact data types, when high orders are
     1163        used, due to weak instabilities of the recursion!
     1164        """
     1165        k = args[0]
     1166        a = args[1]
     1167        x = args[-1]
     1168       
     1169        if k == 0:
     1170            return 1
     1171        elif k == 1:
     1172            return 1-x + a
     1173        else:
     1174            help1 = 1
     1175            help2 = 1-x + a
     1176            if is_Expression(x):
     1177                help3 = 0
     1178                for j in xrange(0,k+1):
     1179                    help3 = help3 + (-1)**j*binomial(k+a,k-j)/factorial(j)*x**j
     1180            else:
     1181                for j in xrange(1,k):
     1182                    help3 = -x*help2 + (2*j+1+a)*help2 - (j+a)*help1
     1183                    help3 = help3/(j+1)
     1184                    help1 = help2
     1185                    help2 = help3
     1186               
     1187            return help3   
    4201188
    421 def gen_legendre_P(n,m,x):
    422     r"""
    423     Returns the generalized (or associated) Legendre function of the
    424     first kind for integers `n > -1, m > -1`.
     1189    def _maxima_init_evaled_(self, *args):
     1190        n = args[0]
     1191        a = args[1]
     1192        x = args[-1]
     1193       
     1194        from sage.functions.all import sqrt
     1195        _init()
     1196        return sage_eval(maxima.eval('gen_laguerre(%s,%s,x)'%(ZZ(n),a)),\
     1197                         locals={'x':x})
     1198
     1199    def _evalf_(self, *args,**kwds):
     1200        """
     1201        Evals generalized laguerre polynomial
     1202        numerically with mpmath.
     1203        EXAMPLES::
     1204            sage: gen_laguerre(3,2,5.*I).n(53)
     1205            -52.5000000000000 - 29.1666666666667*I
     1206        """
     1207        try:
     1208            step_parent = kwds['parent']
     1209        except KeyError:
     1210            step_parent = parent(args[-1])
     1211
     1212        try:
     1213            precision = step_parent.prec()
     1214        except AttributeError:
     1215            precision = RR.prec()
     1216
     1217        from sage.libs.mpmath.all import call as mpcall
     1218        from sage.libs.mpmath.all import laguerre as mplaguerre
     1219
     1220        return mpcall(
     1221            mplaguerre,args[0],args[1],args[-1],prec = precision)
    4251222   
    426     The awkward code for when m is odd and 1 results from the fact that
    427     Maxima is happy with, for example, `(1 - t^2)^3/2`, but
    428     Sage is not. For these cases the function is computed from the
    429     (m-1)-case using one of the recursions satisfied by the Legendre
    430     functions.
     1223    def _eval_special_values_(self,*args):
     1224        """
     1225        Special values known.
     1226        EXAMPLES:
     1227            sage: var('n')
     1228            n
     1229            sage: gen_laguerre(n,1,0)
     1230            binomial(n + 1, n)
     1231        """
     1232       
     1233        if (args[-1] == 0):
     1234            try:
     1235                return binomial(args[0]+args[1],args[0])
     1236            except TypeError:
     1237                pass
     1238
     1239        raise ValueError("Value not found")
     1240
     1241    def _eval_numpy_(self, *args):
     1242        """
     1243        EXAMPLES::
     1244        sage: import numpy
     1245        sage: z = numpy.array([1,2])
     1246        sage: z2 = numpy.array([[1,2],[1,2]])
     1247        sage: z3 = numpy.array([1,2,3.])
     1248        sage: gen_laguerre(1,1,z)
     1249        array([ 1.,  0.])
     1250        sage: gen_laguerre(1,1,z2)
     1251        array([[ 1.,  0.],
     1252               [ 1.,  0.]])
     1253        sage: gen_laguerre(1,1,z3)
     1254        array([ 1.,  0., -1.])
     1255        """
     1256
     1257        return eval_genlaguerre(args[0],args[1],args[-1])
     1258
    4311259   
    432     REFERENCE:
     1260    def _derivative_(self,*args,**kwds):
     1261        """return the derivative of gen_laguerre in
     1262        form of the gen. Laguerre Polynomial.
     1263        EXAMPLES::
     1264        sage: n = var('n')
     1265        sage: derivative(gen_laguerre(3,1,x),x)
     1266        -1/2*x^2 + 4*x - 6
     1267        sage: derivative(gen_laguerre(n,1,x),x)
     1268        -((n + 1)*gen_laguerre(n - 1, 1, x) - n*gen_laguerre(n, 1, x))/x
     1269        """
     1270        diff_param = kwds['diff_param']
     1271        if diff_param == 0:
     1272            raise NotImplementedError(
     1273"Derivative w.r.t. to the index is not supported, \
     1274yet, and perhaps never will be...")
     1275        else:
     1276            return (args[0]*gen_laguerre(args[0],args[1],args[-1])-(args[0]+args[1])*\
     1277                    gen_laguerre(args[0]-1,args[1],args[-1]))/args[-1]
    4331278
    434     - Gradshteyn and Ryzhik 8.706 page 1000.
    435    
    436     EXAMPLES::
    437    
    438         sage: P.<t> = QQ[]
    439         sage: gen_legendre_P(2, 0, t)
    440         3/2*t^2 - 1/2
    441         sage: gen_legendre_P(2, 0, t) == legendre_P(2, t)
    442         True
    443         sage: gen_legendre_P(3, 1, t)
    444         -3/2*sqrt(-t^2 + 1)*(5*t^2 - 1)
    445         sage: gen_legendre_P(4, 3, t)
    446         105*sqrt(-t^2 + 1)*(t^2 - 1)*t
    447         sage: gen_legendre_P(7, 3, I).expand()
    448         -16695*sqrt(2)
    449         sage: gen_legendre_P(4, 1, 2.5)
    450         -583.562373654533*I
    451     """
    452     from sage.functions.all import sqrt
    453     _init()
    454     if m.mod(2).is_zero() or m.is_one():
    455         return sage_eval(maxima.eval('assoc_legendre_p(%s,%s,x)'%(ZZ(n),ZZ(m))), locals={'x':x})
    456     else:
    457         return sqrt(1-x**2)*(((n-m+1)*x*gen_legendre_P(n,m-1,x)-(n+m-1)*gen_legendre_P(n-1,m-1,x))/(1-x**2))
     1279gen_laguerre = Func_gen_laguerre()
    4581280
    459 def gen_legendre_Q(n,m,x):
    460     """
    461     Returns the generalized (or associated) Legendre function of the
    462     second kind for integers `n>-1`, `m>-1`.
    463    
    464     Maxima restricts m = n. Hence the cases m n are computed using the
    465     same recursion used for gen_legendre_P(n,m,x) when m is odd and
    466     1.
    467    
    468     EXAMPLES::
    469    
    470         sage: P.<t> = QQ[]
    471         sage: gen_legendre_Q(2,0,t)
    472         3/4*t^2*log(-(t + 1)/(t - 1)) - 3/2*t - 1/4*log(-(t + 1)/(t - 1))
    473         sage: gen_legendre_Q(2,0,t) - legendre_Q(2, t)
    474         0
    475         sage: gen_legendre_Q(3,1,0.5)
    476         2.49185259170895
    477         sage: gen_legendre_Q(0, 1, x)
    478         -1/sqrt(-x^2 + 1)
    479         sage: gen_legendre_Q(2, 4, x).factor()
    480         48*x/((x - 1)^2*(x + 1)^2)
    481     """
    482     from sage.functions.all import sqrt
    483     if m <= n:
    484         _init()
    485         return sage_eval(maxima.eval('assoc_legendre_q(%s,%s,x)'%(ZZ(n),ZZ(m))), locals={'x':x})
    486     if m == n + 1 or n == 0:
    487         if m.mod(2).is_zero():
    488             denom = (1 - x**2)**(m/2)
    489         else:
    490             denom = sqrt(1 - x**2)*(1 - x**2)**((m-1)/2)   
    491         if m == n + 1:
    492             return (-1)**m*(m-1).factorial()*2**n/denom
    493         else:
    494             return (-1)**m*(m-1).factorial()*((x+1)**m - (x-1)**m)/(2*denom)
    495     else:
    496         return ((n-m+1)*x*gen_legendre_Q(n,m-1,x)-(n+m-1)*gen_legendre_Q(n-1,m-1,x))/sqrt(1-x**2)
    497 
    498 def hermite(n,x):
     1281class Func_hermite(OrthogonalPolynomial):
    4991282    """
    5001283    Returns the Hermite polynomial for integers `n > -1`.
    5011284   
     
    5031286
    5041287    - AS 22.5.40 and 22.5.41, page 779.
    5051288   
     1289    #sage: S.<y> = PolynomialRing(RR)
     1290    #   <- Here a bug of the base class BuiltinFunction occours!
     1291
    5061292    EXAMPLES::
    5071293   
    5081294        sage: x = PolynomialRing(QQ, 'x').gen()
     
    5121298        8*x^3 - 12*x
    5131299        sage: hermite(3,2)
    5141300        40
    515         sage: S.<y> = PolynomialRing(RR)
    516         sage: hermite(3,y)
     1301        sage: y = var('y')
     1302        sage: hermite(3,y).polynomial(PolynomialRing(RR, 'y'))
    5171303        8.00000000000000*y^3 - 12.0000000000000*y
    5181304        sage: R.<x,y> = QQ[]
    5191305        sage: hermite(3,y^2)
    5201306        8*y^6 - 12*y^2
    5211307        sage: w = var('w')
    5221308        sage: hermite(3,2*w)
    523         8*(8*w^2 - 3)*w
     1309        64*w^3 - 24*w
    5241310    """
    525     _init()
    526     return sage_eval(maxima.eval('hermite(%s,x)'%ZZ(n)), locals={'x':x})
     1311   
     1312    def __init__(self):
     1313        OrthogonalPolynomial.__init__(self,"hermite",nargs = 2,
     1314conversions =dict(maxima='hermite',mathematica='HermiteH'))
    5271315
    528 def jacobi_P(n,a,b,x):
     1316    def _clenshaw_method_(self,*args):
     1317        """
     1318        Clenshaw method for hermite polynomial (means use the recursion...)
     1319        See A.S. 227 (p. 782) for details for the recurions
     1320        For the symbolic evaluation, maxima seems to be quite fast.
     1321        The break even point between the recursion and Maxima is about
     1322        n = 25
     1323        """
     1324        k = args[0]
     1325        x = args[1]
     1326           
     1327        if k == 0:
     1328            return 1
     1329        elif k == 1:
     1330            return 2*x
     1331        else:
     1332            help1 = 1
     1333            help2 = 2*x
     1334            if is_Expression(x):
     1335                help3 = 0
     1336                for j in xrange(0,floor(k/2)+1):
     1337                    help3 = help3 + (-1)**j*(2*x)**(k-2*j)/factorial(j)/\
     1338                            factorial(k-2*j)
     1339                help3 = help3*factorial(k)
     1340            else:
     1341                for j in xrange(1,k):
     1342                    help3 = 2*x*help2 - 2*j*help1
     1343                    help1 = help2
     1344                    help2 = help3
     1345                   
     1346            return help3   
     1347   
     1348           
     1349    def _evalf_(self, *args,**kwds):
     1350        """
     1351        Evals hermite
     1352        numerically with mpmath.
     1353        EXAMPLES::
     1354            sage: hermite(3,2.).n(74)
     1355            40.0000000000000000000
     1356        """
     1357        try:
     1358            step_parent = kwds['parent']
     1359        except KeyError:
     1360            step_parent = parent(args[-1])
     1361
     1362        try:
     1363            precision = step_parent.prec()
     1364        except AttributeError:
     1365            precision = RR.prec()
     1366
     1367        from sage.libs.mpmath.all import call as mpcall
     1368        from sage.libs.mpmath.all import hermite as mphermite
     1369
     1370        return mpcall(mphermite,args[0],args[-1],prec = precision)
     1371   
     1372    def _eval_special_values_(self,*args):
     1373        """
     1374        Special values known. A.S. 22.4 (p.777)
     1375        EXAMPLES:
     1376       
     1377            sage: var('n')
     1378            n
     1379            sage: hermite(n,0)
     1380            ((-1)^n + 1)*(-1)^(1/2*n)*factorial(n)/gamma(1/2*n + 1)
     1381        """
     1382
     1383        if (args[-1] == 0):
     1384            return (1+(-1)**args[0])*(-1)**(args[0]/2)*\
     1385                   factorial(args[0])/gamma(args[0]/2+1)
     1386
     1387        raise ValueError("Value not found")
     1388
     1389    def _maxima_init_evaled_(self, *args):
     1390        """
     1391        Old maxima method.
     1392        """
     1393        n = args[0]
     1394        x = args[1]
     1395        return sage_eval(maxima.eval('hermite(%s,x)'%ZZ(n)), locals={'x':x})
     1396
     1397    def _eval_numpy_(self, *args):
     1398        """
     1399        EXAMPLES::
     1400        sage: import numpy
     1401        sage: z = numpy.array([1,2])
     1402        sage: z2 = numpy.array([[1,2],[1,2]])
     1403        sage: z3 = numpy.array([1,2,3.])
     1404        sage: hermite(1,z)
     1405        array([ 2.,  4.])
     1406        sage: hermite(1,z2)
     1407        array([[ 2.,  4.],
     1408               [ 2.,  4.]])
     1409        sage: hermite(1,z3)
     1410        array([ 2.,  4.,  6.])
     1411
     1412        """
     1413
     1414        #we have to convert the input to float first
     1415        return eval_hermite(nparray(args[0],dtype = float),\
     1416                            nparray(args[-1],dtype = float))
     1417
     1418
     1419    def _derivative_(self, *args, **kwds):
     1420        """
     1421        Returns the derivative of the hermite polynomial in form of the chebyshev
     1422        Polynomials of the first and second kind
     1423       
     1424        EXAMPLES::
     1425            sage: var('k')
     1426            k
     1427            sage: derivative(hermite(k,x),x)
     1428            2*k*hermite(k - 1, x)
     1429             
     1430        """
     1431        diff_param = kwds['diff_param']
     1432        if diff_param == 0:
     1433            raise NotImplementedError(
     1434"Derivative w.r.t. to the index is not supported, yet, and perhaps never will be...")
     1435        else:
     1436            return 2*args[0]*hermite(args[0]-1,args[1])
     1437 
     1438hermite = Func_hermite()       
     1439
     1440
     1441class Func_jacobi_P(OrthogonalPolynomial):
    5291442    r"""
    5301443    Returns the Jacobi polynomial `P_n^{(a,b)}(x)` for
    5311444    integers `n > -1` and a and b symbolic or `a > -1`
     
    5461459        sage: jacobi_P(2,1,2,1.2)        # random output of low order bits
    5471460        5.009999999999998
    5481461    """
    549     _init()
    550     return sage_eval(maxima.eval('jacobi_p(%s,%s,%s,x)'%(ZZ(n),a,b)), locals={'x':x})
    5511462
    552 def laguerre(n,x):
     1463    def __init__(self):
     1464        OrthogonalPolynomial.__init__(self,"jacobi_P",nargs = 4,
     1465conversions =dict(maxima='jacobi_p',mathematica='JacobiP'))
     1466   
     1467    def _clenshaw_method_(self,*args):
     1468        """
     1469        Clenshaw method for jacobi_P (means use the recursion,
     1470        or sum formula)
     1471        This is much faster for numerical evaluation than maxima!
     1472        See A.S. 227 (p. 782) for details for the recurions.
     1473        Warning: The clanshaw method for the Jacobi Polynomials
     1474        should only used for exact data types, when high orders are
     1475        used, due to weak instabilities of the recursion!
     1476        """
     1477        k = args[0]
     1478        x = args[-1]
     1479        alpha = args[1]
     1480        beta = args[2]
     1481
     1482        if k == 0:
     1483            return 1
     1484        elif k == 1:
     1485            return (alpha-beta + (alpha+beta+2)*x)/2
     1486        else:
     1487           
     1488            if is_Expression(x) or is_Expression(alpha) or is_Expression(beta):
     1489                #Here we use the sum formula of jacobi_P it seems this is rather
     1490                #optimal for use.
     1491                help1 = gamma(alpha+k+1)/factorial(k)/gamma(alpha+beta+k+1)
     1492                help2 = 0
     1493                for j in xrange(0,k+1):
     1494                    help2 = help2 + binomial(k,j)*gamma(alpha+beta+k+j+1)/\
     1495                            gamma(alpha+j+1)*((x-1)/2)**j
     1496                return help1*help2
     1497            else:
     1498                help1 = 1
     1499                help2 = (alpha-beta + (alpha+beta+2)*x)/2
     1500
     1501                for j in xrange(1,k):
     1502                    a1 = 2*(j+1)*(j+alpha+beta+1)*(2*j+alpha+beta)
     1503                    a2 = (2*j+alpha+beta+1)*(alpha**2-beta**2)
     1504                    a3 = (2*j+alpha+beta)*(2*j+alpha+beta+1)*(2*j+alpha+beta+2)
     1505                    a4 = 2*(j+alpha)*(j+beta)*(2*j+alpha+beta+2)
     1506                    help3 = (a2+a3*x)*help2 - a4*help1
     1507                    help3 = help3/a1
     1508                    help1 = help2
     1509                    help2 = help3
     1510               
     1511            return help3   
     1512
     1513    def _maxima_init_evaled_(self, *args):
     1514        """
     1515        Old maxima method.
     1516        """
     1517        n = args[0]
     1518        a = args[1]
     1519        b = args[2]
     1520        x = args[3]
     1521        return sage_eval(maxima.eval('jacobi_p(%s,%s,%s,x)'%(ZZ(n),a,b)),\
     1522                         locals={'x':x})
     1523
     1524    def _evalf_(self, *args,**kwds):
     1525        """
     1526        Evals jacobi_P
     1527        numerically with mpmath.
     1528        EXAMPLES::
     1529            sage: jacobi_P(10,2,3,3).n(75)
     1530            1.322776620000000000000e8
     1531        """
     1532        try:
     1533            step_parent = kwds['parent']
     1534        except KeyError:
     1535            step_parent = parent(args[-1])
     1536
     1537        try:
     1538            precision = step_parent.prec()
     1539        except AttributeError:
     1540            precision = RR.prec()
     1541
     1542        from sage.libs.mpmath.all import call as mpcall
     1543        from sage.libs.mpmath.all import jacobi as mpjacobi
     1544
     1545        return mpcall(mpjacobi,args[0],args[1],args[2],args[-1],\
     1546                           prec = precision)
     1547
     1548    def _eval_special_values_(self,*args):
     1549        """
     1550        Special values known. A.S. 22.4 (p.777)
     1551        EXAMPLES:
     1552           
     1553            sage: var('n k a y')
     1554            (n, k, a, y)
     1555            sage: jacobi_P(n,k,a,1)
     1556            binomial(k + n, n)
     1557        """
     1558        if args[-1] == 1:
     1559            return binomial(args[0]+args[1],args[0])
     1560       
     1561        if args[-1] == -1:
     1562            return (-1)**args[0]*binomial(args[0]+args[1],args[0])
     1563       
     1564        if args[1] == 0 and args[2] == 0:
     1565            return legendre_P(args[0],args[-1])
     1566
     1567        if args[1] == -0.5 and args[2] == -0.5:
     1568            try:
     1569                return binomial(2*args[0],args[0])*\
     1570                       chebyshev_T(args[0],args[-1])/4**args[0]
     1571            except TypeError:
     1572                pass
     1573
     1574        raise ValueError("Value not found")
     1575
     1576    def _eval_numpy_(self, *args):
     1577        """
     1578        EXAMPLES::
     1579        sage: import numpy
     1580        sage: z = numpy.array([1,2])
     1581        sage: z2 = numpy.array([[1,2],[1,2]])
     1582        sage: z3 = numpy.array([1,2,3.])
     1583        sage: jacobi_P(3,2,5,z)
     1584        array([  10.  ,  183.25])
     1585        sage: jacobi_P(3,2,z,0)
     1586        array([ -5.00000000e-01,  -1.38777878e-16])
     1587        sage: jacobi_P(3,2,1,z2)
     1588        array([[ 10. ,  90.5],
     1589               [ 10. ,  90.5]])
     1590        sage: jacobi_P(3,z3,1,0)
     1591        array([ 0. , -0.5, -1. ])
     1592        """
     1593
     1594        return eval_jacobi(args[0],args[1],args[2],args[-1])
     1595
     1596    def _derivative_(self, *args, **kwds):
     1597        """
     1598        Returns the derivative of jacobi_P in form of jacobi_polynomials
     1599       
     1600        EXAMPLES::
     1601            sage: var('k a b')
     1602            (k, a, b)
     1603            sage: derivative(jacobi_P(k,a,b,x),x)
     1604            -(2*(b + k)*(a + k)*jacobi_P(k - 1, a, b, x) - ((a + b + 2*k)*x - a + b)*k*jacobi_P(k, a, b, x))/((x^2 - 1)*(a + b + 2*k))
     1605            sage: derivative(jacobi_P(2,1,3,x),x)
     1606            14*x - 7/2
     1607            sage: derivative(jacobi_P(k,a,b,x),a)
     1608            Traceback (most recent call last):
     1609            ...
     1610            NotImplementedError: Derivative w.r.t. to the index is not supported, yet, and perhaps never will be...
     1611 
     1612        """
     1613        diff_param = kwds['diff_param']
     1614        if diff_param in [0,1,2]:
     1615            raise NotImplementedError(
     1616"Derivative w.r.t. to the index is not supported, yet, and perhaps never will be...")
     1617        else:
     1618            return (args[0]*(args[1]-args[2]-(2*args[0]+args[1]+args[2])*args[-1])*\
     1619                    jacobi_P(args[0],args[1],args[2],args[-1])+ 2*(args[0]+args[1])*\
     1620                    (args[0]+args[2])*jacobi_P(args[0]-1,args[1],args[2],args[-1]))/\
     1621                    (2*args[0]+args[1]+args[2])/(1-args[-1]**2)
     1622       
     1623       
     1624jacobi_P = Func_jacobi_P()       
     1625
     1626class Func_laguerre(OrthogonalPolynomial):
    5531627    """
    5541628    Returns the Laguerre polynomial for integers `n > -1`.
    5551629   
     
    5671641        sage: laguerre(2,2)
    5681642        -1
    5691643    """
    570     _init()
    571     return sage_eval(maxima.eval('laguerre(%s,x)'%ZZ(n)), locals={'x':x})
     1644    def __init__(self):
     1645        OrthogonalPolynomial.__init__(self,"laguerre",nargs = 2,
     1646conversions =dict(maxima='laguerre',mathematica='LaguerreL'))
     1647       
     1648    def _clenshaw_method_(self,*args):
     1649        """
     1650        Clenshaw method for laguerre polynomial (means use the recursion...)
     1651        This is much faster for numerical evaluation than maxima!
     1652        See A.S. 227 (p. 782) for details for the recurions.
     1653        Warning: The clanshaw method for the laguerre Polynomials
     1654        should only used for exact data types, when high orders are
     1655        used, due to weak instabilities of the recursion!
     1656        """
     1657        k = args[0]
     1658        x = args[-1]
     1659       
     1660        if k == 0:
     1661            return 1
     1662        elif k == 1:
     1663            return 1-x
     1664        else:
     1665            help1 = 1
     1666            help2 = 1-x
     1667            if is_Expression(x):
     1668                help3 = 0
     1669                for j in xrange(0,k+1):
     1670                    help3 = help3 + (-1)**j*binomial(k,k-j)/factorial(j)*x**j
     1671            else:
     1672                for j in xrange(1,k):
     1673                    help3 = -x*help2 + (2*j+1)*help2 - j*help1
     1674                    help3 = help3/(j+1)
     1675                    help1 = help2
     1676                    help2 = help3
     1677               
     1678            return help3   
    5721679
    573 def legendre_P(n,x):
     1680    def _maxima_init_evaled_(self, *args):
     1681        n = args[0]
     1682        x = args[1]
     1683        return sage_eval(maxima.eval('laguerre(%s,x)'%ZZ(n)), locals={'x':x})
     1684
     1685    def _evalf_(self, *args,**kwds):
     1686        """
     1687        Evals laguerre polynomial
     1688        numerically with mpmath.
     1689        EXAMPLES::
     1690            sage: laguerre(3,5.).n(53)
     1691            2.66666666666667
     1692        """
     1693        try:
     1694            step_parent = kwds['parent']
     1695        except KeyError:
     1696            step_parent = parent(args[-1])
     1697
     1698        try:
     1699            precision = step_parent.prec()
     1700        except AttributeError:
     1701            precision = RR.prec()
     1702
     1703        from sage.libs.mpmath.all import call as mpcall
     1704        from sage.libs.mpmath.all import laguerre as mplaguerre
     1705
     1706        return mpcall(mplaguerre,args[0],0,args[-1],prec = precision)
     1707   
     1708    def _eval_special_values_(self,*args):
     1709        """
     1710        Special values known.
     1711        EXAMPLES:
     1712           
     1713            sage: var('n')
     1714            n
     1715            sage: laguerre(n,0)
     1716            1
     1717        """
     1718       
     1719        if (args[-1] == 0):
     1720            try:
     1721                return 1
     1722            except TypeError:
     1723                pass
     1724
     1725        raise ValueError("Value not found")
     1726 
     1727
     1728    def _eval_numpy_(self, *args):
     1729        """
     1730        EXAMPLES::
     1731        sage: import numpy
     1732        sage: z = numpy.array([1,2])
     1733        sage: z2 = numpy.array([[1,2],[1,2]])
     1734        sage: z3 = numpy.array([1,2,3.])
     1735        sage: laguerre(1,z)
     1736        array([ 0., -1.])
     1737        sage: laguerre(1,z2)
     1738        array([[ 0., -1.],
     1739               [ 0., -1.]])
     1740        sage: laguerre(1,z3)
     1741        array([ 0., -1., -2.])
     1742
     1743        """
     1744       
     1745        return eval_laguerre(args[0],args[1])
     1746
     1747    def _derivative_(self,*args,**kwds):
     1748        """return the derivative of laguerre in
     1749        form of the Laguerre Polynomial.
     1750        EXAMPLES::
     1751        sage: n = var('n')
     1752        sage: derivative(laguerre(3,x),x)
     1753        -1/2*x^2 + 3*x - 3
     1754        sage: derivative(laguerre(n,x),x)
     1755        -(n*laguerre(n - 1, x) - n*laguerre(n, x))/x
     1756        """
     1757        diff_param = kwds['diff_param']
     1758        if diff_param == 0:
     1759            raise NotImplementedError(
     1760            "Derivative w.r.t. to the index is not supported, \
     1761             yet, and perhaps never will be...")
     1762        else:
     1763            return (args[0]*laguerre(args[0],args[-1])-args[0]*\
     1764                    laguerre(args[0]-1,args[-1]))/args[1]
     1765
     1766laguerre = Func_laguerre()
     1767
     1768class Func_legendre_P(OrthogonalPolynomial):
    5741769    """
    575     Returns the Legendre polynomial of the first kind for integers
    576     `n > -1`.
     1770    Returns the Legendre polynomial of the first kind..
    5771771   
    5781772    REFERENCE:
    5791773
     
    5881782        1.67750000000000
    5891783        sage: legendre_P(3, 1 + I)
    5901784        7/2*I - 13/2
    591         sage: legendre_P(3, MatrixSpace(ZZ, 2)([1, 2, -4, 7]))
    592         [-179  242]
    593         [-484  547]
     1785        sage: legendre_P(3, MatrixSpace(ZZ, 2)([1, 2, -4, 7]))
     1786        Traceback (most recent call last):
     1787        ...
     1788        TypeError: cannot coerce arguments: no canonical coercion from Full MatrixSpace of 2 by 2 dense matrices over Integer Ring to Symbolic Ring
    5941789        sage: legendre_P(3, GF(11)(5))
    5951790        8
    5961791    """
    597     _init()
    598     return sage_eval(maxima.eval('legendre_p(%s,x)'%ZZ(n)), locals={'x':x})
     1792    def __init__(self):
     1793        OrthogonalPolynomial.__init__(self,"legendre_P",nargs = 2,
     1794conversions =dict(maxima='legendre_p',mathematica='LegendreP'))
    5991795
    600 def legendre_Q(n,x):
     1796    def _clenshaw_method_(self,*args):
     1797        """
     1798        Clenshaw method for legendre_P (means use the recursion...)
     1799        This is much faster for numerical evaluation than maxima!
     1800        See A.S. 227 (p. 782) for details for the recurions.
     1801        Warning: The clanshaw method for the Legendre Polynomials
     1802        should only used for exact data types, when high orders are
     1803        used, due to weak instabilities of the recursion!
     1804        """
     1805        k = args[0]
     1806        x = args[-1]
     1807       
     1808        if k == 0:
     1809            return 1
     1810        elif k == 1:
     1811            return x
     1812        else:
     1813            help1 = 1
     1814            help2 = x
     1815            if is_Expression(x):
     1816                #raise NotImplementedError("Maxima is faster here...")
     1817                help1 = ZZ(2**k) #Workarround because of segmentation fault...
     1818                help3 = 0
     1819                for j in xrange(0,floor(k/2)+1):
     1820                    help3 = help3 + (-1)**j*x**(k-2*j)*binomial(k,j)*\
     1821                            binomial(2*(k-j),k)
     1822               
     1823                help3 = help3/help1
     1824            else:
     1825                for j in xrange(1,k):
     1826                    help3 = (2*j+1)*x*help2 - j*help1
     1827                    help3 = help3/(j+1)
     1828                    help1 = help2
     1829                    help2 = help3
     1830               
     1831            return help3   
     1832
     1833    def _maxima_init_evaled_(self, *args):
     1834        n = args[0]
     1835        x = args[1]
     1836        return sage_eval(maxima.eval('legendre_p(%s,x)'%ZZ(n)),\
     1837                         locals={'x':x})
     1838
     1839    def _evalf_(self, *args,**kwds):
     1840        """
     1841        Evals legendre_P
     1842        numerically with mpmath.
     1843        EXAMPLES::
     1844            sage: legendre_P(10,3).n(75)
     1845            8.097453000000000000000e6
     1846        """
     1847        try:
     1848            step_parent = kwds['parent']
     1849        except KeyError:
     1850            step_parent = parent(args[-1])
     1851
     1852        try:
     1853            precision = step_parent.prec()
     1854        except AttributeError:
     1855            precision = RR.prec()
     1856
     1857        from sage.libs.mpmath.all import call as mpcall
     1858        from sage.libs.mpmath.all import legendre as mplegendre
     1859
     1860        return mpcall(
     1861            mplegendre,args[0],args[-1],prec = precision)
     1862   
     1863
     1864    def _eval_special_values_(self,*args):
     1865        """
     1866        Special values known.
     1867        EXAMPLES:
     1868           
     1869            sage: var('n')
     1870            n
     1871            sage: legendre_P(n,1)
     1872            1
     1873            sage: legendre_P(n,-1)
     1874            (-1)^n
     1875        """
     1876        if args[-1] == 1:
     1877            return 1
     1878       
     1879        if args[-1] == -1:
     1880            return (-1)**args[0]
     1881       
     1882        if (args[-1] == 0):
     1883            try:
     1884                return (1+(-1)**args[0])/2*binomial(args[0],args[0]/2)/4**(args[0]/2)
     1885            except TypeError:
     1886                pass
     1887
     1888        raise ValueError("Value not found")
     1889
     1890    def _eval_numpy_(self, *args):
     1891        """
     1892        EXAMPLES::
     1893        sage: import numpy
     1894        sage: z = numpy.array([1,2])
     1895        sage: z2 = numpy.array([[1,2],[1,2]])
     1896        sage: z3 = numpy.array([1,2,3.])
     1897        sage: legendre_P(1,z)
     1898        array([ 1.,  2.])
     1899        sage: legendre_P(1,z2)
     1900        array([[ 1.,  2.],
     1901               [ 1.,  2.]])
     1902        sage: legendre_P(1,z3)
     1903        array([ 1.,  2.,  3.])
     1904        sage: legendre_P(z3,3)
     1905        array([  3.,  13.,  63.])
     1906
     1907        """
     1908
     1909        return eval_legendre(args[0],args[1])
     1910
     1911    def _derivative_(self,*args,**kwds):
     1912        """return the derivative of legendre_P in
     1913        form of the Legendre Polynomial.
     1914        EXAMPLES::
     1915        sage: n = var('n')
     1916        sage: derivative(legendre_P(n,x),x)
     1917        (n*x*legendre_P(n, x) - n*legendre_P(n - 1, x))/(x^2 - 1)
     1918        sage: derivative(legendre_P(3,x),x)
     1919        15/2*x^2 - 3/2
     1920        """
     1921        diff_param = kwds['diff_param']
     1922        if diff_param == 0:
     1923            raise NotImplementedError(
     1924"Derivative w.r.t. to the index is not supported, yet, and perhaps never will be...")
     1925        else:
     1926            return (args[0]*legendre_P(args[0]-1,args[-1])-args[0]*args[-1]*\
     1927                    legendre_P(args[0],args[-1]))/(1-args[-1]**2)
     1928       
     1929
     1930legendre_P = Func_legendre_P()     
     1931
     1932class Func_legendre_Q(OrthogonalPolynomial):
    6011933    """
    602     Returns the Legendre function of the second kind for integers
    603     `n>-1`.
     1934    Return the chebyshev function of the second kind.
    6041935   
    605     Computed using Maxima.
     1936    REFERENCE:
     1937   
     1938    - A.S. 8 (p. 332)
    6061939   
    6071940    EXAMPLES::
    6081941   
     
    6141947        sage: legendre_Q(4, 2)
    6151948        443/16*I*pi + 443/16*log(3) - 365/12
    6161949        sage: legendre_Q(4, 2.0)
    617         0.00116107583162324 + 86.9828465962674*I
     1950        0.00116107583162041 + 86.9828465962674*I
     1951
    6181952    """
    619     _init()
    620     return sage_eval(maxima.eval('legendre_q(%s,x)'%ZZ(n)), locals={'x':x})
     1953    def __init__(self):
     1954        OrthogonalPolynomial.__init__(self,"legendre_Q",nargs = 2,
     1955conversions =dict(maxima='legendre_q',mathematica='LegendreQ'))
    6211956
    622 def ultraspherical(n,a,x):
    623     """
    624     Returns the ultraspherical (or Gegenbauer) polynomial for integers
    625     `n > -1`.
     1957    def _maxima_init_evaled_(self, *args):
     1958        """
     1959        Maxima seems just fine for legendre Q. So we use it here!
     1960        """
     1961        n = args[0]
     1962        x = args[1]
     1963        return sage_eval(maxima.eval('legendre_q(%s,x)'%ZZ(n)),\
     1964                         locals={'x':x})
     1965
     1966    def _clenshaw_method_(self,*args):
     1967        """
     1968        Clenshaw method for legendre_q (means use the recursion...)
     1969        This is much faster for numerical evaluation than maxima!
     1970        See A.S. 8.5.3 (p. 334) for details for the recurions.
     1971        Warning: The clanshaw method for the Legendre fUNCTIONS
     1972        should only used for exact data types, when high orders are
     1973        used, due to weak instabilities of the recursion!
     1974        """
     1975        raise NotImplementedError("Function not ready yet...")
     1976       
     1977        k = args[0]
     1978        x = args[-1]
     1979       
     1980        if k == 0:
     1981            return ln((1+x)/(1-x))/2
     1982        elif k == 1:
     1983            return x/2*ln((1+x)/(1-x))-1
     1984        else:
     1985            if is_Expression(x):
     1986                raise NotImplementedError("Maxima works fine here!")
     1987                #it seems that the old method just works fine here...
     1988                #raise NotImplementedError("clenshaw does not work well...")
     1989            else:
     1990                help1 = ln((1+x)/(1-x))/2
     1991                help2 = x/2*ln((1+x)/(1-x))-1
     1992
     1993                for j in xrange(1,k):
     1994                    help3 = (2*j+1)*x*help2 - j*help1
     1995                    help3 = help3/(j+1)
     1996                    help1 = help2
     1997                    help2 = help3
     1998               
     1999            return help3   
     2000
     2001    def _eval_special_values_(self,*args):
     2002        """
     2003        Special values known.
     2004        EXAMPLES:
     2005           
     2006            sage: var('n')
     2007            n
     2008            sage: legendre_Q(n,0)
     2009            -1/2*sqrt(pi)*gamma(1/2*n + 1/2)*sin(1/2*pi*n)/gamma(1/2*n + 1)
     2010        """
     2011        if args[-1] == 1:
     2012            return NaN
     2013       
     2014        if args[-1] == -1:
     2015            return NaN
     2016       
     2017        if (args[-1] == 0):
     2018            if is_Expression(args[0]):
     2019                try:
     2020                    return -(sqrt(SR.pi()))/2*sin(SR.pi()/2*args[0])*\
     2021                           gamma((args[0]+1)/2)/gamma(args[0]/2 + 1)
     2022                except TypeError:
     2023                    pass
     2024            else:
     2025                return -(sqrt(math.pi))/2*sin(math.pi/2*args[0])*\
     2026                       gamma((args[0]+1)/2)/gamma(args[0]/2. + 1)
     2027
     2028        raise ValueError("Value not found")
     2029
     2030    def _evalf_(self, *args,**kwds):
     2031        """
     2032        Evals legendre_Q
     2033        numerically with mpmath.
     2034        EXAMPLES::
     2035            sage: legendre_Q(10,3).n(75)
     2036            2.079454941572578263731e-9 + 1.271944942879431601408e7*I
     2037            sage: legendre_Q(3, 0.5).n(53)
     2038            -0.198654771479482
     2039        """
     2040
     2041        try:
     2042            step_parent = kwds['parent']
     2043        except KeyError:
     2044            step_parent = parent(args[-1])
     2045
     2046        try:
     2047            precision = step_parent.prec()
     2048        except AttributeError:
     2049            precision = RR.prec()
     2050       
     2051        from sage.libs.mpmath.all import call as mpcall
     2052        from sage.libs.mpmath.all import legenq as mplegenq
     2053
     2054        return conjugate(
     2055            mpcall(mplegenq,args[0],0,args[-1],prec = precision))
     2056        #it seems that mpmath uses here a different branch of the logarithm
     2057
     2058    def _eval_numpy_(self, *args):
     2059        #TODO: numpy_eval with help of the a newer scipy version!!!!
     2060        #Reason scipy supports stable and fast numerical evaluation
     2061        #of ortho polys, but not the current version!
     2062        #Now this only evaluates the array pointwise, and only the first one...
     2063        #if isinstance(arg[0], numpy.ndarray). This is a hack to provide compability!
     2064        """
     2065        EXAMPLES::
     2066            sage: import numpy
     2067            sage: z = numpy.array([1,2])
     2068            sage: z2 = numpy.array([[1,2],[1,2]])
     2069            sage: z3 = numpy.array([1,2,3.])
     2070            sage: legendre_Q(1,z/5.)
     2071            array([-0.95945349, -0.83054043])
     2072            sage: legendre_Q(1,z2/5.)
     2073            array([[-0.95945349, -0.83054043],
     2074                   [-0.95945349, -0.83054043]])
     2075            sage: legendre_Q(1,z3/5.)
     2076            array([-0.95945349, -0.83054043, -0.58411169])
     2077        """
     2078        #The imports are made here because this isn't optimal anyway
     2079        from numpy import ndarray, zeros
     2080
     2081        if isinstance(args[0],ndarray):
     2082            raise NotImplementedError(
     2083           "Support for numpy array in first argument(s) is not supported yet!")
     2084
     2085        result = zeros(args[-1].shape).tolist()
     2086        if isinstance(args[-1][0],int):
     2087            for k in xrange(len(args[-1])):
     2088                result[k] = legendre_Q(args[0],ZZ(args[-1][k]))
     2089
     2090        if isinstance(args[-1][0],float):
     2091            for k in xrange(len(args[-1])):
     2092                result[k] = legendre_Q(args[0],RR(args[-1][k]))
     2093
     2094        if isinstance(args[-1][0],ndarray):
     2095            for k in xrange(len(args[-1])):
     2096                result[k] = legendre_Q(args[0],args[-1][k])
     2097
     2098        return nparray(result)
     2099
     2100    def _derivative_(self,*args,**kwds):
     2101        """return the derivative of legendre_Q in
     2102        form of the Legendre Function.
     2103        EXAMPLES::
     2104        n = var('n')
     2105        derivative(legendre_Q(n,x),x)
     2106        (n*x*legendre_Q(n, x) - n*legendre_Q(n - 1, x))/(x^2 - 1)
     2107        sage: derivative(legendre_Q(3,x),x)
     2108        5/4*(x - 1)*(1/(x - 1) - (x + 1)/(x - 1)^2)*x^3/(x + 1) + 15/4*x^2*log(-(x + 1)/(x - 1)) - 3/4*(x - 1)*(1/(x - 1) - (x + 1)/(x - 1)^2)*x/(x + 1) - 5*x - 3/4*log(-(x + 1)/(x - 1))
     2109        """
     2110        diff_param = kwds['diff_param']
     2111        if diff_param == 0:
     2112            raise NotImplementedError(
     2113"Derivative w.r.t. to the index is not supported, yet, and perhaps never will be...")
     2114        else:
     2115            return (args[0]*args[-1]*legendre_Q(args[0],args[-1])-args[0]*\
     2116                    legendre_Q(args[0]-1,args[-1]))/(args[-1]**2-1)
     2117
     2118legendre_Q = Func_legendre_Q()
     2119
     2120
     2121class Func_gen_legendre_P(OrthogonalPolynomial):
     2122
     2123    def __init__(self):
     2124        OrthogonalPolynomial.__init__(self,"gen_legendre_P",nargs = 3,
     2125conversions =dict(maxima='assoc_legendre_p',mathematica='LegendreP'))
     2126
     2127    def _evalf_(self, *args,**kwds):
     2128        """
     2129        Evals gen_legendre_P
     2130        numerically with mpmath.
     2131        EXAMPLES::
     2132            sage: gen_legendre_P(10,2,3).n(75)
     2133            -7.194963600000000000000e8
     2134        """
     2135
     2136        try:
     2137            step_parent = kwds['parent']
     2138        except KeyError:
     2139            step_parent = parent(args[-1])
     2140
     2141        try:
     2142            precision = step_parent.prec()
     2143        except AttributeError:
     2144            precision = RR.prec()
     2145
     2146        from sage.libs.mpmath.all import call as mpcall
     2147        from sage.libs.mpmath.all import legenp as mplegenp
     2148
     2149        return mpcall(
     2150            mplegenp,args[0],args[1],args[-1],prec = precision)
     2151
     2152    def _eval_special_values_(self,*args):
     2153        """
     2154        Special values known.
     2155        EXAMPLES:
     2156           
     2157            sage: n, m = var('n m')
     2158            sage: gen_legendre_P(n,m,0)
     2159            2^m*gamma(1/2*m + 1/2*n + 1/2)*cos(1/2*(m + n)*pi)/(sqrt(pi)*gamma(-1/2*m + 1/2*n + 1))
     2160        """
     2161           
     2162        if args[1] == 0:
     2163            return legendre_P(args[0],args[-1])
     2164
     2165        if (args[-1] == 0):
     2166            if is_Expression(args[0]):
     2167                try:
     2168                    return cos(SR.pi()/2*(args[0]+args[1]))/(sqrt(SR.pi()))*\
     2169                           gamma((args[0]+args[1]+1)/2)/\
     2170                           gamma((args[0]-args[1])/2 + 1)*2**(args[1])
     2171                except TypeError:
     2172                    pass
     2173            else:
     2174                return cos(math.pi/2*(args[0]+args[1]))/(sqrt(math.pi))*\
     2175                       gamma((args[0]+args[1]+1)/2)/\
     2176                       gamma((args[0]-args[1])/2. + 1)*2**args[1]
     2177
     2178        raise ValueError("Value not found")
     2179
     2180    def _maxima_init_evaled_(self, *args):
     2181        n = args[0]
     2182        m = args[1]
     2183        x = args[2]
     2184        if is_Expression(n) or is_Expression(m):
     2185            return None
     2186       
     2187        from sage.functions.all import sqrt
    6262188   
    627     Computed using Maxima.
     2189        if m.mod(2).is_zero() or m.is_one():
     2190            return sage_eval(maxima.eval('assoc_legendre_p(%s,%s,x)'\
     2191                                         %(ZZ(n),ZZ(m))), locals={'x':x})
     2192        else:
     2193            return sqrt(1-x**2)*(((n-m+1)*x*gen_legendre_P(n,m-1,x)-\
     2194                                  (n+m-1)*gen_legendre_P(n-1,m-1,x))/(1-x**2))
     2195
     2196
     2197    def _eval_numpy_(self, *args):
     2198        """
     2199        EXAMPLES::
     2200        sage: import numpy
     2201        sage: z = numpy.array([1,2])
     2202        sage: z2 = numpy.array([[1,2],[1,2]])
     2203        sage: z3 = numpy.array([1,2,3.])
     2204        sage: gen_legendre_P(1,1,z/5.)
     2205        array([-0.9797959 , -0.91651514])
     2206        sage: gen_legendre_P(1,1,z2/5.)
     2207        array([[-0.9797959 , -0.91651514],
     2208               [-0.9797959 , -0.91651514]])
     2209        sage: gen_legendre_P(1,1,z3/5.)
     2210        array([-0.9797959 , -0.91651514, -0.8       ])
     2211        """
     2212        #in scipy arguments 0,1 are in reverse order
     2213        return lpmv(args[1],args[0],args[-1])
     2214       
     2215     
     2216    def _derivative_(self,*args,**kwds):
     2217        """return the derivative of gen_legendre_P in
     2218        form of the Legendre Function.
     2219        EXAMPLES::
     2220        sage: n,m = var('n m')
     2221        sage: derivative(gen_legendre_P(n,m,x),x)
     2222        (n*x*gen_legendre_P(n, m, x) - (m + n)*gen_legendre_P(n - 1, m, x))/(x^2 - 1)
     2223        sage: derivative(gen_legendre_P(3,1,x),x)
     2224        -15*sqrt(-x^2 + 1)*x + 3/2*(5*(x - 1)^2 + 10*x - 6)*x/sqrt(-x^2 + 1)
     2225        """
     2226        diff_param = kwds['diff_param']
     2227        if diff_param in [0,1]:
     2228            raise NotImplementedError(
     2229            "Derivative w.r.t. to the index is not supported,\
     2230            yet, and perhaps never will be...")
     2231        else:
     2232            return (args[0]*args[-1]*gen_legendre_P(args[0],args[1],args[-1])\
     2233                    -(args[0]+args[1])*\
     2234                    gen_legendre_P(args[0]-1,args[1],args[-1]))/(args[-1]**2-1)
     2235
     2236gen_legendre_P = Func_gen_legendre_P()
     2237
     2238class Func_gen_legendre_Q(OrthogonalPolynomial):
     2239
     2240    def __init__(self):
     2241        OrthogonalPolynomial.__init__(self,"gen_legendre_Q",nargs = 3,
     2242conversions =dict(maxima='assoc_legendre_q',mathematica='LegendreQ'))
     2243
     2244    def _evalf_(self, *args,**kwds):
     2245        """
     2246        Evals gen_legendre_Q
     2247        numerically with mpmath.
     2248        EXAMPLES::
     2249            sage: gen_legendre_Q(10,2,3).n(75)
     2250            -2.773909528741569374688e-7 - 1.130182239430298584113e9*I
     2251        """
     2252
     2253        try:
     2254            step_parent = kwds['parent']
     2255        except KeyError:
     2256            step_parent = parent(args[-1])
     2257
     2258        try:
     2259            precision = step_parent.prec()
     2260        except AttributeError:
     2261            precision = RR.prec()
     2262
     2263        from sage.libs.mpmath.all import call as mpcall
     2264        from sage.libs.mpmath.all import legenq as mplegenq
     2265
     2266        return mpcall(
     2267            mplegenq,args[0],args[1],args[-1],prec = precision)
     2268
     2269    def _eval_special_values_(self,*args):
     2270        """
     2271        Special values known.
     2272        EXAMPLES:
     2273           
     2274            sage: n, m = var('n m')
     2275            sage: gen_legendre_Q(n,m,0)
     2276            -sqrt(pi)*2^(m - 1)*gamma(1/2*m + 1/2*n + 1/2)*sin(1/2*(m + n)*pi)/gamma(-1/2*m + 1/2*n + 1)
     2277        """
     2278           
     2279        if args[1] == 0:
     2280            return legendre_Q(args[0],args[-1])
     2281
     2282        if (args[-1] == 0):
     2283            if is_Expression(args[0]):
     2284                try:
     2285                    return -(sqrt(SR.pi()))*sin(SR.pi()/2*(args[0]+args[1]))\
     2286                           *gamma((args[0]+args[1]+1)/2)/\
     2287                           gamma((args[0]-args[1])/2 + 1)*2**(args[1]-1)
     2288                except TypeError:
     2289                    pass
     2290            else:
     2291                return -(sqrt(math.pi))/2*sin(math.pi/2*(args[0]+args[1]))\
     2292                       *gamma((args[0]+args[1]+1)/2)/\
     2293                       gamma((args[0]-args[1])/2. + 1)*2**args[1]
     2294
     2295        raise ValueError("Value not found")
     2296
     2297    def _maxima_init_evaled_(self, *args):
     2298        n = args[0]
     2299        m = args[1]
     2300        x = args[2]
     2301        if is_Expression(n) or is_Expression(m):
     2302            return None
     2303       
     2304        from sage.functions.all import sqrt
    6282305   
    629     REFERENCE:
     2306        if m <= n:
     2307            return sage_eval(maxima.eval('assoc_legendre_q(%s,%s,x)'\
     2308                                         %(ZZ(n),ZZ(m))), locals={'x':x})
     2309        if m == n + 1 or n == 0:
     2310            if m.mod(2).is_zero():
     2311                denom = (1 - x**2)**(m/2)
     2312            else:
     2313                denom = sqrt(1 - x**2)*(1 - x**2)**((m-1)/2)   
     2314                if m == n + 1:
     2315                    return (-1)**m*(m-1).factorial()*2**n/denom
     2316                else:
     2317                    return (-1)**m*(m-1).factorial()*((x+1)**m - (x-1)**m)/(2*denom)
     2318        else:
     2319            return ((n-m+1)*x*gen_legendre_Q(n,m-1,x)-\
     2320                    (n+m-1)*gen_legendre_Q(n-1,m-1,x))/sqrt(1-x**2)
    6302321
    631     - AS 22.5.27
    632    
    633     EXAMPLES::
    634    
    635         sage: x = PolynomialRing(QQ, 'x').gen()
    636         sage: ultraspherical(2,3/2,x)
    637         15/2*x^2 - 3/2
    638         sage: ultraspherical(2,1/2,x)
    639         3/2*x^2 - 1/2
    640         sage: ultraspherical(1,1,x)
    641         2*x     
    642         sage: t = PolynomialRing(RationalField(),"t").gen()
    643         sage: gegenbauer(3,2,t)
    644         32*t^3 - 12*t
    645     """
    646     _init()
    647     return sage_eval(maxima.eval('ultraspherical(%s,%s,x)'%(ZZ(n),a)), locals={'x':x})
    6482322
    649 gegenbauer = ultraspherical
     2323    def _eval_numpy_(self, *args):
     2324        #TODO: numpy_eval with help of a new scipy version!!!!
     2325        #Now this only evaluates the array pointwise, and only the first one...
     2326        #if isinstance(arg[0], numpy.ndarray).
     2327        #This is a hack to provide compability with older releases of sage!
     2328        """
     2329        EXAMPLES::
     2330        sage: import numpy
     2331        sage: z = numpy.array([1,2])
     2332        sage: z2 = numpy.array([[1,2],[1,2]])
     2333        sage: z3 = numpy.array([1,2,3.])
     2334        sage: gen_legendre_Q(1,1,z/5.)
     2335        array([-0.40276067, -0.82471644])
     2336        sage: gen_legendre_Q(1,1,z2/5.)
     2337        array([[-0.40276067, -0.82471644],
     2338               [-0.40276067, -0.82471644]])
     2339        sage: gen_legendre_Q(1,1,z3/5.)
     2340        array([-0.40276067, -0.82471644, -1.30451774])
     2341        """
     2342        #imports are included here, because this isn't optimal anyway
     2343        from numpy import ndarray, zeros
     2344
     2345        if isinstance(args[0],ndarray) or isinstance(args[1],ndarray):
     2346            raise NotImplementedError(
     2347           "Support for numpy array in first argument(s) is not supported yet!")
     2348       
     2349        result = zeros(args[-1].shape).tolist()
     2350        if isinstance(args[-1][0],int):
     2351            for k in xrange(len(args[-1])):
     2352                result[k] = gen_legendre_Q(args[0],args[1],ZZ(args[-1][k]))
     2353
     2354        if isinstance(args[-1][0],float):
     2355            for k in xrange(len(args[-1])):
     2356                result[k] = gen_legendre_Q(args[0],args[1],RR(args[-1][k]))
     2357
     2358        if isinstance(args[-1][0],ndarray):
     2359            for k in xrange(len(args[-1])):
     2360                result[k] = gen_legendre_Q(args[0],args[1],args[-1][k])
     2361
     2362        return nparray(result)
     2363     
     2364    def _derivative_(self,*args,**kwds):
     2365        """return the derivative of gen_legendre_Q in
     2366        form of the Legendre Function.
     2367        EXAMPLES::
     2368        sage: n,m = var('n m')
     2369        sage: derivative(gen_legendre_Q(n,m,x),x)
     2370        (n*x*gen_legendre_Q(n, m, x) - (m + n)*gen_legendre_Q(n - 1, m, x))/(x^2 - 1)
     2371        sage: derivative(gen_legendre_Q(0,1,x),x)
     2372        -x/(-x^2 + 1)^(3/2)
     2373        """
     2374        diff_param = kwds['diff_param']
     2375        if diff_param in [0,1]:
     2376            raise NotImplementedError("Derivative w.r.t. to the index \
     2377            is not supported, yet, and perhaps never will be...")
     2378        else:
     2379            return (args[0]*args[-1]*gen_legendre_Q(args[0],args[1],args[-1])\
     2380                    -(args[0]+args[1])*\
     2381                    gen_legendre_Q(args[0]-1,args[1],args[-1]))/(args[-1]**2-1)
     2382
     2383
     2384gen_legendre_Q = Func_gen_legendre_Q()
     2385
     2386
     2387
     2388
     2389
     2390
     2391
  • sage/symbolic/pynac.pyx

    diff -r 7ae7b5175304 -r 5fac168ed1fc sage/symbolic/pynac.pyx
    a b  
    379379    Test custom func::
    380380
    381381        sage: def my_print(self, *args): return "my args are: " + ', '.join(map(repr, args))
     382
    382383        sage: foo = function('foo', nargs=2, print_latex_func=my_print)
    383         sage: for i in range(get_ginac_serial(), get_ginac_serial()+50):
     384        sage: for i in range(get_ginac_serial(), get_ginac_serial()+100):
    384385        ...     if get_sfunction_from_serial(i) == foo: break
    385386
    386387        sage: get_sfunction_from_serial(i) == foo
     
    460461        (x, y, z)
    461462        sage: from sage.symbolic.function import get_sfunction_from_serial
    462463        sage: foo = function('foo', nargs=2)
    463         sage: for i in range(get_ginac_serial(), get_ginac_serial()+50):
     464        sage: for i in range(get_ginac_serial(), get_ginac_serial()+100):
    464465        ...     if get_sfunction_from_serial(i) == foo: break
    465466
    466467        sage: get_sfunction_from_serial(i) == foo
     
    472473
    473474        sage: def my_print(self, *args): return "func_with_args(" + ', '.join(map(repr, args)) +')'
    474475        sage: foo = function('foo', nargs=2, print_func=my_print)
    475         sage: for i in range(get_ginac_serial(), get_ginac_serial()+50):
     476        sage: for i in range(get_ginac_serial(), get_ginac_serial()+100):
    476477        ...     if get_sfunction_from_serial(i) == foo: break
    477478
    478479        sage: get_sfunction_from_serial(i) == foo
     
    522523    Test latex_name::
    523524
    524525        sage: foo = function('foo', nargs=2, latex_name=r'\mathrm{bar}')
    525         sage: for i in range(get_ginac_serial(), get_ginac_serial()+50):
     526        sage: for i in range(get_ginac_serial(), get_ginac_serial()+100):
    526527        ...     if get_sfunction_from_serial(i) == foo: break
    527528
    528529        sage: get_sfunction_from_serial(i) == foo
     
    534535
    535536        sage: def my_print(self, *args): return "func_with_args(" + ', '.join(map(repr, args)) +')'
    536537        sage: foo = function('foo', nargs=2, print_latex_func=my_print)
    537         sage: for i in range(get_ginac_serial(), get_ginac_serial()+50):
     538        sage: for i in range(get_ginac_serial(), get_ginac_serial()+100):
    538539        ...     if get_sfunction_from_serial(i) == foo: break
    539540
    540541        sage: get_sfunction_from_serial(i) == foo
  • sage/symbolic/random_tests.py

    diff -r 7ae7b5175304 -r 5fac168ed1fc sage/symbolic/random_tests.py
    a b  
    1515
    1616        sage: from sage.symbolic.random_tests import _mk_full_functions
    1717        sage: [f for (one,f,arity) in _mk_full_functions()]
    18         [Ei, abs, arccos, arccosh, arccot, arccoth, arccsc, arccsch,
    19         arcsec, arcsech, arcsin, arcsinh, arctan, arctan2, arctanh,
    20         binomial, ceil, conjugate, cos, cosh, cot, coth, csc, csch,
    21         dickman_rho, dilog, dirac_delta, elliptic_e, elliptic_ec,
    22         elliptic_eu, elliptic_f, elliptic_kc, elliptic_pi, erf, exp,
    23         factorial, floor, heaviside, imag_part, integrate,
    24         kronecker_delta, log, polylog, real_part, sec, sech, sgn, sin,
    25         sinh, tan, tanh, unit_step, zeta, zetaderiv]
     18        [Ei, abs, arccos, arccosh, arccot, arccoth, arccsc, arccsch, arcsec,
     19        arcsech, arcsin, arcsinh, arctan, arctan2, arctanh, binomial, ceil,
     20        chebyshev_T, chebyshev_U, conjugate, cos, cosh, cot, coth, csc, csch, dickman_rho,
     21        dilog, dirac_delta, elliptic_e, elliptic_ec, elliptic_eu, elliptic_f, elliptic_kc,
     22        elliptic_pi, erf, exp, factorial, floor, gegenbauer, gen_laguerre, gen_legendre_P,
     23        gen_legendre_Q, heaviside, hermite, imag_part, integrate, jacobi_P, kronecker_delta,
     24        laguerre, legendre_P, legendre_Q, log, polylog, real_part, sec, sech, sgn, sin, sinh, tan,
     25        tanh, unit_step, zeta, zetaderiv]
     26
    2627
    2728    Note that this doctest will fail whenever a Pynac function is added or
    2829    removed.  In that case, it is very likely that the doctests for
     
    234235
    235236        sage: from sage.symbolic.random_tests import *
    236237        sage: set_random_seed(2)
    237         sage: random_expr(50, nvars=3, coeff_generator=CDF.random_element)
    238         (euler_gamma - v3^(-e) + (v2 - factorial(-e/v2))^(((2.85879036573 - 1.18163393202*I)*v2 + (2.85879036573 - 1.18163393202*I)*v3)*pi - 0.247786879678 + 0.931826724898*I)*arccsc((0.891138386848 - 0.0936820840629*I)/v1) - (0.553423153995 - 0.5481180572*I)*v3 + 0.149683576515 - 0.155746451854*I)*v1 + arccsch(pi + e)*elliptic_f(khinchin*v2, 1.4656989704 + 0.863754357069*I)
     238        sage: random_expr(60, nvars=3, coeff_generator=CDF.random_element)
     239        ceil(-(v1 + 0.149683576515 - 0.155746451854*I)*((catalan + v1)*v3 - csc(v3 - 0.23544741659 + 0.838316605873*I))/v1)/(pi*v1 + (1.38519897979 - 1.09891497125*I)/elliptic_ec(jacobi_P(-1/(v3 - arccos(v3 + 0.423315908054 - 0.313381569707*I)), -(0.33181415602 + 0.938333932603*I)*v2, -0.811204447488 + 0.547720758646*I, 1.14668794934 + 0.40478610032*I)) + erf(-v1 + v3 + elliptic_f(-0.418259557401 + 0.0431731633732*I, 1/euler_gamma) + 0.308808895751 - 0.288302424226*I) + 0.390821882125)
    239240        sage: random_expr(5, verbose=True)
    240         About to apply dirac_delta to [1]
    241         About to apply arccsch to [0]
    242         About to apply <built-in function add> to [0, arccsch(0)]
    243         arccsch(0)
     241        About to apply <built-in function div> to [v1, e]
     242        About to apply <built-in function inv> to [v1*e^(-1)]
     243        About to apply <built-in function neg> to [e/v1]
     244        -e/v1
     245
    244246    """
    245247    vars = [(1.0, sage.calculus.calculus.var('v%d' % (n+1))) for n in range(nvars)]
    246248    if ncoeffs is None: