199 Iterative Circle Packing

Problem 199 – Project Euler

Apollonian gasket っていうらしいです.フラクタル.

Apollonian gasket – Wikipedia, the free encyclopedia

問題は3つの円に接する円の半径を求めることだと思いますが…

Descartes’ theorem – Wikipedia, the free encyclopedia

ということで,これが分かれば,あとは再帰的に隙間を埋めてやればよい.

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らしきことをして,計算結果の再利用をしてみた.