connett circles

Like Barry Martin's 'Hopalong' fractal, this dynamical system from John Connett was first published in Scientific American in 1986. This demo is interactive: successively clicking two points specifies a rectangle to zoom into. Doing so, you'll see that the system isn't actually a fractal. Instead of self-similarity, deep zooms reveal peacock-like images.

	{-# LANGUAGE NoMonomorphismRestriction #-}
	import Data.Complex (magnitude, Complex((:+)))
	import Graphics.UI.SDL as SDL
	(xres, yres, cast) = (768, 768, fromIntegral)
	main = withInit [InitVideo] $ do
	w <- setVideoMode xres yres 32 [NoFrame]
	enableEvent SDLMouseMotion False
	setCaption "Connett Circles" "Connett Circles"
	zoom w (-2000, -2000) (2000, 2000)
	zoom w (x,y) (s,t) = do
	let r = area x s y t
	render w r
	run w r (0,0)
	run w (xs,ys) p@(x,y) = do
	delay 32
	e <- pollEvent
	case e of
	KeyUp (Keysym SDLK_ESCAPE _ _) -> return ()
	MouseButtonUp i j _ -> click $ f i j
	_ -> run w (xs,ys) p
	where
	f i j = (xs !! cast i, ys !! cast j)
	click q
	| p == (0,0) = run w (xs,ys) q
	| otherwise = zoom w p q
	render w r = draw w . concat . f $ set r
	where
	f = zipWith zip $ map (zip [0..] . repeat) [0..]
	draw w r = do
	s <- createRGBSurface [] 1 1 32 0 0 0 0
	mapM_ (draw s) r
	SDL.flip w
	where
	rect (x,y) = Just $ Rect x y 1 1
	draw s (p,n) = do fillRect s Nothing $ Pixel $ rgb n
	blitSurface s Nothing w $ rect p
	------------------------------------------------------
	rgb n
	| m `mod` 7 == 0 = 0xFFFFFF
	| m `mod` 6 == 0 = 0x00FFFF
	| m `mod` 5 == 0 = 0xFFFF00
	| m `mod` 4 == 0 = 0xFF0000
	| m `mod` 3 == 0 = 0xFF00
	| m `mod` 2 == 0 = 0xFF
	| otherwise = 0
	where
	f x y = round $ magnitude $ (cast x^2) :+ (cast y^2)
	m = abs n
	set (xs, ys) = [[f x y | x <- xs] | y <- ys]
	where
	f x y = round $ magnitude $ (x^2) :+ (y^2)
	area x x' y y' = ([t, t+a.. s], [u, u-b.. v])
	where
	(a,b) = ((s - t) / dx, (u - v) / dy)
	(s,t) = (max x x', min x x')
	(u,v) = (max y y', min y y')
	(dx, dy) = (cast xres, cast yres)

view raw connett_circles.hs hosted with ❤ by GitHub

martin attractor

This pattern generator, discovered by Barry Martin, was nicknamed 'Hopalong' when Scientific American introduced it in their September '86 issue.

Clicking the window adjusts the viewport position; there is also an alternate version with color and animation.

Also, for a certain Rubyist friend, I wrote another lazily-evaluated, colored and animated implementation in Ruby.

	import Control.Arrow ((***))
	import Graphics.UI.SDL as SDL
	(xres, yres) = (850, 850)
	(a,b,c) = (5, 1, 20)
	main = withInit [InitVideo] $ do
	w <- setVideoMode xres yres 32 [NoFrame]
	enableEvent SDLMouseMotion False
	setCaption "Hopalong Pattern" "Barry Martin"
	render w p ps
	loop w p ps
	where
	p = (xres `div` 2, yres `div` 2)
	ps = take 7000000 points
	loop w p ps = do
	delay 32
	e <- pollEvent
	case e of
	KeyUp (Keysym SDLK_ESCAPE _ _) -> return ()
	MouseButtonUp x y _ -> click $ f x y
	_ -> loop w p ps
	where
	click p = render w p ps >> loop w p ps
	f x y = (fromIntegral x, fromIntegral y)
	render w (x,y) ps = do
	s <- createRGBSurface [SWSurface] 1 1 32 0 0 0 0
	fillRect w (Just $ Rect 0 0 xres yres) (Pixel 0)
	mapM_ (draw s 0xFFFFFF . rect . center) ps
	SDL.flip w
	where
	center = ((xres - x) +) *** (+ (yres - y))
	rect (x,y) = Just $ Rect x y 1 1
	draw s c p = do fillRect s Nothing $ Pixel c
	blitSurface s Nothing w p
	points = map (round *** round) $ iterate f (0,0)
	where
	f (x,y) = (g x y, a - x)
	g x y = y - (signum x * (sqrt . abs $ b*x - c))

view raw barry-martin-hopalong.hs hosted with ❤ by GitHub

lyapunov fractals

It took me some experimentation to figure out how to color this derivation of the logistic map; I'm still not quite sure how the hues should scale as you zoom. But the bi-tonal method shown below works well enough to produce the image at left - click it for more detail.

	{-# LANGUAGE NoMonomorphismRestriction #-}
	import Control.Monad (liftM2)
	import Data.Bits (shiftL)
	import Graphics.UI.SDL as SDL
	(xres, yres, cast) = (400, 400, fromIntegral)
	main = withInit [InitVideo] $ do
	w <- setVideoMode xres yres 32 [NoFrame]
	enableEvent SDLMouseMotion False
	setCaption "Lyapunov" "Lyapunov"
	render w 1024 (3.2, 3.2) (3.8, 3.8)
	render w i (x,y) (s,t) = do
	draw w (zip ps $ area x s y t) i
	run w
	where
	ps = liftM2 (,) [0.. xres] [0.. yres]
	draw w cs i = do
	s <- createRGBSurface [] 1 1 32 0 0 0 0
	mapM_ (draw s) cs
	SDL.flip w
	where
	draw s (p,q) = do
	let rect (x,y) = Just $ Rect x y 1 1
	fillRect s Nothing $ rgb $ lyapunov i q
	blitSurface s Nothing w $ rect p
	run w = do
	delay 32
	e <- pollEvent
	case e of
	KeyUp (Keysym SDLK_ESCAPE _ _) -> return ()
	_ -> run w
	-----------------------------------------------------
	rgb n | n < 0 = Pixel $ f n
	| True = Pixel $ f n `shiftL` 16
	where
	f x = round (min 255 $ abs n * 255)
	lyapunov i p = (1 / cast i) * sum rs
	where
	rs = filter (not . isInfinite) $ f xs
	f = map (log . abs) . zipWith ($) (str p)
	xs = map ((1 -) . (* 2)) $ verhulst i p
	verhulst i p = take i $ go 0.5 $ str p
	where
	go x (f:fs) = x : go (f (x * (1-x))) fs
	str (x,y) = map (*) $ concat $ repeat [x,x,y,x,y]
	area x x' y y' = liftM2 (,) [s, s+a.. t] [u, u+b.. v]
	where
	(a,b) = ((t - s) / m, (v - u) / n)
	(s,t) = (min x x', max x x')
	(u,v) = (min y y', max y y')
	(m,n) = (cast xres, cast yres)

view raw lyapunov-fractals.hs hosted with ❤ by GitHub

the mandelbrot set

What programmer hasn't at some point written an implementation of Benoit Mandelbrot's great discovery, the most famous fractal in the world? Here's my own minimal version, with the simplest possible coloring scheme. To interact with it, just click any two points: the window will zoom in on the rectangle they define.

	import Data.Complex (magnitude, Complex((:+)))
	import Graphics.UI.SDL as SDL
	(xres, yres) = (900, 750)
	main = withInit [InitVideo] $ do
	w <- setVideoMode xres yres 32 [NoFrame]
	enableEvent SDLMouseMotion False
	setCaption "Mandelbrot" "Mandelbrot"
	zoom w 1 (-2, 1) (0.5, -1)
	zoom w i (x,y) (s,t) = do
	let r = area x s y t
	render w r i
	run w (i + 1) r (0,0)
	render w r i = draw w . concat . f $ set i r
	where
	f = zipWith zip $ map (zip [0..] . repeat) [0..]
	draw w cs = do
	s <- createRGBSurface [] 1 1 32 0 0 0 0
	mapM_ (draw s) cs
	SDL.flip w
	where
	rgb n = div (floor $ min n 2 * 512) 2
	rect (x,y) = Just $ Rect x y 1 1
	draw s (p,n) = do fillRect s Nothing $ Pixel $ rgb n
	blitSurface s Nothing w $ rect p
	run w i (xs,ys) p@(x,y) = do
	e <- pollEvent
	case e of
	KeyUp (Keysym SDLK_ESCAPE _ _) -> return ()
	MouseButtonUp i j _ -> click $ f i j
	_ -> run w i (xs,ys) p
	where
	f i j = (xs !! fromIntegral i, ys !! fromIntegral j)
	click q
	| p == (0,0) = run w i (xs,ys) q
	| otherwise = zoom w i p q
	------------------------------------------------------
	area x x' y y' = ([t, t+a.. s], [u, u-b.. v])
	where
	(a,b) = ((s - t) / dx, (u - v) / dy)
	(s,t) = (max x x', min x x')
	(u,v) = (max y y', min y y')
	(dx, dy) = (fromIntegral xres, fromIntegral yres)
	set i (xs, ys) = [[f x y | x <- xs] | y <- ys]
	where
	f x = magnitude . g . (x :+)
	g a = iterate ((+ a) . (^2)) 0 !! (i * 20)

view raw mandelbrot_set_viewer.hs hosted with ❤ by GitHub

a lava lamp

This 'lava lamp' is actually a simple cyclic cellular automaton; the CA rule is courtesy of Jason Rampe. In keeping with the spirit of Haskell, I chose to implement it with a hashmap of points rather than an array. Needless to say, that's not practical; but this is just a demonstration.

	import Control.Arrow ((***))
	import Control.Monad (join, liftM2)
	import Data.List (unfoldr)
	import Data.Maybe (fromJust)
	import qualified Data.HashMap.Strict as M
	import Graphics.UI.SDL as SDL
	import System.Random.Mersenne.Pure64 (newPureMT, randomInt)
	(xres, yres, unit) = (160, 90, 10)
	main = withInit [InitVideo] $ do
	w <- setVideoMode (f xres) (f yres) 32 [NoFrame]
	enableEvent SDLMouseMotion False
	setCaption "Cyclic CA" "Cyclic CA"
	pause w =<< randPattern
	where
	f = (* unit)
	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 cs = do
	drawCells w cs
	SDL.flip w
	e <- pollEvent
	case e of
	KeyUp (Keysym SDLK_ESCAPE _ _) -> return ()
	KeyUp (Keysym SDLK_SPACE _ _) -> pause w cs
	_ -> run w cs'
	where
	cs' = M.mapWithKey (next cs) cs
	drawCells w cs = do
	s <- createRGBSurface [] unit unit 32 0 0 0 0
	mapM_ (draw s) $ M.toList cs
	where
	rgb n = [0xFF0000, 0xAA0000, 0x770000] !! n
	rect (x,y) = Just $ Rect x y unit unit
	scale = rect . join (***) (* unit)
	draw s (p,n) = do fillRect s Nothing $ Pixel $ rgb n
	blitSurface s Nothing w $ scale p
	randPattern = fmap (M.fromList . zip ps) ns
	where
	g = map (`mod` 3) . unfoldr (Just . randomInt)
	ns = (return . g) =<< newPureMT
	ps = liftM2 (,) [0.. xres - 1] [0.. yres - 1]
	next cs p n = if m >= 10 then n' else n
	where
	m = length $ filter (== n') area
	area = map (fromJust . (`M.lookup` cs)) $ moore p
	n' = (n + 1) `mod` 3
	moore (x,y) = filter p ps
	where
	ps = liftM2 (,) [x-2.. x+2] [y-2.. y+2]
	p (x,y) = x >= 0 && x < xres && y >= 0 && y < yres

view raw lava-lamp-automaton.hs hosted with ❤ by GitHub

random undirected graphs

This little program generates random, undirected graphs without loops or isolated vertices (though the graph is not necessarily connected).

I have found it useful to generate random inputs for testing graph algorithms.

	import Data.List (nub, unfoldr)
	import Control.Arrow ((***))
	import Control.Monad (join, replicateM)
	import Graphics.UI.SDL as SDL
	import Graphics.UI.SDL.Primitives (line)
	import System.Random.Shuffle (shuffle')
	import System.Random.Mersenne.Pure64 (newPureMT, randomInt)
	[xres, yres, unit] = map (1 *) [512, 512, 5]
	main = withInit [InitVideo] $ do
	window <- setVideoMode xres yres 32 [NoFrame]
	enableEvent SDLMouseMotion False
	setCaption "Random Graphs" "Random Graphs"
	render window =<< genGraph
	run
	run = do
	delay 128
	e <- pollEvent
	case e of
	KeyUp (Keysym SDLK_ESCAPE _ _) -> return ()
	_ -> run
	render w ls = do
	s <- createRGBSurface [SWSurface] unit unit 32 0 0 0 0
	fillRect s Nothing $ Pixel 0xFFFFFF
	fillRect w (Just $ Rect 0 0 xres yres) (Pixel 0)
	mapM_ (draw w s) ls
	SDL.flip w
	draw w s ((a,b), (x,y)) = do
	blitSurface s Nothing w $ rect (a,b)
	blitSurface s Nothing w $ rect (x,y)
	line w a' b' x' y' $ Pixel 0xFF0000FF
	where
	rect (x,y) = Just $ Rect x y unit unit
	[a', b', x', y'] = map ((+2) . fromIntegral) [a,b,x,y]
	genGraph = do
	gs <- replicateM 2 newPureMT
	ps <- return . nub . take 15 $ rands gs
	qs <- permute ps
	return $ filter p $ zip ps qs
	where
	permute xs = fmap (shuffle' xs $ length xs) newPureMT
	p ((a,b), (x,y)) = (a,b) /= (y,x) && (a,b) /= (x,y)
	rands [g, g'] = map (f . h) $ zip xs ys
	where
	xs = map ((8+) . (`mod` x)) $ unfoldr (Just . randomInt) g
	ys = map ((8+) . (`mod` y)) $ unfoldr (Just . randomInt) g'
	(x,y) = (xres - 16, yres - 16)
	(f,h) = (join (***) (* unit), join (***) (`div` unit))

view raw random-graphs.hs hosted with ❤ by GitHub

sqrt 2, visualized

This code, inspired by my earlier post, animates the digits of √2 in an unsuual way. Each digit advances the curve in the direction given by a numeric keypad, with 5 and 0 both considered (0,0). The thumbnail at left links to a 4000x4000 window on 9856041 digits of the curve, which originates at the image center.

The program makes for an interesting screensaver, but does it follow any pattern, or illustrate any special properties of √2 ? Though I'm no mathematician, from what I've read this seems unlikely. √2 is suspected (but not proved) to be a normal number. This would mean the digits of its expansion (respective to some number base) follow a uniform distribution; no digit would be more likely to appear than any other.

Nonetheless, I'd be curious to see the results of a really long run!

	import Control.Arrow ((***))
	import Control.Monad (join, liftM2, liftM3, when)
	import Data.Bits (shift)
	import Data.Maybe (fromJust)
	import Graphics.UI.SDL as SDL
	(xres, yres) = (1600, 900)
	main = withInit [InitVideo] $ do
	w <- setVideoMode xres yres 32 [NoFrame]
	s <- createRGBSurface [SWSurface] 1 1 32 0 0 0 0
	enableEvent SDLMouseMotion False
	setCaption "Sqrt 2 Curve" "Sqrt 2 Curve"
	run w s (xres `div` 2, yres `div` 2) dirs rgbs [1..]
	run w s p (f:fs) (c:cs) (n:ns) = do
	drawCell w s p c
	when (n `mod` 47 == 0) $ SDL.flip w
	e <- pollEvent
	case e of
	KeyUp (Keysym SDLK_ESCAPE _ _) -> save >> print n
	_ -> run w s (f p) fs cs ns
	where
	save = saveBMP w "out.bmp"
	drawCell w s (x,y) c = draw x y
	where
	rect x y = Just $ Rect x y 1 1
	draw x y = do fillRect s (rect 0 0) (Pixel c)
	blitSurface s (rect 0 0) w $ rect x y
	sqrt2 = let ns = f 2 0 in concatMap (show . (ns !!)) [0..]
	where
	f x r = d : f (100 * (x - (20 * r + d) * d)) (10 * r + d)
	where
	d = head (dropWhile p [0..]) - 1
	p n = (20 * r + n) * n < x
	dirs = map f . expand $ map numPad sqrt2
	where
	f (x,y) = (x+) *** (+y)
	rgbs = expand . cycle . map f $ liftM3 (,,) ns ns ns
	where
	ns = [71, 77.. 255]
	f (r,g,b) = shift r 16 + shift g 8 + b
	numPad = fromJust . (`lookup` pad)
	where
	pad = zip "0528639417" $ head ps : ps
	ps = (join $ liftM2 (,)) [0, 1, -1]
	expand (x:xs) = replicate 8 x ++ expand xs

view raw sqrt-two-visualized.hs hosted with ❤ by GitHub

triangle geometry

Here's a trivial post, mostly just to show off how nice mathematical code can look in Haskell. It illustrates the centroid, incircle, and circumcircle of a random triangle, with everything defined in terms of vertex triples. The viewport centers on the circumcenter, and the triangle's interior is filled using barycentric coordinates.

	import Control.Arrow ((***))
	import Control.Monad (ap, liftM2, void, when)
	import Data.List (unfoldr)
	import Graphics.UI.SDL as SDL
	import Graphics.UI.SDL.Primitives (circle, filledCircle, pixel)
	import System.Random.Mersenne.Pure64 (newPureMT, randomInt)
	(xres, yres) = (1600, 900)
	main = withInit [InitVideo] $ do
	w <- setVideoMode xres yres 32 [NoFrame]
	rng <- newPureMT
	enableEvent SDLMouseMotion False
	setCaption "Triangle" "Triangle"
	render w rng $ liftM2 (,) [0.. xres] [0.. yres]
	SDL.flip w
	pause w
	pause w = do
	delay 128
	e <- pollEvent
	case e of
	KeyUp (Keysym SDLK_ESCAPE _ _) -> return ()
	_ -> pause w
	render w rng ps = do
	let pxl (x,y) = void . pixel w x y $ Pixel 0x0000FFFF
	let (x,y) = centroid tri
	let (p,q) = circumcenter tri
	let (a,b,r) = incircle tri
	let rad = dist (p,q) $ head tri
	let f = round . (+ (fromIntegral xres / 2 - p))
	let g = round . (+ (fromIntegral yres / 2 - q))
	let h = (f . fromIntegral) *** (g . fromIntegral)
	mapM_ (ap (when . within tri) (pxl . h)) ps
	filledCircle w (f x) (g y) 2 $ Pixel 0xFFFF00FF
	filledCircle w (f a) (g b) 2 $ Pixel 0xFF00FF
	circle w (f p) (g q) rad $ Pixel 0xFF0000FF
	circle w (f a) (g b) (round r) $ Pixel 0xFF00FF
	where
	tri = head $ filter p $ triples rng
	p t = let n = area t in n > 15000 && n < 40000
	------------------------------------------------------------
	within [(ax,ay), (bx,by), (cx,cy)] (x', y') = f [a,b,g]
	where
	(x,y) = (fromIntegral x', fromIntegral y')
	a = ((by-cy) * (x-cx) + (cx-bx) * (y-cy)) / d
	b = ((cy-ay) * (x-cx) + (ax-cx) * (y-cy)) / d
	d = (by-cy) * (ax-cx) + (cx-bx) * (ay-cy)
	g = 1 - a - b;
	f = all (> 0)
	circumcenter [(ax,ay), (bx,by), (cx,cy)] = (x/d, y/d)
	where
	x = ((ax-cx) * (ax+cx) + (ay-cy) * (ay+cy)) / 2 * (by-cy)
	- ((bx-cx) * (bx+cx) + (by-cy) * (by+cy)) / 2 * (ay-cy)
	y = ((bx-cx) * (bx+cx) + (by-cy) * (by+cy)) / 2 * (ax-cx)
	- ((ax-cx) * (ax+cx) + (ay-cy) * (ay+cy)) / 2 * (bx-cx)
	d = (ax-cx) * (by-cy) - (bx-cx) * (ay-cy)
	incircle [(x0,y0), (x1,y1), (x2,y2)] = (x/d, y/d, r)
	where
	x = a*x0 + b*x1 + c*x2
	y = a*y0 + b*y1 + c*y2
	a = fromIntegral $ dist (x1,y1) (x2,y2)
	b = fromIntegral $ dist (x2,y2) (x0,y0)
	c = fromIntegral $ dist (x0,y0) (x1,y1)
	d = a + b + c
	r = let s = d/2 in sqrt ((s-a)*(s-b)*(s-c) / s)
	centroid [(ax,ay), (bx,by), (cx,cy)] = (x,y)
	where
	(x,y) = ((ax+bx+cx) / 3, (ay+by+cy) / 3)
	area [(ax,ay), (bx,by), (cx,cy)] = abs $ (a+b+c) / 2
	where
	a = ax * (by - cy)
	b = bx * (cy - ay)
	c = cx * (ay - by)
	dist (p,p') (q,q') = round . sqrt $ (p-q)^2 + (p'-q')^2
	triples rng = let f xs = take 3 xs : f (drop 3 xs) in f ps
	where
	ps = go $ map ((+ 10) . (`mod` 580)) rands
	rands = unfoldr (Just . randomInt) rng
	go (x:y:ns) = (fromIntegral x, fromIntegral y) : go ns

view raw triangle-geometry-demo.hs hosted with ❤ by GitHub

primitive totalistic automata

This code renders any of the 2187 possible 3-colored, 1-dimensional, totalistic cellular automata. I was charmed by these and many other beautiful demonstrations in Stephen Wolfram's notorious compendium, though I regret I can't say the same for its tendentious style.

The program input is an integer representing the intended CA rule in base 3.

	import Control.Monad (when)
	import Data.List (group, unfoldr)
	import Graphics.UI.SDL as SDL
	import System.Environment (getArgs)
	import System.Random.Mersenne.Pure64 (pureMT, randomInt)
	(xres, yres, unit) = (400, 400, 2)
	main = withInit [InitVideo] $ do
	[n] <- getArgs
	w <- setVideoMode xres yres 32 [NoFrame]
	enableEvent SDLMouseMotion False
	setCaption "Elementary CA" "Elementary CA"
	run w $ zip (concat $ cells (read n) seed) xys
	run w (b:bs) = do
	drawCell w b
	e <- pollEvent
	case e of
	KeyUp (Keysym SDLK_ESCAPE _ _) -> return ()
	_ -> run w bs
	drawCell w (b, (x,y)) = do
	f =<< createRGBSurface [SWSurface] unit unit 32 0 0 0 0
	when (x `mod` 157 == 0) $ SDL.flip w
	where
	rect = Just $ Rect x y unit unit
	rgb = maybe 0 id $ lookup b rgbs
	rgbs = [(0, 0xffffff), (1, 0x225577), (2,0)]
	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 (`mod` 3) ns
	where
	ns = unfoldr (Just . randomInt) $ pureMT 71
	cells n = iterate next
	where
	next = map (f . sum) . chunk . pad
	f = maybe 0 id . (`lookup` rule)
	rule = zip [6,5..0] $ tern n
	pad s = [0,0] ++ s ++ [0,0]
	chunk [_,_] = []
	chunk ns = take 3 ns : chunk (tail ns)
	tern n = map (length . f) pows
	where
	pows = map ((:[]) . (3^)) [6,5..0]
	f m | m `elem` group (bits n) = m
	f m | (m ++ m) `elem` group (bits n) = m ++ m
	| otherwise = []
	bits 0 = []
	bits n = let x = f n in x : bits (n-x)
	where
	f n = last $ takeWhile (<= n) $ map (3^) [0..]

view raw primitive-totalistic-automata.hs hosted with ❤ by GitHub

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 _ = []

view raw elementary-cellular-automata.hs hosted with ❤ by GitHub

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 _ = []

view raw rule-30-random-gen.hs hosted with ❤ by GitHub