```{-# LANGUAGE DeriveDataTypeable #-}
{- |
Copyright   :  (c) Claude Heiland-Allen 2010-2011

Maintainer  :  claudiusmaximus@goto10.org
Stability   :  unstable
Portability :  portable

External angles give rise to kneading sequences under the angle doubling
when extended to angled internal addresses they distinguish hyperbolic
components in a concise and meaningful way.

The algorithms are mostly based on Dierk Schleicher's paper
/Internal Addresses Of The Mandelbrot Set And Galois Groups Of Polynomials (version of February 5, 2008)/
<http://arxiv.org/abs/math/9411238v2>.
-}

( Angle, double, wrap, prettyAngle, prettyAngles
) where

import Data.Data (Data())
import Data.Typeable (Typeable())
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)

-- | 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
| f < 0 = 1 + f
| otherwise = f
where
(_, f) = properFraction a :: (Integer, Angle)

-- | Angle doubling map.
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.
= Zero
| One
| Star
deriving (Read, Show, Eq, Ord, Enum, Bounded, Data, Typeable)

-- | Knead character representation.
kneadChar Zero = '0'
kneadChar One  = '1'
kneadChar Star = '*'

-- | Kneading sequences.  Note that the 'Aperiodic' case has an infinite list.
deriving (Read, Show, Eq, Ord, Data, Typeable)

-- | Kneading sequence as a string.  The 'Aperiodic' case is truncated arbitrarily.
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.
| a0 == 0 = StarPeriodic [Star]
| otherwise = fst kneads
where
a0 = wrap a0'
lo =  a0      / 2
hi = (a0 + 1) / 2
ks = (a0, One) : snd kneads
| 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
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)

-- | The period of a kneading sequence, or 'Nothing' when it isn't periodic.
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 a kneading sequence to an infinite list.
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

-- | Internal addresses are a non-empty sequence of strictly increasing
--   integers beginning with '1'.
deriving (Read, Show, Eq, Ord, Data, Typeable)

-- | Internal address as a string.

-- | Construct a valid 'InternalAddress', checking the precondition.
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

-- | Extract the sequence of integers.
internalToList :: InternalAddress -> [Integer]
internalToList (InternalAddress xs) = xs

-- | Construct an 'InternalAddress' from a kneading sequence.
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'per :: Integer -> [Knead] -> [Integer]
address'per p v = takeWhile (<= p) \$ address' v

where
address'' sk vk = sk : address'' sk' vk'
where
sk' = (1 +) . genericLength . takeWhile id . zipWith (==) v . cycle \$ vk
vk' = genericTake sk' (cycle v)

-- | A star-periodic kneading sequence's upper and lower associated
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

-- | The upper associated kneading sequence.
upper = fmap fst . associated

-- | The lower associated kneading sequence.
lower = fmap snd . associated

-- | Angled internal addresses have angles between each integer in an
= Unangled Integer
| Angled Integer Angle AngledInternalAddress
deriving (Read, Show, Eq, Ord, Data, Typeable)

-- | Angled internal address as a 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',
--   and the 'Integer' must be strictly increasing.
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"

-- | Convert an 'AngledInternalAddress' to a list.
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

-- | The angled internal address corresponding to an external angle.
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

-- | Split an angled internal address at the last island.
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 (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 (Unangled p) = p
addressPeriod (Angled _ _ a) = addressPeriod a

-- | Discard angle information from an internal address.
stripAngles = InternalAddress . map fst . angledToList

-- | The pair of external angles whose rays land at the root of the
--   hyperbolic component described by the angled internal address.
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
-}
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
| 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 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 s = case runP pKnead () "" s of
Left _ -> Nothing
Right k -> Just k

-- | Parse a non-aperiodic kneading sequence.
parseKneading s = case runP pKneading () "" s of
Left _ -> Nothing
Right ks -> Just ks

-- | Parse an internal address.
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.
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 " ")