From: "Andre Ratel"
To: "Earl F. Glynn"
Subject: Random number generators
Date: Tuesday, December 12, 2000 12:35 PM
Earl:
The 1988 algorithm:
==============
Last August, we exchanged some emails about the random
number generator found in Random.txt. We weren't sure
about the range of values returned by function Random.
Since this function is rather involved, I decided to use
brute force and to let the computer do all the work. The
idea was (using the notation of L'Ecuyer's 1988 paper)
- to try all possible values for the seed of s1, from 1 to
(m1 - 1) and to look at the range of values coming out
for the next s1
- to proceed similarly for variable s2, using all possible
values for the seed, from 1 to (m2 - 1).
This took almost two months of calculations (full time) on
my 486-66 for the s1 part and only 13 hours, or so, on the
Pentium III-733 of a friend for the s2 part. (No doubt: I
really need to upgrade my system!)
Here's what came out of this:
- for s1 seeds in the range 1..(m1 - 1), the loop
<<
k:= s div q1;
s1:= a1*(s - k*q1) - k*r1;
if s1 < 0 then
s1:= s1 + m1;
>>
outputs values in the range 1..(m1 - 1)
- for s2 seeds in the range 1..(m2 - 1), the loop
<<
k:= s div q2;
s2:= a2*(s - k*q2) - k*r2;
if s2 < 0 then
s2:= s2 + m2;
>>
outputs values in the range 1..(m2 - 1).
In other words, after using all possible seed values, the
minimum value I got for s1 is 1 and its maximum value
is (m1 - 1), and similarly for s2. The ranges of values for
s1 and s2 are then the same as those for their seeds.
Knowing this, it is easy to deduce the range of the other
variables.
Here's then is my implementation of L'Ecuyer's 1988 random
number generator
<<<
function CMLCG1988: extended;
{This function combines two multiplicative linear congruential number
generators to return a pseudo-random floating point value with a
uniform probability distribution on the open interval (0.0, 1.0).
It uses the algorithm found in L'Ecuyer [1988], p. 747, Fig. 3.
Rem: - As mentioned by L'Ecuyer [1988], p. 747 (right column), this
algorithm
"will work as long as the machine can represent
all integers in the range [-2^31 + 85, 2^31 - 85]."
The interval corresponds to [-2147483563, 2147483563] and it
shouldn't cause any problem since Delphi's integer type covers
the range 2147483648..2147483647.
- CMLCG stands for "combined multiplicative linear congruential
generators"
- We use here the same variable names as those found in L'Ecuyer's
paper.}
const
{Constants used to generate the first integer:}
m1 = 2147483563;
a1 = 40014;
q1 = 53668;
r1 = 12211;
{Constants used to generate the second integer:}
m2 = 2147483399;
a2 = 40692;
q2 = 52774;
r2 = 3791;
{Constant used to obtain a value between 0.0 and 1.0:}
h = 1.0/m1;
var
{Used for intermediary results:}
k, Z: Integer;
begin
{Generating the first integer value:}
k:= s1 div q1;
s1:= a1*(s1 - k*q1) - k*r1;
if s1 < 0 then
s1:= s1 + m1;
{rem: This gives an integer in the range
1 <= s1 <= (m1 - 1).}
{Generating the second integer value:}
k:= s2 div q2;
s2:= a2*(s2 - k*q2) - k*r2;
if s2 < 0 then
s2:= s2 + m2;
{rem: This gives an integer in the range
1 <= s2 <= (m2 - 1).}
{Combining the two results:}
Z:= s1 - s2;
if (Z < 1) then
Z:= Z + (m1 - 1);
{rem: Since
-(m2 - 2) <= (s1 - s2) <= (m1 - 2)
this yields an integer in the range
1 <= Z <= (m1 - 1).}
{Returning a floating point value:}
Result:= Z*h;
{rem: This yields
1/m1 <= Result <= (1.0 - 1/m1).}
end; {CMLCG1988}
>>>
The 1991 algorithm:
==============
The work described above was done more or less in vain since,
while reading
Pierre L'Ecuyer and Shue Tezuka,
"Structural properties for two classes of combined generators",
Mathematics of Computation, 57, (1991), 735--746.
I learned that, by changing the numerical values in CMLCG1988, it
is possible ot obtain, according to the authors, a LCG which has a
lattice structure of slightly better quality and with much more noise.
The new numerical values to be used are:
<<<
const
{Constants used to generate the first integer:}
m1 = 2147483647;
a1 = 26756;
q1 = 80261;
r1 = 20331;
{Constants used to generate the second integer:}
m2 = 2145483479;
a2 = 30318;
q2 = 70765;
r2 = 30209;
{Constant used to obtain a value between 0.0 and 1.0:}
h = 1.0/m1;
{Rem: As in L'Ecuyer [1988], p. 744, Eqs (15) and (16), the values
of q1, r1, q2, and r2 were obtained from
q1:= m1 div a1;
r1:= m1 - a1*q1;
q2:= m2 div a2;
r2:= m2 - a2*q2; }
>>>
The 1999 algorithms:
===============
The algorithms I found in
combmrg2.ps
P. L'Ecuyer,
"Good parameter sets for combined multiple recursive random number
generators",
Figs. I, II, III and Table VIII.
are, however, by far the best ones.
rem: This PostScript file is available at
http://www.iro.umontreal.ca/~lecuyer/papers.html
In the article, the code is written in C but it is really easy to translate it:
<<<
var
{Global variables used by random number generators MRG32k3a and MRG32k5a:}
s10, s11, s12, s13, s14, s20, s21, s22, s23, s24: Extended;
{Global variables used by random number generators MRG63k3a:}
is10, is11, is12, is20, is21, is22: Int64;
{rem: In L'Ecuyer Fig. III, these variables are named
s10, s11, s12, s20, s21, s22
Here, we changed their names so that functions MRG32k3a, MRG32k5a,
MRG63k3a could be part of the same unit.}
function MRG32k3a: Extended;
{This function combines 2 multiple recursive number generators of order 3
to return a pseudo-random number with a uniform probability distribution
in the open interval (0.0, 1.0).
rem: - This generator is well behaved in all dimensions up to at
least 45. Its period is
(m1^3 -1)*(m2^3 -1) = 2^191 = 3.139E57
(cf L'Ecuyer, Section 3, p. 11)
- This generator has been designed such that it can never return
0.0 or 1.0 exactly. This will allow one to generate exponential
random variables by using the logarithm of u or (1.0 - u), with
u random.
(cf L'Ecuyer, Section 3, p. 13)}
const
norm = 2.328306549295728e-10;
m1 = 4294967087.0;
m2 = 4294944443.0;
a12 = 1403580.0;
a13n = 810728.0;
a21 = 527612.0;
a23n = 1370589.0;
var
k: Integer;
p1, p2: Extended;
begin
{Component 1:}
p1:= a12 * s11 - a13n * s10;
k:= SciRound(p1/m1);
p1:= p1 - k * m1;
if (p1 < 0.0) then
p1:= p1 + m1;
s10:= s11;
s11:= s12;
s12:= p1;
{Component 2:}
p2:= a21 * s22 - a23n * s20;
k:= SciRound(p2/m2);
p2:= p2 - k * m2;
if (p2 < 0.0) then
p2:= p2 + m2;
s20:= s21;
s21:= s22;
s22:= p2;
{Combination:}
if (p1 <= p2) then
Result:= (p1 - p2 + m1) * norm
else
Result:= (p1 - p2) * norm;
end; {MRG32k3a}
function MRG32k5a: Extended;
{This function combines 2 multiple recursive number generators of order 5
to return a pseudo-random number with a uniform probability distribution
in the open interval (0.0, 1.0).
rem: The period of this generator is
(m1^5 -1)*(m2^5 -1) = 2^319 = 1.068E96
(cf L'Ecuyer, Section 3, p. 13).}
const
norm = 2.3283163396834613E-10;
m1 = 4294949027.0;
m2 = 4294934327.0;
a12 = 1154721.0;
a14 = 1739991.0;
a15n = 1108499.0;
a21 = 1776413.0;
a23 = 865203.0;
a25n = 1641052.0;
var
k: Integer;
p1, p2: Extended;
begin
{Component 1:}
p1:= a12 * s13 - a15n * s10;
if (p1 > 0.0) then
p1:= p1 - a14 * m1;
p1:= p1 + a14 * s11;
k:= SciRound(p1/m1);
p1:= p1 - k * m1;
if (p1 < 0.0) then
p1:= p1 + m1;
s10:= s11;
s11:= s12;
s12:= s13;
s13:= s14;
s14:= p1;
{Component 2:}
p2:= a21 * s24 - a25n * s20;
if (p2 > 0.0) then
p2:= p2 - a23 * m2;
p2:= p2 + a23 * s22;
k:= SciRound(p2/m2);
p2:= p2 - k * m2;
if (p2 < 0.0) then
p2:= p2 + m2;
s20:= s21;
s21:= s22;
s22:= s23;
s23:= s24;
s24:= p2;
{Combination:}
if (p1 <= p2) then
Result:= (p1 - p2 + m1) * norm
else
Result:= (p1 - p2) * norm;
end; {MRG32k5a}
function MRG63k3a: Extended;
{This function combines 2 multiple recursive number generators of order 3
to return a pseudo-random number with a uniform probability distribution
in the open interval (0.0, 1.0). Here, the generator has been inplemented
in 64-bit integer arithmetic.
rem: The period of this generator is
(m1^3 -1)*(m2^3 -1) = 2^377
(cf L'Ecuyer, Section 3, p. 14).}
const
norm = 1.0842021724855052e-19;
m1 = 9223372036854769163;
m2 = 9223372036854754679;
a12 = 1754669720;
q12 = 5256471877;
r12 = 251304723;
a13n = 3182104042;
q13 = 2898513661;
r13 = 394451401;
a21 = 31387477935;
q21 = 293855150;
r21 = 143639429;
a23n = 6199136374;
q23 = 1487847900;
r23 = 985240079;
var
h, p12, p13, p21, p23: Int64;
begin
{Component 1:}
h:= is10 div q13;
p13:= a13n * (is10 - h * q13) - h * r13;
h:= is11 div q12;
p12:= a12 * (is11 - h * q12) - h * r12;
if (p13 < 0) then
p13:= p13 + m1;
if (p12 < 0) then
p12:= p12 + m1 - p13
else
p12:= p12 - p13;
if (p12 < 0) then
p12:= p12 + m1;
is10:= is11;
is11:= is12;
is12:= p12;
{Component 2:}
h:= is20 div q23;
p23:= a23n * (is20 - h * q23) - h * r23;
h:= is22 div q21;
p21:= a21 * (is22 - h * q21) - h * r21;
if (p23 < 0) then
p23:= p23 + m2;
if (p21 < 0) then
p21:= p21 + m2 - p23
else
p21:= p21 - p23;
if (p21 < 0) then
p21:= p21 + m2;
is20:= is21;
is21:= is22;
is22:= p21;
{Combination:}
if (p12 > p21) then
Result:= (p12 - p21) * norm
else
Result:= (p12 - p21 + m1) * norm;
end; {MRG63k3a}
>>>
In this listing, the invoked SciRound is my rounding function:
<<<
function SciRound(x: extended): int64;
{This function rounds a floating point number x to its closest integer
value. It differs from Delphi's Round function for which (according to
the OnLine Help):
"If X is exactly halfway between two whole numbers, the result is
always the even number."
This is called "banker's rounding".
Our function implements scientific rounding: real numbers having 5 as
their first decimal are always rounded towards the integer with the
largest absolute value.
Examples:
x Round(x) SciRound(x)
--------------------------------
-13.6 -14 -14
-13.5 -14 -14
-13.4 -13 -13
-12.6 -13 -13
-12.5 -12 -13 <<<<<<<
-12.4 -12 -12
12.4 12 12
12.5 12 13 <<<<<<<
12.6 13 13
13.4 13 13
13.5 14 14
13.6 14 14
with the differences beween the results marked by <<<<<<< symbols.}
begin
if (Trunc(10*Frac(Abs(x))) >= 5) then
Result:= Trunc(x) + Sgn(x)*1
else
Result:= Trunc(x);
{rem: Trunc(10*Frac(Abs(x))) is the first decimal of x. We need to take
the Abs of x value since functions Frac and Truc keep the sign.}
end; {SciRound}
>>>
MRG32k3a, MRG32k5a, and MRG63k3a are most probably the random
number generators I'll use.
Sincerely.
Andre
=============================================================================
Date: Wed, 20 Dec 2000 06:20:39 -0500
To: "Earl F. Glynn"
From: Andre Ratel
Subject: More on random number generators
Earl:
In my previous email, I should also have mentioned the following
article
clcg4.ps
P. L'Ecuyer and T. H. Andres,
"A Random Number Generator Based on the Combination of Four LCGs'',
Mathematics and Computers in Simulation, 44 (1997), 99--107
It contains a little more than just a portable random number generator
and might be of some interest to people wanting to write more complex
simulation programs.
Here's the abstract:
<<<
A portable package for uniform random number generation is proposed,
based on a backbone generator with a period near 2^121 which is a
combination of four linear congruential generators. The package provides
for multiple (virtual) generators evolving in parallel. Each generator also
has many disjoint subsequences, and solftware tools are provided to reset
the state of any generator to the beginning of the first, previous or current
subsequences. Such facilities are helpful to maintain synchronization for
implementing variance reduction methods in simulation. Computer
implementations are available in the C and Modula-2 languages.
>>>
The Modula-2 files, also available on L'Ecuyer's Web page, are really nice
and read just like Delphi code.
Andre