ruff 0.1 → 0.2
raw patch · 7 files changed
+435/−80 lines, 7 filesdep −wl-pprint-text
Dependencies removed: wl-pprint-text
Files
- Fractal/RUFF/Mandelbrot/Address.hs +182/−32
- Fractal/RUFF/Mandelbrot/Image.hs +38/−16
- Fractal/RUFF/Mandelbrot/Iterate.hs +2/−4
- Fractal/RUFF/Mandelbrot/Nucleus.hs +88/−0
- Fractal/RUFF/Mandelbrot/Ray.hs +80/−14
- Fractal/RUFF/Types/Complex.hs +27/−6
- ruff.cabal +18/−8
Fractal/RUFF/Mandelbrot/Address.hs view
@@ -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
Fractal/RUFF/Mandelbrot/Image.hs view
@@ -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 -}
Fractal/RUFF/Mandelbrot/Iterate.hs view
@@ -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
+ Fractal/RUFF/Mandelbrot/Nucleus.hs view
@@ -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/+-- <http://mrob.com/pub/muency/newtonraphsonmethod.html>.+--+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+-- <http://mrob.com/pub/muency/derivative.html>.+--+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/+-- <http://mrob.com/pub/muency/period.html>.+--+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
Fractal/RUFF/Mandelbrot/Ray.hs view
@@ -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/-<http://www.math.nagoya-u.ac.jp/~kawahira/programs/mandel-exray.pdf>.+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/+-- <http://www.math.nagoya-u.ac.jp/~kawahira/programs/mandel-exray.pdf>.+--+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 []
Fractal/RUFF/Types/Complex.hs view
@@ -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
ruff.cabal view
@@ -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