Problem 198

Problem 198 - Project Euler
ちょっとしたミスでした．基本は途中経過で述べたこと．
以下コード．

exGcd :: (Integral a) => a -> a -> ((a, a), a)
exGcd x y | y == 0    = ((1,0),x)
          | otherwise = let (q,r) = divMod x y
                            ((a,b),d) = exGcd y r
                        in ((b,a-b*q),d)

invMod :: (Integral b) => b -> b -> b
invMod n = fst.fst.exGcd n

farey :: (Integral a) => a -> [(a, a)]
farey n = takeWhile ((<n).fst) fs
    where fs = (0,1):(1,n):zipWith tie fs (tail fs)
          tie (a,b) (c,d) = let k = div (n+b) d
                            in (k*c-a,k*d-b)

left,right :: Int -> Int -> (Int, Int) -> Int
left x b (p,q) | p == 0        = 0
               | x*p <= q      = u1 --p/q <= 1/x
               | min u1 u2 > 0 = min u1 u2
               | otherwise     = 0
    where s0 = (invMod p q) `mod` q
          u1 = (b - 2*q*s0) `div` (2*q^2)
          u2 = (x - 2*s0*(x*p-q)) `div` (2*q*(x*p-q))
right x b (p,q) | x*p >= q  = 0 -- p/q >= 1/x
                | l > 0     = u - l
                | otherwise = u
    where s0 = (- invMod p q) `mod` q
          l = (x - 2*s0*(q-x*p)) `div` (2*q*(q-x*p))
          u = (b - 2*q*s0) `div` (2*q^2)

ambiguous :: Int -> Int -> (Int, Int) -> Int
ambiguous x b (p,q) = left x b (p,q) + right x b (p,q)

p198 :: Int -> Int -> Int
p198 x b = sum.map (ambiguous x b).takeWhile small.farey $ n
    where n = floor.sqrt.fromIntegral $ div b 2
          small (p,q) = x*p < 2*q -- p/q < 2/x ?

main :: IO ()
main = print.p198 100 $ 10^8

計算量は10^6くらい(たぶん)．
まぁ，まぁ，速いんですが，コードに簡潔さが足りない気がする．
フォーラムに素晴しいコードがのってた．