Programming Praxis

In 1988, Stephen K. Park and Keith W. Miller, reacting to the plethora of unsatisfactory random number generators then available, published a linear congruential random number generator that they claimed should be the “minimum standard” for an acceptable random number generator. Twenty-five years later, the situation is only somewhat better, and embarrassingly bad random number generators are still more common than they ought to be. Today’s exercise implements the original minimum standard random number generator in several different forms.

We begin with the original minimum standard random number generator. Given a random integer x_n, the next random integer in a random sequence is given by computing x_n+1 = a x_n (mod m), where a = 7⁵ = 16807 and m = 2³¹ − 1 = 2147483647; as a check on your implementation, if x₀ = 1, then x₁₀₀₀₀ = 1043618065. Park and Miller chose m as the largest Mersenne prime less than 2³²; the smallest primitive root of m is 7, and since 5 is also prime, 7⁵ is also a primitive root, hence their choice of a. Because a is a primitive root of m, all values in the range 1 to m − 1 inclusive will be generated before any repeat, so the random number generator has full period. The multiplier a = 16807 has been shown to have good randomness properties. Subsequent to their original paper, Park and Miller recommended 48271 as an improvement, and some people use 69621, but we’ll continue to use 16807.

The easiest way to implement that is obvious: just multiply a by the current value of x and compute the modulus. But that may cause overflow in the intermediate multiplication, rendering the results incorrect. A trick of Linus Schrage allows that multiplication to be done without overflow: Compute q = ⌊m / a⌋ and r = m (mod a) so that m = a q + r. Then a new x can be computed by hi = ⌊x / q⌋, lo = x (mod q), x = a · lo − r · hi, then adding m to x if x ≤ 0. In the case a = 16807 and m = 2147483647, q = 127773 and r = 2836; if you use a = 48271, then q = 44488 and r = 3399.

That works properly, without overflow, but each step in the sequence requires two divisions, two multiplications, a subtraction, a comparison, and possibly an addition, which can be slow. David Carta found a clever way to make that computation without any divisions: The product a · x is typically a 46-bit number. If we set q to the 31 low-order bits and p to the 15 high-order bits, then the product is q + p · 0x80000000, which is equivalent to q + p · 0x7FFFFFFF + p. But since we are working mod 0x7FFFFFFF, we can ignore the middle term and just take q + p. That 32-bit value might be larger tham m, in which case we can just subtract m to get it in range.

Your task is to implement all three versions of the minimum standard random number generator. When you are finished, you are welcome to read or run a suggested solution, or to post your own solution or discuss the exercise in the comments below.

Pages: 1 2