Wednesday, December 13, 2017

gumowski-mira attractor

Another attractor! The internet doesn't seem to have much information on this one, except that it was developed to model particle trajectories (an image search for 'mira fractals' does turn up some nice results though).

This system seems to only give interesting results when b is close to one. It also shows some continuity when b > 1 is fixed, so you can actually animate it - click the image at right 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))

Monday, May 15, 2017

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. Thanks to Jeff Moser and Sam Trenholme for their clear 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))

Monday, March 6, 2017

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). To see the game properly, be sure the 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

Friday, February 24, 2017

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

Wednesday, February 22, 2017

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)

Wednesday, February 15, 2017

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 = n2 / 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..]