gumowski-mira attractor

Here's an unusual chaotic attractor: the web doesn't seem to have much information on this one, except that it was invented to model particle trajectories in physics. A google image search for 'mira fractals' does turn up some pretty results though.

The system seems to give interesting results only when b is close to one. It behaves less chaotically when b > 1 is fixed, so you can actually animate it - click the thumbnail to see.

	import Control.Arrow ((***))
	import Graphics.UI.SDL as SDL
	(xres, yres) = (200, 200)
	main = withInit [InitVideo] $ do
	w <- setVideoMode xres yres 32 [NoFrame]
	enableEvent SDLMouseMotion False
	setCaption "Gumowski-Mira" "Gumowski-Mira"
	loop w p pss
	where
	p = (xres `div` 2, yres `div` 2)
	pss = map points [0.31, 0.3101.. 0.316]
	loop w p pss = do
	render w p $ map (round *** round) $ head pss
	e <- pollEvent
	case e of
	KeyUp (Keysym SDLK_ESCAPE _ _) -> return ()
	_ -> loop w p $ tail pss
	render w (x,y) ps = do
	s <- createRGBSurface [SWSurface] 2 2 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 a = take 650000 . map (h *** h) $ iterate f (12, 0)
	where
	f (x,y) = let x' = b*y + g x in (x', -x + g x')
	g x = a*x + (2 * (1 - a) * x^2) / (1 + x^2)
	(b,h) = (1.005, (* 4))

view raw gumowski-mira.hs hosted with ❤ by GitHub

128-bit AES electronic codebook

Rijndael (the core of AES) is an algorithm that might actually be shorter and simpler in C, but I was curious to see how natural I could make it look in Haskell. Thanks to Jeff Moser and Sam Trenholme for their clear elucidations.

Note that this code only does ECB mode; it computes rather than hard-codes the S-box; and it could be vulnerable to side-channel attacks. So enjoy reading it, but don't try to make a serious encryption app out of it. That kind of thing is best left to the professionals :)

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

view raw rijndael-AES.hs hosted with ❤ by GitHub

butterfly curve

Temple Fay discovered this complicated curve in 1989. It can be defined either parametrically or as a polar equation; I did it the former way.

One application I thought of for this is object motion in games: I tried it out by writing this little Canvas game, where the comets follow the curve's trajectory. The differences in plot density along the curve create natural-looking comet tails.

	{-# 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

view raw fay_butterfly_curve.hs hosted with ❤ by GitHub

bézier curves

De Casteljau's algorithm is a fast, numerically stable way to rasterize Bézier curves. This code implements an interactive demo: click to append nodes to the curve's defining polygon, and 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

view raw interactive_bezier.hs hosted with ❤ by GitHub

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 the algorithm, noted by A.M. Andrew, sorts the set lexicographically and finds the upper and lower hull chains separately. This 'monotone chain' version is often preferred, since it's easier to do robustly.

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

view raw graham-scan.hs hosted with ❤ by GitHub

more Π obscurantism

Here's an illustration of an approach to Π 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, Π = 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..]

view raw kevin_brown_PI.hs hosted with ❤ by GitHub