I managed to find the source code in sexp.c, and the algorithm (Ahrens & Dieter). I mostly understand the first half of the code - it seems like it finds the coarse location of the returned value, using the fact that $\int_{k\cdot\ln2}^{(k+1)\cdot\ln2}e^{-x}dx=2^{-k-1}$. So if we stopped there, more or less, we would get a distribution which is a step function of the exp: enter image description here

I don't understand the 2nd half.

I guess I agree that $u$ (in this point) is again uniform distributed in $U(0, \ln 2)$ (the $\ln2$ comes from multiplying the result). I more or less understand that: "With probabilities ... $(\ln 2)^i/i!$ consider the minimum of $i$ uniform samples from $(0, \ln2)$". I just don't know how you show this is indeed an exponential distribution ... I know how to find the pdf of $\min u_i$ but here it has different probabilities for each $i$.

Can anyone show this?

Here is the code:

double exp_rand(void)
    /* q[k-1] = sum(log(2)^k / k!)  k=1,..,n, */
    /* The highest n (here 16) is determined by q[n-1] = 1.0 */
    /* within standard precision */
    const static double q[] =
    double a = 0.;
    double u = unif_rand();    /* precaution if u = 0 is ever returned */
    while(u <= 0. || u >= 1.) u = unif_rand();     for (;;) {
        u += u;
        if (u > 1.)
        a += q[0];
        u -= 1.;
        if (u <= q[0])
        return a + u;  # Up to here I understand
        int i = 0;
        double ustar = unif_rand(), umin = ustar;
        do {
        ustar = unif_rand();
        if (umin > ustar)
            umin = ustar;
        } while (u > q[i]);
        return a + umin * q[0];
  • $\begingroup$ The famous software is called R not r. Edited! $\endgroup$ Sep 6, 2022 at 13:25

1 Answer 1


The answer is (as usual!) provided in Devroyes' Non-uniform random variate generation (1986, p.396). The principle for the Ahrens-Dieter (1972) algorithm is a result due to George Marsaglia (1961):

Theorem IX.2.1$\ \ \ $If $U_1,U_2,\ldots$ is a series of iid $\mathcal U(0,1)$ random variables, if $Z$ is an independent positive Poisson $\mathcal P_+(\mu)$ random variable, and if $M$ is an independent Geometric $\mathcal G(1-e^{-\mu})$ random variable, then $$X = \mu(M+\min(U_1,\ldots,U_Z))\sim\mathcal E(1)$$

The proof proceeds as follows (reproducing from Devroye, p.395): \begin{align*} \mathbb P(\mu\min(U_1,\ldots,U_Z)\le x) &=\mathbb E^Z[\mathbb P(\mu\min(U_1,\ldots,U_z)\le x|Z=z)]\\ &=\mathbb E^Z[1-(1-x/\mu)^Z)]\\ &=1-\dfrac{e^{\mu-x}-1}{e^{\mu}-1}\\ &=\dfrac{1-e^{-x}}{1-e^{-\mu}} \end{align*} which is the cdf of the exponential distribution truncated to $(0,\mu)$ (see here). A preliminary result$^1$ due to Von Neumann (Lemma IV.2., p.125 & p.393) is that, when $$Z\sim \mathcal G(1-e^{-\mu})\qquad Y\sim \dfrac{e^{-y}}{1-e^{-\mu}}\mathbb I_{(0,1)}(y)$$then$$\mu(Z-1)+Y\sim\mathcal E(1)$$ This follows from considering the moment generating function (for $t<\min(\mu,1)$) $$\mathbb E[e^{tX}]=\mathbb E[e^{t\mu(Z-1)}]\mathbb E[e^{tY}]= \frac{1-e^{-\mu}}{1-e^{-\mu(1-t)}}\frac{1-e^{-\mu(1-t)}}{(1-e^{-\mu)}(1-t)}=\frac{1}{1-t}$$ which concludes the proof.

Relating to the graph included in the question, this means that $\mu M$ corresponds to the area under the step function, while $\mu\min(U_1,\ldots,U_Z))$ corresponds to the residual area between the step function and the Exponential $\mathcal E(1)$ density, which is therefore independent from $M$.

Marsaglia then derives his Exponential algorithm from Theorem IX.2.1:

  1. Generate a Geometric $M\sim\mathcal G(1-e^{-\mu})$ variable
  2. Generate two Uniform $U$ and $V$
  3. Set $Y=V$ and $Z=1$
  4. While $U>F_\mu(Z)$, increase $Z$ to $Z+1$ and decrease $Y$ to $\min(Y,W)$, where $W$ is Uniform
  5. Return $\mu(M+Y)$ as an Exponential $\mathcal E(1)$ variate.

and Ahrens and Dieter (1972) consists in an optimisation of the above, for instance by choosing$^2$ $\mu=\log(2)$, refining the generation of $M$, and storing the cdf $F(\cdot)$.

A detailed explanation of Ahrens and Dieter (1972) version:

  1. The prerecorded table corresponds to the first terms of the cdf of the Poisson $\mathcal P_+(\log 2)$ distribution (with q[0] equal to $\mu$)
  2. The first loop returns a rescaled Geometric $\mathcal G(1-e^{-\mu})$ variate $\mu M$ as a (using a sequential search inversion method and the property that $F_M(i)=1-2^{-i}$ when $\mu=\log 2$)
  3. The residual u-1 produces$^3$ an independent uniform variate $V$
  4. The case u <= q[0] corresponds to $Z=1$ and avoids running the second loop
  5. The second loop while (u > q[i]) produces the Poisson $\mathcal P_+(\log 2)$ variate $Z$ as i-1and the associated $Y=\min(U_1,\ldots,U_Z)$ as umin
  6. The outcome a+umin*q[0] is indeed $X = \mu(M+\min(U_1,\ldots,U_Z))\sim\mathcal E(1)$

$^1$An ingenious algorithm for generating from the Exponential distribution is derived from this lemma and consists in only producing sequences of Uniform variates $U_0,U_1,\ldots$.

$^2$When $\mu=\log 2$, the geometric random variate corresponds to the number of $0$ before the first $1$ in the binary expansion of $U\sim\mathcal U(0,1)$. Sampling directly these bits proves to be much faster than the first loop until (u > 1.).

$^3$Devroye remarks that "Ahrens and Dieter squeeze the first uniform [0,1] random variate $U$ dry". The efficiency of the method is such that it requires on average $1+\log(2)$ uniforms.

  • $\begingroup$ Obviously... But can you maybe elaborate for those it's not so obvious to them, like me :-) How do you show this? CDF? Forget the $\mu$ for a second, how do you show that $M+min(U_1,...U_Z) \sim Exp(1)$ ? $\endgroup$ May 1, 2021 at 19:43
  • $\begingroup$ ah ok, I didn't know about this bible, thanks! $\endgroup$ May 1, 2021 at 20:43
  • $\begingroup$ I can't find this in the link, and my university doesn't have access to full book. Also, maybe I got something wrong, but my calculations are: $\mathbb E^Z[1-(1-x)^Z)] = 1-\sum_{z=0}^\infty (1-x)^z \frac{e^{-1}}{z!} = 1-e^{-1}e^1e^{-x} = 1-e^{-x}$ $\endgroup$ May 1, 2021 at 21:06
  • $\begingroup$ Oh I got it, Z is from 1 up, without 0. $\endgroup$ May 1, 2021 at 21:29
  • $\begingroup$ Any chance of posting the Lemma IV 2? $\endgroup$ May 1, 2021 at 21:34

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.