Ticket #9706: trac_9706_ortho_polys.patch

File trac_9706_ortho_polys.patch, 79.6 KB (added by maldun, 9 years ago)

Patch for latest version with some code cleanup (no program changes)

  • sage/functions/orthogonal_polys.py

    # HG changeset patch
    # User maldun <domors@gmx.net>
    # Date 1283013861 -7200
    # Node ID 2a39ac68c1dbf8977f43c7cf27de0061b5cb1ca5
    # Parent  5b338f2e484f2065d3d30d47bc204d6e9ed13d12
    Trac 9706: New version of orthogonal polynomials. There is a new base class for
    generation of orthogonal polynomials now, derived from BuiltinFunction to provide
    better symbolic handling. Except for the legendre functions legendre_Q, gen_legendre_P,
    and gen_legendre_Q which are using maxima for symbolic evaluation, all polys now
    use explicit formulas, because this is faster for evaluation, and recursion
    for numerical evaluation of exact data types, because this is faster.
    For numerical evaluation mpmath is used, because it provides arbitrary precision
    
    diff -r 5b338f2e484f -r 2a39ac68c1db 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. Also the orthogonal polys
     16are handled as symbolic expressions, if the input is not known
    1117
    1218
    1319-  The Chebyshev polynomial of the first kind arises as a solution
     
    289295AUTHORS:
    290296
    291297- David Joyner (2006-06)
     298- Stefan Reiterer (2010-)
    292299"""
    293300
    294301#*****************************************************************************
    295302#       Copyright (C) 2006 William Stein <wstein@gmail.com>
    296303#                     2006 David Joyner <wdj@usna.edu>
     304#                     2010 Stefan Reiterer <stefan.reiterer@uni-graz.at>
    297305#
    298306#  Distributed under the terms of the GNU General Public License (GPL)
    299307#
     
    314322from sage.rings.rational_field import RationalField
    315323from sage.rings.real_mpfr import RealField
    316324from sage.misc.sage_eval import sage_eval
    317 from sage.rings.all import QQ, ZZ, CDF, RDF
     325from sage.rings.all import QQ, ZZ, RR, CDF, RDF
     326from sage.symbolic.ring import SR
    318327import sage.rings.commutative_ring as commutative_ring
    319328import sage.rings.ring as ring
    320329from sage.calculus.calculus import maxima
     330from sage.symbolic.function import BuiltinFunction, GinacFunction, is_inexact
     331from sage.symbolic.expression import is_Expression
     332import sage.functions.special
     333from sage.functions.special import MaximaFunction, meval
     334import numpy,scipy
     335import math
     336#import sage.libs.mpmath.all as mpmath #This line does evil things...
     337from sage.functions.other import floor, gamma, factorial, abs, binomial
     338from sage.functions.other import sqrt, conjugate
     339from sage.functions.trig import sin, cos
     340from sage.functions.log import ln
     341import sage.symbolic.expression as expression
     342from sage.structure.parent import Parent
     343from sage.structure.coerce import parent
    321344
    322345_done = False
    323346def _init():
     
    333356        sage: sage.functions.orthogonal_polys._done
    334357        False
    335358
    336     Then after using one of these functions, it changes::
     359    Then after using one of these functions, it changes:: (The value is now
     360    False for chebyshev_T because chebyshev_T uses clenshaw method instead...)
    337361
    338362        sage: from sage.functions.orthogonal_polys import chebyshev_T
    339363        sage: chebyshev_T(2,x)
    340         2*(x - 1)^2 + 4*x - 3
     364        2*x^2 - 1
     365        sage: sage.functions.orthogonal_polys._done
     366        False
     367        sage: legendre_Q(1,x)
     368        1/2*x*log(-(x + 1)/(x - 1)) - 1
    341369        sage: sage.functions.orthogonal_polys._done
    342370        True
    343 
     371   
    344372    Note that because here we use a Pynac variable ``x``,
    345373    the representation of the function is different from
    346374    its actual doctest, where a polynomial indeterminate
     
    355383    maxima.eval('orthopoly_returns_intervals:false;')
    356384    _done = True
    357385
     386#load /home/maldun/sage/sage-4.5.2/devel/sage-ortho/sage/functions/orthogonal_polys.py
    358387
    359 def chebyshev_T(n,x):
     388class OrthogonalPolynomial(BuiltinFunction):
    360389    """
    361     Returns the Chebyshev function of the first kind for integers
    362     `n>-1`.
     390    Base Class for Orthogonal Polynomials. The evaluation as a polynomial
     391    is done via maxima due performance reasons. Therefore the internal name
     392    in maxima maxima_name has to be declared.
     393    Convention: The first argument is always the order of the polynomial,
     394    he last one is always the value x where the polynomial is evaluated.
     395
     396    """
     397    def __init__(self, name, nargs = 2, latex_name = None, conversions = {}):
     398        try:
     399            self._maxima_name = conversions['maxima']
     400        except KeyError:
     401            self._maxima_name = None
     402
     403        BuiltinFunction.__init__(self, name = name,
     404        nargs = nargs, latex_name = latex_name, conversions = conversions)
     405
     406    def _maxima_init_evaled_(self, *args):
     407        """
     408        Returns a string which represents this function evaluated at
     409        *args* in Maxima.
     410        In fact these are thought to be the old wrappers for the orthogonal
     411        polynomials. These are used when the other evaluation methods fail,
     412        or are not fast enough. Experiments showed that for the symbolic
     413        evaluation for larger n maxima is faster, but for small n simply use
     414        of the recursion formulas is faster. A little switch does the trick...
     415
     416        EXAMPLES::
     417
     418            sage: chebyshev_T(3,x)
     419            4*x^3 - 3*x
     420        """
     421        return None
     422
     423    def _clenshaw_method_(self,*args):
     424        """
     425        The Clenshaw method uses the three term recursion of the polynomial,
     426        or explicit formulas instead of maxima to evaluate the polynomial
     427        efficiently, if the x argument is not a symbolic expression.
     428        The name comes from the Clenshaw algorithm for fast evaluation of
     429        polynomialsin chebyshev form.
     430       
     431        comparison: Maxima and Clenshaw algorithm for non-symbolic evaluation:
     432         #sage: time chebyshev_T(50,10) #clenshaw
     433         #CPU times: user 0.00 s, sys: 0.00 s, total: 0.00 s
     434         #Wall time: 0.00 s
     435         #49656746733678312490954442369580252421769338391329426325400124999
     436         #sage: time sage.functions.orthogonal_polys.chebyshev_T(50,10)
     437         #maxima
     438         #CPU times: user 0.00 s, sys: 0.00 s, total: 0.00 s
     439         #Wall time: 0.05 s
     440         #49656746733678312490954442369580252421769338391329426325400124999
     441
     442         #sage: time chebyshev_T(500,10); #clenshaw
     443         #CPU times: user 0.00 s, sys: 0.00 s, total: 0.00 s
     444         #Wall time: 0.00 s
     445         #sage: time sage.functions.orthogonal_polys.chebyshev_T(500,10);
     446         #maxima
     447         #CPU times: user 0.11 s, sys: 0.00 s, total: 0.11 s
     448         #Wall time: 0.77 s
     449        """
     450        raise NotImplementedError(
     451            "No recursive calculation of values implemented (yet)!")
     452
     453    def _eval_special_values_(self,*args):
     454        """
     455        Evals the polynomial explicitly for special values.
     456        EXAMPLES:
     457           
     458            sage: var('n')
     459            n
     460            sage: chebyshev_T(n,-1)
     461            (-1)^n
     462        """
     463        raise ValueError("No special values known!")
     464           
     465
     466    def _eval_(self, *args):
     467        """
     468       
     469        The symbolic evaluation is done with maxima, because the evaluation of
     470        the Polynomial representation seems to be quite clever.
     471        For the fast numerical evaluation an other method should be used...
     472        Therefore I suggest Clenshaw's algorithm, which uses the rekursion!
     473        The function also checks for special values, and if
     474        the order is an integer and in range!
     475       
     476        performance:
     477        #sage: time chebyshev_T(5,x) #maxima
     478        #CPU times: user 0.00 s, sys: 0.00 s, total: 0.00 s
     479        #Wall time: 0.16 s
     480        #16*(x - 1)^5 + 80*(x - 1)^4 + 140*(x - 1)^3 + 100*(x - 1)^2 + 25*x - 24
     481
     482        #sage: time chebyshev_T(5,x) #clenshaw
     483        #CPU times: user 0.01 s, sys: 0.00 s, total: 0.01 s
     484        #Wall time: 0.01 s
     485        #16*x^5 - 20*x^3 + 5*x
     486
     487        #time chebyshev_T(50,x)
     488        #CPU times: user 0.03 s, sys: 0.00 s, total: 0.03 s
     489        #Wall time: 0.04 s
     490        #562949953421312*(x - 1)^50 + 28147497671065600*(x - 1)^49 +....
     491
     492        #time chebyshev_T(100,x);
     493        #CPU times: user 0.08 s, sys: 0.00 s, total: 0.08 s
     494        #Wall time: 0.08 s
     495
     496        EXAMPLES::
     497            sage: chebyshev_T(5,x)
     498            16*x^5 - 20*x^3 + 5*x 
     499            sage: var('n')
     500            n
     501            sage: chebyshev_T(n,-1)
     502            (-1)^n
     503            sage: chebyshev_T(-7,x)
     504            chebyshev_T(-7, x)
     505            sage: chebyshev_T(3/2,x)
     506            chebyshev_T(3/2, x)
     507   
     508        """
     509       
     510        if not is_Expression(args[0]):
     511           
     512            if not is_Expression(args[-1]) and is_inexact(args[-1]):
     513                try:
     514                    import sage.libs.mpmath.all as mpmath
     515                    return self._evalf_(*args)
     516                except AttributeError:
     517                    pass
     518                except mpmath.NoConvergence:
     519                    print "Warning: mpmath returns NoConvergence!"
     520                    print "Switching to clenshaw_method, but it \
     521                    may not be stable!"
     522                except ValueError:
     523                    pass
     524
     525            #A faster check would be nice...
     526            if args[0] != floor(args[0]):
     527                if not is_Expression(args[-1]):
     528                    try:
     529                        return self._evalf_(*args)
     530                    except AttributeError:
     531                        pass
     532                else:
     533                    return None
     534       
     535            if args[0] < 0:
     536                return None
     537           
     538           
     539        try:
     540            return self._eval_special_values_(*args)
     541        except ValueError:
     542            pass
     543
     544   
     545        if not is_Expression(args[0]):
     546           
     547            try:
     548                return self._clenshaw_method_(*args)
     549            except NotImplementedError:
     550                pass
     551       
     552        if self._maxima_name is None:
     553            return None
     554        else:
     555            _init()
     556            try:
     557                #s = maxima(self._maxima_init_evaled_(*args))
     558                #This above is very inefficient! The older
     559                #methods were much faster...
     560                return self._maxima_init_evaled_(*args)
     561            except TypeError:
     562                return None
     563            if self._maxima_name in repr(s):
     564                return None
     565            else:
     566                return s.sage()
     567
     568    #def _eval_numpy_(self, *args):
     569        #TODO: numpy_eval with help of the new scipy version!!!!
     570        #Reason scipy 0.8 supports stable and fast numerical evaluation
     571        #of ortho polys
     572        #Now this only evaluates the array pointwise, and only the first one...
     573        #if isinstance(arg[0], numpy.ndarray)
     574       
     575       
     576       
     577       
     578class Func_chebyshev_T(OrthogonalPolynomial):
     579
     580    """
     581    Class for the Chebyshev polynomial of the first kind.
    363582   
    364583    REFERENCE:
    365584
     
    367586   
    368587    EXAMPLES::
    369588   
    370         sage: x = PolynomialRing(QQ, 'x').gen()
    371         sage: chebyshev_T(2,x)
    372         2*x^2 - 1
     589       sage: chebyshev_T(3,x)
     590       4*x^3 - 3*x
     591       sage: var('k')
     592       k
     593       sage: test = chebyshev_T(k,x)
     594       sage: test
     595       chebyshev_T(k, x)
    373596    """
    374     _init()
    375     return sage_eval(maxima.eval('chebyshev_t(%s,x)'%ZZ(n)), locals={'x':x})
     597   
     598    def __init__(self):
     599        OrthogonalPolynomial.__init__(self,"chebyshev_T",nargs = 2,
     600conversions =dict(maxima='chebyshev_t',mathematica='ChebyshevT'))
     601   
     602    def _eval_special_values_(self,*args):
     603        """
     604        Values known for special values of x.
     605        For details see A.S. 22.4 (p. 777)
     606       
     607        EXAMPLES:
     608           
     609            sage: var('n')
     610            n
     611            sage: chebyshev_T(n,1)
     612            1
     613            sage: chebyshev_T(n,-1)
     614            (-1)^n
     615        """
     616        if args[-1] == 1:
     617            return 1
     618       
     619        if args[-1] == -1:
     620            return (-1)**args[0]
    376621
    377 def chebyshev_U(n,x):
     622        if (args[-1] == 0):
     623            return (1+(-1)**args[0])*(-1)**(args[0]/2)/2
     624
     625        raise ValueError("Value not found!")
     626   
     627    def _evalf_(self, *args,**kwds):
     628        """
     629        Evals chebyshev_T
     630        numerically with mpmath.
     631        EXAMPLES::
     632            sage: chebyshev_T(10,3).n(75)
     633            2.261953700000000000000e7
     634        """
     635        try:
     636            step_parent = kwds['parent']
     637        except KeyError:
     638            step_parent = parent(args[-1])
     639
     640        import sage.libs.mpmath.all as mpmath
     641
     642        try:
     643            precision = step_parent.prec()
     644        except AttributeError:
     645            precision = mpmath.mp.prec
     646
     647        return mpmath.call(mpmath.chebyt,args[0],args[-1],prec = precision)
     648
     649    def _maxima_init_evaled_(self, *args):
     650        n = args[0]
     651        x = args[1]
     652        return sage_eval(maxima.eval('chebyshev_t(%s,x)'%ZZ(n)), locals={'x':x})
     653
     654    def _clenshaw_method_(self,*args):
     655            """
     656            Clenshaw method for chebyshev_T (means use recursions in this case)
     657            This is much faster for numerical evaluation than maxima!
     658            See A.S. 227 (p. 782) for details for the recurions
     659            """
     660
     661            k = args[0]
     662            x = args[1]
     663
     664            if k == 0:
     665                return 1
     666            elif k == 1:
     667                return x
     668            else:
     669                #TODO: When evaluation of Symbolic Expressions works better
     670                #use these explicit formulas instead!
     671                #if -1 <= x <= 1:
     672                #    return cos(k*acos(x))
     673                #elif 1 < x:
     674                #    return cosh(k*acosh(x))
     675                #else: # x < -1
     676                #    return (-1)**(k%2)*cosh(k*acosh(-x))
     677               
     678                help1 = 1
     679                help2 = x
     680                if is_Expression(x):
     681                    #raise NotImplementedError
     682                    help3 = 0
     683                    for j in xrange(0,floor(k/2)+1):
     684                        help3 = \
     685help3 +(-1)**j*(2*x)**(k-2*j)*factorial(k-j-1)/factorial(j)/factorial(k-2*j)
     686                    help3 = help3*k/2
     687                else:
     688                     for j in xrange(0,k-1):
     689                        help3 = 2*x*help2 - help1
     690                        help1 = help2
     691                        help2 = help3
     692               
     693                return help3
     694
     695    def _eval_numpy_(self, *args):
     696        #TODO: numpy_eval with help of the new scipy version!!!!
     697        #Reason scipy 0.8 supports stable and fast numerical evaluation
     698        #of ortho polys
     699        #Now this only evaluates the array pointwise, and only the first one...
     700        #if isinstance(arg[0], numpy.ndarray). This is a hack to provide compability!
     701        """
     702        EXAMPLES::
     703        sage: import numpy
     704        sage: z = numpy.array([1,2])
     705        sage: z2 = numpy.array([[1,2],[1,2]])
     706        sage: z3 = numpy.array([1,2,3.])
     707        sage: chebyshev_T(1,z)
     708        array([1, 2])
     709        sage: chebyshev_T(1,z2)
     710        array([[1, 2],
     711              [1, 2]])
     712        sage: chebyshev_T(1,z3)
     713        array([ 1.,  2.,  3.])
     714
     715        """
     716        if isinstance(args[0],numpy.ndarray):
     717            raise NotImplementedError("Support for numpy array in \
     718            first argument(s) is not supported yet!")
     719
     720        result = numpy.zeros(args[-1].shape).tolist()
     721        if isinstance(args[-1][0],int):
     722            for k in xrange(len(args[-1])):
     723                result[k] = chebyshev_T(args[0],ZZ(args[-1][k]))
     724
     725        if isinstance(args[-1][0],float):
     726            for k in xrange(len(args[-1])):
     727                result[k] = chebyshev_T(args[0],RR(args[-1][k]))
     728
     729        if isinstance(args[-1][0],numpy.ndarray):
     730            for k in xrange(len(args[-1])):
     731                result[k] = chebyshev_T(args[0],args[-1][k])
     732
     733        return numpy.array(result)
     734
     735    def _derivative_(self, *args, **kwds):
     736        """
     737        Returns the derivative of chebyshev_T in form of the chebyshev Polynomial
     738        of the second kind chebyshev_U
     739        EXAMPLES::
     740            sage: var('k')
     741            k
     742            sage: derivative(chebyshev_T(k,x),x)
     743            k*chebyshev_U(k - 1, x)
     744            sage: derivative(chebyshev_T(3,x),x)
     745            12*x^2 - 3
     746            sage: derivative(chebyshev_T(k,x),k)
     747            Traceback (most recent call last):
     748            ...
     749            NotImplementedError: Derivative w.r.t. to the index is not supported, yet, and perhaps never will be...
     750 
     751        """
     752        diff_param = kwds['diff_param']
     753        if diff_param == 0:
     754            raise NotImplementedError(
     755"Derivative w.r.t. to the index is not supported, yet, \
     756and perhaps never will be...")
     757        else:
     758            return args[0]*chebyshev_U(args[0]-1,args[1])
     759
     760   
     761chebyshev_T = Func_chebyshev_T()
     762
     763class Func_chebyshev_U(OrthogonalPolynomial):
     764   
    378765    """
    379     Returns the Chebyshev function of the second kind for integers `n>-1`.
     766    Class for the Chebyshev polynomial of the second kind.
    380767   
    381768    REFERENCE:
    382769
     
    384771   
    385772    EXAMPLES::
    386773   
     774        sage: chebyshev_U(3,x)
     775        8*x^3 - 4*x
     776    """
     777
     778    def __init__(self):
     779        OrthogonalPolynomial.__init__(self,"chebyshev_U",
     780nargs = 2,conversions =dict(maxima='chebyshev_u',mathematica='ChebyshevU'))
     781       
     782    def _clenshaw_method_(self,*args):
     783        """
     784        Clenshaw method for chebyshev_U (means use the recursion...)
     785        This is much faster for numerical evaluation than maxima!
     786        See A.S. 227 (p. 782) for details for the recurions
     787        """
     788        k = args[0]
     789        x = args[1]
     790           
     791        if k == 0:
     792            return 1
     793        elif k == 1:
     794            return 2*x
     795        else:
     796            help1 = 1
     797            help2 = 2*x
     798            if is_Expression(x):
     799                #raise NotImplementedError("Maxima is faster here!")
     800                help3 = 0
     801                for j in xrange(0,floor(k/2)+1):
     802                    help3 = \
     803help3 + (-1)**j*(2*x)**(k-2*j)*factorial(k-j)/factorial(j)/factorial(k-2*j)
     804               
     805            else:
     806                for j in xrange(0,k-1):
     807                    help3 = 2*x*help2 - help1
     808                    help1 = help2
     809                    help2 = help3
     810                   
     811            return help3   
     812   
     813    def _maxima_init_evaled_(self, *args):
     814        """
     815        Uses 
     816        EXAMPLES::
    387817        sage: x = PolynomialRing(QQ, 'x').gen()
    388         sage: chebyshev_U(2,x)
    389         4*x^2 - 1
     818        sage: chebyshev_T(2,x)
     819        2*x^2 - 1
     820        """
     821        n = args[0]
     822        x = args[1]
     823        return sage_eval(maxima.eval('chebyshev_u(%s,x)'%ZZ(n)), locals={'x':x})
     824
     825       
     826    def _evalf_(self, *args,**kwds):
     827        """
     828        Evals chebyshev_U
     829        numerically with mpmath.
     830        EXAMPLES::
     831            sage: chebyshev_U(10,3).n(75)
     832            4.661117900000000000000e7
     833        """
     834        try:
     835            step_parent = kwds['parent']
     836        except KeyError:
     837            step_parent = parent(args[-1])
     838
     839        import sage.libs.mpmath.all as mpmath
     840
     841        try:
     842            precision = step_parent.prec()
     843        except AttributeError:
     844            precision = mpmath.mp.prec
     845
     846        return mpmath.call(mpmath.chebyu,args[0],args[-1],prec = precision)
     847
     848    def _eval_special_values_(self,*args):
     849        """
     850        Special values known. A.S. 22.4 (p.777).
     851        EXAMPLES:
     852       
     853            sage: var('n')
     854            n
     855            sage: chebyshev_U(n,1)
     856            n + 1
     857            sage: chebyshev_U(n,-1)
     858            (n + 1)*(-1)^n
     859        """
     860        if args[-1] == 1:
     861            return (args[0]+1)
     862       
     863        if args[-1] == -1:
     864            return (-1)**args[0]*(args[0]+1)
     865
     866        if (args[-1] == 0):
     867            return (1+(-1)**args[0])*(-1)**(args[0]/2)/2
     868
     869        raise ValueError("Value not found")
     870
     871    def _eval_numpy_(self, *args):
     872        #TODO: numpy_eval with help of the new scipy version!!!!
     873        #Reason scipy 0.8 supports stable and fast numerical evaluation
     874        #of ortho polys
     875        #Now this only evaluates the array pointwise, and only the first one...
     876        #if isinstance(arg[0], numpy.ndarray).
     877        #This is a hack to provide compability!
     878        """
     879        EXAMPLES::
     880        sage: import numpy
     881        sage: z = numpy.array([1,2])
     882        sage: z2 = numpy.array([[1,2],[1,2]])
     883        sage: z3 = numpy.array([1,2,3.])
     884        sage: chebyshev_U(1,z)
     885        array([2, 4])
     886        sage: chebyshev_U(1,z2)
     887        array([[2, 4],
     888              [2, 4]])
     889        sage: chebyshev_U(1,z3)
     890        array([ 2.,  4.,  6.])
     891
     892        """
     893        if isinstance(args[0],numpy.ndarray):
     894            raise NotImplementedError(
     895           "Support for numpy array in first argument(s) is not supported yet!")
     896
     897        result = numpy.zeros(args[-1].shape).tolist()
     898        if isinstance(args[-1][0],int):
     899            for k in xrange(len(args[-1])):
     900                result[k] = chebyshev_U(args[0],ZZ(args[-1][k]))
     901
     902        if isinstance(args[-1][0],float):
     903            for k in xrange(len(args[-1])):
     904                result[k] = chebyshev_U(args[0],RR(args[-1][k]))
     905
     906        if isinstance(args[-1][0],numpy.ndarray):
     907            for k in xrange(len(args[-1])):
     908                result[k] = chebyshev_U(args[0],args[-1][k])
     909
     910        return numpy.array(result)
     911
     912
     913    def _derivative_(self, *args, **kwds):
     914        """
     915        Returns the derivative of chebyshev_U in form of the chebyshev
     916        Polynomials of the first and second kind
     917       
     918        EXAMPLES::
     919            sage: var('k')
     920            k
     921            sage: derivative(chebyshev_U(k,x),x)
     922            ((k + 1)*chebyshev_T(k + 1, x) - x*chebyshev_U(k, x))/(x^2 - 1)
     923            sage: derivative(chebyshev_U(3,x),x)
     924            24*x^2 - 4
     925            sage: derivative(chebyshev_U(k,x),k)
     926            Traceback (most recent call last):
     927            ...
     928            NotImplementedError: Derivative w.r.t. to the index is not supported, yet, and perhaps never will be...
     929 
     930        """
     931        diff_param = kwds['diff_param']
     932        if diff_param == 0:
     933            raise NotImplementedError(
     934"Derivative w.r.t. to the index is not supported, \
     935yet, and perhaps never will be...")
     936        else:
     937            return ((args[0]+1)*chebyshev_T(args[0]+1,args[1])-args[1]*
     938                    chebyshev_U(args[0],args[1]))/(args[1]**2-1)
     939       
     940chebyshev_U = Func_chebyshev_U()
     941
     942class Func_gegenbauer(OrthogonalPolynomial):
    390943    """
    391     _init()
    392     return sage_eval(maxima.eval('chebyshev_u(%s,x)'%ZZ(n)), locals={'x':x})
     944    Returns the ultraspherical (or Gegenbauer) polynomial.
     945   
     946    REFERENCE:
    393947
    394 def gen_laguerre(n,a,x):
     948    - AS 22.5.27
     949   
     950    EXAMPLES::
     951   
     952        sage: x = PolynomialRing(QQ, 'x').gen()
     953        sage: ultraspherical(2,3/2,x)
     954        15/2*x^2 - 3/2
     955        sage: ultraspherical(2,1/2,x)
     956        3/2*x^2 - 1/2
     957        sage: ultraspherical(1,1,x)
     958        2*x     
     959        sage: t = PolynomialRing(RationalField(),"t").gen()
     960        sage: gegenbauer(3,2,t)
     961        32*t^3 - 12*t
     962    """   
     963    def __init__(self):
     964        OrthogonalPolynomial.__init__(self,"gegenbauer",nargs = 3,
     965        conversions =dict(maxima='ultraspherical',mathematica='GegenbauerC'))
     966
     967    def _clenshaw_method_(self,*args):
     968        """
     969        Clenshaw method for gegenbauer poly (means use the recursion...)
     970        This is much faster for numerical evaluation than maxima!
     971        See A.S. 227 (p. 782) for details for the recurions
     972        """
     973        k = args[0]
     974        x = args[-1]
     975        alpha = args[1]
     976
     977        if is_Expression(alpha) or abs(alpha) > numpy.finfo(float).eps:
     978            alpha_zero = False
     979        else:
     980            alpha_zero = True
     981
     982        if k == 0:
     983            return 1
     984        elif k == 1:
     985            if not alpha_zero:
     986                return 2*x*alpha
     987            else:
     988                return 2*x # It seems that maxima evals this the wrong way!
     989                           #(see A.S. 22.4 (p.777))
     990        else:
     991            help1 = 1
     992            if alpha_zero:
     993                help2 = 2*x
     994                gamma_alpha = 1
     995            else:
     996                help2 = 2*x*alpha
     997                gamma_alpha = gamma(alpha)
     998
     999            help3 = 0
     1000            if is_Expression(x):
     1001                for j in xrange(0,floor(k/2)+1):
     1002                    help3 = \
     1003help3 + (-1)**j*gamma(alpha+k-j)/factorial(j)/factorial(k-2*j)*(2*x)**(k-2*j)/gamma_alpha
     1004               
     1005            else:
     1006                for j in xrange(1,k):
     1007                    help3 = 2*(j+alpha)*x*help2 - (j+2*alpha-1)*help1
     1008                    help3 = help3/(j+1)
     1009                    help1 = help2
     1010                    help2 = help3
     1011                   
     1012            return help3   
     1013
     1014    def _maxima_init_evaled_(self, *args):
     1015        """
     1016        Uses 
     1017        EXAMPLES::
     1018        sage: x = PolynomialRing(QQ, 'x').gen()
     1019        sage: ultraspherical(2,1/2,x)
     1020        3/2*x^2 - 1/2
     1021        """
     1022        n = args[0]
     1023        a = args[1]
     1024        x = args[2]
     1025        return sage_eval(maxima.eval('ultraspherical(%s,%s,x)'%(ZZ(n),a)),\
     1026                         locals={'x':x})
     1027
     1028
     1029    def _evalf_(self, *args,**kwds):
     1030        """
     1031        Evals gegenbauer
     1032        numerically with mpmath.
     1033        EXAMPLES::
     1034            sage: gegenbauer(10,2,3.).n(54)
     1035            5.25360702000000e8
     1036        """
     1037        if args[1] == 0:
     1038            raise ValueError("mpmath don't handle alpha = 0 correctly!")
     1039           
     1040        try:
     1041            step_parent = kwds['parent']
     1042        except KeyError:
     1043            step_parent = parent(args[-1])
     1044
     1045        import sage.libs.mpmath.all as mpmath
     1046
     1047        try:
     1048            precision = step_parent.prec()
     1049        except AttributeError:
     1050            precision = mpmath.mp.prec
     1051
     1052        return mpmath.call(
     1053            mpmath.gegenbauer,args[0],args[1],args[-1],prec = precision)
     1054
     1055    def _eval_special_values_(self,*args):
     1056        """
     1057        Special values known. A.S. 22.4 (p.777)
     1058        EXAMPLES:
     1059       
     1060            sage: var('n a')
     1061            (n, a)
     1062            sage: gegenbauer(n,1/2,x)
     1063            legendre_P(n, x)
     1064            sage: gegenbauer(n,0,x)
     1065            1/2*n*chebyshev_T(n, x)
     1066            sage: gegenbauer(n,1,x)
     1067            chebyshev_U(n, x)
     1068            sage: gegenbauer(n,a,1)
     1069            binomial(2*a + n - 1, n)
     1070            sage: gegenbauer(n,a,-1)
     1071            (-1)^n*binomial(2*a + n - 1, n)
     1072            sage: gegenbauer(n,a,0)
     1073            1/2*((-1)^n + 1)*(-1)^(1/2*n)*gamma(a + 1/2*n)/(gamma(a)*gamma(1/2*n))
     1074        """
     1075        if args[1] == 0 and args[0] != 0:
     1076            return args[0]*chebyshev_T(args[0],args[-1])/2
     1077
     1078        if args[1] == 0.5:
     1079            return legendre_P(args[0],args[-1])
     1080
     1081        if args[1] == 1:
     1082            return chebyshev_U(args[0],args[-1])
     1083
     1084        if args[-1] == 1:
     1085            return binomial(args[0] + 2*args[1] - 1,args[0])
     1086       
     1087        if args[-1] == -1:
     1088            return (-1)**args[0]*binomial(args[0] + 2*args[1] - 1,args[0])
     1089
     1090        if (args[-1] == 0):
     1091            return (1+(-1)**args[0])*(-1)**(args[0]/2)/2*\
     1092            gamma(args[1]+args[0]/2)/gamma(args[1])/gamma(args[0]/2)
     1093
     1094        raise ValueError("Value not found")
     1095
     1096    def _eval_numpy_(self, *args):
     1097        #TODO: numpy_eval with help of the new scipy version!!!!
     1098        #Reason scipy 0.8 supports stable and fast numerical evaluation
     1099        #of ortho polys
     1100        #Now this only evaluates the array pointwise, and only the first one...
     1101        #if isinstance(arg[0], numpy.ndarray). This is a hack to provide compability!
     1102        """
     1103        EXAMPLES::
     1104        sage: import numpy
     1105        sage: z = numpy.array([1,2])
     1106        sage: z2 = numpy.array([[1,2],[1,2]])
     1107        sage: z3 = numpy.array([1,2,3.])
     1108        sage: gegenbauer(2,0,z)
     1109        array([1, 7])
     1110        sage: gegenbauer(2,1,z)
     1111        array([ 3, 15])
     1112        sage: gegenbauer(2,1,z2)
     1113        array([[ 3, 15],
     1114              [ 3, 15]])
     1115        sage: gegenbauer(2,1,z3)
     1116        array([  3.,  15.,  35.])
     1117        """
     1118
     1119        if isinstance(args[0],numpy.ndarray) or isinstance(args[1],numpy.ndarray):
     1120            raise NotImplementedError(
     1121                "Support for numpy array in first argument(s) is not supported yet!")
     1122       
     1123        result = numpy.zeros(args[-1].shape).tolist()
     1124        if isinstance(args[-1][0],int):
     1125            for k in xrange(len(args[-1])):
     1126                result[k] = gegenbauer(args[0],args[1],ZZ(args[-1][k]))
     1127
     1128        if isinstance(args[-1][0],float):
     1129            for k in xrange(len(args[-1])):
     1130                result[k] = gegenbauer(args[0],args[1],RR(args[-1][k]))
     1131
     1132        if isinstance(args[-1][0],numpy.ndarray):
     1133            for k in xrange(len(args[-1])):
     1134                result[k] = gegenbauer(args[0],args[1],args[-1][k])
     1135
     1136        return numpy.array(result)
     1137
     1138    def _derivative_(self, *args, **kwds):
     1139        """
     1140        Returns the derivative of chebyshev_U in form of the chebyshev
     1141        Polynomials of the first and second kind
     1142       
     1143        EXAMPLES::
     1144            sage: var('k a')
     1145            (k, a)
     1146            sage: derivative(gegenbauer(k,a,x),x)
     1147            (k*x*gegenbauer(k, a, x) - (2*a + k - 1)*gegenbauer(k - 1, a, x))/(x^2 - 1)
     1148            sage: derivative(gegenbauer(4,3,x),x)
     1149            960*x^3 - 240*x
     1150            sage: derivative(gegenbauer(k,a,x),a)
     1151            Traceback (most recent call last):
     1152            ...
     1153            NotImplementedError: Derivative w.r.t. to the index is not supported, yet, and perhaps never will be...
     1154 
     1155        """
     1156        diff_param = kwds['diff_param']
     1157        if diff_param in [0,1]:
     1158            raise NotImplementedError(
     1159"Derivative w.r.t. to the index is not supported, yet,\
     1160and perhaps never will be...")
     1161        else:
     1162            return (-args[0]*args[-1]*gegenbauer(args[0],args[1],args[2])+\
     1163(args[0] + 2*args[1]-1)*gegenbauer(args[0]-1,args[1],args[2]))/(1-args[-1]**2)
     1164
     1165gegenbauer = Func_gegenbauer()
     1166ultraspherical = Func_gegenbauer()
     1167
     1168class Func_gen_laguerre(OrthogonalPolynomial):
    3951169    """
    3961170    Returns the generalized Laguerre polynomial for integers `n > -1`.
    3971171    Typically, a = 1/2 or a = -1/2.
     
    4141188        sage: gen_laguerre(3,0,x)
    4151189        -1/6*x^3 + 3/2*x^2 - 3*x + 1
    4161190    """
    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})
     1191    def __init__(self):
     1192        OrthogonalPolynomial.__init__(self,"gen_laguerre",nargs = 3,
     1193conversions =dict(maxima='gen_laguerre',mathematica='LaguerreL'))
     1194       
     1195    def _clenshaw_method_(self,*args):
     1196        """
     1197        Clenshaw method for gen_laguerre polynomial (means use the recursion...)
     1198        This is much faster for numerical evaluation than maxima!
     1199        See A.S. 227 (p. 782) for details for the recurions.
     1200        Warning: The clanshaw method for the laguerre Polynomials
     1201        should only used for exact data types, when high orders are
     1202        used, due to weak instabilities of the recursion!
     1203        """
     1204        k = args[0]
     1205        a = args[1]
     1206        x = args[-1]
     1207       
     1208        if k == 0:
     1209            return 1
     1210        elif k == 1:
     1211            return 1-x + a
     1212        else:
     1213            help1 = 1
     1214            help2 = 1-x + a
     1215            if is_Expression(x):
     1216                help3 = 0
     1217                for j in xrange(0,k+1):
     1218                    help3 = help3 + (-1)**j*binomial(k+a,k-j)/factorial(j)*x**j
     1219            else:
     1220                for j in xrange(1,k):
     1221                    help3 = -x*help2 + (2*j+1+a)*help2 - (j+a)*help1
     1222                    help3 = help3/(j+1)
     1223                    help1 = help2
     1224                    help2 = help3
     1225               
     1226            return help3   
    4201227
    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`.
     1228    def _maxima_init_evaled_(self, *args):
     1229        n = args[0]
     1230        a = args[1]
     1231        x = args[-1]
     1232       
     1233        from sage.functions.all import sqrt
     1234        _init()
     1235        return sage_eval(maxima.eval('gen_laguerre(%s,%s,x)'%(ZZ(n),a)),\
     1236                         locals={'x':x})
     1237
     1238    def _evalf_(self, *args,**kwds):
     1239        """
     1240        Evals laguerre polynomial
     1241        numerically with mpmath.
     1242        EXAMPLES::
     1243            sage: laguerre(3,5.).n(53)
     1244            2.66666666666667
     1245        """
     1246        try:
     1247            step_parent = kwds['parent']
     1248        except KeyError:
     1249            step_parent = parent(args[-1])
     1250
     1251        import sage.libs.mpmath.all as mpmath
     1252
     1253        try:
     1254            precision = step_parent.prec()
     1255        except AttributeError:
     1256            precision = mpmath.mp.prec
     1257
     1258        return mpmath.call(
     1259            mpmath.laguerre,args[0],args[1],args[-1],prec = precision)
    4251260   
    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.
     1261    def _eval_special_values_(self,*args):
     1262        """
     1263        Special values known.
     1264        EXAMPLES:
     1265           
     1266            sage: var('n')
     1267            n
     1268            sage: gen_laguerre(n,1,0)
     1269            binomial(n + 1, n)
     1270        """
     1271       
     1272        if (args[-1] == 0):
     1273            try:
     1274                return binomial(args[0]+args[1],args[0])
     1275            except TypeError:
     1276                pass
     1277
     1278        raise ValueError("Value not found")
     1279
     1280    def _eval_numpy_(self, *args):
     1281        #TODO: numpy_eval with help of the new scipy version!!!!
     1282        #Reason scipy 0.8 supports stable and fast numerical evaluation
     1283        #of ortho polys
     1284        #Now this only evaluates the array pointwise, and only the first one...
     1285        #if isinstance(arg[0], numpy.ndarray). This is a hack to provide compability!
     1286        """
     1287        EXAMPLES::
     1288        sage: import numpy
     1289        sage: z = numpy.array([1,2])
     1290        sage: z2 = numpy.array([[1,2],[1,2]])
     1291        sage: z3 = numpy.array([1,2,3.])
     1292        sage: gen_laguerre(1,1,z)
     1293        array([1, 0])
     1294        sage: gen_laguerre(1,1,z2)
     1295        array([[1, 0],
     1296              [1, 0]])
     1297        sage: gen_laguerre(1,1,z3)
     1298        array([ 1.,  0., -1.])
     1299        """
     1300
     1301        if isinstance(args[0],numpy.ndarray) or isinstance(args[1],numpy.ndarray):
     1302            raise NotImplementedError(
     1303           "Support for numpy array in first argument(s) is not supported yet!")
     1304       
     1305        result = numpy.zeros(args[-1].shape).tolist()
     1306        if isinstance(args[-1][0],int):
     1307            for k in xrange(len(args[-1])):
     1308                result[k] = gen_laguerre(args[0],args[1],ZZ(args[-1][k]))
     1309
     1310        if isinstance(args[-1][0],float):
     1311            for k in xrange(len(args[-1])):
     1312                result[k] = gen_laguerre(args[0],args[1],RR(args[-1][k]))
     1313
     1314        if isinstance(args[-1][0],numpy.ndarray):
     1315            for k in xrange(len(args[-1])):
     1316                result[k] = gen_laguerre(args[0],args[1],args[-1][k])
     1317
     1318        return numpy.array(result)
     1319
    4311320   
    432     REFERENCE:
    4331321
    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))
     1322    def _derivative_(self,*args,**kwds):
     1323        """return the derivative of gen_laguerre in
     1324        form of the gen. Laguerre Polynomial.
     1325        EXAMPLES::
     1326        sage: n = var('n')
     1327        sage: derivative(gen_laguerre(3,1,x),x)
     1328        -1/2*x^2 + 4*x - 6
     1329        sage: derivative(gen_laguerre(n,1,x),x)
     1330        -((n + 1)*gen_laguerre(n - 1, 1, x) - n*gen_laguerre(n, 1, x))/x
     1331        """
     1332        diff_param = kwds['diff_param']
     1333        if diff_param == 0:
     1334            raise NotImplementedError(
     1335"Derivative w.r.t. to the index is not supported, \
     1336yet, and perhaps never will be...")
     1337        else:
     1338            return (args[0]*gen_laguerre(args[0],args[1],args[-1])-(args[0]+args[1])*\
     1339                    gen_laguerre(args[0]-1,args[1],args[-1]))/args[-1]
    4581340
    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)
     1341gen_laguerre = Func_gen_laguerre()
    4971342
    498 def hermite(n,x):
     1343class Func_hermite(OrthogonalPolynomial):
    4991344    """
    5001345    Returns the Hermite polynomial for integers `n > -1`.
    5011346   
     
    5031348
    5041349    - AS 22.5.40 and 22.5.41, page 779.
    5051350   
     1351    #sage: S.<y> = PolynomialRing(RR)
     1352    #   <- Here a bug of the base class BuiltinFunction occours!
     1353
    5061354    EXAMPLES::
    5071355   
    5081356        sage: x = PolynomialRing(QQ, 'x').gen()
     
    5121360        8*x^3 - 12*x
    5131361        sage: hermite(3,2)
    5141362        40
    515         sage: S.<y> = PolynomialRing(RR)
    516         sage: hermite(3,y)
     1363        sage: y = var('y')
     1364        sage: hermite(3,y).polynomial(PolynomialRing(RR, 'y'))
    5171365        8.00000000000000*y^3 - 12.0000000000000*y
    5181366        sage: R.<x,y> = QQ[]
    5191367        sage: hermite(3,y^2)
    5201368        8*y^6 - 12*y^2
    5211369        sage: w = var('w')
    5221370        sage: hermite(3,2*w)
    523         8*(8*w^2 - 3)*w
     1371        64*w^3 - 24*w
    5241372    """
    525     _init()
    526     return sage_eval(maxima.eval('hermite(%s,x)'%ZZ(n)), locals={'x':x})
     1373   
     1374    def __init__(self):
     1375        OrthogonalPolynomial.__init__(self,"hermite",nargs = 2,
     1376conversions =dict(maxima='hermite',mathematica='HermiteH'))
    5271377
    528 def jacobi_P(n,a,b,x):
     1378    def _clenshaw_method_(self,*args):
     1379        """
     1380        Clenshaw method for hermite polynomial (means use the recursion...)
     1381        See A.S. 227 (p. 782) for details for the recurions
     1382        For the symbolic evaluation, maxima seems to be quite fast.
     1383        The break even point between the recursion and Maxima is about
     1384        n = 25
     1385        """
     1386        k = args[0]
     1387        x = args[1]
     1388           
     1389        if k == 0:
     1390            return 1
     1391        elif k == 1:
     1392            return 2*x
     1393        else:
     1394            help1 = 1
     1395            help2 = 2*x
     1396            if is_Expression(x):
     1397                help3 = 0
     1398                for j in xrange(0,floor(k/2)+1):
     1399                    help3 = help3 + (-1)**j*(2*x)**(k-2*j)/factorial(j)/\
     1400                            factorial(k-2*j)
     1401                help3 = help3*factorial(k)
     1402            else:
     1403                for j in xrange(1,k):
     1404                    help3 = 2*x*help2 - 2*j*help1
     1405                    help1 = help2
     1406                    help2 = help3
     1407                   
     1408            return help3   
     1409   
     1410           
     1411    def _evalf_(self, *args,**kwds):
     1412        """
     1413        Evals hermite
     1414        numerically with mpmath.
     1415        EXAMPLES::
     1416            sage: hermite(3,2.).n(74)
     1417            40.0000000000000000000
     1418        """
     1419        try:
     1420            step_parent = kwds['parent']
     1421        except KeyError:
     1422            step_parent = parent(args[-1])
     1423
     1424        import sage.libs.mpmath.all as mpmath
     1425
     1426        try:
     1427            precision = step_parent.prec()
     1428        except AttributeError:
     1429            precision = mpmath.mp.prec
     1430
     1431        return mpmath.call(mpmath.hermite,args[0],args[-1],prec = precision)
     1432   
     1433    def _eval_special_values_(self,*args):
     1434        """
     1435        Special values known. A.S. 22.4 (p.777)
     1436        EXAMPLES:
     1437       
     1438            sage: var('n')
     1439            n
     1440            sage: hermite(n,0)
     1441            ((-1)^n + 1)*(-1)^(1/2*n)*factorial(n)/gamma(1/2*n + 1)
     1442        """
     1443
     1444        if (args[-1] == 0):
     1445            return (1+(-1)**args[0])*(-1)**(args[0]/2)*\
     1446                   factorial(args[0])/gamma(args[0]/2+1)
     1447
     1448        raise ValueError("Value not found")
     1449
     1450    def _maxima_init_evaled_(self, *args):
     1451        """
     1452        Old maxima method.
     1453        """
     1454        n = args[0]
     1455        x = args[1]
     1456        return sage_eval(maxima.eval('hermite(%s,x)'%ZZ(n)), locals={'x':x})
     1457
     1458    def _eval_numpy_(self, *args):
     1459        #TODO: numpy_eval with help of the new scipy version!!!!
     1460        #Reason scipy 0.8 supports stable and fast numerical evaluation
     1461        #of ortho polys
     1462        #Now this only evaluates the array pointwise, and only the first one...
     1463        #if isinstance(arg[0], numpy.ndarray).
     1464        #This is a hack to provide compability!
     1465        """
     1466        EXAMPLES::
     1467        sage: import numpy
     1468        sage: z = numpy.array([1,2])
     1469        sage: z2 = numpy.array([[1,2],[1,2]])
     1470        sage: z3 = numpy.array([1,2,3.])
     1471        sage: hermite(1,z)
     1472        array([2, 4])
     1473        sage: hermite(1,z2)
     1474        array([[2, 4],
     1475              [2, 4]])
     1476        sage: hermite(1,z3)
     1477        array([ 2.,  4.,  6.])
     1478
     1479        """
     1480        if isinstance(args[0],numpy.ndarray):
     1481            raise NotImplementedError(
     1482          "Support for numpy array in first argument(s) is not supported yet!")
     1483
     1484        result = numpy.zeros(args[-1].shape).tolist()
     1485        if isinstance(args[-1][0],int):
     1486            for k in xrange(len(args[-1])):
     1487                result[k] = hermite(args[0],ZZ(args[-1][k]))
     1488
     1489        if isinstance(args[-1][0],float):
     1490            for k in xrange(len(args[-1])):
     1491                result[k] = hermite(args[0],RR(args[-1][k]))
     1492
     1493        if isinstance(args[-1][0],numpy.ndarray):
     1494            for k in xrange(len(args[-1])):
     1495                result[k] = hermite(args[0],args[-1][k])
     1496
     1497        return numpy.array(result)
     1498
     1499
     1500    def _derivative_(self, *args, **kwds):
     1501        """
     1502        Returns the derivative of the hermite polynomial in form of the chebyshev
     1503        Polynomials of the first and second kind
     1504       
     1505        EXAMPLES::
     1506            sage: var('k')
     1507            k
     1508            sage: derivative(hermite(k,x),x)
     1509            2*k*hermite(k - 1, x)
     1510             
     1511        """
     1512        diff_param = kwds['diff_param']
     1513        if diff_param == 0:
     1514            raise NotImplementedError(
     1515"Derivative w.r.t. to the index is not supported, yet, and perhaps never will be...")
     1516        else:
     1517            return 2*args[0]*hermite(args[0]-1,args[1])
     1518 
     1519hermite = Func_hermite()       
     1520
     1521
     1522class Func_jacobi_P(OrthogonalPolynomial):
    5291523    r"""
    5301524    Returns the Jacobi polynomial `P_n^{(a,b)}(x)` for
    5311525    integers `n > -1` and a and b symbolic or `a > -1`
     
    5461540        sage: jacobi_P(2,1,2,1.2)        # random output of low order bits
    5471541        5.009999999999998
    5481542    """
    549     _init()
    550     return sage_eval(maxima.eval('jacobi_p(%s,%s,%s,x)'%(ZZ(n),a,b)), locals={'x':x})
    5511543
    552 def laguerre(n,x):
     1544    def __init__(self):
     1545        OrthogonalPolynomial.__init__(self,"jacobi_P",nargs = 4,
     1546conversions =dict(maxima='jacobi_p',mathematica='JacobiP'))
     1547   
     1548    def _clenshaw_method_(self,*args):
     1549        """
     1550        Clenshaw method for jacobi_P (means use the recursion,
     1551        or sum formula)
     1552        This is much faster for numerical evaluation than maxima!
     1553        See A.S. 227 (p. 782) for details for the recurions.
     1554        Warning: The clanshaw method for the Jacobi Polynomials
     1555        should only used for exact data types, when high orders are
     1556        used, due to weak instabilities of the recursion!
     1557        """
     1558        k = args[0]
     1559        x = args[-1]
     1560        alpha = args[1]
     1561        beta = args[2]
     1562
     1563        if k == 0:
     1564            return 1
     1565        elif k == 1:
     1566            return (alpha-beta + (alpha+beta+2)*x)/2
     1567        else:
     1568           
     1569            if is_Expression(x) or is_Expression(alpha) or is_Expression(beta):
     1570                #Here we use the sum formula of jacobi_P it seems this is rather
     1571                #optimal for use.
     1572                help1 = gamma(alpha+k+1)/factorial(k)/gamma(alpha+beta+k+1)
     1573                help2 = 0
     1574                for j in xrange(0,k+1):
     1575                    help2 = help2 + binomial(k,j)*gamma(alpha+beta+k+j+1)/\
     1576                            gamma(alpha+j+1)*((x-1)/2)**j
     1577                return help1*help2
     1578            else:
     1579                help1 = 1
     1580                help2 = (alpha-beta + (alpha+beta+2)*x)/2
     1581
     1582                for j in xrange(1,k):
     1583                    a1 = 2*(j+1)*(j+alpha+beta+1)*(2*j+alpha+beta)
     1584                    a2 = (2*j+alpha+beta+1)*(alpha**2-beta**2)
     1585                    a3 = (2*j+alpha+beta)*(2*j+alpha+beta+1)*(2*j+alpha+beta+2)
     1586                    a4 = 2*(j+alpha)*(j+beta)*(2*j+alpha+beta+2)
     1587                    help3 = (a2+a3*x)*help2 - a4*help1
     1588                    help3 = help3/a1
     1589                    help1 = help2
     1590                    help2 = help3
     1591               
     1592            return help3   
     1593
     1594    def _maxima_init_evaled_(self, *args):
     1595        """
     1596        Old maxima method.
     1597        """
     1598        n = args[0]
     1599        a = args[1]
     1600        b = args[2]
     1601        x = args[3]
     1602        return sage_eval(maxima.eval('jacobi_p(%s,%s,%s,x)'%(ZZ(n),a,b)),\
     1603                         locals={'x':x})
     1604
     1605    def _evalf_(self, *args,**kwds):
     1606        """
     1607        Evals jacobi_P
     1608        numerically with mpmath.
     1609        EXAMPLES::
     1610            sage: jacobi_P(10,2,3,3).n(75)
     1611            1.322776620000000000000e8
     1612        """
     1613        try:
     1614            step_parent = kwds['parent']
     1615        except KeyError:
     1616            step_parent = parent(args[-1])
     1617
     1618        import sage.libs.mpmath.all as mpmath
     1619
     1620        try:
     1621            precision = step_parent.prec()
     1622        except AttributeError:
     1623            precision = mpmath.mp.prec
     1624
     1625        return mpmath.call(mpmath.jacobi,args[0],args[1],args[2],args[-1],\
     1626                           prec = precision)
     1627
     1628    def _eval_special_values_(self,*args):
     1629        """
     1630        Special values known. A.S. 22.4 (p.777)
     1631        EXAMPLES:
     1632           
     1633            sage: var('n')
     1634            n
     1635            sage: legendre_P(n,1)
     1636            1
     1637            sage: legendre_P(n,-1)
     1638            (-1)^n
     1639        """
     1640        if args[-1] == 1:
     1641            return binomial(args[0]+args[1],args[0])
     1642       
     1643        if args[-1] == -1:
     1644            return (-1)**args[0]*binomial(args[0]+args[1],args[0])
     1645       
     1646        if args[1] == 0 and args[2] == 0:
     1647            return legendre_P(args[0],args[-1])
     1648
     1649        if args[1] == -0.5 and args[2] == -0.5:
     1650            try:
     1651                return binomial(2*args[0],args[0])*\
     1652                       chebyshev_T(args[0],args[-1])/4**args[0]
     1653            except TypeError:
     1654                pass
     1655
     1656        raise ValueError("Value not found")
     1657
     1658    def _eval_numpy_(self, *args):
     1659        #TODO: numpy_eval with help of the new scipy version!!!!
     1660        #Reason scipy 0.8 supports stable and fast numerical evaluation
     1661        #of ortho polys
     1662        #Now this only evaluates the array pointwise, and only the first one...
     1663        #if isinstance(arg[0], numpy.ndarray).
     1664        #This is a hack to provide compability!
     1665        """
     1666        EXAMPLES::
     1667        sage: import numpy
     1668        sage: z = numpy.array([1,2])
     1669        sage: z2 = numpy.array([[1,2],[1,2]])
     1670        sage: z3 = numpy.array([1,2,3.])
     1671        sage: gen_laguerre(1,1,z)
     1672        array([1, 0])
     1673        sage: gen_laguerre(1,1,z2)
     1674        array([[1, 0],
     1675              [1, 0]])
     1676        sage: gen_laguerre(1,1,z3)
     1677        array([ 1.,  0., -1.])
     1678        """
     1679
     1680        if isinstance(args[0],numpy.ndarray) or \
     1681               isinstance(args[1],numpy.ndarray) or \
     1682               isinstance(args[2],numpy.ndarray):
     1683            raise NotImplementedError(
     1684           "Support for numpy array in first argument(s) is not supported yet!")
     1685       
     1686        result = numpy.zeros(args[-1].shape).tolist()
     1687        if isinstance(args[-1][0],int):
     1688            for k in xrange(len(args[-1])):
     1689                result[k] = jacobi_P(args[0],args[1],args[2],ZZ(args[-1][k]))
     1690
     1691        if isinstance(args[-1][0],float):
     1692            for k in xrange(len(args[-1])):
     1693                result[k] = jacobi_P(args[0],args[1],args[2],RR(args[-1][k]))
     1694
     1695        if isinstance(args[-1][0],numpy.ndarray):
     1696            for k in xrange(len(args[-1])):
     1697                result[k] = jacobi_P(args[0],args[1],args[2],args[-1][k])
     1698
     1699        return numpy.array(result)
     1700
     1701    def _derivative_(self, *args, **kwds):
     1702        """
     1703        Returns the derivative of jacobi_P in form of jacobi_polynomials
     1704       
     1705        EXAMPLES::
     1706            sage: var('k a b')
     1707            (k, a, b)
     1708            sage: derivative(jacobi_P(k,a,b,x),x)
     1709            -(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))
     1710            sage: derivative(jacobi_P(2,1,3,x),x)
     1711            14*x - 7/2
     1712            sage: derivative(jacobi_P(k,a,b,x),a)
     1713            Traceback (most recent call last):
     1714            ...
     1715            NotImplementedError: Derivative w.r.t. to the index is not supported, yet, and perhaps never will be...
     1716 
     1717        """
     1718        diff_param = kwds['diff_param']
     1719        if diff_param in [0,1,2]:
     1720            raise NotImplementedError(
     1721"Derivative w.r.t. to the index is not supported, yet, and perhaps never will be...")
     1722        else:
     1723            return (args[0]*(args[1]-args[2]-(2*args[0]+args[1]+args[2])*args[-1])*\
     1724                    jacobi_P(args[0],args[1],args[2],args[-1])+ 2*(args[0]+args[1])*\
     1725                    (args[0]+args[2])*jacobi_P(args[0]-1,args[1],args[2],args[-1]))/\
     1726                    (2*args[0]+args[1]+args[2])/(1-args[-1]**2)
     1727       
     1728       
     1729jacobi_P = Func_jacobi_P()       
     1730
     1731class Func_laguerre(OrthogonalPolynomial):
    5531732    """
    5541733    Returns the Laguerre polynomial for integers `n > -1`.
    5551734   
     
    5671746        sage: laguerre(2,2)
    5681747        -1
    5691748    """
    570     _init()
    571     return sage_eval(maxima.eval('laguerre(%s,x)'%ZZ(n)), locals={'x':x})
     1749    def __init__(self):
     1750        OrthogonalPolynomial.__init__(self,"laguerre",nargs = 2,
     1751conversions =dict(maxima='laguerre',mathematica='LaguerreL'))
     1752       
     1753    def _clenshaw_method_(self,*args):
     1754        """
     1755        Clenshaw method for laguerre polynomial (means use the recursion...)
     1756        This is much faster for numerical evaluation than maxima!
     1757        See A.S. 227 (p. 782) for details for the recurions.
     1758        Warning: The clanshaw method for the laguerre Polynomials
     1759        should only used for exact data types, when high orders are
     1760        used, due to weak instabilities of the recursion!
     1761        """
     1762        k = args[0]
     1763        x = args[-1]
     1764       
     1765        if k == 0:
     1766            return 1
     1767        elif k == 1:
     1768            return 1-x
     1769        else:
     1770            help1 = 1
     1771            help2 = 1-x
     1772            if is_Expression(x):
     1773                help3 = 0
     1774                for j in xrange(0,k+1):
     1775                    help3 = help3 + (-1)**j*binomial(k,k-j)/factorial(j)*x**j
     1776            else:
     1777                for j in xrange(1,k):
     1778                    help3 = -x*help2 + (2*j+1)*help2 - j*help1
     1779                    help3 = help3/(j+1)
     1780                    help1 = help2
     1781                    help2 = help3
     1782               
     1783            return help3   
    5721784
    573 def legendre_P(n,x):
     1785    def _maxima_init_evaled_(self, *args):
     1786        n = args[0]
     1787        x = args[1]
     1788        return sage_eval(maxima.eval('laguerre(%s,x)'%ZZ(n)), locals={'x':x})
     1789
     1790    def _evalf_(self, *args,**kwds):
     1791        """
     1792        Evals laguerre polynomial
     1793        numerically with mpmath.
     1794        EXAMPLES::
     1795            sage: laguerre(3,5.).n(53)
     1796            2.66666666666667
     1797        """
     1798        try:
     1799            step_parent = kwds['parent']
     1800        except KeyError:
     1801            step_parent = parent(args[-1])
     1802
     1803        import sage.libs.mpmath.all as mpmath
     1804
     1805        try:
     1806            precision = step_parent.prec()
     1807        except AttributeError:
     1808            precision = mpmath.mp.prec
     1809
     1810        return mpmath.call(mpmath.laguerre,args[0],0,args[-1],prec = precision)
     1811   
     1812    def _eval_special_values_(self,*args):
     1813        """
     1814        Special values known.
     1815        EXAMPLES:
     1816           
     1817            sage: var('n')
     1818            n
     1819            sage: laguerre(n,0)
     1820            1
     1821        """
     1822       
     1823        if (args[-1] == 0):
     1824            try:
     1825                return 1
     1826            except TypeError:
     1827                pass
     1828
     1829        raise ValueError("Value not found")
     1830 
     1831
     1832    def _eval_numpy_(self, *args):
     1833        #TODO: numpy_eval with help of the new scipy version!!!!
     1834        #Reason scipy 0.8 supports stable and fast numerical evaluation
     1835        #of ortho polys
     1836        #Now this only evaluates the array pointwise, and only the first one...
     1837        #if isinstance(arg[0], numpy.ndarray).
     1838        #This is a hack to provide compability!
     1839        """
     1840        EXAMPLES::
     1841        sage: import numpy
     1842        sage: z = numpy.array([1,2])
     1843        sage: z2 = numpy.array([[1,2],[1,2]])
     1844        sage: z3 = numpy.array([1,2,3.])
     1845        sage: laguerre(1,z)
     1846        array([ 0, -1])
     1847        sage: laguerre(1,z2)
     1848        array([[ 0, -1],
     1849               [ 0, -1]])
     1850        sage: laguerre(1,z3)
     1851        array([ 0., -1., -2.])
     1852
     1853        """
     1854        if isinstance(args[0],numpy.ndarray):
     1855            raise NotImplementedError(
     1856           "Support for numpy array in first argument(s) is not supported yet!")
     1857
     1858        result = numpy.zeros(args[-1].shape).tolist()
     1859        if isinstance(args[-1][0],int):
     1860            for k in xrange(len(args[-1])):
     1861                result[k] = laguerre(args[0],ZZ(args[-1][k]))
     1862
     1863        if isinstance(args[-1][0],float):
     1864            for k in xrange(len(args[-1])):
     1865                result[k] = laguerre(args[0],RR(args[-1][k]))
     1866
     1867        if isinstance(args[-1][0],numpy.ndarray):
     1868            for k in xrange(len(args[-1])):
     1869                result[k] = laguerre(args[0],args[-1][k])
     1870
     1871        return numpy.array(result)
     1872
     1873    def _derivative_(self,*args,**kwds):
     1874        """return the derivative of laguerre in
     1875        form of the Laguerre Polynomial.
     1876        EXAMPLES::
     1877        sage: n = var('n')
     1878        sage: derivative(laguerre(3,x),x)
     1879        -1/2*x^2 + 3*x - 3
     1880        sage: derivative(laguerre(n,x),x)
     1881        -(n*laguerre(n - 1, x) - n*laguerre(n, x))/x
     1882        """
     1883        diff_param = kwds['diff_param']
     1884        if diff_param == 0:
     1885            raise NotImplementedError(
     1886            "Derivative w.r.t. to the index is not supported, \
     1887             yet, and perhaps never will be...")
     1888        else:
     1889            return (args[0]*laguerre(args[0],args[-1])-args[0]*\
     1890                    laguerre(args[0]-1,args[-1]))/args[1]
     1891
     1892laguerre = Func_laguerre()
     1893
     1894class Func_legendre_P(OrthogonalPolynomial):
    5741895    """
    575     Returns the Legendre polynomial of the first kind for integers
    576     `n > -1`.
     1896    Returns the Legendre polynomial of the first kind..
    5771897   
    5781898    REFERENCE:
    5791899
     
    5881908        1.67750000000000
    5891909        sage: legendre_P(3, 1 + I)
    5901910        7/2*I - 13/2
    591         sage: legendre_P(3, MatrixSpace(ZZ, 2)([1, 2, -4, 7]))
    592         [-179  242]
    593         [-484  547]
     1911        sage: legendre_P(3, MatrixSpace(ZZ, 2)([1, 2, -4, 7]))
     1912        Traceback (most recent call last):
     1913        ...
     1914        TypeError: cannot coerce arguments: no canonical coercion from Full MatrixSpace of 2 by 2 dense matrices over Integer Ring to Symbolic Ring
    5941915        sage: legendre_P(3, GF(11)(5))
    5951916        8
    5961917    """
    597     _init()
    598     return sage_eval(maxima.eval('legendre_p(%s,x)'%ZZ(n)), locals={'x':x})
     1918    def __init__(self):
     1919        OrthogonalPolynomial.__init__(self,"legendre_P",nargs = 2,
     1920conversions =dict(maxima='legendre_p',mathematica='LegendreP'))
    5991921
    600 def legendre_Q(n,x):
     1922    def _clenshaw_method_(self,*args):
     1923        """
     1924        Clenshaw method for legendre_P (means use the recursion...)
     1925        This is much faster for numerical evaluation than maxima!
     1926        See A.S. 227 (p. 782) for details for the recurions.
     1927        Warning: The clanshaw method for the Legendre Polynomials
     1928        should only used for exact data types, when high orders are
     1929        used, due to weak instabilities of the recursion!
     1930        """
     1931        k = args[0]
     1932        x = args[-1]
     1933       
     1934        if k == 0:
     1935            return 1
     1936        elif k == 1:
     1937            return x
     1938        else:
     1939            help1 = 1
     1940            help2 = x
     1941            if is_Expression(x):
     1942                #raise NotImplementedError("Maxima is faster here...")
     1943                help1 = ZZ(2**k) #Workarround because of segmentation fault...
     1944                help3 = 0
     1945                for j in xrange(0,floor(k/2)+1):
     1946                    help3 = help3 + (-1)**j*x**(k-2*j)*binomial(k,j)*\
     1947                            binomial(2*(k-j),k)
     1948               
     1949                help3 = help3/help1
     1950            else:
     1951                for j in xrange(1,k):
     1952                    help3 = (2*j+1)*x*help2 - j*help1
     1953                    help3 = help3/(j+1)
     1954                    help1 = help2
     1955                    help2 = help3
     1956               
     1957            return help3   
     1958
     1959    def _maxima_init_evaled_(self, *args):
     1960        n = args[0]
     1961        x = args[1]
     1962        return sage_eval(maxima.eval('legendre_p(%s,x)'%ZZ(n)),\
     1963                         locals={'x':x})
     1964
     1965    def _evalf_(self, *args,**kwds):
     1966        """
     1967        Evals legendre_P
     1968        numerically with mpmath.
     1969        EXAMPLES::
     1970            sage: legendre_P(10,3).n(75)
     1971            8.097453000000000000000e6
     1972        """
     1973        try:
     1974            step_parent = kwds['parent']
     1975        except KeyError:
     1976            step_parent = parent(args[-1])
     1977
     1978        import sage.libs.mpmath.all as mpmath
     1979
     1980        try:
     1981            precision = step_parent.prec()
     1982        except AttributeError:
     1983            precision = mpmath.mp.prec
     1984
     1985        return mpmath.call(
     1986            mpmath.legendre,args[0],args[-1],prec = precision)
     1987   
     1988
     1989    def _eval_special_values_(self,*args):
     1990        """
     1991        Special values known.
     1992        EXAMPLES:
     1993           
     1994            sage: var('n')
     1995            n
     1996            sage: legendre_P(n,1)
     1997            1
     1998            sage: legendre_P(n,-1)
     1999            (-1)^n
     2000        """
     2001        if args[-1] == 1:
     2002            return 1
     2003       
     2004        if args[-1] == -1:
     2005            return (-1)**args[0]
     2006       
     2007        if (args[-1] == 0):
     2008            try:
     2009                return (1+(-1)**args[0])/2*binomial(args[0],args[0]/2)/4**(args[0]/2)
     2010            except TypeError:
     2011                pass
     2012
     2013        raise ValueError("Value not found")
     2014
     2015    def _eval_numpy_(self, *args):
     2016        #TODO: numpy_eval with help of the new scipy version!!!!
     2017        #Reason scipy 0.8 supports stable and fast numerical evaluation
     2018        #of ortho polys
     2019        #Now this only evaluates the array pointwise, and only the first one...
     2020        #if isinstance(arg[0], numpy.ndarray).
     2021        #This is a hack to provide compability!
     2022        """
     2023        EXAMPLES::
     2024        sage: import numpy
     2025        sage: z = numpy.array([1,2])
     2026        sage: z2 = numpy.array([[1,2],[1,2]])
     2027        sage: z3 = numpy.array([1,2,3.])
     2028        sage: legendre_P(1,z)
     2029        array([1, 2])
     2030        sage: legendre_P(1,z2)
     2031        array([[1, 2],
     2032              [1, 2]])
     2033        sage: legendre_P(1,z3)
     2034        array([ 1.,  2.,  3.])
     2035
     2036        """
     2037        if isinstance(args[0],numpy.ndarray):
     2038            raise NotImplementedError(
     2039           "Support for numpy array in first argument(s) is not supported yet!")
     2040
     2041        result = numpy.zeros(args[-1].shape).tolist()
     2042        if isinstance(args[-1][0],int):
     2043            for k in xrange(len(args[-1])):
     2044                result[k] = legendre_P(args[0],ZZ(args[-1][k]))
     2045
     2046        if isinstance(args[-1][0],float):
     2047            for k in xrange(len(args[-1])):
     2048                result[k] = legendre_P(args[0],RR(args[-1][k]))
     2049
     2050        if isinstance(args[-1][0],numpy.ndarray):
     2051            for k in xrange(len(args[-1])):
     2052                result[k] = legendre_P(args[0],args[-1][k])
     2053
     2054        return numpy.array(result)
     2055
     2056    def _derivative_(self,*args,**kwds):
     2057        """return the derivative of legendre_P in
     2058        form of the Legendre Polynomial.
     2059        EXAMPLES::
     2060        sage: n = var('n')
     2061        sage: derivative(legendre_P(n,x),x)
     2062        (n*x*legendre_P(n, x) - n*legendre_P(n - 1, x))/(x^2 - 1)
     2063        sage: derivative(legendre_P(3,x),x)
     2064        15/2*x^2 - 3/2
     2065        """
     2066        diff_param = kwds['diff_param']
     2067        if diff_param == 0:
     2068            raise NotImplementedError(
     2069"Derivative w.r.t. to the index is not supported, yet, and perhaps never will be...")
     2070        else:
     2071            return (args[0]*legendre_P(args[0]-1,args[-1])-args[0]*args[-1]*\
     2072                    legendre_P(args[0],args[-1]))/(1-args[-1]**2)
     2073       
     2074
     2075legendre_P = Func_legendre_P()     
     2076
     2077class Func_legendre_Q(OrthogonalPolynomial):
    6012078    """
    602     Returns the Legendre function of the second kind for integers
    603     `n>-1`.
     2079    Return the chebyshev function of the second kind.
    6042080   
    605     Computed using Maxima.
     2081    REFERENCE:
     2082   
     2083    - A.S. 8 (p. 332)
    6062084   
    6072085    EXAMPLES::
    6082086   
     
    6142092        sage: legendre_Q(4, 2)
    6152093        443/16*I*pi + 443/16*log(3) - 365/12
    6162094        sage: legendre_Q(4, 2.0)
    617         0.00116107583162324 + 86.9828465962674*I
     2095        0.00116107583162041 + 86.9828465962674*I
     2096
    6182097    """
    619     _init()
    620     return sage_eval(maxima.eval('legendre_q(%s,x)'%ZZ(n)), locals={'x':x})
     2098    def __init__(self):
     2099        OrthogonalPolynomial.__init__(self,"legendre_Q",nargs = 2,
     2100conversions =dict(maxima='legendre_q',mathematica='LegendreQ'))
    6212101
    622 def ultraspherical(n,a,x):
    623     """
    624     Returns the ultraspherical (or Gegenbauer) polynomial for integers
    625     `n > -1`.
     2102    def _maxima_init_evaled_(self, *args):
     2103        """
     2104        Maxima seems just fine for legendre Q. So we use it here!
     2105        """
     2106        n = args[0]
     2107        x = args[1]
     2108        return sage_eval(maxima.eval('legendre_q(%s,x)'%ZZ(n)),\
     2109                         locals={'x':x})
     2110
     2111    def _clenshaw_method_(self,*args):
     2112        """
     2113        Clenshaw method for legendre_q (means use the recursion...)
     2114        This is much faster for numerical evaluation than maxima!
     2115        See A.S. 8.5.3 (p. 334) for details for the recurions.
     2116        Warning: The clanshaw method for the Legendre fUNCTIONS
     2117        should only used for exact data types, when high orders are
     2118        used, due to weak instabilities of the recursion!
     2119        """
     2120        raise NotImplementedError("Function not ready yet...")
     2121       
     2122        k = args[0]
     2123        x = args[-1]
     2124       
     2125        if k == 0:
     2126            return ln((1+x)/(1-x))/2
     2127        elif k == 1:
     2128            return x/2*ln((1+x)/(1-x))-1
     2129        else:
     2130            if is_Expression(x):
     2131                raise NotImplementedError("Maxima works fine here!")
     2132                #it seems that the old method just works fine here...
     2133                #raise NotImplementedError("clenshaw does not work well...")
     2134            else:
     2135                help1 = ln((1+x)/(1-x))/2
     2136                help2 = x/2*ln((1+x)/(1-x))-1
     2137
     2138                for j in xrange(1,k):
     2139                    help3 = (2*j+1)*x*help2 - j*help1
     2140                    help3 = help3/(j+1)
     2141                    help1 = help2
     2142                    help2 = help3
     2143               
     2144            return help3   
     2145
     2146    def _eval_special_values_(self,*args):
     2147        """
     2148        Special values known.
     2149        EXAMPLES:
     2150           
     2151            sage: var('n')
     2152            n
     2153            sage: legendre_Q(n,0)
     2154            -1/2*sqrt(pi)*gamma(1/2*n + 1/2)*sin(1/2*pi*n)/gamma(1/2*n + 1)
     2155        """
     2156        if args[-1] == 1:
     2157            return NaN
     2158       
     2159        if args[-1] == -1:
     2160            return NaN
     2161       
     2162        if (args[-1] == 0):
     2163            if is_Expression(args[0]):
     2164                try:
     2165                    return -(sqrt(SR.pi()))/2*sin(SR.pi()/2*args[0])*\
     2166                           gamma((args[0]+1)/2)/gamma(args[0]/2 + 1)
     2167                except TypeError:
     2168                    pass
     2169            else:
     2170                return -(sqrt(math.pi))/2*sin(math.pi/2*args[0])*\
     2171                       gamma((args[0]+1)/2)/gamma(args[0]/2. + 1)
     2172
     2173        raise ValueError("Value not found")
     2174
     2175    def _evalf_(self, *args,**kwds):
     2176        """
     2177        Evals legendre_Q
     2178        numerically with mpmath.
     2179        EXAMPLES::
     2180            sage: legendre_Q(10,3).n(75)
     2181            2.079454941572578263731e-9 + 1.271944942879431601408e7*I
     2182        """
     2183
     2184        try:
     2185            step_parent = kwds['parent']
     2186        except KeyError:
     2187            step_parent = parent(args[-1])
     2188
     2189        import sage.libs.mpmath.all as mpmath
     2190
     2191        try:
     2192            precision = step_parent.prec()
     2193        except AttributeError:
     2194            precision = mpmath.mp.prec
     2195
     2196        return conjugate(
     2197            mpmath.call(mpmath.legenq,args[0],0,args[-1],prec = precision))
     2198        #it seems that mpmath uses here a different branch of the logarithm
     2199
     2200    def _eval_numpy_(self, *args):
     2201        #TODO: numpy_eval with help of the new scipy version!!!!
     2202        #Reason scipy 0.8 supports stable and fast numerical evaluation
     2203        #of ortho polys
     2204        #Now this only evaluates the array pointwise, and only the first one...
     2205        #if isinstance(arg[0], numpy.ndarray). This is a hack to provide compability!
     2206        """
     2207        EXAMPLES::
     2208        sage: import numpy
     2209        sage: z = numpy.array([1,2])
     2210        sage: z2 = numpy.array([[1,2],[1,2]])
     2211        sage: z3 = numpy.array([1,2,3.])
     2212        sage: legendre_Q(1,z/5.)
     2213        array([-0.95945349, -0.83054043])
     2214        sage: legendre_Q(1,z2/5.)
     2215        array([[-0.95945349, -0.83054043],
     2216               [-0.95945349, -0.83054043]])
     2217        sage: legendre_Q(1,z3/5.)
     2218        array([-0.95945349, -0.83054043, -0.58411169])
     2219
     2220        """
     2221        if isinstance(args[0],numpy.ndarray):
     2222            raise NotImplementedError(
     2223           "Support for numpy array in first argument(s) is not supported yet!")
     2224
     2225        result = numpy.zeros(args[-1].shape).tolist()
     2226        if isinstance(args[-1][0],int):
     2227            for k in xrange(len(args[-1])):
     2228                result[k] = legendre_Q(args[0],ZZ(args[-1][k]))
     2229
     2230        if isinstance(args[-1][0],float):
     2231            for k in xrange(len(args[-1])):
     2232                result[k] = legendre_Q(args[0],RR(args[-1][k]))
     2233
     2234        if isinstance(args[-1][0],numpy.ndarray):
     2235            for k in xrange(len(args[-1])):
     2236                result[k] = legendre_Q(args[0],args[-1][k])
     2237
     2238        return numpy.array(result)
     2239
     2240    def _derivative_(self,*args,**kwds):
     2241        """return the derivative of legendre_Q in
     2242        form of the Legendre Function.
     2243        EXAMPLES::
     2244        n = var('n')
     2245        derivative(legendre_Q(n,x),x)
     2246        (n*x*legendre_Q(n, x) - n*legendre_Q(n - 1, x))/(x^2 - 1)
     2247        sage: derivative(legendre_Q(3,x),x)
     2248        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))
     2249        """
     2250        diff_param = kwds['diff_param']
     2251        if diff_param == 0:
     2252            raise NotImplementedError(
     2253"Derivative w.r.t. to the index is not supported, yet, and perhaps never will be...")
     2254        else:
     2255            return (args[0]*args[-1]*legendre_Q(args[0],args[-1])-args[0]*\
     2256                    legendre_Q(args[0]-1,args[-1]))/(args[-1]**2-1)
     2257
     2258legendre_Q = Func_legendre_Q()
     2259
     2260
     2261class Func_gen_legendre_P(OrthogonalPolynomial):
     2262
     2263    def __init__(self):
     2264        OrthogonalPolynomial.__init__(self,"gen_legendre_P",nargs = 3,
     2265conversions =dict(maxima='assoc_legendre_p',mathematica='LegendreP'))
     2266
     2267    def _evalf_(self, *args,**kwds):
     2268        """
     2269        Evals gen_legendre_P
     2270        numerically with mpmath.
     2271        EXAMPLES::
     2272            sage: gen_legendre_P(10,2,3).n(75)
     2273            -7.194963600000000000000e8
     2274        """
     2275
     2276        try:
     2277            step_parent = kwds['parent']
     2278        except KeyError:
     2279            step_parent = parent(args[-1])
     2280
     2281        import sage.libs.mpmath.all as mpmath
     2282
     2283        try:
     2284            precision = step_parent.prec()
     2285        except AttributeError:
     2286            precision = mpmath.mp.prec
     2287
     2288        return mpmath.call(
     2289            mpmath.legenp,args[0],args[1],args[-1],prec = precision)
     2290
     2291    def _eval_special_values_(self,*args):
     2292        """
     2293        Special values known.
     2294        EXAMPLES:
     2295           
     2296            sage: n, m = var('n m')
     2297            sage: gen_legendre_P(n,m,0)
     2298            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))
     2299        """
     2300           
     2301        if args[1] == 0:
     2302            return legendre_P(args[0],args[-1])
     2303
     2304        if (args[-1] == 0):
     2305            if is_Expression(args[0]):
     2306                try:
     2307                    return cos(SR.pi()/2*(args[0]+args[1]))/(sqrt(SR.pi()))*\
     2308                           gamma((args[0]+args[1]+1)/2)/\
     2309                           gamma((args[0]-args[1])/2 + 1)*2**(args[1])
     2310                except TypeError:
     2311                    pass
     2312            else:
     2313                return cos(math.pi/2*(args[0]+args[1]))/(sqrt(math.pi))*\
     2314                       gamma((args[0]+args[1]+1)/2)/\
     2315                       gamma((args[0]-args[1])/2. + 1)*2**args[1]
     2316
     2317        raise ValueError("Value not found")
     2318
     2319    def _maxima_init_evaled_(self, *args):
     2320        n = args[0]
     2321        m = args[1]
     2322        x = args[2]
     2323        if is_Expression(n) or is_Expression(m):
     2324            return None
     2325       
     2326        from sage.functions.all import sqrt
    6262327   
    627     Computed using Maxima.
     2328        if m.mod(2).is_zero() or m.is_one():
     2329            return sage_eval(maxima.eval('assoc_legendre_p(%s,%s,x)'\
     2330                                         %(ZZ(n),ZZ(m))), locals={'x':x})
     2331        else:
     2332            return sqrt(1-x**2)*(((n-m+1)*x*gen_legendre_P(n,m-1,x)-\
     2333                                  (n+m-1)*gen_legendre_P(n-1,m-1,x))/(1-x**2))
     2334
     2335
     2336    def _eval_numpy_(self, *args):
     2337        #TODO: numpy_eval with help of the new scipy version!!!!
     2338        #Reason scipy 0.8 supports stable and fast numerical evaluation
     2339        #of ortho polys
     2340        #Now this only evaluates the array pointwise, and only the first one...
     2341        #if isinstance(arg[0], numpy.ndarray).
     2342        #This is a hack to provide compability!
     2343        """
     2344        EXAMPLES::
     2345        sage: import numpy
     2346        sage: z = numpy.array([1,2])
     2347        sage: z2 = numpy.array([[1,2],[1,2]])
     2348        sage: z3 = numpy.array([1,2,3.])
     2349        sage: gen_legendre_P(1,1,z/5.)
     2350        array([-0.9797959 , -0.91651514])
     2351        sage: gen_legendre_P(1,1,z2/5.)
     2352        array([[-0.9797959 , -0.91651514],
     2353               [-0.9797959 , -0.91651514]])
     2354        sage: gen_legendre_P(1,1,z3/5.)
     2355        array([-0.9797959 , -0.91651514, -0.8       ])
     2356        """
     2357
     2358        if isinstance(args[0],numpy.ndarray) or isinstance(args[1],numpy.ndarray):
     2359            raise NotImplementedError(
     2360           "Support for numpy array in first argument(s) is not supported yet!")
     2361       
     2362        result = numpy.zeros(args[-1].shape).tolist()
     2363        if isinstance(args[-1][0],int):
     2364            for k in xrange(len(args[-1])):
     2365                result[k] = gen_legendre_P(args[0],args[1],ZZ(args[-1][k]))
     2366
     2367        if isinstance(args[-1][0],float):
     2368            for k in xrange(len(args[-1])):
     2369                result[k] = gen_legendre_P(args[0],args[1],RR(args[-1][k]))
     2370
     2371        if isinstance(args[-1][0],numpy.ndarray):
     2372            for k in xrange(len(args[-1])):
     2373                result[k] = gen_legendre_P(args[0],args[1],args[-1][k])
     2374
     2375        return numpy.array(result)
     2376     
     2377    def _derivative_(self,*args,**kwds):
     2378        """return the derivative of gen_legendre_P in
     2379        form of the Legendre Function.
     2380        EXAMPLES::
     2381        sage: n,m = var('n m')
     2382        sage: derivative(gen_legendre_P(n,m,x),x)
     2383        (n*x*gen_legendre_P(n, m, x) - (m + n)*gen_legendre_P(n - 1, m, x))/(x^2 - 1)
     2384        sage: derivative(gen_legendre_P(3,1,x),x)
     2385        -15*sqrt(-x^2 + 1)*x + 3/2*(5*(x - 1)^2 + 10*x - 6)*x/sqrt(-x^2 + 1)
     2386        """
     2387        diff_param = kwds['diff_param']
     2388        if diff_param in [0,1]:
     2389            raise NotImplementedError(
     2390            "Derivative w.r.t. to the index is not supported,\
     2391            yet, and perhaps never will be...")
     2392        else:
     2393            return (args[0]*args[-1]*gen_legendre_P(args[0],args[1],args[-1])\
     2394                    -(args[0]+args[1])*\
     2395                    gen_legendre_P(args[0]-1,args[1],args[-1]))/(args[-1]**2-1)
     2396
     2397gen_legendre_P = Func_gen_legendre_P()
     2398
     2399class Func_gen_legendre_Q(OrthogonalPolynomial):
     2400
     2401    def __init__(self):
     2402        OrthogonalPolynomial.__init__(self,"gen_legendre_Q",nargs = 3,
     2403conversions =dict(maxima='assoc_legendre_q',mathematica='LegendreQ'))
     2404
     2405    def _evalf_(self, *args,**kwds):
     2406        """
     2407        Evals gen_legendre_Q
     2408        numerically with mpmath.
     2409        EXAMPLES::
     2410            sage: gen_legendre_Q(10,2,3).n(75)
     2411            -2.773909528741569374688e-7 - 1.130182239430298584113e9*I
     2412        """
     2413
     2414        try:
     2415            step_parent = kwds['parent']
     2416        except KeyError:
     2417            step_parent = parent(args[-1])
     2418
     2419        import sage.libs.mpmath.all as mpmath
     2420
     2421        try:
     2422            precision = step_parent.prec()
     2423        except AttributeError:
     2424            precision = mpmath.mp.prec
     2425
     2426        return mpmath.call(
     2427            mpmath.legenq,args[0],args[1],args[-1],prec = precision)
     2428
     2429    def _eval_special_values_(self,*args):
     2430        """
     2431        Special values known.
     2432        EXAMPLES:
     2433           
     2434            sage: n, m = var('n m')
     2435            sage: gen_legendre_Q(n,m,0)
     2436            -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)
     2437        """
     2438           
     2439        if args[1] == 0:
     2440            return legendre_Q(args[0],args[-1])
     2441
     2442        if (args[-1] == 0):
     2443            if is_Expression(args[0]):
     2444                try:
     2445                    return -(sqrt(SR.pi()))*sin(SR.pi()/2*(args[0]+args[1]))\
     2446                           *gamma((args[0]+args[1]+1)/2)/\
     2447                           gamma((args[0]-args[1])/2 + 1)*2**(args[1]-1)
     2448                except TypeError:
     2449                    pass
     2450            else:
     2451                return -(sqrt(math.pi))/2*sin(math.pi/2*(args[0]+args[1]))\
     2452                       *gamma((args[0]+args[1]+1)/2)/\
     2453                       gamma((args[0]-args[1])/2. + 1)*2**args[1]
     2454
     2455        raise ValueError("Value not found")
     2456
     2457    def _maxima_init_evaled_(self, *args):
     2458        n = args[0]
     2459        m = args[1]
     2460        x = args[2]
     2461        if is_Expression(n) or is_Expression(m):
     2462            return None
     2463       
     2464        from sage.functions.all import sqrt
    6282465   
    629     REFERENCE:
     2466        if m <= n:
     2467            return sage_eval(maxima.eval('assoc_legendre_q(%s,%s,x)'\
     2468                                         %(ZZ(n),ZZ(m))), locals={'x':x})
     2469        if m == n + 1 or n == 0:
     2470            if m.mod(2).is_zero():
     2471                denom = (1 - x**2)**(m/2)
     2472            else:
     2473                denom = sqrt(1 - x**2)*(1 - x**2)**((m-1)/2)   
     2474                if m == n + 1:
     2475                    return (-1)**m*(m-1).factorial()*2**n/denom
     2476                else:
     2477                    return (-1)**m*(m-1).factorial()*((x+1)**m - (x-1)**m)/(2*denom)
     2478        else:
     2479            return ((n-m+1)*x*gen_legendre_Q(n,m-1,x)-\
     2480                    (n+m-1)*gen_legendre_Q(n-1,m-1,x))/sqrt(1-x**2)
    6302481
    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})
    6482482
    649 gegenbauer = ultraspherical
     2483    def _eval_numpy_(self, *args):
     2484        #TODO: numpy_eval with help of the new scipy version!!!!
     2485        #Reason scipy 0.8 supports stable and fast numerical evaluation
     2486        #of ortho polys
     2487        #Now this only evaluates the array pointwise, and only the first one...
     2488        #if isinstance(arg[0], numpy.ndarray).
     2489        #This is a hack to provide compability!
     2490        """
     2491        EXAMPLES::
     2492        sage: import numpy
     2493        sage: z = numpy.array([1,2])
     2494        sage: z2 = numpy.array([[1,2],[1,2]])
     2495        sage: z3 = numpy.array([1,2,3.])
     2496        sage: gen_legendre_Q(1,1,z/5.)
     2497        array([-0.40276067, -0.82471644])
     2498        sage: gen_legendre_Q(1,1,z2/5.)
     2499        array([[-0.40276067, -0.82471644],
     2500               [-0.40276067, -0.82471644]])
     2501        sage: gen_legendre_Q(1,1,z3/5.)
     2502        array([-0.40276067, -0.82471644, -1.30451774])
     2503        """
     2504
     2505        if isinstance(args[0],numpy.ndarray) or isinstance(args[1],numpy.ndarray):
     2506            raise NotImplementedError(
     2507           "Support for numpy array in first argument(s) is not supported yet!")
     2508       
     2509        result = numpy.zeros(args[-1].shape).tolist()
     2510        if isinstance(args[-1][0],int):
     2511            for k in xrange(len(args[-1])):
     2512                result[k] = gen_legendre_Q(args[0],args[1],ZZ(args[-1][k]))
     2513
     2514        if isinstance(args[-1][0],float):
     2515            for k in xrange(len(args[-1])):
     2516                result[k] = gen_legendre_Q(args[0],args[1],RR(args[-1][k]))
     2517
     2518        if isinstance(args[-1][0],numpy.ndarray):
     2519            for k in xrange(len(args[-1])):
     2520                result[k] = gen_legendre_Q(args[0],args[1],args[-1][k])
     2521
     2522        return numpy.array(result)
     2523     
     2524    def _derivative_(self,*args,**kwds):
     2525        """return the derivative of gen_legendre_Q in
     2526        form of the Legendre Function.
     2527        EXAMPLES::
     2528        sage: n,m = var('n m')
     2529        sage: derivative(gen_legendre_Q(n,m,x),x)
     2530        (n*x*gen_legendre_Q(n, m, x) - (m + n)*gen_legendre_Q(n - 1, m, x))/(x^2 - 1)
     2531        sage: derivative(gen_legendre_Q(0,1,x),x)
     2532        -x/(-x^2 + 1)^(3/2)
     2533        """
     2534        diff_param = kwds['diff_param']
     2535        if diff_param in [0,1]:
     2536            raise NotImplementedError("Derivative w.r.t. to the index \
     2537            is not supported, yet, and perhaps never will be...")
     2538        else:
     2539            return (args[0]*args[-1]*gen_legendre_Q(args[0],args[1],args[-1])\
     2540                    -(args[0]+args[1])*\
     2541                    gen_legendre_Q(args[0]-1,args[1],args[-1]))/(args[-1]**2-1)
     2542
     2543
     2544gen_legendre_Q = Func_gen_legendre_Q()
     2545
     2546
     2547
     2548
     2549
     2550
     2551
     2552
     2553
     2554
     2555
     2556