mrpy.base.special.gammainc¶
-
mrpy.base.special.
gammainc
(z, x)¶ Upper incomplete gamma function.
Note
This is exactly as defined by mpmath, but modified to take
array_like
arguments, and return results with type float.Notes
gammainc(z, a=0, b=inf)
computes the (generalized) incomplete gamma function with integration limits [a, b]:\[\Gamma(z,a,b) = \int_a^b t^{z-1} e^{-t} \, dt\]The generalized incomplete gamma function reduces to the following special cases when one or both endpoints are fixed:
- Gamma(z,0,infty) is the standard (“complete”)
gamma function, Gamma(z) (available directly
as the mpmath function
gamma()
) - Gamma(z,a,infty) is the “upper” incomplete gamma function, Gamma(z,a)
- Gamma(z,0,b) is the “lower” incomplete gamma function, gamma(z,b).
Of course, we have Gamma(z,0,x) + Gamma(z,x,infty) = Gamma(z) for all z and x.
Note however that some authors reverse the order of the arguments when defining the lower and upper incomplete gamma function, so one should be careful to get the correct definition.
If also given the keyword argument
regularized=True
,gammainc()
computes the “regularized” incomplete gamma function\[P(z,a,b) = \frac{\Gamma(z,a,b)}{\Gamma(z)}.\]Examples
We can compare with numerical quadrature to verify that
gammainc()
computes the integral in the definition:>>> from mpmath import * >>> mp.dps = 25; mp.pretty = True >>> gammainc(2+3j, 4, 10) (0.00977212668627705160602312 - 0.0770637306312989892451977j) >>> quad(lambda t: t**(2+3j-1) * exp(-t), [4, 10]) (0.00977212668627705160602312 - 0.0770637306312989892451977j)
Argument symmetries follow directly from the integral definition:
>>> gammainc(3, 4, 5) + gammainc(3, 5, 4) 0.0 >>> gammainc(3,0,2) + gammainc(3,2,4); gammainc(3,0,4) 1.523793388892911312363331 1.523793388892911312363331 >>> findroot(lambda z: gammainc(2,z,3), 1) 3.0
Evaluation for arbitrarily large arguments:
>>> gammainc(10, 100) 4.083660630910611272288592e-26 >>> gammainc(10, 10000000000000000) 5.290402449901174752972486e-4342944819032375 >>> gammainc(3+4j, 1000000+1000000j) (-1.257913707524362408877881e-434284 + 2.556691003883483531962095e-434284j)
Evaluation of a generalized incomplete gamma function automatically chooses the representation that gives a more accurate result, depending on which parameter is larger:
>>> gammainc(10000000, 3) - gammainc(10000000, 2) # Bad 0.0 >>> gammainc(10000000, 2, 3) # Good 1.755146243738946045873491e+4771204 >>> gammainc(2, 0, 100000001) - gammainc(2, 0, 100000000) # Bad 0.0 >>> gammainc(2, 100000000, 100000001) # Good 4.078258353474186729184421e-43429441
The incomplete gamma functions satisfy simple recurrence relations:
>>> mp.dps = 25 >>> z, a = mpf(3.5), mpf(2) >>> gammainc(z+1, a); z*gammainc(z,a) + a**z*exp(-a) 10.60130296933533459267329 10.60130296933533459267329 >>> gammainc(z+1,0,a); z*gammainc(z,0,a) - a**z*exp(-a) 1.030425427232114336470932 1.030425427232114336470932
Evaluation at integers and poles:
>>> gammainc(-3, -4, -5) (-0.2214577048967798566234192 + 0.0j) >>> gammainc(-3, 0, 5) +inf
If z is an integer, the recurrence reduces the incomplete gamma function to P(a) exp(-a) + Q(b) exp(-b) where P and Q are polynomials:
>>> gammainc(1, 2); exp(-2) 0.1353352832366126918939995 0.1353352832366126918939995 >>> mp.dps = 50 >>> identify(gammainc(6, 1, 2), ['exp(-1)', 'exp(-2)']) '(326*exp(-1) + (-872)*exp(-2))'
The incomplete gamma functions reduce to functions such as the exponential integral Ei and the error function for special arguments:
>>> mp.dps = 25 >>> gammainc(0, 4); -ei(-4) 0.00377935240984890647887486 0.00377935240984890647887486 >>> gammainc(0.5, 0, 2); sqrt(pi)*erf(sqrt(2)) 1.691806732945198336509541 1.691806732945198336509541
- Gamma(z,0,infty) is the standard (“complete”)
gamma function, Gamma(z) (available directly
as the mpmath function