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