Take the 2-minute tour ×
Computational Science Stack Exchange is a question and answer site for scientists using computers to solve scientific problems. It's 100% free, no registration required.

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?

share|improve this question

1 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.

share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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