# Calculation among Gamma functions

I was calculating gamma functions in Mathematica while it does not give me an agreed answer.

By definition, $$\Gamma[\alpha]=\int_0^\infty t^{\alpha-1}e^{-t}dt$$, $$\Gamma[\alpha,z]=\int_z^\infty t^{\alpha-1}e^{-t}dt$$ and $$\Gamma[\alpha,z_1,z_2]=\int_{z_1}^{z_2}t^{\alpha-1}e^{-t}dt$$. So we should have $$\Gamma[\alpha]-\Gamma[\alpha,z]$$ should be equal to $$\Gamma[\alpha,0,z]$$. But this is not the case. I simply let $$\alpha=200$$ and $$z$$ from 1 to 5, he numerical values do not fit; see below.

 α = 200.;
Table[{Gamma[α] - Gamma[α, z], Gamma[α, 0, z]}, {z, 1, 5}]

{{0.*10^359, 0.00184859}, {0.*10^359, 1.0983*10^57}, {0.*10^359,6.71225*10^91},
{0.*10^359, 2.41279*10^116}, {0.*10^359,  2.14999*10^135}}


Why is that?

## 1 Answer

This is precisely the reason Mathematica supplies the three-argument form of Gamma[] (and GammaRegularized[] as well); in computing Gamma[α] - Gamma[α, z] for large α, you are effectively subtracting two large numbers to get a (comparatively) tiny result. (Similar remarks hold for computing Cosh[x] - Sinh[x] just to evaluate Exp[-x]; attempts to perform this for sufficiently large x are doomed.)

Thus, whenever one has to compute differences of incomplete gamma functions (as in this recent example), think about whether it can be reformulated in terms of three-argument Gamma[].

See this related thread on Erf[] as well.