## George Marsaglia’s Random Number Generators

### October 5, 2010

Marsaglia used 32-bit unsigned integers, so all the arithmetic is modulo 2^{32}, 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 2^{8}. 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.

[…] 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 […]

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

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

[…] 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 looked). […]