Simple algorithm for generating Poisson distribution I found a very simple algorithm that draws values from a Poisson distribution from this project.
The algorithm's code in Java is:
public final int poisson(double a) {
        double limit = Math.exp(-a), prod = nextDouble();
        int n;
        for (n = 0; prod >= limit; n++)
            prod *= nextDouble();
        return n;
        }

nextDouble() is a function from the Random package in Java that returns a uniformly distributed random double, for example 0.885598042879084.
I can't understand how this creates a Poisson distribution. 
Can someone explain?
 A: It is related to the Poisson process: suppose $a$ is fixed and $N$ is the number of independent Exponential (mean 1) RV's to be added until the sum exceeds $a$. In this case, $N \sim Poisson(a)$. 
In the code above, everything is anti-logged. For example, instead of adding Exponentials, they multiply Uniforms, due to the relation $Y = - \ln(U)$ is an Exponential (mean 1) RV whenever $U$ is Uniform[0,1].
A: Suppose that we have $k$ independent and identically distributed exponential random variables $X_1, \dotsc, X_k$ with parameter $\lambda$. If we define a counting process $\{N(t)\}_{t \geq 0}$ such that $S_k := X_1 + \dotsb + X_k$ is the occurring time of the $k$-th event, then this is a Poisson process with rate $\lambda$. This means that $\mathbb{P}\{N(t) \leq k\} = \mathbb{P}\{S_k \geq t\}$, or more informally that you can describe the distribution of $N(t)$ through the distribution of $S_k$.
Now consider that if $X$ is an exponential random variable with parameter $\lambda$, then $Y := F_X(X) = 1 - e^{-\lambda X}$ is a standard uniform random variable (here $F_X$ is the cumulative distribution function of $X$). Note also that
$$
F_X^{-1}(Y) = -\frac{1}{\lambda}\log(1 - Y)
$$
is an exponential random variable with parameter $\lambda$. Since also $U := 1 - Y$ is a standard uniform random variable, the random variable
$$
X' = -\frac{1}{\lambda}\log U
$$
is again an exponential random variable with parameter $\lambda$.
We want to count (at most) $k$ events in a given time interval $[0, t]$ (that we can suppose without loss of generality to be $[0, 1]$, i.e. we arbitrarily fix $t = 1$), that is we want $N(1) \leq k$. By what we observed at the beginning, this corresponds to the condition $S_k \geq 1$, which, in its turn, corresponds to the following:
$$
\sum_{i = 1}^k X_i \geq 1 \iff \sum_{i = 1}^k -\frac{1}{\lambda}\log U_i \geq 1 \iff -\frac{1}{\lambda}\log\bigg( \prod_{i = 1}^k U_i \bigg) \geq 1 \iff
\prod_{i = 1}^k U_i \leq e^{-\lambda},
$$
which is precisely what is described in Knuth's algorithm.
A: From wikipedia:

The Poisson distribution is a discrete probability distribution that expresses the probability of a given number of events occurring in a fixed interval of time and/or space if these events occur with a known average rate and independently of the time since the last event.

I think that you could get something similar by adding together random values until you hit a maximum, as the random function has a "known average rate" and is basically "independent" of the last function call. It looks to me like multiplying it like so is just a transformation to simplify the code.
A: Yes, the previous answer is correct. This is simply a version of Knuth's algorithm.  It gets slower the larger the parameter lambda.  
