Address.hs

{-# LANGUAGE DeriveDataTypeable #-}
{- |
Module      :  Fractal.RUFF.Mandelbrot.Address
Copyright   :  (c) Claude Heiland-Allen 2010-2011
License     :  BSD3

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

External angles give rise to kneading sequences under the angle doubling
map.  Internal addresses encode kneading sequences in human-readable form,
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>.
-}

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

-- | 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.
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 = '*'

-- | Kneading sequences.  Note that the 'Aperiodic' case has an infinite list,
--   which the 'Pretty' instance truncates arbitrarily.
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)

-- | The kneading sequence for an external angle.
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)

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

-- | Internal addresses are a non-empty sequence of strictly increasing
--   integers beginning with '1'.
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)

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

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

-- | The upper associated kneading sequence.
upper :: Kneading -> Maybe Kneading
upper = fmap fst . associated

-- | The lower associated kneading sequence.
lower :: Kneading -> Maybe Kneading
lower = fmap fst . associated

-- | Angled internal addresses have angles between each integer in an
--   internal address.
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

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

-- | Discard angle information from an internal address.
stripAngles :: AngledInternalAddress -> InternalAddress
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 :: 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 " ")