/
from sympy.core import Add, S, C, sympify, oo, pi
from sympy.core.function import Function, ArgumentIndexError
from zeta_functions import zeta
from error_functions import erf
from sympy.core import Dummy,Rational
from sympy.functions.elementary.exponential import log
from sympy.functions.elementary.integers import floor
from sympy.functions.elementary.miscellaneous import sqrt
from sympy.functions.combinatorial.numbers import bernoulli
from sympy.functions.combinatorial.factorials import rf
###############################################################################
############################ COMPLETE GAMMA FUNCTION ##########################
###############################################################################
[docs]class gamma(Function):
"""The gamma function returns a function which passes through the integral
values of the factorial function, i.e. though defined in the complex plane,
when n is an integer, `gamma(n) = (n - 1)!`
Reference:
http://en.wikipedia.org/wiki/Gamma_function
"""
nargs = 1
unbranched = True
def fdiff(self, argindex=1):
if argindex == 1:
return gamma(self.args[0])*polygamma(0, self.args[0])
else:
raise ArgumentIndexError(self, argindex)
@classmethod
def eval(cls, arg):
if arg.is_Number:
if arg is S.NaN:
return S.NaN
elif arg is S.Infinity:
return S.Infinity
elif arg.is_Integer:
if arg.is_positive:
return C.factorial(arg-1)
else:
return S.ComplexInfinity
elif arg.is_Rational:
if arg.q == 2:
n = abs(arg.p) // arg.q
if arg.is_positive:
k, coeff = n, S.One
else:
n = k = n + 1
if n & 1 == 0:
coeff = S.One
else:
coeff = S.NegativeOne
for i in range(3, 2*k, 2):
coeff *= i
if arg.is_positive:
return coeff*sqrt(S.Pi) / 2**n
else:
return 2**n*sqrt(S.Pi) / coeff
def _eval_expand_func(self, **hints):
arg = self.args[0]
if arg.is_Rational:
if abs(arg.p) > arg.q:
x = Dummy('x')
n = arg.p // arg.q
p = arg.p - n*arg.q
return gamma(x + n)._eval_expand_func().subs(x, Rational(p, arg.q))
if arg.is_Add:
coeff, tail = arg.as_coeff_add()
if coeff and coeff.q != 1:
intpart = floor(coeff)
tail = (coeff - intpart,) + tail
coeff = intpart
tail = arg._new_rawargs(*tail, **dict(reeval=False))
return gamma(tail)*C.RisingFactorial(tail, coeff)
return self.func(*self.args)
def _eval_is_real(self):
return self.args[0].is_real
def _eval_rewrite_as_tractable(self, z):
return C.exp(loggamma(z))
def _eval_nseries(self, x, n, logx):
x0 = self.args[0].limit(x, 0)
if not (x0.is_Integer and x0 <= 0):
return super(gamma, self)._eval_nseries(x, n, logx)
t = self.args[0] - x0
return (gamma(t + 1)/rf(self.args[0], -x0 + 1))._eval_nseries(x, n, logx)
###############################################################################
################## LOWER and UPPER INCOMPLETE GAMMA FUNCTIONS #################
###############################################################################
[docs]class lowergamma(Function):
r"""
The lower incomplete gamma function.
It can be defined as the meromorphic continuation of
.. math ::
\gamma(s, x) = \int_0^x t^{s-1} e^{-t} \mathrm{d}t.
This can be shown to be the same as
.. math ::
\gamma(s, x) = \frac{x^s}{s} {}_1F_1\left({s \atop s+1} \middle| -x\right),
where :math:`{}_1F_1` is the (confluent) hypergeometric function.
See Also
========
gamma, uppergamma
sympy.functions.special.hyper.hyper
Examples
========
>>> from sympy import lowergamma, S
>>> from sympy.abc import s, x
>>> lowergamma(s, x)
lowergamma(s, x)
>>> lowergamma(3, x)
-x**2*exp(-x) - 2*x*exp(-x) + 2 - 2*exp(-x)
>>> lowergamma(-S(1)/2, x)
-2*sqrt(pi)*erf(sqrt(x)) - 2*exp(-x)/sqrt(x)
References
==========
- Abramowitz, Milton; Stegun, Irene A., eds. (1965), Chapter 6, Section 5,
Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical
Tables
- http://en.wikipedia.org/wiki/Incomplete_gamma_function
"""
nargs = 2
def fdiff(self, argindex=2):
from sympy import meijerg, unpolarify
if argindex == 2:
a, z = self.args
return C.exp(-unpolarify(z))*z**(a-1)
elif argindex == 1:
a, z = self.args
return gamma(a)*digamma(a) - log(z)*uppergamma(a, z) \
+ meijerg([], [1, 1], [0, 0, a], [], z)
else:
raise ArgumentIndexError(self, argindex)
@classmethod
def eval(cls, a, x):
# For lack of a better place, we use this one to extract branching
# information. The following can be
# found in the literature (c/f references given above), albeit scattered:
# 1) For fixed x != 0, lowergamma(s, x) is an entire function of s
# 2) For fixed positive integers s, lowergamma(s, x) is an entire
# function of x.
# 3) For fixed non-positive integers s,
# lowergamma(s, exp(I*2*pi*n)*x) =
# 2*pi*I*n*(-1)**(-s)/factorial(-s) + lowergamma(s, x)
# (this follows from lowergamma(s, x).diff(x) = x**(s-1)*exp(-x)).
# 4) For fixed non-integral s,
# lowergamma(s, x) = x**s*gamma(s)*lowergamma_unbranched(s, x),
# where lowergamma_unbranched(s, x) is an entire function (in fact
# of both s and x), i.e.
# lowergamma(s, exp(2*I*pi*n)*x) = exp(2*pi*I*n*a)*lowergamma(a, x)
from sympy import unpolarify, I, factorial, exp
nx, n = x.extract_branch_factor()
if a.is_integer and a > 0:
nx = unpolarify(x)
if nx != x:
return lowergamma(a, nx)
elif a.is_integer and a <= 0:
if n != 0:
return 2*pi*I*n*(-1)**(-a)/factorial(-a) + lowergamma(a, nx)
elif n != 0:
return exp(2*pi*I*n*a)*lowergamma(a, nx)
# Special values.
if a.is_Number:
# TODO this should be non-recursive
if a is S.One:
return S.One - C.exp(-x)
elif a is S.Half:
return sqrt(pi)*erf(sqrt(x))
elif a.is_Integer or (2*a).is_Integer:
b = a - 1
if b.is_positive:
return b*cls(b, x) - x**b * C.exp(-x)
if not a.is_Integer:
return (cls(a + 1, x) + x**a * C.exp(-x))/a
def _eval_evalf(self, prec):
from sympy.mpmath import mp
from sympy import Expr
a = self.args[0]._to_mpmath(prec)
z = self.args[1]._to_mpmath(prec)
oprec = mp.prec
mp.prec = prec
res = mp.gammainc(a, 0, z)
mp.prec = oprec
return Expr._from_mpmath(res, prec)
def _eval_rewrite_as_uppergamma(self, s, x):
return gamma(s) - uppergamma(s, x)
def _eval_rewrite_as_expint(self, s, x):
from sympy import expint
if s.is_integer and s.is_nonpositive:
return self
return self.rewrite(uppergamma).rewrite(expint)
[docs]class uppergamma(Function):
r"""
Upper incomplete gamma function
It can be defined as the meromorphic continuation of
.. math ::
\Gamma(s, x) = \int_x^\infty t^{s-1} e^{-t} \mathrm{d}t
= \Gamma(s) - \gamma(s, x).
This can be shown to be the same as
.. math ::
\Gamma(s, x) = \Gamma(s)
- \frac{x^s}{s} {}_1F_1\left({s \atop s+1} \middle| -x\right),
where :math:`{}_1F_1` is the (confluent) hypergeometric function.
The upper incomplete gamma function is also essentially equivalent to the
generalized exponential integral.
Examples
========
>>> from sympy import uppergamma, S
>>> from sympy.abc import s, x
>>> uppergamma(s, x)
uppergamma(s, x)
>>> uppergamma(3, x)
x**2*exp(-x) + 2*x*exp(-x) + 2*exp(-x)
>>> uppergamma(-S(1)/2, x)
-2*sqrt(pi)*(-erf(sqrt(x)) + 1) + 2*exp(-x)/sqrt(x)
>>> uppergamma(-2, x)
expint(3, x)/x**2
See Also
========
gamma, lowergamma
sympy.functions.special.hyper.hyper
References
==========
- Abramowitz, Milton; Stegun, Irene A., eds. (1965), Chapter 6, Section 5,
Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical
Tables
- http://en.wikipedia.org/wiki/Incomplete_gamma_function
"""
nargs = 2
def fdiff(self, argindex=2):
from sympy import meijerg, unpolarify
if argindex == 2:
a, z = self.args
return -C.exp(-unpolarify(z))*z**(a-1)
elif argindex == 1:
a, z = self.args
return uppergamma(a, z)*log(z) + meijerg([], [1, 1], [0, 0, a], [], z)
else:
raise ArgumentIndexError(self, argindex)
def _eval_evalf(self, prec):
from sympy.mpmath import mp
from sympy import Expr
a = self.args[0]._to_mpmath(prec)
z = self.args[1]._to_mpmath(prec)
oprec = mp.prec
mp.prec = prec
res = mp.gammainc(a, z, mp.inf)
mp.prec = oprec
return Expr._from_mpmath(res, prec)
@classmethod
def eval(cls, a, z):
from sympy import unpolarify, I, factorial, exp, expint
if z.is_Number:
if z is S.NaN:
return S.NaN
elif z is S.Infinity:
return S.Zero
elif z is S.Zero:
return gamma(a)
# We extract branching information here. C/f lowergamma.
nx, n = z.extract_branch_factor()
if a.is_integer and a > 0:
nx = unpolarify(z)
if z != nx:
return uppergamma(a, nx)
elif a.is_integer and a <= 0:
if n != 0:
return -2*pi*I*n*(-1)**(-a)/factorial(-a) + uppergamma(a, nx)
elif n != 0:
return gamma(a)*(1 - exp(2*pi*I*n*a)) + exp(2*pi*I*n*a)*uppergamma(a, nx)
# Special values.
if a.is_Number:
# TODO this should be non-recursive
if a is S.One:
return C.exp(-z)
elif a is S.Half:
return sqrt(pi)*(1 - erf(sqrt(z))) # TODO could use erfc...
elif a.is_Integer or (2*a).is_Integer:
b = a - 1
if b.is_positive:
return b*cls(b, z) + z**b * C.exp(-z)
elif b.is_Integer:
return expint(-b, z)*unpolarify(z)**(b + 1)
if not a.is_Integer:
return (cls(a + 1, z) - z**a * C.exp(-z))/a
def _eval_rewrite_as_lowergamma(self, s, x):
return gamma(s) - lowergamma(s, x)
def _eval_rewrite_as_expint(self, s, x):
from sympy import expint
return expint(1 - s, x)*x**s
###############################################################################
########################### GAMMA RELATED FUNCTIONS ###########################
###############################################################################
[docs]class polygamma(Function):
"""The function `polygamma(n, z)` returns `log(gamma(z)).diff(n + 1)`
See Also
========
gamma, digamma, trigamma
Reference:
http://en.wikipedia.org/wiki/Polygamma_function
"""
nargs = 2
def fdiff(self, argindex=2):
if argindex == 2:
n, z = self.args[:2]
return polygamma(n+1, z)
else:
raise ArgumentIndexError(self, argindex)
def _eval_is_positive(self):
if self.args[1].is_positive and self.args[0] > 0:
return self.args[0].is_odd
def _eval_is_negative(self):
if self.args[1].is_positive and self.args[0] > 0:
return self.args[0].is_even
def _eval_is_real(self):
return self.args[0].is_real
def _eval_aseries(self, n, args0, x, logx):
if args0[1] != oo or not \
(self.args[0].is_Integer and self.args[0].is_nonnegative):
return super(polygamma, self)._eval_aseries(n, args0, x, logx)
z = self.args[1]
N = self.args[0]
if N == 0:
# digamma function series
# Abramowitz & Stegun, p. 259, 6.3.18
r = log(z) - 1/(2*z)
o = None
if n < 2:
o = C.Order(1/z, x)
else:
m = C.ceiling((n+1)//2)
l = [bernoulli(2*k) / (2*k*z**(2*k)) for k in range(1, m)]
r -= Add(*l)
o = C.Order(1/z**(2*m), x)
return r._eval_nseries(x, n, logx) + o
else:
# proper polygamma function
# Abramowitz & Stegun, p. 260, 6.4.10
# We return terms to order higher than O(x**n) on purpose
# -- otherwise we would not be able to return any terms for
# quite a long time!
fac = gamma(N)
e0 = fac + N*fac/(2*z)
m = C.ceiling((n+1)//2)
for k in range(1, m):
fac = fac*(2*k+N-1)*(2*k+N-2) / ((2*k)*(2*k-1))
e0 += bernoulli(2*k)*fac/z**(2*k)
o = C.Order(1/z**(2*m), x)
if n == 0:
o = C.Order(1/z, x)
elif n == 1:
o = C.Order(1/z**2, x)
r = e0._eval_nseries(z, n, logx) + o
return -1 * (-1/z)**N * r
@classmethod
def eval(cls, n, z):
n, z = map(sympify, (n, z))
from sympy import unpolarify
if n.is_integer:
if n.is_nonnegative:
nz = unpolarify(z)
if z != nz:
return polygamma(n, nz)
if n == -1:
return loggamma(z)
else:
if z.is_Number:
if z is S.NaN:
return S.NaN
elif z is S.Infinity:
if n.is_Number:
if n is S.Zero:
return S.Infinity
else:
return S.Zero
elif z.is_Integer:
if z.is_nonpositive:
return S.ComplexInfinity
else:
if n is S.Zero:
return -S.EulerGamma + C.harmonic(z-1, 1)
elif n.is_odd:
return (-1)**(n+1)*C.factorial(n)*zeta(n+1, z)
if n == 0 and z.is_Rational:
# TODO actually *any* n/m can be done, but that is messy
lookup = {S(1)/2: -2*log(2) - S.EulerGamma,
S(1)/3: -S.Pi/2/sqrt(3) - 3*log(3)/2 - S.EulerGamma,
S(1)/4: -S.Pi/2 - 3*log(2) - S.EulerGamma,
S(3)/4: -3*log(2) - S.EulerGamma + S.Pi/2,
S(2)/3: -3*log(3)/2 + S.Pi/2/sqrt(3) - S.EulerGamma}
if z > 0:
n = floor(z)
z0 = z - n
if z0 in lookup:
return lookup[z0] + Add(*[1/(z0 + k) for k in range(n)])
elif z < 0:
n = floor(1 - z)
z0 = z + n
if z0 in lookup:
return lookup[z0] - Add(*[1/(z0 - 1 - k) for k in range(n)])
# TODO n == 1 also can do some rational z
def _eval_expand_func(self, **hints):
n, z = self.args
if n.is_Integer and n.is_nonnegative:
if z.is_Add:
coeff = z.args[0]
if coeff.is_Integer:
e = -(n + 1)
if coeff > 0:
tail = Add(*[C.Pow(z - i, e) for i in xrange(1, int(coeff) + 1)])
else:
tail = -Add(*[C.Pow(z + i, e) for i in xrange(0, int(-coeff))])
return polygamma(n, z - coeff) + (-1)**n*C.factorial(n)*tail
elif z.is_Mul:
coeff, z = z.as_two_terms()
if coeff.is_Integer and coeff.is_positive:
tail = [ polygamma(n, z + C.Rational(i, coeff)) for i in xrange(0, int(coeff)) ]
if n == 0:
return Add(*tail)/coeff + log(coeff)
else:
return Add(*tail)/coeff**(n+1)
z *= coeff
return polygamma(n, z)
def _eval_rewrite_as_zeta(self, n, z):
return (-1)**(n+1)*C.factorial(n)*zeta(n+1, z-1)
[docs]class loggamma(Function):
"""
The loggamma function is `ln(gamma(x))`.
References
==========
http://mathworld.wolfram.com/LogGammaFunction.html
"""
nargs = 1
def _eval_aseries(self, n, args0, x, logx):
if args0[0] != oo:
return super(loggamma, self)._eval_aseries(n, args0, x, logx)
z = self.args[0]
m = min(n, C.ceiling((n+S(1))/2))
r = log(z)*(z-S(1)/2) - z + log(2*pi)/2
l = [bernoulli(2*k) / (2*k*(2*k-1)*z**(2*k-1)) for k in range(1, m)]
o = None
if m == 0:
o = C.Order(1, x)
else:
o = C.Order(1/z**(2*m-1), x)
# It is very inefficient to first add the order and then do the nseries
return (r + Add(*l))._eval_nseries(x, n, logx) + o
def _eval_rewrite_as_intractable(self, z):
return log(gamma(z))
def _eval_is_real(self):
return self.args[0].is_real
def fdiff(self, argindex=1):
if argindex == 1:
return polygamma(0, self.args[0])
else:
raise ArgumentIndexError(self, argindex)
[docs]def digamma(x):
"""
The digamma function is the logarithmic derivative of the gamma function.
In this case, `digamma(x) = polygamma(0, x)`.
See Also
========
gamma, trigamma, polygamma
"""
return polygamma(0, x)
[docs]def trigamma(x):
"""
The trigamma function is the second of the polygamma functions.
In this case, `trigamma(x) = polygamma(1, x)`.
See Also
========
gamma, digamma, polygamma
"""
return polygamma(1, x)
[docs]def beta(x, y):
"""
Euler Beta function
``beta(x, y) == gamma(x)*gamma(y) / gamma(x+y)``
"""
return gamma(x)*gamma(y) / gamma(x+y)