128-bit AES electronic codebook

Rijndael is an algorithm that might actually be more succinct in C than Haskell, but I always wanted to learn the details of AES, and writing it in C wouldn't be nearly as fun. My thanks to Jeff Moser and Sam Trenholme for their excellent elucidations.

Note that this implementation is ECB mode, it doesn't include any decryption code, it computes rather than hard-codes the S-box, and it's probably vulnerable to side-channel attacks - so of course it's neither intended nor safe for production use.

{-# LANGUAGE NoMonomorphismRestriction #-}

import Control.Applicative (liftA2)
import Data.Bits (xor, shiftL, shiftR, (.|.), (.&.))
import Data.List (transpose, sortBy, foldl')
import Data.Ord (comparing)
import Data.Word (Word8)

encrypt input key = last ks `g` sRows (h t)
 where
  t  = foldl1 (g . f) $ init (k : tail ks)
  f  = transpose . map mix . transpose . sRows . h
  g  = zipWith $ zipWith xor
  h  = map $ map sub
  k  = input `g` head ks
  ks = expand key

mix [a,b,c,d] = [a', b', c', d']
 where
  a' = w ⊕ d ⊕ c ⊕ x ⊕ b
  b' = x ⊕ a ⊕ d ⊕ y ⊕ c
  c' = y ⊕ b ⊕ a ⊕ z ⊕ d
  d' = z ⊕ c ⊕ b ⊕ w ⊕ a 
  [w,x,y,z] = map fg [a,b,c,d]

fg b = b''
 where
  b'  = shiftL b 1
  b'' = ((b .&. 0x80) == 0x80) ? (b' ⊕ 0x1B, b')

sRows (w:ws) = w : zipWith f ws [1,2,3]
 where
  f w i = take 4 $ drop i $ cycle w

-----------------------------------------------------------
expand k = scanl f k [1, 2, 4, 8, 16, 32, 64, 128, 27, 54]
 where
  f n w = xpndE (transpose n) . xpndC . xpndB . xpndA $ xpnd0 w n

xpndE n [a,b,c,_] = transpose [a, b, c, zipWith xor c $ last n]

xpndC [a,b,c,d] = [a, b, zipWith xor b c, d]
xpndB [a,b,c,d] = [a, zipWith xor a b, c, d]
xpndA [a,b,c,d] = zipWith xor a d : [b,c,d]

xpnd0 rc ws = take 3 tW ++ [w']
 where
  w' = zipWith xor (map sub w) [rc, 0, 0, 0]
  tW = transpose ws
  w  = take 4 $ tail $ cycle $ last tW

----------------------------------------------------
sub w = get sbox (fromIntegral lo) $ fromIntegral hi
 where
  (hi, lo) = nibs w

nibs w    = (shiftR (w .&. 0xF0) 4, w .&. 0x0F)
(⊕)      = xor
p ? (a,b) = if p then a else b; infix 2 ?

get wss x y = (wss !! y) !! x

----------------------------------------------------
sbox = grid 16 $ map snd $ sortBy (comparing fst) $ sbx 1 1 []

sbx :: Word8 -> Word8 -> [(Word8, Word8)] -> [(Word8, Word8)]
sbx p q ws
  | length ws == 255 = (0, 0x63) : ws
  | otherwise = sbx p' r $ (p', xf ⊕ 0x63) : ws
 where
  p' = p ⊕ shiftL p 1 ⊕ ((p .&. 0x80 /= 0) ? (0x1B, 0))
  q1 = foldl' (liftA2 (.) xor shiftL) q [1, 2, 4]
  r  = q1 ⊕ ((q1 .&. 0x80 /= 0) ? (0x09, 0))
  xf = r  ⊕ rotl8 r 1 ⊕ rotl8 r 2 ⊕ rotl8 r 3 ⊕ rotl8 r 4

grid _ [] = []
grid n xs = take n xs : grid n (drop n xs)

rotl8 w n = (w `shiftL` n) .|. (w `shiftR` (8 - n))

butterfly curve

Temple Fay discovered this delightfully suggestive curve in 1989. It can be defined either parametrically or as a polar equation; I use the former method here.

I've also written a vector graphics game demo which uses this curve, and the Bézier curve methods in the below post, for object motion. The differences in plot density along the curve create natural-looking trails (e.g. for comets, or exhaust) without requiring any extra logic. (If you want to check out the game, be sure to use Chrome...)

{-# LANGUAGE NoMonomorphismRestriction #-}

import Data.Bits ((.|.), shiftL)
import Control.Arrow ((***), (&&&))
import Graphics.UI.SDL as SDL
import Graphics.UI.SDL.Primitives (pixel)

(xres, yres, zz) = (1050, 1050, round *** round)

main = withInit [InitVideo] $ do
  w <- setVideoMode xres yres 32 [NoFrame]
  enableEvent SDLMouseMotion False
  setCaption "Butterfly Curve" "Butterfly Curve"
  fillRect w (Just $ Rect 0 0 xres yres) $ Pixel 0
  plot w $ center $ scale curve
  loop w []
 where
  scale  = map ((150 *) *** (* 150))
  center = map (((xres / 2) +) *** (+ ((yres + 190) / 2)))

curve = map (f &&& (negate . g)) ts
 where
  ts  = [-999, -998.99.. 999]
  f t = sin t * (e ** cos t - 2 * cos (4*t) - h (t / 12))
  g t = cos t * (e ** cos t - 2 * cos (4*t) - h (t / 12))
  h   = foldr1 (.) $ replicate 5 sin
  e   = exp 1

plot w = mapM_ (f . zz)
 where
  f (x,y) = pixel w (fromIntegral x) (fromIntegral y) $ Pixel rgb
   where
    rgb   = rgb'  .|. shiftL (255 `div` (max x y `div` min x y)) 8
    rgb'  = rgb'' .|. shiftL (255 `div` (xres `div` x)) 16
    rgb'' = 0xFF  .|. shiftL (255 `div` (yres `div` y)) 24

loop w ps = do
  delay 128
  event <- pollEvent
  case event of
   KeyUp (Keysym SDLK_ESCAPE _ _) -> return ()
   _                              -> loop w ps

bézier curves

De Casteljau's algorithm is a fast, numerically stable way to rasterize Bézier curves (I was disappointed to see Wikipedia already gives succinct Haskell for it!). Clicking appends nodes to the curve's defining polygon; you can also drag any node to alter the curve.

{-# LANGUAGE NoMonomorphismRestriction #-}

import Control.Monad (when)
import Control.Arrow ((***))
import Data.List (findIndex)
import Graphics.UI.SDL as SDL hiding (init)
import Graphics.UI.SDL.Primitives (circle, line)

dim = 400

main = withInit [InitVideo] $ do
  w <- setVideoMode dim dim 32 []
  enableEvent SDLMouseMotion False
  setCaption "Bézier Curves" "Bézier Curves"
  loop w []

plot w ps = do
  fillRect w (Just $ Rect 0 0 dim dim) $ Pixel 0xFF222255
  mapM_ (f 2 0xFFFFFFFF . zz) [head ps, last ps]
  when b $ mapM_ (f 3 0x888888FF . zz) controls
  when b $ mapM_ (f 2 0xBBBBBBFF . zz) controls
 where
  f r c (x,y)   = circle w x y r $ Pixel c
  (b, controls) = (length ps > 2, tail $ init ps)

limn w [_] = SDL.flip w
limn w ((a,b):(x,y):ps) = do
  line w a b x y $ Pixel 0xFFFFFFFF
  limn w $ (x,y) : ps

loop w ps = do
  delay 128
  event <- pollEvent
  case event of
   KeyUp (Keysym SDLK_ESCAPE _ _) -> return ()
   MouseButtonDown x y _          -> click x y
   _                              -> loop w ps
 where
  click x y = let p = rr (x,y) in
   case findIndex ((10 >) . dist p) ps of
    Just i  -> drag w i ps
    Nothing -> do
      let ps' = p : ps
      plot w ps' >> SDL.flip w
      when (length ps' > 2) $ render w ps'
      loop w ps'

drag w i ps = do
  delay 16
  (x,y,_) <- getMouseState
  event   <- pollEvent
  let ps'  = swap i ps $ rr (x,y)
  plot w ps'
  when (length ps' > 2) $ render w ps'
  case event of
   MouseButtonUp x y _ -> loop w ps'
   _                   -> drag w i ps'

render w ps = limn w $ map zz curve
 where
  curve = map (casteljau ps) [0, 0.001.. 1]

casteljau [p] t = p
casteljau ps  t = casteljau ps' t
 where
  ps'             = zipWith (g t) ps $ tail ps
  g t (a,b) (c,d) = (f t a c, f t b d)
  f t a b         = (1 - t) * a + t * b

swap 0 ps p = p : tail ps
swap i ps p = take i ps ++ p : drop (i+1) ps

dist (a,b) (c,d) = sqrt $ (a-c)^2 + (b-d)^2

rr = fromIntegral *** fromIntegral
zz = round *** round

convex hulls

This code demonstrates the Graham scan, an O(n log n) method for finding the convex hull (smallest enclosing polygon) of a planar point set. It's a great example of exploiting order: it works by sorting the set by angle about a known extreme point, allowing the hull points to be found in linear time.

A variation on this algorithm, noted by A.M Andrew, sorts the set lexicographically and finds the upper and lower hull chains separately. This 'monotone chain' technique is often preferred, since it's easier to do robustly.

I've also written a couple JavaScript examples, the Graham scan and a slower (but constant-space and adaptable to higher dimensions) method called the Jarvis march.

{-# LANGUAGE NoMonomorphismRestriction #-}

import Control.Arrow ((***))
import Data.List (maximumBy, delete, sort, sortBy, unfoldr)
import Data.Ord (comparing)
import Graphics.UI.SDL as SDL
import Graphics.UI.SDL.Primitives (filledCircle, line)
import System.Random.Mersenne.Pure64 (newPureMT, randomDouble)

res = 250

main = withInit [InitVideo] $ do
  w  <- setVideoMode res res 32 []
  ps <- randPoints
  enableEvent SDLMouseMotion False
  setCaption "Graham Scan" "Graham Scan"
  fillRect w (Just $ Rect 0 0 res res) $ Pixel 0
  limn w $ map (round *** round) $ hull ps
  plot w ps
  pause

plot w ps = do
  mapM_ (f . (round *** round)) ps
  SDL.flip w
 where
  f (x,y) = filledCircle w x y 1 $ Pixel 0xFFFFFFFF

limn w ps = f $ ps ++ [head ps]
 where
  f [_] = return ()
  f ((a,b):(x,y):ps) = do
    line w a b x y $ Pixel 0xFF0000FF
    f $ (x,y) : ps

pause = do
  delay 128
  e <- pollEvent
  case e of
   KeyUp (Keysym SDLK_ESCAPE _ _) -> return ()
   _                              -> pause

hull qs = go (drop 2 ps) $ reverse $ take 2 ps
 where
  o  = bottomRightP qs
  ps = o : sortBy (ccw o) (delete o qs)

go [] qs = qs

go (p:ps) s@(a:b:qs)
  | ccw a b p /= GT = go (p:ps) $ b:qs
  | otherwise       = go ps $ p:s

ccw (ax, ay) (bx, by) (cx, cy)
  | d < 0 = LT
  | d > 0 = GT
  | True  = EQ
 where
  d = (bx - ax) * (cy - ay) - (by - ay) * (cx - ax)

bottomRightP = maximumBy (comparing snd) . sort

randPoints = fmap f newPureMT
 where
  f = uncurry zip . splitAt 20 . g
  g = map (* res) . unfoldr (Just . randomDouble)

more Π obscurantism

Here's an illustration of an approach to pi pointed out by Kevin Brown. Define f(n) as the nearest greater or equal multiple of n-1, then of n-2, etc (yielding OEIS sequence 2491). Then, inverting a result found by Duane Broline and Daniel Loeb, pi = n² / f(n).

But as you can see from the comment, the series converges very slowly!

main = print $ e**2 / f e e 1
 where
  e = 900000  -- yields 3.1416003...

f n 1 _ = n
f n k l = f n' (k - 1) $ n' / k
 where
  n' = head $ dropWhile (< n) $ map (* k) [l..]

segment intersection search

Here's a simple method I learned from Gareth Rees for finding the intersection (or parallel / collinear status) of two line segments; Rees credits Ronald Goldman, though I'd imagine the technique goes further back.

import Control.Arrow ((***))
import Control.Monad (liftM2)
import Data.Int (Int16)
import Data.List (unfoldr)
import Graphics.UI.SDL as SDL
import Graphics.UI.SDL.Primitives (line, circle)
import System.Random.Mersenne.Pure64 (pureMT, randomInt)

(xres, yres) = (1600, 900)
(red, white) = (0xFF0000FF, 0xFFFFFFFF)

main = withInit [InitVideo] $ do
  w <- setVideoMode xres yres 32 []
  enableEvent SDLMouseMotion False
  setCaption "Line Intersection" "Line Intersection"
  run w . take 15 $ filter f segments
 where
  f pq = dist pq < 1000
 
run w ss = do
  delay 256
  drawSegs w ss
  let ss' = map (vectorize . cst) ss
  drawPoints w $ liftM2 crossing ss' ss'
  SDL.flip w
  e <- pollEvent
  case e of
   KeyUp (Keysym SDLK_ESCAPE _ _) -> return ()
   _                              -> run w ss
  
drawSegs w = mapM_ f
 where
  f ((a,b), (c,d)) = line w a b c d $ Pixel red

drawPoints w = mapM_ (f . (round *** round))
 where
  f (x,y) = circle w x y 5 $ Pixel white 

-------------------------------------------------------

segments = h $ zip (f xres 13) $ f yres 23
 where
  f m      = map (fromIntegral . (`mod` m)) . g . pureMT
  g        = unfoldr $ fmap randomInt . Just
  h (p:ps) = (p, head ps) : h (tail ps)

dist (p,q) = sqrt $ (a-c)^2 + (b-d)^2
 where
  [a,b,c,d] = map f [fst p, snd p, fst q, snd q]
  f         = fromIntegral

crossing (p,r) (q,s)
  | f    = p .+. (t .^. r)  -- exactly one intersection
  | True = (-5,-5)          -- parallel, colinear, etc.
 where
  t   = (q .-. p) .*. s / (r .*. s)
  u   = (q .-. p) .*. r / (r .*. s)
  f   = (r .*. s) /= 0 && h t && h u
  h n = 0 <= n && n <= 1

(a,b) .+. (c,d) = (a+c, b+d)
(a,b) .-. (c,d) = (a-c, b-d)
n     .^. (x,y) = (n*x, n*y)   -- vector scaling
(a,b) .*. (c,d) =  a*d - b*c   -- 2D cross product

vectorize ((a,b), (c,d)) = ((a,b), (c-a, d-b))

cst (p,q) = (f p, f q)
 where
  f = fromIntegral *** fromIntegral

today in oblique approaches

If your standard library offers complex numbers and you're not in any hurry, you can't ask for much simpler (or more obscure) ways to compute pi than this.

import Data.Complex (magnitude, Complex((:+)))
 
pi = length $ takeWhile (< 2) $ map magnitude zs
 where
  zs = iterate (\z -> z^2 + ((-0.75) :+ 0.0000001)) 0

ikeda map

Continuing the theme of strange attractors, here's the well-known one embedded in the Ikeda map. The thumbnail at right shows the central 'vortex' of the attractor, and links to a larger viewport.

In these images, I've plotted the real and imaginary components along the x and y axes respectively. But the more popular way to visualize this attractor adds an extra parameter to the system and is expressed in trigonometric functions. Such adaptation of the code below yields these results.

{-# LANGUAGE NoMonomorphismRestriction #-}
 
import Data.Complex
import Graphics.UI.SDL as SDL
  
(xRes, yRes) = (1366, 768)
[a, b, k, p] = [0.85, 0.9, 0.4, 7.7]
 
main = withInit [InitVideo] $ do
  w <- setVideoMode xRes yRes 32 [NoFrame]
  s <- createRGBSurface [] 1 1 32 0 0 0 0
  fillRect s Nothing $ Pixel 0xFFFFFF
  enableEvent SDLMouseMotion False
  setCaption "Ikeda" "Ikeda"
  render w s $ ikeda $ 0 :+ 0
  SDL.flip w
  run w
 
run w = do
  e <- pollEvent
  delay 64
  case e of
   KeyUp (Keysym SDLK_ESCAPE _ _) -> return ()
   _                              -> run w
 
render w s = mapM_ draw
 where
  rect (x,y) = Just $ Rect (round x) (round y) 1 1
  g (x,y)    = (x + xRes / 24, y + yRes / 2.25)
  draw z     = blitSurface s Nothing w $ rect $ g p
   where
    p = (1050 * realPart z, 1050 * imagPart z)
 
ikeda = take 1500000 . filter g . iterate f
 where
  f z      = a + b * z * exp(i * (k-p) * (1 + abs z^2))
  i        = 0 :+ 1
  g (r:+i) = -35 < r && r < 35 && -18 < i && i < 18

hénon attractor

French astronomer Michel Hénon reported on this chaotic dynamical system in 1976.

{-# LANGUAGE NoMonomorphismRestriction #-}

import Graphics.UI.SDL as SDL
 
(xRes, yRes) = (1366, 768)
(a, b)       = (1.4, 0.3)

main = withInit [InitVideo] $ do
  w <- setVideoMode xRes yRes 32 [NoFrame]
  s <- createRGBSurface [] 1 1 32 0 0 0 0
  fillRect s Nothing $ Pixel 0xFFFFFF
  enableEvent SDLMouseMotion False
  setCaption "Hénon" "Hénon"
  render w s henon
  SDL.flip w
  run w

run w = do
  e <- pollEvent
  delay 64
  case e of
   KeyUp (Keysym SDLK_ESCAPE _ _) -> return ()
   _                              -> run w

render w s = mapM_ draw
 where
  rect (x,y) = Just $ Rect (round x) (round y) 1 1
  g (x,y)    = (500*x + xRes / 2, 900*y + yRes / 2)
  draw (x,y) = blitSurface s Nothing w $ rect $ g (x,y)

henon = take 99999 $ filter g $ iterate f (0,0)
 where
  f (x,y) = (1 - a*x^2 + y, b*x)
  g (x,y) = -1.5 < x && x < 1.5 && -0.45 < y && y < 0.45

langton's ant

For a round of code golf, I wrote this spare implementation of Chris Langton's remarkably simple universal computer. If you want amenities like pause, random starting pattern or even quit, check out this more complete version.

import Control.Arrow ((***))
import Control.Monad (when, join)
import Data.Set (insert, delete, member, empty)
import Graphics.UI.SDL as SDL

(xres, yres, sq) = (1366, 768, 4)

main = withInit [InitVideo] $ do
  w <- setVideoMode xres yres 32 [NoFrame]
  enableEvent SDLMouseMotion False
  setCaption "Langton's Ant" "Langton's Ant"
  run w [1..] p (0,1) empty
 where
  p = (xres `div` 2 `div` sq, yres `div` 2 `div` sq)

run w (n:ns) p v ps = do
  when (n `mod` 7 == 0) render
  run w ns (move p $ g v) (g v) $ f p ps
 where
  b          = member p ps
  (f,g)      = if b then (delete, fl) else (insert, fr)
  fr (x,y)   = if x == 0 then (-y,x) else (y,x)
  fl (x,y)   = if x == 0 then (y,x)  else (y,-x)
  move (x,y) = (x+) *** (+y)
  render     = do
    fillRect w (Just $ Rect 0 0 xres yres) $ Pixel 0
    mapM_ (draw w . join (***) (* sq)) ps
    SDL.flip w
 
draw w p = f p =<< g [SWSurface] sq sq 32 0 0 0 0
 where
  rect x y  = Just $ Rect x y sq sq
  g         = createRGBSurface
  f (x,y) s = do fillRect s (rect 0 0) $ Pixel 0xFFFFFF
                 blitSurface s (rect 0 0) w $ rect x y

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

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)

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)

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