module Fractal.RUFF.Mandelbrot.Address
( Angle, tune, prettyAngle, prettyAngles, angles
, BinAngle, binary, unbinary, btune, prettyBinAngle, parseBinAngle, bperiod, bpreperiod, bangles
, 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 Prelude hiding (Rational)
import Safe (headNote)
import Data.Data (Data())
import Data.Typeable (Typeable())
import Control.Monad (guard)
import Control.Monad.Identity (Identity())
import Data.Char (digitToInt)
import Data.Bits (shiftL, shiftR, (.&.), (.|.))
import Data.List (elemIndex, foldl', sort, sortBy, nub)
import Data.Maybe (mapMaybe)
import Data.Ord (comparing)
import Data.Strict.Tuple (Pair((:!:)))
import Fractal.RUFF.Types.Ratio (Q(..), Rational)
import Text.Parsec (ParsecT(), choice, digit, eof, many, many1, runParser, sepEndBy, string, try)
ceiling' :: (Q r, Z r ~ Integer) => r -> Int -> Integer
ceiling' x y = ((numerator x `shiftL` y) numerator x + denominator x 1) `div` denominator x
floor' :: (Q r, Z r ~ Integer) => r -> Int -> Integer
floor' x y = ((numerator x `shiftL` y) numerator x) `div` denominator x
angles :: Angle -> [Angle]
angles = map unbinary . bangles . binary
bangles :: BinAngle -> [BinAngle]
bangles = rays
where
periodic :: Int -> BinAngle -> Maybe BinAngle
periodic = preperiodic 0
preperiodic :: Int -> Int -> BinAngle -> Maybe BinAngle
preperiodic preperiod' period' (pre, per) =
let (pre', per') = splitAt preperiod' (pre ++ cycle per)
in check (pre', take period' per')
where
check t@(_, p)
| null p = Nothing
| binary (unbinary t) == t = Just t
| otherwise = Nothing
rays :: BinAngle -> [BinAngle]
rays t
| pp == 0 = case fmap (map binary . sort . (\(a,b) -> [a,b])) . (externalAngles =<<) . angledInternalAddress . unbinary $ t of
Just xs -> xs
Nothing -> []
| pp > 0 = case kneading (unbinary t) of
PrePeriodic _ kper -> case p `divMod` length kper of
(n, m)
| m /= 0 -> error $ "rays Preperiodic: " ++ show (p, length kper, n, m)
| n > 1 -> headNote "rays PrePeriodic" . dropWhile ((n /=) . length) . iterate (rays' n pp p) $ [t]
| n == 1 -> rays'' pp p t
| otherwise -> error $ "rays " ++ show pp ++ " " ++ show p ++ " " ++ show n
k -> error $ "rays " ++ show pp ++ " " ++ show p ++ " " ++ show k
| otherwise = error $ "rays " ++ show pp ++ " " ++ show p ++ " " ++ show t
where
pp = bpreperiod t
p = bperiod t
rays' :: Int -> Int -> Int -> [BinAngle] -> [BinAngle]
rays' n pp p ts
| not (null ts)
= sortBy (comparing unbinary)
. take (n `min` (length ts + 2))
. nub
. mapMaybe (preperiodic pp p)
. concat
. mapMaybe
( fmap (map binary . (\(a,b) -> [a,b]))
. (externalAngles =<<)
. (angledInternalAddress =<<)
. fmap unbinary
)
$ [ periodic m t | m <- [2 * pp + p ..], t <- ts ]
| otherwise = error "rays' null ts"
rays'' :: Int -> Int -> BinAngle -> [BinAngle]
rays'' pp p t
= sortBy (comparing unbinary)
. nub
. mapMaybe (fmap (binary . unbinary) . preperiodic pp p)
. concat
. mapMaybe
( fmap (map binary . (\(a,b) -> [a,b]))
. (externalAngles =<<)
. (angledInternalAddress =<<)
. fmap unbinary
)
$ [ periodic m t | m <- [2 * (pp + p) .. 3 * (pp + p)] ]
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
type BinAngle = ([Bool], [Bool])
prettyBinAngle :: BinAngle -> String
prettyBinAngle (pre, per) = "." ++ map b pre ++ "(" ++ map b per ++ ")"
where
b False = '0'
b True = '1'
parseBinAngle :: String -> Maybe BinAngle
parseBinAngle s =
case s of
'.':s1 -> case break ('('==) s1 of
(pre, '(':s2) -> case break (')'==) s2 of
(per, ")") -> case all (`elem`"01") (pre ++ per) && not (null per) of
True -> Just (map b pre, map b per)
_ -> Nothing
_ -> Nothing
_ -> Nothing
_ -> Nothing
where
b '0' = False
b '1' = True
b c = error $ "parseBinAngle.b " ++ [c]
bpreperiod :: BinAngle -> Int
bpreperiod (pre, _) = length pre
bperiod :: BinAngle -> Int
bperiod (_, per) = length per
unbinary :: BinAngle -> Angle
unbinary (pre, per)
| n == 0 = bits pre % (1 `shiftL` m)
| otherwise = (bits pre * ((1 `shiftL` n) 1) + bits per) % (((1 `shiftL` n) 1) `shiftL` m)
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 = binary' . wrap
where
binary' a
| a == zero = ([], [False])
| even (denominator a) =
let ~(pre, per) = binary' (double a)
in ((a >= half) : pre, per)
| otherwise = ([], (a >= half) : binary'' (doubleOdd a))
where
binary'' a'
| a' == a = []
| otherwise = (a' >= half) : binary'' (doubleOdd a')
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 == zero = StarPeriodic [Star]
| otherwise = case span (even . denominator . fst) . kneading' $ a0 of
(pre, ak1@(a1,_):aks) -> case takeWhile ((a1 /=) . fst) aks of
aks' ->
let per = map snd $ ak1 : aks'
in case (null pre, last per) of
(True, Star) -> StarPeriodic per
(True, _) -> Periodic (canonical per)
(False, _) -> PrePeriodic (map snd pre) (canonical per)
ppp -> error $ "kneading: " ++ show a0' ++ " " ++ show ppp
where
a0 = wrap a0'
(lo, hi) = preimages a0
kneading' :: Angle -> [(Angle, Knead)]
kneading' a
| even (denominator a) = (a, knead a) : kneading' (double a)
| otherwise = kneading'' a
kneading'' :: Angle -> [(Angle, Knead)]
kneading'' a = (a, knead a) : kneading'' (doubleOdd a)
knead a
| a == lo = Star
| a == hi = Star
| lo < a && a < hi = One
| hi < a || a < lo = Zero
| otherwise = error $ "knead " ++ show a ++ " " ++ show lo ++ " " ++ show hi
canonical ks = headNote "kneading canonical" ([ ks' | m <- [1..n], n `mod` m == 0, let ks' = take m ks, ks == take n (cycle ks') ])
where n = length ks
period :: Kneading -> Maybe Int
period (StarPeriodic k) = Just (length k)
period (Periodic k) = Just (length k)
period _ = Nothing
rho :: Kneading -> Int -> Int
rho v = rho'
where
rho' r
| r >= 1 && r `mod` n /= 0 = ((1 + r) +) . length . takeWhile id . zipWith (==) (unwrap v) . drop r $ (unwrap v)
| otherwise = rho' (r + 1)
Just n = period 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 [Int]
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 :: [Int] -> 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 -> [Int]
internalToList (InternalAddress xs) = xs
internalAddress :: Kneading -> Maybe InternalAddress
internalAddress (StarPeriodic [Star]) = Just (InternalAddress [1])
internalAddress (StarPeriodic v@(One:_)) = Just . InternalAddress . address'per (length v) $ v
internalAddress (Periodic v@(One:_)) = Just . InternalAddress . address'per (length v) $ v
internalAddress k@(Aperiodic (One:_)) = Just . InternalAddress . address'inf . unwrap $ k
internalAddress _ = Nothing
address'inf :: [Knead] -> [Int]
address'inf v = address' v
address'per :: Int -> [Knead] -> [Int]
address'per p v = takeWhile (<= p) $ address' v
address' :: [Knead] -> [Int]
address' v = address'' 1 [One]
where
address'' sk vk = sk : address'' sk' vk'
where
sk' = (1 +) . length . takeWhile id . zipWith (==) v . cycle $ vk
vk' = take sk' (cycle v)
associated :: Kneading -> Maybe (Kneading, Kneading)
associated (StarPeriodic k) = Just (Periodic a, Periodic abar)
where
n = length k
divisors = [ m | m <- [1 .. n], n `mod` m == 0 ]
abar = headNote "associated abar" . filter (and . zipWith (==) a' . cycle) . map (`take` 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 Int
| Angled Int Angle AngledInternalAddress
deriving (Read, Show, Eq, Ord, Data, Typeable)
prettyAngledInternalAddress :: AngledInternalAddress -> String
prettyAngledInternalAddress (Unangled n) = show n
prettyAngledInternalAddress (Angled n r a)
| r /= half = show n ++ " " ++ show (numerator r) ++ "/" ++ show (denominator r) ++ " " ++ prettyAngledInternalAddress a
| otherwise = show n ++ " " ++ prettyAngledInternalAddress a
angledFromList :: [(Int, 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 && zero < r && r < one = Angled n r `fmap` fromList' n xs
fromList' _ _ = Nothing
unsafeAngledFromList :: [(Int, Maybe Angle)] -> AngledInternalAddress
unsafeAngledFromList = fromList' 0
where
fromList' x [(n, Nothing)] | n > x = Unangled n
fromList' x ((n, Just r) : xs) | n > x && zero < r && r < one = Angled n r (fromList' n xs)
fromList' _ _ = error "Fractal.Mandelbrot.Address.unsafeAngledFromList"
angledToList :: AngledInternalAddress -> [(Int, Maybe Angle)]
angledToList (Unangled n) = [(n, Nothing)]
angledToList (Angled n r a) = (n, Just r) : angledToList a
denominators :: InternalAddress -> Kneading -> [Int]
denominators a v = denominators' (internalToList a)
where
denominators' (s0:ss@(s1:_)) =
let rr = r s0 s1
in (((s1 rr) `div` s0) + if (s0 ==) . headNote "denominators" . dropWhile (< 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 -> [Int] -> [Int]
numerators r a qs = zipWith num (internalToList a) qs
where
num s q = length . filter (<= r) . map (rs !!) $ [0 .. q 2]
where
rs = iterate (\t -> foldr (.) id (replicate s (if even (denominator t) then double else doubleOdd)) $ t) (wrap 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 (\a b -> fromIntegral a % fromIntegral b) 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 == fromIntegral (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 * fromIntegral (denominator r)) rs)
joinAddress (Angled p r a) rs = Angled p r (joinAddress a rs)
addressPeriod :: AngledInternalAddress -> Int
addressPeriod (Unangled p) = p
addressPeriod (Angled _ _ a) = addressPeriod a
stripAngles :: AngledInternalAddress -> InternalAddress
stripAngles = InternalAddress . map fst . angledToList
externalAngles :: AngledInternalAddress -> Maybe (Angle, Angle)
externalAngles = externalAngles' 1 (zero, one)
externalAngles' :: Int -> (Angle, Angle) -> AngledInternalAddress -> Maybe (Angle, Angle)
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 (zero, one) (fromIntegral den)
nums = [ num' | num' <- [ 1.. den 1 ], let r' = num' % den :: Angle, denominator r' == den ]
nws, nnums :: Int
nws = length ws
nnums = length nums
guard (nws == nnums)
i <- elemIndex num nums
(l,h) <- safeIndex ws i
externalAngles' (p * fromIntegral den) (if p > 1 then (tune l lohi, tune h lohi) else (l, h)) a
wakees :: (Angle, Angle) -> Int -> [(Angle, Angle)]
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 = (1 `shiftL` n) 1
in [ r
| (l :!: h) <- gs
, num <- [ ceiling' l n .. floor' h n ]
, fullperiod n num
, let r = num % den
, l < r, r < h
]
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
fullperiod bs = \n -> and [ (((n `shiftR` b) .|. (n `shiftL` (bs b))) .&. mask) /= n | b <- factors ]
where
factors = [ b | b <- [ bs 1, bs 2 .. 1 ], bs `mod` b == 0 ]
mask = (1 `shiftL` bs) 1
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"
safeIndex :: [a] -> Int -> Maybe a
safeIndex [] _ = Nothing
safeIndex (x:xs) i
| i < 0 = Nothing
| i > 0 = safeIndex xs (i 1)
| otherwise = Just x
parseAngle :: String -> Maybe Angle
parseAngle s = case runParser pFraction () "" s of
Left _ -> Nothing
Right f -> Just (unFraction f)
parseAngles :: String -> Maybe [Angle]
parseAngles s = case runParser (many pFraction) () "" s of
Left _ -> Nothing
Right fs -> Just (map unFraction fs)
parseKnead :: String -> Maybe Knead
parseKnead s = case runParser pKnead () "" s of
Left _ -> Nothing
Right k -> Just k
parseKneading :: String -> Maybe Kneading
parseKneading s = case runParser pKneading () "" s of
Left _ -> Nothing
Right ks -> Just ks
parseInternalAddress :: String -> Maybe InternalAddress
parseInternalAddress s = case runParser (many pNumber) () "" s of
Left _ -> Nothing
Right ns -> internalFromList (map (fromIntegral . unNumber) ns)
parseAngledInternalAddress :: String -> Maybe AngledInternalAddress
parseAngledInternalAddress s = case runParser 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 (fromIntegral p)
accum _ [Number n] = return $ Unangled (fromIntegral n)
accum _ (Number n : ts@(Number _ : _)) = do
a <- accum n ts
return $ Angled (fromIntegral n) (1%2) a
accum _ (Number n : Fraction t b : ts) = do
a <- accum (n * b) ts
return $ Angled (fromIntegral n) (t%b) a
accum p (Fraction t b : ts) = do
a <- accum (p * b) ts
return $ Angled (fromIntegral 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