Ticket #4102: trac_symbolic_bessel.metaclass.2.patch

File trac_symbolic_bessel.metaclass.2.patch, 30.6 KB (added by benjaminfjones, 9 years ago)

work in progress making Bessel functions symbolic

  • sage/functions/all.py

    # HG changeset patch
    # User Benjamin Jones <benjaminfjones@gmail.com>
    # Date 1353303890 28800
    # Node ID a15b2a0d86d0c48ca8c33de6689b31a49148385c
    # Parent  d06cf4b2215d37d3a87a58f65ac53234502dd471
    Trac 4102: implement symbolic Bessel functions
    
    diff --git a/sage/functions/all.py b/sage/functions/all.py
    a b  
    88                   arctan2, atan2)
    99
    1010from hyperbolic import ( tanh, sinh, cosh, coth, sech, csch,
    11                          asinh, acosh, atanh, acoth, asech, acsch, 
     11                         asinh, acosh, atanh, acoth, asech, acsch,
    1212                         arcsinh, arccosh, arctanh, arccoth, arcsech, arccsch )
    1313
    1414reciprocal_trig_functions = {'sec': cos, 'csc': sin, 'cot': tan, 'sech': cosh, 'csch': sinh, 'coth': tanh}
     
    2626
    2727from transcendental import (zeta, zetaderiv, zeta_symmetric, dickman_rho)
    2828
    29 from special import (bessel_I, bessel_J, bessel_K, bessel_Y,
    30                      hypergeometric_U, Bessel,
     29from special import (Bessel, bessel_I, bessel_J, bessel_K, bessel_Y,
     30                     hypergeometric_U,
    3131                     spherical_bessel_J, spherical_bessel_Y,
    3232                     spherical_hankel1, spherical_hankel2,
    3333                     spherical_harmonic, jacobi,
     
    3636                     elliptic_f, elliptic_ec, elliptic_eu,
    3737                     elliptic_kc, elliptic_pi, elliptic_j,
    3838                     airy_ai, airy_bi)
    39                        
     39
    4040from orthogonal_polys import (chebyshev_T,
    4141                              chebyshev_U,
    4242                              gen_laguerre,
  • sage/functions/special.py

    diff --git a/sage/functions/special.py b/sage/functions/special.py
    a b  
    384384from sage.plot.plot import plot
    385385from sage.rings.real_mpfr import RealField
    386386from sage.rings.complex_field import ComplexField
    387 from sage.misc.sage_eval import sage_eval
    388387from sage.rings.all import ZZ, RR, RDF
    389388from sage.functions.other import real, imag, log_gamma
    390 from sage.symbolic.function import BuiltinFunction
     389from sage.symbolic.function import BuiltinFunction, is_inexact
     390from sage.symbolic.expression import Expression
     391from sage.structure.coerce import parent
     392import sage.structure.element
    391393from sage.calculus.calculus import maxima
     394from sage.libs.mpmath import utils as mpmath_utils
     395from sage.misc.sage_eval import sage_eval
    392396
    393397_done = False
    394398def _init():
     
    555559
    556560    EXAMPLES::
    557561
    558         sage: n(bessel_J(3,10,"maxima"))
    559         0.0583793793051...
     562        sage: elliptic_e(0.5, 0.1)
     563        0.498011394499
    560564        sage: spherical_hankel2(2,i)
    561565        -e
    562566    """
     
    570574
    571575            TESTS::
    572576
    573                 sage: n(bessel_J(3,10,"maxima"))
    574                 0.0583793793051...
     577                sage: elliptic_e(0.5, 0.1)
     578                0.498011394499
    575579                sage: spherical_hankel2(2,x)
    576580                (-I*x^2 - 3*x + 3*I)*e^(-I*x)/x^3
    577581            """
     
    653657   return RDF(meval("airy_bi(%s)"%RDF(x)))
    654658
    655659
    656 def bessel_I(nu,z,algorithm = "pari",prec=53):
    657     r"""
    658     Implements the "I-Bessel function", or "modified Bessel function,
    659     1st kind", with index (or "order") nu and argument z.
    660    
    661     INPUT:
    662    
    663    
    664     -  ``nu`` - a real (or complex, for pari) number
    665    
    666     -  ``z`` - a real (positive) algorithm - "pari" or
    667        "maxima" or "scipy" prec - real precision (for PARI only)
    668    
    669    
    670     DEFINITION::
    671    
    672             Maxima:
    673                              inf
    674                             ====   - nu - 2 k  nu + 2 k
    675                             \     2          z
    676                              >    -------------------
    677                             /     k! Gamma(nu + k + 1)
    678                             ====
    679                             k = 0
    680        
    681             PARI:
    682            
    683                              inf
    684                             ====   - 2 k  2 k
    685                             \     2      z    Gamma(nu + 1)
    686                              >    -----------------------
    687                             /       k! Gamma(nu + k + 1)
    688                             ====
    689                             k = 0
    690        
    691            
    692    
    693     Sometimes ``bessel_I(nu,z)`` is denoted
    694     ``I_nu(z)`` in the literature.
    695    
    696     .. warning::
     660def Bessel(*args, **kwds):
     661    """
     662    A class factory that produces symbolic I, J, K, and Y Bessel functions. There
     663    are several ways to call this function:
    697664
    698        In Maxima (the manual says) i0 is deprecated but
    699        ``bessel_i(0,*)`` is broken. (Was fixed in recent CVS patch
    700        though.)
    701    
    702     EXAMPLES::
    703    
    704         sage: bessel_I(1,1,"pari",500)
    705         0.565159103992485027207696027609863307328899621621092009480294489479255640964371134092664997766814410064677886055526302676857637684917179812041131208121
    706         sage: bessel_I(1,1)
    707         0.565159103992485
    708         sage: bessel_I(2,1.1,"maxima") 
    709         0.16708949925104...
    710         sage: bessel_I(0,1.1,"maxima")
    711         1.32616018371265...
    712         sage: bessel_I(0,1,"maxima")   
    713         1.2660658777520...
    714         sage: bessel_I(1,1,"scipy")
    715         0.565159103992...
     665        - ``Bessel(order, type)``
     666        - ``Bessel(order)`` -- type defaults to 'J'
     667        - ``Bessel(order, typ=T)``
     668        - ``Bessel(typ=T)`` -- order is unspecified, this is a 2-parameter function
     669        - ``Bessel()`` -- order is unspecified, type is 'J'
    716670
    717     Check whether the return value is real whenever the argument is real (#10251)::
    718    
    719         sage: bessel_I(5, 1.5, algorithm='scipy') in RR
    720         True
    721        
    722     """
    723     if algorithm=="pari":
    724         from sage.libs.pari.all import pari
    725         try:
    726             R = RealField(prec)
    727             nu = R(nu)
    728             z = R(z)
    729         except TypeError:
    730             C = ComplexField(prec)
    731             nu = C(nu)
    732             z = C(z)
    733             K = C
    734         K = z.parent()
    735         return K(pari(nu).besseli(z, precision=prec))
    736     elif algorithm=="scipy":
    737         if prec != 53:
    738             raise ValueError, "for the scipy algorithm the precision must be 53"
    739         import scipy.special
    740         ans = str(scipy.special.iv(float(nu),complex(real(z),imag(z))))
    741         ans = ans.replace("(","")
    742         ans = ans.replace(")","")
    743         ans = ans.replace("j","*I")
    744         ans = sage_eval(ans)
    745         return real(ans) if z in RR else ans # Return real value when arg is real
    746     elif algorithm == "maxima":
    747         if prec != 53:
    748             raise ValueError, "for the maxima algorithm the precision must be 53"
    749         return sage_eval(maxima.eval("bessel_i(%s,%s)"%(float(nu),float(z))))
    750     else:
    751         raise ValueError, "unknown algorithm '%s'"%algorithm
    752        
    753 def bessel_J(nu,z,algorithm="pari",prec=53):
    754     r"""
    755     Return value of the "J-Bessel function", or "Bessel function, 1st
    756     kind", with index (or "order") nu and argument z.
    757    
    758     ::
    759    
    760             Defn:
    761             Maxima:
    762                              inf
    763                             ====          - nu - 2 k  nu + 2 k
    764                             \     (-1)^k 2           z
    765                              >    -------------------------
    766                             /        k! Gamma(nu + k + 1)
    767                             ====
    768                             k = 0
    769        
    770             PARI:
    771            
    772                              inf
    773                             ====          - 2k    2k
    774                             \     (-1)^k 2      z    Gamma(nu + 1)
    775                              >    ----------------------------
    776                             /         k! Gamma(nu + k + 1)
    777                             ====
    778                             k = 0
    779            
    780    
    781     Sometimes bessel_J(nu,z) is denoted J_nu(z) in the literature.
    782    
    783     .. warning::
     671    where `order` can be any integer and T must be one of the strings 'I', 'J',
     672    'K', or 'Y'.
    784673
    785        Inaccurate for small values of z.
    786    
    787     EXAMPLES::
    788    
    789         sage: bessel_J(2,1.1)
    790         0.136564153956658
    791         sage: bessel_J(0,1.1)
    792         0.719622018527511
    793         sage: bessel_J(0,1)
    794         0.765197686557967
    795         sage: bessel_J(0,0)
    796         1.00000000000000
    797         sage: bessel_J(0.1,0.1)
    798         0.777264368097005
    799    
    800     We check consistency of PARI and Maxima::
    801    
    802         sage: n(bessel_J(3,10,"maxima"))
    803         0.0583793793051...
    804         sage: n(bessel_J(3,10,"pari")) 
    805         0.0583793793051868
    806         sage: bessel_J(3,10,"scipy")
    807         0.0583793793052...
     674    See the EXAMPLES below.
    808675
    809     Check whether the return value is real whenever the argument is real (#10251)::                                                                                                                                                           
    810         sage: bessel_J(5, 1.5, algorithm='scipy') in RR                                                                     
    811         True
    812     """
    813    
    814     if algorithm=="pari":
    815         from sage.libs.pari.all import pari
    816         try:
    817             R = RealField(prec)
    818             nu = R(nu)
    819             z = R(z)
    820         except TypeError:
    821             C = ComplexField(prec)
    822             nu = C(nu)
    823             z = C(z)
    824             K = C
    825         if nu == 0:
    826             nu = ZZ(0)
    827         K = z.parent()
    828         return K(pari(nu).besselj(z, precision=prec))
    829     elif algorithm=="scipy":
    830         if prec != 53:
    831             raise ValueError, "for the scipy algorithm the precision must be 53"
    832         import scipy.special
    833         ans = str(scipy.special.jv(float(nu),complex(real(z),imag(z))))
    834         ans = ans.replace("(","")
    835         ans = ans.replace(")","")
    836         ans = ans.replace("j","*I")
    837         ans = sage_eval(ans)
    838         return real(ans) if z in RR else ans
    839     elif algorithm == "maxima":
    840         if prec != 53:
    841             raise ValueError, "for the maxima algorithm the precision must be 53"
    842         return maxima_function("bessel_j")(nu, z)
    843     else:
    844         raise ValueError, "unknown algorithm '%s'"%algorithm
     676    EXAMPLES:
    845677
    846 def bessel_K(nu,z,algorithm="pari",prec=53):
    847     r"""
    848     Implements the "K-Bessel function", or "modified Bessel function,
    849     2nd kind", with index (or "order") nu and argument z. Defn::
    850    
    851                     pi*(bessel_I(-nu, z) - bessel_I(nu, z))
    852                    ----------------------------------------
    853                                 2*sin(pi*nu)
    854            
    855    
    856     if nu is not an integer and by taking a limit otherwise.
    857    
    858     Sometimes bessel_K(nu,z) is denoted K_nu(z) in the literature. In
    859     PARI, nu can be complex and z must be real and positive.
    860    
    861     EXAMPLES::
    862    
    863         sage: bessel_K(3,2,"scipy")
    864         0.64738539094...
    865         sage: bessel_K(3,2)
    866         0.64738539094...
    867         sage: bessel_K(1,1)
    868         0.60190723019...
    869         sage: bessel_K(1,1,"pari",10)
    870         0.60
    871         sage: bessel_K(1,1,"pari",100)
    872         0.60190723019723457473754000154
     678    Construction of Bessel functions with various orders and types::
    873679
    874     TESTS::
     680        sage: Bessel()
     681        bessel_J
     682        sage: Bessel(1)
     683        bessel_{1,J}
     684        sage: Bessel(1, 'Y')
     685        bessel_{1,Y}
     686        sage: Bessel(-2, 'Y')
     687        bessel_{-2,Y}
     688        sage: Bessel(typ='K')
     689        bessel_K
     690        sage: Bessel(0, typ='I')
     691        bessel_{0,I}
    875692
    876         sage: bessel_K(2,1.1, algorithm="maxima")
    877         Traceback (most recent call last):
    878         ...
    879         NotImplementedError: The K-Bessel function is only implemented for the pari and scipy algorithms
     693        sage: f = Bessel(4, 'Y')
     694        sage: f.get_type()
     695        'Y'
     696        sage: f.get_order()
     697        4
     698        sage: f.get_system()
     699        'mpmath'
    880700
    881         Check whether the return value is real whenever the argument is real (#10251)::
     701    Evaluation::
    882702
    883         sage: bessel_K(5, 1.5, algorithm='scipy') in RR
     703        sage: f = Bessel(1)
     704        sage: f(3.0)
     705        0.339058958525936
     706        sage: f(3).n(digits=50)
     707        0.33905895852593645892551459720647889697308041819801
     708
     709        sage: g = Bessel(typ='J')
     710        sage: g(1,3)
     711        bessel_J(1, 3)
     712        sage: g(2, 3+I).n()
     713        0.634160370148554 + 0.0253384000032695*I
     714        sage: abs(numerical_integral(1/pi*cos(3*sin(x)), 0.0, pi)[0] - Bessel(0, 'J')(3.0)) < 1e-15
    884715        True
    885716
     717    Symbolic calculus::
     718
     719        sage: f(x) = Bessel(0, 'J')(x)
     720        sage: derivative(f, x)
     721        x |--> bessel_{-1,J}(x)
     722        sage: derivative(f, x, x)
     723        x |--> bessel_{-1,J}(x)/x + bessel_{-2,J}(x)
     724
     725    Conversion to other systems::
     726
     727        sage: x,y = var('x,y')
     728        sage: f = maxima(Bessel(typ='K')(x,y))
     729        sage: f.derivative('x')
     730        %pi*csc(%pi*x)*('diff(bessel_i(-x,y),x,1)-'diff(bessel_i(x,y),x,1))/2-%pi*bessel_k(x,y)*cot(%pi*x)
     731        sage: f.derivative('y')
     732        -(bessel_k(x+1,y)+bessel_k(x-1,y))/2
     733
     734    Plotting::
     735
     736        sage: f(x) = Bessel(0)(x); f
     737        x |--> bessel_{0,J}(x)
     738        sage: plot(f, (x, 1, 10))
     739
     740        sage: plot([ Bessel(i, 'J') for i in range(5) ], 2, 10)
     741
     742        sage: G = Graphics()
     743        sage: for i in range(5): G += plot(Bessel(i), 0, 4*pi, rgbcolor=hue(sin(pi*i/10)))
     744        sage: show(G)
     745
     746    A recreation of Abramowitz and Stegun Figure 9.1::
     747
     748        sage: G  = plot(Bessel(0, 'J'), 0, 15, color='black')
     749        sage: G += plot(Bessel(0, 'Y'), 0, 15, color='black')
     750        sage: G += plot(Bessel(1, 'J'), 0, 15, color='black', linestyle='dotted')
     751        sage: G += plot(Bessel(1, 'Y'), 0, 15, color='black', linestyle='dotted')
     752        sage: show(G, ymin=-1, ymax=1)
     753
    886754    """
    887     if algorithm=="scipy":
    888         if prec != 53:
    889             raise ValueError, "for the scipy algorithm the precision must be 53"
    890         import scipy.special
    891         ans = str(scipy.special.kv(float(nu),float(z)))
    892         ans = ans.replace("(","")
    893         ans = ans.replace(")","")
    894         ans = ans.replace("j","*I")
    895         ans = sage_eval(ans)
    896         return real(ans) if z in RR else ans
    897     elif algorithm == 'pari':
    898         from sage.libs.pari.all import pari
    899         try:
    900             R = RealField(prec)
    901             nu = R(nu)
    902             z = R(z)
    903         except TypeError:
    904             C = ComplexField(prec)
    905             nu = C(nu)
    906             z = C(z)
    907             K = C
    908         K = z.parent()
    909         return K(pari(nu).besselk(z, precision=prec))
    910     elif algorithm == 'maxima':
    911         raise NotImplementedError, "The K-Bessel function is only implemented for the pari and scipy algorithms"
     755
     756    # conversions dictionary for the Bessel function
     757    _conversions_dict = { 'I': { 'maxima': 'bessel_i', 'mathematica': 'BesselI' }
     758                        , 'J': { 'maxima': 'bessel_j', 'mathematica': 'BesselJ' }
     759                        , 'K': { 'maxima': 'bessel_k', 'mathematica': 'BesselK' }
     760                        , 'Y': { 'maxima': 'bessel_y', 'mathematica': 'BesselY' } }
     761
     762    class_attrs = {} # storage for attributes of the class to be created
     763
     764    # determine the order and type of function from the arguments and keywords
     765    if 'typ' in kwds:
     766        class_attrs['_type'] = kwds['typ']
    912767    else:
    913         raise ValueError, "unknown algorithm '%s'"%algorithm
    914    
     768        class_attrs['_type'] = 'J'
    915769
    916 def bessel_Y(nu,z,algorithm="maxima", prec=53):
    917     r"""
    918     Implements the "Y-Bessel function", or "Bessel function of the 2nd
    919     kind", with index (or "order") nu and argument z.
    920    
    921     .. note::
     770    if len(args) == 0: # no order specified
     771        class_attrs['_order'] = None
     772        class_attrs['_nargs'] = 2
     773        class_attrs['_name']= "bessel_%s" % class_attrs['_type']
     774        class_attrs['_latex_name'] = r"bessel\_%s" % class_attrs['_type']
     775    elif len(args) == 1: # order is specified
     776        class_attrs['_order'] = args[0]
     777        class_attrs['_nargs'] = 1
     778        class_attrs['_name'] = "bessel_{%s,%s}" % (class_attrs['_order'], class_attrs['_type'])
     779        class_attrs['_latex_name'] = r"bessel\_\{%s,%s\}" % (class_attrs['_order'], class_attrs['_type'])
     780    elif len(args) == 2: # both order and type are positional arguments
     781        class_attrs['_order'] = args[0]
     782        class_attrs['_type'] = args[1]
     783        class_attrs['_nargs'] = 1
     784        class_attrs['_name'] = "bessel_{%s,%s}" % (class_attrs['_order'], class_attrs['_type'])
     785        class_attrs['_latex_name'] = r"bessel\_\{%s,%s\}" % (class_attrs['_order'], class_attrs['_type'])
     786    else:
     787        raise TypeError("at most two positional arguments may be specified, "
     788                       +"see the docstring for Bessel")
    922789
    923        Currently only prec=53 is supported.
    924    
    925     Defn::
    926    
    927                     cos(pi n)*bessel_J(nu, z) - bessel_J(-nu, z)
    928                    -------------------------------------------------
    929                                      sin(nu*pi)
    930    
    931     if nu is not an integer and by taking a limit otherwise.
    932    
    933     Sometimes bessel_Y(n,z) is denoted Y_n(z) in the literature.
    934    
    935     This is computed using Maxima by default.
    936    
    937     EXAMPLES::
    938    
    939         sage: bessel_Y(2,1.1,"scipy")
    940         -1.4314714939...
    941         sage: bessel_Y(2,1.1)   
    942         -1.4314714939590...
    943         sage: bessel_Y(3.001,2.1)
    944         -1.0299574976424...
     790    # check the function type
     791    if not (class_attrs['_type'] in ['I', 'J', 'K', 'Y']):
     792        raise ValueError, "type must be one of I, J, K, Y"
    945793
    946     TESTS::
     794    # record the numerical evaluation system
     795    if 'algorithm' in kwds:
     796        class_attrs['_system'] = kwds['algorithm']
     797    else:
     798        class_attrs['_system'] = 'mpmath'
    947799
    948         sage: bessel_Y(2,1.1, algorithm="pari")
    949         Traceback (most recent call last):
    950         ...
    951         NotImplementedError: The Y-Bessel function is only implemented for the maxima and scipy algorithms
    952     """
    953     if algorithm=="scipy":
    954         if prec != 53:
    955             raise ValueError, "for the scipy algorithm the precision must be 53"
    956         import scipy.special
    957         ans = str(scipy.special.yv(float(nu),complex(real(z),imag(z))))
    958         ans = ans.replace("(","")
    959         ans = ans.replace(")","")
    960         ans = ans.replace("j","*I")
    961         ans = sage_eval(ans)
    962         return real(ans) if z in RR else ans
    963     elif algorithm == "maxima":
    964         if prec != 53:
    965             raise ValueError, "for the maxima algorithm the precision must be 53"
    966         return RR(maxima.eval("bessel_y(%s,%s)"%(float(nu),float(z))))
    967     elif algorithm == "pari":
    968         raise NotImplementedError, "The Y-Bessel function is only implemented for the maxima and scipy algorithms"
     800    # determine the numerical evaluation method
     801    if class_attrs['_system'] == 'mpmath':
     802        import mpmath
     803        mpmath_bessel_functions = { 'I': mpmath.besseli, 'J': mpmath.besselj,
     804                                    'K': mpmath.besselk, 'Y': mpmath.bessely }
     805        mpf = mpmath_bessel_functions[class_attrs['_type']]
     806        if class_attrs['_nargs'] == 1:
     807            def f(self, z, **kwds):
     808                return mpmath_utils.call(mpf, self._order, z, **kwds)
     809        else:
     810            def f(self, n, z, **kwds):
     811                return mpmath_utils.call(mpf, n, z, **kwds)
     812        class_attrs['_evalf_'] = f
     813    elif class_attrs['_system'] == 'maxima':
     814        raise NotImplementedError("maxima evaluation is not implemented yet")
     815    elif class_attrs['_system'] == 'scipy':
     816        raise NotImplementedError("scipy evaluation is not implemented yet")
    969817    else:
    970         raise ValueError, "unknown algorithm '%s'"%algorithm
    971    
    972 class Bessel():
    973     """
    974     A class implementing the I, J, K, and Y Bessel functions.
    975    
    976     EXAMPLES::
    977    
    978         sage: g = Bessel(2); g
    979         J_{2}
    980         sage: print g
    981         J-Bessel function of order 2
    982         sage: g.plot(0,10)
     818        raise NotImplementedError("unknown algorithm specified")
    983819
    984     ::
     820    # conversions are only defined when the index is
     821    # unspecified since the Bessel functions in Maxima take two parameters:
     822    # the index and the argument.
     823    if class_attrs['_nargs'] == 2:
     824        class_attrs['_conversions'] = _conversions_dict[class_attrs['_type']]
     825    else:
     826        class_attrs['_conversions'] = {}
    985827
    986         sage: Bessel(2, typ='I')(pi)
    987         2.61849485263445
    988         sage: Bessel(2, typ='J')(pi)
    989         0.485433932631509
    990         sage: Bessel(2, typ='K')(pi)
    991         0.0510986902537926
    992         sage: Bessel(2, typ='Y')(pi)
    993         -0.0999007139289404
    994     """
    995     def __init__(self, nu, typ = "J", algorithm = None, prec = 53):
     828    def __init__(self, attrs):
    996829        """
    997         Initializes new instance of the Bessel class.
     830        Initializes new instance of a Bessel class.
    998831
    999832        INPUT:
    1000833
    1001          - ``typ`` -- (default: J) the type of Bessel function: 'I', 'J', 'K'
    1002            or 'Y'.
    1003 
    1004          - ``algorithm`` -- (default: maxima for type Y, pari for other types)
    1005            algorithm to use to compute the Bessel function: 'pari', 'maxima' or
    1006            'scipy'.  Note that type K is not implemented in Maxima and type Y
    1007            is not implemented in PARI.
    1008 
    1009          - ``prec`` -- (default: 53) precision in bits of the Bessel function.
    1010            Only supported for the PARI algorithm.
     834            - ``attrs`` - a dictionary containing self attributes to set during
     835                          initialization. This should include at least the keys:
     836                          '_type', '_order', '_nargs', '_name', '_latex_name',
     837                          '_system', and '_nargs'.
    1011838
    1012839        EXAMPLES::
    1013840
    1014841            sage: g = Bessel(2); g
    1015             J_{2}
    1016             sage: Bessel(1,'I')
    1017             I_{1}
    1018             sage: Bessel(6, prec=120)(pi)
    1019             0.014545966982505560573660369604001804
    1020             sage: Bessel(6, algorithm="pari")(pi)
    1021             0.0145459669825056
     842            bessel_{2,J}
     843            sage: Bessel(0, 'Y')
     844            bessel_{0,Y}
     845        """
     846        for k in ['_type', '_order', '_nargs', '_name', '_latex_name', '_system',
     847                  '_conversions']:
     848            setattr(self, k, attrs[k])
    1022849
    1023         For the Bessel J-function, Maxima returns a symbolic result.  For
    1024         types I and Y, we always get a numeric result::
     850        BuiltinFunction.__init__(self, self._name, nargs=self._nargs,
     851                                 latex_name=self._latex_name,
     852                                 conversions=self._conversions)
    1025853
    1026             sage: b = Bessel(6, algorithm="maxima")(pi); b
    1027             bessel_j(6, pi)
    1028             sage: b.n(53)
    1029             0.0145459669825056
    1030             sage: Bessel(6, typ='I', algorithm="maxima")(pi)
    1031             0.0294619840059568
    1032             sage: Bessel(6, typ='Y', algorithm="maxima")(pi)
    1033             -4.33932818939038
     854    def _eval_(self, *args):
     855        """
     856        EXAMPLES::
    1034857
    1035         SciPy usually gives less precise results::
     858            sage: Bessel(0, 'J')(1)
     859            bessel_{0,J}(1)
     860            sage: Bessel(0, 'J')(1.0)
     861            0.765197686557966
     862            sage: Bessel(1, 'Y')(x)
     863            bessel_{1,Y}(x)
     864        """
     865        if (len(args) == 1 and
     866           (not isinstance(args[0], Expression) and is_inexact(args[0]))):
     867                return self._evalf_(args[0], parent=parent(args[0]))
     868        elif (len(args) == 2 and
     869             (not isinstance(args[0], Expression) and
     870             not isinstance(args[1], Expression) and
     871             (is_inexact(args[0]) or is_inexact(args[1])))):
     872                coercion_model = sage.structure.element.get_coercion_model()
     873                n, z = coercion_model.canonical_coercion(*args)
     874                return self._evalf_(n, z, parent=parent(n))
     875        return None # leaves the expression unevaluated
    1036876
    1037             sage: Bessel(6, algorithm="scipy")(pi)
    1038             0.0145459669825000...
     877    def _derivative_(self, z, diff_param=None):
     878        r"""
     879        The derivative of any of the Bessel functions of order `n` and argument `z`
     880        is:
    1039881
    1040         TESTS::
     882        .. math::
    1041883
    1042             sage: Bessel(1,'Z')
    1043             Traceback (most recent call last):
    1044             ...
    1045             ValueError: typ must be one of I, J, K, Y
     884            \left\( \frac{1}{z} \frac{d}{dz} \right\)^k \left\( z^n J_n(z) \right\) = z^{n-k} J_{n-k}(z)
     885
     886        Specializing to `k = 1` we get:
     887
     888        ..math::
     889
     890            \frac{d}{dz} (z^n J_n(z)) = z^n J_{n-1}(z)
     891
     892        It follows that,
     893
     894        ..math::
     895
     896            \frac{d}{dz} J_n(z) = \frac{1}{z^n} \left\(z^n J_{n-1}(z) - n z^{n-1} J_n(z) \right\)
     897
     898        The same formula holds for 'J' replaced by 'Y', 'K', or 'I'.
     899
     900        EXAMPLES::
     901
     902            sage: f(x) = Bessel(10)(x)
     903            sage: derivative(f, x)
     904            x |--> -10*bessel_{10,J}(x)/x + bessel_{9,J}(x)
     905            sage: f(x) = Bessel(10)(Bessel(-2)(x))
     906            sage: derivative(f, x)
     907            x |--> -(10*bessel_{10,J}(bessel_{-2,J}(x))/bessel_{-2,J}(x) - bessel_{9,J}(bessel_{-2,J}(x)))*(2*bessel_{-2,J}(x)/x + bessel_{-3,J}(x))
    1046908        """
    1047         if not (typ in ['I', 'J', 'K', 'Y']):
    1048             raise ValueError, "typ must be one of I, J, K, Y"
     909        n = self._order
     910        if n is None:
     911            raise NotImplementedError("derivative of two parameter Bessel function is not implemented")
     912        return Bessel(n-1, self._type)(z) - n/z*Bessel(n, self._type)(z)
    1049913
    1050         # Did the user ask for the default algorithm?
    1051         if algorithm is None:
    1052             if typ == 'Y':
    1053                 algorithm = 'maxima'
    1054             else:
    1055                 algorithm = 'pari'
    1056 
    1057         self._system = algorithm
    1058         self._order = nu
    1059         self._type = typ
    1060         prec = int(prec)
    1061         if prec < 0:
    1062             raise ValueError, "prec must be a positive integer"
    1063         self._prec = int(prec)
    1064 
    1065     def __str__(self):
    1066         """
    1067         Returns a string representation of this Bessel object.
    1068 
    1069         TEST::
    1070 
    1071             sage: a = Bessel(1,'I')
    1072             sage: str(a)
    1073             'I-Bessel function of order 1'
    1074         """
    1075         return self.type()+"-Bessel function of order "+str(self.order())
    1076    
    1077     def __repr__(self):
    1078         """
    1079         Returns a string representation of this Bessel object.
    1080 
    1081         TESTS::
    1082 
    1083             sage: Bessel(1,'I')
    1084             I_{1}
    1085         """
    1086         return self.type()+"_{"+str(self.order())+"}"
    1087    
    1088     def type(self):
     914    def get_type(self):
    1089915        """
    1090916        Returns the type of this Bessel object.
    1091917
    1092         TEST::
     918        EXAMPLES::
    1093919
    1094920            sage: a = Bessel(3,'K')
    1095             sage: a.type()
     921            sage: a.get_type()
    1096922            'K'
    1097923        """
    1098924        return self._type
    1099    
    1100     def prec(self):
    1101         """
    1102         Returns the precision (in number of bits) used to represent this
    1103         Bessel function.
    1104925
    1105         TESTS::
    1106 
    1107             sage: a = Bessel(3,'K')
    1108             sage: a.prec()
    1109             53
    1110             sage: B = Bessel(20,prec=100); B
    1111             J_{20}
    1112             sage: B.prec()
    1113             100
    1114         """
    1115         return self._prec
    1116 
    1117     def order(self):
     926    def get_order(self):
    1118927        """
    1119928        Returns the order of this Bessel function.
    1120929
    1121         TEST::
     930        EXAMPLES::
    1122931
    1123932            sage: a = Bessel(3,'K')
    1124             sage: a.order()
     933            sage: a.get_order()
    1125934            3
    1126935        """
    1127936        return self._order
    1128937
    1129     def system(self):
     938    def get_system(self):
    1130939        """
    1131         Returns the package used, e.g. Maxima, PARI, or SciPy, to compute with
     940        Returns the package used, e.g. 'mpmath', to compute with
    1132941        this Bessel function.
    1133942
    1134         TESTS::
     943        EXAMPLES::
    1135944
    1136             sage: Bessel(20,algorithm='maxima').system()
    1137             'maxima'
    1138             sage: Bessel(20,prec=100).system()
    1139             'pari'
     945            sage: Bessel(20).get_system()
     946            'mpmath'
    1140947        """
    1141948        return self._system
    1142949
    1143     def __call__(self,z):
    1144         """
    1145         Implements evaluation of all the Bessel functions directly
    1146         from the Bessel class. This essentially allows a Bessel object to
    1147         behave like a function that can be invoked.
     950    # populate the new classes attributes
     951    class_attrs.update({'__init__': __init__,
     952                        '_eval_': _eval_,
     953                        '_derivative_': _derivative_,
     954                        'get_type': get_type,
     955                        'get_system': get_system,
     956                        'get_order': get_order})
    1148957
    1149         TESTS::
     958    # last step: create the new subclass of BuiltinFunction and return its
     959    # instantiation.
     960    if class_attrs['_order'] is None:
     961        class_name = "Bessel_Function_%s" % (class_attrs['_type'],)
     962    else:
     963        class_name = "Bessel_Function_%d_%s" % (class_attrs['_order'], class_attrs['_type'])
     964    new_class = type(class_name, (BuiltinFunction,), class_attrs)
     965    return new_class(class_attrs)
    1150966
    1151             sage: Bessel(3,'K')(5.0)
    1152             0.00829176841523093
    1153             sage: Bessel(20,algorithm='maxima')(5.0)
    1154             2.77033005213e-11
    1155             sage: Bessel(20,prec=100)(5.0101010101010101)
    1156             2.8809188227195382093062257967e-11
    1157             sage: B = Bessel(2,'Y',algorithm='scipy',prec=50)
    1158             sage: B(2.0)
    1159             Traceback (most recent call last):
    1160             ...
    1161             ValueError: for the scipy algorithm the precision must be 53
    1162         """
    1163         nu = self.order()
    1164         t = self.type()
    1165         s = self.system()
    1166         p = self.prec()
    1167         if t == "I":
    1168             return bessel_I(nu,z,algorithm=s,prec=p)
    1169         if t == "J":
    1170             return bessel_J(nu,z,algorithm=s,prec=p)
    1171         if t == "K":
    1172             return bessel_K(nu,z,algorithm=s,prec=p)
    1173         if t == "Y":
    1174             return bessel_Y(nu,z,algorithm=s,prec=p)
    1175        
    1176     def plot(self,a,b):
    1177         """
    1178         Enables easy plotting of all the Bessel functions directly
    1179         from the Bessel class.
     967# Construct top level Bessel functions of the 4 types
     968bessel_I = Bessel(typ='I')
     969bessel_J = Bessel(typ='J')
     970bessel_K = Bessel(typ='K')
     971bessel_Y = Bessel(typ='Y')
    1180972
    1181         TESTS::
    1182 
    1183             sage: plot(Bessel(2),3,4)
    1184             sage: Bessel(2).plot(3,4)
    1185             sage: P = Bessel(2,'I').plot(1,5)
    1186             sage: P += Bessel(2,'J').plot(1,5)
    1187             sage: P += Bessel(2,'K').plot(1,5)
    1188             sage: P += Bessel(2,'Y').plot(1,5)
    1189             sage: show(P)
    1190         """
    1191         nu = self.order()
    1192         s = self.system()
    1193         t = self.type()
    1194         if t == "I":
    1195             f = lambda z: bessel_I(nu,z,s) 
    1196             P = plot(f,a,b)
    1197         if t == "J":
    1198             f = lambda z: bessel_J(nu,z,s)
    1199             P = plot(f,a,b)
    1200         if t == "K":
    1201             f = lambda z: bessel_K(nu,z,s)
    1202             P = plot(f,a,b)
    1203         if t == "Y":
    1204             f = lambda z: bessel_Y(nu,z,s)
    1205             P = plot(f,a,b)
    1206         return P
    1207    
    1208973def hypergeometric_U(alpha,beta,x,algorithm="pari",prec=53):
    1209974    r"""
    1210975    Default is a wrap of PARI's hyperu(alpha,beta,x) function.