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-24an-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.

About these ads

Pages: 1 2

4 Responses to “GB_FLIP”

  1. Graham said

    I’ve been trying to come up with a Python version, but the global variables are giving me the run around.

  2. Mike said

    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!"
    
    
  3. [...] 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. [...]

  4. Ying Yin said

    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;
    	}
    }
    

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

%d bloggers like this: