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


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.

| cite | improve this answer | |
  • $\begingroup$ Many thanks. I was unable to see the simple solution. It is, however, unclear why such a substitution is used. Perhaps using LGAMMA function is numerically better than using GAMMA function itself. Again, many thanks. $\endgroup$ – rk85 Sep 2 '17 at 13:56
  • $\begingroup$ Just to add. The reference to Applied Statistics does not contain any description, but the original old Fotran code: lib.stat.cmu.edu/apstat/147. $\endgroup$ – rk85 Sep 2 '17 at 14:08

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.