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