{-# LANGUAGE FlexibleContexts      #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE Rank2Types            #-}
{-# LANGUAGE ViewPatterns          #-}
{-# OPTIONS_GHC -fno-warn-missing-methods #-}

-----------------------------------------------------------------------------
-- |
-- Module      :  Diagrams.TwoD.Apollonian
-- Copyright   :  (c) 2011, 2016 Brent Yorgey
-- License     :  BSD-style (see LICENSE)
-- Maintainer  :  byorgey@cis.upenn.edu
--
-- Generation of Apollonian gaskets.  Any three mutually tangent
-- circles uniquely determine exactly two others which are mutually
-- tangent to all three.  This process can be repeated, generating a
-- fractal circle packing.
--
-- See J. Lagarias, C. Mallows, and A. Wilks, \"Beyond the Descartes
-- circle theorem\", /Amer. Math. Monthly/ 109 (2002), 338--361.
-- <http://arxiv.org/abs/math/0101066>.
--
-- A few examples:
--
-- > import Diagrams.TwoD.Apollonian
-- > apollonian1 = apollonianGasket 0.01 2 2 2
--
-- <<diagrams/src_Diagrams_TwoD_Apollonian_apollonian1.svg#diagram=apollonian1&width=400>>
--
-- > import Diagrams.TwoD.Apollonian
-- > apollonian2 = apollonianGasket 0.01 2 3 3
--
-- <<diagrams/src_Diagrams_TwoD_Apollonian_apollonian2.svg#diagram=apollonian2&width=400>>
--
-- > import Diagrams.TwoD.Apollonian
-- > apollonian3 = apollonianGasket 0.01 2 4 7
--
-- <<diagrams/src_Diagrams_TwoD_Apollonian_apollonian3.svg#diagram=apollonian3&width=400>>
--
-----------------------------------------------------------------------------

module Diagrams.TwoD.Apollonian
       ( -- * Circles

         Circle(..), mkCircle, center, radius

         -- * Descartes' Theorem
       , descartes, other, initialConfig

         -- * Apollonian gasket generation

       , apollonian

         -- ** Kissing sets

       , KissingSet(..), kissingSets, flipSelected, selectOthers

         -- ** Apollonian trees

       , apollonianTrees, apollonianTree

         -- * Diagram generation

       , drawCircle
       , drawGasket
       , apollonianGasket

       ) where

import           Data.Complex
import qualified Data.Foldable    as F
import           Data.Maybe       (catMaybes)
import           Data.Tree

import           Diagrams.Prelude hiding (center, radius)

import           Control.Arrow    (second, (&&&))

------------------------------------------------------------
--  Circles
------------------------------------------------------------

-- | Representation for circles that lets us quickly compute an
--   Apollonian gasket.
data Circle n = Circle
  { forall n. Circle n -> n
bend :: n
    -- ^ The bend is the reciprocal of signed
    --   radius: a negative radius means the
    --   outside and inside of the circle are
    --   switched.  The bends of any four mutually
    --   tangent circles satisfy Descartes'
    --   Theorem.
  , forall n. Circle n -> Complex n
cb   :: Complex n
    -- ^ /Product/ of bend and center represented
    --   as a complex number.  Amazingly, these
    --   products also satisfy the equation of
    --   Descartes' Theorem.
  }
  deriving (Circle n -> Circle n -> Bool
(Circle n -> Circle n -> Bool)
-> (Circle n -> Circle n -> Bool) -> Eq (Circle n)
forall n. Eq n => Circle n -> Circle n -> Bool
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: forall n. Eq n => Circle n -> Circle n -> Bool
== :: Circle n -> Circle n -> Bool
$c/= :: forall n. Eq n => Circle n -> Circle n -> Bool
/= :: Circle n -> Circle n -> Bool
Eq, Int -> Circle n -> ShowS
[Circle n] -> ShowS
Circle n -> String
(Int -> Circle n -> ShowS)
-> (Circle n -> String) -> ([Circle n] -> ShowS) -> Show (Circle n)
forall n. Show n => Int -> Circle n -> ShowS
forall n. Show n => [Circle n] -> ShowS
forall n. Show n => Circle n -> String
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: forall n. Show n => Int -> Circle n -> ShowS
showsPrec :: Int -> Circle n -> ShowS
$cshow :: forall n. Show n => Circle n -> String
show :: Circle n -> String
$cshowList :: forall n. Show n => [Circle n] -> ShowS
showList :: [Circle n] -> ShowS
Show)

-- | Create a @Circle@ given a signed radius and a location for its center.
mkCircle :: Fractional n =>
            n -- ^ signed radius
         -> P2 n     -- ^ center
         -> Circle n
mkCircle :: forall n. Fractional n => n -> P2 n -> Circle n
mkCircle n
r (P2 n -> (n, n)
forall n. P2 n -> (n, n)
unp2 -> (n
x,n
y)) = n -> Complex n -> Circle n
forall n. n -> Complex n -> Circle n
Circle (n
1n -> n -> n
forall a. Fractional a => a -> a -> a
/n
r) (n
bn -> n -> n
forall a. Num a => a -> a -> a
*n
x n -> n -> Complex n
forall a. a -> a -> Complex a
:+ n
bn -> n -> n
forall a. Num a => a -> a -> a
*n
y)
  where b :: n
b = n
1n -> n -> n
forall a. Fractional a => a -> a -> a
/n
r

-- | Get the center of a circle.
center :: Fractional n => Circle n -> P2 n
center :: forall n. Fractional n => Circle n -> P2 n
center (Circle n
b (n
cbx :+ n
cby)) = (n, n) -> P2 n
forall n. (n, n) -> P2 n
p2 (n
cbx n -> n -> n
forall a. Fractional a => a -> a -> a
/ n
b, n
cby n -> n -> n
forall a. Fractional a => a -> a -> a
/ n
b)

-- | Get the (unsigned) radius of a circle.
radius :: Fractional n => Circle n -> n
radius :: forall n. Fractional n => Circle n -> n
radius = n -> n
forall a. Num a => a -> a
abs (n -> n) -> (Circle n -> n) -> Circle n -> n
forall b c a. (b -> c) -> (a -> b) -> a -> c
. n -> n
forall a. Fractional a => a -> a
recip (n -> n) -> (Circle n -> n) -> Circle n -> n
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Circle n -> n
forall n. Circle n -> n
bend

liftF :: RealFloat n => (forall a. Floating a => a -> a) -> Circle n -> Circle n
liftF :: forall n.
RealFloat n =>
(forall a. Floating a => a -> a) -> Circle n -> Circle n
liftF forall a. Floating a => a -> a
f (Circle n
b Complex n
c) = n -> Complex n -> Circle n
forall n. n -> Complex n -> Circle n
Circle (n -> n
forall a. Floating a => a -> a
f n
b) (Complex n -> Complex n
forall a. Floating a => a -> a
f Complex n
c)

liftF2 :: RealFloat n => (forall a. Floating a => a -> a -> a) ->
          Circle n -> Circle n -> Circle n
liftF2 :: forall n.
RealFloat n =>
(forall a. Floating a => a -> a -> a)
-> Circle n -> Circle n -> Circle n
liftF2 forall a. Floating a => a -> a -> a
f (Circle n
b1 Complex n
cb1) (Circle n
b2 Complex n
cb2) = n -> Complex n -> Circle n
forall n. n -> Complex n -> Circle n
Circle (n -> n -> n
forall a. Floating a => a -> a -> a
f n
b1 n
b2) (Complex n -> Complex n -> Complex n
forall a. Floating a => a -> a -> a
f Complex n
cb1 Complex n
cb2)

instance RealFloat n => Num (Circle n) where
  + :: Circle n -> Circle n -> Circle n
(+)           = (forall a. Floating a => a -> a -> a)
-> Circle n -> Circle n -> Circle n
forall n.
RealFloat n =>
(forall a. Floating a => a -> a -> a)
-> Circle n -> Circle n -> Circle n
liftF2 a -> a -> a
forall a. Floating a => a -> a -> a
forall a. Num a => a -> a -> a
(+)
  (-)           = (forall a. Floating a => a -> a -> a)
-> Circle n -> Circle n -> Circle n
forall n.
RealFloat n =>
(forall a. Floating a => a -> a -> a)
-> Circle n -> Circle n -> Circle n
liftF2 (-)
  * :: Circle n -> Circle n -> Circle n
(*)           = (forall a. Floating a => a -> a -> a)
-> Circle n -> Circle n -> Circle n
forall n.
RealFloat n =>
(forall a. Floating a => a -> a -> a)
-> Circle n -> Circle n -> Circle n
liftF2 a -> a -> a
forall a. Floating a => a -> a -> a
forall a. Num a => a -> a -> a
(*)
  negate :: Circle n -> Circle n
negate        = (forall a. Floating a => a -> a) -> Circle n -> Circle n
forall n.
RealFloat n =>
(forall a. Floating a => a -> a) -> Circle n -> Circle n
liftF a -> a
forall a. Floating a => a -> a
forall a. Num a => a -> a
negate
  abs :: Circle n -> Circle n
abs           = (forall a. Floating a => a -> a) -> Circle n -> Circle n
forall n.
RealFloat n =>
(forall a. Floating a => a -> a) -> Circle n -> Circle n
liftF a -> a
forall a. Floating a => a -> a
forall a. Num a => a -> a
abs
  fromInteger :: Integer -> Circle n
fromInteger Integer
n = n -> Complex n -> Circle n
forall n. n -> Complex n -> Circle n
Circle (Integer -> n
forall a. Num a => Integer -> a
fromInteger Integer
n) (Integer -> Complex n
forall a. Num a => Integer -> a
fromInteger Integer
n)

instance RealFloat n => Fractional (Circle n) where
  / :: Circle n -> Circle n -> Circle n
(/)   = (forall a. Floating a => a -> a -> a)
-> Circle n -> Circle n -> Circle n
forall n.
RealFloat n =>
(forall a. Floating a => a -> a -> a)
-> Circle n -> Circle n -> Circle n
liftF2 a -> a -> a
forall a. Floating a => a -> a -> a
forall a. Fractional a => a -> a -> a
(/)
  recip :: Circle n -> Circle n
recip = (forall a. Floating a => a -> a) -> Circle n -> Circle n
forall n.
RealFloat n =>
(forall a. Floating a => a -> a) -> Circle n -> Circle n
liftF a -> a
forall a. Floating a => a -> a
forall a. Fractional a => a -> a
recip

-- | The @Num@, @Fractional@, and @Floating@ instances for @Circle@
--   (all simply lifted elementwise over @Circle@'s fields) let us use
--   Descartes' Theorem directly on circles.
instance RealFloat n => Floating (Circle n) where
  sqrt :: Circle n -> Circle n
sqrt = (forall a. Floating a => a -> a) -> Circle n -> Circle n
forall n.
RealFloat n =>
(forall a. Floating a => a -> a) -> Circle n -> Circle n
liftF a -> a
forall a. Floating a => a -> a
sqrt

------------------------------------------------------------
--  Descartes' Theorem
------------------------------------------------------------

-- XXX generalize these for higher dimensions?

-- | Descartes' Theorem states that if @b1@, @b2@, @b3@ and @b4@ are
--   the bends of four mutually tangent circles, then
--
--   @
--     b1^2 + b2^2 + b3^2 + b4^2 = 1/2 * (b1 + b2 + b3 + b4)^2.
--   @
--
--   Surprisingly, if we replace each of the @bi@ with the /product/
--   of @bi@ and the center of the corresponding circle (represented
--   as a complex number), the equation continues to hold! (See the
--   paper referenced at the top of the module.)
--
--   @descartes [b1,b2,b3]@ solves for @b4@, returning both solutions.
--   Notably, @descartes@ works for any instance of @Floating@, which
--   includes both @Double@ (for bends), @Complex Double@ (for
--   bend/center product), and @Circle@ (for both at once).
descartes :: Floating n => [n] -> [n]
descartes :: forall n. Floating n => [n] -> [n]
descartes [n
b1,n
b2,n
b3] = [n
r n -> n -> n
forall a. Num a => a -> a -> a
+ n
s, -n
r n -> n -> n
forall a. Num a => a -> a -> a
+ n
s]
  where r :: n
r = n
2 n -> n -> n
forall a. Num a => a -> a -> a
* n -> n
forall a. Floating a => a -> a
sqrt (n
b1n -> n -> n
forall a. Num a => a -> a -> a
*n
b2 n -> n -> n
forall a. Num a => a -> a -> a
+ n
b1n -> n -> n
forall a. Num a => a -> a -> a
*n
b3 n -> n -> n
forall a. Num a => a -> a -> a
+ n
b2n -> n -> n
forall a. Num a => a -> a -> a
*n
b3)
        s :: n
s = n
b1n -> n -> n
forall a. Num a => a -> a -> a
+n
b2n -> n -> n
forall a. Num a => a -> a -> a
+n
b3
descartes [n]
_ = String -> [n]
forall a. HasCallStack => String -> a
error String
"descartes must be called on a list of length 3"

-- | If we have /four/ mutually tangent circles we can choose one of
--   them to replace; the remaining three determine exactly one other
--   circle which is mutually tangent.  However, in this situation
--   there is no need to apply 'descartes' again, since the two
--   solutions @b4@ and @b4'@ satisfy
--
--   @
--     b4 + b4' = 2 * (b1 + b2 + b3)
--   @
--
--   Hence, to replace @b4@ with its dual, we need only sum the other
--   three, multiply by two, and subtract @b4@.  Again, this works for
--   bends as well as bend/center products.
other :: Num n => [n] -> n -> n
other :: forall n. Num n => [n] -> n -> n
other [n]
xs n
x = n
2 n -> n -> n
forall a. Num a => a -> a -> a
* [n] -> n
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [n]
xs n -> n -> n
forall a. Num a => a -> a -> a
- n
x

-- | Generate an initial configuration of four mutually tangent
--   circles, given just the signed bends of three of them.
initialConfig :: RealFloat n => n -> n -> n -> [Circle n]
initialConfig :: forall n. RealFloat n => n -> n -> n -> [Circle n]
initialConfig n
b1 n
b2 n
b3 = [Circle n]
cs [Circle n] -> [Circle n] -> [Circle n]
forall a. [a] -> [a] -> [a]
++ [Circle n
c4]
  where cs :: [Circle n]
cs     = [n -> Complex n -> Circle n
forall n. n -> Complex n -> Circle n
Circle n
b1 Complex n
0, n -> Complex n -> Circle n
forall n. n -> Complex n -> Circle n
Circle n
b2 ((n
b2n -> n -> n
forall a. Fractional a => a -> a -> a
/n
b1 n -> n -> n
forall a. Num a => a -> a -> a
+ n
1) n -> n -> Complex n
forall a. a -> a -> Complex a
:+ n
0), n -> Complex n -> Circle n
forall n. n -> Complex n -> Circle n
Circle n
b3 Complex n
cb3]
        a :: n
a      = n
1n -> n -> n
forall a. Fractional a => a -> a -> a
/n
b1 n -> n -> n
forall a. Num a => a -> a -> a
+ n
1n -> n -> n
forall a. Fractional a => a -> a -> a
/n
b2
        b :: n
b      = n
1n -> n -> n
forall a. Fractional a => a -> a -> a
/n
b1 n -> n -> n
forall a. Num a => a -> a -> a
+ n
1n -> n -> n
forall a. Fractional a => a -> a -> a
/n
b3
        c :: n
c      = n
1n -> n -> n
forall a. Fractional a => a -> a -> a
/n
b2 n -> n -> n
forall a. Num a => a -> a -> a
+ n
1n -> n -> n
forall a. Fractional a => a -> a -> a
/n
b3
        x :: n
x      = (n
bn -> n -> n
forall a. Num a => a -> a -> a
*n
b n -> n -> n
forall a. Num a => a -> a -> a
+ n
an -> n -> n
forall a. Num a => a -> a -> a
*n
a n -> n -> n
forall a. Num a => a -> a -> a
- n
cn -> n -> n
forall a. Num a => a -> a -> a
*n
c)n -> n -> n
forall a. Fractional a => a -> a -> a
/(n
2n -> n -> n
forall a. Num a => a -> a -> a
*n
a)
        y :: n
y      = n -> n
forall a. Floating a => a -> a
sqrt (n
bn -> n -> n
forall a. Num a => a -> a -> a
*n
b n -> n -> n
forall a. Num a => a -> a -> a
- n
xn -> n -> n
forall a. Num a => a -> a -> a
*n
x)
        cb3 :: Complex n
cb3    = n
b3n -> n -> n
forall a. Num a => a -> a -> a
*n
x n -> n -> Complex n
forall a. a -> a -> Complex a
:+ n
b3n -> n -> n
forall a. Num a => a -> a -> a
*n
y
        [Circle n
c4,Circle n
_] = [Circle n] -> [Circle n]
forall n. Floating n => [n] -> [n]
descartes [Circle n]
cs

------------------------------------------------------------
--  Gasket generation
------------------------------------------------------------

select :: [a] -> [(a, [a])]
select :: forall a. [a] -> [(a, [a])]
select [] = []
select (a
x:[a]
xs) = (a
x,[a]
xs) (a, [a]) -> [(a, [a])] -> [(a, [a])]
forall a. a -> [a] -> [a]
: (((a, [a]) -> (a, [a])) -> [(a, [a])] -> [(a, [a])]
forall a b. (a -> b) -> [a] -> [b]
map (((a, [a]) -> (a, [a])) -> [(a, [a])] -> [(a, [a])])
-> (([a] -> [a]) -> (a, [a]) -> (a, [a]))
-> ([a] -> [a])
-> [(a, [a])]
-> [(a, [a])]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. ([a] -> [a]) -> (a, [a]) -> (a, [a])
forall b c d. (b -> c) -> (d, b) -> (d, c)
forall (a :: * -> * -> *) b c d.
Arrow a =>
a b c -> a (d, b) (d, c)
second) (a
xa -> [a] -> [a]
forall a. a -> [a] -> [a]
:) ([a] -> [(a, [a])]
forall a. [a] -> [(a, [a])]
select [a]
xs)

-- | The basic idea of a kissing set is supposed to represent a set of
--   four mutually tangent circles with one selected, though in fact
--   it is more general than that: it represents any set of objects
--   with one distinguished object selected.
data KissingSet n = KS { forall n. KissingSet n -> n
selected :: n, forall n. KissingSet n -> [n]
others :: [n] }
  deriving (Int -> KissingSet n -> ShowS
[KissingSet n] -> ShowS
KissingSet n -> String
(Int -> KissingSet n -> ShowS)
-> (KissingSet n -> String)
-> ([KissingSet n] -> ShowS)
-> Show (KissingSet n)
forall n. Show n => Int -> KissingSet n -> ShowS
forall n. Show n => [KissingSet n] -> ShowS
forall n. Show n => KissingSet n -> String
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: forall n. Show n => Int -> KissingSet n -> ShowS
showsPrec :: Int -> KissingSet n -> ShowS
$cshow :: forall n. Show n => KissingSet n -> String
show :: KissingSet n -> String
$cshowList :: forall n. Show n => [KissingSet n] -> ShowS
showList :: [KissingSet n] -> ShowS
Show)

-- | Generate all possible kissing sets from a set of objects by
--   selecting each object in turn.
kissingSets :: [n] -> [KissingSet n]
kissingSets :: forall n. [n] -> [KissingSet n]
kissingSets = ((n, [n]) -> KissingSet n) -> [(n, [n])] -> [KissingSet n]
forall a b. (a -> b) -> [a] -> [b]
map ((n -> [n] -> KissingSet n) -> (n, [n]) -> KissingSet n
forall a b c. (a -> b -> c) -> (a, b) -> c
uncurry n -> [n] -> KissingSet n
forall n. n -> [n] -> KissingSet n
KS) ([(n, [n])] -> [KissingSet n])
-> ([n] -> [(n, [n])]) -> [n] -> [KissingSet n]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [n] -> [(n, [n])]
forall a. [a] -> [(a, [a])]
select

-- | \"Flip\" the selected circle to the 'other' circle mutually tangent
--   to the other three.  The new circle remains selected.
flipSelected :: Num n => KissingSet n -> KissingSet n
flipSelected :: forall n. Num n => KissingSet n -> KissingSet n
flipSelected (KS n
c [n]
cs) = n -> [n] -> KissingSet n
forall n. n -> [n] -> KissingSet n
KS ([n] -> n -> n
forall n. Num n => [n] -> n -> n
other [n]
cs n
c) [n]
cs

-- | Make the selected circle unselected, and select each of the
--   others, generating a new kissing set for each.
selectOthers :: KissingSet n -> [KissingSet n]
selectOthers :: forall n. KissingSet n -> [KissingSet n]
selectOthers (KS n
c [n]
cs) = [ n -> [n] -> KissingSet n
forall n. n -> [n] -> KissingSet n
KS n
c' (n
cn -> [n] -> [n]
forall a. a -> [a] -> [a]
:[n]
cs') | (n
c',[n]
cs') <- [n] -> [(n, [n])]
forall a. [a] -> [(a, [a])]
select [n]
cs ]

-- | Given a threshold radius and a list of /four/ mutually tangent
--   circles, generate the Apollonian gasket containing those circles.
--   Stop the recursion when encountering a circle with an (unsigned)
--   radius smaller than the threshold.
apollonian :: RealFloat n => n -> [Circle n] -> [Circle n]
apollonian :: forall n. RealFloat n => n -> [Circle n] -> [Circle n]
apollonian n
thresh [Circle n]
cs
  = ([Circle n]
cs[Circle n] -> [Circle n] -> [Circle n]
forall a. [a] -> [a] -> [a]
++)
  ([Circle n] -> [Circle n])
-> ([Circle n] -> [Circle n]) -> [Circle n] -> [Circle n]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [[Circle n]] -> [Circle n]
forall (t :: * -> *) a. Foldable t => t [a] -> [a]
concat
  ([[Circle n]] -> [Circle n])
-> ([Circle n] -> [[Circle n]]) -> [Circle n] -> [Circle n]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Tree (KissingSet (Circle n)) -> [Circle n])
-> [Tree (KissingSet (Circle n))] -> [[Circle n]]
forall a b. (a -> b) -> [a] -> [b]
map ([Circle n]
-> (Tree (Circle n) -> [Circle n])
-> Maybe (Tree (Circle n))
-> [Circle n]
forall b a. b -> (a -> b) -> Maybe a -> b
maybe [] Tree (Circle n) -> [Circle n]
forall a. Tree a -> [a]
flatten (Maybe (Tree (Circle n)) -> [Circle n])
-> (Tree (KissingSet (Circle n)) -> Maybe (Tree (Circle n)))
-> Tree (KissingSet (Circle n))
-> [Circle n]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Circle n -> Bool) -> Tree (Circle n) -> Maybe (Tree (Circle n))
forall a. (a -> Bool) -> Tree a -> Maybe (Tree a)
prune Circle n -> Bool
p (Tree (Circle n) -> Maybe (Tree (Circle n)))
-> (Tree (KissingSet (Circle n)) -> Tree (Circle n))
-> Tree (KissingSet (Circle n))
-> Maybe (Tree (Circle n))
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (KissingSet (Circle n) -> Circle n)
-> Tree (KissingSet (Circle n)) -> Tree (Circle n)
forall a b. (a -> b) -> Tree a -> Tree b
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap KissingSet (Circle n) -> Circle n
forall n. KissingSet n -> n
selected)
  ([Tree (KissingSet (Circle n))] -> [[Circle n]])
-> ([Circle n] -> [Tree (KissingSet (Circle n))])
-> [Circle n]
-> [[Circle n]]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [Circle n] -> [Tree (KissingSet (Circle n))]
forall n.
RealFloat n =>
[Circle n] -> [Tree (KissingSet (Circle n))]
apollonianTrees
  ([Circle n] -> [Circle n]) -> [Circle n] -> [Circle n]
forall a b. (a -> b) -> a -> b
$ [Circle n]
cs
  where
    p :: Circle n -> Bool
p Circle n
c = Circle n -> n
forall n. Fractional n => Circle n -> n
radius Circle n
c n -> n -> Bool
forall a. Ord a => a -> a -> Bool
>= n
thresh

-- | Given a set of /four/ mutually tangent circles, generate the
--   infinite Apollonian tree rooted at the given set, represented as
--   a list of four subtrees.  Each node in the tree is a kissing set
--   with one circle selected which has just been flipped.  The three
--   children of a node represent the kissing sets obtained by
--   selecting each of the other three circles and flipping them.  The
--   initial roots of the four trees are chosen by selecting and
--   flipping each of the circles in the starting set. This
--   representation has the property that each circle in the
--   Apollonian gasket is the selected circle in exactly one node
--   (except that the initial four circles never appear as the
--   selected circle in any node).
apollonianTrees :: RealFloat n => [Circle n] -> [Tree (KissingSet (Circle n))]
apollonianTrees :: forall n.
RealFloat n =>
[Circle n] -> [Tree (KissingSet (Circle n))]
apollonianTrees = (KissingSet (Circle n) -> Tree (KissingSet (Circle n)))
-> [KissingSet (Circle n)] -> [Tree (KissingSet (Circle n))]
forall a b. (a -> b) -> [a] -> [b]
map (KissingSet (Circle n) -> Tree (KissingSet (Circle n))
forall n.
RealFloat n =>
KissingSet (Circle n) -> Tree (KissingSet (Circle n))
apollonianTree (KissingSet (Circle n) -> Tree (KissingSet (Circle n)))
-> (KissingSet (Circle n) -> KissingSet (Circle n))
-> KissingSet (Circle n)
-> Tree (KissingSet (Circle n))
forall b c a. (b -> c) -> (a -> b) -> a -> c
. KissingSet (Circle n) -> KissingSet (Circle n)
forall n. Num n => KissingSet n -> KissingSet n
flipSelected) ([KissingSet (Circle n)] -> [Tree (KissingSet (Circle n))])
-> ([Circle n] -> [KissingSet (Circle n)])
-> [Circle n]
-> [Tree (KissingSet (Circle n))]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [Circle n] -> [KissingSet (Circle n)]
forall n. [n] -> [KissingSet n]
kissingSets

-- | Generate a single Apollonian tree from a root kissing set.  See
--   the documentation for 'apollonianTrees' for an explanation.
apollonianTree :: RealFloat n => KissingSet (Circle n) -> Tree (KissingSet (Circle n))
apollonianTree :: forall n.
RealFloat n =>
KissingSet (Circle n) -> Tree (KissingSet (Circle n))
apollonianTree = (KissingSet (Circle n)
 -> (KissingSet (Circle n), [KissingSet (Circle n)]))
-> KissingSet (Circle n) -> Tree (KissingSet (Circle n))
forall b a. (b -> (a, [b])) -> b -> Tree a
unfoldTree (KissingSet (Circle n) -> KissingSet (Circle n)
forall a. a -> a
id (KissingSet (Circle n) -> KissingSet (Circle n))
-> (KissingSet (Circle n) -> [KissingSet (Circle n)])
-> KissingSet (Circle n)
-> (KissingSet (Circle n), [KissingSet (Circle n)])
forall b c c'. (b -> c) -> (b -> c') -> b -> (c, c')
forall (a :: * -> * -> *) b c c'.
Arrow a =>
a b c -> a b c' -> a b (c, c')
&&& ((KissingSet (Circle n) -> KissingSet (Circle n))
-> [KissingSet (Circle n)] -> [KissingSet (Circle n)]
forall a b. (a -> b) -> [a] -> [b]
map KissingSet (Circle n) -> KissingSet (Circle n)
forall n. Num n => KissingSet n -> KissingSet n
flipSelected ([KissingSet (Circle n)] -> [KissingSet (Circle n)])
-> (KissingSet (Circle n) -> [KissingSet (Circle n)])
-> KissingSet (Circle n)
-> [KissingSet (Circle n)]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. KissingSet (Circle n) -> [KissingSet (Circle n)]
forall n. KissingSet n -> [KissingSet n]
selectOthers))

-- | Prune a tree at the shallowest points where the predicate is not
--   satisfied.
prune :: (a -> Bool) -> Tree a -> Maybe (Tree a)
prune :: forall a. (a -> Bool) -> Tree a -> Maybe (Tree a)
prune a -> Bool
p (Node a
a [Tree a]
ts)
  | Bool -> Bool
not (a -> Bool
p a
a) = Maybe (Tree a)
forall a. Maybe a
Nothing
  | Bool
otherwise = Tree a -> Maybe (Tree a)
forall a. a -> Maybe a
Just (Tree a -> Maybe (Tree a)) -> Tree a -> Maybe (Tree a)
forall a b. (a -> b) -> a -> b
$ a -> [Tree a] -> Tree a
forall a. a -> [Tree a] -> Tree a
Node a
a ([Maybe (Tree a)] -> [Tree a]
forall a. [Maybe a] -> [a]
catMaybes ((Tree a -> Maybe (Tree a)) -> [Tree a] -> [Maybe (Tree a)]
forall a b. (a -> b) -> [a] -> [b]
map ((a -> Bool) -> Tree a -> Maybe (Tree a)
forall a. (a -> Bool) -> Tree a -> Maybe (Tree a)
prune a -> Bool
p) [Tree a]
ts))

------------------------------------------------------------
--  Diagram generation
------------------------------------------------------------

-- | Draw a circle.
drawCircle :: (Renderable (Path V2 n) b, TypeableFloat n) =>
              Circle n -> QDiagram b V2 n Any
drawCircle :: forall n b.
(Renderable (Path V2 n) b, TypeableFloat n) =>
Circle n -> QDiagram b V2 n Any
drawCircle Circle n
c = n -> QDiagram b V2 n Any
forall t n.
(TrailLike t, V t ~ V2, N t ~ n, Transformable t) =>
n -> t
circle (Circle n -> n
forall n. Fractional n => Circle n -> n
radius Circle n
c) QDiagram b V2 n Any
-> (QDiagram b V2 n Any -> QDiagram b V2 n Any)
-> QDiagram b V2 n Any
forall a b. a -> (a -> b) -> b
# Point V2 n -> QDiagram b V2 n Any -> QDiagram b V2 n Any
forall (v :: * -> *) n t.
(InSpace v n t, HasOrigin t) =>
Point v n -> t -> t
moveTo (Circle n -> Point V2 n
forall n. Fractional n => Circle n -> P2 n
center Circle n
c)
                                 # fcA transparent

-- | Draw a generated gasket, using a line width 0.003 times the
--   radius of the largest circle.
drawGasket :: (Renderable (Path V2 n) b, TypeableFloat n) =>
              [Circle n] -> QDiagram b V2 n Any
drawGasket :: forall n b.
(Renderable (Path V2 n) b, TypeableFloat n) =>
[Circle n] -> QDiagram b V2 n Any
drawGasket [Circle n]
cs = (Circle n -> QDiagram b V2 n Any)
-> [Circle n] -> QDiagram b V2 n Any
forall m a. Monoid m => (a -> m) -> [a] -> m
forall (t :: * -> *) m a.
(Foldable t, Monoid m) =>
(a -> m) -> t a -> m
F.foldMap Circle n -> QDiagram b V2 n Any
forall n b.
(Renderable (Path V2 n) b, TypeableFloat n) =>
Circle n -> QDiagram b V2 n Any
drawCircle [Circle n]
cs

-- | Draw an Apollonian gasket: the first argument is the threshold;
--   the recursion will stop upon reaching circles with radii less than
--   it. The next three arguments are bends of three circles.
apollonianGasket :: (Renderable (Path V2 n) b, TypeableFloat n)
                 => n -> n -> n -> n -> QDiagram b V2 n Any
apollonianGasket :: forall n b.
(Renderable (Path V2 n) b, TypeableFloat n) =>
n -> n -> n -> n -> QDiagram b V2 n Any
apollonianGasket n
thresh n
b1 n
b2 n
b3 = [Circle n] -> QDiagram b V2 n Any
forall n b.
(Renderable (Path V2 n) b, TypeableFloat n) =>
[Circle n] -> QDiagram b V2 n Any
drawGasket ([Circle n] -> QDiagram b V2 n Any)
-> ([Circle n] -> [Circle n]) -> [Circle n] -> QDiagram b V2 n Any
forall b c a. (b -> c) -> (a -> b) -> a -> c
. n -> [Circle n] -> [Circle n]
forall n. RealFloat n => n -> [Circle n] -> [Circle n]
apollonian n
thresh ([Circle n] -> QDiagram b V2 n Any)
-> [Circle n] -> QDiagram b V2 n Any
forall a b. (a -> b) -> a -> b
$ (n -> n -> n -> [Circle n]
forall n. RealFloat n => n -> n -> n -> [Circle n]
initialConfig n
b1 n
b2 n
b3)

------------------------------------------------------------
-- Some notes on floating-point error
-- (only for the intrepid)
------------------------------------------------------------
{-

-- code from Gerald Gutierrez, personal communication

module Main where

import           Data.Complex
import           Diagrams.Backend.SVG.CmdLine
import           Diagrams.Prelude

-- ------^---------^---------^---------^---------^---------^---------^--------

data Circle = Circle Double (Complex Double) deriving (Show)

descartes a b c
  = (s + r, s - r)
  where
    s = a + b + c
    r = 2 * sqrt ((a * b) + (b * c) + (c * a))

descartesDual a b c d
  = 2 * (a + b + c) - d

soddies (Circle k1 b1) (Circle k2 b2) (Circle k3 b3)
  = ( Circle k4 b4
    , Circle k5 b5 )
  where
    (k4, k5) = descartes k1 k2 k3
    (b4, b5) = descartes b1 b2 b3

soddiesDual (Circle k1 b1) (Circle k2 b2) (Circle k3 b3) (Circle k4 b4)
  = Circle (descartesDual k1 k2 k3 k4) (descartesDual b1 b2 b3 b4)

mutuallyTangentCirclesFromTriangle z1 z2 z3
  = ( Circle k1 (z1 * (k1 :+ 0))
    , Circle k2 (z2 * (k2 :+ 0))
    , Circle k3 (z3 * (k3 :+ 0)) )
  where
    a = magnitude (z2 - z3)
    b = magnitude (z3 - z1)
    c = magnitude (z1 - z2)
    s = (a + b + c) / 2
    k1 = 1 / (s - a)
    k2 = 1 / (s - b)
    k3 = 1 / (s - c)

main :: IO ()
main = mainWith picture

pic = mainWith picture

mkCircle :: Circle -> Diagram B
mkCircle (Circle k b)
  = circle (1 / k) # moveTo (p2 (realPart z, imagPart z))
  where
    z = b / (k :+ 0)

picture :: Diagram B
picture
  = mkCircle c1 <> mkCircle c2 <> mkCircle c3 <> mkCircle c4 <> mkCircle c5
  where
    (c1, c2, c3) = mutuallyTangentCirclesFromTriangle z1 z2 z3
    (c4, c5) = soddies c1 c2 c3

    -- z1 = 0 :+ 0
    -- z2 = 3 :+ 0
    -- z3 = 0 :+ 4

    -- z1 = (-0.546) :+ (-0.755)
    -- z2 = ( 0.341) :+ (-0.755)
    -- z3 = (-0.250) :+ ( 0.428)

ofsBad = 0.15397 -- doesn't work
ofsGood = 0.15398 -- works

ofs = ofsGood

z1 = ofs + ((-0.546) :+ (-0.755))
z2 = ofs + (( 0.341) :+ (-0.755))
z3 = ofs + ((-0.250) :+ ( 0.428))


------------------------------------------------------------

Email to Gerald Gutierrez, 30 Sep 2016:

I got a chance to sit down and think hard about your code today, and
I think I have figured out what the problem is.

If you look at the outputs of 'soddies c1 c2 c3' using the two
different values for 'ofs', you will see that they are *almost* the
same, except that the complex numbers have been switched (though in
fact their real components are almost the same so you only notice the
difference in the imaginary components).  This clearly causes
incorrect results since the y-coordinates represented by the imaginary
components are now getting scaled by different bend values.  So the
resulting circles are of the right size and have (close to) the
correct x-coordinates, but their y-coordinates are wrong.

The problem is that in the 'soddies' function, you make independent
calls to 'descartes k1 k2 k3' (to solve for the bends) and 'descartes
b1 b2 b3' (to solve for the bend-center products), and then *assume*
that they return their solutions *in the same order*, so that you can
match up the first values to make one circle and the second values to
make another circle.  I would guess that this is *usually* true but
apparently not when certain values are hovering right around a branch
cut of sqrt, or something like that.

I think my code in Diagrams.TwoD.Apollonian suffers from the same
problem.  Ultimately I think it is due to the inherent inaccuracy of
floating-point numbers; there is nothing wrong with the formula
itself.  It would be nice to figure out how to correct for this, but I
am not sure how off the top of my head.  The interesting thing is that
switching the bend-center products does not seem to violate the
generalized Descartes' Theorem at all --- so it would seem that it
does not actually completely characterize mutually tangent circles,
that is, there exist sets of circles satisfying the theorem which are
not in fact mutually tangent.

-}