IMO, the most significant improvement in the latest version of the code is that it computes at most 18 factorials altogether (9 in the onesided case) per call to fishers_exact. After the first term in the sum, all subsequent terms are computed as multiplicative adjustments to the first term:
$delta *= ( ( $aa * $dd )/( ++$bb * ++$cc ) );
I have not seen any significant loss of accuracy from this in my tests.
To put in perspective the effect of using Gosper's approximations only when Math::Pari::factorial bombs, instead of using a hard cutoff of 10000 for the switchover, the relative difference (i.e. difference divided by average) between M::P::factorial(10_000) and Gosper's approximation to 10,000! is smaller than 1e10, and it gets only better for larger arguments. Therefore Gosper's is more than adequate in this range, AFAIC.
My main reason for removing the hard cutoff was just to simplify the code for posting. In my own production code I don't wait for M::P::factorial to fail before switching over to Gosper's, because for large arguments the benefits in terms of accuracy from using M::P::factorial are minuscule but the performance hit is substantial. For example, for 10,000! (relative error for Gosper's: < 1e10):
Rate pari gospers
pari 604/s  94%
gospers 9657/s 1500% 
and for 100,000! (relative error for Gosper's: < 1e12):
Rate pari gospers
pari 49.8/s  100%
gospers 10001/s 19996% 
Also, note that, as for the previous version of the code, this one returns a value very close to 1 for fishers_exact( 989, 9400, 43300, 2400 ).
And thanks for the headsup about the stray code in Fisher's Exact Test. Fixed.
