4

I'm estimating the ratio of two gamma functions. Both estimates for gamma(x) and gamma(y) are very large ( > 10^300) but the ratio of the two should be fairly small. 

from scipy.special import gamma
gamma(x) / gamma(y)

Unfortunately, there is a point where gamma(x) is too large and the scipy returns an inf value. Is there a way to stop this, increase the threshold, or alternatively calculate this ratio?

Thanks

user3439329
  • 851
  • 4
  • 10
  • 24

1 Answers1

4

Assuming x and y are real and positive you can use the following identity:

a / b = e^(ln(a) - ln(b))

I can suggest using gmpy2 for the arbitrary precision calculation. It has gmpy2.lgamma, which returns the logarithm of the gamma function, and you can use gmpy2.exp to convert from the logarithm form to your desired ratio.

orlp
  • 112,504
  • 36
  • 218
  • 315
  • 1
    `scipy.special` has `gammaln` (http://docs.scipy.org/doc/scipy/reference/generated/scipy.special.gammaln.html), so `exp(gammaln(a) - gammaln(b))` should work. – Warren Weckesser Sep 10 '15 at 17:56
  • @WarrenWeckesser No, `scipy.special`'s data type is a `double`, so it can't handle the big numbers the OP requires. – orlp Sep 10 '15 at 18:16
  • Right, that's the point of using `gammaln`. The OP says `gamma(x)` is greater than 1e300, but doesn't say how much greater. Suppose `x` is, say, 425, and `y` is 420. Both `gamma(x)` and `gamma(y)` are much greater than 1e300 (and far exceed the maximum 64 bit floating point value), but `exp(gammaln(x) - gammaln(y))` is easily computed using 64 bit floating point to be `1.3883e13`. How far one can go with 64 bit floating point depends on how big `x` and `y` really are, and on how big the ratio of `gamma(x)/gamma(y)` is. – Warren Weckesser Sep 10 '15 at 19:06
  • @WarrenWeckesser Oh, I assumed `x` was on the order of 10^300. Misread. If your `x` and `y` fit in a double, feel free to use the scipy function. – orlp Sep 10 '15 at 19:09