3 L'Ecuyer's two-sequence pseudorandom generator with a 32 cell
4 Bays-Durham shuffle on the back end.
6 When updating the individual generators, an intermediate value as
7 large as 64 bits is required. Many modern C compilers provide a
8 "long long" 64-bit integer type. If this type is available,
9 compile this file with LONGLONG defined and it will be used,
10 generally improving execution speed. If the compiler does not
11 implement "long long", compile without LONGLONG and Schrage's
12 algorithm will be used to perform the computation using only
15 For additional details, see L'Ecuyer's paper:
17 L'Ecuyer, P. "Communications of the ACM", Vol. 31, p. 742 (1988).
19 and that of Bays and Durham:
21 Bays, Carter, and S.D. Durham. "ACM Transactions on
22 Mathematical Software", Vol. 2, p. 59 (1976).
24 Schrage's multiple-precision modular algorithm is described in:
26 Schrage, P. "ACM Transactions on Mathematical Software",
27 Vol. 5, p. 132 (1979).
37 #define int64 long long
43 /* L'Ecuyer's recommended multiplier and modulus for the two
44 multiplicative congruential generators. Even though the
45 values fit in 32 bits, we declare them as int64 so that the
46 arithmetic in calculating the next value will be automatically
47 done in int64 without need for casting when LONGLONG is defined.
48 In a non-LONGLONG build int64 is defined as long, so these values
49 participate correctly in the Schrage algorithm computation. */
51 #define mul1 ((int64) 40014)
52 #define mod1 ((int64) 2147483563)
53 #define mul2 ((int64) 40692)
54 #define mod2 ((int64) 2147483399)
56 #define shuffleSize 32 /* Shuffle table size */
57 #define warmup 19 /* Number of initial warmup results to "burn" */
59 static int32 gen1, gen2, state;
60 static int32 shuffle[shuffleSize]; /* Bays-Durham shuffle table */
62 /* Update generator which, using either long long arithmetic or
63 Schrage's algorithm depending on whether LONGLONG is defined.
64 The definition for Schrage's algorithm relies, for efficiency,
65 on the C compiler performing compile-time constant arithmetic
66 for the quotient and remainder of the modulus and multiplier.
67 If you're porting this to a language which lacks that feature,
68 you'll want to predefine the quotient and remainder for each
69 generator and use the explicit values. */
72 #define updgen(which) gen##which = (int32) ((gen##which * mul##which) % mod##which)
74 #define updgen(which) { \
75 int32 t = gen##which / (mod##which / mul##which); \
77 gen##which = mul##which * (gen##which - (t * (mod##which / mul##which))) - \
78 t * ((mod##which % mul##which)); \
79 if (gen##which < 0) { \
80 gen##which += mod##which; \
85 /* LEsetSeed -- Set seed for generator. Subsequent values will be based
86 on the given nonzero seed. */
88 void LEsetSeed(int32 seed)
94 gen1 = gen2 = (int32) (seed & 0x7FFFFFFFL);
96 /* "Warm up" the generator for a number of rounds to eliminate
97 any residual inflence of the seed. */
99 for (i = 0; i < warmup; i++) {
103 /* Fill the shuffle table with values. */
105 for (i = 0; i < shuffleSize; i++) {
107 shuffle[(shuffleSize - 1) - i] = gen1;
112 /* LEnextByte -- Get next byte from generator. */
114 unsigned char LEnextByte(void)
121 /* Extract shuffle table index from most significant part
122 of the previous result. */
124 i = state / (1 + (((int32) mod1) - 1) / shuffleSize);
126 /* New state is sum of generators modulo one of their moduli. */
128 state = (int32) (((shuffle[i]) + ((unsigned int32) gen2)) % mod1);
130 /* Replace value in shuffle table with generator 1 result. */
134 return (unsigned char) (state / (1 + (((int32) mod1) - 1) / 256));