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