Generating random numbers with a Poisson distribution (how-to)

The naive cumulative probability method, though the easiest method to think of, is not necessarily the simplest or fastest method, and the cumulative version of the coefficients has numerical stability issues related to repeated floating-point multiplication and division, which can lead to a distorted distribution. In fact, if generating random numbers is faster than floating point multiplication and division, the following method, developed by Knuth, will be faster:

const int PoissonRandomNumber(const double lambda)
{
  int k=0;
  const double target=exp(-lambda);
  double p=UniformRandomNumber();
  while (p>target)
  {
    p*=UniformRandomNumber();
    k+=1;
  }
  return k;
}

Although the code's operation is easy to understand, that it actually produces a Poisson distribution is somewhat harder. One explanation is to note that, if u is a uniformly-distributed random number, -ln(u) has an exponential distribution, which is the distribution that describes the time between events in a Poisson process; the Poisson distribution describes the number of events that occur in a specific amount of time.