Rule 30 RNG

April 29, 2011

We begin with some global variables. Size is the length of a row; research shows that an odd length is slightly better than an even length, and we want at least 200 cells in the row. Result is the index of the central cell that will be the source of the output bits. State is a vector that holds the current row of the generator:

(define size (* 7 43))
(define result (quotient size 2))
(define state (make-vector size 0))

We begin with the function that steps from one state to the next, returning the result bit and updating the state as a side effect. We implement Rule 30 by adding the bit-values of the three neighboring cells as a binary number, then looking up the new state. The two if statements implement the wrap-around, and the rule30 vector holds the bits of a binary 30 in reverse order:

(define rule30 #(0 1 1 1 1 0 0 0))

(define (step)
  (let ((next (make-vector size 0)))
    (do ((i 0 (+ i 1))) ((= i size))
      (vector-set! next i
        (vector-ref rule30
          (+ (* 4 (if (zero? i)
                      (vector-ref state (- size 1))
                      (vector-ref state (- i 1))))
             (* 2 (vector-ref state i))
             (if (= i (- size 1))
                 (vector-ref state 0)
                 (vector-ref state (+ i 1)))))))
    (set! state next)
    (vector-ref state result)))

If you are clever, you can operate on bits in a block as long as the longest integer type available, which is much faster than operating on each bit separately. But that requires special handling at the ends of the blocks, and isn’t really feasible in Scheme, even with the bitwise operators of R6RS Scheme, because there are no native bit-vectors. Thus we use the bit-at-a-time method.

Step makes it is easy to get successive states, accumulating them until the desired number is reached:

(define (get n)
  (let loop ((n n) (xs '()))
    (if (zero? n)
        (undigits (reverse xs) 2)
        (loop (- n 1) (cons (step) xs)))))

There is no pre-defined way to initialize the state, so we invent our own. We’ll take a pass phrase, like Programming Praxis, expand it as a base-256 number modulo the block size, and copy it into each initial block of the state. Init-state takes one or more numbers, each the bits for a block, and cycles if there aren’t enough blocks given. Note that we run the generator for a little while (state size times the number of blocks) to get it “warmed up” for use. The block size must be a divisor of the state size:

(define block 43) ; must be a divisor of size

(define (init-state . xs)
  (set! state (make-vector size 0))
  (do ((i (- block 1) (+ i block))
       (xs (cycle xs) (cdr xs)))
      ((< size i))
    (do ((j i (- j 1))
         (ds (reverse (digits (car xs) 2)) (cdr ds)))
        ((null? ds))
      (vector-set! state j (car ds))))
  (do ((i 0 (+ i 1))) ((< (* size size (/ block)) i))
    (step)))

(define (pass n str)
  (let loop ((ss (string->list str)) (sum 0))
    (if (null? ss) sum
      (loop (cdr ss)
            (modulo (+ sum (* sum 256)
                       (char->integer (car ss)))
                    (expt 2 n))))))

(init-state (pass block "Programming Praxis"))

In practice, you would probably want the block size to be either 32 or 64, depending on the native word size for your machine; 256 makes a good state size. On the other hand, there is some evidence that odd state sizes make better random number generators. We choose a block size of 43 and state size of 7 × 43 = 301. For testing, you can call (init-state 0) to clear the state, then (vector-set! state result 1) turns on a single bit in the center of the state, as Wolfram did in his book; then the next 32 bits are 1, 0, 1, 1, 1, 0, 0, 1, 1, 0, 0, 0, 1, 0, 1, 1, 0, 0, 1, 0, 0, 1, 1, 1, 0, 1, 0, 1, 1, 1, 0, and 0 (Sloane’s A051023), which you can see by saying (get 32), which will return 3112904540.

We used cycle, digits and undigits from the Standard Prelude. You can run the program at http://programmingpraxis.codepad.org/zoGntTfu.

About these ads

Pages: 1 2

9 Responses to “Rule 30 RNG”

  1. [...] today’s Programming Praxis exercise, our goal is to write a random number generator based on the Rule 30 [...]

  2. My Haskell solution (see http://bonsaicode.wordpress.com/2011/04/29/programming-praxis-rule-30-rng/ for a version with comments):

    import Control.Monad.State
    import Data.Bits
    import Data.Word
    
    size :: Int
    size = 7 * 43
    
    fromBits :: (Bits a, Integral a) => [Bool] -> a
    fromBits = foldl (\a x -> shift a 1 .|. fromIntegral (fromEnum x)) 0
    
    step :: Monad m => Int -> StateT Integer m Bool
    step r = modify (\s -> let b i = testBit s (mod i size) in
        fromBits [ testBit r $ fromBits [b (n+1), b n, b (n - 1)]
                 | n <- [size - 1, size - 2..0]]) >> gets (`testBit` div size 2)
    
    randomBits :: Monad m => Int -> Int -> StateT Integer m Word
    randomBits n r = return . fromBits =<< replicateM n (step r)
    
    main :: IO ()
    main = evalStateT (do r <- randomBits 32 30
                          liftIO $ print (r == 3112904540)
                          liftIO . print =<< replicateM 100 (randomBits 8 30)
                      ) $ bit (div size 2)
    
  3. arturasl said

    Don’t really know if I understand algorithm, but here is my attempt in python:

    #!/usr/bin/python3
    
    from math import log, floor
    from functools import reduce
    
    lst_lookup = [1, 1, 1, 0, 0, 0, 0, 1]
    
    (state, iterations) = 123456789, 5
    size = floor(log(state, 2)) + 1
    main_index = size // 2
    
    def gen_states2(n):
        cur_state = state
        for i in range(0, n):
            print([int((cur_state & (1 << i)) != 0) for i in range(size - 1, -1, -1)])
            yield cur_state
    
            cur_state = (cur_state << 1) | ((cur_state & 1) << (size + 1)) | (cur_state & (1 << (size - 1))) >> (size - 1)
            cur_state = reduce(lambda x, y: x + (lst_lookup[(cur_state >> y) & 7] << y), range(0, size))
    
    print(reduce(lambda x, y: x * 2 + y, [((i >> main_index) & 1) for i in gen_states2(iterations)], 0))
    

    Not the most efficient, but I wanted to look at python’s bitwise capabilities and to play with all that funky stuff Graham uses all the time :D
    Basically I use simple number as a state (each binary digit represents cell) and then use gen_state2 to generate all the preceding states.

    [1, 1, 1, 0, 1, 0, 1, 1, 0, 1, 1, 1, 1, 0, 0, 1, 1, 0, 1, 0, 0, 0, 1, 0, 1, 0, 1]
    [1, 1, 0, 0, 1, 0, 0, 0, 0, 0, 1, 1, 0, 0, 1, 0, 0, 0, 1, 0, 1, 1, 1, 0, 1, 0, 0]
    [0, 0, 0, 1, 1, 0, 1, 1, 1, 1, 0, 0, 0, 1, 1, 0, 1, 1, 1, 0, 0, 1, 0, 0, 1, 0, 0]
    [1, 1, 1, 0, 0, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 0, 0, 1, 0, 0, 1, 1, 0, 1, 1, 0, 0]
    [0, 1, 0, 0, 1, 1, 1, 0, 0, 0, 1, 0, 0, 0, 1, 1, 1, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0]
    4
    
  4. Mike said

    Implemented a Rule30 class. An instance of the class can be called like a function to return the next random number. The class also supports the iteration protocol.

    The basic strategy is to store the state of the rng in an integer. left and right are the state rotated 1 bit to the right and left, respectively, so the three bits that make up the next state “line up”. Then a boolean expression involving left, state and right is used to calculate the next state.

    The LSBs of the state are accumulated to determine the random number

    from random import randint
    
    class Rule30:
        def __init__(self, seed=None, nbits=32):
            self.mask  = 2**(2*nbits) - 1
            self.nbits = nbits
            self.shift = 2*nbits - 1
    
            self.state  = seed or randint(1,self.mask)
    
        def __call__(self):
            val = 0
    
            for n in range(self.nbits):
                left  = (self.state >> 1 | self.state << self.shift) & self.mask
                right = (self.state << 1 | self.state >> self.shift) & self.mask
    
                self.state = ~left&self.state | ~self.state&(left^right)
    
                val = val<<1 | self.state&1
    
            return val
    
        def __iter__(self):
            return self
    
        next = __call__
    
    # examples
    from itertools import islice
    
    # using an instance as an iterable/generator
    for n in islice(Rule30(),10):
        print n
    
    # calling an instance
    rng = Rule30(3495872309)
    rng()
    
  5. David said

    Forth version. We take a 30 bit number, prep it by shifting once left and wrapping the bits. Each iteration of the loop, use table to look up the value of current bit based on three bits in original word.

    \ 30 bit random
    
    create state-trans
        1 C,  1 C,  1 C,  0 C,  0 C,  0 C,  0 C,  1 C,
    
    : prep  ( n -- n )  \ takes 30 bit number, wraps around low bit & high bit
        $3FFFFFFF and 1 lshift
        dup $00000002 and IF $80000000 or THEN
        dup $40000000 and IF $00000001 or THEN ;
    
    : rng30 ( n -- n )
        0 locals| result |
        prep 7
        30 0 DO
            2dup and i rshift state-trans + c@
            i lshift result or  to result
            2*
         LOOP
         2drop result ;    
    
    create seed  counter $3FFFFFFF and ,
    : random-30
        seed @  dup rng30  seed ! ;
    
    : .30bin 
        base @ >r  2 base !
        0 <# 30 0 do # loop #> type space
        r> base ! ;
    

    Session:

    include rng30  ok
    random-30 .30bin 000000111001010111100101100100  ok
    random-30 .30bin 111111010011010011001100001101  ok
    random-30 .30bin 111110010100010100010001110000  ok
    random-30 .30bin 011100110101110101110110100111  ok
    random-30 .30bin 001001000100100100100000101010  ok
    random-30 .30bin 111011011101101101101111101010  ok
    random-30 .30bin 010000001000000000000111001010  ok
    random-30 .30bin 110111111011111111111010011010  ok
    random-30 .30bin 000011110001111111110010100010  ok
    random-30 .30bin 111101100110111111100110101110  ok
    
  6. Graham said

    @arturasl, you flatter me. As for me, here’s my Python solution.
    Perhaps not as elegant as others, but I wanted to play with Python’s bitarrays for speed.

  7. Graham said

    @arturasl: I just try to make my Python solutions get anywhere near the great functional programming
    we see in the “official” solution and in the Haskell ones @Remco posts :-)

  8. JP said

    That was pretty neat. I’m not sure how valuable it is as a RNG, as there are easier to code RNGs that seem more random (heck, this post from here at Programming Praxis even) but as a thought experiment it’s worthwhile. And it gave me a chance to build on the cellular automaton code I’d previously built in JavaScript and transfer that to Racket.

    In any case, here’s my Racket version: Rule 30 RNG

  9. [...] number generator in specific: a Rule 30¬†psuedo-random number generator¬†(PRNG). (Here’s the motivating post from Programming [...]

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: