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?