```{-# 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
, stripAngles
, parse
) where

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

-- | Angle as a fraction of a turn, usually in [0, 1).
type Angle = Rational

-- | 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)

-- | Elements of kneading sequences.
= Zero
| One
| Star
deriving (Read, Show, Eq, Ord, Enum, Bounded, Data, Typeable)

prettyList = pretty . map kneadChar

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

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)

-- | The kneading sequence for an external angle.
| a0 == 0 = StarPeriodic [Star]
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, [])
| 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)

pretty (InternalAddress [x]) = pretty x
pretty (InternalAddress (x:ys)) = pretty x <> char ' ' <> pretty (InternalAddress ys)

-- | 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.

where
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 fst . associated

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

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

-- | 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 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.
let r = wrap r0
let d = denominators i k
n = numerators r i d
return . unsafeAngledFromList . zip (internalToList i) . (++ [Nothing]) . map Just . zipWith (%) n \$ d

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

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 an angled internal address, accepting some unambiguous
--   abbreviations.
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 = 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 " ")
```