gammainc

The incomplete gamma function

Definition

There are actually four types of incomplete gamma functions:

• lower (which is the default one) $$P(a,x) = \frac{1}{\Gamma(a)}\int_{0}^{x}e^{-t}t^{a-1}$$
• upper $$Q(a,x) = \frac{1}{\Gamma(a)}\int_{x}^{+\infty}t^{a-1}e^{-t}dt$$
• scaled lower and scaled upper which return the corresponding incomplete gamma function scaled by $$\Gamma(a+1)\frac{e^x}{x^a}$$

Computation

To compute the incomplete gamma function, there is a total of seven strategies depending on the position of the arguments in the $x-a$ plane.

1. $x=0,a=0$. In this case, the output is given as follows:
   • for "lower" and "scaledlower" is $1$
   • for "upper" and "scaledupper" is $0$
2. $x=0,a\neq0$. In this case:
   • for "lower" is $0$
   • for "upper" and "scalelower" is $1$
   • for "scaledupper" is $+\infty$
3. $x\neq0,a=0$. In this case:
   • for "lower" is $1$
   • for "scalelower" is $e^x$
   • for "upper" and "scaledupper" is $0$
4. $x\neq0,a=1$. In this case:
   • for "lower" is $1-e^x$
   • for "upper" is $e^{-x}$
   • for "scaledlower" is $\frac{e^x-1}{x}$
   • for "scaledupper" is $\frac{1}{x}$
5. $0< x\leq36,a\in\{2,3,4,\ldots,17,18\}$. In this case we used the following expansion (which holds true only for $n\in\mathbb{N}$): denote $\gamma(a,x):=\Gamma(a)P(a,x)$, then $$\gamma(n,x) = (n-1)!\left( 1-e^{-x}\sum_{k=0}^{n-1} \frac{x^k}{k!}\right).$$
6. $x+0.25< a $ or $x<0$, or $x$ not real or $|x|<1$ or $a<5$. Here we used the expansion $$\gamma(a,x) = e^{-x} x^a \sum_{k=0}^{+\infty} \frac{\Gamma(a)}{\Gamma(a+1+k)}x^k.$$
7. for all the remaining cases was used a continued fraction expansion. Call $\Gamma(a,x) = \Gamma(a) Q(a,x)$, the expansion is $$\Gamma(a,x) = e^{-x}x^a\left( \frac{1}{x+1-a-}\frac{1(1-a)}{x+3-a-}\frac{2(2-a)}{x+5-a-}\cdots \right).$$ This case is handled by an external .cc function.

My contribution

Actually I dind't participate from the start in the (re)making of this function: most of the work done on it is due to Nir and Marco (see the discussion on the bug tracker).
My contribution was, before the GSoC, to adapt the codes in such a way to make the type of the output (and the tolerances inside the algorithms) coherent on the type of the input (single, int or double). Then I corrected few small bugs and made a patch helped by Carne during the OctConf in Geneve.
During this first week of GSoC I fixed the input validation and added tests to this function, finding some problem in the implementation. Most of them depended on the fact that the continued fractions were used in cases in which they are not optimal. To solve these problems I changed the conditions to use the series expansion instead of the continued fractions.
Now gammainc works properly also with complex argument (for the $x$ variable), while Matlab doesn't accept non-real input.
You can find my work on my repository on bitbucket here, bookmark "gammainc". I will work on this same bookmark the next two weeks while I will work on gammaincinv.