ruff (empty) → 0.1
raw patch · 9 files changed
+844/−0 lines, 9 filesdep +arraydep +basedep +mtlsetup-changed
Dependencies added: array, base, mtl, parsec, wl-pprint-text
Files
- Fractal/RUFF/Mandelbrot/Address.hs +407/−0
- Fractal/RUFF/Mandelbrot/Image.hs +118/−0
- Fractal/RUFF/Mandelbrot/Iterate.hs +116/−0
- Fractal/RUFF/Mandelbrot/Ray.hs +46/−0
- Fractal/RUFF/Types/Complex.hs +78/−0
- Fractal/RUFF/Types/Tuple.hs +23/−0
- LICENSE +30/−0
- Setup.hs +2/−0
- ruff.cabal +24/−0
+ Fractal/RUFF/Mandelbrot/Address.hs view
@@ -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)/+<http://arxiv.org/abs/math/9411238v2>.+-}++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 " ")
+ Fractal/RUFF/Mandelbrot/Image.hs view
@@ -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
+ Fractal/RUFF/Mandelbrot/Iterate.hs view
@@ -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
+ Fractal/RUFF/Mandelbrot/Ray.hs view
@@ -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/+<http://www.math.nagoya-u.ac.jp/~kawahira/programs/mandel-exray.pdf>.+-}++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)
+ Fractal/RUFF/Types/Complex.hs view
@@ -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)
+ Fractal/RUFF/Types/Tuple.hs view
@@ -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)
+ LICENSE view
@@ -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.
+ Setup.hs view
@@ -0,0 +1,2 @@+import Distribution.Simple+main = defaultMain
+ ruff.cabal view
@@ -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