Computing Poisson probabilities

We propose an algorithm to compute the set of individual (nonnegligible) Poisson probabilities, rigorously bound truncation error, and guarantee no overflow or underflow. Work and space requirements are modest, both proportional to the square root of the Poisson parameter. Our algorithm appears numerically stable. We know no other algorithm with all these (good) features. Our algorithm speeds generation of truncated Poisson variates and the computation of expected terminal reward in continuous-time, uniformizable Markov chains. More generally, our algorithm can be used to evaluate formulas involving Poisson probabilities.