incgammabeta.h: bug in Gamma (fixed in V3.04)

[Bill Press]

In versions 3.02 and earlier, there is a bug in Gamma. For (unusually) large values of a and x, and with x << a, so that gammp should be 0, gammq should be 1, incorrect values (that swap 1 and 0) can be returned.

The bug is in gammpapprox(). To fix the bug, the line

return (psig?(ans>0.0? 1.0-ans:-ans):(ans>=0.0? ans:1.0+ans));

should be replaced by

return (psig?(x>a1? 1.0-ans:-ans):(x>a1? ans:1.0+ans));

This bug is also responsible for Poissondist::cdf(n) throwing exceptions in unusual cases, see here

[emilpohl]

Hi!

I'm using NR3 hosted by the International Centre for Theoretical Physics. I've noticed that the code:

return (psig?(ans>0.0? 1.0-ans:-ans):(ans>=0.0? ans:1.0+ans));

wasn't changed after the Bug report. Is there a reason for this?

best regards,

Carlos

[Bill Press]

Good call, Carlos. In general, we only change the code distributions when we change the version numbers (often corresponding to new printings of the book) when many bugs are fixed at one time. We don't change the distributions after each individual bug.
Cheers,
Bill P.