Wednesday, December 23, 2015

connett circles

Like Barry Martin's 'Hopalong' fractal, this dynamical system from John Connett was first published in Scientific American in 1986. However, it's not actually a fractal: when zooming in (by clicking two points to specify a rectangle), viewers will see the pattern is not infinitely detailed. Instead, peacock-like images appear.
{-# 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)

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

Update: 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))

Tuesday, December 22, 2015

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 right.
{-# 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)

Tuesday, December 15, 2015

the mandelbrot set

Benoit's famous fractal. Click any two points to specify a zoom rectangle.
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)

Monday, November 2, 2015

a lava lamp

This 'lava lamp' is actually a simple cyclic cellular automaton. The rule is courtesy of Jason Rampe.
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

Wednesday, August 19, 2015

random undirected graphs

This little program generates random, undirected graphs without loops or isolated vertices - though the graph is not necessarily connected.
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))

Wednesday, July 8, 2015

sqrt 2, visualized

This program, inspired by my earlier post, illustrates the digits of sqrt 2. Each digit advances the curve eight pixels in the direction given by a numeric keypad, with 5 and 0 both considered (0,0). The thumbnail at right links to a 4000x4000 window on 9856041 digits of the curve, which originates at the image center.
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

Sunday, June 28, 2015

triangle geometry

Illustrates the centroid, incircle, and circumcircle of a random triangle; the viewport centers on the circumcenter. All functions are defined in terms of vertex triples. The triangle 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

Monday, March 30, 2015

primitive totalistic automata

This code renders any of the 2187 possible 3-colored, 1-dimensional totalistic CAs; its input is an integer representing the 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..]

Wednesday, February 18, 2015

elementary cellular automata

Takes the rule number for an elementary CA as input. The seed is provided by 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 right shows the infamous 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

The center column of the rule 30 elementary CA: a chaotic bit stream.

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'd imagine that robust implementations (Mathematica uses this PRNG, among others) re-seed the automaton after every so many bits.

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