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, (<>))
type Angle = Rational
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)
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 = '*'
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)
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)
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)
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 fst . associated
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
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
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
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 :: 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 " ")