Problem 198 – Project Euler

ちょっとしたミスでした．基本は途中経過で述べたこと．

以下コード． 

exGcd :: (Integral a) => a -> a -> ((a, a), a)
exGcd x y | y ==     = ((1,),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 = (,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 ==         = 
| x*p <= q      = u1 --p/q <= 1/x
| min u1 u2 >  = min u1 u2
| otherwise     = 
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  =  -- p/q >= 1/x
| l >      = 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くらい(たぶん)．

まぁ，まぁ，速いんですが，コードに簡潔さが足りない気がする．

フォーラムに素晴しいコードがのってた．