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

## Friday, June 3, 2016

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

## Monday, March 28, 2016

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

## Wednesday, February 17, 2016

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

## Saturday, February 13, 2016

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

## Tuesday, January 26, 2016

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

Labels:
cellular automata,
chaos,
computational universality

Subscribe to:
Posts (Atom)