Project Euler Solution 64: Odd period square roots

Code & Zahlen

Project Euler Problem 64: Odd period square roots teaches us about square root expansion as continued fractions.

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 {a3+ \dots}}} $$

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+ \dots}}}} $$

  • $a_0=4, \frac 1 {\sqrt{23}-4}=\frac {\sqrt{23}+4} 7=1+\frac {\sqrt{23}-3} 7$
  • $a_1=1, \frac 7 {\sqrt{23}-3}=\frac {7(\sqrt{23}+3)} {14}=3+\frac {\sqrt{23}-3} 2$
  • $a_2=3, \frac 2 {\sqrt{23}-3}=\frac {2(\sqrt{23}+3)} {14}=1+\frac {\sqrt{23}-4} 7$
  • $a_3=1, \frac 7 {\sqrt{23}-4}=\frac {7(\sqrt{23}+4)} 7=8+\sqrt{23}-4$
  • $a_4=8, \frac 1 {\sqrt{23}-4}=\frac {\sqrt{23}+4} 7=1+\frac {\sqrt{23}-3} 7$
  • $a_5=1, \frac 7 {\sqrt{23}-3}=\frac {7 (\sqrt{23}+3)} {14}=3+\frac {\sqrt{23}-3} 2$
  • $a_6=3, \frac 2 {\sqrt{23}-3}=\frac {2(\sqrt{23}+3)} {14}=1+\frac {\sqrt{23}-4} 7$
  • $a_7=1, \frac 7 {\sqrt{23}-4}=\frac {7(\sqrt{23}+4)} {7}=8+\sqrt{23}-4$

It can be seen that 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:

  • $\sqrt{2}=[1;(2)]$, period=1
  • $\sqrt{3}=[1;(1,2)]$, period=2
  • $\sqrt{5}=[2;(4)]$, period=1
  • $\sqrt{6}=[2;(2,4)]$, period=2
  • $\sqrt{7}=[2;(1,1,1,4)]$, period=4
  • $\sqrt{8}=[2;(1,4)]$, period=2
  • $\sqrt{10}=[3;(6)]$, period=1
  • $\sqrt{11}=[3;(3,6)]$, period=2
  • $\sqrt{12}=[3;(2,6)]$, period=2
  • $\sqrt{13}=[3;(1,1,1,1,6)]$, period=5

Exactly four continued fractions, for $N \le 13$, have an odd period.

How many continued fractions for $N \le 10\,000$ have an odd period?

We can use the given examples to write a test first:

def test_expansion() -> None:
    expand_root(2) == ([1], [2])
    expand_root(3) == ([1], [1, 2])
    expand_root(5) == ([2], [4])
    expand_root(6) == ([2], [2, 4])
    expand_root(7) == ([2], [1, 1, 1, 4])
    expand_root(8) == ([2], [1, 4])
    expand_root(10) == ([3], [6])
    expand_root(11) == ([3], [3, 6])
    expand_root(12) == ([3], [2, 6])
    expand_root(13) == ([3], [1, 1, 1, 1, 6])
    expand_root(23) == ([4], [1, 3, 1, 8])

The first part are the coefficients before the periodic part, the second list is the periodic part.

Now we need to think of actually building this function. Looking at the example, we can see how to proceed. We start with a number $N$. The $\sqrt N$ usually isn't an integer. We define an integer number $f$ such that $f \leq \sqrt N < f+1$ simply by using the floor function, $f := \lfloor \sqrt x \rfloor$. In the first step. The first step is writing $\sqrt N$ as $f + \sqrt N - f$. We know that $\sqrt N - f$ has to be smaller than 1 by definition. This way we can already split off $f$ and have the coefficient $a_0$ from the sequence already.

The next task is to further expand $\frac{1}{\sqrt N - f}$ into the a form $\frac{\sqrt N + f}{d}$. The denominator $d$ is easily found as $N - f^2$ using the third binomial equation. In general we will start from a form $\frac{b}{\sqrt N - c}$ and want to rewrite this as $b \frac{\sqrt N + c}{d}$.

The next step is the most difficult one. We need to pull out a full integer such that the fraction becomes smaller than one. For this we use that $f < \sqrt N < f + 1$ and therefore we can just see what we can do with the fraction $b \frac{f + c}{d}$. We first cancel $b$ and $d$ and then split the fraction into an integer $a_n$ and the new remaining part $\frac{\sqrt N - c'}{d'}$. In the next step we identify $c = - c'$ and $b = d'$ and do it again.

For the cancellation of the fraction we use greatest_common_denominator from Solution 33: Digit cancelling fractions. Together with the algorithm we then get this code, which isn't particularly pretty:

def expand_root(number: int) -> tuple[list[int], list[int]]:
    floor = int(math.sqrt(number))
    if floor**2 == number:
        return [floor], []
    results = [floor]
    states = [(1, floor)]
    c = floor
    b = 1
    while True:
        # print()
        assert c > 0
        # print(f"{b}/(sqrt({number}) - {c})")
        d = number - c**2
        gcd = greatest_common_denominator(b, d)
        # print(f"{b} (sqrt({number}) + {c})/{d}")
        b //= gcd
        d //= gcd
        # print(f"{b} (sqrt({number}) + {c})/{d}")
        split = (floor + c) // d
        a = split * b
        c -= split * d
        # print(f"{a} + {b} (sqrt({number}) + {c})/{d}")
        c = -c
        b = d
        state = (b, c)
        results.append(a)
        if state in states:
            break
        states.append(state)
    i = states.index(state) + 1
    return results[:i], results[i:]

It tracks the states of $b$ and $c$ such that it can detect a period. Once it has found a recurring element in the list of states, it can split off the prefix and the periodic part.

The solution to the problem then is just checking the length of the period for a range of numbers. Perfect squares will have a period part of zero, so they don't count towards the solution anyway and we don't need to exclude them.

def solution() -> int:
    result = 0
    for number in range(2, 10_000):
        beginning, period = expand_root(number)
        if len(period) % 2 == 1:
            result += 1
    return result

This produces the result within 265 ms, so that seems like the correct way to do it.