32 lines
1.1 KiB
Haskell
32 lines
1.1 KiB
Haskell
-- For how many N <= 10000 does sqrt(N) have an odd period when written as a continued fraction?
|
|
import Euler
|
|
|
|
-- Term r c === r*sqrt(n)+c
|
|
data Term = Term (Ratio Integer) (Ratio Integer) deriving (Eq, Show)
|
|
|
|
-- 1/(r*sqrt(n) + c)
|
|
-- = (r*sqrt(n) - c) / (r^2*n - c^2)
|
|
-- = (r/(r^2*n - c^2))*sqrt(n) - c/(r^2*n - c^2)
|
|
termRecip :: Term -> Integer -> Term
|
|
termRecip (Term r c) n = Term (r / denom) (-c / denom)
|
|
where denom = r^2 * (fromIntegral n) - c^2
|
|
|
|
-- sqrt(n) = rem0
|
|
-- = a0 + 1/rem1
|
|
-- = a0 + 1/(a1 + 1/rem2)
|
|
-- = ...
|
|
-- = a0 + 1/(a1 + 1/(a2 + 1/(a3 + ...)))
|
|
--
|
|
-- rem0 = sqrt(n), a0 = floor(rem0)
|
|
-- rem1 = 1/(rem0 - a0), a1 = floor(rem1)
|
|
-- rem2 = 1/(rem1 - a1), a2 = floor(rem2)
|
|
--
|
|
toFraction :: Integer -> ([Integer], [Integer])
|
|
toFraction n = let sqrtN = sqrt $ fromIntegral n in
|
|
if floor sqrtN ^ 2 == n then ([floor sqrtN], [])
|
|
else flip unfoldPeriodicState (Term 1 0) $ state $ \(Term r c) ->
|
|
let a = floor (fromRational r * sqrtN + fromRational c) :: Integer in
|
|
(Just a, termRecip (Term r (c - fromIntegral a)) n)
|
|
|
|
main = print $ length $ filter (odd . length . snd . toFraction) [1..10000]
|