import Control.Arrow ((***)) import Control.Monad (join, void, when) import Data.Set (delete, empty, insert, member, toList) import Graphics.UI.SDL as SDL (xres, yres, sq) = (1600, 900, 2) main = withInit [InitVideo] $ do w <- setVideoMode xres yres 32 [NoFrame] enableEvent SDLMouseMotion False setCaption "Turmite" "Turmite" run w [1..] 0 p (0,1) empty where p = (xres `div` 2 `div` sq, yres `div` 2 `div` sq) run w (n:ns) s p v ps = do when (n `mod` 23 == 0) $ render w ps e <- pollEvent case e of KeyUp (Keysym SDLK_ESCAPE _ _) -> print n >> printout KeyUp (Keysym SDLK_SPACE _ _) -> pause w s p v ps _ -> continue where continue = go w ns s p v ps printout = void $ saveBMP w "output.bmp" pause w s p v ps = do delay 128 e <- pollEvent case e of KeyUp (Keysym SDLK_ESCAPE _ _) -> return () KeyUp (Keysym SDLK_SPACE _ _) -> run w [1..] s p v ps _ -> pause w s p v ps render w ps = do fillRect w (Just $ Rect 0 0 xres yres) $ Pixel 0 mapM_ (draw w . join (***) (* sq)) $ toList 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 go w ns s p v ps | s == 0 && not b = run w ns 0 (f p $ l' v) (l' v) qs | s == 0 && True = run w ns 1 (f p $ r' v) (r' v) qs | s == 1 && not b = run w ns 0 (f p $ r' v) (r' v) rs | otherwise = run w ns 1 (f p $ l' v) (l' v) rs where b = member p ps (qs, rs) = (insert p ps, delete p ps) f (x, y) = (x +) *** (+ y) r' (0, y) = (-y, 0) r' (x, y) = ( y, x) l' (0, y) = ( y, 0) l' (x, y) = ( y, -x)

# normal order

## Tuesday, March 13, 2018

### generalizing langton's ant

## Wednesday, December 13, 2017

### gumowski-mira attractor

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

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

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

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

*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

*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

^{2}/ 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..]

## Friday, June 3, 2016

### segment intersection search

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

## Monday, March 28, 2016

### today in oblique approaches

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

## Wednesday, February 17, 2016

### ikeda map

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

## Saturday, February 13, 2016

### hénon attractor

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

## Tuesday, January 26, 2016

### langton's ant

*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

## Wednesday, December 23, 2015

### connett circles

*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

*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

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