continuedFraction :: Integral a => a -> a -> a -> [a]
continuedFraction root num denom = whole : (continuedFraction root newdenom outdenom)
where whole = floor (((sqrt $ fromIntegral root) + (fromIntegral num)) / (fromIntegral denom))
newdenom = (denom * whole) - num
outdenom = (root - newdenom*newdenom) `div` denom
evaluateConvergent :: Integral a => [a] -> (a, a)
evaluateConvergent terms = foldr (\ term (x, y) -> (y + term * x, x)) (1, 0) terms
getConvergents :: Integral a => a -> [(a, a)]
getConvergents root = [evaluateConvergent $ take x $ continuedFraction root 0 1 | x <- [1..]]
getX :: Integral a => a -> a
getX d = fst $ head $ filter (\ (x, y) -> x*x - d*y*y == 1) $ getConvergents d
isPerfectSquare :: Integral a => a -> Bool
isPerfectSquare val = val == squareroot * squareroot
where squareroot = floor $ sqrt $ fromIntegral val
argmax :: Ord b => (a -> b) -> [a] -> a
argmax func list = fst $ foldl1 (\ (x1, y1) (x2, y2) -> if y1 > y2 then (x1, y1) else (x2, y2)) $ zip list (map func list)
main = print $ argmax getX (filter (not . isPerfectSquare) [2..1000])
And there you have it. Solution runs in about 0.4 seconds on my laptop.
No comments:
Post a Comment