My question is about R's built-in exponential random number generator, the function rexp(). When trying to generate exponentially distributed random numbers, many textbooks recommend the inverse transform method as outlined in this Wikipedia page. I am aware that there are other methods to accomplish this task. In particular, R's source code uses the algorithm outlined in a paper by Ahrens & Dieter (1972).

I have convinced myself that the Ahrens-Dieter (AD) method is correct. Still, I do not see the benefit of using their method compared to the inverse transform (IT) method. AD is not only more complex to implement than IT. There does not seem to be a speed benefit either. Here is my R code to benchmark both methods followed by the results.

invTrans <- function(n)
    -log(runif(n))
print("For the inverse transform:")
print(system.time(invTrans(1e8)))
print("For the Ahrens-Dieter algorithm:")
print(system.time(rexp(1e8)))

Results:

[1] "For the inverse transform:"

user system elapsed

4.227 0.266 4.597

[1] "For the Ahrens-Dieter algorithm:"

user system elapsed

4.919 0.265 5.213

Comparing the code for the two methods, AD draws at least two uniform random numbers (with the C function unif_rand()) to get one exponential random number. IT only needs one uniform random number. Presumably the R core team decided against implementing IT because it assumed that taking the logarithm may be slower than generating more uniform random numbers. I understand that the speed of taking logarithms can be machine-dependent, but at least for me the opposite is true. Perhaps there are issues around IT’s numerical precision having to do with the singularity of the logarithm at 0? But then, the R source code sexp.c reveals that the implementation of AD also loses some numerical precision because the following portion of C code removes the leading bits from the uniform random number u.

double u = unif_rand();
while(u <= 0. || u >= 1.) u = unif_rand();
for (;;) {
    u += u;
    if (u > 1.)
        break;
    a += q[0];
}
u -= 1.;

u is later recycled as a uniform random number in the remainder of sexp.c. So far, it appears as if

• IT is easier to code,
• IT is faster, and
• both IT and AD possibly lose numerical accuracy.

I would really appreciate if someone could explain why R still implements AD as the only available option for rexp().