Problem 165

http://projecteuler.net/index.php?section=problems&id=165

この問題は結構苦労した。

線分交差列挙問題

といったら平面走査法とおもい、実装しようと試みた（前から実装してみたいとは思っていた）が…

無理

例外処理が予想以上に大変だ。

処理しなくてはいけないと思われる例外

走査線と平行な線分
端点での交差
3線分以上の交差

線分を蓄えるデータ構造（平衡木など）が必要になるが、

この要素（線分）の順序関係が動的に変化していくことと例外の存在が非常に厄介。

どこからか「そんな時は記号摂動だよ」という声が聞こえてくるがそんな気力はない。

というかそもそも平面走査法の計算量は交差数をkとすると O ((n+k)log n)　（実は O(nlog n + k)があるらしいが複雑らしい)。

つまり、線分が密なとき（k が n^2 程度 )は大して効率よくないよ。

この問題は [0,499]x[0,499]のなかに5000本の線分があるんだから密に決まってるよ。

というわけで O (n^2) のナイーブな解法を採用。

計算終わらない -> 有理数方から実数型
答えあわない -> 共有点の個数ではなく　異なる　共有点の個数だった。
なんか遅い -> 重複要素を取り除くのに quick sort みたいな　nub を使った (できたリストの順序は非保存、quick sort はもらってきた）

import Data.Maybe (mapMaybe)
import Data.Function (on)
import Control.Arrow ((***))
type Point = (Integer,Integer)
type Segment = (Point,Point,Integer,Integer,Integer)
type Intersection = (Double,Double)
blumBlumShub :: [Integer]
blumBlumShub = map (flip mod 500).tail.iterate sqMod $ 290797
where sqMod s = mod (s*s) 50515093
takeSegs :: Int -> [Integer] -> [Segment]
takeSegs  _ = []
takeSegs (n+1) (a:b:c:d:es) = segment (a,b) (c,d) : takeSegs n es
segment :: Point -> Point -> Segment
segment p@(x1,y1) q@(x2,y2) = (p,q,a,b,a*x1+b*y1)
where a = y1-y2
b = x2-x1
innerSeg (p1,q1,_,_,_) (x,y) = closure (fst p1,fst q1) x && closure (snd p1,snd q1) y && not endpoint
where closure (a,b) c = fromIntegral (min a b) <= c && c <= fromIntegral (max a b)
endpoint = same p1 (x,y) || same q1 (x,y)
same (a,b) (c,d) = fromIntegral a == c && fromIntegral b == d
overlapp (a,b) (c,d) = min a b < max c d || max a b > min b c
intersection :: Segment -> Segment -> Maybe Intersection
intersection l1@(p1,q1,a1,b1,d1) l2@(p2,q2,a2,b2,d2)
| overlappX && overlappY &&  det /=  &&
innerSeg l1 (x,y)  && innerSeg l2 (x,y) = Just (x,y)
| otherwise = Nothing
where overlappX = on overlapp (fst***fst) (p1,q1) (p2,q2)
overlappY = on overlapp (snd***snd) (p1,q1) (p2,q2)
det = fromIntegral $ a1*b2-a2*b1
x = fromIntegral (b2*d1-b1*d2) / det
y = fromIntegral (a1*d2-a2*d1) / det
choose2 :: [a] -> [(a, a)]
choose2 [] = []
choose2 (x:xs) = map ((,) x) xs ++ choose2 xs
main = print.length.quicknub.intersections $ 5000
intersections n = mapMaybe (uncurry intersection).choose2 $ takeSegs n blumBlumShub
-- quick sort like nub
quicknub :: Ord a => [a] -> [a]
quicknub xs = qnub compare xs []
qnub :: (a -> a -> Ordering) -> [a] -> [a] -> [a]
qnub _   []     r = r
qnub _   [x]    r = x:r
qnub cmp (x:xs) r = qpart cmp x xs [] [] r
qpart :: (a -> a -> Ordering) -> a -> [a] -> [a] -> [a] -> [a] -> [a]
qpart cmp x [] rlt rge r =
qnub cmp rlt (x:qnub cmp rge r)
qpart cmp x (y:ys) rlt rge r =
case cmp x y of
LT -> qpart cmp x ys rlt (y:rge) r
EQ -> qpart cmp x ys rlt rge r
GT -> qpart cmp x ys (y:rlt) rge r

メモリ足りなくなったり、苦労したんだけど、遅いし、きれいなコードじゃない。

コンテンツ