Problem 199

199 Iterative Circle Packing
Problem 199 - Project Euler
Apollonian gasket っていうらしいです．フラクタル．
Apollonian gasket - Wikipedia
問題は3つの円に接する円の半径を求めることだと思いますが…
Descartes' theorem - Wikipedia
ということで，これが分かれば，あとは再帰的に隙間を埋めてやればよい．

import Data.Map (toList, fromListWith)
import Control.Arrow ((***))
import Data.List (sort)

type State = ((Double, Double, Double), Int)

next :: State -> ((Double, Int), [State])
next ((a,b,c),n) = ((r,n),zip s.repeat $ n)
    where r = 1 / (curv (1/a) (1/b) (1/c))
          s = [(r,b,c),(r,a,c),(r,a,b)]

curv :: (Floating a) => a -> a -> a -> a
curv a b c = a+b+c + 2*sqrt (a*b+b*c+c*a)

gather :: Ord k => [(k, Int)] -> [(k, Int)]
gather = toList.fromListWith (+)

step :: ([(Double, Int)], [State]) -> ([(Double, Int)], [State])
step (_,ss) = (gather***gather.concat).unzip.map next $ ss

p199 :: Int -> Double
p199 m = (pi-cover) / pi
    where area (r,n) = r*r*pi*fromIntegral n
          cover = sum.sort.map area.concat.take (m+1).map fst.iterate step $ ini
          ini = ([(r1,3)],[(s1,1),(s2,3)])
          s1 = (r1, r1, r1)
          s2 = (-1, r1, r1)
          r1 = 2*sqrt 3 - 3

一応DPらしきことをして，計算結果の再利用をしてみた．