ISO 28640:2010

INTERNATIONAL ISO
STANDARD 28640
First edition
2010-03-15
Random variate generation methods
Méthodes de génération de nombres pseudo-aléatoires

Reference number
©
ISO 2010
PDF disclaimer
This PDF file may contain embedded typefaces. In accordance with Adobe's licensing policy, this file may be printed or viewed but
shall not be edited unless the typefaces which are embedded are licensed to and installed on the computer performing the editing. In
downloading this file, parties accept therein the responsibility of not infringing Adobe's licensing policy. The ISO Central Secretariat
accepts no liability in this area.
Adobe is a trademark of Adobe Systems Incorporated.
Details of the software products used to create this PDF file can be found in the General Info relative to the file; the PDF-creation
parameters were optimized for printing. Every care has been taken to ensure that the file is suitable for use by ISO member bodies. In
the unlikely event that a problem relating to it is found, please inform the Central Secretariat at the address given below.

©  ISO 2010
All rights reserved. Unless otherwise specified, no part of this publication may be reproduced or utilized in any form or by any means,
electronic or mechanical, including photocopying and microfilm, without permission in writing from either ISO at the address below or
ISO's member body in the country of the requester.
ISO copyright office
Case postale 56 • CH-1211 Geneva 20
Tel. + 41 22 749 01 11
Fax + 41 22 749 09 47
E-mail copyright@iso.org
Web www.iso.org
Published in Switzerland
ii © ISO 2010 – All rights reserved

Contents Page
Foreword .iv
Introduction.v
1 Scope.1
2 Normative references.1
3 Terms and definitions .1
4 Symbols and mathematical binary operations.2
4.1 Symbols.2
4.2 Mathematical binary operations .2
5 Uniformly distributed pseudo-random numbers.3
5.1 General .3
5.2 M-sequence method definition .3
5.3 Pentanomial GFSR method .4
5.4 Combined Tausworthe method.4
5.5 Mersenne Twister method .5
6 Generation of random numbers from various distributions.6
6.1 Introduction.6
6.2 Uniform distribution .6
6.3 Standard beta distribution.7
6.4 Triangular distribution .8
6.5 General exponential distribution with location and scale parameters .8
6.6 Normal distribution .9
6.7 Gamma distribution.9
6.8 Weibull distribution .11
6.9 Lognormal distribution .11
6.10 Logistic distribution .11
6.11 Multivariate normal random variate generation .12
6.12 Binomial distribution.12
6.13 Poisson distribution.14
6.14 Discrete uniform distribution .14
Annex A (informative) Table of physical random numbers .16
A.1 Table of random numbers .16
A.2 Method of physical random number generation.17
Annex B (informative) Algorithm for pseudo-random number generation.18
B.1 Program code for the trinomial GFSR method.18
B.2 Program code for the pentanomial GFSR method.22
B.3 Program code for the combined Tausworthe method.28
B.4 Program code for the Mersenne Twister method .32
B.5 Linear congruential method .38
B.6 Reference examples.52
Bibliography.53

Foreword
ISO (the International Organization for Standardization) is a worldwide federation of national standards bodies
(ISO member bodies). The work of preparing International Standards is normally carried out through ISO
technical committees. Each member body interested in a subject for which a technical committee has been
established has the right to be represented on that committee. International organizations, governmental and
non-governmental, in liaison with ISO, also take part in the work. ISO collaborates closely with the
International Electrotechnical Commission (IEC) on all matters of electrotechnical standardization.
International Standards are drafted in accordance with the rules given in the ISO/IEC Directives, Part 2.
The main task of technical committees is to prepare International Standards. Draft International Standards
adopted by the technical committees are circulated to the member bodies for voting. Publication as an
International Standard requires approval by at least 75 % of the member bodies casting a vote.
Attention is drawn to the possibility that some of the elements of this document may be the subject of patent
rights. ISO shall not be held responsible for identifying any or all such patent rights.
ISO 28640 was prepared by Technical Committee ISO/TC 69, Applications of statistical methods.
This is the first edition.
iv © ISO 2010 – All rights reserved

Introduction
This International Standard specifies typical algorithms by which the users can regard the generated
numerical sequences as if they were real random variates.
Nowadays most statisticians, scientists and engineers have enough computer power at their disposal to carry
out large computer simulations, and it is important that these be based on sound pseudo-random generators.
This International Standard has been developed to help ensure that randomization, where needed, is carried
out correctly and efficiently.
Six uses of randomization can be identified in statistical standardization:
⎯ selection of a random sample;
⎯ analysis of sample data;
⎯ development of standards;
⎯ checking theoretical results;
⎯ demonstrating that a proposed procedure has the properties claimed of it;
⎯ resolving uncertainty in the statistical literature.

INTERNATIONAL STANDARD ISO 28640:2010(E)

Random variate generation methods
1 Scope
This International Standard specifies methods for generating uniform and non-uniform random variates for
Monte Carlo simulation purposes. Cryptographic random number generation methods are not included. This
International Standard is applicable, inter alia, by
⎯ researchers, industrial engineers or experts in operations management, who use statistical simulation,
⎯ statisticians who need randomization related to SQC methods, statistical design of experiments or sample
surveys,
⎯ applied mathematicians who plan complex optimization procedures that require the use of Monte Carlo
methods, and
⎯ software engineers who implement algorithms for random variate generation.
2 Normative references
The following referenced documents are indispensable for the application of this document. For dated
references, only the edition cited applies. For undated references, the latest edition of the referenced
document (including any amendments) applies.
ISO/IEC 2382-1, Information technology — Vocabulary — Part 1: Fundamental terms
ISO 3534-1, Statistics — Vocabulary and symbols — Part 1: General statistical terms and terms used in
probability
ISO 3534-2, Statistics — Vocabulary and symbols — Part 2: Applied statistics
3 Terms and definitions
For the purposes of this document, the terms and definitions given in ISO/IEC 2382-1, ISO 3534-1 and
ISO 3534-2 apply, except where redefined below.
3.1
random variate
random number
number as the realization of a specific random variable
NOTE 1 The term “random number” is often used for uniformly distributed random variate.
NOTE 2 Random numbers provided as a sequence are called a “random number sequence”.
3.2
pseudo-random number
random number (3.1) generated by an algorithm, that appears to be random
NOTE If there is no fear of misunderstanding, a pseudo-random number may simply be called a “random number”.
3.3
physical random number
random number (3.1) generated by a physical mechanism
3.4
binary random number sequence
random number (3.1) sequence consisting of zeros and ones
3.5
seed
initialization value required for pseudo-random number generation
4 Symbols and mathematical binary operations
4.1 Symbols
For the purposes of this document, the symbols given in the normative references as the latest versions of
ISO/IEC 2382-1, ISO 3534-1 and ISO 3534-2 apply, except where redefined below.
The symbols and abbreviations specifically used in this International Standard are as follows:
X integer type uniform random number
U standard uniform random number
Z normal random variate
n suffix of random number sequence
4.2 Mathematical binary operations
The mathematical binary operations specifically used in this International Standard are as follows:
mod(m; k) residue from dividing integer m by k
m ⊕ k bitwise exclusive logical disjunction of binary integers m and k
EXAMPLE 1 1 ⊕ 1 = 0
0 ⊕ 1 = 1
1 ⊕ 0 = 1
0 ⊕ 0 = 0
1010 ⊕ 1100 = 0110
2 © ISO 2010 – All rights reserved

m ∧ k bitwise logical conjunction of binary integers m and k
EXAMPLE 2 1 ∧ 1 = 1
0 ∧ 1 = 0
1 ∧ 0 = 0
0 ∧ 0 = 0
1010 ∧ 1100 = 1000
m := k replaces value m by k
m >> k k-bit right shift of binary integer m
m << k k-bit left shift of binary integer m
5 Uniformly distributed pseudo-random numbers
5.1 General
This clause provides algorithms for generating uniformly distributed pseudo-random numbers based on
M-sequence methods (see 5.2).
Annex A introduces the concept of physically generated random numbers for information.
Annex B includes C and full Basic codes for all the recommended algorithms for information. Although the
linear congruential method is not recommended for complex Monte Carlo simulations, it is also included in
Annex B for information.
5.2 M-sequence method definition
a) Let p be a natural number, and c , c , ., c be specified to be 0 or 1, and define the recurrence
1 2 p − 1
formula
x = c x + c x + . + c x + x (mod 2)  (n = 1, 2, 3, .)

n + p p − 1 n + p − 1 p − 2 n + p − 2 1 n + 1 n
b) The least positive integer N such that x = x for all n is called the period of the sequence. This
n + N n
p
sequence is called an M-sequence in cases where its period is 2 − 1.
c) The polynomial
p p − 1
t + c t + . + c t + 1
p − 1 1
is called the characteristic polynomial of the above-mentioned recurrence formula.
NOTE 1 A necessary and sufficient condition for the above-mentioned recurrence formula to generate an M-sequence
is that at least one of the seeds x , x , ., x is not zero.
1 2 p
NOTE 2 The letter M of the M-sequence originates from the English word “maximum”, which means the largest. The
p
period of any sequence generated by the above recurrence formula cannot exceed 2 − 1. Therefore, if there is a series
p
that has a period of 2 − 1, it is the series that has the largest period.
NOTE 3 When this method is used, either one of the polynomials listed in Table 1 or another primitive polynomial listed
in the literature is chosen as the characteristic polynomial and its coefficients are used to define the recurrence formula
in a).
5.3 Pentanomial GFSR method
This method uses a characteristic polynomial of 5 terms, and it generates binary integer sequences of w bits
by the following recurrence formula. This algorithm is called the GFSR or “generalized feedback shift register”
random number generator.
X = X ⊕ X ⊕ X ⊕ X  (n = 1, 2, 3, .)
n + p n + q1 n + q2 n + q3 n
The parameters are (p, q1, q2, q3, w) and X , ., X are initially given as seeds. Examples of parameters p, q ,
1 p 1
p
q , q giving the largest period 2 − 1 are indicated in Table 1.
2 3
Table 1 — Pentanomial characteristic polynomials
p q q q
1 2 3
89 20 40 69
107 31 57 82
127 22 63 83
521 86 197 447
607 167 307 461
1 279 339 630 988
2 203 585 1 197 1 656
2 281 577 1 109 1 709
3 217 809 1 621 2 381
4 253 1 093 2 254 3 297
4 423 1 171 2 273 3 299
9 689 2 799 5 463 7 712
NOTE q , q , q represent exponents of the non-zero terms of the
1 2 3
characteristic polynomial.
5.4 Combined Tausworthe method
Let x , x , x , … be an M-sequence generated by the recurrence relationship:
0 1 2
x = x + x (mod 2)  (n = 0, 1, 2, …)
n + p n + q n
Using this M-sequence, a w-bit integer sequence called a simple Tausworthe sequence with parameters
(p, q, t) is obtained as follows:
X = x x …x  (n = 0, 1, 2, …)
n nt nt + 1 nt + w − 1
where
p
t is a natural number which is coprime to the period 2 − 1 of the M-sequence;
w is the word length which does not exceed p.
p
The period of this sequence is also 2 − 1.
NOTE 1 Two integers are said to be coprime, or relatively prime, when they have no common divisors other than unity.
4 © ISO 2010 – All rights reserved

EXAMPLE If a primitive polynomial t + t + 1 is chosen, set p = 4, and q = 1 in the above recurrence relationship. If
the seeds (x , x , x , x ) = (1,1,1,1) are given to the recurrence, then the M-sequence obtained by the recurrence will be
0 1 2 3
1,1,1,1, 0,0,0,1, 0,0,1,1, 0,1,0,1, 1,1,1,0, … , and the period of the sequence is 2 − 1 = 15. Taking, for example, t = 4
which is coprime to 15, and w = 4, the simple Tausworthe sequence {X } with parameters (4, 1, 4) is obtained as follows:
n
X = x x x x = 1111 (= 15)
0 0 1 2 3
X = x x x x = 0001 (= 1)
1 4 5 6 7
X = x x x x = 0011 (= 3)
2 8 9 10 11
X = x x x x = 0101 (= 5)
3 12 13 14 0
X = x x x x = 1110 (= 14)
4 1 2 3 4
X = x x x x = 0010 (= 2)
5 5 6 7 8
.....
The simple Tausworthe sequence obtained in this way will be, in decimal notation, 15, 1, 3, 5, 14, 2, 6, 11, 12,
4, 13, 7, 8, 9, 10, 15, 1, 3, … , and its period is 2 − 1 = 15.
(j)
Suppose now that there is a multiple, say J, of simple Tausworthe sequences {X }, j = 1, 2, ., J with the
n
same word length w. The combined Tausworthe method is a technique that generates a sequence of
pseudo-random numbers {X } as the bitwise exclusive logical disjunction in the binary representation of these
n
J sequences.
(1) (2) (J)
X = X ⊕ X ⊕ … ⊕ X  (n = 0, 1, 2, …)
n n n n
The parameters and the seeds of the combined Tausworthe sequence are combinations of the parameters
and the seeds of each simple Tausworthe sequence. If the periods of the J simple Tausworthe sequences are
coprime, then the period of the combined Tausworthe sequence is the product of the periods of the J
sequences.
NOTE 2 This method can generate sequences with good multidimensional equidistribution characteristics. The
algorithm taus88_31( ) given in Annex A generates a sequence of 31-bit integers by combining three simple Tausworthe
generators with parameters (p, q, t) = (31, 13, 12), (29, 2, 4), and (28, 3, 17), respectively. The period length of the
31 29 28 88
combined sequence is (2 − 1)(2 − 1)(2 − 1), i.e. about 2 . Many other combinations are suggested in References [7]
and [8] in the Bibliography.
5.5 Mersenne Twister method
Let X be a binary integer of w bits. Then, the Mersenne Twister method generates a sequence of binary
n
integer pseudo-random numbers of w bits according to the following recurrence formula with integer constants
p, q, r and a binary integer a of w bits.
f l (r)
X = X ⊕ (X |X ) A ,  (n = 1, 2, 3, .)
n + p n + q n n+1
f l (r) f l
where (X |X ) represents a binary integer that is obtained by a concatenation of X and X , the first
n n+1 n n + 1
w − r bits of X and the last r bits of X in this order. A is a w × w 0-1 matrix, which is determined by a, and
n n + 1
the product XA is given by the following formula.
X >> 1 (when the last bit of X = 0)
XA = (X >> 1) ⊕ a (when the last bit of X = 1)
Here, X is regarded as a w dimensional 0-1 vector.
pw−r
NOTE The necessary amount of memory for this computation is p words, the period becomes 2 − 1, and the
efficiency is better than that of the GFSR methods described previously. To improve the randomness of the first w − r bits,
the following series of conversions can be applied to X .
n
y := X
n
y := y ⊕ (y >> u)
y := y ⊕ [(y << s) ∧ b]
y := y ⊕ [(y << t) ∧ c]
y := y ⊕ (y >> l)
where b, c are constant bits masks to improve the randomness of the first w − r bits. The parameters of this
algorithm are (p, q, r, w, a, u, s, t, l, b, c). The seeds are X , ., X and the first w − r bits of X .
2 q + 1 1
The final value of y is the pseudo-random number.
6 Generation of random numbers from various distributions
6.1 Introduction
Methods of generating random numbers Y from various distributions by using uniform random numbers X of
integer type, are described below.
The distribution function is denoted by F(y). If it is a continuous distribution, its probability density function is
denoted by f(y), and if it is a discrete distribution, its probability mass function is denoted by p(y).
6.2 Uniform distribution
6.2.1 Standard uniform distribution
6.2.1.1 Probability density function
⎧1, 0uuy 1
fy() =
⎨
0, otherwise
⎩
6.2.1.2 Random variate generation method
If the maximum value of uniform random number X of integer type is m − 1, the following formula should be
used to generate standard uniform random numbers.
X
U =
m
EXAMPLE For any w-bit integer sequences generated by the method described in 5.2 through 5.5, m = 2.
NOTE 1 Because X takes on discrete values, the values of U are also discrete.
NOTE 2 The value of U never becomes 1,0. The value of U becomes 0,0 only when X = 0. In the case of M-sequence
random numbers, any generation method may cause this phenomenon.
NOTE 3 Random numbers from the standard uniform distribution are called standard uniform random numbers, and
are represented by U , U , … They are assumed to be independent of each other.
1 2
6 © ISO 2010 – All rights reserved

6.2.2 General uniform distribution
6.2.2.1 Probability density function
⎧1/ba, uuy a +b
fy() =
⎨
0, otherwise
⎩
where b > 0.
6.2.2.2 Random variate generation method
If the standard uniform random number U is generated by the method specified in 6.2.1.2, then the general
uniform random number should be generated by the following formula:
Yb=+U a
6.3 Standard beta distribution
6.3.1 Probability density function
d −1
⎧ c−1
yy1−
()
⎪
01uuy
,
fy() =
⎨
Β cd,
()
⎪
otherwise
0,
⎩
d −1
c−1
where Β cd,1=−x xdx is the beta function and the parameters c and d are greater than 0.
() ( )
∫
6.3.2 Random variate generation method by Jöhnk
If the standard uniform random numbers U and U are independently generated by the method specified in
1 2
6.2.1, then the standard beta random number Y should be generated by the following procedures.
1/cd1/ 1/ c
 
If YU=+U is less than or equal to 1, set YU= /Y ; otherwise, generate two standard uniform
12 1
random numbers again until the inequality is satisfied.
6.3.3 Random variate generation method by Cheng
If the standard uniform random numbers U and U are independently generated by the method specified in
1 2
6.2.1, then the standard beta random number Y should be generated by the following procedures.
[Set-up]
a) Let
⎧ < 1
min()cd, ,      if min(c,d)
⎪
q =
⎨
2(cd−+c d)
,     otherwise
⎪
cd+− 2
⎩
[Generation]
b) Let
1 U
VW==,ecxp(V)
qU1−
c) If
⎛⎞cd+
()cd++ln (c+q)V−ln4Wln(UU)
⎜⎟ 12
dW+
⎝⎠
then employ
W
Y = ; and stop.
dW+
d) Generate U , U , and go to b).
1 2
U
1 ⎛⎞cd+ W
Vc=+()dln +(c+q)V−ln4Wln(UU)
⎜⎟ 12
q 1−+U dW dW+
⎝⎠
Jöhnk's method is recommended when max(c, d) u 1; otherwise, Cheng's method is recommended.
NOTE General beta random variates with the support [a, a + b] will be obtained by a linear transformation similar to
the one described in 6.2.2.2.
6.4 Triangular distribution
6.4.1 Probability density function
⎧
ba−−y
⎪
, ab−+uuy a b
fy() = 2
⎨
b
⎪
0,             otherwise
⎩
where b > 0.
6.4.2 Random variate generation method
If the standard uniform random numbers U and U are independently generated by the method specified in
1 2
6.2.1, then the triangular random number Y should be generated by Ya=+b()U +U −1 .
6.5 General exponential distribution with location and scale parameters
6.5.1 Probability density function
⎧
exp{}−−()ya/b , yWa
⎪
fy() = b
⎨
⎪
⎩0,              ya<
where a and b are the location and scale parameters of the exponential distribution, respectively.
6.5.2 Random variate generation method
If the standard uniform random number U is generated by the method specified in 6.2.1, then the general
exponential random number should be generated by
Y = − b ln(U) + a
8 © ISO 2010 – All rights reserved

6.6 Normal distribution
6.6.1 Probability density function
⎧⎫
11 2
fz=−exp z− µ ,−∞ () ()
⎨⎬
2πσ
⎩⎭2σ
where µ and σ are the mean and standard deviation of the normal distribution, respectively.
NOTE The symbol Z is used for a normal random variate.
6.6.2 The Box-Müller method
If the standard uniform random numbers U and U are independently generated by the method specified in
1 2
6.2.1, then two independent normal random numbers Z , Z will be generated by the following procedures:
1 2
Z =+µσ −2ln 1−UUcos(2π )
()
11 2
Z 2=+µσ −2ln 1−UUsin(2π )
()
NOTE 1 Since U is not continuous, Z , Z cannot be normally distributed in a strict sense. For example, using this
1 1 2
procedure, the upper bound of the absolute value of the pseudo-standardized normal variates is
−1
32 31
M=−2ln()mm=2In ; thus, when m = 2 , M ≈ 6,660 4, and when m = 2 − 1, M ≈ 6,555 5. However, since the
−10
probability that absolute values of random variables from a true standard normal distribution exceed M is only about 10 ,
this will rarely be a problem in practice.
NOTE 2 When generating U , U , by a linear congruential method sequentially, U and U are not independent of each
1 2 1 2
other, so the tail of the distribution of the generated Z and Z can depart considerably from the true normal distribution.
1 2
6.7 Gamma distribution
6.7.1 Probability density function
⎧ 1
c−1
()ya−−/b exp(y−a)/b ,  yWa
{} { }
⎪
fy() = bcΓ()
⎨
⎪
se
⎩0,                         otherwi
where a, b, c are the location, scale and shape parameters of the distribution, respectively.
6.7.2 Random variate generation methods
6.7.2.1 General
Algorithms are given for three special cases depending on the shape parameter value c as follows.
6.7.2.2 Case of c = k (k: integer)
Using independent uniform random numbers U , U , . , U , the transformation formula
1 2 k
Ya=−b ln{(1−U )(1−U ).(1−U )}
12 k
should be used.
NOTE The chi-squared distribution with even degrees of freedom can be generated by this method when a = 0 and
b = 2.
6.7.2.3 Case of c = k + 1/2 (k: integer)
Using a standard normal random number Z and uniform random number U , U , . , U , the transformation
0 1 2 k
formula
⎡⎤
Ya=+bZ /2l−n (1−U )(1−U ).(1−U )
{}
012 k
⎣⎦
should be used. In the case where k = 0, the logarithm term disappears.
NOTE The chi-squared distribution with odd degrees of freedom can be generated by this method when a = 0 and
b = 2.
6.7.2.4 Approximate generation method when c > 1/3
a) Set rc=−1/3, , s= rt=r−rln()r, p=1/(3 s) and qr=−3 .
b) Generate standard normal random number Z.
c) If Z < q, then go to b).
d) Set Yp=+()Z s, 2V=Z/, and generate U.
e) If (Y − r) /Y − V u U, then employ Y := a + bY and stop.
f) Set W = Y − r ln(Y) − t − V.
g) If W u U, then employ Y := a + bY and stop.
h) If W > −In(1,0 − U), then go to b).
NOTE This method is based on the Wilson-Hilferty transformation of chi-square distributions to an approximate
standard normal distribution. The accuracy of this approximation depends on the parameter values of c, and the
dependency is rather complicated. A very rough idea is as follows: the absolute difference between the percentage point
of the approximate distribution and the exact distribution is always less than 0,2.
6.7.2.5 Exact generation method when c > 1/2, by Cheng
a) Set q = c – ln 4 and rc=+ 2c−1.
b) Generate standard uniform random numbers U and U .
1 2
U
1 2
c) Set Vc==ln(),W c exp(U),Z=U U and R = q + rV − W.
1−U
d) If R W 4,5Z − (1 + In 4,5) then employ Y = a + bW and stop.
e) If R W InZ then employ Y = a + bW and stop.
f) Go to b).
10 © ISO 2010 – All rights reserved

p
⎛⎞⎛ ⎞
UU
pc=1/2 −1, lq=−cn4 r=+c 2c−1 q+prln −c W4,5UU −(1+ln4,5)
⎜⎟⎜ ⎟
()12
11−−UU
⎝⎠11⎝ ⎠
p
⎛⎞
U
Ya=+bc
⎜⎟
1−U
⎝⎠1
6.8 Weibull distribution
6.8.1 Probability distribution function
c
⎧
⎧⎫
ya−
⎪⎪⎛⎞
⎪1e−−xp , yaW
⎪⎨ ⎬
⎜⎟
Fy() = b
⎝⎠
⎨
⎪⎪
⎩⎭
⎪
0,               ya<
⎪
⎩
where a, b and c are the location, scale and shape parameters of the distribution, respectively.
6.8.2 Random variate generation method
If the standard uniform random numbers U are generated by the method specified in 6.2.1, then general
Weibull random numbers are generated by the following formula:
1/c
Y = a – b{ln(1 – U)}
6.9 Lognormal distribution
6.9.1 Probability density function
⎧
⎧⎫
11 ya−
⎪⎪⎛⎞
⎪
exp − , yaW
⎨⎬
⎜⎟
fy =
() ⎨
2 b
2{π−ya b} ⎝⎠
()
⎪⎪
⎩⎭
⎪
⎩
0,                      ya<
where a and b are the location and scale parameters of the normal distribution, respectively.
6.9.2 Random variate generation method
Using standard normal random numbers Z,
Ya=+ exp()bZ
generates lognormal random numbers.
6.10 Logistic distribution
6.10.1 Probability function
Fy()=−, ∞ 1e+−xp{}()ya− /b
where a and b are the location and scale parameters of the distribution, respectively.
6.10.2 Random variate generation method
If standard uniform random numbers U are generated by the method specified in 6.2.1, then logistic random
numbers are generated by
U
⎛⎞
Ya=+b ln
⎜⎟
1−U
⎝⎠
6.11 Multivariate normal random variate generation
Random numbers Y , Y , ., Y from an n-dimensional normal distribution, with mean values µ , µ , ., µ and
1 2 n 1 2 n
variances and covariances σ (1 u i, j u n) are generated as follows using mutually independent standard
ij
normal random numbers Z , ., Z .
1 n
Y = µ + a Z
1 1 11 1
Y = µ + a Z + a Z
2 2 21 1 22 2
...
Y = µ + a Z + a Z + . + a Z
n n n1 1 n2 2 nn n
where a , ., a are constants that are calculated before start of random number generation from variances
11 nn
and covariances by following Cholesky factorization procedures, as given below.
NOTE σ (1 u i, j u n), σ means the variance of Y .
ij ii i
a) Set a = σ , a = σ /a (2 u i u n)
11 11 i1 i1 11
b) For j = 2, ., n set
j−1
⎛⎞ 2
⎜⎟
aa=−σ
jj jj ∑ jk
⎜⎟
k=1
⎝⎠
and
j−1
⎛⎞
⎜⎟
aa=−σ a /(aj+1uuin)
ij ij ik jk jj
∑
⎜⎟
k =1
⎝⎠
6.12 Binomial distribution
6.12.1 Probability mass function
When some event occurs with probability p at each trial, the probability that this event occurs y times in n trials
is given by the following formula:
⎛⎞n
yn−y
p(yp)=−(1p) ,y= 0,1, .,n
⎜⎟
y
⎝⎠
where 0 < p < 1.
12 © ISO 2010 – All rights reserved

6.12.2 Random variate generation methods
6.12.2.1 General
The following methods should be used for generating random numbers Y from this distribution.
6.12.2.2 Direct method
Generate n standard uniform random numbers U, and let Y be the number of these values of U that are less
than p.
6.12.2.3 Inverse function method
Compute the distribution function as follows:
y
n
⎛⎞kn−k
F(yp)=−(1p) ,y= 0,1, .,n
⎜⎟
∑
k
⎝⎠
k=0
Whenever a random number is required, generate a standard uniform random number U, and let Y be the
minimum y that satisfies U u F(y).
6.12.2.4 Alias method
First, n + 1 parameters v , v , ., v and n + 1 other parameters a , a , ., a , are calculated by the following
0 1 n 0 1 n
procedures.
a) v = (n + 1) p(y), y = 0, 1, …, n.
y
b) Let G be the set of suffices y that satisfies v W 1 and S be the set of suffices y that satisfies v < 1.
y y
c) While S is not empty, repeat the following operations from 1) to 4).
1) Select any element i from G and any element j from S.
2) Set a = i and v = v − (1 − v ).
j i i j
3) If v < 1, delete element i from G and move it to S.
i
4) Delete element j from S.
If the preparations mentioned above are completed, a binomial random number Y will be obtained by
performing the following operations d) to f).
d) Generate a standard uniform random number U, and set V = (n + 1)U.
e) Set k = (integer part of V) and u = V − k.
f) If u u v , set Y = k; otherwise, set Y = a .
k k
6.13 Poisson distribution
6.13.1 Probability mass function
The probability mass function of a Poisson distribution with mean µ is defined as follows.
y
µ
py()=−exp( µ) , 0y=, 1, 2, .,
y!
where µ > 0.
6.13.2 Method of using a relationship with an exponential distribution
Generate standard uniform random numbers U , U , ., and let Y be the maximum n that satisfies the
1 2
following inequality:
−−ln (1UU)(1− ).(1−U )< µ
{ }
12 n
6.13.3 Alias method
First, select a constant n for which the probability that Y > n is negligibly small, for example, the integer part of
µµ+ 6 could be specified to be n. Then use the procedure 6.12.2.4 alias method of the binomial distribution;
however, this time the probability mass function of the Poisson distribution shall be used for p(y).
NOTE This method is efficient when µ is of medium size (about 10 to 100).
6.14 Discrete uniform distribution
For generating discrete uniform random variates from M to N, a binary r bit random number sequence
generated by the method standardized in 5.1 is converted by the following procedures, where N − M + 1 is
r
assumed to be not greater than 2 .
a) Find the natural number k that satisfies the following inequality:
kk−1
21+−uuNM+12
NOTE 1 Such k is the minimum natural number that is equal to or greater than log (N − M + 1)
.
6 7
EXAMPLE 1 When N − M + 1 = 100, k = 7 because 2 + 1 = 65 u 100 u 2 = 128.
b) Add 1 to the binary integer that is constructed from the first k bits of a random number, and convert the
result to a decimal number.
k − 1 k − 2 k − 3 k − 4
NOTE 2 A k bit binary number Z Z Z Z .Z converts to a decimal number 2 Z + 2 Z + 2 Z + 2 Z
1 2 3 4 k 1 2 3 4
+ . + Z .
k
EXAMPLE 2 When the upper 7 bits are 1 011 001, the corresponding decimal number is 64 + 16 + 8 + 1 = 89,
and the decimal random number becomes 89.
c) The required decimal random number is the converted decimal number plus M − 1 by skipping numbers
greater than N.
r
NOTE 3 When N − M + 1 is more than 2 the required decimal random number can be obtained by concatenating
two or more binary random numbers to get one binary random number.
14 © ISO 2010 – All rights reserved

NOTE 4 When using the linear congruential method for generating pseudo-random numbers, k shall not be
specified to be equal to r.
Further, when N − M + 1 is a decimal k digit natural number, and k is not too large, say k is less than 20,
the method given in 5.2 can be used. The procedure is as follows.
d) Generate a decimal random number sequence of k digits by using procedure 5.2.
e) From the random number sequence which is generated by d) above, remove the numbers greater than N.
The sequence thus obtained is the required decimal random number sequence.
Annex A
(informative)
Table of physical random numbers
A.1 Table of random numbers
Physically generated random numbers have no functional relationship like pseudo-random numbers, and no
periodicity. Table A.1 shows a physically generated random number sequence obtained as measured values
of a property of a random physical system.
Table A.1 — Physical random number table
1 93 90 60 02 17 25 89 42 27 41 64 45 08 02 70 42 49 41 55 98
2 34 19 39 65 54 32 14 02 06 84 43 65 97 97 65 05 40 55 65 06
3 27 88 28 07 16 05 18 96 81 69 53 34 79 84 83 44 07 12 00 38
4 95 16 61 89 77 47 14 14 40 87 12 40 15 18 54 89 72 88 59 67
5 50 45 95 10 48 25 29 74 63 48 44 06 18 67 19 90 52 44 05 85
6 11 72 79 70 41 08 85 77 03 32 46 28 83 22 48 61 93 19 98 60
7 19 31 85 29 48 89 59 53 99 46 72 29 49 06 58 65 69 06 87 09
8 14 58 90 27 73 67 17 08 43 78 71 32 21 97 02 25 27 22 81 74
9 28 04 62 77 82 73 00 73 83 17 27 79 37 13 76 29 90 07 36 47
10 37 43 04 36 86 72 63 43 21 06 10 35 13 61 01 98 23 67 45 21
11 74 47 22 71 36 15 67 41 77 67 40 00 67 24 00 08 98 27 98 56
12 48 85 81 89 45 27 98 41 77 78 24 26 98 03 14 25 73 84 48 28
13 55 81 09 70 17 78 18 54 62 06 50 64 90 30 15 78 60 63 54 56
14 22 18 73 19 32 54 05 18 36 45 87 23 42 43 91 63 50 95 69 09
15 78 29 64 22 97 95 94 54 64 28 34 34 88 98 14 21 38 45 37 87
16 97 51 38 62 95 83 45 12 72 28 70 23 67 04 28 55 20 20 96 57
17 42 91 81 16 52 44 71 99 68 55 16 32 83 27 03 44 93 81 69 58
18 07 84 27 76 18 24 95 78 67 33 45 68 38 56 64 51 10 79 15 46
19 60 31 55 42 68 53 27 82 67 68 73 09 98 45 72 02 87 79 32 84
20 47 10 36 20 10 48 09 72 35 94 12 94 78 29 14 80 77 27 05 67
21 73 63 78 70 96 12 40 36 80 49 23 29 26 69 01 13 39 71 33 17
22 70 65 19 86 11 30 16 23 21 55 04 72 30 01 22 53 24 13 40 63
23 86 37 79 75 97 29 19 00 30 01 22 89 11 84 55 08 40 91 26 61
24 28 00 93 29 59 54 71 77 75 24 10 65 69 15 66 90 47 90 48 80
25 40 74 69 14 01 78 36 13 06 30 79 04 03 28 87 59 85 93 25 73

16 © ISO 2010 – All rights reserved

A.2 Method of physical random number generation
The method by which the physical random numbers in Table A.1 were generated is described below. The
source of the numbers is the noise generated by a diode. In a diode, the noise signal is large because, by
amplification using the electron avalanche effect, stabilized noise is easily obtained. For this reason, it is used
1)
very often as a noise source. For the element, NC2401 of Noisecom in the United States was used. This
element has a noise source and an amplifier built-in, and its band width is 1 GHz, while its amplitude is
160 mVrms.
The methods of digitalizing the noise signal are
a) analogue/digital conversion,
b) regarding noise as a pulse sequence, by counting pulses per unit time,
c) regarding noise as a pulse sequence, by measuring pulse interval.
2)
For example, consider the use of a DAS-4102 converter produced by the Keithley Instruments, Inc. for
analogue/digital conversion. This equipment has a resolution of 8 bits with a sampling period of 64 MHz
maximum. Data was sampled at 1 MHz to produce the attached table. Measuring was done with resolution
ability 3,91 mV/digit and only the lowest bit was used as a random number source.
Because the analogue/digital conversion equipment has errors that are peculiar to the equipment, the
frequency distribution of values after conversion is not uniform. To obtain a more uniform distribution, 2 bits
were generated from the same random number source, and
(0,1) → Random number (Rn) = 0
(1,0) → Random number (Rn) = 1
(0,0), (1,1) → abandoned
Random numbers in Table A.1 were generated according to the above scheme. If the probabilities of (0,1) and
(1,0) are equal to each other, the random number is uniformly distributed. Because the intervals between
successive measurements are as short as 1 ms, the characteristics of the equipment would scarcely change
in this time interval. Therefore, (0,1) and (1,0) are considered to conform to the same probability distribution.
An alternative method of correcting is by formerly measuring the probability distribution of the characteristics,
but, because this distribution varies from equipment to equipment, this method was not employed. Further, for
safety, 32 bits were gathered in one group, and exclusive or was done with pseudo-random numbers using
the Mersenne Twister (routine name genrand) which is described in 5.5. The Mersenne Twister was initialized
by the routine init_genrand(s), s = 19660809. If the original random number sequences are required, they can
be regenerated by operating exclusive or again using the Mersenne Twister.
Table A.1 is a decimal random number sequence generated by the above-mentioned method, taking the
upper 4 bits of the 32-bit random number sequence. If the value of this is not less than 0 and not more than 9,
the value is employed as the random number value; however, if the value of this is 10 or more, it is
abandoned and the next random number is generated.

1) NC2401 is the trade name of a product supplied by Noisecom. This information is given for the convenience of users
of this document and does not constitute an endorsement by ISO of the product named.
2) DAS-4102 is the trade name of a product supplied by Keithley Instruments, Inc. This information is given for the
convenience of users of this document and does not constitute an endorsement by ISO of the product named.
Annex B
(informative)
Algorithm for pseudo-random number generation
B.1 Program code for the trinomial GFSR method
The following C program given below, which is in accordance with ISO/IEC 9899, is an example with
1 279
parameters (p, q, w) = (1 279, 418, 32) and period (2 − 1). When the function gfsr( ) is called, it generates
an integer between 0 and (2 − 1) inclusive. When the function gfsr_31( ) is called, it generates an integer
between 0 and (2 − 1) inclusive. Before calling the functions gfsr( ) and gfsr_31( ), initialization is necessary
by calling init_gfsr(s) once.The initialization function init_gfsr(s) executes initialization under the condition that
an unsigned 32-bit integer [integer between 0 and (2 − 1)] is used as the seed. The generated sequence
can be used to provide 39 independent series, each of which has negligible auto-correlation, is 39-distributed
(uniformly distributed in a 39-dimensional hyper-cube) with 32-bit precision, and its auto-correlation function
1 274
assumes values close to zero up to phase shift 2 .
To obtain different pseudo-random number series, change the seed s given to init_gfsr(s). Only the constants
p, q, w in the program may be changed. The value of w will be a power of 2 within the word length of the
machine. The value of w will generally be 32 or 64, according to the machine. For example, if the word length
of the machine is 64, the constant w in the program is set to 64, and then gfsr( ) generates integers between 0
64 63
and (2 − 1) inclusive, while gfsr_31( ) generates integers between 0 and (2 − 1) inclusive.
In this program, the length of type “unsigned long” is presumed to be not less than 32 bits.
/*************************************************
C code : Trinomial GFSR
*************************************************/

#define P 1279
#define Q 418
#define W 32 /* W should be a power of 2 */

static unsigned long state [P] ;
static int state_i ;
void init_gfsr (unsigned long s)
{
int i, j, k;
static unsigned long x [P] ;
s &= 0xffffffffUL;
for (i=0 ; i

x [i] = s>>31 ;
18 © ISO 2010 – All rights reserved

s = 1664525UL * s + 1UL ;
s &= 0xffffffffUL ;
}
for (k=0, i=0 ; i

state [i] = 0UL ;
for (j=0 ; j state [i] <<= 1 ;
state [i] |= x [k]
...