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: %f\n" i $ solve n cs

eps :: Double
eps = 1.0E-7

solve :: t -> [[Double]] -> Double
solve n cs = binSearch (feasible n cs') 0 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 0 _ = [[]]
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 0 $ 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

しかし，遅いぞ．これ．