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.
Programming Praxis 
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()[…] built several random number generators: [1], [2], [3], [4], [5], [6], [7], [8], [9] (I didn’t realize it was so many until I went back and looked). In today’s […]