Question
All square roots are periodic when written as continued fractions and can be written in the form:
\sqrt{N} = a_0 + \frac{1}{a_1 + \frac{1}{a_2 + \frac{1}{a_3 + ...}}}
For example, let us consider \sqrt{23}:
\sqrt{23} = 4 + \sqrt{23} - 4 = 4 + \frac{1}{\frac{1}{\sqrt{23} - 4}} = 4 + \frac{1}{1 + \frac{\sqrt{23} - 3}{7}}
If we continue we would get the following expansion:
\sqrt{23} = 4 + \frac{1}{1 + \frac{1}{3 + \frac{1}{1 + \frac{1}{8 + ...}}}}
The sequence is repeating. For conciseness, we use the notation \sqrt{23} = [4;(1,3,1,8)], to indicate that the block (1,3,1,8) repeats indefinitely.
The first ten continued fraction representations of (irrational) square roots are:
\begin{aligned} \sqrt{2} &= [1;(2)], \, \mathrm{period}=1 \\ \sqrt{3} &= [1;(1,2)], \, \mathrm{period}=2 \\ \sqrt{5} &= [2;(4)], \, \mathrm{period}=1 \\ \sqrt{6} &= [2;(2,4)], \, \mathrm{period}=2 \\ \sqrt{7} &= [2;(1,1,1,4)], \, \mathrm{period}=4 \\ \sqrt{8} &= [2;(1,4)], \, \mathrm{period}=2 \\ \sqrt{10} &= [3;(6)], \, \mathrm{period}=1 \\ \sqrt{11} &= [3;(3,6)], \, \mathrm{period}=2 \\ \sqrt{12} &= [3;(2,6)], \, \mathrm{period}=2 \\ \sqrt{13} &= [3;(1,1,1,1,6)], \, \mathrm{period}=5 \\ \end{aligned}
Exactly four continued fractions, for N \leq 13, have an odd period.
How many continued fractions for N \leq 10000 have an odd period?
Haskell
cycles :: Eq a => [a] -> [[a]]
cycles xs = map fst $ dropWhile (\(a, b) -> a ++ a /= b) $ zip cs (tail cs)
where cs = [left | n <- [1..], let (left, right) = splitAt n xs, left == take n right]
isSquare :: Int -> Bool
isSquare n = root == fromIntegral (round root)
where root = sqrt (fromIntegral n)
isPrime :: Int -> Bool
isPrime n | n < 1 = False
| otherwise = not $ or [n `rem` x == 0 | x <- [2..floor $ sqrt $ fromIntegral n]]
expansion :: Int -> [Int]
expansion s | isSquare s = []
| isSquare (s-1) = [fromIntegral $ 2 * (floor $ sqrt $ fromIntegral (s-1))]
| otherwise = head $ dropWhile (all (== first)) cs
where cs = cycles $ map a [1..]
first = (head . head) cs
m = (map m' [0..] !!)
m' 0 = 0
m' n = (d (n-1))*(a (n-1)) - (m (n-1))
d = (map d' [0..] !!)
d' 0 = 1
d' n = (s - (m n)^2) `quot` (d (n-1))
a = (map a' [0..] !!)
a' 0 = floor $ sqrt $ fromIntegral s
a' n = floor $ (fromIntegral ((a 0) + (m n))) / (fromIntegral (d n))
main :: IO ()
main = print $ length $ filter (\x -> odd $ length $ expansion x) [2..10000]$ ghc -O2 -o continued-fractions continued-fractions.hs
$ time ./continued-fractions
real 0m4.609s
user 0m4.593s
sys 0m0.016s