Bill Press
11-28-2008, 01:17 PM
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)
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)