Problem 192 – Project Euler

最良近似有理数とでも呼ぶのでしょうか.

探索範囲と√nの形から,どうやら連分数と関係がありそうだ,と思っていたら,

Continued fraction – Wikipedia, the free encyclopedia

あとはSemiconvergentsの単調性を利用.

import Data.List (())
import Data.Ratio ((%),denominator)
sqrt' :: Integer -> Integer
sqrt' = floor.sqrt.fromIntegral
next :: (Integral a) => a -> a -> (a, (a, a)) -> (a, (a, a))
next x s (a,(p,q)) = let p' = a*q-p
q' = (x-p'^2) `div` q
in ((s+p')`div`q',(p',q'))
fractions :: Integer -> [Integer]
fractions x = map fst.iterate (next x.sqrt' $ x) $ (sqrt' x,(,1))
convergents :: Integer -> [(Integer, Integer)]
convergents x = cs
where cs = (1,):zipWith3 f (fractions x) cs ((,1):cs)
f a (n1,d1) (n2,d2) = (a*n1+n2,a*d1+d2)
bestDenominator :: Integer -> Integer -> Integer
bestDenominator b x
| abs ((n1%d1)^2 - x') < abs (s^2 - x') = d1
| otherwise = denominator s
where x' = fromIntegral x
((n1,d1):(n2,d2):_) = reverse.takeWhile ((<=b).snd).convergents $ x
s = (n1*((b-d2)`div`d1)+n2) % (d1*((b-d2)`div`d1)+d2)
p192 :: Integer -> Integer -> Integer
p192 b n = sum.map (bestDenominator b) $ [2..n]  (map (^2) [2..sqrt' n])
main :: IO ()
main = print.p192 (10^12) $ 100000