Project Euler Problem 64 Solution

Question

All square roots are periodic when written as continued fractions and can be written in the form:

N=a0+1a1+1a2+1a3+...\displaystyle \sqrt{N} = a_0 + \frac{1}{a_1 + \frac{1}{a_2 + \frac{1}{a_3 + ...}}}

For example, let us consider 23\sqrt{23}:

23=4+234=4+11234=4+11+2337\displaystyle \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:

23=4+11+13+11+18+...\displaystyle \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 23=[4;(1,3,1,8)]\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:

2=[1;(2)],period=13=[1;(1,2)],period=25=[2;(4)],period=16=[2;(2,4)],period=27=[2;(1,1,1,4)],period=48=[2;(1,4)],period=210=[3;(6)],period=111=[3;(3,6)],period=212=[3;(2,6)],period=213=[3;(1,1,1,1,6)],period=5\displaystyle \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 N13N \leq 13, have an odd period.

How many continued fractions for N10000N \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.136s
user   0m4.084s
sys    0m0.008s