Ticket #2516: hyper.patch

File hyper.patch, 40.1 KB (added by fredrik.johansson, 11 years ago)

tentative implementation of generalized hypergeometric functions

  • sage/functions/all.py

    # HG changeset patch
    # User fredrik@scv
    # Date 1280114065 -7200
    # Node ID 7d20ef3a988c8054f41b4179735eb765d9a4980c
    # Parent  2ecc6b33a4469a76abc0133bc48aeea31c87e70b
    tentative implementation of hypergeometric functions
    
    diff -r 2ecc6b33a446 -r 7d20ef3a988c sage/functions/all.py
    a b  
    11from piecewise import piecewise, Piecewise
    22
    3 from trig import ( sin, cos, sec, csc, cot, tan, 
     3from trig import ( sin, cos, sec, csc, cot, tan,
    44                   asin, acos, atan,
    55                   acot, acsc, asec,
    6                    arcsin, arccos, arctan, 
     6                   arcsin, arccos, arctan,
    77                   arccot, arccsc, arcsec,
    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 (exponential_integral_1,
    2828                            zeta, zetaderiv, zeta_symmetric,
    29                             Li, Ei, 
     29                            Li, Ei,
    3030                            dickman_rho)
    3131
    3232from special import (bessel_I, bessel_J, bessel_K, bessel_Y,
     
    3939                     elliptic_f, elliptic_ec, elliptic_eu,
    4040                     elliptic_kc, elliptic_pi, elliptic_j,
    4141                     airy_ai, airy_bi)
    42                        
     42
    4343from orthogonal_polys import (chebyshev_T,
    4444                              chebyshev_U,
    4545                              gen_laguerre,
     
    6060from wigner import (wigner_3j, clebsch_gordan, racah, wigner_6j,
    6161                    wigner_9j, gaunt)
    6262
    63 from generalized import (dirac_delta, heaviside, unit_step, sgn, sign, 
     63from generalized import (dirac_delta, heaviside, unit_step, sgn, sign,
    6464                         kronecker_delta)
    6565
    6666from min_max import max_symbolic, min_symbolic
     67
     68from hypergeometric import (FormalPFQ, PFQ, hyp,
     69     bessel_i, bessel_j, bessel_k, bessel_y)
  • new file sage/functions/hypergeometric.py

    diff -r 2ecc6b33a446 -r 7d20ef3a988c sage/functions/hypergeometric.py
    - +  
     1"""
     2This module implements manipulation of infinite hypergeometric series
     3represented in standard parametric form (as $\,_pFq$ functions).
     4
     5Also included are some hypergeometric-type special functions
     6(currently, only Bessel functions).
     7
     8"""
     9
     10from sage.rings.all import Integer, ZZ, QQ, Infinity, binomial, rising_factorial
     11from sage.rings.all import factorial as fac
     12from sage.functions.other import sqrt, gamma, real_part
     13from sage.functions.log import exp, log
     14from sage.functions.trig import cos, sin
     15from sage.functions.hyperbolic import cosh, sinh
     16from sage.functions.other import erf
     17from sage.symbolic.constants import pi
     18from sage.symbolic.all import I
     19from sage.symbolic.function import BuiltinFunction
     20from sage.symbolic.ring import SR
     21
     22from sage.misc.latex import latex
     23
     24
     25def rational_param_as_tuple(x):
     26    """
     27    Utility function for converting rational PFQ parameters to
     28    tuples (which mpmath handles more efficiently).
     29
     30        sage: from sage.functions.hypergeometric import rational_param_as_tuple as r
     31        sage: r(1/2)
     32        (1, 2)
     33        sage: r(3)
     34        3
     35        sage: r(pi)
     36        pi
     37
     38    """
     39    try:
     40        x = x.pyobject()
     41    except (TypeError, AttributeError):
     42        pass
     43    try:
     44        if x.parent() is QQ:
     45            p = int(x.numer())
     46            q = int(x.denom())
     47            return p, q
     48    except (TypeError, AttributeError):
     49        pass
     50    return x
     51
     52def expand_hyper(x):
     53    """
     54    Expand compatible functions as hypergeometric series.
     55
     56    EXAMPLES ::
     57
     58        sage: from sage.functions.hypergeometric import expand_hyper
     59        sage: expand_hyper(bessel_j(2,x))
     60        1/8*x^2*PFQ(0, 1, 3, -1/4*x^2)
     61
     62    """
     63    try:
     64        op = x.operator()
     65    except AttributeError:
     66        f = SR(x)
     67        op = x.operator()
     68    if hasattr(op, "_expand_hyper"):
     69        return op._expand_hyper(*x.operands())
     70    if op is None:
     71        return x
     72    return reduce(op, [expand_hyper(t) for t in x.operands()])
     73
     74class ForceSR(BuiltinFunction):
     75    """
     76    Force inputs to SR
     77
     78    XXX: this does not work as intended because of the way
     79    BuiltinFunction.__call__ works.
     80    """
     81    def _eval_(self, *args):
     82        def __force_SR(x):
     83            # use py_scalar_parent?
     84            if isinstance(x, (int, long)):
     85                return SR(Integer(x))
     86            else:
     87                return SR(x)
     88        args = [__force_SR(a) for a in args]
     89        value = self._eval_symbolic_(*args)
     90        if value is None:
     91            return None
     92        g = __force_SR(value)
     93        return g
     94
     95
     96class FormalPFQ:
     97    r"""
     98    Represents a (formal) generalized infinite hypergeometric series.
     99    """
     100
     101    def __init__(self, a, b, z):
     102        """
     103        EXAMPLES ::
     104
     105            sage: FormalPFQ([],[],1)
     106            0F0([]; []; 1)
     107            sage: FormalPFQ([],[1],1)
     108            0F1([]; [1]; 1)
     109            sage: FormalPFQ([2,3],[1],1)
     110            2F1([2, 3]; [1]; 1)
     111
     112        """
     113        self.a = list(a)
     114        self.b = list(b)
     115        self.z = z
     116
     117    def __eq__(self, other):
     118        """
     119        TESTS ::
     120
     121            sage: FormalPFQ([],[],1) == FormalPFQ([],[],1)
     122            True
     123            sage: FormalPFQ([2],[],1) == FormalPFQ([],[],1)
     124            False
     125
     126        """
     127        if not isinstance(other, FormalPFQ):
     128            return False
     129        return self.a == other.a and self.b == other.b and self.z == other.z
     130
     131    def __ne__(self, other):
     132        """
     133        TESTS ::
     134
     135            sage: FormalPFQ([],[],1) != FormalPFQ([],[],1)
     136            False
     137        """
     138        return not self.__eq__(other)
     139
     140    def __hash__(self):
     141        """
     142        TESTS ::
     143
     144            sage: a = FormalPFQ([1],[2],3)
     145            sage: b = FormalPFQ([1],[2],3)
     146            sage: hash(a) == hash(b)
     147            True
     148
     149        """
     150        return hash((tuple(self.a), tuple(self.b), self.z))
     151
     152    @property
     153    def p(self):
     154        """
     155        Return the numerator degree of self.
     156
     157        EXAMPLES ::
     158
     159            sage: FormalPFQ([],[],1).p
     160            0
     161            sage: FormalPFQ([1],[],1).p
     162            1
     163            sage: FormalPFQ([1,1],[],1).p
     164            2
     165
     166        """
     167        return len(self.a)
     168
     169    @property
     170    def q(self):
     171        """
     172        Return the denominator degree of self.
     173
     174        EXAMPLES ::
     175
     176            sage: FormalPFQ([],[],1).q
     177            0
     178            sage: FormalPFQ([],[1],1).q
     179            1
     180            sage: FormalPFQ([],[1,1],1).q
     181            2
     182
     183        """
     184        return len(self.b)
     185
     186    def __repr__(self):
     187        """
     188        TESTS ::
     189
     190            sage: FormalPFQ([1],[2],3)
     191            1F1([1]; [2]; 3)
     192
     193        """
     194        return "%iF%i(%s; %s; %s)" % (self.p, self.q, self.a, self.b, self.z)
     195
     196    def to_PFQ(self):
     197        """
     198        Convert to a PFQ (symbolic function) instance.
     199
     200        EXAMPLES ::
     201
     202            sage: FormalPFQ([1,2],[3,4],5).to_PFQ()
     203            PFQ(2, 2, 1, 2, 3, 4, 5)
     204
     205        """
     206        return PFQ(*([self.p, self.q] + self.a + self.b + [self.z]))
     207
     208    def _print_latex_(self):
     209        r"""
     210        TESTS ::
     211
     212            sage: latex(FormalPFQ([1,1],[2],-1))
     213            \texttt{2F1([1, 1]; [2]; -1)}
     214
     215        """
     216        p = self.p
     217        q = self.q
     218        a = ",".join(latex(c) for c in self.a)
     219        b = ",".join(latex(c) for c in self.b)
     220        z = latex(self.z)
     221        return r"\,_%iF_%i\left(\begin{matrix} %s \\ %s \end{matrix} ; %s \right)" % (p, q, a, b, z)
     222
     223    def n(self, **kwargs):
     224        """
     225        TESTS ::
     226
     227            sage: FormalPFQ([1,1],[2],-1).n()
     228            0.693147180559945
     229
     230        """
     231        from sage.libs.mpmath.all import hyper, call
     232        a = [rational_param_as_tuple(c) for c in self.a]
     233        b = [rational_param_as_tuple(c) for c in self.b]
     234        return call(hyper, a, b, self.z, **kwargs)
     235
     236    def sort_parameters(self):
     237        """
     238        Sort parameters in a canonical order.
     239
     240        EXAMPLES::
     241
     242            sage: FormalPFQ([2,1,3],[5,4],1/2).sort_parameters()
     243            3F2([1, 2, 3]; [4, 5]; 1/2)
     244
     245        """
     246        return FormalPFQ(sorted(self.a), sorted(self.b), self.z)
     247
     248    def eliminate_parameters(self):
     249        """
     250        Eliminate repeated parameters.
     251
     252        EXAMPLES::
     253
     254            sage: FormalPFQ([1,1,2,5],[5,1,4],1/2).eliminate_parameters()
     255            2F1([1, 2]; [4]; 1/2)
     256
     257        """
     258        aa = self.a[:]
     259        bb = self.b[:]
     260        p = self.p
     261        q = self.q
     262        i = 0
     263        while i < q and aa:
     264            b = bb[i]
     265            if b in aa:
     266                aa.remove(b)
     267                bb.remove(b)
     268                p -= 1
     269                q -= 1
     270            else:
     271                i += 1
     272        if (p,q) != (self.p,self.q):
     273            return FormalPFQ(aa, bb, self.z)
     274        return self
     275
     276    def is_termwise_finite(self):
     277        """
     278        Determine whether all terms of self are finite. Any infinite
     279        terms or ambiguous terms beyond the first zero, if one exists,
     280        are ignored.
     281
     282        Ambiguous cases (where a term is the product of both zero
     283        and an infinity) are not considered finite.
     284
     285            sage: FormalPFQ([2],[3,4],5).is_termwise_finite()
     286            True
     287            sage: FormalPFQ([2],[-3,4],5).is_termwise_finite()
     288            False
     289            sage: FormalPFQ([-2],[-3,4],5).is_termwise_finite()
     290            True
     291            sage: FormalPFQ([-3],[-3,4],5).is_termwise_finite()  # ambiguous
     292            False
     293
     294            sage: FormalPFQ([0],[-1],5).is_termwise_finite()
     295            True
     296            sage: FormalPFQ([0],[0],5).is_termwise_finite()   # ambiguous
     297            False
     298            sage: FormalPFQ([1],[2],Infinity).is_termwise_finite()
     299            False
     300            sage: FormalPFQ([0],[0],Infinity).is_termwise_finite()  # ambiguous
     301            False
     302            sage: FormalPFQ([0],[],Infinity).is_termwise_finite()   # ambiguous
     303            False
     304
     305        """
     306        if self.z == 0:
     307            return 0 not in self.b
     308        if abs(self.z) == Infinity:
     309            return False
     310        if abs(self.z) == Infinity:
     311            return False
     312        for b in self.b:
     313            if b in ZZ and b <= 0:
     314                if any((a in ZZ) and (b < a <= 0) for a in self.a):
     315                    continue
     316                return False
     317        return True
     318
     319    def is_terminating(self):
     320        """
     321        Determine whether the series represented by self terminates
     322        after a finite number of terms, i.e. whether any of the
     323        numerator parameters are nonnegative integers (with no
     324        preceding nonnegative denominator parameters), or z = 0.
     325
     326        If terminating, the series represents a polynomial of self.z
     327
     328        EXAMPLES ::
     329
     330            sage: z = var('z')
     331            sage: FormalPFQ([1,2],[3,4],z).is_terminating()
     332            False
     333            sage: FormalPFQ([1,-2],[3,4],z).is_terminating()
     334            True
     335            sage: FormalPFQ([1,-2],[],z).is_terminating()
     336            True
     337            sage: FormalPFQ([2,3],[],0).is_terminating()
     338            True
     339
     340        """
     341        if self.z == 0:
     342            return True
     343        for a in self.a:
     344            if (a in ZZ) and (a <= 0):
     345                return self.is_termwise_finite()
     346        return False
     347
     348    def is_absolutely_convergent(self):
     349        """
     350        Determine whether self converges absolutely as an infinite series.
     351        False is returned if not all terms are finite.
     352
     353        EXAMPLES ::
     354
     355        Degree giving infinite radius of convergence::
     356
     357            sage: FormalPFQ([2,3],[4,5],6).is_absolutely_convergent()
     358            True
     359            sage: FormalPFQ([2,3],[-4,5],6).is_absolutely_convergent()  # undefined
     360            False
     361            sage: FormalPFQ([2,3],[-4,5],Infinity).is_absolutely_convergent()  # undefined
     362            False
     363
     364        Ordinary geometric series (unit radius of convergence)::
     365
     366            sage: FormalPFQ([1],[],1/2).is_absolutely_convergent()
     367            True
     368            sage: FormalPFQ([1],[],2).is_absolutely_convergent()
     369            False
     370            sage: FormalPFQ([1],[],1).is_absolutely_convergent()
     371            False
     372            sage: FormalPFQ([1],[],-1).is_absolutely_convergent()
     373            False
     374            sage: FormalPFQ([1],[],-1).n()    # Sum still exists
     375            0.500000000000000
     376
     377        Degree $p = q+1$ (unit radius of convergence)::
     378
     379            sage: FormalPFQ([2,3],[4],6).is_absolutely_convergent()
     380            False
     381            sage: FormalPFQ([2,3],[4],1).is_absolutely_convergent()
     382            False
     383            sage: FormalPFQ([2,3],[5],1).is_absolutely_convergent()
     384            False
     385            sage: FormalPFQ([2,3],[6],1).is_absolutely_convergent()
     386            True
     387            sage: FormalPFQ([-2,3],[4],5).is_absolutely_convergent()
     388            True
     389            sage: FormalPFQ([2,-3],[4],5).is_absolutely_convergent()
     390            True
     391            sage: FormalPFQ([2,-3],[-4],5).is_absolutely_convergent()
     392            True
     393            sage: FormalPFQ([2,-3],[-1],5).is_absolutely_convergent()
     394            False
     395
     396        Degree giving zero radius of convergence::
     397
     398            sage: FormalPFQ([1,2,3],[4],2).is_absolutely_convergent()
     399            False
     400            sage: FormalPFQ([1,2,3],[4],1/2).is_absolutely_convergent()
     401            False
     402            sage: FormalPFQ([1,2,3],[4],0).is_absolutely_convergent()
     403            True
     404            sage: FormalPFQ([1,2,-3],[4],1/2).is_absolutely_convergent()  # polynomial
     405            True
     406
     407        """
     408        if not self.is_termwise_finite():
     409            return False
     410        if self.p <= self.q: return True
     411        if self.is_terminating(): return True
     412        if self.p == self.q+1:
     413            if abs(self.z) < 1:
     414                return True
     415            if abs(self.z) == 1:
     416                if real_part(sum(self.b) - sum(self.a)) > 0:
     417                    return True
     418        return False
     419
     420    def terms(self, n=None):
     421        """
     422        Generate the terms of self (optionally only n terms).
     423
     424        EXAMPLES ::
     425
     426            sage: list(FormalPFQ([-2,1],[3,4],x).terms())
     427            [1, -1/6*x, 1/120*x^2]
     428            sage: list(FormalPFQ([-2,1],[3,4],x).terms(2))
     429            [1, -1/6*x]
     430            sage: list(FormalPFQ([-2,1],[3,4],x).terms(0))
     431            []
     432
     433        """
     434        if n is None:
     435            n = Infinity
     436        t = Integer(1)
     437        k = 1
     438        while k <= n:
     439            yield t
     440            for a in self.a: t *= (a+k-1)
     441            for b in self.b: t /= (b+k-1)
     442            t *= self.z
     443            if t == 0:
     444                break
     445            t /= k
     446            k += 1
     447
     448    def add_terms(self, n=None):
     449        """
     450        Add the first n terms of self. If n=None, self must be
     451        terminating and all terms will be expanded.
     452
     453        EXAMPLES::
     454
     455            sage: z = var('z')
     456            sage: FormalPFQ([],[],z).add_terms(0)
     457            0
     458            sage: FormalPFQ([],[],z).add_terms(1)
     459            1
     460            sage: FormalPFQ([],[],z).add_terms(2)
     461            z + 1
     462            sage: FormalPFQ([],[],z).add_terms(3)
     463            1/2*z^2 + z + 1
     464
     465            sage: FormalPFQ([],[],z).add_terms()
     466            Traceback (most recent call last):
     467                ...
     468            ValueError: self must be terminating
     469
     470            sage: FormalPFQ([-2],[],z).add_terms()
     471            z^2 - 2*z + 1
     472            sage: FormalPFQ([-2],[],z).add_terms(3)
     473            z^2 - 2*z + 1
     474            sage: FormalPFQ([-2],[],z).add_terms(6)
     475            z^2 - 2*z + 1
     476            sage: FormalPFQ([-2],[],z).add_terms(2)
     477            -2*z + 1
     478
     479            sage: FormalPFQ([],[],z).add_terms(6)
     480            1/120*z^5 + 1/24*z^4 + 1/6*z^3 + 1/2*z^2 + z + 1
     481            sage: FormalPFQ([1],[],z).add_terms(6)
     482            z^5 + z^4 + z^3 + z^2 + z + 1
     483            sage: FormalPFQ([],[1/2],-z^2/4).add_terms(6)
     484            -1/3628800*z^10 + 1/40320*z^8 - 1/720*z^6 + 1/24*z^4 - 1/2*z^2 + 1
     485
     486            sage: FormalPFQ([1,2],[3],1/3).add_terms(6).n()
     487            1.29788359788360
     488            sage: FormalPFQ([1,2],[3],1/3).n()
     489            1.29837194594696
     490
     491        """
     492        if n is None:
     493            if not self.is_terminating():
     494                raise ValueError("self must be terminating")
     495        return sum(self.terms(n))
     496
     497    def deflate(self):
     498        """
     499        Rewrite as a linear combination of functions of strictly lower degree
     500        by eliminating all parameters a[i] and b[j] such that a[i] = b[i] + m
     501        for nonnegative integer m.
     502
     503        EXAMPLES ::
     504
     505            sage: x = FormalPFQ([5],[4],3)
     506            sage: y = x.deflate()
     507            sage: y
     508            (7/4)*0F0([]; []; 3)
     509            sage: x.n(); y.n()
     510            35.1496896155784
     511            35.1496896155784
     512
     513            sage: x = FormalPFQ([6,1],[3,4,5],10)
     514            sage: y = x.deflate()
     515            sage: y
     516            (1)*1F2([1]; [4, 5]; 10) + (1/2)*1F2([2]; [5, 6]; 10) + (1/12)*1F2([3]; [6, 7]; 10) + (1/252)*1F2([4]; [7, 8]; 10)
     517            sage: x.n(); y.n()
     518            2.87893612686782
     519            2.87893612686782
     520
     521            sage: x = FormalPFQ([6,7],[3,4,5],10)
     522            sage: y = x.deflate()
     523            sage: y
     524            (1)*0F1([]; [5]; 10) + (5)*0F1([]; [6]; 10) + (19/3)*0F1([]; [7]; 10) + (181/63)*0F1([]; [8]; 10) + (265/504)*0F1([]; [9]; 10) + (25/648)*0F1([]; [10]; 10) + (25/27216)*0F1([]; [11]; 10)
     525            sage: x.n(); y.n()
     526            63.0734110716969
     527            63.0734110716969
     528
     529        """
     530        self = self.eliminate_parameters()
     531        for i,a in enumerate(self.a):
     532            for j,b in enumerate(self.b):
     533                m = a-b
     534                if m in ZZ and m > 0:
     535                    aa = self.a[:i] + self.a[i+1:]
     536                    bb = self.b[:j] + self.b[j+1:]
     537                    terms = []
     538                    for k in xrange(m+1):
     539                        # TODO: could rewrite prefactors as recurrence
     540                        term = binomial(m,k)
     541                        for c in aa: term *= rising_factorial(c,k)
     542                        for c in bb: term /= rising_factorial(c,k)
     543                        term *= self.z**k
     544                        term /= rising_factorial(a-m, k)
     545                        F = FormalPFQ([c+k for c in aa], [c+k for c in bb], self.z)
     546                        Fterms = F.deflate().terms
     547                        terms += [(term*termG, G) for (termG, G) in Fterms]
     548                    return FormalPFQLinearCombination(terms)
     549        return FormalPFQLinearCombination([(1,self)])
     550
     551    def closed_form(self):
     552        """
     553        Try to evaluate self in closed form using elementary
     554        (and other simple) functions.
     555
     556        It may be necessary to call .deflate() first to
     557        find some closed forms.
     558
     559        EXAMPLES ::
     560
     561            sage: z = var('z')
     562            sage: FormalPFQ([1],[],1+z).closed_form()
     563            -1/z
     564            sage: FormalPFQ([],[],1+z).closed_form()
     565            e^(z + 1)
     566            sage: FormalPFQ([],[1/2],4).closed_form()
     567            cosh(4)
     568            sage: FormalPFQ([],[3/2],4).closed_form()
     569            1/4*sinh(4)
     570            sage: FormalPFQ([],[5/2],4).closed_form()
     571            -3/64*sinh(4) + 3/16*cosh(4)
     572            sage: FormalPFQ([],[-3/2],4).closed_form()
     573            -4*sinh(4) + 19/3*cosh(4)
     574            sage: FormalPFQ([-3,1],[var('a')],z).closed_form()
     575            -6*z^3/((a + 1)*(a + 2)*a) + 6*z^2/((a + 1)*a) - 3*z/a + 1
     576
     577        """
     578        if self.is_terminating():
     579            return self.add_terms()
     580
     581        # necessary?
     582        new = self.sort_parameters()
     583        new = self.eliminate_parameters()
     584        a, b, z, p, q = new.a, new.b, new.z, new.p, new.q
     585
     586        if z == 0:
     587            return Integer(1)
     588        if p == q == 0:
     589            return exp(z)
     590        if p == 1 and q == 0:
     591            return (1-z)**(-a[0])
     592
     593        if p == 0 and q == 1:
     594            # TODO: make this require only linear time
     595            def _0f1(b,z):
     596                F12 = cosh(2*sqrt(z))
     597                F32 = sinh(2*sqrt(z))/(2*sqrt(z))
     598                if 2*b == 1: return F12
     599                if 2*b == 3: return F32
     600                if 2*b > 3:
     601                    return (b-2)*(b-1)/z * (_0f1(b-2,z) - _0f1(b-1,z))
     602                if 2*b < 1:
     603                    return _0f1(b+1,z) + z/(b*(b+1)) * _0f1(b+2,z)
     604                raise ValueError
     605            # Can evaluate 0F1 in terms of elementary functions when
     606            # the parameter is a half-integer
     607            if 2*b[0] in ZZ and b[0] not in ZZ:
     608                return _0f1(b[0], z)
     609
     610        # Confluent hypergeometric function
     611        if p == 1 and q == 1:
     612            aa, bb = a[0], b[0]
     613            if aa*2 == 1 and bb*2 == 3:
     614                t = sqrt(-z)
     615                return sqrt(pi)/2*erf(t)/t
     616            if a == 1 and b == 2:
     617                return (exp(z)-1)/z
     618            n, m = aa, bb
     619            if n in ZZ and m in ZZ and m > 0 and n > 0:
     620                rf = rising_factorial
     621                if m <= n:
     622                    return exp(z)*sum(rf(m-n,k)*(-z)**k/fac(k)/rf(m,k) for k in xrange(n-m+1))
     623                else:
     624                    T = sum(rf(n-m+1,k)*z**k/(fac(k)*rf(2-m,k)) for k in xrange(m-n))
     625                    U = sum(rf(1-n,k)*(-z)**k/(fac(k)*rf(2-m,k)) for k in xrange(n))
     626                    return fac(m-2)*rf(1-m,n)*z**(1-m)/fac(n-1)*(T - exp(z)*U)
     627
     628        if p == 2 and q == 1:
     629            R12 = QQ('1/2')
     630            R32 = QQ('3/2')
     631            def _2f1(a,b,c,z):
     632                """
     633                Evaluation of 2F1(a,b,c,z), assuming a,b,c positive
     634                integers or half-integers
     635                """
     636                if b == c:
     637                    return (1-z)**(-a)
     638                if a == c:
     639                    return (1-z)**(-b)
     640                if a == 0 or b == 0:
     641                    return Integer(1)
     642                if a > b:
     643                    a, b = b, a
     644                if b >= 2:
     645                    F1 = _2f1(a,b-1,c,z)
     646                    F2 = _2f1(a,b-2,c,z)
     647                    q = (b-1)*(z-1)
     648                    return ((c-2*b+2+(b-a-1)*z)*F1 + (b-c-1)*F2) / q
     649                if c > 2:
     650                    # how to handle this case?
     651                    if a-c+1 == 0 or b-c+1 == 0:
     652                        raise NotImplementedError
     653                    F1 = _2f1(a,b,c-1,z)
     654                    F2 = _2f1(a,b,c-2,z)
     655                    r1 = (c-1)*(2-c-(a+b-2*c+3)*z)
     656                    r2 = (c-1)*(c-2)*(1-z)
     657                    q = (a-c+1)*(b-c+1)*z
     658                    return (r1*F1 + r2*F2) / q
     659                if (a,b,c) == (R12,1,2): return (2-2*sqrt(1-z))/z
     660                if (a,b,c) == (1,1,2): return -log(1-z)/z
     661                if (a,b,c) == (1,R32,R12): return (1+z)/(1-z)**2
     662                if (a,b,c) == (1,R32,2): return 2*(1/sqrt(1-z)-1)/z
     663                if (a,b,c) == (R32,2,R12): return (1+3*z)/(1-z)**3
     664                if (a,b,c) == (R32,2,1): return (2+z)/(2*(sqrt(1-z)*(1-z)**2))
     665                if (a,b,c) == (2,2,1): return (1+z)/(1-z)**3
     666                raise NotImplementedError
     667            aa, bb = a; cc, = b
     668            if z == 1:
     669                return gamma(cc)*gamma(cc-aa-bb)/gamma(cc-aa)/gamma(cc-bb)
     670            if (aa*2) in ZZ and (bb*2) in ZZ and (cc*2) in ZZ and \
     671                aa > 0 and bb > 0 and cc > 0:
     672                try:
     673                    return _2f1(aa,bb,cc,z)
     674                except NotImplementedError:
     675                    pass
     676
     677        return self
     678
     679
     680class FormalPFQLinearCombination:
     681    """
     682    Holds a weighted sum of FormalPFQs.
     683    """
     684
     685    def __init__(self, terms=[], collect=True):
     686        """
     687            sage: from sage.functions.hypergeometric import FormalPFQLinearCombination
     688            sage: FormalPFQLinearCombination([(3,FormalPFQ([],[1],x))])
     689            (3)*0F1([]; [1]; x)
     690
     691        """
     692        if collect:
     693            unique = []
     694            counts = []
     695            for c,f in terms:
     696                if f in unique:
     697                    counts[unique.index(f)] += c
     698                else:
     699                    unique.append(f)
     700                    counts.append(c)
     701            terms = zip(counts, unique)
     702        self.terms = list(terms)
     703
     704    def __repr__(self):
     705        return " + ".join(("(%s)*%s" % t) for t in self.terms)
     706
     707    def n(self, **kwargs):
     708        return sum(L.n(**kwargs)*R.n(**kwargs) for (L,R) in self.terms)
     709
     710
     711class Function_PFQ(BuiltinFunction):
     712    r"""
     713    Gives the generalized hypergeometric function
     714
     715    .. math ::
     716
     717        \,_pF_q\left(\begin{matrix} a_1,a_2,\ldots,a_p \\
     718            b_1,b_2,\ldots,b_q \end{matrix} ; z \right) =
     719            \sum_{k=0}^{\infty} \frac{(a_1)_k \cdots (a_p)_k}
     720                {(b_1)_k \cdots (b_q)_k} \frac{z^k}{k!}
     721
     722    where $(c)_k = a (a+1) \cdots (a+k-1)$.
     723
     724    If a representation for the function value in terms of elementary
     725    functions can be found, that representation is returned.
     726
     727    Note: for technical reasons, PFQ() must currently be called
     728    with the single list of parameters as PFQ(p,q,a1,...,ap,b1,...,bq,z).
     729    The alias hyp([a1,...,ap],[b1,...,bq],z) can be used instead.
     730
     731    EXAMPLES ::
     732
     733        sage: var('a b c z')
     734        (a, b, c, z)
     735        sage: hyp([a,b],[c],z)
     736        PFQ(2, 1, a, b, c, z)
     737        sage: latex(hyp([a,b],[c],z))
     738        \,_2F_1\left(\begin{matrix} a,b \\ c \end{matrix} ; z \right)
     739        sage: hyp([-3,1/3],[-4],z)
     740        7/162*z^3 + 1/9*z^2 + 1/4*z + 1
     741        sage: hyp([1,2],[3],0)
     742        1
     743        sage: hyp([a,b],[c],0)
     744        1
     745        sage: hyp([],[],z)
     746        e^z
     747        sage: hyp([a],[],z)
     748        (-z + 1)^(-a)
     749        sage: hyp([1,1,2],[1,1],z)
     750        (z - 1)^(-2)
     751        sage: hyp([2,3],[1],x)
     752        -1/(x - 1)^3 + 3*x/(x - 1)^4
     753        sage: hyp([1/2],[3/2],-5)
     754        1/10*sqrt(pi)*sqrt(5)*erf(sqrt(5))
     755        sage: hyp([2],[5],3)
     756        4
     757        sage: hyp([2],[5],5)
     758        48/625*e^5 + 612/625
     759        sage: hyp([1/2,7/2],[3/2],z)
     760        1/sqrt(-z + 1) + 2/3*z/(-z + 1)^(3/2) + 1/5*z^2/(-z + 1)^(5/2)
     761        sage: hyp([1/2,1],[2],z)
     762        -2*(sqrt(-z + 1) - 1)/z
     763        sage: hyp([1,1],[2],z)
     764        -log(-z + 1)/z
     765        sage: hyp([1,1],[3],z)
     766        -2*((z - 1)*log(-z + 1)/z - 1)/z
     767        sage: hyp([1],[5],x).series(x,5)
     768        1 + 1/5*x + 1/30*x^2 + 1/210*x^3 + 1/1680*x^4 + Order(x^5)
     769        sage: hyp([1,2,3],[4,5,6],1/2).n()
     770        1.02573619590134
     771        sage: hyp([1,2,3],[4,5,6],1/2).n(digits=30)
     772        1.02573619590133865036584139535
     773        sage: hyp([5-3*I],[3/2,2+I,sqrt(2)],4+I).n()
     774        5.52605111678803 - 7.86331357527540*I
     775
     776
     777    """
     778
     779    def __init__(self):
     780        """
     781        EXAMPLES ::
     782
     783            sage: var('a b c z')
     784            (a, b, c, z)
     785            sage: PFQ(2,1,a,b,c,z)
     786            PFQ(2, 1, a, b, c, z)
     787
     788        """
     789        BuiltinFunction.__init__(self, 'PFQ', nargs=0,
     790            latex_name=None, evalf_params_first=False)
     791
     792    def _parse(self, args):
     793        """
     794        Internal hack to convert from symbolic function parameters
     795        list to (p,q,a,b,z) data.
     796
     797        TESTS ::
     798
     799            sage: from sage.functions.hypergeometric import Function_PFQ
     800            sage: Function_PFQ._parse(PFQ,[2,1,3,4,5,6])
     801            (2, 1, [3, 4], [5], 6)
     802
     803        """
     804        assert len(args) >= 3
     805        p = int(args[0])
     806        q = int(args[1])
     807        assert len(args) == p+q+3
     808        z = args[-1]
     809        a = args[2:p+2]
     810        b = args[p+2:p+q+2]
     811        return p, q, a, b, z
     812
     813    def _eval_(self, *args):
     814        """
     815        Try to simplify parameters or rewrite in closed form; otherwise
     816        return unevaluated.
     817
     818        EXAMPLES ::
     819
     820            sage: hyp([2],[1],2)
     821            3*e^2
     822            sage: hyp([2],[1/3],2)
     823            PFQ(1, 1, 2, 1/3, 2)
     824
     825        """
     826        p, q, a, b, z = self._parse(args)
     827        f_input = FormalPFQ(a, b, z)
     828        f = f_input.sort_parameters()
     829        if not f.is_termwise_finite():
     830            raise ValueError("the function is undefined with the given parameters and/or argument")
     831        if f.is_terminating():
     832            return f.add_terms()
     833        terms = f.deflate().terms
     834        s = 0
     835        for coeff, pfq in terms:
     836            g = pfq.closed_form()
     837            if isinstance(g, FormalPFQ):
     838                # No simplification was found, so may as well return
     839                # original PFQ
     840                if f == f_input:
     841                    return None
     842                else:
     843                    return f.to_PFQ()
     844            s += coeff * g
     845        return s
     846
     847    def _derivative_(self, *args, **kwargs):
     848        """
     849        EXAMPLES ::
     850
     851            sage: hyp([1/3,2/3],[5],x^2).diff(x)
     852            4/45*x*PFQ(2, 1, 4/3, 5/3, 6, x^2)
     853            sage: hyp([1,2],[x],2).diff(x)
     854            Traceback (most recent call last):
     855                ...
     856            NotImplementedError: derivative of hypergeometric function with respect to parameters
     857
     858        """
     859        diff_param = kwargs['diff_param']
     860        p, q, a, b, z = self._parse(args)
     861        if diff_param != len(args)-1:
     862            raise NotImplementedError("derivative of hypergeometric function with respect to parameters")
     863        t = Integer(1)
     864        for c in a: t *= c
     865        for c in b: t /= c
     866        return t * self(*([p,q] + [c+1 for c in a] + [c+1 for c in b] + [z]))
     867
     868    def _print_latex_(self, *args):
     869        """
     870        EXAMPLES ::
     871
     872            sage: var('a b c z')
     873            (a, b, c, z)
     874            sage: latex(hyp([a,b],[c],z))
     875            \,_2F_1\left(\begin{matrix} a,b \\ c \end{matrix} ; z \right)
     876
     877        """
     878        p, q, a, b, z = self._parse(args)
     879        return FormalPFQ(a,b,z)._print_latex_()
     880
     881    def _evalf_(self, *args, **kwargs):
     882        """
     883        EXAMPLES ::
     884
     885            sage: hyp([1,2,3],[4,5,6],1/2).n()
     886            1.02573619590134
     887            sage: hyp([1,2,3],[4,5,6],1/2).n(digits=25)
     888            1.025736195901338650365841
     889
     890        """
     891        parent = kwargs['parent']
     892        p, q, a, b, z = self._parse(args)
     893        return FormalPFQ(a,b,z).n(**kwargs)
     894
     895
     896
     897PFQ = Function_PFQ()
     898
     899
     900def hyp(a, b, z):
     901    """
     902    hyp([a1,...,ap],[b1,...,bq],z) gives the hypergeometric function
     903    $\,_pF_q$, and is equivalent to calling PFQ(p,q,a1,...,ap,b1,...,bq,z).
     904
     905    EXAMPLES ::
     906
     907        sage: hyp([],[5],1)
     908        PFQ(0, 1, 5, 1)
     909
     910    """
     911    return FormalPFQ(a,b,z).to_PFQ()
     912
     913
     914
     915
     916class BesselFunction(ForceSR):
     917    """
     918    Base class for symbolic Bessel functions.
     919    """
     920
     921    def __init__(self):
     922        """
     923
     924            sage: bessel_j(1,1)
     925            bessel_j(1, 1)
     926
     927        """
     928        BuiltinFunction.__init__(self,
     929            ('bessel_%s' % self.char.lower()), nargs=2)
     930
     931    def _print_latex_(self, nu, z):
     932        r"""
     933
     934            sage: latex(bessel_j(1,1))
     935            J_{1}\left(1\right)
     936
     937        """
     938        return r"%s_{%s}\left(%s\right)" % (self.char, latex(nu), latex(z))
     939
     940    def _evalf_(self, nu, z, **kwargs):
     941        """
     942
     943            sage: bessel_j(1,1).n()
     944            0.440050585744933
     945
     946        """
     947        import sage.libs.mpmath.all as a
     948        return a.call(getattr(a, "bessel"+self.char.lower()), nu, z, **kwargs)
     949
     950    def _series_(self, nu, z, var, **kwargs):
     951        """
     952
     953            sage: z = var('z')
     954            sage: bessel_j(1,z).series(z,5)
     955            1/2*z + (-1/16)*z^3 + Order(z^5)
     956
     957        """
     958        # XXX
     959        assert kwargs.get('at') == 0
     960        assert kwargs.get('options') == 0
     961        return self._expand_hyper(nu, z).series(var,kwargs['order'])
     962
     963
     964class Function_bessel_j(BesselFunction):
     965    r"""
     966    Gives the Bessel function of the first kind $J_{\nu}(z)$,
     967    with order (parameter) $\nu$ and argument $z$.
     968
     969    EXAMPLES::
     970
     971    Symbolic input and symmetry transformations::
     972
     973        sage: var('nu z')
     974        (nu, z)
     975        sage: bessel_j(nu, z)
     976        bessel_j(nu, z)
     977        sage: bessel_j(0, 1)
     978        bessel_j(0, 1)
     979        sage: bessel_j(nu,0)
     980        bessel_j(nu, 0)
     981        sage: bessel_j(0,z)
     982        bessel_j(0, z)
     983        sage: bessel_j(-3, 4)
     984        -bessel_j(3, 4)
     985        sage: bessel_j(3, -4)
     986        -bessel_j(3, 4)
     987        sage: bessel_j(-3, -4)
     988        bessel_j(3, 4)
     989        sage: bessel_j(-4, 5)
     990        bessel_j(4, 5)
     991
     992    Exact values and limits::
     993
     994        sage: bessel_j(0, 0)
     995        1
     996        sage: bessel_j(1, 0)
     997        0
     998        sage: bessel_j(2, 0)
     999        0
     1000        sage: bessel_j(-2, 0)
     1001        0
     1002        sage: bessel_j(nu, Infinity)
     1003        0
     1004        sage: bessel_j(nu, -Infinity)
     1005        0
     1006        sage: bessel_j(1/2,pi)
     1007        0
     1008        sage: bessel_j(1/2,-8*pi)
     1009        0
     1010
     1011    Numerical evaluation::
     1012
     1013        sage: bessel_j(0,1).n()
     1014        0.765197686557966
     1015        sage: bessel_j(0,1).n(digits=30)
     1016        0.765197686557966551449717526103
     1017        sage: bessel_j(-5, 125000000).n()
     1018        0.0000711791694771539
     1019        sage: bessel_j(2+3*I,4-5*I).n(digits=30)
     1020        1643.61099442312390496727290972 - 848.086110331983438476479306274*I
     1021
     1022    Symbolic and explicit series expansions::
     1023
     1024        sage: bessel_j(1,z).series(z,5)
     1025        1/2*z + (-1/16)*z^3 + Order(z^5)
     1026        sage: from sage.functions.hypergeometric import expand_hyper
     1027        sage: expand_hyper(bessel_j(nu,z))
     1028        (1/2*z)^nu*PFQ(0, 1, nu + 1, -1/4*z^2)/gamma(nu + 1)
     1029        sage: expand_hyper(bessel_j(0,z)).series(z,5)
     1030        1 + (-1/4)*z^2 + 1/64*z^4 + Order(z^5)
     1031        sage: expand_hyper(bessel_j(1,z)).series(z,7)
     1032        1/2*z + (-1/16)*z^3 + 1/384*z^5 + Order(z^7)
     1033
     1034    """
     1035
     1036    char = 'J'
     1037
     1038    def _eval_symbolic_(self, nu, z):
     1039        """
     1040        Evaluate or simplify for symbolic input.
     1041
     1042        EXAMPLES ::
     1043
     1044            sage: bessel_j(1, 0)
     1045            0
     1046
     1047        """
     1048        try:
     1049            n = nu._integer_()
     1050            if n < 0:
     1051                return (-1)**abs(n) * self(-n,z)
     1052            if z < -z:
     1053                return (-1)**abs(n) * self(n,-z)
     1054        except TypeError:
     1055            pass
     1056        if z == 0:
     1057            r = real_part(nu)
     1058            if r > 0:
     1059                return z
     1060            if r < 0:
     1061                return Infinity
     1062            if nu == 0:
     1063                return 1
     1064        if z == Infinity or z == -Infinity:
     1065            return 0
     1066        if 2*nu == 1:
     1067            if (z / pi)._is_integer():
     1068                return 0
     1069        return None
     1070
     1071    def _expand_hyper(self, nu, z):
     1072        """
     1073        Expand as a hypergeometric series.
     1074
     1075            sage: from sage.functions.hypergeometric import expand_hyper
     1076            sage: var('z')
     1077            z
     1078            sage: expand_hyper(bessel_j(1,z))
     1079            1/2*z*PFQ(0, 1, 2, -1/4*z^2)
     1080
     1081        """
     1082        return (z/2)**nu / gamma(nu+1) * hyp([],[nu+1],-(z/2)**2)
     1083
     1084
     1085class Function_bessel_y(BesselFunction):
     1086    """
     1087    Gives the Bessel function of the second kind $Y_{\nu}(z)$,
     1088    with order (parameter) $\nu$ and argument $z$.
     1089    """
     1090
     1091    char = 'Y'
     1092
     1093    def _eval_symbolic_(self, nu, z):
     1094        """
     1095        Evaluate or simplify for symbolic input.
     1096
     1097        EXAMPLES ::
     1098
     1099            sage: bessel_y(1,1)
     1100            bessel_y(1, 1)
     1101
     1102        """
     1103        return None
     1104
     1105    def _expand_hyper(self, nu, z):
     1106        """
     1107        Expand as a hypergeometric series (this only makes sense when
     1108        nu is not an integer).
     1109
     1110            sage: from sage.functions.hypergeometric import expand_hyper
     1111            sage: var('nu z')
     1112            (nu, z)
     1113            sage: expand_hyper(bessel_y(nu,z))
     1114            (pi*(1/2*z)^nu*nu*PFQ(0, 1, nu + 1, -1/4*z^2)/gamma(nu + 1) - (1/2*z)^(-nu)*PFQ(0, 1, -nu + 1, -1/4*z^2)/gamma(-nu + 1))/(pi*nu)
     1115            sage: expand_hyper(bessel_y(1/2,z))
     1116            -2*sqrt(1/2)*cosh(sqrt(-z^2))/(sqrt(pi)*sqrt(z))
     1117            sage: expand_hyper(bessel_y(1,z))
     1118            Traceback (most recent call last):
     1119               ...
     1120            ZeroDivisionError: Symbolic division by zero
     1121
     1122        """
     1123        f = (bessel_j(nu,z)*cos(nu*pi) - bessel_j(-nu,z)) / sin(nu*pi)
     1124        return expand_hyper(f)
     1125
     1126
     1127class Function_bessel_i(BesselFunction):
     1128    r"""
     1129    Gives the modified Bessel function of the first kind $I_{\nu}(z)$,
     1130    with order (parameter) $\nu$ and argument $z$.
     1131
     1132    EXAMPLES::
     1133
     1134    Symbolic input and symmetry transformations::
     1135
     1136        sage: var('nu z')
     1137        (nu, z)
     1138        sage: bessel_i(nu, z)
     1139        bessel_i(nu, z)
     1140        sage: bessel_i(0,1)
     1141        bessel_i(0, 1)
     1142        sage: bessel_i(nu,0)
     1143        bessel_i(nu, 0)
     1144        sage: bessel_i(0,z)
     1145        bessel_i(0, z)
     1146        sage: bessel_i(-3,4)
     1147        bessel_i(3, 4)
     1148        sage: bessel_i(3,-4)
     1149        -bessel_i(3, 4)
     1150        sage: bessel_i(-3,-4)
     1151        -bessel_i(3, 4)
     1152        sage: bessel_i(-4,5)
     1153        bessel_i(4, 5)
     1154
     1155    Exact values and limits::
     1156
     1157        sage: bessel_i(0, 0)
     1158        1
     1159        sage: bessel_i(1, 0)
     1160        0
     1161        sage: bessel_i(2, 0)
     1162        0
     1163        sage: bessel_i(-2, 0)
     1164        0
     1165        sage: bessel_i(1, Infinity)
     1166        +Infinity
     1167        sage: bessel_i(1, -Infinity)
     1168        -Infinity
     1169        sage: bessel_i(2, Infinity)
     1170        +Infinity
     1171        sage: bessel_i(2, -Infinity)
     1172        +Infinity
     1173        sage: bessel_i(1/2,I*pi)
     1174        0
     1175        sage: bessel_i(1/2,-8*I*pi)
     1176        0
     1177
     1178    Numerical evaluation::
     1179
     1180        sage: bessel_j(0,1).n()
     1181        0.765197686557966
     1182        sage: bessel_j(0,1).n(digits=30)
     1183        0.765197686557966551449717526103
     1184        sage: bessel_j(-5, 125000000).n()
     1185        0.0000711791694771539
     1186        sage: bessel_i(2+3*I, 4-5*I).n(digits=30)
     1187        -6.69838203392827123593220944996 + 13.1635256869717070482549506439*I
     1188
     1189    Symbolic and explicit series expansions::
     1190
     1191        sage: from sage.functions.hypergeometric import expand_hyper
     1192        sage: expand_hyper(bessel_i(nu,z))
     1193        (1/2*z)^nu*PFQ(0, 1, nu + 1, 1/4*z^2)/gamma(nu + 1)
     1194        sage: expand_hyper(bessel_i(0,z)).series(z,5)
     1195        1 + 1/4*z^2 + 1/64*z^4 + Order(z^5)
     1196        sage: expand_hyper(bessel_i(1,z)).series(z,7)
     1197        1/2*z + 1/16*z^3 + 1/384*z^5 + Order(z^7)
     1198
     1199    """
     1200
     1201    char = 'I'
     1202
     1203    def _eval_symbolic_(self, nu, z):
     1204        """
     1205        Evaluate or simplify for symbolic input.
     1206
     1207        EXAMPLES ::
     1208
     1209            sage: bessel_i(0,0)
     1210            1
     1211
     1212        """
     1213        try:
     1214            n = nu._integer_()
     1215            if n < 0:
     1216                return self(-n,z)
     1217            if z < -z:
     1218                return (-1)**abs(n) * self(n,-z)
     1219        except TypeError:
     1220            pass
     1221        if z == 0:
     1222            r = real_part(nu)
     1223            if r > 0:
     1224                return z
     1225            if r < 0:
     1226                return Infinity
     1227            if nu == 0:
     1228                return 1
     1229        if z == Infinity:
     1230            return Infinity
     1231        if z == -Infinity:
     1232            return (-1)**abs(nu) * Infinity
     1233        if 2*nu == 1:
     1234            if (z / (pi*I))._is_integer():
     1235                return 0
     1236        return None
     1237
     1238    def _expand_hyper(self, nu, z):
     1239        """
     1240        Expand as a hypergeometric series.
     1241
     1242            sage: from sage.functions.hypergeometric import expand_hyper
     1243            sage: var('z')
     1244            z
     1245            sage: expand_hyper(bessel_i(1,z))
     1246            1/2*z*PFQ(0, 1, 2, 1/4*z^2)
     1247
     1248        """
     1249        return (z/2)**nu / gamma(nu+1) * hyp([],[nu+1],(z/2)**2)
     1250
     1251class Function_bessel_k(BesselFunction):
     1252    r"""
     1253    Gives the modified Bessel function of the second kind $K_{\nu}(z)$,
     1254    with order (parameter) $\nu$ and argument $z$.
     1255    """
     1256
     1257    char = 'K'
     1258
     1259    def _eval_symbolic_(self, nu, z):
     1260        """
     1261        Evaluate or simplify for symbolic input.
     1262
     1263        EXAMPLES ::
     1264
     1265            sage: bessel_k(1,1)
     1266            bessel_k(1, 1)
     1267
     1268        """
     1269        return None
     1270
     1271
     1272
     1273bessel_j = Function_bessel_j()
     1274bessel_y = Function_bessel_y()
     1275bessel_i = Function_bessel_i()
     1276bessel_k = Function_bessel_k()