Shuffle Box

January 17, 2014

We begin with Knuth’s random number generator. Note that Scheme integers have no fixed size, so the modulo operation must be made explicit. Here is Knuth’s 32-bit random number generator following the same conventions as minstd from the previous exercise:

(define knuth
  (let* ((a 69069) (c 1234567) (m 4294967296)
         (seed (time-second (current-time))))
    (lambda args
      (when (pair? args) (set! seed (modulo (car args) m)))
      (set! seed (modulo (+ (* a seed) c) m))
      (/ seed m))))

Adding a shuffle box isn’t hard. The array of length k is stored in variable box, local function init takes a seed and returns a box, and local function next updates the random sequence. Seed and an element of box are always updated at the same time:

(define rand ; knuth random number generator with shuffle box
  (let* ((a 69069) (c 1234567) (m 4294967296) (k 32) ; 32-bit
         ; (a 6364136223846793005) (c 1442695040888963407)
         ; (m 18446744073709551616) (k 256) ; 64-bit
         (seed (time-second (current-time)))
         (next (lambda ()
           (set! seed (modulo (+ (* a seed) c) m)) seed))
         (init (lambda (seed) (let ((box (make-vector k)))
           (do ((j 0 (+ j 1))) ((= j k) box)
             (vector-set! box j (next))))))
         (box (init seed)))
    (lambda args
      (when (pair? args)
        (set! seed (modulo (car args) m)) (set! box (init seed)))
      (let* ((j (quotient (* k seed) m)) (n (vector-ref box j)))
        (set! seed (next)) (vector-set! box j seed) (/ n m)))))

Adding a shuffle box to the Park-Miller random number generator requires only a change to the next function; everything else stays the same. The function admits Schrage multiplication if your language worries about overflow, or you can use Carta’s algorithm, though the code below does neither:

(define minstd ; minimum standard rng with shuffle box
  (let* ((a 16807) (m 2147483647) (k 32)
         (seed (time-second (current-time)))
         (next (lambda ()
           (set! seed (modulo (* a seed) m)) seed))
         (init (lambda (seed) (let ((box (make-vector k)))
           (do ((j 0 (+ j 1))) ((= j k) box)
             (vector-set! box j (next))))))
         (box (init seed)))
    (lambda args
      (when (pair? args)
        (set! seed (modulo (car args) m)) (set! box (init seed)))
      (let* ((j (quotient (* k seed) m)) (n (vector-ref box j)))
        (set! seed (next)) (vector-set! box j seed) (/ n m)))))

Both these functions are fine random number generators with a reasonable period and good spectral properties, simple and straight forward to implement and use, returning random fractions from zero to one, suitable for any non-cryptographic purpose, including simulation. Of the two, I have a small preference for Knuth’s, because it regularizes zero, has a somewhat larger period, and a 64-bit version is available.

Often, instead of a random fraction, what is needed is a random integer on a specific range. A random fraction can be converted to an integer on the range lo inclusive to hi exclusive like this:

(define (randint . args)
  (let ((lo (if (pair? (cdr args)) (car args) 0))
        (hi (if (pair? (cdr args)) (cadr args) (car args))))
    (+ lo (floor (* (rand) (- hi lo))))))

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

About these ads

Pages: 1 2

One Response to “Shuffle Box”

  1. Paul said

    In Python. I.s.o. creating dedicated shuffle boxes for Park Miller and Knuth prng, I created a shuffle box, that takes any prng.

    from __future__ import division
    import time
    T0 = time. clock()
    
    def get_seed(seed=None):
        if seed is None:
            return (int(time.clock() * 2147483647 * 2147483647))
        return seed
        
    def park_miller(seed=None):
        """  yields floats in range (0, 1)
        """
        M = 2147483647
        A = 16807
        seed = get_seed(seed)
        while 1:
            seed = (A * seed) % M
            yield seed / M
            
    def knuth(seed=None, bits=32):
        """   yields floats in range [0, 1)"""
        a, c, m = ((69069, 1234567, 2 ** 32) if bits == 32 else
            (6364136223846793005, 1442695040888963407, 2 ** 64))
        seed  = get_seed(seed)
        while 1:
            seed = (seed * a + c) % m
            yield seed / m
    
    def prng_shuffle(prng, nboxes=55):
        """generator that yields a uniform random number in range[0,1)
            prng: function that returns a uniform random number in range[0,1)
        """
        boxes = [prng() for _ in range(nboxes)]
        index = int(prng() * nboxes)
        while 1:
            result = boxes[index]
            boxes[index] = prng()
            index = int(result * nboxes)
            yield result
            
    #example
    rand = prng_shuffle(knuth(bits=32).next).next
    print rand()
    

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 627 other followers

%d bloggers like this: