Which algorithm used by the rnorm function

Question

Which algorithm is using by the rnorm function by default to generate standard-normally distributed random numbers?

Answer 1

See ?RNGkind. The default is an inversion algorithm:

normal.kind can be "Kinderman-Ramage", "Buggy Kinderman-Ramage" (not for set.seed), "Ahrens-Dieter", "Box-Muller", "Inversion" (the default), or "user-supplied". (For inversion, see the reference in qnorm.) The Kinderman-Ramage generator used in versions prior to 1.7.1 (now called "Buggy") had several approximation errors and should only be used for reproduction of old results. The "Box-Muller" generator is stateful as pairs of normals are generated and returned sequentially. The state is reset whenever it is selected (even if it is the current normal generator) and when kind is changed.

You can change the algorithm by

RNGkind(normal.kind = "Box-Muller")

You can find what is currently set by looking at RNGkind()[2].

Answer 2

The other answer is sufficient, but left me with some more questions; in particular, I didn't see anywhere in the documentation* what on earth the "Inversion" algorithm is, so I dived into the source code, which also gives academic references to the papers originating the other possible algorithms, to figure out what exactly is being done.

    case INVERSION:
#define BIG 134217728 /* 2^27 */
    /* unif_rand() alone is not of high enough precision */
    u1 = unif_rand();
    u1 = (int)(BIG*u1) + unif_rand();
    return qnorm5(u1/BIG, 0.0, 1.0, 1, 0);

So it seems at base the default "Inversion" algorithm generates a high precision floating point number (looks like 53 bits, or the mantissa size for 64-bit floating numbers), then sends it to the qnorm5 function which is a CDF function for the normal distribution.

As to how the qnorm5 function works (given there is no closed form for the Normal CDF nor inverse CDF), I haven't had much luck cracking what seems to be the source code here, but they do give further academic references, namely Beasley, J. D. and S. G. Springer (1977) and Wichura, M.J. (1988); the former being typically used for small quantiles of the CDF and the latter for large (z>7 or so).

It may also be interesting to note that (as of this writing) this algorithm appears to be shared by the Julia language, which also shares the qnorm5 code used by R.

*To be fair, in retrospect, Wichura is mentioned in ?qnorm which is referenced above. Still it's worthwhile to spell things out in this thread, I think.