module Fractal.RUFF.Mandelbrot.Address
( 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())
import Data.Typeable (Typeable())
import Control.Monad (guard)
import Control.Monad.Identity (Identity())
import Data.Char (digitToInt)
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)
type Angle = Rational
prettyAngle :: Angle -> String
prettyAngle a = show (numerator a) ++ " / " ++ show (denominator a)
prettyAngles :: [Angle] -> String
prettyAngles [] = ""
prettyAngles [a] = show (numerator a) ++ "/" ++ show (denominator a)
prettyAngles (a:as) = show (numerator a) ++ "/" ++ show (denominator a) ++ " " ++ prettyAngles as
wrap :: Angle -> Angle
wrap a
| f < 0 = 1 + f
| otherwise = f
where
(_, f) = properFraction a :: (Integer, Angle)
double :: Angle -> Angle
double a = wrap (2 * a)
type BinAngle = ([Bool], [Bool])
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
bits :: [Bool] -> Integer
bits = foldl' (\ a b -> 2 * a + if b then 1 else 0) 0
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)
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"
tune :: Angle -> (Angle, Angle) -> Angle
tune t (t0, t1) = unbinary $ btune (binary t) (binary t0, binary t1)
data Knead
= Zero
| One
| Star
deriving (Read, Show, Eq, Ord, Enum, Bounded, Data, Typeable)
kneadChar :: Knead -> Char
kneadChar Zero = '0'
kneadChar One = '1'
kneadChar Star = '*'
data Kneading
= Aperiodic [Knead]
| PrePeriodic [Knead] [Knead]
| StarPeriodic [Knead]
| Periodic [Knead]
deriving (Read, Show, Eq, Ord, Data, Typeable)
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 ++ ")"
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)
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 :: 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
data InternalAddress = InternalAddress [Integer]
deriving (Read, Show, Eq, Ord, Data, Typeable)
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)
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
internalToList :: InternalAddress -> [Integer]
internalToList (InternalAddress xs) = xs
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)
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
upper :: Kneading -> Maybe Kneading
upper = fmap fst . associated
lower :: Kneading -> Maybe Kneading
lower = fmap snd . associated
data AngledInternalAddress
= Unangled Integer
| Angled Integer Angle AngledInternalAddress
deriving (Read, Show, Eq, Ord, Data, Typeable)
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
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"
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
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
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 []
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)
addressPeriod :: AngledInternalAddress -> Integer
addressPeriod (Unangled p) = p
addressPeriod (Angled _ _ a) = addressPeriod a
stripAngles :: AngledInternalAddress -> InternalAddress
stripAngles = InternalAddress . map fst . angledToList
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
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)]
| 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
parseAngle :: String -> Maybe Angle
parseAngle s = case runP pFraction () "" s of
Left _ -> Nothing
Right f -> Just (unFraction f)
parseAngles :: String -> Maybe [Angle]
parseAngles s = case runP (many pFraction) () "" s of
Left _ -> Nothing
Right fs -> Just (map unFraction fs)
parseKnead :: String -> Maybe Knead
parseKnead s = case runP pKnead () "" s of
Left _ -> Nothing
Right k -> Just k
parseKneading :: String -> Maybe Kneading
parseKneading s = case runP pKneading () "" s of
Left _ -> Nothing
Right ks -> Just ks
parseInternalAddress :: String -> Maybe InternalAddress
parseInternalAddress s = case runP (many pNumber) () "" s of
Left _ -> Nothing
Right ns -> internalFromList (map unNumber ns)
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
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 " ")
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