# Why doesn't R use Inverse Transform Sampling to sample from the Exponential Distribution?

I was reading this question about the algorithm that R uses to sample from the Exponential($$\lambda$$) distribution. It looks like R uses the Ahrens-Dieter algorithm to sample from the exponential distribution. The reference for it is a paper from 1972 that makes a surprising claim that this algorithm is a little faster (in terms of microseconds, not in terms of Big O computational cycles) than just delivering $$\frac{-\ln(u)}{\lambda}$$ (where U~Uniform(0,1); as the Inverse Transform would have us do).

I find this very surprising. First of all, I am not sure the results from 1972 carry much weight in 2022. Secondly, why the need for an algorithm in the first place as opposed to just delivering $$\frac{-\ln(u)}{\lambda}$$ (as the Inverse Transform would have us do)? How slow can calculating $$\frac{-\ln(u)}{\lambda}$$ be? Maybe it was slow in 1972. Or maybe there is another reason?

• Are you sure it's actually using the 1972 algorithm? There are later algorithms as well:, for example, "AHRENS, J. g. AND DIETER, U. 1988. Efficient table-free sampling methods for the exponential, Cauchy, and normal distributions. Commun. ACM 31, 11 (Nov.), 1330-1337" and a paper by Hamilton (Algorithm 780: exponential pseudorandom distribution in ACM Transactions on Mathematical Software, Volume 24, Issue 1 March 1998 pp 102–1060) which claims to have fixed a bug in the 1988 paper. Sep 27, 2022 at 21:22
• In R I just typed ?rexp ... then if I scroll down I see: rexp uses Ahrens, J. H. and Dieter, U. (1972). Computer methods for sampling from the exponential and normal distributions. Communications of the ACM, 15, 873–882. That said maybe the core R team referenced the wrong article here; hard to say. But at least that is what the R documentation says. Sep 27, 2022 at 21:24
• Hmmm... you might want to check the timings in the Hamilton paper (with some comments in the text about compilers and chip architecture) just for your interest. Link here: dl.acm.org/doi/epdf/10.1145/285861.285866 - page 105. Sep 27, 2022 at 21:26
• The questions about historical computing equipment are interesting but off topic. Why not just ask R for the answer? On my system, the log-uniform method is 60+% slower. with(data.frame(n = 1e8, lambda = pi), { print(system.time(x <- rexp(n, 1/lambda))); print(system.time(x <- -log(runif(n)) * lambda)) })
– whuber
Sep 27, 2022 at 21:43
• @ColorStatistics if it helps others, the manual page you are looking at is at stat.ethz.ch/R-manual/R-devel/library/stats/html/… Sep 28, 2022 at 13:27

I have written rather a lot of code to test this on two machines, the first of which is a Windows machine that has an Intel(R) Core(TM) i7-9700K CPU @ 3.60GHz processor. I am using a recent MinGW compiler via the CodeBlocks IDE, with optimizations set at O3 and expensive_optimizations enabled, but no chip-specific optimizations.

I have run the same code on an Intel(R) Xeon CPU E5-2650L v4 @ 1.70GHz processor, Ubuntu OS, and GCC 7.5 compiler, with the same compiler options.

I am using the pcg32 random number generator to generate the required uniform variates inside the code helpfully provided by @Alex. pcg32 is extremely fast - my implementation, copied from who knows where, takes only 18% longer than the C library rand() function while returning an unsigned integer between 0 and 4294967296$$= 2^{32}$$, whereas rand() returns a signed integer between 0 and 32767 (much poorer granularity), and has excellent properties. See https://www.pcg-random.org/ for more.

Fast Windows machine: Generating 10 million variates using the pcg32 version of the code that R uses took 568,095 microseconds, including the overhead induced by the for loop.

Generating 10 million variates using the inverse probability transform took 527,346 microseconds, including the overhead induced by the for loop. This is roughly 93% of the time that the A-D algorithm uses. Roughly 40% of the time for the inverse probability transform algorithm appears to be loop overhead and the uniform RNG.

Slower Linux machine: the A-D algorithm took 666,358 microseconds, the inverse probability transform took only 555,192 microseconds. This is roughly 83% of the time that the A-D algorithm uses.

These results certainly validate the OP's suspicion that things may have changed since the 1970s. Regardless of the algorithm choice, being able to generate 15-20 million exponential random numbers in one second on one thread is certainly a nice capability to have!

One interesting finding is that the runtimes don't scale with the CPU speed; they are a little greater on the Linux box than on the Windows box, but not nearly the 2x+ difference in speed. Of course, the compilers are different, and GHz is by no means the sole influence on runtime.

I'd be happy to post the code, but there's about 130 lines of it, including some comments and blank lines. Thoughts?

• I'm pretty sure a key difference between machines in the 60's and 70's and those we have today is transcendental functions like exp, log, etc. were not implemented in microcode on the chip. Fortran (on which R was built) was the language of choice for scientific computing partly because it included a numeric function library which included optimized, high-accuracy implementations of these functions in software. (Expert Fortran programmers understood this and distinguished such "internal" functions from "externals.") Avoiding a call to a logarithm might have been a big savings.
– whuber
Sep 28, 2022 at 13:35
• @whuber: very insightful; it is starting to make sense why they went the route they did Sep 28, 2022 at 16:01
• @whuber: I'd be very interested in reading more about what you wrote in this comment. If you can think of any resources I could consult to that end, please let me know. Sep 28, 2022 at 21:39
• @Color Wikipedia has an article on the IBM System 360 architecture. It's not very detailed, though. I found a list of machine instructions at en.wikibooks.org/wiki/360_Assembly/360_Instructions. It doesn't appear to include anything more complicated than a floating point square root. The ACM has a copy of the machine manual Principles of Operation: Google "IBM system/360 principles of operation ACM". (One often got started in computing by reading manuals like this.)
– whuber
Sep 28, 2022 at 21:51
• ibm.com/ibm/history/exhibits/mainframe/mainframe_profiles.html is a nice resource concerning IBM mainframe architectures and capabilities from the early 1950's through the 1980's.
– whuber
Sep 28, 2022 at 22:02

Just for the review sake, the source code for this function in R (written in the C language):

#include "nmath.h"

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[] = {
0.6931471805599453,
0.9333736875190459,
0.9888777961838675,
0.9984959252914960,
0.9998292811061389,
0.9999833164100727,
0.9999985691438767,
0.9999998906925558,
0.9999999924734159,
0.9999999995283275,
0.9999999999728814,
0.9999999999985598,
0.9999999999999289,
0.9999999999999968,
0.9999999999999999,
1.0000000000000000
};

double a = 0.0;
double u = unif_rand();    /* precaution if u = 0 is ever returned */
while(u <= 0.0 || u >= 1.0) u = unif_rand();
for (;;) {
u += u;
if (u > 1.0) break;
a += q[0];
}
u -= 1.0;

if (u <= q[0]) return a + u;

int i = 0;
double ustar = unif_rand(), umin = ustar;
do {
ustar = unif_rand();
if (umin > ustar) umin = ustar;
i++;
} while (u > q[i]);
return a + umin * q[0];
}
$$$$
`