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.

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

%d bloggers like this: