George Marsaglia’s Random Number Generators

October 5, 2010

Marsaglia used 32-bit unsigned integers, so all the arithmetic is modulo 232, which C does by simply overflowing the register but we must do with an explicit call to modulo; note though that arithmetic on the c variable is modulo 28. Bit operations are expressed more compactly in C than Scheme, and also seem to work a “bit” faster, if you’ll pardon my pun. The let makes the state variables of the random number generators private, rather than global as in the C program:

(define mwc #f)
(define shr3 #f)
(define cong #f)
(define fib #f)
(define kiss #f)
(define lfib4 #f)
(define swb #f)
(define uni #f)
(define vni #f)
(define settable #f)

(let ((z 362436069) (w 521288629) (jsr 123456789)
      (jcong 380116160) (a 224466889) (b 7584631)
      (t (make-vector 256 0)) (x 0) (y 0) (c 0))

  (define (mod8 n) (modulo n 256))
  (define (mod32 n) (modulo n 4294967296))
  (define (ref i) (vector-ref t (mod8 i)))

  (set! mwc (lambda ()
    (set! z (mod32 (+ (* 36969 (logand z 65535)) (ash z -16))))
    (set! w (mod32 (+ (* 18000 (logand w 65535)) (ash w -16))))
    (mod32 (+ (ash z 16) w))))

  (set! shr3 (lambda ()
    (set! jsr (mod32 (logxor jsr (ash jsr 17))))
    (set! jsr (mod32 (logxor jsr (ash jsr -13))))
    (set! jsr (mod32 (logxor jsr (ash jsr 5)))) jsr))

  (set! cong (lambda ()
    (set! jcong (mod32 (+ (* 69069 jcong) 1234567))) jcong))

  (set! fib (lambda ()
    (set! b (mod32 (+ a b))) (set! a (mod32 (- b a))) a))

  (set! kiss (lambda ()
    (mod32 (+ (logxor (mwc) (cong)) (shr3)))))

  (set! lfib4 (lambda ()
    (set! c (mod8 (+ c 1)))
    (vector-set! t c (mod32 (+ (ref c) (ref (+ c 58))
      (ref (+ c 119)) (ref (+ c 178))))) (ref c)))

  (set! swb (lambda ()
    (set! c (mod8 (+ c 1)))
    (let ((bro (if (< x y) 1 0)))
      (set! x (mod32 (ref (+ c 34))))
      (set! y (mod32 (+ (ref (+ c 19)) bro)))
      (vector-set! t c (mod32 (- x y)))
      (vector-ref t c))))

  (set! uni (lambda ()
    (* (kiss) 2.328306e-10)))

  (set! vni (lambda ()
    (* (- (kiss) 2147483648) 4.6566133-10)))

  (set! settable (lambda (i1 i2 i3 i4 i5 i6)
    (set! z i1) (set! w i2) (set! jsr i3) (set! jcong i4)
    (set! a i5) (set! b i6) (set! x 0) (set! y 0) (set! c 0)
    (do ((i 0 (+ i 1))) ((= i 256))
      (vector-set! t i (kiss))))))

The random number generators can be a lot faster if your Scheme system provides 32-bit unsigned integers and the corresponding bit operations on them, either natively or by some kind of foreign function interface. If your Scheme system doesn’t provide them, you will need logand, logxor and ash from the Standard Prelude; you will also need ipow, which is called by ash. And you will need the assert macro from the Standard Prelude for our version of the test program:

(define (test-rng)
  (let ((k 0))
    (settable 12345 65435 34221 12345 9983651 95746118)
    (do ((i 0 (+ i 1))) ((= i 1e6) (assert k 1064612766)) (set! k (lfib4)))
    (do ((i 0 (+ i 1))) ((= i 1e6) (assert k 627749721)) (set! k (swb)))
    (do ((i 0 (+ i 1))) ((= i 1e6) (assert k 1372460312)) (set! k (kiss)))
    (do ((i 0 (+ i 1))) ((= i 1e6) (assert k 1529210297)) (set! k (cong)))
    (do ((i 0 (+ i 1))) ((= i 1e6) (assert k 2642725982)) (set! k (shr3)))
    (do ((i 0 (+ i 1))) ((= i 1e6) (assert k 904977562)) (set! k (mwc)))
    (do ((i 0 (+ i 1))) ((= i 1e6) (assert k 3519793928)) (set! k (fib)))))

You can run the program at http://programmingpraxis.codepad.org/sf8Z4pJP.

About these ads

Pages: 1 2

3 Responses to “George Marsaglia’s Random Number Generators”

  1. [...] have examined several different random number generators in past exercises, including Donald Knuth’s lagged-fibonacci generator that is used in the Standard Prelude. We [...]

  2. Note that the SHR3 generator is flawed: it doesn’t actually have a period of 2^32-1 but several different cycles, some with quite short periods. Probably a typo: should instead be:

    #define SHR3  (jsr^=(jsr<>17), jsr^=(jsr<<5))

    That is, the 13 and 17 have been swapped. This generator has the period of 2^32-1 as intended.

    Beware also that for MWC generator you have to avoid certain "bad" seeds.

    I've implemented some of Marsaglia's generators in C and Python. Also with proper checking of valid seeds. See:

    https://bitbucket.org/cmcqueen1975/simplerandom/wiki/Home

  3. Sorry the previous comment had its #define line messed up by blog formatting, and I don’t know how to fix it.

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

%d bloggers like this: