diff --git a/Fractal/RUFF/Mandelbrot/Address.hs b/Fractal/RUFF/Mandelbrot/Address.hs --- a/Fractal/RUFF/Mandelbrot/Address.hs +++ b/Fractal/RUFF/Mandelbrot/Address.hs @@ -19,12 +19,13 @@ -} 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 + ( Angle, double, wrap, prettyAngle, prettyAngles + , Knead(..), kneadChar + , Kneading(..), prettyKneading, kneading, period, unwrap, associated, upper, lower + , InternalAddress(..), prettyInternalAddress, internalAddress, internalFromList, internalToList + , AngledInternalAddress(..), prettyAngledInternalAddress, angledInternalAddress, angledFromList, angledToList + , externalAngles, stripAngles, splitAddress, joinAddress, addressPeriod + , parseAngle, parseAngles, parseKnead, parseKneading, parseInternalAddress, parseAngledInternalAddress ) where import Data.Data (Data()) @@ -32,15 +33,25 @@ import Control.Monad (guard) import Control.Monad.Identity (Identity()) import Data.Char (digitToInt) -import Data.List (genericDrop, genericIndex, genericLength, genericReplicate, genericSplitAt, genericTake) +import Data.Bits (testBit) +import Data.List (genericDrop, genericIndex, genericLength, genericReplicate, genericSplitAt, genericTake, foldl') 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 +-- | Convert to human readable form. +prettyAngle :: Angle -> String +prettyAngle a = show (numerator a) ++ " / " ++ show (denominator a) + +-- | Convert to human readable form. +prettyAngles :: [Angle] -> String +prettyAngles [] = "" +prettyAngles [a] = show (numerator a) ++ "/" ++ show (denominator a) +prettyAngles (a:as) = show (numerator a) ++ "/" ++ show (denominator a) ++ " " ++ prettyAngles as + -- | Wrap an angle into [0, 1). wrap :: Angle -> Angle wrap a @@ -53,6 +64,50 @@ double :: Angle -> Angle double a = wrap (2 * a) +-- | Binary representation of a (pre-)periodic angle. +type BinAngle = ([Bool], [Bool]) + +-- | Convert an angle from binary representation. +unbinary :: BinAngle -> Angle +unbinary (pre, per) + | n == 0 = bits pre % (2 ^ m) + | otherwise = (bits pre % (2 ^ m)) + (bits per % (2 ^ m * (2 ^ n - 1))) + where + m = length pre + n = length per + +-- | Convert a list of bits to an integer. +bits :: [Bool] -> Integer +bits = foldl' (\ a b -> 2 * a + if b then 1 else 0) 0 + +-- | Convert an angle to binary representation. +binary :: Angle -> BinAngle +binary a + | a == 0 = ([], []) + | even (denominator a) = + let (pre, per) = binary (double a) + b = a >= 1/2 + in (b:pre, per) + | otherwise = + let (t, p) = head . dropWhile ((1 /=) . denominator . fst) . map (\q -> (a * (2^q - 1), q)) $ [ (1 :: Int) ..] + s = numerator t + n = fromIntegral p + per = [ s `testBit` i | i <- [n - 1, n - 2 .. 0] ] + in ([], per) + +-- | Tuning transformation for binary represented periodic angles. +-- Probably only valid for angle pairs presenting ray pairs. +btune :: BinAngle -> (BinAngle, BinAngle) -> BinAngle +btune (tpre, tper) (([], per0), ([], per1)) = (concatMap f tpre, concatMap f tper) + where + f False = per0 + f True = per1 +btune _ _ = error "btune: can't handle pre-periods" + +-- | Tuning transformation for angles. +tune :: Angle -> (Angle, Angle) -> Angle +tune t (t0, t1) = unbinary $ btune (binary t) (binary t0, binary t1) + -- | Elements of kneading sequences. data Knead = Zero @@ -60,17 +115,13 @@ | Star deriving (Read, Show, Eq, Ord, Enum, Bounded, Data, Typeable) -instance Pretty Knead where - pretty = char . kneadChar - prettyList = pretty . map kneadChar - +-- | Knead character representation. 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. +-- | Kneading sequences. Note that the 'Aperiodic' case has an infinite list. data Kneading = Aperiodic [Knead] | PrePeriodic [Knead] [Knead] @@ -78,11 +129,12 @@ | 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) +-- | Kneading sequence as a string. The 'Aperiodic' case is truncated arbitrarily. +prettyKneading :: Kneading -> String +prettyKneading (Aperiodic ks) = map kneadChar (take 17 ks) ++ "..." +prettyKneading (PrePeriodic us vs) = map kneadChar us ++ "(" ++ map kneadChar vs ++ ")" +prettyKneading (StarPeriodic vs) = "(" ++ map kneadChar vs ++ ")" +prettyKneading (Periodic vs) = "(" ++ map kneadChar vs ++ ")" -- | The kneading sequence for an external angle. kneading :: Angle -> Kneading @@ -143,10 +195,11 @@ 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) +-- | Internal address as a string. +prettyInternalAddress :: InternalAddress -> String +prettyInternalAddress (InternalAddress []) = error "Fractal.Mandelbrot.Address.InternalAddress.pretty" +prettyInternalAddress (InternalAddress [x]) = show x +prettyInternalAddress (InternalAddress (x:ys)) = show x ++ " " ++ prettyInternalAddress (InternalAddress ys) -- | Construct a valid 'InternalAddress', checking the precondition. internalFromList :: [Integer] -> Maybe InternalAddress @@ -202,7 +255,7 @@ -- | The lower associated kneading sequence. lower :: Kneading -> Maybe Kneading -lower = fmap fst . associated +lower = fmap snd . associated -- | Angled internal addresses have angles between each integer in an -- internal address. @@ -211,11 +264,12 @@ | 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 +-- | Angled internal address as a string. +prettyAngledInternalAddress :: AngledInternalAddress -> String +prettyAngledInternalAddress (Unangled n) = show n +prettyAngledInternalAddress (Angled n r a) + | r /= 1/2 = show n ++ " " ++ show (numerator r) ++ "/" ++ show (denominator r) ++ " " ++ prettyAngledInternalAddress a + | otherwise = show n ++ " " ++ prettyAngledInternalAddress a -- | Builds a valid 'AngledInternalAddress' from a list, checking the -- precondition that only the last 'Maybe Angle' should be 'Nothing', @@ -268,6 +322,36 @@ n = numerators r i d return . unsafeAngledFromList . zip (internalToList i) . (++ [Nothing]) . map Just . zipWith (%) n $ d +-- | Split an angled internal address at the last island. +splitAddress :: AngledInternalAddress -> (AngledInternalAddress, [Angle]) +splitAddress a = + let (ps0, rs0) = unzip $ angledToList a + ps1 = reverse ps0 + rs1 = reverse (Nothing : init rs0) + prs1 = zip ps1 rs1 + f ((p, Just r):qrs@((q, _):_)) acc + | p == denominator r * q = f qrs (r : acc) + f prs acc = g prs acc + g prs acc = + let (ps2, rs2) = unzip prs + ps3 = reverse ps2 + rs3 = reverse (Nothing : init rs2) + prs3 = zip ps3 rs3 + aa = unsafeAngledFromList prs3 + in (aa, acc) + in f prs1 [] + +-- | The inverse of 'splitAddress'. +joinAddress :: AngledInternalAddress -> [Angle] -> AngledInternalAddress +joinAddress (Unangled p) [] = Unangled p +joinAddress (Unangled p) (r:rs) = Angled p r (joinAddress (Unangled $ p * denominator r) rs) +joinAddress (Angled p r a) rs = Angled p r (joinAddress a rs) + +-- | The period of an angled internal address. +addressPeriod :: AngledInternalAddress -> Integer +addressPeriod (Unangled p) = p +addressPeriod (Angled _ _ a) = addressPeriod a + -- | Discard angle information from an internal address. stripAngles :: AngledInternalAddress -> InternalAddress stripAngles = InternalAddress . map fst . angledToList @@ -288,6 +372,7 @@ [lh] -> externalAngles' p lh a0 _ -> Nothing | otherwise = do +{- let num = numerator r den = denominator r q = p * den @@ -300,11 +385,23 @@ i <- genericElemIndex num nums lh <- safeGenericIndex ws (i :: Integer) externalAngles' q lh a - +-} + let num = numerator r + den = denominator r + ws = wakees (0, 1) den + 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 + (l,h) <- safeGenericIndex ws (i :: Integer) + externalAngles' (p * den) (if p > 1 then (tune l lohi, tune h lohi) else (l, h)) a wakees :: (Rational, Rational) -> Integer -> [(Rational, Rational)] wakees (lo, hi) q = let gaps (l, h) n | n == 0 = [(l, h)] +-- | h - l < 1 % (2 ^ n - 1) = [(l, h)] | n > 0 = let gs = gaps (l, h) (n - 1) cs = candidates n gs in accumulate cs gs @@ -346,15 +443,54 @@ | i > 0 = safeGenericIndex xs (i - 1) | otherwise = Just x +-- | Parse an angle. +parseAngle :: String -> Maybe Angle +parseAngle s = case runP pFraction () "" s of + Left _ -> Nothing + Right f -> Just (unFraction f) + +-- | Parse a list of angles. +parseAngles :: String -> Maybe [Angle] +parseAngles s = case runP (many pFraction) () "" s of + Left _ -> Nothing + Right fs -> Just (map unFraction fs) + +-- | Parse a kneading element. +parseKnead :: String -> Maybe Knead +parseKnead s = case runP pKnead () "" s of + Left _ -> Nothing + Right k -> Just k + +-- | Parse a non-aperiodic kneading sequence. +parseKneading :: String -> Maybe Kneading +parseKneading s = case runP pKneading () "" s of + Left _ -> Nothing + Right ks -> Just ks + +-- | Parse an internal address. +parseInternalAddress :: String -> Maybe InternalAddress +parseInternalAddress s = case runP (many pNumber) () "" s of + Left _ -> Nothing + Right ns -> internalFromList (map unNumber ns) + -- | Parse an angled internal address, accepting some unambiguous -- abbreviations. -parse :: String -> Maybe AngledInternalAddress -parse s = case runP parser () "" s of +parseAngledInternalAddress :: String -> Maybe AngledInternalAddress +parseAngledInternalAddress s = case runP parser () "" s of Left _ -> Nothing Right a -> Just a data Token = Number Integer | Fraction Integer Integer +unFraction :: Token -> Angle +unFraction (Fraction t b) = t % b +unFraction _ = error "Fractal.Mandelbrot.Address.unFraction" + +unNumber :: Token -> Integer +unNumber (Number n) = n +unNumber _ = error "Fractal.Mandelbrot.Address.unNumber" + + type Parse t = ParsecT String () Identity t parser :: Parse AngledInternalAddress @@ -405,3 +541,17 @@ pOptionalSpace :: Parse [String] pOptionalSpace = many (string " ") + +pKnead :: Parse Knead +pKnead = choice [ string "0" >> return Zero, string "1" >> return One, string "*" >> return Star ] + +pKneading :: Parse Kneading +pKneading = do + pre <- many pKnead + _ <- string "(" + per <- many1 pKnead + _ <- string ")" + return $ case (null pre, last per) of + (False, _) -> PrePeriodic pre per + (True, Star) -> StarPeriodic per + _ -> Periodic per diff --git a/Fractal/RUFF/Mandelbrot/Image.hs b/Fractal/RUFF/Mandelbrot/Image.hs --- a/Fractal/RUFF/Mandelbrot/Image.hs +++ b/Fractal/RUFF/Mandelbrot/Image.hs @@ -8,20 +8,25 @@ Stability : unstable Portability : portable -Generic functions to render images. +Generic (slow) functions to render images. -} module Fractal.RUFF.Mandelbrot.Image ( simpleImage, complexImage, imageLoop, coordinates, ascii, unicode + , Channel(..), Coordinates, border ) 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 Data.Array.Unboxed (UArray, (!), bounds, range, amap, ixmap) -import Fractal.RUFF.Types.Complex (Complex((:+))) +import Data.Ix (Ix) +import Data.Data (Data) +import Data.Typeable (Typeable) + +import Fractal.RUFF.Types.Complex (Complex((:+)), magnitude) 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) @@ -29,15 +34,14 @@ -- 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 -} +-- > putStr . unicode $ simpleImage (coordinates 100 100 ((-1.861):+0) (0.001)) 1000000000 +simpleImage :: (Ord r, Floating r) => Coordinates r {- ^ coordinates -} -> Int {- ^ max iterations -} -> UArray (Int, Int) Bool {- ^ image -} {-# INLINABLE simpleImage #-} -simpleImage width height c0 r0 n0 = runSTUArray $ do +simpleImage (bs, cs) 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 @@ -47,23 +51,30 @@ -- | 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 -} +-- +-- > putStr . unicode . border $ complexImage (coordinates 100 100 ((-1.861):+0) (0.001)) 1000000000 +complexImage :: (Ord r, Real r, Floating r) => Coordinates r {-^ coordinates -} -> Int {- ^ max iterations -} -> UArray (Int, Int, Channel) Float {- ^ image -} {-# INLINABLE complexImage #-} -complexImage !width !height !c0 !r0 !n0 = runSTUArray $ do +complexImage (((jlo,ilo),(jhi,ihi)), cs) !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 + bs = ((jlo,ilo,minBound), (jhi,ihi,maxBound)) + (_, cx0):(_, cx1):_ = cs + pixelSpacing = magnitude (cx1 - cx0) 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) + writeArray a (j, i, EscapeTime) (realToFrac et) + writeArray a (j, i, DistanceEstimate') (realToFrac (de / pixelSpacing)) + writeArray a (j, i, FinalAngle) (realToFrac fa) modifySTRef' s (+ 1) out _ _ _ = return () +-- | Channels in an image. +data Channel = EscapeTime {- ^ continuous dwell -} | DistanceEstimate' {- ^ normalized to pixel spacing -} | FinalAngle {- ^ in [-pi,pi] -} + deriving (Eq, Ord, Enum, Bounded, Ix, Read, Show, Data, Typeable) + -- | 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 #-} @@ -75,8 +86,11 @@ o <- readSTRef s if null is || (f && o == 0) || n > n0 then return a else loop (f || o > 0) (n + m) (m * 2) is' +-- | Image bounds and coordinates. +type Coordinates r = (((Int,Int),(Int,Int)), [(Tuple2 Int Int, Complex r)]) + -- | 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) -} +coordinates :: (Ord r, Floating r) => Int {- ^ width -} -> Int {- ^ height -} -> Complex r {- ^ center -} -> r {- ^ size -} -> Coordinates r {-# INLINABLE coordinates #-} coordinates !width !height !(c0r :+ c0i) !r0 = (bs, cs) where @@ -84,13 +98,21 @@ cs = [ (Tuple2 j i, c) | (j,i) <- range bs , let y = (fromIntegral j - h) / h - , let x = (fromIntegral i - w) / w + , let x = (fromIntegral i - w) / h , 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 distance estimate image to a near-boundary bit array. +-- The input image must have a DistanceEstimate' channel. +border :: UArray (Int, Int, Channel) Float {- ^ image -} -> UArray (Int, Int) Bool +border a = amap (\x -> x > 0 && x < 1) . ixmap bs (\(j, i) -> (j, i, DistanceEstimate')) $ a + where + ((jlo, ilo, _), (jhi, ihi, _)) = bounds a + bs = ((jlo, ilo), (jhi, ihi)) -- | Convert a bit array to ascii graphics. ascii :: UArray (Int, Int) Bool {- ^ image -} -> String {- ^ ascii -} diff --git a/Fractal/RUFF/Mandelbrot/Iterate.hs b/Fractal/RUFF/Mandelbrot/Iterate.hs --- a/Fractal/RUFF/Mandelbrot/Iterate.hs +++ b/Fractal/RUFF/Mandelbrot/Iterate.hs @@ -8,8 +8,7 @@ Stability : unstable Portability : portable -Generic functions to iterate points. - +Generic (slow) functions to iterate points. -} module Fractal.RUFF.Mandelbrot.Iterate @@ -24,7 +23,6 @@ import Data.Data (Data) import Data.Typeable (Typeable) import Prelude hiding (iterate) - import Fractal.RUFF.Types.Complex (Complex(..), phase) -- | Iteration mode. @@ -36,7 +34,7 @@ = 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) + deriving (Read, Show, Eq, Data, Typeable) -- | Iteration initial state. initial :: Num r => Mode -> u -> Complex r -> Iterate u r diff --git a/Fractal/RUFF/Mandelbrot/Nucleus.hs b/Fractal/RUFF/Mandelbrot/Nucleus.hs new file mode 100644 --- /dev/null +++ b/Fractal/RUFF/Mandelbrot/Nucleus.hs @@ -0,0 +1,88 @@ +{-# LANGUAGE BangPatterns #-} +{- | +Module : Fractal.RUFF.Mandelbrot.Nucleus +Copyright : (c) Claude Heiland-Allen 2011 +License : BSD3 + +Maintainer : claudiusmaximus@goto10.org +Stability : unstable +Portability : portable + +Mu-atom period, nucleus and bond point finding. +-} +module Fractal.RUFF.Mandelbrot.Nucleus (findPeriod, findNucleus, findBond, findInternal) where + +import Data.List (genericIndex) +import Data.Maybe (listToMaybe) +import Fractal.RUFF.Types.Complex (Complex((:+)), mkPolar, magnitude2) + +-- | Given the period and approximate location, successively refine +-- this estimate to a nucleus. +-- +-- The algorithm is based on Robert Munafo's page +-- /Newton-Raphson method/ +-- . +-- +findNucleus :: (Floating r, Fractional r) => Integer {- ^ period -} -> Complex r {- ^ estimate -} -> [Complex r] +findNucleus p g = iterate go g + where + go !c = + let step (!z, !d) = (z * z + c, 2 * z * d + 1) + (zn, dn) = iterate step (0, 0) `genericIndex` p + in c - zn / dn + +-- | Given the period and nucleus, find succesive refinements to the +-- bond point at a given internal angle. +-- +-- The algorithm is based on ideas from +-- . +-- +findBond :: (Floating r, Fractional r) => Integer {- ^ period -} -> Complex r {- ^ nucleus -} -> r {- ^ angle -} -> [Complex r] +findBond p c0 a0 = findInternal p c0 1 a0 + +-- | Given the period and nucleus, find an interior point at a given internal +-- angle and radius in (0,1]. +-- +findInternal :: (Floating r, Fractional r) => Integer {- ^ period -} -> Complex r {- ^ nucleus -} -> r {- ^ radius -} -> r {- ^ angle -} -> [Complex r] +findInternal p c0 r0 a0 = snd `map` iterate go (c0, c0) + where + b0 = mkPolar r0 (2 * pi * a0) + go (!z1, !c1) = + let step (!a, !b, !c, !d, !e) = + ( a * a + c1 + , 2 * a * b + , 2 * (b * b + a * c) + , 2 * a * d + 1 + , 2 * (a * e + b * d) + ) + (an, bn, cn, dn, en) = iterate step (z1, 1, 0, 0, 0) `genericIndex` p + y0 = z1 - an + y1 = b0 - bn + bn1 = bn - 1 + m = bn1 * en - dn * cn + d0 = (y0 * en - dn * y1) / m + d1 = (bn1 * y1 - y0 * cn) / m + in (z1 + d0, c1 + d1) + +-- | Find the period of the lowest period nucleus inside a square. +-- +-- The algorithm is based on Robert Munafo's page, +-- /Finding the Period of a mu-Atom/ +-- . +-- +findPeriod :: (Floating r, Ord r) => Integer {- ^ maximum period -} -> r {- ^ radius -} -> Complex r {- ^ center -} -> Maybe Integer +findPeriod m r c = + let cs = [ c + (r:+r), c + (r:+(-r)), c + ((-r):+(-r)), c + ((-r):+r) ] + zs = iterate (zipWith (\cc z -> z * z + cc) cs) [0,0,0,0] + -- kludge = if r > 0 then 1 else 2 -- fixes space leak from long literal list (CAF?) + in fmap fst . listToMaybe . dropWhile (not . straddlesOrigin . snd) . takeWhile (all ((< 65536) . magnitude2) . snd) . zip [{-kludge + 0 - kludge-} 0 .. m ] $ zs + +straddlesOrigin :: (Ord r, Num r) => [Complex r] -> Bool +straddlesOrigin ps = odd . length . filter id . zipWith positiveReal ps $ (drop 1 ps ++ take 1 ps) + +positiveReal :: (Ord r, Num r) => Complex r -> Complex r -> Bool +positiveReal (u:+v) (x:+y) + | v < 0 && y < 0 = False + | v > 0 && y > 0 = False + | (u * (y - v) - v * (x - u)) * (y - v) > 0 = True + | otherwise = False diff --git a/Fractal/RUFF/Mandelbrot/Ray.hs b/Fractal/RUFF/Mandelbrot/Ray.hs --- a/Fractal/RUFF/Mandelbrot/Ray.hs +++ b/Fractal/RUFF/Mandelbrot/Ray.hs @@ -1,4 +1,4 @@ -{-# LANGUAGE BangPatterns, DeriveDataTypeable #-} +{-# LANGUAGE BangPatterns #-} {- | Module : Fractal.RUFF.Mandelbrot.Ray Copyright : (c) Claude Heiland-Allen 2011 @@ -9,16 +9,15 @@ 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/ -. +the circle at infinity to parameters near the boundary of the Mandelbrot +Set. Conversely, parameters near the boundary of the Mandelbrot Set can +be traced outwards to compute external angles. -} +module Fractal.RUFF.Mandelbrot.Ray (externalRay, externalRayOut) where -module Fractal.RUFF.Mandelbrot.Ray (externalRay) where +import Data.Maybe (fromMaybe) -import Fractal.RUFF.Types.Complex (Complex(), magnitude, mkPolar) +import Fractal.RUFF.Types.Complex (Complex, magnitude2, magnitude, phase, mkPolar) import Fractal.RUFF.Mandelbrot.Address (Angle, double) -- | Compute the external ray for an external angle with a given @@ -26,21 +25,88 @@ -- -- > 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)) +-- The algorithm is based on Tomoki Kawahira's paper +-- /An algorithm to draw external rays of the Mandelbrot set/ +-- . +-- +externalRay :: (Ord r, Floating r) => r {- ^ accuracy -} -> Int {- ^ sharpness -} -> r {- ^ radius -} -> Angle {- ^ external angle -} -> [Complex r] +externalRay accuracy sharpness radius angle = map fst3 . iterate step $ (mkPolar radius (2 * pi * fromRational angle), accuracy * radius, (0, 0)) where + fst3 (x, _, _) = x -- 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)) + step (!c, !epsilon, (!k0, !j0)) + | j > sharpness = step (c, epsilon, (k0 + 1, 0)) + | otherwise = + let c' = n c + epsilon' = accuracy * magnitude (c' - c) + in (c', epsilon', (k0, j0 + 1)) where + epsilon2 = epsilon * epsilon 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) + n !z = let d = (cc - t) / dd in if not (magnitude2 d > epsilon2) 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) + +-- | Compute the external ray outwards from a given parameter value. +-- If the result @rs@ satisfies: +-- +-- > c = last rs +-- > magnitude c > radius +-- +-- then the external angle is given by @t@: +-- +-- > a = phase c / (2 * pi) +-- > t = a - fromIntegral (floor a) +-- +externalRayOut :: (Ord r, Floating r, RealFrac r) + => Int {- ^ iterations -} + -> r {- ^ epsilon -} + -> r {- ^ accuracy -} + -> Int {- ^ sharpness -} + -> r {- ^ radius -} + -> Complex r {- ^ parameter -} + -> [Complex r] +externalRayOut maxIters epsilon accuracy sharpness radius = go (epsilon * epsilon) + where + radius2 = radius * radius + iter !c !n !z + | magnitude2 z > radius2 = Just (n, z) + | n > maxIters = Nothing + | otherwise = iter c (n + 1) (z * z + c) + iterd !c !z !dz !m + | m == 0 = (z, dz) + | otherwise = iterd c (z * z + c) (2 * z * dz + 1) (m - 1) + go !epsilon2 !c = (c :) . fromMaybe [] $ do + (n, z) <- iter c 0 0 + let d = fromIntegral n - logBase 2 (log (magnitude2 z) / log radius2) + d' = d - 1 / fromIntegral sharpness + m = ceiling d' + r = radius ** (2 ** (fromIntegral m - d')) + a = phase z / (2 * pi) + t = a - fromIntegral (floor a :: Int) + k0 = mkPolar r (phase z) + k1 = mkPolar r (pi * t ) + k2 = mkPolar r (pi * (t + 1)) + step !k !c0 = let (f, df) = iterd c0 0 0 m + dc = (f - k) / df + c0' = c0 - dc + in c0 : if not (magnitude2 dc > epsilon2) then [] else step k c0' + steps k = step k c + if m == n + then do + return $ let c' = last $ steps k0 in go (accuracy * magnitude2 (c' - c)) c' + else if m > 0 then do + let (c1, c2) = last (steps k1 `zip` steps k2) + (n1, _) <- iter c1 0 0 + (n2, _) <- iter c2 0 0 + let (c', n') + | magnitude2 (c1 - c) < magnitude2 (c2 - c) = (c1, n1) + | otherwise = (c2, n2) + return $ if n' == m then go (accuracy * magnitude2 (c' - c)) c' else [] + else return [] diff --git a/Fractal/RUFF/Types/Complex.hs b/Fractal/RUFF/Types/Complex.hs --- a/Fractal/RUFF/Types/Complex.hs +++ b/Fractal/RUFF/Types/Complex.hs @@ -14,7 +14,7 @@ module Fractal.RUFF.Types.Complex ( Complex((:+)), cis, mkPolar , realPart, imagPart, conjugate - , magnitude, phase, polar + , magnitude2, magnitude, phase, polar ) where import Data.Data (Data) @@ -22,7 +22,7 @@ -- | Complex number type without the 'RealFloat' constraint. data Complex r = !r :+ !r - deriving (Read, Show, Eq, Ord, Data, Typeable) + deriving (Read, Show, Eq, Data, Typeable) instance Num r => Num (Complex r) where (x :+ y) + (u :+ v) = (x + u) :+ (y + v) @@ -37,6 +37,23 @@ (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 +instance (Ord r, Floating r) => Floating (Complex r) where + pi = pi :+ 0 + exp (x :+ y) = mkPolar (exp x) y + log z = let (r, t) = polar z in log r :+ t + sin (x :+ y) = (sin x * cosh y) :+ (cos x * sinh y) + cos (x :+ y) = (cos x * cosh y) :+ negate (sin x * sinh y) + tan z = sin z / cos z + asin z = negate i * log (i * z + sqrt (1 - z*z)) where i = 0:+1 + acos z = negate i * log (z + sqrt (z*z - 1)) where i = 0:+1 + atan z = 1/2 * i * log ((1 - iz)/(1 + iz)) where i = 0:+1 ; iz = i * z + sinh z = (exp z - exp (-z)) / 2 + cosh z = (exp z + exp (-z)) / 2 + tanh z = let ez2 = exp (2 * z) in (ez2 - 1) / (ez2 + 1) + asinh z = log (z + sqrt (z*z + 1)) + acosh z = log (z + sqrt (z*z - 1)) + atanh z = 1/2 * log ((1 + z) / (1 - z)) + -- | Extract the real part. realPart :: Complex r -> r realPart (r :+ _) = r @@ -49,6 +66,14 @@ conjugate :: Num r => Complex r -> Complex r conjugate (r :+ i) = r :+ negate i +-- | Complex magnitude squared. +magnitude2 :: Num r => Complex r -> r +magnitude2 (r :+ i) = r * r + i * i + +-- | Complex magnitude. +magnitude :: Floating r => Complex r -> r +magnitude = sqrt . magnitude2 + -- | Complex phase. phase :: (Ord r, Floating r) => Complex r -> r phase (r :+ i) @@ -60,10 +85,6 @@ | 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 diff --git a/ruff.cabal b/ruff.cabal --- a/ruff.cabal +++ b/ruff.cabal @@ -1,10 +1,11 @@ Name: ruff -Version: 0.1 +Version: 0.2 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. + A library for analysis and exploration of fractals, providing + angled internal addresses, external ray tracing, nucleus and + bond point finding, and iterations for images of the Mandelbrot + Set. Homepage: https://gitorious.org/ruff License: BSD3 @@ -18,7 +19,16 @@ 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 + Exposed-modules: Fractal.RUFF.Mandelbrot.Address + Fractal.RUFF.Mandelbrot.Image + Fractal.RUFF.Mandelbrot.Iterate + Fractal.RUFF.Mandelbrot.Nucleus + 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 + GHC-Options: -Wall -O2 + GHC-Prof-Options: -prof -auto-all -caf-all