numerical evaluation of quotients of gamma functions



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