GB_FLIP
May 25, 2010
The Stanford GraphBase of Donald Knuth provides a portable, high-quality random number generator called GB_FLIP. Based on the lagged fibonacci generator an = (an-24 − an-55) mod m, GB_FLIP provides values on the range zero inclusive to 231 exclusive, has a period of 285 − 230, and the low-order bits of the generated numbers are just as random as the high-order bits. You can read Knuth’s original version of the program at http://tex.loria.fr/sgb/gb_flip.pdf.
Your task is to write your own version of GB_FLIP. 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.
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 […]