Google Code Jam 2009 Round2 D (メモ)

書いてみた.

上位の方のソースを見ると2分探索で解いていた.

まね.

import Control.Monad
import Text.Printf
import Debug.Trace
-- main part
data Circle = C { x::Double, y::Double, r::Double}
main = do (c:_) <- getList :: IO [Int]
forM_ [1..c] $ i -> do
(n:_) <- getList :: IO [Int]
cs <- replicateM n getList :: IO [[Double]]
printf "Case #%d: %fn" i $ solve n cs
eps :: Double
eps = 1.0E-7
solve :: t -> [[Double]] -> Double
solve n cs = binSearch (feasible n cs')  eps
where cs' = map ([x,y,r] -> C x y r) cs
-- min x s.t. p x = True
--   where p is monotone (i.e x < y => p x < p y)
--         p l = False
binSearch :: (Fractional a, Ord a) => (a -> Bool) -> a -> a -> a
binSearch p l e = binpart (l + s) (l + 2 * s)
where s = until (x -> p (l + 2 * x)) (* 2) 1
binpart a b | b - a < e = m
| p m       = trace ("searching: " ++ show m) $ binpart a m
| otherwise = trace ("searching: " ++ show m) $ binpart m b
where m = (a + b) / 2
feasible :: t -> [Circle] -> Double -> Bool
feasible n cs s = any cover.choose 2 $ candidate n cs'
where cs' = map ((C x y r) -> C x y (s-r)) cs
cover [p,q] = all (c -> inCircle p c || inCircle q c) cs'
inCircle (u,v) c = (u - x c)^2 + (v - y c)^2 <= (r c)^2 + eps
candidate :: t -> [Circle] -> [(Double, Double)]
candidate n cs = cs' ++ concat [intersection c1 c2| [c1,c2] <- choose 2 cs]
where cs' = map ((C x y r) -> (x,y)) cs
choose :: (Integral t) => t -> [a] -> [[a]]
choose _ [] = []
choose  _ = [[]]
choose (n+1) (x:xs) = map (x:) (choose n xs) ++ choose (n+1) xs
intersection :: Circle -> Circle -> [(Double, Double)]
intersection c1 c2
| d > r c1 + r c2 +  eps      = [] -- too far each other
| d < abs (r c1 - r c2) + eps = [] -- one contain another
| d < eps                     = [] -- same center
| abs a > eps                 = ipart a b c (x c1) (y c1) (r c1)
| otherwise                   = map swap $ ipart b a c (y c1) (x c1) (r c1)
where a = 2 * (x c1 - x c2)
b = 2 * (y c1 - y c2)
c = (r c1)^2 - (x c1)^2 - (y c1)^2 - ( (r c2)^2 - (x c2)^2 - (y c2)^2 )
d = sqrt $ (a/2)^2 + (b/2)^2
swap (x,y) = (y,x)
-- solve a*x+b*y+c=0, (x-u)^2+(y-v)^2=r^2
ipart :: Double -> Double -> Double -> Double -> Double -> Double -> [(Double, Double)]
ipart a b c u v r = [(x1,y1), (x2,y2)]
where a' =  a^2 + b^2
b' = 2 * (a*b*u + b*c - a^2*v)
c' = c^2 + a^2*u^2 + 2*a*c*u + a^2*v^2 - a^2*r^2
d' = max  $ b'^2 - 4*a'*c'
y1 = (-b' + sqrt d') / (2*a')
x1 = (-b*y1 - c) / a
y2 = (-b' - sqrt d') / (2*a')
x2 = (-b*y2 - c) / a
-- output and input function
getList :: Read a => IO [a]
getList = liftM (map read.words) getLine

しかし,遅いぞ.これ.