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