diff --git a/Address.hs b/Address.hs new file mode 100644 --- /dev/null +++ b/Address.hs @@ -0,0 +1,299 @@ +{- + + gmndl -- Mandelbrot Set explorer + Copyright (C) 2010,2011,2014 Claude Heiland-Allen + + This program is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License along + with this program; if not, write to the Free Software Foundation, Inc., + 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + +-} + +module Address(Address(..), angledInternalAddress, externalAngles, rayEnd, parameter, parse, pretty) where + +import Prelude hiding (isNaN) + +import Control.Monad (guard) +import Control.Monad.Identity (Identity()) +import Data.Char (digitToInt) +import Data.List (genericDrop, genericLength, genericTake, unfoldr) +import Data.Maybe (listToMaybe) +import Data.Ratio ((%), numerator, denominator) +import Data.Vec (NearZero()) +import Text.Parsec (ParsecT(), choice, digit, eof, many, many1, sepBy, string, try) +import Text.Parsec.Prim (runP) + +import Complex (Complex((:+)), mkPolar, Turbo) +import MuAtom (refineNucleus) + +isNaN x = not (x == x) + +double :: Rational -> Rational +double angle = wrap (2 * angle) + +wrap :: Rational -> Rational +wrap angle + | frac < 0 = frac + 1 + | otherwise = frac + where + _i :: Integer + (_i, frac) = properFraction angle + +data Knead = Zero | One | Star + deriving (Eq, Show) + +knead :: Rational -> [Knead] +knead angle + | angle == 0 || angle == 1 = [Star] + | otherwise = (++[Star]) . takeWhile (/= Star) . map k . iterate double $ angle + where + k a + | a `elem` [ angle / 2 , (angle + 1) / 2 ] = Star + | angle / 2 < a && a < (angle + 1) / 2 = One + | a < angle / 2 || (angle + 1) / 2 < a = Zero + +period :: Rational -> Integer +period angle = genericLength (knead angle) + +internalAddress :: [Knead] -> [Integer] +internalAddress v = iA 1 [Star] + where + iA sk vk + | sk == genericLength v = [ ] + | otherwise = sk' : iA sk' vk' + where + sk' = (+) 1 . genericLength . takeWhile id $ zipWith (==) (cycle v) (cycle vk) + vk' = genericTake sk' v + +orbit :: Eq a => (a -> Maybe a) -> a -> [a] +orbit f x = x : unfoldr (fmap both . f) x + where + both z = (z, z) + +rho :: [Knead] -> Integer -> Maybe Integer +rho v r = listToMaybe . concat $ zipWith3 f [r+1 .. 1000000] (zipWith (flip const) v (genericDrop r (cycle v))) v + where + f k a b + | a /= b = [k] + | otherwise = [] + +denominators :: [Knead] -> [Integer] +denominators v = zipWith f a (tail a) + where + a = internalAddress v + f sk sk1 + | sk `elem` orbit (rho v) r = (sk1 - r) `div` sk + 1 + | otherwise = (sk1 - r) `div` sk + 2 + where + r | sk1 `mod` sk == 0 = sk + | otherwise = sk1 `mod` sk + +numerators :: Rational -> [Integer] -> [Integer] -> [Integer] +numerators angle = zipWith f + where + f qk sk = genericLength . filter (<= angle) $ [ wrap $ 2^(i * sk) * angle | i <- [0 .. qk - 2] ] + +data Address = P Integer | S Integer Rational Address + deriving (Eq, Ord, Show) + +angledInternalAddress :: Rational -> Address +angledInternalAddress angle = foldr (\(s, pq) a -> S s pq a) (P (last ss)) (zip ss rs) + where + rs = zipWith (%) ns ds + ns = numerators angle ds ss + ds = denominators ks + ss = internalAddress ks + ks = knead angle + +externalAngles :: Address -> Maybe (Rational, Rational) +externalAngles = externalAngles' 1 (0, 1) + +externalAngles' :: Integer -> (Rational, Rational) -> Address -> Maybe (Rational, Rational) +externalAngles' p0 lohi a0@(P p) + | p0 /= p = case wakees lohi p of + [lh] -> externalAngles' p lh a0 + _ -> Nothing + | otherwise = Just lohi +externalAngles' p0 lohi a0@(S p r a) + | p0 /= p = case wakees lohi p of + [lh] -> externalAngles' p lh a0 + _ -> Nothing + | otherwise = do + let num = numerator r + den = denominator r + q = p * den + ws = wakees lohi q + nums = [ num' | num' <- [ 1.. den - 1 ], let r' = num' % den, denominator r' == den ] + nws, nnums :: Integer + nws = genericLength ws + nnums = genericLength nums + guard (nws == nnums) + i <- genericElemIndex num nums + lh <- safeGenericIndex ws (i :: Integer) + externalAngles' q lh a + +wakees :: (Rational, Rational) -> Integer -> [(Rational, Rational)] +wakees (lo, hi) q = + let gaps (l, h) n + | n == 0 = [(l, h)] + | n > 0 = let gs = gaps (l, h) (n - 1) + cs = candidates n gs + in accumulate cs gs + candidates n gs = + let den = 2 ^ n - 1 + in [ r + | (l, h) <- gs + , num <- [ ceiling (l * fromInteger den) + .. floor (h * fromInteger den) ] + , let r = num % den + , l < r, r < h + , period r == n + ] + accumulate [] ws = ws + accumulate (l : h : lhs) ws = + let (ls, ms@((ml, _):_)) = break (l `inside`) ws + (_s, (_, rh):rs) = break (h `inside`) ms + in ls ++ [(ml, l)] ++ accumulate lhs ((h, rh) : rs) + inside x (l, h) = l < x && x < h + in chunk2 . candidates q . gaps (lo, hi) $ (q - 1) + +chunk2 :: [t] -> [(t, t)] +chunk2 [] = [] +chunk2 (x:y:zs) = (x, y) : chunk2 zs + +genericElemIndex :: (Eq a, Integral b) => a -> [a] -> Maybe b +genericElemIndex _ [] = Nothing +genericElemIndex e (f:fs) + | e == f = Just 0 + | otherwise = (1 +) `fmap` genericElemIndex e fs + +safeGenericIndex :: Integral b => [a] -> b -> Maybe a +safeGenericIndex [] _ = Nothing +safeGenericIndex (x:xs) i + | i < 0 = Nothing + | i > 0 = safeGenericIndex xs (i - 1) + | otherwise = Just x + +safeLast :: [a] -> Maybe a +safeLast [] = Nothing +safeLast xs = Just (last xs) + +radius :: (Real r, Floating r) => r +radius = 2 ** 24 + +sharpness :: Int +sharpness = 4 + +limit :: Int +limit = 64 + +distance :: Int +distance = 64 + +ray :: (Real r, Floating r, Turbo r) => Rational -> [Complex r] +ray angle = map fst . iterate (step angle) $ (mkPolar radius (2 * pi * fromRational angle), (0, 0)) + +step :: (Real r, Floating r, Turbo r) => Rational -> (Complex r, (Int, Int)) -> (Complex r, (Int, Int)) +step angle (c, (k0, j0)) + | j > sharpness = step angle (c, (k0 + 1, 0)) + | otherwise = (c', (k0, j0 + 1)) + where + k = k0 + 1 + j = j0 + 1 + m = (k - 1) * sharpness + j + r = radius ** ((1/2) ** (fromIntegral m / fromIntegral sharpness)) + t = mkPolar (r ** (2 ** fromIntegral k0)) ((2 ** fromIntegral k0) * 2 * pi * fromRational angle) + c' = iterate n c !! limit + n z = z - (cc - t) / dd + where + (cc, dd) = ncnd k + ncnd 1 = (z, 1) + ncnd i = let (nc, nd) = ncnd (i - 1) in (nc * nc + z, 2 * nc * nd + 1) + +rayEnd :: (Real r, Floating r, Turbo r) => Rational -> Maybe (Complex r) +rayEnd = safeLast . takeWhile (\(r:+i) -> not (isNaN r || isNaN i)) . take (sharpness * distance) . ray + +parameter :: (NearZero r, Real r, Floating r, Turbo r) => Address -> Maybe (r, r, r) +parameter a = do + (lo, hi) <- externalAngles a + c1 <- rayEnd lo + c2 <- rayEnd hi + let c = 0.5 * (c1 + c2) + return $ refineNucleus (addressPeriod a) c + +addressPeriod :: Address -> Integer +addressPeriod (P p) = p +addressPeriod (S _ _ a) = addressPeriod a + +parse :: String -> Maybe Address +parse s = case runP parser () "" s of + Left _ -> Nothing + Right a -> Just a + +data Token = Number Integer | Fraction Integer Integer + +type Parse t = ParsecT String () Identity t + +parser :: Parse Address +parser = do + ts <- pTokens + accum 1 ts + where + accum p [] = return $ P p + accum _ [Number n] = return $ P n + accum _ (Number n : ts@(Number _ : _)) = do + a <- accum n ts + return $ S n (1%2) a + accum _ (Number n : Fraction t b : ts) = do + a <- accum (n * b) ts + return $ S n (t%b) a + accum p (Fraction t b : ts) = do + a <- accum (p * b) ts + return $ S p (t % b) a + +pTokens :: Parse [Token] +pTokens = do + _ <- pOptionalSpace + ts <- pToken `sepBy` pSpace + _ <- pOptionalSpace + eof + return ts + +pToken :: Parse Token +pToken = choice [ try pFraction, pNumber ] + +pFraction :: Parse Token +pFraction = do + Number top <- pNumber + _ <- pOptionalSpace + _ <- string "/" + _ <- pOptionalSpace + Number bottom <- pNumber + guard $ top < bottom + return $ Fraction top bottom + +pNumber :: Parse Token +pNumber = do + n <- foldl (\x y -> 10 * x + y) 0 `fmap` map (toInteger . digitToInt) `fmap` many1 digit + guard $ 0 < n + return $ Number n + +pSpace :: Parse [String] +pSpace = many1 (string " ") + +pOptionalSpace :: Parse [String] +pOptionalSpace = many (string " ") + +pretty :: Address -> String +pretty (P p) = show p +pretty (S p r a) = show p ++ " " ++ show (numerator r) ++ "/" ++ show (denominator r) ++ " " ++ pretty a diff --git a/Calculate.hs b/Calculate.hs new file mode 100644 --- /dev/null +++ b/Calculate.hs @@ -0,0 +1,286 @@ +{-# LANGUAGE BangPatterns #-} + +{- + + gmndl -- Mandelbrot Set explorer + Copyright (C) 2010,2011,2014 Claude Heiland-Allen + + This program is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License along + with this program; if not, write to the Free Software Foundation, Inc., + 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + +-} + +module Calculate (convert, renderer) where + +-- simple helpers +import Control.Monad (when) + +-- concurrent renderer with capability-specific scheduling +import Control.Concurrent (MVar, newEmptyMVar, takeMVar, tryTakeMVar, tryPutMVar, threadDelay, forkIO, killThread) +import GHC.Conc (forkOn, numCapabilities) + +-- each worker uses a mutable unboxed array of Bool to know which pixels +-- it has already started to render, to avoid pointless work duplication +import Data.Array.IO (IOUArray, newArray, readArray, writeArray, inRange) + +-- each worker thread keeps a queue of pixels that it needs to render or +-- to continue rendering later +import Data.PriorityQueue (PriorityQueue, newPriorityQueue, enqueue, enqueueBatch, dequeue) + +-- poking bytes into memory is dirty, but it's quick and allows use of +-- other fast functions like memset and easy integration with OpenGL +import Foreign (Word8) + +-- higher precision arithmetic using libqd +import Numeric.QD.DoubleDouble (DoubleDouble()) +import Numeric.QD.QuadDouble (QuadDouble()) + +import Complex (Complex((:+)), Turbo(sqr, twice), convert) + +-- some type aliases to shorten things +type B = Word8 +type N = Int +type R = Double + +{- +-- colour space conversion from HSV [0..1] to RGB [0..1] +-- HSV looks quite 'chemical' to my eyes, need to investigate something +-- better to make it feel more 'natural' +hsv2rgb :: R -> R -> R -> (R, R, R) +hsv2rgb !h !s !v + | s == 0 = (v, v, v) + | h == 1 = hsv2rgb 0 s v + | otherwise = + let !i = floor (h * 6) `mod` 6 :: N + !f = (h * 6) - fromIntegral i + !p = v * (1 - s) + !q = v * (1 - s * f) + !t = v * (1 - s * (1 - f)) + in case i of + 0 -> (v, t, p) + 1 -> (q, v, p) + 2 -> (p, v, t) + 3 -> (p, q, v) + 4 -> (t, p, v) + 5 -> (v, p, q) + _ -> (0, 0, 0) +-} + +hsv2rgb' :: R -> R -> R -> (R, R, R) +hsv2rgb' !h !s !l = + let !a = 2 * pi * h + !ca = cos a + !sa = sin a + !y = l / 2 + !u = s / 2 * ca + !v = s / 2 * sa + !r = y + 1.407 * u + !g = y - 0.677 * u - 0.236 * v + !b = y + 1.848 * v + in (r, g, b) + +-- compute RGB [0..255] bytes from the results of the complex iterations +-- don't need very high precision for this, as spatial aliasing will be +-- much more of a problem in intricate regions of the fractal +colour :: Complex Double -> Complex Double -> N -> (B, B, B) +colour !(zr:+zi) !(dzr:+dzi) !n = + let -- micro-optimization - there is no log2 function + !il2 = 1 / log 2 + !zd2 = sqr zr + sqr zi + !dzd2 = sqr dzr + sqr dzi + -- normalized escape time + !d = (fromIntegral n :: R) - log (log zd2 / log escapeR2) * il2 + !dwell = fromIntegral (floor d :: N) + -- final angle of the iterate + !finala = atan2 zi zr + -- distance estimate + !de = (log zd2 * il2) * sqrt (zd2 / dzd2) + !dscale = -log de * il2 + -- HSV is based on escape time, distance estimate, and angle + !hue = log ((- log de * il2) `max` 1) * il2 / 16 + log d * il2 + !saturation = 0 `max` (log d * il2 / 8) `min` 1 + !value = 0 `max` (1 - dscale / 256) `min` 1 + !h = hue - fromIntegral (floor hue :: N) + -- adjust saturation to give concentric striped pattern + !k = dwell / 2 + !satf = if k - fromIntegral (floor k :: N) >= (0.5 :: R) then 0.9 else 1 + -- adjust value to give tiled pattern + !valf = if finala < 0 then 0.9 else 1 + -- convert to RGB + (!r, !g, !b) = hsv2rgb' h (satf * saturation) (valf * value) + -- convert to bytes + !rr = floor $ 0 `max` (255 * r) `min` 255 + !gg = floor $ 0 `max` (255 * g) `min` 255 + !bb = floor $ 0 `max` (255 * b) `min` 255 + in (rr, gg, bb) + +-- a Job stores a pixel undergoing iterations +data Job c = Job !N !N !(Complex c) !(Complex c) !(Complex c) !N + +-- the priority of a Job is how many iterations have been computed: +-- so 'fresher' pixels drop to the front of the queue in the hope of +-- avoiding too much work iterating pixels that will never escape +priority :: Job c -> N +priority !(Job _ _ _ _ _ n) = n + +-- add a job to a work queue, taking care not to duplicate work +-- there is no race condition here as each worker has its own queue +addJob :: (Real c, Floating c, Turbo c) => N -> N -> Complex c -> c -> PriorityQueue IO (Job c) -> IOUArray (N,N) Bool -> N -> N -> IO () +addJob !w !h !c !zradius' todo sync !i !j = do + already <- readArray sync (j, i) + when (not already) $ do + writeArray sync (j, i) True + enqueue todo $! Job i j (coords w h c zradius' i j) 0 0 0 + +-- spawns a new batch of workers to render an image +-- returns an action that stops the rendering +renderer' :: (Turbo c, Real c, Floating c) => MVar () -> ((N,N),(N,N)) -> (N -> N -> B -> B -> B -> IO ()) -> Complex c -> c -> IO (IO ()) +renderer' done rng output !c !zradius' = do + wdog <- newEmptyMVar + workerts <- mapM (\w -> forkOn w $ worker wdog rng c zradius' output w) [ 0 .. workers - 1 ] + watcher <- forkIO $ do + () <- takeMVar wdog + let loop = do + threadDelay 10000000 -- 10 seconds + m <- tryTakeMVar wdog + case m of + Nothing -> mapM_ killThread workerts >> tryPutMVar done () >> return () + Just () -> loop + loop + return $ killThread watcher >> mapM_ killThread workerts + +-- compute the Complex 'c' coordinate for a pixel in the image +coords :: (Real c, Floating c, Turbo c) => N -> N -> Complex c -> c -> N -> N -> Complex c +coords !w !h !c !zradius' !i !j = c + ( (fromIntegral (i - w`div`2) * k) + :+(fromIntegral (h`div`2 - j) * k)) + where !k = zradius' / (fromIntegral $ (w `div` 2) `min` (h `div` 2)) + +-- the worker thread enqueues its border and starts computing iterations +worker :: (Turbo c, Real c, Floating c) => MVar () -> ((N,N),(N,N)) -> Complex c -> c -> (N -> N -> B -> B -> B -> IO ()) -> N -> IO () +worker wdog rng@((y0,x0),(y1,x1)) !c !zradius' output !me = do + sync <- newArray rng False + queue <- newPriorityQueue priority + let addJ = addJob w h c zradius' queue sync + js = filter mine (border w h) + w = x1 - x0 + 1 + h = y1 - y0 + 1 + mapM_ (flip (writeArray sync) True) js + enqueueBatch queue (map (\(j,i) -> Job i j (coords w h c zradius' i j) 0 0 0) js) + compute wdog rng addJ output queue + where mine (_, i) = i `mod` workers == me -- another dependency on spread + +-- the compute engine pulls pixels from the queue until there are no +-- more, and calculates a batch of iterations for each +compute :: (Turbo c, Real c, Floating c) => MVar () -> ((N,N),(N,N)) -> (N -> N -> IO ()) -> (N -> N -> B -> B -> B -> IO ()) -> PriorityQueue IO (Job c) -> IO () +compute wdog rng addJ output queue = do + mjob <- dequeue queue + case mjob of + Just (Job i j c z dz n) -> do + let -- called when the pixel escapes + done' !(zr:+zi) !(dzr:+dzi) !n' = {-# SCC "done" #-} do + _ <- tryPutMVar wdog () + let (r, g, b) = colour (convert zr :+ convert zi) (convert dzr :+ convert dzi) n' + output i j r g b + -- a wavefront of computation spreads to neighbouring pixels + sequence_ + [ addJ x y + | u <- spreadX + , v <- spreadY + , let x = i + u + , let y = j + v + , inRange rng (y, x) + ] + -- called when the pixel doesn't escape yet + todo' !z' !dz' !n' = {-# SCC "todo" #-} {- output i j 255 0 0 >> -} enqueue queue $! Job i j c z' dz' n' + calculate c limit z dz n done' todo' + compute wdog rng addJ output queue + Nothing -> return () -- no pixels left to render, so finish quietly + +-- the raw z->z^2+c calculation engine +-- also computes the derivative for distance estimation calculations +-- this function is crucial for speed, too much allocation will slooow +-- everything down severely +calculate :: (Turbo c, Real c, Floating c) => Complex c -> N -> Complex c -> Complex c -> N -> (Complex c -> Complex c -> N -> IO ()) -> (Complex c -> Complex c -> N -> IO ()) -> IO () +calculate !c !m0 !z0 !dz0 !n0 done todo = go m0 z0 dz0 n0 + where + go !m !z@(zr:+zi) !dz !n + | not (sqr zr + sqr zi < er2) = done z dz n + | m <= 0 = todo z dz n + | otherwise = go (m - 1) (sqr z + c) (let !zdz = z * dz in twice zdz + 1) (n + 1) + !er2 = convert escapeR2 + +-- dispatch to different instances of renderer depending on required precision +-- if zoom is low, single precision Float is ok, but as soon as pixel spacing +-- gets really small, it's necessary to increase it +-- it's probably not even worth using Float - worth benchmarking this and +-- also the DD and QD types (which cause a massively noticeable slowdown) +renderer :: (Real c, Floating c) => MVar () -> ((N,N),(N,N)) -> (N -> N -> B -> B -> B -> IO ()) -> Complex c -> c -> IO (IO ()) +renderer done rng output !c !zradius' + | zoom' < 20 = {-# SCC "rF" #-} renderer' done rng output (f c :: Complex Float ) (g zradius') + | zoom' < 50 = {-# SCC "rD" #-} renderer' done rng output (f c :: Complex Double ) (g zradius') + | zoom' < 100 = {-# SCC "rDD" #-} renderer' done rng output (f c :: Complex DoubleDouble) (g zradius') + | otherwise = {-# SCC "rQD" #-} renderer' done rng output (f c :: Complex QuadDouble ) (g zradius') + where f !(cr :+ ci) = convert cr :+ convert ci + g !x = convert x + zoom' = - logBase 2 (zradius' / (fromIntegral $ w `min` h)) + ((x0,y0), (x1, y1)) = rng + w = x1 - x0 + 1 + h = y1 - y0 + 1 + +-- start rendering pixels from the edge of the image +-- the Mandelbrot Set and its complement are both simply-connected +-- discounting spatial aliasing any point inside the boundary that is +-- in the complement is 'reachable' from a point on the boundary that +-- is also in the complement - probably some heavy math involved to +-- prove this though +-- note: this implicitly depends on the spread values below - it's +-- necessary for each interlaced subimage (one per worker) to have +-- at least a one pixel deep border +border :: N -> N -> [(N, N)] +border !w !h = concat $ + [ [ (j, i) | i <- [ 0 .. w - 1 ], j <- [ 0 ] ] + , [ (j, i) | j <- [ 0 .. h - 1 ], i <- [ 0 .. workers - 1 ] ] + , [ (j, i) | j <- [ 0 .. h - 1 ], i <- [ w - workers .. w - 1 ] ] + , [ (j, i) | i <- [ 0 .. w - 1 ], j <- [ h - 1 ] ] + ] + +-- which neighbours to activate once a pixel has escaped +-- there are essentially two choices, with x<->y swapped +-- choose greater X spread because images are often wider than tall +-- other schemes wherein the spread is split in both directions +-- might benefit appearance with large worker count, but too complicated +spreadX, spreadY :: [ N ] +spreadX = [ -workers, 0, workers ] +spreadY = [ -1, 0, 1 ] + +-- number of worker threads +-- use as many worker threads as capabilities, with the workers +-- distributed 1-1 onto capabilities to maximize CPU utilization +workers :: N +workers = numCapabilities + +-- iteration limit per pixel +-- at most this many iterations are performed on each pixel before it +-- is shunted to the back of the work queue +-- this should be tuneable to balance display updates against overheads +limit :: N +limit = (2^(13::N)-1) + +-- escape radius for fractal iteration calculations +-- once the complex iterate exceeds this, it's never coming back +-- theoretically escapeR = 2 would work +-- but higher values like this give a significantly smoother picture +escapeR, escapeR2 :: R +escapeR = 65536 +escapeR2 = escapeR * escapeR diff --git a/Complex.hs b/Complex.hs new file mode 100644 --- /dev/null +++ b/Complex.hs @@ -0,0 +1,152 @@ +{-# LANGUAGE BangPatterns, FlexibleContexts #-} + +{- + + gmndl -- Mandelbrot Set explorer + Copyright (C) 2010,2011,2014 Claude Heiland-Allen + + This program is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License along + with this program; if not, write to the Free Software Foundation, Inc., + 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + +-} + +module Complex where + +import Prelude hiding (atan2) + +import Foreign.C (CDouble) + +-- higher precision arithmetic using libqd +import Numeric.QD.DoubleDouble (DoubleDouble(DoubleDouble)) +import Numeric.QD.QuadDouble (QuadDouble(QuadDouble)) +import qualified Numeric.QD.DoubleDouble as DD +import qualified Numeric.QD.QuadDouble as QD +import Numeric.AD.Mode.Reverse (Reverse) +import Data.Reflection (Reifies) +import Numeric.AD.Internal.Reverse (Tape) + +-- ugly! but the default realToFrac :: (C)Double -> (C)Double is slooow +import Unsafe.Coerce (unsafeCoerce) + +-- don't look! this is really really ugly, and should be benchmarked +-- to see how really necessary it is, or at least made into a type class +convert :: (Real a, Fractional b) => a -> b +convert = realToFrac +{-# NOINLINE convert #-} +convertDouble2CDouble :: Double -> CDouble +convertDouble2CDouble !x = unsafeCoerce x +convertCDouble2Double :: CDouble -> Double +convertCDouble2Double !x = unsafeCoerce x +convertDouble2DoubleDouble :: Double -> DoubleDouble +convertDouble2DoubleDouble !x = convertCDouble2DoubleDouble . convertDouble2CDouble $ x +convertCDouble2DoubleDouble :: CDouble -> DoubleDouble +convertCDouble2DoubleDouble !x = DoubleDouble x 0 +convertDoubleDouble2Double :: DoubleDouble -> Double +convertDoubleDouble2Double !(DoubleDouble x _) = convertCDouble2Double x +convertDoubleDouble2CDouble :: DoubleDouble -> CDouble +convertDoubleDouble2CDouble !(DoubleDouble x _) = x +{-# RULES "convert/Double2CDouble" convert = convertDouble2CDouble #-} +{-# RULES "convert/CDouble2Double" convert = convertCDouble2Double #-} +{-# RULES "convert/Double2DoubleDouble" convert = convertDouble2DoubleDouble #-} +{-# RULES "convert/CDouble2DoubleDouble" convert = convertCDouble2DoubleDouble #-} +{-# RULES "convert/DoubleDouble2Double" convert = convertDoubleDouble2Double #-} +{-# RULES "convert/DoubleDouble2CDouble" convert = convertDoubleDouble2CDouble #-} + +{- +-- this is ugly too: can't use Data.Complex because the qd bindings do +-- not implement some low-level functions properly, leading to obscure +-- crashes inside various Data.Complex functions... +data Complex c = {-# UNPACK #-} !c :+ {-# UNPACK #-} !c deriving (Read, Show, Eq) + +-- complex number arithmetic, with extra strictness and cost-centres +instance Num c => Num (Complex c) where + (!(a :+ b)) + (!(c :+ d)) = {-# SCC "C+" #-} ((a + c) :+ (b + d)) + (!(a :+ b)) - (!(c :+ d)) = {-# SCC "C-" #-} ((a - c) :+ (b - d)) + (!(a :+ b)) * (!(c :+ d)) = {-# SCC "C*" #-} ((a * c - b * d) :+ (a * d + b * c)) + negate !(a :+ b) = (-a) :+ (-b) + abs x = error $ "Complex.abs: " ++ show x + signum x = error $ "Complex.signum: " ++ show x + fromInteger !x = fromInteger x :+ 0 +-} + +-- an extra class for some operations that can be made faster for things +-- like DoubleDouble: probably should have given this a better name +class Num c => Turbo c where + sqr :: c -> c + sqr !x = x * x + twice :: c -> c + twice !x = x + x + +-- the default methods are fine for simple primitive types... +instance Turbo Float where +instance Turbo Double where +instance Turbo CDouble where + +-- ...and complex numbers +instance (Real c, Floating c, Turbo c) => Turbo (Complex c) where + sqr !(r :+ i) = (sqr r - sqr i) :+ (twice (r * i)) + twice !(r :+ i) = (twice r) :+ (twice i) + +-- use the specific implementations for the higher precision types +instance Turbo DoubleDouble where + sqr !x = DD.sqr x + twice !(DoubleDouble a b) = DoubleDouble (twice a) (twice b) + +instance Turbo QuadDouble where + sqr !x = QD.sqr x + twice !(QuadDouble a b c d) = QuadDouble (twice a) (twice b) (twice c) (twice d) + +instance (Reifies s Tape, Num r) => Turbo (Reverse s r) where + + +data Complex r = !r :+ !r + deriving (Read, Show, Eq) + +instance (Real r, Floating r, Turbo r) => Num (Complex r) where + (a :+ b) + (x :+ y) = (a + x) :+ (b + y) + (a :+ b) - (x :+ y) = (a - x) :+ (b - y) + (a :+ b) * (x :+ y) = (a * x - b * y) :+ (a * y + b * x) + negate (a :+ b) = negate a :+ negate b + abs c = magnitude c :+ 0 + signum = normalize + fromInteger n = fromInteger n :+ 0 + +instance (Real r, Floating r, Turbo r) => Fractional (Complex r) where + (a:+b) / (c:+d) = ((a * c + b * d)/m2) :+ ((b * c - a * d)/m2) where m2 = sqr c + sqr d + fromRational r = fromRational r :+ 0 + +magnitude :: (Real c, Floating c, Turbo c) => Complex c -> c +magnitude (re:+im) = sqrt $ sqr re + sqr im + +cis :: (Real c, Floating c) => c -> Complex c +cis a = cos a :+ sin a + +mkPolar :: (Real c, Floating c) => c -> c -> Complex c +mkPolar r a = (r * cos a) :+ (r * sin a) + +phase :: (Real c, Floating c) => Complex c -> c +phase (re:+im) = atan2 im re + +normalize :: (Real c, Floating c, Turbo c) => Complex c -> Complex c +normalize z@(re:+im) = let m = magnitude z in (re / m) :+ (im / m) + +atan2 :: (Real c, Floating c) => c -> c -> c +atan2 y x + | x > 0 = atan (y/x) + | x == 0 && y > 0 = pi/2 + | x < 0 && y > 0 = pi + atan (y/x) + | x <= 0 && y < 0 = -atan2 (-y) x + | y == 0 && x < 0 = pi -- must be after the previous test on zero y + | x == 0 && y == 0 = y -- must be after the other double zero tests + | otherwise = x + y -- x or y is a NaN, return a NaN (via +) diff --git a/Image.hs b/Image.hs new file mode 100644 --- /dev/null +++ b/Image.hs @@ -0,0 +1,135 @@ +{-# LANGUAGE ForeignFunctionInterface #-} + +{- + + gmndl -- Mandelbrot Set explorer + Copyright (C) 2010,2011,2014 Claude Heiland-Allen + + This program is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License along + with this program; if not, write to the Free Software Foundation, Inc., + 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + +-} + +module Image (Image(), new, clear, plot, upload, with, draw, putPPM, hPutPPM) where + +-- poking bytes into memory is dirty, but it's quick and allows use of +-- other fast functions like memset and easy integration with OpenGL +import Foreign (castPtr, mallocBytes, nullPtr, plusPtr, pokeByteOff, Ptr, Word8) +import Foreign.C (CInt(..), CSize(..)) +import qualified Graphics.Rendering.OpenGL as GL +import Graphics.Rendering.OpenGL (($=)) +import Graphics.Rendering.OpenGL.Raw (glGenerateMipmap, gl_TEXTURE_2D) +import System.IO (Handle, hPutBuf, hPutStr, stdout) + +data Image + = Image + { pixels :: Ptr Word8 + , width :: Int + , height :: Int + , size :: Int + , texture :: GL.TextureObject + } + +-- these images are RGB only +channels :: Int +channels = 3 + +-- create a new image +-- note: this needs an OpenGL context +new :: Int -> Int -> IO Image +new w h = do + p <- mallocBytes $ w * h * channels + [t] <- GL.genObjectNames 1 + let s = roundUp2 (w `max` h) + i = Image{ pixels = p, width = w, height = h, size = s, texture = t } + dim = GL.TextureSize2D (fromIntegral s) (fromIntegral s) + img = GL.PixelData GL.RGB GL.UnsignedByte nullPtr + with i $ do + GL.texImage2D GL.Texture2D GL.NoProxy 0 GL.RGB' dim 0 img + GL.textureFilter GL.Texture2D $= ((GL.Linear', Just GL.Linear'), GL.Linear') + GL.textureWrapMode GL.Texture2D GL.S $= (GL.Repeated, GL.ClampToEdge) + GL.textureWrapMode GL.Texture2D GL.T $= (GL.Repeated, GL.ClampToEdge) + clear i + upload i + return i + +-- clear with white +clear :: Image -> IO () +clear i = do + let bytes = width i * height i * channels + _ <- memset (castPtr (pixels i)) 255 (fromIntegral bytes) + return () + +-- plot a pixel +plot :: Image -> Int -> Int -> Word8 -> Word8 -> Word8 -> IO () +plot i x y r g b = do + let p = pixels i `plusPtr` ((y * width i + x) * channels) + pokeByteOff p 0 r + pokeByteOff p 1 g + pokeByteOff p 2 b + +-- upload an image to its texture +upload :: Image -> IO () +upload i = do + let pos = GL.TexturePosition2D 0 0 + dim = GL.TextureSize2D (fromIntegral $ width i) (fromIntegral $ height i) + img = GL.PixelData GL.RGB GL.UnsignedByte (pixels i) + with i $ do + GL.texSubImage2D GL.Texture2D 0 pos dim img + glGenerateMipmap gl_TEXTURE_2D + +-- use a texture +-- FIXME TODO preserve the previous texture binding instead of clearing +with :: Image -> IO a -> IO a +with i act = do + GL.textureBinding GL.Texture2D $= Just (texture i) + r <- act + GL.textureBinding GL.Texture2D $= Nothing + return r + +-- draw textured unit quad +draw :: Image -> IO () +draw i = do + let v :: GL.GLfloat -> GL.GLfloat -> GL.GLfloat -> GL.GLfloat -> IO () + v tx ty vx vy = GL.texCoord (GL.TexCoord2 tx ty) >> GL.vertex (GL.Vertex2 vx vy) + sx = fromIntegral (width i) / fromIntegral (size i) + sy = fromIntegral (height i) / fromIntegral (size i) + with i $ GL.renderPrimitive GL.Quads $ do + v 0 sy 0 0 + v 0 0 0 1 + v sx 0 1 1 + v sx sy 1 0 + +-- save as PPM +putPPM :: Image -> IO () +putPPM = hPutPPM stdout + +hPutPPM :: Handle -> Image -> IO () +hPutPPM h i = do + hPutStr h ("P6\n" ++ show (width i) ++ " " ++ show (height i) ++ " 255\n") + hPutBuf h (pixels i) (width i * height i * channels) + +-- round up to nearest power of two +-- this will probably explode when n gets large, but it's only used +-- for OpenGL texture dimensions so you'll run out of memory first +roundUp2 :: Int -> Int +roundUp2 n = head . dropWhile (< n) . iterate (2*) $ 1 + +-- import standard C library memset for clearing images efficiently +-- previous implementation used pokeArray ... (replicate ...) ... +-- which had a nasty habit of keeping the list around in memory +foreign import ccall unsafe "string.h memset" + c_memset :: Ptr Word8 -> CInt -> CSize -> IO (Ptr Word8) +memset :: Ptr Word8 -> Word8 -> CSize -> IO (Ptr Word8) +memset p w s = c_memset p (fromIntegral w) s diff --git a/MuAtom.hs b/MuAtom.hs --- a/MuAtom.hs +++ b/MuAtom.hs @@ -1,9 +1,9 @@ -{-# LANGUAGE BangPatterns, RecordWildCards #-} +{-# LANGUAGE BangPatterns, RecordWildCards, Rank2Types, FlexibleContexts #-} {- gmndl -- Mandelbrot Set explorer - Copyright (C) 2010 Claude Heiland-Allen + Copyright (C) 2010,2011,2014 Claude Heiland-Allen This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by @@ -21,16 +21,18 @@ -} -module MuAtom (muAtom) where +module MuAtom (muAtom, refineNucleus) where -import Data.Complex (Complex((:+)), cis, magnitude, mkPolar, phase) -import Data.Ratio (denominator) -import Numeric.GSL.Root (RootMethod(Hybrids)) -import qualified Numeric.GSL.Root as GSL +import Data.Ratio (numerator, denominator) +import Data.List (genericIndex, genericSplitAt) +import Data.Vec (NearZero()) +import Numeric.QD (QuadDouble()) +import Roots (root2, root4, lift, FF) +import Complex (Complex((:+)), magnitude, phase, cis, mkPolar, normalize, Turbo) type N = Integer type Q = Rational -type R = Double +type R = QuadDouble type C = Complex R -- a Mandelbrot Set mu-atom @@ -43,70 +45,80 @@ -- finding bond points -f :: N -> C -> C -> C -f !n !z !c - | n == 0 = z - | otherwise = let w = f (n - 1) z c in w * w + c - -df :: N -> C -> C -> C -df !n !z !c = 2 ^ n * product [ f i z c | i <- [0 .. n - 1] ] +fdf :: (Integral i, Num c) => i -> c -> c -> (c, c) +fdf !n !z !c = let (fzs, fz:_) = genericSplitAt n $ iterate (\w -> w * w + c) z + in (fz, 2 ^ n * product fzs) -- [ f i z c | i <- [0 .. n - 1] ] -bondIter :: N -> C -> [R] -> [R] -bondIter !n !b [x0, x1, x2, x3] = - let z = x0:+x1 - c = x2:+x3 - y0 :+ y1 = f n z c - z - y2 :+ y3 = df n z c - b +bondIter :: (Real r, Floating r) => Integer -> Complex r -> FF [] [] r +bondIter !n !(br:+bi) [x0, x1, x2, x3] = + let !z = x0:+x1 + !c = x2:+x3 + !b = lift br :+ lift bi + (!fz, !dfz) = fdf n z c + !(y0 :+ y1) = fz - z -- f n z c - z + !(y2 :+ y3) = dfz - b -- df n z c - (lift br :+ lift bi) in [y0, y1, y2, y3] bondIter _ _ _ = error "MuAtom.bondIter: internal error" -- finding nucleus -l :: N -> C -> C -l !n !c - | n == 0 = 0 - | otherwise = let z = l (n - 1) c in z * z + c +l :: (Integral i, Num c) => i -> c -> c +l !n !c = (`genericIndex` n) . iterate (\z -> z * z + c) $ 0 -nucleusIter :: N -> [R] -> [R] +nucleusIter :: (Real r, Floating r) => Integer -> FF [] [] r nucleusIter !n [x0, x1] = - let c = x0:+x1 - y0 :+ y1 = l n c + let !c = x0 :+ x1 + !(y0 :+ y1) = l n c in [y0, y1] nucleusIter _ _ = error "MuAtom.nucleusIter: internal error" + +refineNucleus :: (NearZero r, Real r, Floating r, Turbo r) => Integer -> Complex r -> (r, r, r) +refineNucleus p root@(gr :+ gi) = + let eps = 1e-20 -- FIXME + [cr, ci] = root2 eps (nucleusIter p) [gr, gi] + [_, _, b0r, b0i] = root4 eps (bondIter p ( 1)) [cr, ci, cr, ci] + [_, _, b1r, b1i] = root4 eps (bondIter p (-1)) [cr, ci, cr, ci] + bond0 = b0r :+ b0i + bond1 = b1r :+ b1i + r = magnitude (bond1 - bond0) + in (cr, ci, r) + -- finding descendants muChild :: Atom -> Q -> Atom -muChild Atom{..} address = +muChild !Atom{..} !address = let -- some properties of the parent and its relation to the child - size = magnitude (root - nucleus) - angle = 2 * pi * realToFrac address - bondAngle = cis angle - child = period * denominator address - solve = GSL.root Hybrids 1e-10 10000 + !size = magnitude (root - nucleus) + !address' = fromIntegral (numerator address) / fromIntegral (denominator address) + !angle = 2 * pi * address' + !(bar :+ bai) = cis angle + !bondAngle = bar :+ bai + !child = period * denominator address -- perturb from the stable nucleus to help ensure convergence to the bond point - _initial@(ir :+ ii) = nucleus + mkPolar (size / 2) (phase (root - nucleus) + angle) - ([_, _, br, bi], _) = solve (bondIter period bondAngle) [ ir, ii, ir, ii ] - bondPoint = br :+ bi + !_initial@(ir :+ ii) = nucleus + mkPolar (size / 2) (phase (root - nucleus) + angle) + [_, _, br, bi] = {-# SCC "bond" #-} root4 eps (bondIter period bondAngle) [ ir, ii, ir, ii ] + !bondPoint = br :+ bi -- estimate where the nucleus will be - radiusEstimate - | cardioid = size / m2 * sin (pi * realToFrac address) + !radiusEstimate + | cardioid = size / m2 * sin (pi * address') | otherwise = size / m2 where m2 = fromIntegral (denominator address) ^ (2 :: N) - deltaEstimate = bondPoint - nucleus - _guess@(gr :+ gi) = bondPoint + mkPolar radiusEstimate (phase deltaEstimate) + !deltaEstimate = bondPoint - nucleus + !_guess@(gr :+ gi) = bondPoint + (radiusEstimate :+ 0) * normalize deltaEstimate -- refine the nucleus estimate - ([cr, ci], _) = solve (nucleusIter child) [gr, gi] - childNucleus = cr :+ ci + [cr, ci] = {-# SCC "nucleus" #-} root2 eps (nucleusIter child) [gr, gi] + !childNucleus = cr :+ ci + eps = radiusEstimate / 10000000 in Atom childNucleus child bondPoint False muChildren :: Atom -> [Q] -> [Atom] -muChildren a [] = [a] -muChildren a (q:qs) = let b = muChild a q in a : muChildren b qs +muChildren !a [] = [a] +muChildren !a (q:qs) = let b = muChild a q in a : muChildren b qs -- interface to the outside world -muAtom :: [Rational] -> (Double, Double, Double, Integer) +muAtom :: [Q] -> (R, R, R, N) muAtom qs = let Atom{..} = last $ muChildren continent qs r :+ i = nucleus diff --git a/Roots.hs b/Roots.hs new file mode 100644 --- /dev/null +++ b/Roots.hs @@ -0,0 +1,78 @@ +{-# LANGUAGE BangPatterns, Rank2Types, ScopedTypeVariables, FlexibleContexts, NoMonomorphismRestriction #-} + +{- + + gmndl -- Mandelbrot Set explorer + Copyright (C) 2010,2011,2014 Claude Heiland-Allen + + This program is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License along + with this program; if not, write to the Free Software Foundation, Inc., + 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + +-} + +module Roots (root2, root4, FF, lift) where + +import Prelude hiding (zipWith) +import Data.Maybe (fromJust) +import Data.Functor ((<$>)) +import Data.Vec (toList, solve, fromList, matFromLists, zipWith, Vec2, Mat22, Vec4, Mat44, NearZero(nearZero)) +import Data.Reflection (Reifies) +import Numeric.AD (jacobian', auto) +import Numeric.AD.Mode.Reverse (Reverse) +import Numeric.AD.Internal.Reverse (Tape) +import Numeric.QD (QuadDouble()) + +lift = auto +type FF f g a = forall s. Reifies s Tape => f (Reverse s a) -> g (Reverse s a) + +root2' :: forall r . (Fractional r, NearZero r, Ord r) => r -> FF [] [] r -> Vec2 r -> Vec2 r +root2' eps f !x = go x + where + jf = jacobian' f + go x0 = + let (ys, js) = unzip $ jf (toList x0) + y = fromList (negate <$> ys) :: Vec2 r + j = matFromLists js :: Mat22 r + dx = fromJust $ solve j y + x1 = zipWith (+) x0 dx + in if all (not . (> eps)) (abs <$> ys) + then x0 + else if all (not . (> eps)) (abs <$> toList dx) + then x1 + else go x1 + +root2 :: (Fractional r, NearZero r, Ord r) => r -> FF [] [] r -> [r] -> [r] +root2 eps f = toList . root2' eps f . fromList + +root4' :: forall r . (Fractional r, NearZero r, Ord r) => r -> FF [] [] r -> Vec4 r -> Vec4 r +root4' eps f !x = go x + where + jf = jacobian' f + go x0 = + let (ys, js) = unzip $ jf (toList x0) + y = fromList (negate <$> ys) :: Vec4 r + j = matFromLists js :: Mat44 r + dx = fromJust $ solve j y + x1 = zipWith (+) x0 dx + in if all (not . (> eps)) (abs <$> ys) + then x0 + else if all (not . (> eps)) (abs <$> toList dx) + then x1 + else go x1 + +root4 :: (Fractional r, NearZero r, Ord r) => r -> FF [] [] r -> [r] -> [r] +root4 eps f = toList . root4' eps f . fromList + +instance NearZero QuadDouble where + nearZero x = not (abs x > 1e-60) -- NearZero Double has 1e-14 diff --git a/gmndl.cabal b/gmndl.cabal --- a/gmndl.cabal +++ b/gmndl.cabal @@ -1,5 +1,5 @@ Name: gmndl -Version: 0.3 +Version: 0.4.0.2 Synopsis: Mandelbrot Set explorer using GTK Description: @@ -11,30 +11,47 @@ . Left-click to zoom in, right-click to zoom out. The status bar shows where you are, and the entry field - takes a space-separated list of fractions strictly - between 0 and 1, try for example: + takes an angled internal address, try for example: . - @1\/2 2\/3 1\/4 3\/5@ + @1\/3 1\/2 5 7@ -Cabal-version: >=1.4 +Cabal-version: >=1.6 License: GPL-2 License-file: LICENSE Author: Claude Heiland-Allen -Maintainer: claudiusmaximus@goto10.org +Maintainer: claude@mathr.co.uk Category: Graphics Build-type: Simple Executable gmndl Main-is: gmndl.hs - Other-modules: MuAtom + Other-modules: Address + Calculate + Complex + Image + MuAtom + Roots Build-depends: base >= 4 && < 5, - array >= 0.3 && < 0.4, - gtk >= 0.11 && < 0.13, - gtkglext >= 0.11 && < 0.13, - hmatrix >= 0.10 && < 0.11, - mtl, - OpenGL >= 2.4 && < 2.5, - priority-queue >= 0.2.1 && < 0.3, - qd >= 0.2 && < 0.3 - GHC-options: -O2 -Wall -threaded -fno-excess-precision -funbox-strict-fields - GHC-prof-options: -O2 -Wall -threaded -fno-excess-precision -funbox-strict-fields -prof -auto-all -caf-all + array >= 0.5 && < 0.6, + gtk >= 0.12 && < 0.13, + gtkglext >= 0.12 && < 0.13, + mtl >= 2.2 && < 2.3, + OpenGL >= 2.9 && < 2.10, + OpenGLRaw >= 1.5 && < 1.6, + parsec >= 3.1 && < 3.2, + priority-queue >= 0.2 && < 0.3, + qd >= 1.0 && < 1.1, + ad >= 4.2 && < 4.3, + reflection >= 1.5 && < 1.6, + Vec >= 1.0 && < 1.1 + GHC-options: -O2 -Wall -threaded -fno-excess-precision -funbox-strict-fields -rtsopts + GHC-prof-options: -O2 -Wall -threaded -fno-excess-precision -funbox-strict-fields -rtsopts -prof -auto-all -caf-all + +source-repository head + type: git + location: git://gitorious.org/maximus/gmndl.git + +source-repository this + type: git + location: git://gitorious.org/maximus/gmndl.git + tag: v0.4.0.2 diff --git a/gmndl.hs b/gmndl.hs --- a/gmndl.hs +++ b/gmndl.hs @@ -1,9 +1,7 @@ -{-# LANGUAGE BangPatterns, ForeignFunctionInterface #-} - {- gmndl -- Mandelbrot Set explorer - Copyright (C) 2010 Claude Heiland-Allen + Copyright (C) 2010,2011,2014 Claude Heiland-Allen This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by @@ -23,445 +21,209 @@ module Main (main) where +import Prelude hiding (isNaN) + +import Control.Concurrent (forkIO, newEmptyMVar, takeMVar, putMVar) +import Control.Monad (forever) + -- some simple helpers -import Control.Monad (when) import Data.List (isPrefixOf) --- concurrent renderer with capability-specific scheduling -import Control.Concurrent (forkIO, killThread) -import GHC.Conc (forkOnIO, numCapabilities) - -- the dependency on mtl is just for this! import Control.Monad.Trans (liftIO) --- each worker uses a mutable unboxed array of Bool to know which pixels --- it has already started to render, to avoid pointless work duplication -import Data.Array.IO (IOUArray, newArray, readArray, writeArray, inRange) - -- the main program thread needs to store some thread-local state import Data.IORef (newIORef, readIORef, writeIORef) --- each worker thread keeps a queue of pixels that it needs to render or --- to continue rendering later -import Data.PriorityQueue (PriorityQueue, newPriorityQueue, enqueue, enqueueBatch, dequeue) - --- poking bytes into memory is dirty, but it's quick and allows use of --- other fast functions like memset and easy integration with OpenGL -import Foreign (castPtr, mallocBytes, nullPtr, plusPtr, pokeByteOff, Ptr, Word8) -import Foreign.C (CDouble, CInt, CSize) - -- build the interface with GTK to allow more fancy controls later import Graphics.UI.Gtk -- use OpenGL to display frequently update images on a textured quad import Graphics.UI.Gtk.OpenGL import qualified Graphics.Rendering.OpenGL as GL -import Graphics.Rendering.OpenGL (($=), GLfloat) +import Graphics.Rendering.OpenGL (($=)) --- higher precision arithmetic using libqd -import Numeric.QD.DoubleDouble (DoubleDouble(DoubleDouble)) -import Numeric.QD.QuadDouble (QuadDouble(QuadDouble)) -import qualified Numeric.QD.DoubleDouble as DD -import qualified Numeric.QD.QuadDouble as QD -import Numeric.QD.FPU.Raw (fpu_fix_start) +-- need a hack to ensure correct qd operation +import Numeric.QD.QuadDouble (QuadDouble()) +import Foreign (nullPtr) --- ugly! but the default realToFrac :: (C)Double -> (C)Double is slooow -import Unsafe.Coerce (unsafeCoerce) +import Complex (Complex((:+))) -- mu-atom properties -import MuAtom (muAtom) - --- some type aliases to shorten things -type B = Word8 -type N = Int -type R = Double - --- don't look! this is really really ugly, and should be benchmarked --- to see how really necessary it is, or at least made into a type class -convert :: (Real a, Fractional b) => a -> b -convert = realToFrac -convertDouble2CDouble :: Double -> CDouble -convertDouble2CDouble !x = unsafeCoerce x -convertCDouble2Double :: CDouble -> Double -convertCDouble2Double !x = unsafeCoerce x -convertDouble2DoubleDouble :: Double -> DoubleDouble -convertDouble2DoubleDouble !x = convertCDouble2DoubleDouble . convertDouble2CDouble $ x -convertCDouble2DoubleDouble :: CDouble -> DoubleDouble -convertCDouble2DoubleDouble !x = DoubleDouble x 0 -convertDoubleDouble2Double :: DoubleDouble -> Double -convertDoubleDouble2Double !(DoubleDouble x _) = convertCDouble2Double x -convertDoubleDouble2CDouble :: DoubleDouble -> CDouble -convertDoubleDouble2CDouble !(DoubleDouble x _) = x -{-# RULES "convert/Double2CDouble" convert = convertDouble2CDouble #-} -{-# RULES "convert/CDouble2Double" convert = convertCDouble2Double #-} -{-# RULES "convert/Double2DoubleDouble" convert = convertDouble2DoubleDouble #-} -{-# RULES "convert/CDouble2DoubleDouble" convert = convertCDouble2DoubleDouble #-} -{-# RULES "convert/DoubleDouble2Double" convert = convertDoubleDouble2Double #-} -{-# RULES "convert/DoubleDouble2CDouble" convert = convertDoubleDouble2CDouble #-} - --- this is ugly too: can't use Data.Complex because the qd bindings do --- not implement some low-level functions properly, leading to obscure --- crashes inside various Data.Complex functions... -data Complex c = {-# UNPACK #-} !c :+ {-# UNPACK #-} !c deriving (Read, Show, Eq) - --- complex number arithmetic, with extra strictness and cost-centres -instance Num c => Num (Complex c) where - (!(a :+ b)) + (!(c :+ d)) = {-# SCC "C+" #-} ((a + c) :+ (b + d)) - (!(a :+ b)) - (!(c :+ d)) = {-# SCC "C-" #-} ((a - c) :+ (b - d)) - (!(a :+ b)) * (!(c :+ d)) = {-# SCC "C*" #-} ((a * c - b * d) :+ (a * d + b * c)) - negate !(a :+ b) = (-a) :+ (-b) - abs x = error $ "Complex.abs: " ++ show x - signum x = error $ "Complex.signum: " ++ show x - fromInteger !x = fromInteger x :+ 0 - --- an extra class for some operations that can be made faster for things --- like DoubleDouble: probably should have given this a better name -class Num c => Turbo c where - sqr :: c -> c - sqr !x = x * x - twice :: c -> c - twice !x = x + x - --- the default methods are fine for simple primitive types... -instance Turbo Float where -instance Turbo Double where -instance Turbo CDouble where - --- ...and complex numbers -instance Turbo c => Turbo (Complex c) where - sqr !(r :+ i) = (sqr r - sqr i) :+ (twice (r * i)) - twice !(r :+ i) = (twice r) :+ (twice i) - --- use the specific implementations for the higher precision types -instance Turbo DoubleDouble where - sqr !x = DD.sqr x - twice !(DoubleDouble a b) = DoubleDouble (twice a) (twice b) - -instance Turbo QuadDouble where - sqr !x = QD.sqr x - twice !(QuadDouble a b c d) = QuadDouble (twice a) (twice b) (twice c) (twice d) - --- colour space conversion from HSV [0..1] to RGB [0..1] --- HSV looks quite 'chemical' to my eyes, need to investigate something --- better to make it feel more 'natural' -hsv2rgb :: R -> R -> R -> (R, R, R) -hsv2rgb !h !s !v - | s == 0 = (v, v, v) - | h == 1 = hsv2rgb 0 s v - | otherwise = - let !i = floor (h * 6) `mod` 6 :: N - !f = (h * 6) - fromIntegral i - !p = v * (1 - s) - !q = v * (1 - s * f) - !t = v * (1 - s * (1 - f)) - in case i of - 0 -> (v, t, p) - 1 -> (q, v, p) - 2 -> (p, v, t) - 3 -> (p, q, v) - 4 -> (t, p, v) - 5 -> (v, p, q) - _ -> (0, 0, 0) - --- compute RGB [0..255] bytes from the results of the complex iterations --- don't need very high precision for this, as spatial aliasing will be --- much more of a problem in intricate regions of the fractal -colour :: Complex Double -> Complex Double -> N -> (B, B, B) -colour !(zr:+zi) !(dzr:+dzi) !n = - let -- micro-optimization - there is no log2 function - !il2 = 1 / log 2 - !zd2 = sqr zr + sqr zi - !dzd2 = sqr dzr + sqr dzi - -- normalized escape time - !d = (fromIntegral n :: R) - log (log zd2 / log escapeR2) * il2 - !dwell = fromIntegral (floor d :: N) - -- final angle of the iterate - !finala = atan2 zi zr - -- distance estimate - !de = (log zd2 * il2) * sqrt (zd2 / dzd2) - !dscale = log de * il2 + 32 - -- HSV is based on escape time, distance estimate, and angle - !hue = log d * il2 / 3 - !saturation = 0 `max` (log d * il2 / 8) `min` 1 - !value = 0 `max` (1 - dscale / 48) `min` 1 - !h = hue - fromIntegral (floor hue :: N) - -- adjust saturation to give concentric striped pattern - !k = dwell / 2 - !satf = if k - fromIntegral (floor k :: N) >= (0.5 :: R) then 0.9 else 1 - -- adjust value to give tiled pattern - !valf = if finala < 0 then 0.9 else 1 - -- convert to RGB - (!r, !g, !b) = hsv2rgb h (satf * saturation) (valf * value) - -- convert to bytes - !rr = floor $ 0 `max` 255 * r `min` 255 - !gg = floor $ 0 `max` 255 * g `min` 255 - !bb = floor $ 0 `max` 255 * b `min` 255 - in (rr, gg, bb) - --- a Job stores a pixel undergoing iterations -data Job c = Job !N !N !(Complex c) !(Complex c) !(Complex c) !N - --- the priority of a Job is how many iterations have been computed: --- so 'fresher' pixels drop to the front of the queue in the hope of --- avoiding too much work iterating pixels that will never escape -priority :: Job c -> N -priority !(Job _ _ _ _ _ n) = n - --- add a job to a work queue, taking care not to duplicate work --- there is no race condition here as each worker has its own queue -addJob :: RealFloat c => N -> N -> Complex c -> N -> PriorityQueue IO (Job c) -> IOUArray (N,N) Bool -> N -> N -> IO () -addJob !w !h !c !zoom todo sync !i !j = do - already <- readArray sync (j, i) - when (not already) $ do - writeArray sync (j, i) True - enqueue todo $! Job i j (coords w h c zoom i j) 0 0 0 - --- spawns a new batch of workers to render an image --- returns an action that stops the rendering -renderer :: (Turbo c, RealFloat c) => ((N,N),(N,N)) -> (N -> N -> B -> B -> B -> IO ()) -> Complex c -> N -> IO (IO ()) -renderer rng output !c !zoom = do - workerts <- mapM (\w -> forkOnIO w $ worker rng c zoom output w) [ 0 .. workers - 1 ] - return $ do - mapM_ killThread workerts - --- compute the Complex 'c' coordinate for a pixel in the image -coords :: RealFloat c => N -> N -> Complex c -> N -> N -> N -> Complex c -coords !w !h !c !zoom !i !j = c + ( (fromIntegral (i - w`div`2) * k) - :+(fromIntegral (h`div`2 - j) * k)) - where !k = convert (1/2^^zoom :: Double) - --- start rendering pixels from the edge of the image --- the Mandelbrot Set and its complement are both simply-connected --- discounting spatial aliasing any point inside the boundary that is --- in the complement is 'reachable' from a point on the boundary that --- is also in the complement - probably some heavy math involved to --- prove this though --- note: this implicitly depends on the spread values below - it's --- necessary for each interlaced subimage (one per worker) to have --- at least a one pixel deep border -border :: N -> N -> [(N, N)] -border !w !h = concat $ - [ [ (j, i) | i <- [ 0 .. w - 1 ], j <- [ 0 .. workers - 1 ] ] - , [ (j, i) | j <- [ 0 .. h - 1 ], i <- [ 0 .. workers - 1 ] ] - , [ (j, i) | j <- [ 0 .. h - 1 ], i <- [ w - workers .. w - 1 ] ] - , [ (j, i) | i <- [ 0 .. w - 1 ], j <- [ h - workers .. h - 1 ] ] - ] - --- the worker thread enqueues its border and starts computing iterations -worker :: (Turbo c, RealFloat c) => ((N,N),(N,N)) -> Complex c -> N -> (N -> N -> B -> B -> B -> IO ()) -> N -> IO () -worker rng@((y0,x0),(y1,x1)) !c !zoom output !me = do - sync <- newArray rng False - queue <- newPriorityQueue priority - let addJ = addJob w h c zoom queue sync - js = filter mine (border w h) - w = x1 - x0 + 1 - h = y1 - y0 + 1 - mapM_ (flip (writeArray sync) True) js - enqueueBatch queue (map (\(j,i) -> Job i j (coords w h c zoom i j) 0 0 0) js) - compute rng addJ output queue - where mine (j, _) = j `mod` workers == me -- another dependency on spread - --- the compute engine pulls pixels from the queue until there are no --- more, and calculates a batch of iterations for each -compute :: (Turbo c, RealFloat c) => ((N,N),(N,N)) -> (N -> N -> IO ()) -> (N -> N -> B -> B -> B -> IO ()) -> PriorityQueue IO (Job c) -> IO () -compute rng addJ output queue = do - mjob <- dequeue queue - case mjob of - Just (Job i j c z dz n) -> do - let -- called when the pixel escapes - done' !(zr:+zi) !(dzr:+dzi) !n' = {-# SCC "done" #-} do - let (r, g, b) = colour (convert zr :+ convert zi) (convert dzr :+ convert dzi) n' - output i j r g b - -- a wavefront of computation spreads to neighbouring pixels - sequence_ - [ addJ x y - | u <- spreadX - , v <- spreadY - , let x = i + u - , let y = j + v - , inRange rng (y, x) - ] - -- called when the pixel doesn't escape yet - todo' !z' !dz' !n' = {-# SCC "todo" #-} enqueue queue $! Job i j c z' dz' n' - calculate c limit z dz n done' todo' - compute rng addJ output queue - Nothing -> return () -- no pixels left to render, so finish quietly +import Calculate +import qualified Image +import Address (parse, parameter) --- the raw z->z^2+c calculation engine --- also computes the derivative for distance estimation calculations --- this function is crucial for speed, too much allocation will slooow --- everything down severely -calculate :: (Turbo c, RealFloat c) => Complex c -> N -> Complex c -> Complex c -> N -> (Complex c -> Complex c -> N -> IO ()) -> (Complex c -> Complex c -> N -> IO ()) -> IO () -calculate !c !m0 !z0 !dz0 !n0 done todo = go m0 z0 dz0 n0 - where - go !m !z@(zr:+zi) !dz !n - | not (sqr zr + sqr zi < er2) = done z dz n - | m <= 0 = todo z dz n - | otherwise = go (m - 1) (sqr z + c) (let !zdz = z * dz in twice zdz + 1) (n + 1) - !er2 = convert escapeR2 +isNaN x = not (x == x) --- dispatch to different instances of renderer depending on required precision --- if zoom is low, single precision Float is ok, but as soon as pixel spacing --- gets really small, it's necessary to increase it --- it's probably not even worth using Float - worth benchmarking this and --- also the DD and QD types (which cause a massively noticeable slowdown) -renderer' :: Real c => ((N,N),(N,N)) -> (N -> N -> B -> B -> B -> IO ()) -> Complex c -> N -> IO (IO ()) -renderer' rng output !c !zoom - | zoom < 20 = {-# SCC "rF" #-} renderer rng output (f c :: Complex Float ) zoom - | zoom < 50 = {-# SCC "rD" #-} renderer rng output (f c :: Complex Double ) zoom - | zoom < 100 = {-# SCC "rDD" #-} renderer rng output (f c :: Complex DoubleDouble) zoom - | otherwise = {-# SCC "rQD" #-} renderer rng output (f c :: Complex QuadDouble ) zoom - where f !(cr :+ ci) = convert cr :+ convert ci +-- the state we need for everything +data GMNDL + = Invalid + | GMNDL + { center :: Complex QuadDouble + , zradius :: QuadDouble + , image :: Image.Image + , stop :: IO () + } -- command line arguments: currently only initial window dimensions -data Args = Args{ aWidth :: N, aHeight :: N } +data Args = Args{ aWidth :: Int, aHeight :: Int, aOversample :: Int, aRe :: QuadDouble, aIm :: QuadDouble, aZr :: QuadDouble } -- and the defaults are suitable for PAL DVD rendering, if that should -- come to pass in the future defaultArgs :: Args -defaultArgs = Args{ aWidth = 788, aHeight = 576 } +defaultArgs = Args{ aWidth = 788, aHeight = 576, aOversample = 1, aRe = 0, aIm = 0, aZr = 2 } -- braindead argument parser: latest argument takes priority -- probably should use Monoid instances for this stuff combineArgs :: Args -> String -> Args combineArgs a0 s | "--width=" `isPrefixOf` s = a0{ aWidth = read $ "--width=" `dropPrefix` s } + | "-w=" `isPrefixOf` s = a0{ aWidth = read $ "-w=" `dropPrefix` s } | "--height=" `isPrefixOf` s = a0{ aHeight = read $ "--height=" `dropPrefix` s } - | "-w=" `isPrefixOf` s = a0{ aWidth = read $ "-w=" `dropPrefix` s } - | "-h=" `isPrefixOf` s = a0{ aHeight = read $ "-h=" `dropPrefix` s } + | "-h=" `isPrefixOf` s = a0{ aHeight = read $ "-h=" `dropPrefix` s } + | "--aa=" `isPrefixOf` s = a0{ aOversample = read $ "--aa=" `dropPrefix` s } + | "--re=" `isPrefixOf` s = a0{ aRe = read $ "--re=" `dropPrefix` s } + | "--im=" `isPrefixOf` s = a0{ aIm = read $ "--im=" `dropPrefix` s } + | "--zr=" `isPrefixOf` s = a0{ aZr = read $ "--zr=" `dropPrefix` s } | otherwise = a0 -- this is a bit silly, especially with the duplicated string literals.. dropPrefix :: String -> String -> String dropPrefix p s = drop (length p) s --- round up to nearest power of two --- this will probably explode when n gets large, but it's only used --- for OpenGL texture dimensions so you'll run out of memory first -roundUp2 :: N -> N -roundUp2 n = head . dropWhile (< n) . iterate (2*) $ 1 - -- the main program! main :: IO () main = do args <- foldl combineArgs defaultArgs `fmap` unsafeInitGUIForThreadedRTS let width = aWidth args height = aHeight args - size = roundUp2 (width `max` height) - rng = ((0, 0), (height - 1, width - 1)) + oversample = aOversample args + rng = ((0, 0), (oversample * height - 1, oversample * width - 1)) _ <- initGL glconfig <- glConfigNew [ GLModeRGBA, GLModeDouble ] canvas <- glDrawingAreaNew glconfig widgetSetSizeRequest canvas width height - -- allocate some image bytes and clear with white - imgdata <- mallocBytes $ width * height * 3 - let clear = memset (castPtr imgdata) 255 (fromIntegral $ height * width * 3) >> return () - clear - let output x y r g b = do - let p = imgdata `plusPtr` ((y * width + x) * 3) - pokeByteOff p 0 r - pokeByteOff p 1 g - pokeByteOff p 2 b window <- windowNew eventb <- eventBoxNew vbox <- vBoxNew False 0 status <- vBoxNew False 0 - statusRe <- labelNew Nothing - statusIm <- labelNew Nothing - statusZo <- labelNew Nothing + statusRe <- entryNew + statusIm <- entryNew + statusZr <- entryNew ratios <- entryNew boxPackStart vbox eventb PackGrow 0 boxPackStart vbox status PackGrow 0 boxPackStart vbox ratios PackGrow 0 boxPackStart status statusRe PackGrow 0 boxPackStart status statusIm PackGrow 0 - boxPackStart status statusZo PackGrow 0 - let updateStatus re im zo = do -- this updates the status bar - labelSetText statusRe (reshow $ show re) - labelSetText statusIm (reshow $ show im) - labelSetText statusZo (show zo) + boxPackStart status statusZr PackGrow 0 + let -- update the status bar + updateStatus re im zr = do + entrySetText statusRe (show re) + entrySetText statusIm (show im) + entrySetText statusZr (show zr) set window [ containerBorderWidth := 0, containerChild := vbox, windowResizable := False ] set eventb [ containerBorderWidth := 0, containerChild := canvas ] - -- mouse motion events not sent by default as there can a flood - widgetAddEvents eventb [PointerMotionMask] - -- dirty hack to set FPU control words as recommended by libqd docs - -- because it relies on 64bit doubles and some FPU use 80bits inside - mapM_ (flip forkOnIO $ fpu_fix_start nullPtr) [ 0 .. numCapabilities - 1 ] - -- start the renderer for the first time - stop0 <- renderer' rng output c0 zoom0 - -- save the initial state - sR <- newIORef (c0, zoom0, stop0) - -- when the mouse moves, update the status bar with the current coords - _ <- eventb `on` motionNotifyEvent $ {-# SCC "cbMo" #-} tryEvent $ do - (x, y) <- eventCoordinates - liftIO $ do - (c, zoom, _stop) <- readIORef sR - let _c'@(re' :+ im') = c + ((convert x :+ convert (-y)) - (fromIntegral width :+ fromIntegral (-height)) * (0.5 :+ 0)) * ((1/2^^zoom) :+ 0) - updateStatus re' im' zoom - -- when the mouse button is pressed, center and zoom in + -- initial state is invalid because... + sR <- newIORef Invalid + done <- newEmptyMVar + let -- restart the renderer + restart :: IO () + restart = do + g <- readIORef sR + stop g + Image.clear (image g) + stop' <- renderer done rng (Image.plot (image g)) (center g) (zradius g) + writeIORef sR $! g{ stop = stop' } + let re :+ im = center g + updateStatus re im (zradius g) + -- ...need to initialize OpenGL stuff etc in this callback + _ <- onRealize canvas $ {-# SCC "cbRz" #-} withGLDrawingArea canvas $ \_ -> do + GL.matrixMode $= GL.Projection + GL.loadIdentity + GL.ortho 0.0 1.0 0.0 1.0 (-1.0) 1.0 + GL.drawBuffer $= GL.BackBuffers + GL.texture GL.Texture2D $= GL.Enabled + i <- Image.new (oversample * width) (oversample * height) + writeIORef sR $! GMNDL{ image = i, center = aRe args :+ aIm args, zradius = aZr args, stop = return () } + restart + -- when the mouse button is pressed, center and zoom in or out _ <- eventb `on` buttonPressEvent $ {-# SCC "cbEv" #-} tryEvent $ do b <- eventButton (x, y) <- eventCoordinates liftIO $ do - (c, zoom, stop) <- readIORef sR - stop - clear - let c'@(re' :+ im') = c + ((convert x :+ convert (-y)) - (fromIntegral width :+ fromIntegral (-height)) * (0.5 :+ 0)) * ((1/2^^zoom) :+ 0) - zoom' = zoom + delta - delta | b == LeftButton = 1 - | b == RightButton = -1 - | otherwise = 0 - stop' <- renderer' rng output c' zoom' - writeIORef sR $! (c', zoom', stop') - updateStatus re' im' zoom' + g <- readIORef sR + let w2 = fromIntegral width / 2 + h2 = fromIntegral height / 2 + p = convert x :+ convert (-y) + s = (zradius g / (w2 `min` h2)) :+ 0 + c = center g + (p - (w2 :+ (-h2))) * s + zradius' = zradius g * delta + delta | b == LeftButton = 0.5 + | b == RightButton = 2 + | otherwise = 1 + writeIORef sR $! g{ center = c, zradius = zradius' } + restart + -- when typing in the coordinate boxes, zoom to the new place + _ <- statusRe `onEntryActivate` do + s <- entryGetText statusRe + liftIO $ do + g <- readIORef sR + case safeRead s of + Just re -> do + let _ :+ im = center g + writeIORef sR $! g{ center = re :+ im } + restart + Nothing -> return () + _ <- statusIm `onEntryActivate` do + s <- entryGetText statusIm + liftIO $ do + g <- readIORef sR + case safeRead s of + Just im -> do + let re :+ _ = center g + writeIORef sR $! g{ center = re :+ im } + restart + Nothing -> return () + _ <- statusZr `onEntryActivate` do + s <- entryGetText statusZr + liftIO $ do + g <- readIORef sR + case safeRead s of + Just zradius' -> do + writeIORef sR $! g{ zradius = zradius' } + restart + Nothing -> return () -- when pressing return in the ratios list, zoom to that mu-atom + muQueue <- newEmptyMVar + _ <- forkIO . forever $ do + qs <- takeMVar muQueue + case parameter =<< parse qs of + Nothing -> postGUISync $ do + _ <- ratios `widgetSetSensitive` True + return () + Just (cr, ci, radius) -> do + let zradius' = radius * 3 + cr `seq` ci `seq` zradius' `seq` postGUISync $ do + g <- readIORef sR + if isNaN cr || isNaN ci + then writeIORef sR $! g{ center = c0, zradius = zradius0 } + else writeIORef sR $! g{ center = cr :+ ci, zradius = zradius' } + _ <- ratios `widgetSetSensitive` True + restart _ <- ratios `onEntryActivate` do s <- entryGetText ratios - case rationalize s of - Nothing -> return () - Just qs -> do - _ <- ratios `widgetSetSensitive` False - _ <- forkIO $ do - let (cr, ci, radius, _period) = muAtom qs - zoom' = floor $ (logBase 2 . fromIntegral $ width `min` height) - (logBase 2 radius) - 2 - c'@(cr' :+ ci') = convert cr :+ convert ci - cr `seq` ci `seq` zoom' `seq` postGUISync $ do - (_, _, stop) <- readIORef sR - stop - clear - stop' <- renderer' rng output c' zoom' - writeIORef sR $! (c', zoom', stop') - _ <- ratios `widgetSetSensitive` True - updateStatus cr' ci' zoom' - return () - -- need to set up OpenGL stuff in callback just because that's how... - _ <- onRealize canvas $ {-# SCC "cbRz" #-} withGLDrawingArea canvas $ \_ -> do - GL.clearColor $= (GL.Color4 0.0 0.0 0.0 0.0) - GL.matrixMode $= GL.Projection - GL.loadIdentity - GL.ortho 0.0 1.0 0.0 1.0 (-1.0) 1.0 - GL.drawBuffer $= GL.BackBuffers - -- create a new texture and pre-allocate it to a square 2^n size for speed - [tex] <- GL.genObjectNames 1 - GL.texture GL.Texture2D $= GL.Enabled - GL.textureBinding GL.Texture2D $= Just tex - GL.texImage2D Nothing GL.NoProxy 0 GL.RGB' (GL.TextureSize2D (fromIntegral size) (fromIntegral size)) 0 (GL.PixelData GL.RGB GL.UnsignedByte nullPtr) - GL.textureFilter GL.Texture2D $= ((GL.Nearest, Nothing), GL.Nearest) - GL.textureWrapMode GL.Texture2D GL.S $= (GL.Repeated, GL.ClampToEdge) - GL.textureWrapMode GL.Texture2D GL.T $= (GL.Repeated, GL.ClampToEdge) + _ <- ratios `widgetSetSensitive` False + g <- readIORef sR + stop g + putMVar muQueue s -- time to draw the image: upload to the texture and draw a quad _ <- onExpose canvas $ {-# SCC "cbEx" #-} \_ -> do withGLDrawingArea canvas $ \glwindow -> do - let v :: GLfloat -> GLfloat -> GLfloat -> GLfloat -> IO () - v tx ty vx vy = GL.texCoord (GL.TexCoord2 tx ty) >> GL.vertex (GL.Vertex2 vx vy) - w = fromIntegral width - h = fromIntegral height - sx = fromIntegral width / fromIntegral size - sy = fromIntegral height / fromIntegral size - GL.clear [ GL.ColorBuffer ] - GL.texSubImage2D Nothing 0 (GL.TexturePosition2D 0 0) (GL.TextureSize2D w h) (GL.PixelData GL.RGB GL.UnsignedByte imgdata) - GL.renderPrimitive GL.Quads $ do - v 0 sy 0 0 >> v 0 0 0 1 >> v sx 0 1 1 >> v sx sy 1 0 + GMNDL{ image = i } <- readIORef sR + Image.upload i + Image.draw i glDrawableSwapBuffers glwindow return True -- need an exit strategy @@ -472,28 +234,6 @@ widgetShowAll window mainGUI --- which neighbours to activate once a pixel has escaped --- there are essentially two choices, with x<->y swapped --- choose greater X spread because images are often wider than tall --- other schemes wherein the spread is split in both directions --- might benefit appearance with large worker count, but too complicated -spreadX, spreadY :: [ N ] -spreadX = [ -workers, 0, workers ] -spreadY = [ -1, 0, 1 ] - --- number of worker threads --- use as many worker threads as capabilities, with the workers --- distributed 1-1 onto capabilities to maximize CPU utilization -workers :: N -workers = numCapabilities - --- iteration limit per pixel --- at most this many iterations are performed on each pixel before it --- is shunted to the back of the work queue --- this should be tuneable to balance display updates against overheads -limit :: N -limit = (2^(11::N)-1) - -- initial center coordinates -- using the maximum precision available from the start for this makes -- sure that nothing weird happens when precision gets close to the edge @@ -501,52 +241,11 @@ c0 = 0 -- initial zoom level --- neighbouring pixel are 2^(-zoom) units apart -- the initial zoom level should probably depend on initial image size -zoom0 :: N -zoom0 = 6 - --- escape radius for fractal iteration calculations --- once the complex iterate exceeds this, it's never coming back --- theoretically escapeR = 2 would work --- but higher values like this give a significantly smoother picture -escapeR, escapeR2 :: R -escapeR = 65536 -escapeR2 = escapeR * escapeR - --- import standard C library memset for clearing images efficiently --- previous implementation used pokeArray ... (replicate ...) ... --- which had a nasty habit of keeping the list around in memory -foreign import ccall unsafe "string.h memset" c_memset :: Ptr Word8 -> CInt -> CSize -> IO (Ptr Word8) -memset :: Ptr Word8 -> Word8 -> CSize -> IO (Ptr Word8) -memset p w s = c_memset p (fromIntegral w) s - --- convert scientific notation like "-4.12345600000000e-03" --- into human-readable numbers like "-0.004123456" --- do it by string manipulation as QuadDouble has a lot of precision --- and the show instance for QuadDouble gives the problematic form... -reshow :: String -> String -reshow s = - let (front, 'e':rear) = break (=='e') s - (sign, mantissa) = case front of - '-':m -> (True, m) - m -> (False, m) - (big, '.':small) = break (=='.') mantissa - expo = read (dropWhile (=='+') rear) + length big - digits = if expo < 0 then "0." ++ replicate (-expo) '0' ++ big ++ small - else let (x,y) = splitAt expo (big ++ tail small) - in case x of - [] -> "0." ++ y - x' -> x' ++ "." ++ y - in reverse . dropWhile (=='.') . dropWhile (=='0') . reverse . (if sign then ('-':) else id) $ digits +zradius0 :: QuadDouble +zradius0 = 2 --- convert a string of space separated fractions into rationals --- moreover ensure that they are all in 0 Maybe [Rational] -rationalize s = - let f '/' = '%' - f c = c - r t = case reads t of - [(q, "")] | 0 < q && q < 1 -> Just q - _ -> Nothing - in mapM r . words . map f $ s +safeRead :: Read a => String -> Maybe a +safeRead s = case reads s of + [(a, "")] -> Just a + _ -> Nothing