§1  GB FLIP  INTRODUCTION 1

1. Introduction. This is GB FLIP, the module used by GraphBase programs to generate random numbers.

To use the routines in this file, first call the function gb init rand (seed ). Subsequent uses of the macro gb next rand ( ) will then return pseudo-random integers between 0 and 2^{31} − 1, inclusive.

GraphBase programs are designed to produce identical results on almost all existing computers and operating systems. An improved version of the portable subtractive method recommended in Seminumerical Algorithms, Section 3.6, is used to generate random numbers in the routines below. The period length of the generated numbers is at least 2^{55} − 1, and it is in fact plausibly conjectured to be 2^{85} − 2^{30} for all but at most one choice of the seed value. The low-order bits of the generated numbers are just as random as the high-order bits.

2. Changes might be needed when these routines are ported to different systems, because the programs have been written to be most efficient on binary computers that use two’s complement notation. Almost all modern computers are based on two’s complement arithmetic, but if you have a nonconformist machine you might have to revise the code in sections that are listed under ‘system dependencies’ in the index.

A validation program is provided so that installers can tell if GB FLIP is working properly. To make the test, simply run . ⟨ 2 ⟩ ≡ #include #include /∗ all users of GB FLIP should do this ∗/

int main( )

{ long j; gb init rand (−314159_{L}); if (gb next rand ( ) ≠ 119318998) {

fprintf (stderr , ␣ ␣ ␣ ␣ ); return −1;

}

for (j = 1; j ≤ 133; j^{++}) gb next rand ( ); if (gb unif rand (^{#}55555555_{L}) ≠ 748103812) {

fprintf (stderr , ␣ ␣ ␣ ␣ ); return −2;

}

fprintf (stderr , ␣ ␣ ␣ ␣ ␣ ␣ ); return 0;

}

3. The C code for GB FLIP doesn’t have a main routine; it’s just a bunch of subroutines to be incorporated into programs at a higher level via the system loading routine. Here is the general outline of :

⟨ Private declarations 4 ⟩

⟨ External declarations 5 ⟩

⟨ External functions 7 ⟩

2  THE SUBTRACTIVE METHOD  GB FLIP §4

4. The subtractive method. If m is any even number and if the numbers a_{0}, a_{1}, …, a_{54} are not all even, then the numbers generated by the recurrence

a_{n} = (a_{n−55} − a_{n−24}) mod m

have a period length of at least 2^{55} − 1, because the residues a_{n} mod 2 have a period of this length. Furthermore, the numbers 24 and 55 in this recurrence are sufficiently large that deficiencies in randomness due to the simplicity of the recurrence are negligible in most applications.

Here we take m = 2^{31} so that we get the full set of nonnegative numbers on a 32-bit computer. The recurrence is computed by maintaining an array of 55 values, A[1] … A[55]. We also set A[0] = −1 to act as a sentinel. ⟨ Private declarations 4 ⟩ ≡

static long A[56] = {−1}; /∗ pseudo-random values ∗/

This code is used in section 3.

5. Every external variable should be declared twice in this file: once for GB FLIP itself (the “real” declaration for storage allocation purposes), and once in (for cross-references by GB FLIP users).

The pointer variable gb fptr should not be mentioned explicitly by user routines. It is made public only for efficiency, so that the gb next rand macro can access the private A table. ⟨ External declarations 5 ⟩ ≡

long ∗gb fptr = A; /∗ the next A value to be exported ∗/

This code is used in section 3.

6. The numbers generated by gb next rand ( ) seem to be satisfactory for most purposes, but they do fail a stringent test called the “birthday spacings test,” devised by George Marsaglia. [See, for example, Statistics and Probability Letters 8 (1990), 35–39.] One way to get numbers that pass the birthday test is to discard half of the values, for example by changing ‘gb flip cycle( )’ to ‘(gb flip cycle( ), gb flip cycle( ))’ in the definition of gb next rand ( ). Users who wish to make such a change should define their own substitute macro.

Incidentally, we hope that optimizing compilers are smart enough to do the right thing with gb next rand . #define gb next rand () (∗gb fptr ≥ 0 ? ∗gb fptr^{−} : gb flip cycle( )) ⟨ 6 ⟩ ≡ #define gb next rand ( ) (∗gb fptr ≥ 0 ? ∗gb fptr^{−} : gb flip cycle( ))

extern long ∗gb fptr ; /∗ the next A value to be used ∗/

extern long gb flip cycle( ); /∗ compute 55 more pseudo-random numbers ∗/

See also sections 11 and 13.

§7  GB FLIP  THE SUBTRACTIVE METHOD 3

7. The user is not supposed to call gb flip cycle directly either. It is a routine invoked by the macro gb next rand ( ) when gb fptr points to the negative value in A[0].

The purpose of gb flip cycle is to do 55 more steps of the basic recurrence, at high speed, and to reset gb fptr .

The nonnegative remainder of (x − y) mod 2^{31} is computed here by doing a logical-and with the constant #7fffffff. This technique doesn’t work on computers that do not perform two’s complement arithmetic. An alternative for such machines is to add the value 2^{30} twice to (x − y), when (x − y) turns out to be negative. Careful calculations are essential because the GraphBase results must be identical on all computer systems.

The sequence of random numbers returned by successive calls of gb next rand ( ) isn’t really a_{n}, a_{n+1}, …, as defined by the basic recurrence above. Blocks of 55 consecutive values are essentially being “flipped” or “reflected”—output in reverse order—because gb next rand ( ) makes the value of gb fptr decrease instead of increase. But such flips don’t make the results any less random. #define mod diff (x, y) (((x) − (y)) &^{#}7fffffff) /∗ difference modulo 2^{31} ∗/ ⟨ External functions 7 ⟩ ≡

long gb flip cycle( )

{ register long ∗ii , ∗jj ; for (ii = &A[1], jj = &A[32]; jj ≤ &A[55]; ii^{++}, jj^{++}) ∗ii = mod diff (∗ii , ∗jj ); for (jj = &A[1]; ii ≤ &A[55]; ii^{++}, jj^{++}) ∗ii = mod diff (∗ii , ∗jj ); gb fptr = &A[54]; return A[55];

}

See also sections 8 and 12. This code is used in section 3.

4  INITIALIZATION  GB FLIP §8

8. Initialization. To get everything going, we use a scheme like that recommended in Seminumerical Algorithms, but revised so that the least significant bits of the starting values depend on the entire seed, not just on the seed’s least significant bits.

Notice that we jump around in the array by increments of 21, a number that is relatively prime to 55. Repeated skipping by steps of 21 mod 55 keeps the values we’re computing spread out as far from each other as possible in the array, since 21, 34, and 55 are consecutive Fibonacci numbers (see the discussion of Fibonacci hashing in Section 6.4 of Sorting and Searching). Our initialization mechanism would be rather poor if we didn’t do something like that to disperse the values (see Seminumerical Algorithms, exercise 3.2.2–2). ⟨ External functions 7 ⟩ +≡

void gb init rand (seed )

long seed ; { register long i;

register long prev = seed , next = 1;

seed = prev = mod diff (prev , 0); /∗ strip off the sign ∗/ A[55] = prev ; for (i = 21; i; i = (i + 21) % 55) {

A[i] = next ;

⟨ Compute a new next value, based on next , prev , and seed 9 ⟩; prev = A[i];

}

⟨ Get the array values “warmed up” 10 ⟩; }

9. Incidentally, if fails, the person debugging these routines will want to know some of the intermediate numbers computed during initialization. The first nontrivial values calculated by gb init rand are A[42] = 2147326568, A[8] = 1073977445, and A[29] = 536517481. Once you get those right, the rest should be easy.

An early version of this routine simply said ‘seed ≫ 1’ instead of making seed shift cyclically. This method had an interesting flaw: When the original seed was a number of the form 4s + 1, the first 54 elements A[1], …, A[54] were set to exactly the same values as when seed was 4s + 2. Therefore one out of every four seed values was effectively being wasted. ⟨ Compute a new next value, based on next , prev , and seed 9 ⟩ ≡

next = mod diff (prev , next );

if (seed & 1) seed = ^{#}40000000 + (seed ≫ 1);

else seed ≫= 1; /∗ cyclic shift right 1 ∗/

next = mod diff (next , seed );

This code is used in section 8.

§10  GB FLIP  INITIALIZATION 5

10. After the first 55 values have been computed as a function of seed , they aren’t random enough for us to start using them right away. For example, we have set A[21] = 1 in order to ensure that at least one starting value is an odd number. But once the sequence a_{n} gets going far enough from its roots, the initial transients become imperceptible. Therefore we call gb flip cycle five times, effectively skipping past the first 275 elements of the sequence; this has the desired effect. It also initializes gb fptr .

Note: It is possible to express the least significant bit of the generated numbers as a linear combination mod 2 of the 31 bits of seed and of the constant 1. For example, the first generated number turns out to be odd if and only if

s_{24} + s_{23} + s_{22} + s_{21} + s_{19} + s_{18} + s_{15} + s_{14} + s_{13} + s_{11} + s_{10} + s_{8} + s_{7} + s_{6} + s_{2} + s_{1} + s_{0}

is odd, when seed = (s_{31} … s_{1}s_{0})_{2}. We can represent this linear combination conveniently by the hexadecimal number ^{#}01ecedc7; the stands for s_{24} and the final stands for s_{2} + s_{1} + s_{0}. The first ten least-significant

bits turn out to be respectively ^{#}01ecedc7, ^{#}dbbdc362, ^{#}400e0b06, ^{#}0eb73780, ^{#}da0d66ae, ^{#}002b63bc,

#adb801ed, ^{#}8077bbbc, ^{#}803d9db5, and ^{#}401a0eda in this notation (using the sign bit to indicate cases when 1 must be added to the sum).

We must admit that these ten 32-bit patterns do not look at all random; the number of ’s, ’s, and ’s is unusually high. (Before the “warmup cycles,” the patterns are even more regular.) This phenomenon eventually disappears, however, as the sequence proceeds; and it does not seem to imply any serious deficiency in practice, even at the beginning of the sequence, once we’ve done the warmup exercises. ⟨ Get the array values “warmed up” 10 ⟩ ≡

(void) gb flip cycle( );

(void) gb flip cycle( );

(void) gb flip cycle( );

(void) gb flip cycle( );

(void) gb flip cycle( );

This code is used in section 8.

11. ⟨ 6 ⟩ +≡ extern void gb init rand ( );

6  UNIFORM INTEGERS  GB FLIP §12

12. Uniform integers. Here is a simple routine that produces a uniform integer between 0 and m − 1, inclusive, when m is any positive integer less than 2^{31}. It avoids the bias toward small values that would occur if we simply calculated gb next rand ( ) % m. (The bias is insignificant when m is small, but it can be serious when m is large. For example, if m ≈ 2^{32}/3, the simple remainder algorithm would give an answer less than m/2 about 2/3 of the time.)

This routine consumes fewer than two random numbers, on the average, for any fixed m.

In the program (main), this routine should compute t = m, then it should reject the values r = 2081307921, 1621414801, and 1469108743 before returning the answer 748103812. #define two to the 31 ((unsigned long) ^{#}80000000) ⟨ External functions 7 ⟩ +≡

long gb unif rand (m)

long m; { register unsigned long t = two to the 31 − (two to the 31 % m);

register long r;

do {

r = gb next rand ( ); } while (t ≤ (unsigned long) r); return r % m;

}

13. ⟨ 6 ⟩ +≡ extern long gb unif rand ( );

§14  GB FLIP  INDEX 7

14. Index. Here is a list that shows where the identifiers of this program are defined and used.

A: \underline{4}. fprintf : 2. gb flip cycle: \underline{6}, \underline{7}, 10. gb fptr : \underline{5}, \underline{6}, 7, 10. gb init rand : 1, 2, \underline{8}, 9, \underline{11}. gb next rand : 1, 2, 5, \underline{6}, 7, 12. gb unif rand : 2, \underline{12}, \underline{13}. i: \underline{8}. ii : \underline{7}. j: \underline{2}. jj : \underline{7}. m: \underline{12}. main: \underline{2}, 12. mod diff : \underline{7}, 8, 9. next : \underline{8}, 9. prev : \underline{8}, 9. r: \underline{12}. seed : 1, \underline{8}, 9, 10. stderr : 2. system dependencies: 7. t: \underline{12}. two to the 31 : \underline{12}.

⟨ Compute a new next value, based on next , prev , and seed 9 ⟩ Used in section 8. ⟨ External declarations 5 ⟩ Used in section 3.

⟨ External functions 7, 8, 12 ⟩ Used in section 3.

⟨ Get the array values “warmed up” 10 ⟩ Used in section 8. ⟨ Private declarations 4 ⟩ Used in section 3.

⟨ 6, 11, 13 ⟩

⟨ 2 ⟩

May 18, 1999 at 02:58

GB FLIP

Section Page

Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1 The subtractive method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 2 Initialization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 4 Uniform integers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 6 Index . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 7

\begin{matrix}c \\ ◯ \end{matrix} 1993 Stanford University

This file may be freely copied and distributed, provided that no changes whatsoever are made. All users are asked to help keep the Stanford GraphBase files consistent and “uncorrupted,” identical everywhere in the world. Changes are permissible only if the modified file is given a new name, different from the names of existing files in the Stanford GraphBase, and only if the modified file is clearly identified as not being part of that GraphBase. (The CWEB system has a “change file” facility by which users can easily make minor alterations without modifying the master source files in any way. Everybody is supposed to use change files instead of changing the files.) The author has tried his best to produce correct and useful programs, in order to help promote computer science research, but no warranty of any kind should be assumed. Preliminary work on the Stanford GraphBase project was supported in part by National Science Foundation grant CCR-86-10181.