# Ticket #2516: hyper.patch

File hyper.patch, 40.1 KB (added by fredrik.johansson, 12 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 from piecewise import piecewise, Piecewise from trig import ( sin, cos, sec, csc, cot, tan, from trig import ( sin, cos, sec, csc, cot, tan, asin, acos, atan, acot, acsc, asec, arcsin, arccos, arctan, arcsin, arccos, arctan, arccot, arccsc, arcsec, arctan2, atan2) from hyperbolic import ( tanh, sinh, cosh, coth, sech, csch, asinh, acosh, atanh, acoth, asech, acsch, asinh, acosh, atanh, acoth, asech, acsch, arcsinh, arccosh, arctanh, arccoth, arcsech, arccsch ) reciprocal_trig_functions = {'sec': cos, 'csc': sin, 'cot': tan, 'sech': cosh, 'csch': sinh, 'coth': tanh} from transcendental import (exponential_integral_1, zeta, zetaderiv, zeta_symmetric, Li, Ei, Li, Ei, dickman_rho) from special import (bessel_I, bessel_J, bessel_K, bessel_Y, elliptic_f, elliptic_ec, elliptic_eu, elliptic_kc, elliptic_pi, elliptic_j, airy_ai, airy_bi) from orthogonal_polys import (chebyshev_T, chebyshev_U, gen_laguerre, from wigner import (wigner_3j, clebsch_gordan, racah, wigner_6j, wigner_9j, gaunt) from generalized import (dirac_delta, heaviside, unit_step, sgn, sign, from generalized import (dirac_delta, heaviside, unit_step, sgn, sign, kronecker_delta) from min_max import max_symbolic, min_symbolic from hypergeometric import (FormalPFQ, PFQ, hyp, bessel_i, bessel_j, bessel_k, bessel_y)
• ## new file sage/functions/hypergeometric.py

diff -r 2ecc6b33a446 -r 7d20ef3a988c sage/functions/hypergeometric.py
 - """ This module implements manipulation of infinite hypergeometric series represented in standard parametric form (as $\,_pFq$ functions). Also included are some hypergeometric-type special functions (currently, only Bessel functions). """ from sage.rings.all import Integer, ZZ, QQ, Infinity, binomial, rising_factorial from sage.rings.all import factorial as fac from sage.functions.other import sqrt, gamma, real_part from sage.functions.log import exp, log from sage.functions.trig import cos, sin from sage.functions.hyperbolic import cosh, sinh from sage.functions.other import erf from sage.symbolic.constants import pi from sage.symbolic.all import I from sage.symbolic.function import BuiltinFunction from sage.symbolic.ring import SR from sage.misc.latex import latex def rational_param_as_tuple(x): """ Utility function for converting rational PFQ parameters to tuples (which mpmath handles more efficiently). sage: from sage.functions.hypergeometric import rational_param_as_tuple as r sage: r(1/2) (1, 2) sage: r(3) 3 sage: r(pi) pi """ try: x = x.pyobject() except (TypeError, AttributeError): pass try: if x.parent() is QQ: p = int(x.numer()) q = int(x.denom()) return p, q except (TypeError, AttributeError): pass return x def expand_hyper(x): """ Expand compatible functions as hypergeometric series. EXAMPLES :: sage: from sage.functions.hypergeometric import expand_hyper sage: expand_hyper(bessel_j(2,x)) 1/8*x^2*PFQ(0, 1, 3, -1/4*x^2) """ try: op = x.operator() except AttributeError: f = SR(x) op = x.operator() if hasattr(op, "_expand_hyper"): return op._expand_hyper(*x.operands()) if op is None: return x return reduce(op, [expand_hyper(t) for t in x.operands()]) class ForceSR(BuiltinFunction): """ Force inputs to SR XXX: this does not work as intended because of the way BuiltinFunction.__call__ works. """ def _eval_(self, *args): def __force_SR(x): # use py_scalar_parent? if isinstance(x, (int, long)): return SR(Integer(x)) else: return SR(x) args = [__force_SR(a) for a in args] value = self._eval_symbolic_(*args) if value is None: return None g = __force_SR(value) return g class FormalPFQ: r""" Represents a (formal) generalized infinite hypergeometric series. """ def __init__(self, a, b, z): """ EXAMPLES :: sage: FormalPFQ([],[],1) 0F0([]; []; 1) sage: FormalPFQ([],[1],1) 0F1([]; [1]; 1) sage: FormalPFQ([2,3],[1],1) 2F1([2, 3]; [1]; 1) """ self.a = list(a) self.b = list(b) self.z = z def __eq__(self, other): """ TESTS :: sage: FormalPFQ([],[],1) == FormalPFQ([],[],1) True sage: FormalPFQ([2],[],1) == FormalPFQ([],[],1) False """ if not isinstance(other, FormalPFQ): return False return self.a == other.a and self.b == other.b and self.z == other.z def __ne__(self, other): """ TESTS :: sage: FormalPFQ([],[],1) != FormalPFQ([],[],1) False """ return not self.__eq__(other) def __hash__(self): """ TESTS :: sage: a = FormalPFQ([1],[2],3) sage: b = FormalPFQ([1],[2],3) sage: hash(a) == hash(b) True """ return hash((tuple(self.a), tuple(self.b), self.z)) @property def p(self): """ Return the numerator degree of self. EXAMPLES :: sage: FormalPFQ([],[],1).p 0 sage: FormalPFQ([1],[],1).p 1 sage: FormalPFQ([1,1],[],1).p 2 """ return len(self.a) @property def q(self): """ Return the denominator degree of self. EXAMPLES :: sage: FormalPFQ([],[],1).q 0 sage: FormalPFQ([],[1],1).q 1 sage: FormalPFQ([],[1,1],1).q 2 """ return len(self.b) def __repr__(self): """ TESTS :: sage: FormalPFQ([1],[2],3) 1F1([1]; [2]; 3) """ return "%iF%i(%s; %s; %s)" % (self.p, self.q, self.a, self.b, self.z) def to_PFQ(self): """ Convert to a PFQ (symbolic function) instance. EXAMPLES :: sage: FormalPFQ([1,2],[3,4],5).to_PFQ() PFQ(2, 2, 1, 2, 3, 4, 5) """ return PFQ(*([self.p, self.q] + self.a + self.b + [self.z])) def _print_latex_(self): r""" TESTS :: sage: latex(FormalPFQ([1,1],[2],-1)) \texttt{2F1([1, 1]; [2]; -1)} """ p = self.p q = self.q a = ",".join(latex(c) for c in self.a) b = ",".join(latex(c) for c in self.b) z = latex(self.z) return r"\,_%iF_%i\left(\begin{matrix} %s \\ %s \end{matrix} ; %s \right)" % (p, q, a, b, z) def n(self, **kwargs): """ TESTS :: sage: FormalPFQ([1,1],[2],-1).n() 0.693147180559945 """ from sage.libs.mpmath.all import hyper, call a = [rational_param_as_tuple(c) for c in self.a] b = [rational_param_as_tuple(c) for c in self.b] return call(hyper, a, b, self.z, **kwargs) def sort_parameters(self): """ Sort parameters in a canonical order. EXAMPLES:: sage: FormalPFQ([2,1,3],[5,4],1/2).sort_parameters() 3F2([1, 2, 3]; [4, 5]; 1/2) """ return FormalPFQ(sorted(self.a), sorted(self.b), self.z) def eliminate_parameters(self): """ Eliminate repeated parameters. EXAMPLES:: sage: FormalPFQ([1,1,2,5],[5,1,4],1/2).eliminate_parameters() 2F1([1, 2]; [4]; 1/2) """ aa = self.a[:] bb = self.b[:] p = self.p q = self.q i = 0 while i < q and aa: b = bb[i] if b in aa: aa.remove(b) bb.remove(b) p -= 1 q -= 1 else: i += 1 if (p,q) != (self.p,self.q): return FormalPFQ(aa, bb, self.z) return self def is_termwise_finite(self): """ Determine whether all terms of self are finite. Any infinite terms or ambiguous terms beyond the first zero, if one exists, are ignored. Ambiguous cases (where a term is the product of both zero and an infinity) are not considered finite. sage: FormalPFQ([2],[3,4],5).is_termwise_finite() True sage: FormalPFQ([2],[-3,4],5).is_termwise_finite() False sage: FormalPFQ([-2],[-3,4],5).is_termwise_finite() True sage: FormalPFQ([-3],[-3,4],5).is_termwise_finite()  # ambiguous False sage: FormalPFQ([0],[-1],5).is_termwise_finite() True sage: FormalPFQ([0],[0],5).is_termwise_finite()   # ambiguous False sage: FormalPFQ([1],[2],Infinity).is_termwise_finite() False sage: FormalPFQ([0],[0],Infinity).is_termwise_finite()  # ambiguous False sage: FormalPFQ([0],[],Infinity).is_termwise_finite()   # ambiguous False """ if self.z == 0: return 0 not in self.b if abs(self.z) == Infinity: return False if abs(self.z) == Infinity: return False for b in self.b: if b in ZZ and b <= 0: if any((a in ZZ) and (b < a <= 0) for a in self.a): continue return False return True def is_terminating(self): """ Determine whether the series represented by self terminates after a finite number of terms, i.e. whether any of the numerator parameters are nonnegative integers (with no preceding nonnegative denominator parameters), or z = 0. If terminating, the series represents a polynomial of self.z EXAMPLES :: sage: z = var('z') sage: FormalPFQ([1,2],[3,4],z).is_terminating() False sage: FormalPFQ([1,-2],[3,4],z).is_terminating() True sage: FormalPFQ([1,-2],[],z).is_terminating() True sage: FormalPFQ([2,3],[],0).is_terminating() True """ if self.z == 0: return True for a in self.a: if (a in ZZ) and (a <= 0): return self.is_termwise_finite() return False def is_absolutely_convergent(self): """ Determine whether self converges absolutely as an infinite series. False is returned if not all terms are finite. EXAMPLES :: Degree giving infinite radius of convergence:: sage: FormalPFQ([2,3],[4,5],6).is_absolutely_convergent() True sage: FormalPFQ([2,3],[-4,5],6).is_absolutely_convergent()  # undefined False sage: FormalPFQ([2,3],[-4,5],Infinity).is_absolutely_convergent()  # undefined False Ordinary geometric series (unit radius of convergence):: sage: FormalPFQ([1],[],1/2).is_absolutely_convergent() True sage: FormalPFQ([1],[],2).is_absolutely_convergent() False sage: FormalPFQ([1],[],1).is_absolutely_convergent() False sage: FormalPFQ([1],[],-1).is_absolutely_convergent() False sage: FormalPFQ([1],[],-1).n()    # Sum still exists 0.500000000000000 Degree $p = q+1$ (unit radius of convergence):: sage: FormalPFQ([2,3],[4],6).is_absolutely_convergent() False sage: FormalPFQ([2,3],[4],1).is_absolutely_convergent() False sage: FormalPFQ([2,3],[5],1).is_absolutely_convergent() False sage: FormalPFQ([2,3],[6],1).is_absolutely_convergent() True sage: FormalPFQ([-2,3],[4],5).is_absolutely_convergent() True sage: FormalPFQ([2,-3],[4],5).is_absolutely_convergent() True sage: FormalPFQ([2,-3],[-4],5).is_absolutely_convergent() True sage: FormalPFQ([2,-3],[-1],5).is_absolutely_convergent() False Degree giving zero radius of convergence:: sage: FormalPFQ([1,2,3],[4],2).is_absolutely_convergent() False sage: FormalPFQ([1,2,3],[4],1/2).is_absolutely_convergent() False sage: FormalPFQ([1,2,3],[4],0).is_absolutely_convergent() True sage: FormalPFQ([1,2,-3],[4],1/2).is_absolutely_convergent()  # polynomial True """ if not self.is_termwise_finite(): return False if self.p <= self.q: return True if self.is_terminating(): return True if self.p == self.q+1: if abs(self.z) < 1: return True if abs(self.z) == 1: if real_part(sum(self.b) - sum(self.a)) > 0: return True return False def terms(self, n=None): """ Generate the terms of self (optionally only n terms). EXAMPLES :: sage: list(FormalPFQ([-2,1],[3,4],x).terms()) [1, -1/6*x, 1/120*x^2] sage: list(FormalPFQ([-2,1],[3,4],x).terms(2)) [1, -1/6*x] sage: list(FormalPFQ([-2,1],[3,4],x).terms(0)) [] """ if n is None: n = Infinity t = Integer(1) k = 1 while k <= n: yield t for a in self.a: t *= (a+k-1) for b in self.b: t /= (b+k-1) t *= self.z if t == 0: break t /= k k += 1 def add_terms(self, n=None): """ Add the first n terms of self. If n=None, self must be terminating and all terms will be expanded. EXAMPLES:: sage: z = var('z') sage: FormalPFQ([],[],z).add_terms(0) 0 sage: FormalPFQ([],[],z).add_terms(1) 1 sage: FormalPFQ([],[],z).add_terms(2) z + 1 sage: FormalPFQ([],[],z).add_terms(3) 1/2*z^2 + z + 1 sage: FormalPFQ([],[],z).add_terms() Traceback (most recent call last): ... ValueError: self must be terminating sage: FormalPFQ([-2],[],z).add_terms() z^2 - 2*z + 1 sage: FormalPFQ([-2],[],z).add_terms(3) z^2 - 2*z + 1 sage: FormalPFQ([-2],[],z).add_terms(6) z^2 - 2*z + 1 sage: FormalPFQ([-2],[],z).add_terms(2) -2*z + 1 sage: FormalPFQ([],[],z).add_terms(6) 1/120*z^5 + 1/24*z^4 + 1/6*z^3 + 1/2*z^2 + z + 1 sage: FormalPFQ([1],[],z).add_terms(6) z^5 + z^4 + z^3 + z^2 + z + 1 sage: FormalPFQ([],[1/2],-z^2/4).add_terms(6) -1/3628800*z^10 + 1/40320*z^8 - 1/720*z^6 + 1/24*z^4 - 1/2*z^2 + 1 sage: FormalPFQ([1,2],[3],1/3).add_terms(6).n() 1.29788359788360 sage: FormalPFQ([1,2],[3],1/3).n() 1.29837194594696 """ if n is None: if not self.is_terminating(): raise ValueError("self must be terminating") return sum(self.terms(n)) def deflate(self): """ Rewrite as a linear combination of functions of strictly lower degree by eliminating all parameters a[i] and b[j] such that a[i] = b[i] + m for nonnegative integer m. EXAMPLES :: sage: x = FormalPFQ([5],[4],3) sage: y = x.deflate() sage: y (7/4)*0F0([]; []; 3) sage: x.n(); y.n() 35.1496896155784 35.1496896155784 sage: x = FormalPFQ([6,1],[3,4,5],10) sage: y = x.deflate() sage: y (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) sage: x.n(); y.n() 2.87893612686782 2.87893612686782 sage: x = FormalPFQ([6,7],[3,4,5],10) sage: y = x.deflate() sage: y (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) sage: x.n(); y.n() 63.0734110716969 63.0734110716969 """ self = self.eliminate_parameters() for i,a in enumerate(self.a): for j,b in enumerate(self.b): m = a-b if m in ZZ and m > 0: aa = self.a[:i] + self.a[i+1:] bb = self.b[:j] + self.b[j+1:] terms = [] for k in xrange(m+1): # TODO: could rewrite prefactors as recurrence term = binomial(m,k) for c in aa: term *= rising_factorial(c,k) for c in bb: term /= rising_factorial(c,k) term *= self.z**k term /= rising_factorial(a-m, k) F = FormalPFQ([c+k for c in aa], [c+k for c in bb], self.z) Fterms = F.deflate().terms terms += [(term*termG, G) for (termG, G) in Fterms] return FormalPFQLinearCombination(terms) return FormalPFQLinearCombination([(1,self)]) def closed_form(self): """ Try to evaluate self in closed form using elementary (and other simple) functions. It may be necessary to call .deflate() first to find some closed forms. EXAMPLES :: sage: z = var('z') sage: FormalPFQ([1],[],1+z).closed_form() -1/z sage: FormalPFQ([],[],1+z).closed_form() e^(z + 1) sage: FormalPFQ([],[1/2],4).closed_form() cosh(4) sage: FormalPFQ([],[3/2],4).closed_form() 1/4*sinh(4) sage: FormalPFQ([],[5/2],4).closed_form() -3/64*sinh(4) + 3/16*cosh(4) sage: FormalPFQ([],[-3/2],4).closed_form() -4*sinh(4) + 19/3*cosh(4) sage: FormalPFQ([-3,1],[var('a')],z).closed_form() -6*z^3/((a + 1)*(a + 2)*a) + 6*z^2/((a + 1)*a) - 3*z/a + 1 """ if self.is_terminating(): return self.add_terms() # necessary? new = self.sort_parameters() new = self.eliminate_parameters() a, b, z, p, q = new.a, new.b, new.z, new.p, new.q if z == 0: return Integer(1) if p == q == 0: return exp(z) if p == 1 and q == 0: return (1-z)**(-a[0]) if p == 0 and q == 1: # TODO: make this require only linear time def _0f1(b,z): F12 = cosh(2*sqrt(z)) F32 = sinh(2*sqrt(z))/(2*sqrt(z)) if 2*b == 1: return F12 if 2*b == 3: return F32 if 2*b > 3: return (b-2)*(b-1)/z * (_0f1(b-2,z) - _0f1(b-1,z)) if 2*b < 1: return _0f1(b+1,z) + z/(b*(b+1)) * _0f1(b+2,z) raise ValueError # Can evaluate 0F1 in terms of elementary functions when # the parameter is a half-integer if 2*b[0] in ZZ and b[0] not in ZZ: return _0f1(b[0], z) # Confluent hypergeometric function if p == 1 and q == 1: aa, bb = a[0], b[0] if aa*2 == 1 and bb*2 == 3: t = sqrt(-z) return sqrt(pi)/2*erf(t)/t if a == 1 and b == 2: return (exp(z)-1)/z n, m = aa, bb if n in ZZ and m in ZZ and m > 0 and n > 0: rf = rising_factorial if m <= n: return exp(z)*sum(rf(m-n,k)*(-z)**k/fac(k)/rf(m,k) for k in xrange(n-m+1)) else: T = sum(rf(n-m+1,k)*z**k/(fac(k)*rf(2-m,k)) for k in xrange(m-n)) U = sum(rf(1-n,k)*(-z)**k/(fac(k)*rf(2-m,k)) for k in xrange(n)) return fac(m-2)*rf(1-m,n)*z**(1-m)/fac(n-1)*(T - exp(z)*U) if p == 2 and q == 1: R12 = QQ('1/2') R32 = QQ('3/2') def _2f1(a,b,c,z): """ Evaluation of 2F1(a,b,c,z), assuming a,b,c positive integers or half-integers """ if b == c: return (1-z)**(-a) if a == c: return (1-z)**(-b) if a == 0 or b == 0: return Integer(1) if a > b: a, b = b, a if b >= 2: F1 = _2f1(a,b-1,c,z) F2 = _2f1(a,b-2,c,z) q = (b-1)*(z-1) return ((c-2*b+2+(b-a-1)*z)*F1 + (b-c-1)*F2) / q if c > 2: # how to handle this case? if a-c+1 == 0 or b-c+1 == 0: raise NotImplementedError F1 = _2f1(a,b,c-1,z) F2 = _2f1(a,b,c-2,z) r1 = (c-1)*(2-c-(a+b-2*c+3)*z) r2 = (c-1)*(c-2)*(1-z) q = (a-c+1)*(b-c+1)*z return (r1*F1 + r2*F2) / q if (a,b,c) == (R12,1,2): return (2-2*sqrt(1-z))/z if (a,b,c) == (1,1,2): return -log(1-z)/z if (a,b,c) == (1,R32,R12): return (1+z)/(1-z)**2 if (a,b,c) == (1,R32,2): return 2*(1/sqrt(1-z)-1)/z if (a,b,c) == (R32,2,R12): return (1+3*z)/(1-z)**3 if (a,b,c) == (R32,2,1): return (2+z)/(2*(sqrt(1-z)*(1-z)**2)) if (a,b,c) == (2,2,1): return (1+z)/(1-z)**3 raise NotImplementedError aa, bb = a; cc, = b if z == 1: return gamma(cc)*gamma(cc-aa-bb)/gamma(cc-aa)/gamma(cc-bb) if (aa*2) in ZZ and (bb*2) in ZZ and (cc*2) in ZZ and \ aa > 0 and bb > 0 and cc > 0: try: return _2f1(aa,bb,cc,z) except NotImplementedError: pass return self class FormalPFQLinearCombination: """ Holds a weighted sum of FormalPFQs. """ def __init__(self, terms=[], collect=True): """ sage: from sage.functions.hypergeometric import FormalPFQLinearCombination sage: FormalPFQLinearCombination([(3,FormalPFQ([],[1],x))]) (3)*0F1([]; [1]; x) """ if collect: unique = [] counts = [] for c,f in terms: if f in unique: counts[unique.index(f)] += c else: unique.append(f) counts.append(c) terms = zip(counts, unique) self.terms = list(terms) def __repr__(self): return " + ".join(("(%s)*%s" % t) for t in self.terms) def n(self, **kwargs): return sum(L.n(**kwargs)*R.n(**kwargs) for (L,R) in self.terms) class Function_PFQ(BuiltinFunction): r""" Gives the generalized hypergeometric function .. math :: \,_pF_q\left(\begin{matrix} a_1,a_2,\ldots,a_p \\ b_1,b_2,\ldots,b_q \end{matrix} ; z \right) = \sum_{k=0}^{\infty} \frac{(a_1)_k \cdots (a_p)_k} {(b_1)_k \cdots (b_q)_k} \frac{z^k}{k!} where $(c)_k = a (a+1) \cdots (a+k-1)$. If a representation for the function value in terms of elementary functions can be found, that representation is returned. Note: for technical reasons, PFQ() must currently be called with the single list of parameters as PFQ(p,q,a1,...,ap,b1,...,bq,z). The alias hyp([a1,...,ap],[b1,...,bq],z) can be used instead. EXAMPLES :: sage: var('a b c z') (a, b, c, z) sage: hyp([a,b],[c],z) PFQ(2, 1, a, b, c, z) sage: latex(hyp([a,b],[c],z)) \,_2F_1\left(\begin{matrix} a,b \\ c \end{matrix} ; z \right) sage: hyp([-3,1/3],[-4],z) 7/162*z^3 + 1/9*z^2 + 1/4*z + 1 sage: hyp([1,2],[3],0) 1 sage: hyp([a,b],[c],0) 1 sage: hyp([],[],z) e^z sage: hyp([a],[],z) (-z + 1)^(-a) sage: hyp([1,1,2],[1,1],z) (z - 1)^(-2) sage: hyp([2,3],[1],x) -1/(x - 1)^3 + 3*x/(x - 1)^4 sage: hyp([1/2],[3/2],-5) 1/10*sqrt(pi)*sqrt(5)*erf(sqrt(5)) sage: hyp([2],[5],3) 4 sage: hyp([2],[5],5) 48/625*e^5 + 612/625 sage: hyp([1/2,7/2],[3/2],z) 1/sqrt(-z + 1) + 2/3*z/(-z + 1)^(3/2) + 1/5*z^2/(-z + 1)^(5/2) sage: hyp([1/2,1],[2],z) -2*(sqrt(-z + 1) - 1)/z sage: hyp([1,1],[2],z) -log(-z + 1)/z sage: hyp([1,1],[3],z) -2*((z - 1)*log(-z + 1)/z - 1)/z sage: hyp([1],[5],x).series(x,5) 1 + 1/5*x + 1/30*x^2 + 1/210*x^3 + 1/1680*x^4 + Order(x^5) sage: hyp([1,2,3],[4,5,6],1/2).n() 1.02573619590134 sage: hyp([1,2,3],[4,5,6],1/2).n(digits=30) 1.02573619590133865036584139535 sage: hyp([5-3*I],[3/2,2+I,sqrt(2)],4+I).n() 5.52605111678803 - 7.86331357527540*I """ def __init__(self): """ EXAMPLES :: sage: var('a b c z') (a, b, c, z) sage: PFQ(2,1,a,b,c,z) PFQ(2, 1, a, b, c, z) """ BuiltinFunction.__init__(self, 'PFQ', nargs=0, latex_name=None, evalf_params_first=False) def _parse(self, args): """ Internal hack to convert from symbolic function parameters list to (p,q,a,b,z) data. TESTS :: sage: from sage.functions.hypergeometric import Function_PFQ sage: Function_PFQ._parse(PFQ,[2,1,3,4,5,6]) (2, 1, [3, 4], [5], 6) """ assert len(args) >= 3 p = int(args[0]) q = int(args[1]) assert len(args) == p+q+3 z = args[-1] a = args[2:p+2] b = args[p+2:p+q+2] return p, q, a, b, z def _eval_(self, *args): """ Try to simplify parameters or rewrite in closed form; otherwise return unevaluated. EXAMPLES :: sage: hyp([2],[1],2) 3*e^2 sage: hyp([2],[1/3],2) PFQ(1, 1, 2, 1/3, 2) """ p, q, a, b, z = self._parse(args) f_input = FormalPFQ(a, b, z) f = f_input.sort_parameters() if not f.is_termwise_finite(): raise ValueError("the function is undefined with the given parameters and/or argument") if f.is_terminating(): return f.add_terms() terms = f.deflate().terms s = 0 for coeff, pfq in terms: g = pfq.closed_form() if isinstance(g, FormalPFQ): # No simplification was found, so may as well return # original PFQ if f == f_input: return None else: return f.to_PFQ() s += coeff * g return s def _derivative_(self, *args, **kwargs): """ EXAMPLES :: sage: hyp([1/3,2/3],[5],x^2).diff(x) 4/45*x*PFQ(2, 1, 4/3, 5/3, 6, x^2) sage: hyp([1,2],[x],2).diff(x) Traceback (most recent call last): ... NotImplementedError: derivative of hypergeometric function with respect to parameters """ diff_param = kwargs['diff_param'] p, q, a, b, z = self._parse(args) if diff_param != len(args)-1: raise NotImplementedError("derivative of hypergeometric function with respect to parameters") t = Integer(1) for c in a: t *= c for c in b: t /= c return t * self(*([p,q] + [c+1 for c in a] + [c+1 for c in b] + [z])) def _print_latex_(self, *args): """ EXAMPLES :: sage: var('a b c z') (a, b, c, z) sage: latex(hyp([a,b],[c],z)) \,_2F_1\left(\begin{matrix} a,b \\ c \end{matrix} ; z \right) """ p, q, a, b, z = self._parse(args) return FormalPFQ(a,b,z)._print_latex_() def _evalf_(self, *args, **kwargs): """ EXAMPLES :: sage: hyp([1,2,3],[4,5,6],1/2).n() 1.02573619590134 sage: hyp([1,2,3],[4,5,6],1/2).n(digits=25) 1.025736195901338650365841 """ parent = kwargs['parent'] p, q, a, b, z = self._parse(args) return FormalPFQ(a,b,z).n(**kwargs) PFQ = Function_PFQ() def hyp(a, b, z): """ hyp([a1,...,ap],[b1,...,bq],z) gives the hypergeometric function $\,_pF_q$, and is equivalent to calling PFQ(p,q,a1,...,ap,b1,...,bq,z). EXAMPLES :: sage: hyp([],[5],1) PFQ(0, 1, 5, 1) """ return FormalPFQ(a,b,z).to_PFQ() class BesselFunction(ForceSR): """ Base class for symbolic Bessel functions. """ def __init__(self): """ sage: bessel_j(1,1) bessel_j(1, 1) """ BuiltinFunction.__init__(self, ('bessel_%s' % self.char.lower()), nargs=2) def _print_latex_(self, nu, z): r""" sage: latex(bessel_j(1,1)) J_{1}\left(1\right) """ return r"%s_{%s}\left(%s\right)" % (self.char, latex(nu), latex(z)) def _evalf_(self, nu, z, **kwargs): """ sage: bessel_j(1,1).n() 0.440050585744933 """ import sage.libs.mpmath.all as a return a.call(getattr(a, "bessel"+self.char.lower()), nu, z, **kwargs) def _series_(self, nu, z, var, **kwargs): """ sage: z = var('z') sage: bessel_j(1,z).series(z,5) 1/2*z + (-1/16)*z^3 + Order(z^5) """ # XXX assert kwargs.get('at') == 0 assert kwargs.get('options') == 0 return self._expand_hyper(nu, z).series(var,kwargs['order']) class Function_bessel_j(BesselFunction): r""" Gives the Bessel function of the first kind $J_{\nu}(z)$, with order (parameter) $\nu$ and argument $z$. EXAMPLES:: Symbolic input and symmetry transformations:: sage: var('nu z') (nu, z) sage: bessel_j(nu, z) bessel_j(nu, z) sage: bessel_j(0, 1) bessel_j(0, 1) sage: bessel_j(nu,0) bessel_j(nu, 0) sage: bessel_j(0,z) bessel_j(0, z) sage: bessel_j(-3, 4) -bessel_j(3, 4) sage: bessel_j(3, -4) -bessel_j(3, 4) sage: bessel_j(-3, -4) bessel_j(3, 4) sage: bessel_j(-4, 5) bessel_j(4, 5) Exact values and limits:: sage: bessel_j(0, 0) 1 sage: bessel_j(1, 0) 0 sage: bessel_j(2, 0) 0 sage: bessel_j(-2, 0) 0 sage: bessel_j(nu, Infinity) 0 sage: bessel_j(nu, -Infinity) 0 sage: bessel_j(1/2,pi) 0 sage: bessel_j(1/2,-8*pi) 0 Numerical evaluation:: sage: bessel_j(0,1).n() 0.765197686557966 sage: bessel_j(0,1).n(digits=30) 0.765197686557966551449717526103 sage: bessel_j(-5, 125000000).n() 0.0000711791694771539 sage: bessel_j(2+3*I,4-5*I).n(digits=30) 1643.61099442312390496727290972 - 848.086110331983438476479306274*I Symbolic and explicit series expansions:: sage: bessel_j(1,z).series(z,5) 1/2*z + (-1/16)*z^3 + Order(z^5) sage: from sage.functions.hypergeometric import expand_hyper sage: expand_hyper(bessel_j(nu,z)) (1/2*z)^nu*PFQ(0, 1, nu + 1, -1/4*z^2)/gamma(nu + 1) sage: expand_hyper(bessel_j(0,z)).series(z,5) 1 + (-1/4)*z^2 + 1/64*z^4 + Order(z^5) sage: expand_hyper(bessel_j(1,z)).series(z,7) 1/2*z + (-1/16)*z^3 + 1/384*z^5 + Order(z^7) """ char = 'J' def _eval_symbolic_(self, nu, z): """ Evaluate or simplify for symbolic input. EXAMPLES :: sage: bessel_j(1, 0) 0 """ try: n = nu._integer_() if n < 0: return (-1)**abs(n) * self(-n,z) if z < -z: return (-1)**abs(n) * self(n,-z) except TypeError: pass if z == 0: r = real_part(nu) if r > 0: return z if r < 0: return Infinity if nu == 0: return 1 if z == Infinity or z == -Infinity: return 0 if 2*nu == 1: if (z / pi)._is_integer(): return 0 return None def _expand_hyper(self, nu, z): """ Expand as a hypergeometric series. sage: from sage.functions.hypergeometric import expand_hyper sage: var('z') z sage: expand_hyper(bessel_j(1,z)) 1/2*z*PFQ(0, 1, 2, -1/4*z^2) """ return (z/2)**nu / gamma(nu+1) * hyp([],[nu+1],-(z/2)**2) class Function_bessel_y(BesselFunction): """ Gives the Bessel function of the second kind $Y_{\nu}(z)$, with order (parameter) $\nu$ and argument $z$. """ char = 'Y' def _eval_symbolic_(self, nu, z): """ Evaluate or simplify for symbolic input. EXAMPLES :: sage: bessel_y(1,1) bessel_y(1, 1) """ return None def _expand_hyper(self, nu, z): """ Expand as a hypergeometric series (this only makes sense when nu is not an integer). sage: from sage.functions.hypergeometric import expand_hyper sage: var('nu z') (nu, z) sage: expand_hyper(bessel_y(nu,z)) (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) sage: expand_hyper(bessel_y(1/2,z)) -2*sqrt(1/2)*cosh(sqrt(-z^2))/(sqrt(pi)*sqrt(z)) sage: expand_hyper(bessel_y(1,z)) Traceback (most recent call last): ... ZeroDivisionError: Symbolic division by zero """ f = (bessel_j(nu,z)*cos(nu*pi) - bessel_j(-nu,z)) / sin(nu*pi) return expand_hyper(f) class Function_bessel_i(BesselFunction): r""" Gives the modified Bessel function of the first kind $I_{\nu}(z)$, with order (parameter) $\nu$ and argument $z$. EXAMPLES:: Symbolic input and symmetry transformations:: sage: var('nu z') (nu, z) sage: bessel_i(nu, z) bessel_i(nu, z) sage: bessel_i(0,1) bessel_i(0, 1) sage: bessel_i(nu,0) bessel_i(nu, 0) sage: bessel_i(0,z) bessel_i(0, z) sage: bessel_i(-3,4) bessel_i(3, 4) sage: bessel_i(3,-4) -bessel_i(3, 4) sage: bessel_i(-3,-4) -bessel_i(3, 4) sage: bessel_i(-4,5) bessel_i(4, 5) Exact values and limits:: sage: bessel_i(0, 0) 1 sage: bessel_i(1, 0) 0 sage: bessel_i(2, 0) 0 sage: bessel_i(-2, 0) 0 sage: bessel_i(1, Infinity) +Infinity sage: bessel_i(1, -Infinity) -Infinity sage: bessel_i(2, Infinity) +Infinity sage: bessel_i(2, -Infinity) +Infinity sage: bessel_i(1/2,I*pi) 0 sage: bessel_i(1/2,-8*I*pi) 0 Numerical evaluation:: sage: bessel_j(0,1).n() 0.765197686557966 sage: bessel_j(0,1).n(digits=30) 0.765197686557966551449717526103 sage: bessel_j(-5, 125000000).n() 0.0000711791694771539 sage: bessel_i(2+3*I, 4-5*I).n(digits=30) -6.69838203392827123593220944996 + 13.1635256869717070482549506439*I Symbolic and explicit series expansions:: sage: from sage.functions.hypergeometric import expand_hyper sage: expand_hyper(bessel_i(nu,z)) (1/2*z)^nu*PFQ(0, 1, nu + 1, 1/4*z^2)/gamma(nu + 1) sage: expand_hyper(bessel_i(0,z)).series(z,5) 1 + 1/4*z^2 + 1/64*z^4 + Order(z^5) sage: expand_hyper(bessel_i(1,z)).series(z,7) 1/2*z + 1/16*z^3 + 1/384*z^5 + Order(z^7) """ char = 'I' def _eval_symbolic_(self, nu, z): """ Evaluate or simplify for symbolic input. EXAMPLES :: sage: bessel_i(0,0) 1 """ try: n = nu._integer_() if n < 0: return self(-n,z) if z < -z: return (-1)**abs(n) * self(n,-z) except TypeError: pass if z == 0: r = real_part(nu) if r > 0: return z if r < 0: return Infinity if nu == 0: return 1 if z == Infinity: return Infinity if z == -Infinity: return (-1)**abs(nu) * Infinity if 2*nu == 1: if (z / (pi*I))._is_integer(): return 0 return None def _expand_hyper(self, nu, z): """ Expand as a hypergeometric series. sage: from sage.functions.hypergeometric import expand_hyper sage: var('z') z sage: expand_hyper(bessel_i(1,z)) 1/2*z*PFQ(0, 1, 2, 1/4*z^2) """ return (z/2)**nu / gamma(nu+1) * hyp([],[nu+1],(z/2)**2) class Function_bessel_k(BesselFunction): r""" Gives the modified Bessel function of the second kind $K_{\nu}(z)$, with order (parameter) $\nu$ and argument $z$. """ char = 'K' def _eval_symbolic_(self, nu, z): """ Evaluate or simplify for symbolic input. EXAMPLES :: sage: bessel_k(1,1) bessel_k(1, 1) """ return None bessel_j = Function_bessel_j() bessel_y = Function_bessel_y() bessel_i = Function_bessel_i() bessel_k = Function_bessel_k()