Incomplete gamma function algorithm

Question

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

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

We coded this in C++. There is one thing I don't understand in this algorithm. In one place to compute the following expression:

it is substituted by:

Obviously this is the same, but why it is done like this? Is computing exp of lgamma more efficient than computing tgamma function (both lgamma and tgamma are available in C++)?

Answer 1

Does computing exp of lgamma is more efficient than computing tgamma function (both lgamma and tgamma are available in C++)?

computing lgamma is more efficient because it is basically n*log(n) behavior. So typically you have good approximation id you're trying to compute lgamma(x)/x function.

also, keep in mind, that lgamma is often used because it is part of expression which could be computed while tgamma, which is growing as fast as factorial, is not. So it is safe to compute whole expression f(x) as exp(log(f(x))), and if f(x) has product of tgamma, then log(f(x)) will have to sum/substract lgamma.

Good way to avoid overflow, basically

Answer 2

You will find proper implementations of the Gamma for c++ here: http://www.boost.org/doc/libs/1_64_0/libs/math/doc/html/math_toolkit/sf_gamma