The work horse for dpois, dgamma, and dchisq is dpois_raw. It has problems with accuracy when the result is small. For example: > options(digits=20) > dpois(1400,1000) [1] 1.466773344196698e-33 http://www.wolframalpha.com/input/?i=1000^1400*Exp[-1000]%2F1400! 1.466773344196568303089528189573057152277244031720195885... × 10^-33 I.e., only about 13 correct digits. The basic problem is that the function works with log-probability. Basically it does something*exp(-stirlerr(x)-bd0(x,lambda)) stirlerr is about 1/(12x) and thus harmless, but bd0 can get big. In this case it's about 71 and the rounding error on its value is blown up by exp. It's not entirely clear to me how to solve this completely, but there are things that can be done that will be faster and more accurate for certain parameter ranges. Assuming the non-log case: x<=170 && lambda<690: pow(lambda,x)*exp(-lambda)/fact(x) This assumes an accurate fact(double) function. An accurate gamma function can be used if also x>=1 (so that doing x+1 doesn't affect accuracy much). Note: pow(lambda,x) can overflow. If that happens, bail out and try the next formula: x<=170 && lambda<2*690: (f < 1 ? f*f/fact(x) : f/fact(x)*f) where f=pow(lambda,x/2)*exp(-lambda/2) Again, pow can overflow and that needs to be detected. The above can be used for the log case too, except when the products underflow. These ranges are too restrictive to help with the dpois(1400,1000) case, but they will improve things like dpois(40,7.5).

Created attachment 1552 [details] Extended bd0 code I think this problem can be solved with the attached code. The code was developed for Gnumeric, but should work more or less directly in R. The code implements ebd0, an extended version of bd0, that returns its result in two parts. The basic approach is to scale x/M to become close to 1 using a factor from a fixed set of factors whose logs are known to higher-than-double precision. I used ~100 bits precision which is probably somewhat higher than needed. Someone ought to check the values. The result is split in whole and fractional part in a simple way. One could use Dekker addition for even higher accuracy. With this, the main case of dpois_raw just becomes ebd0 (x, lambda, &yh, &yl); yl += stirlerr(x); f2 = M_2PI * x; return give_log ? -yl - yh - 0.5 * log (f2) : exp (-yl) * exp (-yh) / sqrt (f2); Note, that dbinom_raw and dt appear to suffer from the same problem.