Wednesday, February 18, 2015

elementary cellular automata

This code takes the rule number for an elementary cellular automaton as input, and then runs the CA from a random seed, rendering the result. The seed comprises 600 random bits taken from rule 30. Its length, when accounting for scale, is greater than the window width; this helps keep pathological edge effects outside the visible frame. The screenshot at left shows an execution of the famous rule 110.
import Control.Monad (when)
import Graphics.UI.SDL as SDL
import System.Environment (getArgs)
(xres, yres, unit) = (800, 600, 2)
main = withInit [InitVideo] $ do
[n] <- getArgs
w <- setVideoMode xres yres 32 [NoFrame]
enableEvent SDLMouseMotion False
setCaption "Elementary CA" "Elementary CA"
pause w $ zip (concat $ cells (read n) seed) xys
pause w cs = do
delay 128
e <- pollEvent
case e of
KeyUp (Keysym SDLK_ESCAPE _ _) -> return ()
KeyUp (Keysym SDLK_SPACE _ _) -> run w cs
_ -> pause w cs
run w (b:bs) = do
drawCell w b
e <- pollEvent
case e of
KeyUp (Keysym SDLK_ESCAPE _ _) -> return ()
KeyUp (Keysym SDLK_SPACE _ _) -> pause w $ b:bs
_ -> run w bs
drawCell w (b, (x,y)) = do
f =<< createRGBSurface [SWSurface] unit unit 32 0 0 0 0
when (x `mod` 73 == 0) $ SDL.flip w
where
rect = Just $ Rect x y unit unit
rgb = if b == 1 then 0 else 0xFFFFFF
f c = do fillRect c Nothing $ Pixel rgb
blitSurface c Nothing w rect
xys = concat . zipWith zip xs $ map repeat [0, unit..]
where
xs = offset $ iterate f [-299, -298.. 300]
offset = map . map $ (+ (xres `div` 2)) . (* unit)
f (n:ns) = enumFromTo (n-1) $ last ns + 1
-----------------------------------------------------
seed = take 600 . map mid . tail $ cells 30 [1]
where
mid xs = xs !! (length xs `div` 2)
cells n = iterate next
where
next = map f . chunk . pad
f = maybe 0 id . (`lookup` rule n)
pad s = [0,0] ++ s ++ [0,0]
rule = zip ns . bin
where
ns = map (drop 5 . bin) [7, 6..0]
bin n = map (fromEnum . (`elem` bits n)) pows
where
pows = reverse $ map (2^) [0..7]
bits 0 = []
bits n = let x = f n in x : bits (n-x)
where
f n = last $ takeWhile (<= n) $ map (2^) [0..]
chunk (a:b:c:ns) = [a,b,c] : chunk (b:c:ns)
chunk _ = []

Monday, February 16, 2015

wolfram's random generator

Despite its fundamental simplicity, the rule 30 elementary CA is conjectured to generate a purely random sequence along its center column. It reputably excels at many statistical tests of randomness, and Mathematica includes it as a choice of RNG.

For my own conviction, it was enough to find that the expression sum (take 80 $ rands 5000) `div` 80 produces 2525.

Since each random bit requires computing a full cell generation, this algorithm quickly slows down. I assume more pragmatic implementations solve this by perhaps re-seeding the automaton after some number of iterations.

rands n = map ((`mod` n) . dec) $ f cells
where
cells = map mid . tail $ iterate next [1]
mid xs = xs !! (length xs `div` 2)
f xs = take 32 xs : f (drop 32 xs)
next = map f . chunk . pad
where
f = maybe 0 id . (`lookup` rule 30)
pad s = [0,0] ++ s ++ [0,0]
rule = zip ns . bin
where
ns = map (drop 5 . bin) [7, 6..0]
dec = sum . zipWith (*) ns
where
ns = map (2^) [31, 30.. 0]
bin n = map (fromEnum . (`elem` bits n)) pows
where
pows = reverse $ map (2^) [0..7]
bits 0 = []
bits n = let x = f n in x : bits (n-x)
where
f n = last $ takeWhile (<= n) $ map (2^) [0..]
chunk (a:b:c:ns) = [a,b,c] : chunk (b:c:ns)
chunk _ = []