# Incomplete gamma function algorithm

There is a very concise algorithm for computing incomplete gamma function:

https://people.sc.fsu.edu/~jburkardt/f_src/asa147/asa147.html

However, when I look extensively over the Internet or trying to find any references I do not find a rationale for this algorithm. For example why we use log gamma function into it? Would be grateful for any comments which explain the algorithm construction and the mathematical basis for it.

There isn't any magic in it, and the source code does contain references: one to a standard book on special functions (a rough Internet equivalent would be this), another one to an article in "Applied Statistics", but the algorithm uses a simple and well-known series expansion (essentially the one here). It computes $$\frac{\gamma(p,x)}{\Gamma(p)}=x^p e^{-x}\sum^\infty_{k=0}\frac{x^k}{\Gamma(k+p+1)},$$ starting from the first term $$\frac{x^p e^{-x}}{\Gamma(p+1)}=\exp(p\ln x-\ln\Gamma(p+1)-x),$$ that's where the function LGAMMA is used.