numerical evaluation of quotients of gamma functions
Subject: numerical evaluation of quotients of gamma functions
From: Barton Willis
Date: Fri, 1 Nov 2013 17:39:54 +0000
I need to numerically evaluate quotients of products of gamma functions. In order of preference, I would like to
(1) avoid overflows, (2) have good accuracy, and (3) be fast (much less important). Any suggestions?
Simplifications such as gamma(x + 5) / gamma(x) --> pochhmmer(x,5) are pretty safe, especially if x is rational. But
for gamma(1/3 + 100) / gamma(1/3), rational arithmetic might be too spendy.
(%i42) gamma(1/3 + 100) / gamma(1/3);
(%o42) 832434801314587533188150459656[144 digits]465280000000000000000000000000/515377520732011331036461129765621272702107522001
(%i43) float(%);
(%o43) 1.6151942369008391*10^156
(%i44) gamma(1.0/3 + 100) / gamma(1.0/3);
(%o44) 1.6151942369008382*10^156
Also for what I need it would be OK to only consider quotients of the form gamma(x) * gamma(y) / (gamma(p) * gamma(q)).
There is the fun paper
http://projecteuclid.org/DPubS/Repository/1.0/Disseminate?view=body&id=pdf_1&handle=euclid.pjm/1102613160
But that requires some (possibly) messy logic in deciding when to switch to the asymptotic formula--I was hoping for something
simple.
--Barton