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 (http://www.nr.com/forum/showthread.php?t=1658)

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

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.