Minimum Standard Random Number Generator

January 14, 2014

Here is a basic version of the minimum standard random number generator, which updates a global variable seed and returns the new seed.

(define (minstd0) (set! seed (modulo (* 16807 seed) 2147483647)) seed)

A better design hides the seed inside a closure and returns a random fraction between zero and one, exclusive. It is also better to set the seed randomly, then allow the user to reset the seed if it is desired to repeat the same random sequence. A common way to set the seed uses the system clock:

(define minstd1
  (let ((a 16807) (m 2147483647) (seed (time-second (current-time))))
    (lambda args (if (pair? args) (set! seed (modulo (car args) m)))
      (set! seed (modulo (* a seed) m)) (/ seed m))))

The expression (time-second (current-time)) gets the number of seconds since the Unix epoch; it works in Chez Scheme, but other Scheme systems have other conventions for getting the time. When called as (minstd1 seed), this resets the seed before returning a random fraction.

Using Schrage multiplication makes the function only a little bit more complicated:

(define minstd2
  (let ((a 16807) (m 2147483647) (q 127773) (r 2836)
        (seed (time-second (current-time))))
    (lambda args (if (pair? args) (set! seed (modulo (car args) m)))
      (let* ((lo (modulo seed q)) (hi (quotient seed q))
             (test (- (* a lo) (* r hi))))
        (set! seed (if (positive? test) test (+ test m)))
        (/ seed m)))))

Although the point of Schrage multiplication is to prevent intermediate overflow, this doesn’t prevent intermediate overflow in Scheme. Most Scheme systems use the two least-significant bits of any value as tag bits, so the largest integer that can be encoded natively is 230, which is too small. This will work in other languages, however.

Implementation of Carta’s version of the function uses the bit operations from the Standard Prelude. An optimized version of this function is available at http://www.firstpr.com.au/dsp/rand31/.

(define minstd3
  (let ((a 16807) (m 2147483647) (seed (time-second (current-time))))
    (lambda args (if (pair? args) (set! seed (modulo (car args) m)))
      (let* ((lo (* a (logand seed #xFFFF)))
             (hi (* a (ash seed -16)))
             (lo (+ lo (ash (logand hi #x7FFF) 16)))
             (lo (+ lo (ash hi -15)))
             (lo (if (>= lo #x7FFFFFFF) (- lo #x7FFFFFFF) lo)))
        (set! seed lo) (/ seed m)))))

This actually looks better in C, assuming seed is a global unsigned long integer:

long unsigned int minstd()
{
    long unsigned int lo, hi;

    lo = 16807 * (seed & 0xFFFF);
    hi = 16807 * (seed >> 16);

    lo += (hi & 0x7FFF) <> 15;

    if (lo >= 0x7FFFFFFF)
        lo -= 0x7FFFFFFF;

    return ( seed = (long) lo );
}

We can test the four implementations by computing the 10001st value starting from 1; we stop at 9999 because one iteration is consumed by the initial setting of the seed and one iteration is performed in the final call after the loop completes:

> (define seed 1)
> (minstd0)
16807
> (do ((n 1 (+ n 1))) ((= 9999 n) (minstd0)) (minstd0))
1043618065
> (minstd1 1)
16807/2147483647
> (do ((n 1 (+ n 1))) ((= 9999 n) (minstd1)) (minstd1))
1043618065/2147483647
> (minstd2 1)
16807/2147483647
> (do ((n 1 (+ n 1))) ((= 9999 n) (minstd2)) (minstd2))
1043618065/2147483647
> (minstd3 1)
16807/2147483647
> (do ((n 1 (+ n 1))) ((= 9999 n) (minstd3)) (minstd3))
1043618065/2147483647

You can run the program at http://programmingpraxis.codepad.org/vqL1b0fL.

About these ads

Pages: 1 2

10 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 […]

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

Follow

Get every new post delivered to your Inbox.

Join 609 other followers

%d bloggers like this: