## Minimum Standard Random Number Generator

### January 14, 2014

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 xn, the next random integer in a random sequence is given by computing xn+1 = a xn (mod m), where a = 75 = 16807 and m = 231 − 1 = 2147483647; as a check on your implementation, if x0 = 1, then x10000 = 1043618065. Park and Miller chose m as the largest Mersenne prime less than 232; the smallest primitive root of m is 7, and since 5 is also prime, 75 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 · lor · 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

### 11 Responses to “Minimum Standard Random Number Generator”

1. Paul said

In Python. As in Python integer overflow cannot happen, I only implemented the simple version.

```from __future__ import division

import time

def minst(seed=None):
M = 2147483647
if seed is None:
seed = (int(time.clock() * 1e15)) % M
while 1:
seed = (16807 * seed) % M
yield seed / M
```
2. Felipe Ceotto said

I’m having trouble with the results of the suggested solution: when I run minstd0 three times I get results 16807, 282475249 and 1622650073 but when I do manual calculations, the third result for me is: 1622652946. And I can verify that by getting (M * 221) + 1622652946 = (4745938856997 + 1622652946) = 4747561509943 = (16807 * 282475249).

So, which one is correct? Am I missing something in the calculations? Obviously, due to this difference I can’t get the same X(10001) result as suggested in the post.

Thanks!

3. programmingpraxis said

1622650073 is correct. See A096550.

4. Daniel said

I believe that if x0 = 1, then x10000 is 1043618065. It currently says that x10001 is 1043618065.

Here’s a solution using Julia. I’m not too confident that I did nextInt3 as intended, but it got the right answer in my tests.

```function nextInt1(x, a = 16807, m = 2147483647)
(a * x) % m;
end

function nextInt2(x, a = 16807, m = 2147483647)
q = div(m, a);
r = m % a;
lo = x % q;
hi = div(x, q);
x = a*lo - r*hi;
if (x <= 0)
x += m;
end
x;
end

function nextInt3(x, a = 16807, m = 2147483647)
ax = a*x;
q = ax & (1 << 31 - 1);
p = ax >> 31;
s = p + q;
if (s >= m)
s -= m;
end
s;
end

```
5. treeowl said

Don’t you mean “ensuring that r<a” rather than “r<q”? The former is the result of the usual division algorithm; the latter, it seems, may force r to be negative.

6. treeowl said

I just checked the paper, and it does seem that r<a is intended. More significantly, you swapped the definitions of lo and hi, so the math doesn’t work out right.

7. treeowl said

It’s not terribly pretty, but here it is. Prettified on codepad

```module Main where
import           Data.Bits
import           Data.Int

a32=7^5::Int32

a64::Int64
a64=7^5

m32=2^31-1::Int32

m64::Int64
m64= 2^31-1

(q32,r32)=m32 `divMod` a32

getRandom1:: Int64 -> Int64
getRandom1 xold = (a64*xold) `mod` m64

low31=0x7FFFFFFF::Int64

getRandom3 :: Int64 -> Int64
getRandom3 xold = if newmodded Int32

getRandom2 xold = if res >0 then res else res +m32
where
(hi,lo)=xold `divMod` q32
res = a32*lo-r32*hi

randoms1 = iterate getRandom1

randoms2 = iterate getRandom2

randoms3 = iterate getRandom3

numtest = 1000::Int

main :: IO ()
main = do
let ones = take numtest . randoms1
let twos = take numtest . map (fromInteger.fromIntegral) . randoms2
let threes = take numtest . randoms3

print (ones 1 == threes 1)
print (twos 1==threes 1)
```
8. treeowl said

Since that was such a mess, I’ll post a cleaned up version. I’m hoping PP will delete the first one. To actually compile this, you will have to split it into two files with appropriate names.

```{-# START_FILE Lincon.hs #-}
module Lincon (getRandom1,getRandom2,getRandom3,randoms1,randoms2,randoms3) where
import           Data.Bits
import           Data.Int

a::Integral x => x
a=7^5
m::Integral x => x
m=2^31-1

a32=a::Int32
a64=a::Int64
m32=m::Int32
m64=m::Int64

getRandom1:: Int64 -> Int64
getRandom1 xold = (a64*xold) `mod` m64

getRandom3 :: Int64 -> Int64
getRandom3 xold = if newmodded Int32
getRandom2 xold = if res >0 then res else res +m32
where
(q,r)=m `divMod` a
(hi,lo)=xold `divMod` q
res = a32*lo-r*hi

randoms1 = iterate getRandom1
randoms2 = iterate getRandom2
randoms3 = iterate getRandom3

{-# START_FILE Main.hs #-}
module Main (main) where
import Lincon

numtest = 100000::Int

main :: IO ()
main = do
let ones = take numtest . randoms1
let twos = take numtest . map fromIntegral . randoms2
let threes = take numtest . randoms3
print (ones 1 == threes 1)
print (twos 1==threes 1)
```
9. […] working on solutions to the latest Programming Praxis puzzles—the first on implementing the minimal standard random number generator and the second on implementing a […]

10. […] statement Read the problem statement here at Programming Praxis. Be warned that the blog post mixes up the values for lo and hi. It took me […]

11. […] built several random number generators: , , , , , , , ,  (I didn’t realize it was so many until I went back and looked). In today’s […]