GB_FLIP
May 25, 2010
Here is our version:
(define two31 #x80000000)
(define a (make-vector 56 -1))
(define fptr #f)
(define (mod-diff x y)
(modulo (- x y) two31))
; use this if you have logand
; (define (mod-diff x y)
; (logand (- x y) #x7FFFFFFF))
(define (flip-cycle)
(do ((ii 1 (+ ii 1)) (jj 32 (+ jj 1))) ((< 55 jj))
(vector-set! a ii (mod-diff (vector-ref a ii) (vector-ref a jj))))
(do ((ii 25 (+ ii 1)) (jj 1 (+ jj 1))) ((< 55 ii))
(vector-set! a ii (mod-diff (vector-ref a ii) (vector-ref a jj))))
(set! fptr 54)
(vector-ref a 55))
(define (init-rand seed)
(let* ((seed (mod-diff seed 0)) (prev seed) (next 1))
(vector-set! a 55 prev)
(do ((i 21 (modulo (+ i 21) 55))) ((zero? i))
(vector-set! a i next)
(set! next (mod-diff prev next))
(set! seed (+ (quotient seed 2) (if (odd? seed) #x40000000 0)))
(set! next (mod-diff next seed))
(set! prev (vector-ref a i)))
(flip-cycle) (flip-cycle) (flip-cycle) (flip-cycle) (flip-cycle)))
(define (next-rand)
(if (negative? (vector-ref a fptr)) (flip-cycle)
(let ((next (vector-ref a fptr)))
(set! fptr (- fptr 1)) next)))
(define (unif-rand m)
(let ((t (- two31 (modulo two31 m))))
(let loop ((r (next-rand)))
(if (<= t r) (loop (next-rand)) (modulo r m)))))
Standard Scheme doesn’t provide a logand function, so we gave an alternative to Knuth’s mod-diff macro. You should use the logand version if possible, since bit operations are faster than division.
Knuth provides a test function:
(define (test-flip)
(init-rand -314159)
(when (not (= (next-rand) 119318998))
(error 'test-flip "Failure on the first try!"))
(do ((j 1 (+ j 1))) ((< 133 j)) (next-rand))
(when (not (= (unif-rand #x55555555) 748103812))
(error 'test-flip "Failure on the second try!"))
(display "OK, the gb-flip routines seem to work!") (newline))
You can run the program at http://programmingpraxis.codepad.org/nBm7d0wF. We’ll install GB_FLIP in our Standard Prelude, replacing the current random number generator, but using our own interface.
Programming Praxis 
I’ve been trying to come up with a Python version, but the global variables are giving me the run around.
Here’s my Python version. It’s a fairly direct translation of Knuth’s C-code to Python.
from itertools import count, takewhile A = [-1L]*56 fptr = 0 def next_rand(): global A, fptr if A[fptr] >= 0: val = A[fptr] fptr -= 1 else: val = flip_cycle() return val def mod_diff(x, y): return (x - y) & 0x7fffffffL def flip_cycle(): global A, fptr for ii in range(1,56): A[ii] = mod_diff( A[ii], A[ (ii + 30)%55 + 1 ] ) fptr = 54 return A[ 55 ] def init_rand( seed ): global A prev = mod_diff( long( seed ), 0 ) seed = prev A[55] = prev next_ = 1 for i in takewhile( bool, ( 21*n%55 for n in count( 1 ) ) ): A[i] = next_ next_ = mod_diff( prev, next_ ) if seed&1: seed = 0x40000000L | (seed>>1) else: seed >>= 1 next_ = mod_diff( next_, seed ) prev = A[i] for j in range(5): flip_cycle() def unif_rand( m ): r = t = 0x80000000L - (0x80000000L % m) while r >= t: r = next_rand() return r%m def test(): init_rand( -314159 ) if next_rand() != 119318998L: print "Failure on first try!" else: for j in range(133): _ = next_rand() if unif_rand( 0x55555555L ) != 748103812L: print "Failure on second test!" else: print "Ok, the routines seem to work!"[…] random number generators in past exercises, including Donald Knuth’s lagged-fibonacci generator that is used in the Standard Prelude. We also looked at cellular automata in a previous exercise. […]
Java implementation.
package edu.mit.yingyin.util; import java.util.Arrays; /** * Random number generator based on Knuth's GB_FLIP, the module used by * GraphBase programs to generate random numbers. */ public class Random { /** * The index of the next pseudo-random value in randArray to be used. */ private int index = 0; /** * Pseudo-random values. */ private long[] randArray = new long[56]; public Random() { Arrays.fill(randArray, -1); } /** * @return the next random number. */ public long next() { long val; if (randArray[index] >= 0) { val = randArray[index]; index--; } else { val = flipCycle(); } return val; } /** * Returns a uniform integer between 0 and n - 1, inclusive. * @param n a positive integer less than 2^31. * @return */ public long next(long n) { long r = 0x80000000L - (0x80000000L % n); long t = r; while (r >= t) r = next(); return r % n; } /** * Initialize the random number generator with a seed. * @param seed */ public void setSeed(long seed) { // Strips off the sign. long prev = modDiff(seed, 0); seed = prev; randArray[55] = prev; long next = 1L; for (int i = 21; i > 0; i = (i + 21) % 55) { randArray[i] = next; next = modDiff(prev, next); if ((seed & 1L) > 0) seed = 0x40000000L + (seed >> 1); else seed >>= 1; // Cyclic shift right 1. next = modDiff(next, seed); prev = randArray[i]; } for (int i = 0; i < 5; i++) flipCycle(); } /** * Computes 55 more pseudo-random numbers. * @return the next pseudo-random number. */ private long flipCycle() { for (int i = 1; i < 56; i++) { randArray[i] = modDiff(randArray[i], randArray[(i + 30) % 55 + 1]); } index = 54; return randArray[55]; } /** * Computes difference modulo 2^31. * @param x * @param y * @return */ private long modDiff(long x, long y) { return (x - y) & 0x7fffffffL; } }[…] 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 […]