Sampling Exactly from the Normal Distribution
CHARLES F. F. KARNEY, SRI International
An algorithm for sampling exactly from the normal distribution is given. The algorithm reads some number
of uniformly distributed random digits in a given base and generates an initial portion of the representation
of a normal deviate in the same base. Thereafter, uniform. random digits are copied directly into the rep-
resentation of the normal deviate. Thus, in contrast to existing methods, it is possible to generate normal
deviates exactly rounded to any precision with a mean cost that scales linearly in the precision. The method
performs no extended precision arithmetic, calls no transcendental functions, and uses no floating point
arithmetic whatsoever; it uses only simple integer operations. It can easily be adapted to sample exactly
from the discrete normal distribution whose parameters are rational numbers.
Categories and Subject Descriptors: G.3 [Probability and Statistics]: Random Number Generation
General Terms: Algorithms
Additional Key Words and Phrases: Random deviates, normal distribution, exact sampling
ACM Reference Format:
Charles F. F. Karney. 2016. Sampling exactly from the normal distribution. ACM Trans. Math. Softw. 42, 1,
Article 3 (January 2016), 14 pages.
DOI: http://dx.doi.org/10.1145/2710016
1. INTRODUCTION
Random variables with a normal density,
arewidelyusedinMonteCarlosimulations.Overthepast60years,scoresofalgorithms
for generating such normal deviates have been published [Thomas et al. 2007]. In this
article, I give another algorithm, Algorithm N, with the distinguishing feature that,
given a source of uniformly distributed random digits in some base b, it generates exact
normal deviates. In order to make the meaning of “exact” precise, consider Table I,
which shows the operation of the algorithm using b = 10.
The first column shows the random digits used by the algorithm, which, in this
example, are taken from successive lines of the table of random digits produced by the
RAND Corporation [1955], beginning at line 9077. Referring to the first line of this
table, the algorithm completes after reading the 7 random decimal digits, 9148686,
from the source and produces +1.6 as the initial portion of the decimal representation
of a normal deviate. At this point, random digits can be copied directly from the input
to the output, indicated by the ellipses (...) in the table; thus +1.6... represents a
uniform. random sample in the range (1.6,1.7). The next digits of the random sequence
Author’s address: C. F. F. Karney, SRI International, 201 Washington Rd, Princeton, NJ 08543-5300;
email: .
Permission to make digital or hard copies of part or all of this work for personal or classroom use is granted
without fee provided that copies are not made or distributed for profit or commercial advantage and that
copiesshowthisnoticeonthefirstpageorinitialscreenofadisplayalongwiththefullcitation.Copyrightsfor
components of this work owned by others than ACM must be honored. Abstracting with credit is permitted.
To copy otherwise, to republish, to post on servers, to redistribute to lists, or to use any component of this
work in other works requires prior specific permission and/or a fee. Permissions may be requested from
Publications Dept., ACM, Inc., 2 Penn Plaza, Suite 701, New York, NY 10121-0701 USA, fax +1 (212)
869-0481, or .
c©2016 ACM 0098-3500/2016/01-ART3 $15.00
DOI: http://dx.doi.org/10.1145/2710016
ACM Transactions on Mathematical Software, Vol. 42, No. 1, Article 3, Publication date: January 2016.
3:2 C. F. F. Karney
Table I. Sample Input and Output for Algorithm N with b = 10
input u-rand rounded
9148686 | 685171... +1.6... +1.668517(+)
2708 | 5545979... +0... +0.554598(−)
501446297 | 43871... +1.42... +1.424387(+)
065130319777860 | 96289... −0.76... −0.769629(−)
2736 | 0659086... +0... +0.065909(−)
Note: The input consists of uniformly distributed random decimal
digits. The algorithm reads the random digits before the vertical
bar and produces a normal deviate as a u-rand (given in the second
column), which is an initial portion of the decimal representation
of the normal deviate. Thereafter, random digits are copied directly
from the input (the digits after the vertical bar) into the decimal
fraction of the u-rand. The third column shows the result of adding
enough digits to allow the deviate to be rounded to 6 decimal digits;
the parenthetical sign indicates whether the magnitude of true
deviation is greater (+) or smaller (−) than the rounded result.
are 685171...,allowing the normal deviate to be exactly rounded to 6 decimal digits as
1.668517. I call the intermediate result, +1.6...,a “u-rand.” This can be thought of as
a partially sampled uniform. deviate. However, in conjunction with a source of random
digits, it is better to think of it as a compact representation of an arbitrary precision
random deviate. (The results in Table I are not “typical,” because the starting line in
the table of random digits was specifically chosen to limit the number of random digits
used.)
It is clear from this example that the method can be used to generate deviates that
satisfy the conditions of “ideal approximation” [Monahan 1985], that is, that the al-
gorithm is equivalent to sampling a real number from the normal distribution and
rounding it to the nearest representable floating point number. Furthermore, for ap-
plications requiring high-precision normal deviates, the new algorithm offers perfect
scaling: there is an amortized constant cost to producing the initial portion of the nor-
mal deviate; but, thereafter, the digits can be added to the result at a rate limited only
by the cost of producing and copying the random digits. Other sampling methods are
frequently referred to as “exact,” for example, the polar method [Box and Muller 1958]
and the ratio method [Kinderman and Monahan 1977]; but these are merely “accu-
rate to round off,” which, in practice, means only that the accuracy is commensurate
with the precision of the floating point number system. It’s possible to convert such
algorithms to obtain correctly rounded deviates; but this inevitably involves the use
of extended precision arithmetic. I will show, in Section 5, that Algorithm N performs
substantially better.
It is not immediately obvious that such an algorithm for exact sampling is possible.
However, in the early years of the era of modern computing, von Neumann [1951] pre-
sented a remarkably simple algorithm for sampling from the exponential density, in
which “the machine has in effect computed a logarithm by performing only discrimina-
tionsontherelativemagnitudeofnumbersin(0,1).”KnuthandYao[1976]showedthat
the algorithm can easily be adapted to generate exponential deviates that are exact.
The resulting method was extensively analyzed by Flajolet and Saheb [1986]. Sev-
eral authors have generalized von Neumann’s algorithm [Forsythe 1972; Ahrens and
Dieter 1973; Brent 1974; Monahan 1979]. However, these efforts entail using ordinary
floating point arithmetic; thus, the methods do not generate exact deviates.
In this article, I show that von Neumann’s algorithm can be extended to sample
exactly from the unit normal and discrete normal distributions. Although the resulting
algorithms are unlikely to displace existing methods for most applications, a nearly
ACM Transactions on Mathematical Software, Vol. 42, No. 1, Article 3, Publication date: January 2016.
Sampling Exactly from the Normal Distribution 3:3
optimal method is provided for generating normal deviates at arbitrary precision. In
addition, the ability to sample exactly from the discrete normal has applications to
cryptography because the security of cryptographic systems requires that any random
sampling be very accurate. Finally, the method is of theoretical interest as an example
of an algorithm in which exact transcendental results can be achieved with simple
integer arithmetic.
Implementations of the algorithms given in this article are available in ExRandom,
a small “header only” library for C++11, available at http://exrandom.sf.net/, and in
version 3.2.0 of MPFR [2016].
2. VON NEUMANN’S ALGORITHM
IbeginbyreviewingvonNeummann’salgorithm,becausethisisthebasisofthemethod
for sampling from the normal distribution.
Algorithm V (von Neumann). Samples E from the exponential distribution e