diff --git a/Fractal/RUFF/Mandelbrot/Address.hs b/Fractal/RUFF/Mandelbrot/Address.hs new file mode 100644 --- /dev/null +++ b/Fractal/RUFF/Mandelbrot/Address.hs @@ -0,0 +1,407 @@ +{-# LANGUAGE DeriveDataTypeable #-} +{- | +Module : Fractal.RUFF.Mandelbrot.Address +Copyright : (c) Claude Heiland-Allen 2010-2011 +License : BSD3 + +Maintainer : claudiusmaximus@goto10.org +Stability : unstable +Portability : portable + +External angles give rise to kneading sequences under the angle doubling +map. Internal addresses encode kneading sequences in human-readable form, +when extended to angled internal addresses they distinguish hyperbolic +components in a concise and meaningful way. + +The algorithms are mostly based on Dierk Schleicher's paper +/Internal Addresses Of The Mandelbrot Set And Galois Groups Of Polynomials (version of February 5, 2008)/ +. +-} + +module Fractal.RUFF.Mandelbrot.Address + ( Angle, double, wrap + , Knead(..), Kneading(..), kneading, period, unwrap + , InternalAddress(..), internalAddress, associated, upper, lower, internalFromList, internalToList + , AngledInternalAddress(..), angledInternalAddress, angledFromList, angledToList, externalAngles + , stripAngles + , parse + ) where + +import Data.Data (Data()) +import Data.Typeable (Typeable()) +import Control.Monad (guard) +import Control.Monad.Identity (Identity()) +import Data.Char (digitToInt) +import Data.List (genericDrop, genericIndex, genericLength, genericReplicate, genericSplitAt, genericTake) +import Data.Maybe (isJust, listToMaybe) +import Data.Ratio ((%), numerator, denominator) +import Text.Parsec (ParsecT(), choice, digit, eof, many, many1, runP, sepEndBy, string, try) +import Text.PrettyPrint.Leijen.Text (Pretty, pretty, prettyList, char, parens, (<>)) + +-- | Angle as a fraction of a turn, usually in [0, 1). +type Angle = Rational + +-- | Wrap an angle into [0, 1). +wrap :: Angle -> Angle +wrap a + | f < 0 = 1 + f + | otherwise = f + where + (_, f) = properFraction a :: (Integer, Angle) + +-- | Angle doubling map. +double :: Angle -> Angle +double a = wrap (2 * a) + +-- | Elements of kneading sequences. +data Knead + = Zero + | One + | Star + deriving (Read, Show, Eq, Ord, Enum, Bounded, Data, Typeable) + +instance Pretty Knead where + pretty = char . kneadChar + prettyList = pretty . map kneadChar + +kneadChar :: Knead -> Char +kneadChar Zero = '0' +kneadChar One = '1' +kneadChar Star = '*' + +-- | Kneading sequences. Note that the 'Aperiodic' case has an infinite list, +-- which the 'Pretty' instance truncates arbitrarily. +data Kneading + = Aperiodic [Knead] + | PrePeriodic [Knead] [Knead] + | StarPeriodic [Knead] + | Periodic [Knead] + deriving (Read, Show, Eq, Ord, Data, Typeable) + +instance Pretty Kneading where + pretty (Aperiodic ks) = pretty . (++ "ยทยทยท") . map kneadChar . take 17 $ ks + pretty (PrePeriodic us vs) = pretty us <> parens (pretty vs) + pretty (StarPeriodic vs) = parens (pretty vs) + pretty (Periodic vs) = parens (pretty vs) + +-- | The kneading sequence for an external angle. +kneading :: Angle -> Kneading +kneading a0' + | a0 == 0 = StarPeriodic [Star] + | otherwise = fst kneads + where + a0 = wrap a0' + lo = a0 / 2 + hi = (a0 + 1) / 2 + kneads = kneading' 1 (double a0) + ks = (a0, One) : snd kneads + kneading' :: Integer -> Angle -> (Kneading, [(Angle, Knead)]) + kneading' n a + | isJust i = case i of + Just 0 -> case last qs of + Star -> (StarPeriodic qs, []) + _ -> (Periodic qs, []) + Just j -> let (p, q) = genericSplitAt j qs + in (PrePeriodic p q, []) + _ -> error "Fractal.Mandelbrot.Address.kneading (isJust -> Nothing?)" + | a == lo = ((a, Star):) `mapP` k + | a == hi = ((a, Star):) `mapP` k + | lo < a && a < hi = ((a, One ):) `mapP` k + | hi < a || a < lo = ((a, Zero):) `mapP` k + | otherwise = error "Fractal.Mandelbrot.Address.kneading (unmatched?)" + where + k = kneading' (n+1) (double a) + ps = genericTake n ks + qs = map snd ps + i = fmap fst . listToMaybe . filter ((a ==) . fst . snd) . zip [(0 :: Integer) ..] $ ps + mapP f ~(x, y) = (x, f y) + +-- | The period of a kneading sequence, or 'Nothing' when it isn't periodic. +period :: Kneading -> Maybe Integer +period (StarPeriodic k) = Just (genericLength k) +period (Periodic k) = Just (genericLength k) +period _ = Nothing + +rho :: Kneading -> Integer -> Integer +rho v r | r >= 1 && fmap (r`mod`) (period v) /= Just 0 = ((1 + r) +) . genericLength . takeWhile id . zipWith (==) vs . genericDrop r $ vs + | otherwise = rho v (r + 1) + where + vs = unwrap v + +-- | Unwrap a kneading sequence to an infinite list. +unwrap :: Kneading -> [Knead] +unwrap (Aperiodic vs) = vs +unwrap (PrePeriodic us vs) = us ++ cycle vs +unwrap (StarPeriodic vs) = cycle vs +unwrap (Periodic vs) = cycle vs + +orbit :: (a -> a) -> a -> [a] +orbit = iterate + +-- | Internal addresses are a non-empty sequence of strictly increasing +-- integers beginning with '1'. +data InternalAddress = InternalAddress [Integer] + deriving (Read, Show, Eq, Ord, Data, Typeable) + +instance Pretty InternalAddress where + pretty (InternalAddress []) = error "Fractal.Mandelbrot.Address.InternalAddress.pretty" + pretty (InternalAddress [x]) = pretty x + pretty (InternalAddress (x:ys)) = pretty x <> char ' ' <> pretty (InternalAddress ys) + +-- | Construct a valid 'InternalAddress', checking the precondition. +internalFromList :: [Integer] -> Maybe InternalAddress +internalFromList x0s@(1:_) = InternalAddress `fmap` fromList' 0 x0s + where + fromList' n [x] | x > n = Just [x] + fromList' n (x:xs) | x > n = (x:) `fmap` fromList' x xs + fromList' _ _ = Nothing +internalFromList _ = Nothing + +-- | Extract the sequence of integers. +internalToList :: InternalAddress -> [Integer] +internalToList (InternalAddress xs) = xs + +-- | Construct an 'InternalAddress' from a kneading sequence. +internalAddress :: Kneading -> Maybe InternalAddress +internalAddress (StarPeriodic [Star]) = Just (InternalAddress [1]) +internalAddress (StarPeriodic v@(One:_)) = Just . InternalAddress . address'per (genericLength v) $ v +internalAddress (Periodic v@(One:_)) = Just . InternalAddress . address'per (genericLength v) $ v +internalAddress k@(Aperiodic (One:_)) = Just . InternalAddress . address'inf . unwrap $ k +internalAddress _ = Nothing + +address'inf :: [Knead] -> [Integer] +address'inf v = address' v + +address'per :: Integer -> [Knead] -> [Integer] +address'per p v = takeWhile (<= p) $ address' v + +address' :: [Knead] -> [Integer] +address' v = address'' 1 [One] + where + address'' sk vk = sk : address'' sk' vk' + where + sk' = (1 +) . genericLength . takeWhile id . zipWith (==) v . cycle $ vk + vk' = genericTake sk' (cycle v) + +-- | A star-periodic kneading sequence's upper and lower associated +-- kneading sequences. +associated :: Kneading -> Maybe (Kneading, Kneading) +associated (StarPeriodic k) = Just (Periodic a, Periodic abar) + where + n = genericLength k + divisors = [ m | m <- [1 .. n], n `mod` m == 0 ] + abar = head . filter (and . zipWith (==) a' . cycle) . map (`genericTake` a') $ divisors + (a, a') = if ((n `elem`) . internalToList) `fmap` internalAddress (Periodic a1) == Just True then (a1, a2) else (a2, a1) + a1 = map (\s -> case s of Star -> Zero ; t -> t) k + a2 = map (\s -> case s of Star -> One ; t -> t) k +associated _ = Nothing + +-- | The upper associated kneading sequence. +upper :: Kneading -> Maybe Kneading +upper = fmap fst . associated + +-- | The lower associated kneading sequence. +lower :: Kneading -> Maybe Kneading +lower = fmap fst . associated + +-- | Angled internal addresses have angles between each integer in an +-- internal address. +data AngledInternalAddress + = Unangled Integer + | Angled Integer Angle AngledInternalAddress + deriving (Read, Show, Eq, Ord, Data, Typeable) + +instance Pretty AngledInternalAddress where + pretty (Unangled n) = pretty n + pretty (Angled n r a) + | r /= 1/2 = pretty n <> char ' ' <> pretty (numerator r) <> char '/' <> pretty (denominator r) <> char ' ' <> pretty a + | otherwise = pretty n <> char ' ' <> pretty a + +-- | Builds a valid 'AngledInternalAddress' from a list, checking the +-- precondition that only the last 'Maybe Angle' should be 'Nothing', +-- and the 'Integer' must be strictly increasing. +angledFromList :: [(Integer, Maybe Angle)] -> Maybe AngledInternalAddress +angledFromList = fromList' 0 + where + fromList' x [(n, Nothing)] | n > x = Just (Unangled n) + fromList' x ((n, Just r) : xs) | n > x && 0 < r && r < 1 = Angled n r `fmap` fromList' n xs + fromList' _ _ = Nothing + +unsafeAngledFromList :: [(Integer, Maybe Angle)] -> AngledInternalAddress +unsafeAngledFromList = fromList' 0 + where + fromList' x [(n, Nothing)] | n > x = Unangled n + fromList' x ((n, Just r) : xs) | n > x && 0 < r && r < 1 = Angled n r (fromList' n xs) + fromList' _ _ = error "Fractal.Mandelbrot.Address.unsafeAngledFromList" + +-- | Convert an 'AngledInternalAddress' to a list. +angledToList :: AngledInternalAddress -> [(Integer, Maybe Angle)] +angledToList (Unangled n) = [(n, Nothing)] +angledToList (Angled n r a) = (n, Just r) : angledToList a + +denominators :: InternalAddress -> Kneading -> [Integer] +denominators a v = denominators' (internalToList a) + where + denominators' (s0:ss@(s1:_)) = + let rr = r s0 s1 + in (((s1 - rr) `div` s0) + if s0 `elem` takeWhile (<= s0) (orbit p rr) then 1 else 2) : denominators' ss + denominators' _ = [] + r s s' = case s' `mod` s of + 0 -> s + t -> t + p = rho v + +numerators :: Angle -> InternalAddress -> [Integer] -> [Integer] +numerators r a qs = zipWith num (internalToList a) qs + where + num s q = genericLength . filter (<= r) . map (genericIndex rs) $ [0 .. q - 2] + where + rs = iterate (foldr (.) id . genericReplicate s $ double) r + +-- | The angled internal address corresponding to an external angle. +angledInternalAddress :: Angle -> Maybe AngledInternalAddress +angledInternalAddress r0 = do + let r = wrap r0 + k = kneading r + i <- internalAddress k + let d = denominators i k + n = numerators r i d + return . unsafeAngledFromList . zip (internalToList i) . (++ [Nothing]) . map Just . zipWith (%) n $ d + +-- | Discard angle information from an internal address. +stripAngles :: AngledInternalAddress -> InternalAddress +stripAngles = InternalAddress . map fst . angledToList + +-- | The pair of external angles whose rays land at the root of the +-- hyperbolic component described by the angled internal address. +externalAngles :: AngledInternalAddress -> Maybe (Rational, Rational) +externalAngles = externalAngles' 1 (0, 1) + +externalAngles' :: Integer -> (Rational, Rational) -> AngledInternalAddress -> Maybe (Rational, Rational) +externalAngles' p0 lohi a0@(Unangled p) + | p0 /= p = case wakees lohi p of + [lh] -> externalAngles' p lh a0 + _ -> Nothing + | otherwise = Just lohi +externalAngles' p0 lohi a0@(Angled 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 + | otherwise = error "Fractal.Mandelbrot.Address.gaps !(n >= 0)" + 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 (kneading r) == Just 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) + accumulate _ _ = error "Fractal.Mandelbrot.Address.gaps !even" + 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 +chunk2 _ = error "Fractal.Mandelbrot.Address.chunk2 !even" + +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 + +-- | Parse an angled internal address, accepting some unambiguous +-- abbreviations. +parse :: String -> Maybe AngledInternalAddress +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 AngledInternalAddress +parser = do + ts <- pTokens + accum 1 ts + where + accum p [] = return $ Unangled p + accum _ [Number n] = return $ Unangled n + accum _ (Number n : ts@(Number _ : _)) = do + a <- accum n ts + return $ Angled n (1%2) a + accum _ (Number n : Fraction t b : ts) = do + a <- accum (n * b) ts + return $ Angled n (t%b) a + accum p (Fraction t b : ts) = do + a <- accum (p * b) ts + return $ Angled p (t % b) a + +pTokens :: Parse [Token] +pTokens = do + _ <- pOptionalSpace + ts <- pToken `sepEndBy` pSpace + 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 " ") diff --git a/Fractal/RUFF/Mandelbrot/Image.hs b/Fractal/RUFF/Mandelbrot/Image.hs new file mode 100644 --- /dev/null +++ b/Fractal/RUFF/Mandelbrot/Image.hs @@ -0,0 +1,118 @@ +{-# LANGUAGE BangPatterns, DeriveDataTypeable #-} +{- | +Module : Fractal.RUFF.Mandelbrot.Image +Copyright : (c) Claude Heiland-Allen 2011 +License : BSD3 + +Maintainer : claudiusmaximus@goto10.org +Stability : unstable +Portability : portable + +Generic functions to render images. + +-} + +module Fractal.RUFF.Mandelbrot.Image + ( simpleImage, complexImage, imageLoop, coordinates, ascii, unicode + ) where + +import Control.Monad.ST (ST) +import Data.Array.ST (newArray, writeArray, runSTUArray) +import Data.STRef (STRef, newSTRef, readSTRef, writeSTRef) +import Data.Array.Unboxed (UArray, (!), bounds, range) + +import Fractal.RUFF.Types.Complex (Complex((:+))) +import Fractal.RUFF.Types.Tuple (Tuple2(Tuple2)) +import Fractal.RUFF.Mandelbrot.Iterate (iterates, initial, Mode(Simple, DistanceEstimate), Iterate(), Output(OutSimple, OutDistanceEstimate), escapeTime, distanceEstimate, finalAngle, outUser) + +-- | Render an image with the 'Simple' algorithm. The iteration count is +-- doubled until the image is good enough, or the fixed maximum iteration +-- count is reached. +-- +-- > putStr . unicode $ simpleImage 100 100 ((-1.861):+0) (0.001) 1000000000 +simpleImage :: (Ord r, Floating r) => Int {- ^ width -} -> Int {- ^ height -} -> Complex r {- ^ center -} -> r {- ^ radius -} -> Int {- ^ max iterations -} -> UArray (Int, Int) Bool {- ^ image -} +{-# INLINABLE simpleImage #-} +simpleImage width height c0 r0 n0 = runSTUArray $ do + a <- newArray bs True + s <- newSTRef (0 :: Int) + imageLoop s a n0 0 False 64 i0s (out s a) + where + (bs, cs) = coordinates width height c0 r0 + i0s = map (uncurry $ initial Simple) cs + out s a (OutSimple{ outUser = Tuple2 j i }) = do + writeArray a (j, i) False + modifySTRef' s (+ 1) + out _ _ _ = return () + +-- | Render an image with the 'DistanceEstimate' algorithm. The iteration count is +-- doubled until the image is good enough, or the fixed maximum iteration +-- count is reached. The output values are converted to 'Float'. +complexImage :: (Ord r, Real r, Floating r) => Int {- ^ width -} -> Int {- ^ height -} -> Complex r {- ^ center -} -> r {- ^ radius -} -> Int {- ^ max iterations -} -> UArray (Int, Int, Int) Float {- ^ image -} +{-# INLINABLE complexImage #-} +complexImage !width !height !c0 !r0 !n0 = runSTUArray $ do + a <- newArray bs (-1) + s <- newSTRef (0 :: Int) + imageLoop s a n0 0 False 64 i0s (out s a) + where + bs = ((jlo,ilo,0), (jhi,ihi,2)) + (((jlo,ilo),(jhi,ihi)), cs) = coordinates width height c0 r0 + i0s = map (uncurry $ initial DistanceEstimate) cs + out !s !a (OutDistanceEstimate{ escapeTime = et, distanceEstimate = de, finalAngle = fa, outUser = Tuple2 j i }) = {-# SCC "complexImage.out" #-} do + writeArray a (j, i, 0) (realToFrac et) + writeArray a (j, i, 1) (realToFrac de) + writeArray a (j, i, 2) (realToFrac fa) + modifySTRef' s (+ 1) + out _ _ _ = return () + +-- | Image rendering loop. +imageLoop :: (Ord r, Floating r) => STRef s Int {- ^ escapees -} -> a {- ^ output array -} -> Int {- ^ max iterations -} -> Int {- ^ iterations -} -> Bool {- ^ prior escapees -} -> Int {- ^ iterations this phase -} -> [Iterate u r] {- ^ iterates -} -> (Output u r -> ST s ()) {- ^ output callback -} -> ST s a {- ^ output array as given -} +{-# INLINABLE imageLoop #-} +imageLoop s a !n0 !n1 !f1 !m1 is1 out = loop f1 n1 m1 is1 + where + loop !f !n !m is = do + writeSTRef s 0 + is' <- iterates m is out + o <- readSTRef s + if null is || (f && o == 0) || n > n0 then return a else loop (f || o > 0) (n + m) (m * 2) is' + +-- | The parameter plane coordinates for an image, with bounds. +coordinates :: (Ord r, Floating r) => Int {- ^ width -} -> Int {- ^ height -} -> Complex r {- ^ center -} -> r {- ^ radius -} -> (((Int,Int),(Int,Int)), [(Tuple2 Int Int, Complex r)]) {- ^ (bounds, coords) -} +{-# INLINABLE coordinates #-} +coordinates !width !height !(c0r :+ c0i) !r0 = (bs, cs) + where + bs = ((0, 0), (height - 1, width - 1)) + cs = [ (Tuple2 j i, c) + | (j,i) <- range bs + , let y = (fromIntegral j - h) / h + , let x = (fromIntegral i - w) / w + , let ci = c0i + r0 * y + , let cr = c0r + r0 * x + , let c = cr :+ ci + ] + w = fromIntegral $ width `div` 2 + h = fromIntegral $ height `div` 2 + +-- | Convert a bit array to ascii graphics. +ascii :: UArray (Int, Int) Bool {- ^ image -} -> String {- ^ ascii -} +ascii a = unlines . map concat $ [ [ b (a ! (j, i)) | i <- [ ilo .. ihi ] ] | j <- [ jhi, jhi - 1 .. jlo ] ] + where + ((jlo, ilo), (jhi, ihi)) = bounds a + b False = " " + b True = "##" + +-- | Convert a bit array to unicode block graphics. +unicode :: UArray (Int, Int) Bool {- ^ image -} -> String {- ^ unicode -} +unicode a = unlines [ [ b (a ! (j, i)) (a ! (j - 1, i)) | i <- [ ilo .. ihi ] ] | j <- [ jhi, jhi - 2 .. jlo ] ] + where + ((jlo, ilo), (jhi, ihi)) = bounds a + b False False = ' ' + b True False = '\x2580' + b False True = '\x2584' + b True True = '\x2588' + +-- | Strict version of 'modifySTRef'. +modifySTRef' :: STRef s a -> (a -> a) -> ST s () +{-# INLINABLE modifySTRef' #-} +modifySTRef' s f = do + x <- readSTRef s + writeSTRef s $! f x diff --git a/Fractal/RUFF/Mandelbrot/Iterate.hs b/Fractal/RUFF/Mandelbrot/Iterate.hs new file mode 100644 --- /dev/null +++ b/Fractal/RUFF/Mandelbrot/Iterate.hs @@ -0,0 +1,116 @@ +{-# LANGUAGE BangPatterns, DeriveDataTypeable #-} +{- | +Module : Fractal.RUFF.Mandelbrot.Iterate +Copyright : (c) Claude Heiland-Allen 2011 +License : BSD3 + +Maintainer : claudiusmaximus@goto10.org +Stability : unstable +Portability : portable + +Generic functions to iterate points. + +-} + +module Fractal.RUFF.Mandelbrot.Iterate + ( Mode (..) + , Iterate (..) + , Output (..) + , initial + , iterate + , iterates + ) where + +import Data.Data (Data) +import Data.Typeable (Typeable) +import Prelude hiding (iterate) + +import Fractal.RUFF.Types.Complex (Complex(..), phase) + +-- | Iteration mode. +data Mode = Simple | EscapeTime | DistanceEstimate + deriving (Read, Show, Eq, Ord, Enum, Bounded, Data, Typeable) + +-- | Iteration state. +data Iterate u r + = IterSimple{ itc, itz :: !(Complex r), iterUser :: !u } + | IterEscapeTime{ itc, itz :: !(Complex r), itn :: !Int, iterUser :: !u } + | IterDistanceEstimate{ itc, itz, itdz :: !(Complex r), itn :: !Int, iterUser :: !u } + deriving (Read, Show, Eq, Ord, Data, Typeable) + +-- | Iteration initial state. +initial :: Num r => Mode -> u -> Complex r -> Iterate u r +{-# INLINABLE initial #-} +initial Simple u c = IterSimple + { itc = c, itz = 0 :+ 0, iterUser = u } +initial EscapeTime u c = IterEscapeTime + { itc = c, itz = 0 :+ 0, itn = 0, iterUser = u } +initial DistanceEstimate u c = IterDistanceEstimate + { itc = c, itz = 0 :+ 0, itdz = 0 :+ 0, itn = 0, iterUser = u } + +-- | Iteration output. +data Output u r + = OutSimple{ outUser :: !u } + | OutEscapeTime{ escapeTime, finalAngle :: !r, outUser :: !u } + | OutDistanceEstimate{ escapeTime, finalAngle, distanceEstimate :: !r, outUser :: !u } + deriving (Read, Show, Eq, Ord, Data, Typeable) + +-- | Iteration engine. +iterate :: (Ord r, Floating r) => Int -> Iterate u r -> Either (Iterate u r) (Output u r) +{-# INLINABLE iterate #-} +iterate n i@(IterSimple{ itc = cr :+ ci, itz = z0, iterUser = u }) = go 0 z0 + where + go !m !z@(zr :+ zi) + | m < n = let !zrr = zr * zr + !zii = zi * zi + !zri = zr * zi + !e = zrr + zii > 4 + in if e then Right (OutSimple{ outUser = u}) + else go (m + 1) ((zrr - zii + cr) :+ (2 * zri + ci)) + | otherwise = Left (i{ itz = z }) +iterate n i@(IterEscapeTime{ itc = cr :+ ci, itz = z0, itn = n0, iterUser = u }) = go 0 z0 + where + er = 65536 + er2 = er * er + log2 = log 2 + go !m !z@(zr :+ zi) + | m < n = let !zrr = zr * zr + !zii = zi * zi + !zri = zr * zi + !zz = zrr + zii + !e = zz > er2 + in if e then Right (OutEscapeTime + { escapeTime = fromIntegral (n0 + m) + (log (log er) - log (log zz / 2)) / log2 + , finalAngle = phase z + , outUser = u}) + else go (m + 1) ((zrr - zii + cr) :+ (2 * zri + ci)) + | otherwise = Left (i{ itz = z, itn = n0 + n }) +iterate !n !i@(IterDistanceEstimate{ itc = cr :+ ci, itz = z0, itdz = dz0, itn = n0, iterUser = u }) = go 0 z0 dz0 + where + er = 65536 + er2 = er * er + log2 = log 2 + go !m !z@(zr :+ zi) !dz@(dzr :+ dzi) + | m < n = let !zrr = zr * zr + !zii = zi * zi + !zri = zr * zi + !zz = zrr + zii + !e = zz > er2 + !zdzr = zr * dzr - zi * dzi + !zdzi = zr * dzi + zi * dzr + !dzdz = dzr * dzr + dzi * dzi + in if e then Right (OutDistanceEstimate + { escapeTime = fromIntegral (n0 + m) + (log (log er) - log (log zz / 2)) / log2 + , finalAngle = phase z + , distanceEstimate = log zz * sqrt (zz / dzdz) + , outUser = u}) + else go (m + 1) ((zrr - zii + cr) :+ (2 * zri + ci)) ((2 * zdzr + 1) :+ (2 * zdzi)) + | otherwise = Left (i{ itz = z, itdz = dz, itn = n0 + n }) + +-- | Iterate over a list. +iterates :: (Functor m, Monad m, Ord r, Floating r) => Int -> [Iterate u r] -> (Output u r -> m ()) -> m [Iterate u r] +{-# INLINABLE iterates #-} +iterates _ [] _ = return [] +iterates n (x:xs) out = case iterate n x of + Right o -> out o >> iterates n xs out + Left y -> (y:) `fmap` iterates n xs out diff --git a/Fractal/RUFF/Mandelbrot/Ray.hs b/Fractal/RUFF/Mandelbrot/Ray.hs new file mode 100644 --- /dev/null +++ b/Fractal/RUFF/Mandelbrot/Ray.hs @@ -0,0 +1,46 @@ +{-# LANGUAGE BangPatterns, DeriveDataTypeable #-} +{- | +Module : Fractal.RUFF.Mandelbrot.Ray +Copyright : (c) Claude Heiland-Allen 2011 +License : BSD3 + +Maintainer : claudiusmaximus@goto10.org +Stability : unstable +Portability : portable + +External angles define external rays which can be traced back from +the circle at infinity to the boundary of the Mandelbrot Set. + +The algorithm is based on Tomoki Kawahira's paper +/An algorithm to draw external rays of the Mandelbrot set/ +. +-} + +module Fractal.RUFF.Mandelbrot.Ray (externalRay) where + +import Fractal.RUFF.Types.Complex (Complex(), magnitude, mkPolar) +import Fractal.RUFF.Mandelbrot.Address (Angle, double) + +-- | Compute the external ray for an external angle with a given +-- accuracy, sharpness and starting radius. For example: +-- +-- > externalRay 1e-10 8 (2**24) (1/3) +-- +externalRay :: (Ord r, Floating r) => r -> Int -> r -> Angle -> [Complex r] +externalRay epsilon sharpness radius angle = map fst . iterate step $ (mkPolar radius (2 * pi * fromRational angle), (0, 0)) + where + -- step :: (NearZero r, Floating r) => (Complex r, (Int, Int)) -> (Complex r, (Int, Int)) + step (!c, (!k0, !j0)) + | j > sharpness = step (c, (k0 + 1, 0)) + | otherwise = (n 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 * pi * fromRational (iterate double angle !! k0)) + n !z = let d = (cc - t) / dd in if not (magnitude d > epsilon) then z else n (z - d) + 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) diff --git a/Fractal/RUFF/Types/Complex.hs b/Fractal/RUFF/Types/Complex.hs new file mode 100644 --- /dev/null +++ b/Fractal/RUFF/Types/Complex.hs @@ -0,0 +1,78 @@ +{-# LANGUAGE DeriveDataTypeable #-} +{- | +Module : Fractal.RUFF.Types.Complex +Copyright : (c) Claude Heiland-Allen 2011 +License : BSD3 + +Maintainer : claudiusmaximus@goto10.org +Stability : unstable +Portability : portable + +Complex numbers without the 'RealFloat' constraint. +-} + +module Fractal.RUFF.Types.Complex + ( Complex((:+)), cis, mkPolar + , realPart, imagPart, conjugate + , magnitude, phase, polar + ) where + +import Data.Data (Data) +import Data.Typeable (Typeable) + +-- | Complex number type without the 'RealFloat' constraint. +data Complex r = !r :+ !r + deriving (Read, Show, Eq, Ord, Data, Typeable) + +instance Num r => Num (Complex r) where + (x :+ y) + (u :+ v) = (x + u) :+ (y + v) + (x :+ y) - (u :+ v) = (x - u) :+ (y - v) + (x :+ y) * (u :+ v) = (x * u - y * v) :+ (x * v + y * u) + negate (x :+ y) = negate x :+ negate y + abs = error "Fractal.Types.Complex.Num.abs" + signum = error "Fractal.Types.Complex.Num.signum" + fromInteger n = fromInteger n :+ 0 + +instance Fractional r => Fractional (Complex r) where + (x :+ y) / (u :+ v) = ((x * u + y * v) / d) :+ ((y * u - x * v) / d) where d = u * u + v * v + fromRational r = fromRational r :+ 0 + +-- | Extract the real part. +realPart :: Complex r -> r +realPart (r :+ _) = r + +-- | Extract the imaginary part. +imagPart :: Complex r -> r +imagPart (_ :+ i) = i + +-- | Complex conjugate. +conjugate :: Num r => Complex r -> Complex r +conjugate (r :+ i) = r :+ negate i + +-- | Complex phase. +phase :: (Ord r, Floating r) => Complex r -> r +phase (r :+ i) + | r > 0 && i > 0 = atan ( i / r) + | r > 0 && i < 0 = - atan (abs i / r) + | r < 0 && i > 0 = pi - atan ( i / abs r) + | r < 0 && i < 0 = atan (abs i / abs r) - pi + | i > 0 = pi / 2 + | i < 0 = - pi / 2 + | r < 0 = pi + | otherwise = 0 + +-- | Complex magnitude. +magnitude :: Floating r => Complex r -> r +magnitude (r :+ i) = sqrt $ r * r + i * i + +-- | Complex number with the given magnitude and phase. +mkPolar :: Floating r => r -> r -> Complex r +mkPolar r t = (r * cos t) :+ (r * sin t) + +-- | Complex number with magnitude 1 and the given phase. +cis :: Floating r => r -> Complex r +cis t = cos t :+ sin t + +-- | Convert to polar form. +polar :: (Ord r, Floating r) => Complex r -> (r, r) +polar z = (magnitude z, phase z) diff --git a/Fractal/RUFF/Types/Tuple.hs b/Fractal/RUFF/Types/Tuple.hs new file mode 100644 --- /dev/null +++ b/Fractal/RUFF/Types/Tuple.hs @@ -0,0 +1,23 @@ +{-# LANGUAGE DeriveDataTypeable #-} +{- | +Module : Fractal.RUFF.Types.Tuple +Copyright : (c) Claude Heiland-Allen 2011 +License : BSD3 + +Maintainer : claudiusmaximus@goto10.org +Stability : unstable +Portability : portable + +Strict tuples. +-} + +module Fractal.RUFF.Types.Tuple + ( Tuple2(..) + ) where + +import Data.Data (Data) +import Data.Typeable (Typeable) + +-- | Strict 'Tuple2' type. +data Tuple2 l r = Tuple2 !l !r + deriving (Read, Show, Eq, Ord, Data, Typeable) diff --git a/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -0,0 +1,30 @@ +Copyright (c)2011, Claude Heiland-Allen + +All rights reserved. + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions are met: + + * Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. + + * Redistributions in binary form must reproduce the above + copyright notice, this list of conditions and the following + disclaimer in the documentation and/or other materials provided + with the distribution. + + * Neither the name of Claude Heiland-Allen nor the names of other + contributors may be used to endorse or promote products derived + from this software without specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR +A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT +OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, +SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT +LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, +DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY +THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT +(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE +OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. diff --git a/Setup.hs b/Setup.hs new file mode 100644 --- /dev/null +++ b/Setup.hs @@ -0,0 +1,2 @@ +import Distribution.Simple +main = defaultMain diff --git a/ruff.cabal b/ruff.cabal new file mode 100644 --- /dev/null +++ b/ruff.cabal @@ -0,0 +1,24 @@ +Name: ruff +Version: 0.1 +Synopsis: relatively useful fractal functions +Description: + A library for analysis and exploration of fractals. This initial + version provides angled internal addresses, external ray tracing, + and iterations for images of the Mandelbrot Set. + +Homepage: https://gitorious.org/ruff +License: BSD3 +License-file: LICENSE +Author: Claude Heiland-Allen +Maintainer: claudiusmaximus@goto10.org +Copyright: (c) 2011 Claude Heiland-Allen +Category: Math +Build-type: Simple + +Cabal-version: >=1.2 + +Library + Exposed-modules: Fractal.RUFF.Mandelbrot.Address, Fractal.RUFF.Mandelbrot.Iterate, Fractal.RUFF.Mandelbrot.Image, Fractal.RUFF.Mandelbrot.Ray, Fractal.RUFF.Types.Complex, Fractal.RUFF.Types.Tuple + Build-depends: base >= 3 && < 6, array >= 0.3 && < 0.4, mtl >= 2 && < 3, parsec >= 3.1 && < 3.2, wl-pprint-text >= 1 && < 2 + GHC-Options: -Wall -O2 -funbox-strict-fields + GHC-Prof-Options: -prof -auto-all -caf-all