Using Log Gamma function to avoid overflow

Question

I have to do some numerical calculus using gamma functions. I am using the tgamma incluided in the C++ cmath library. The problem is that the arguments inside the gamma functions are positive and big, so I always get a nan.

I know that I can use a log-gamma function, lgamma in C++, to avoid numerical overflow during the calculation. The problem is that I don't know how to apply this in my case.

For example, with something like

$$X=\dfrac{\Gamma(a+b)\Gamma(a-b)}{\Gamma(c)}$$

where $a$, $b$ and $c$ are integers largue enough to produce an overflow.

I know that I could use a log to obtain $\log(X)=\log(\Gamma(a+b)) + \log(\Gamma(a-b)) -\log(\Gamma(c))$, obtain the result from that operation and make the exponential to get the correct result without overflow.

But, what about if I have something like this?

$$X=\sum_i^N p_i\cdot\dfrac{\Gamma(i+b)\Gamma(i-b)}{\Gamma(c)} + \sum_i^N\sum_{j<i} q_i\cdot \dfrac{\Gamma(i+b)\Gamma(j-b)}{\Gamma(c)}$$

In this case $p_i$ and $q_i$ are real numbers.

Now I have a sum, so the $\log$ properties are useless. Of course, I could simply substitute all the gammas with log-gammas, but, in that case, how could I retrieve my original result? Can I do something to avoid the overflow?

Answer 1

You can avoid overflow by rewriting the sum appropriately. Say you have a sum $$ X = \sum_{i=1}^{n} x_i, $$ where some $x_i$'s are positive and large enough to cause overflow, but $\log x_i$ are of reasonable magnitude. If you find the greatest $\log x_M \geq \log x_i$, then you can rewrite the sum as $$ \log X = \log \sum_i e^{\log x_i} = \log x_M + \log \sum_i e^{\log x_i - \log x_M}. $$ By doing this, only logs are used in the computation. In the remaining sum inside the log we have $\log x_i \leq \log x_M$, so $e^{\log x_i - \log x_M} \leq 1$, which means the value of the remaining sum is at most $n$.

Obviously this doesn't work if $x_i$'s are allowed to sometimes be negative - in that case you might get catastrophic cancellation if they can be negative and relatively large compared to the final value, although there is no problem if they are negative and sufficiently small.