Project Euler Solution 66: Diophantine equation
Project Euler Problem 66: Diophantine equation gives us a simple looking quadratic equation which turns out to have quite interesting structural properties.
Consider quadratic Diophantine equations of the form:
$$ x^2 – Dy^2 = 1 $$
For example, when $D=13$, the minimal solution in $x$ is $649^2 – 13\cdot 180^2 = 1$.
It can be assumed that there are no solutions in positive integers when $D$ is square.
By finding minimal solutions in $x$ for $D = {2, 3, 5, 6, 7}$, we obtain the following:
- $3^2 – 2\cdot2^2 = 1$
- $2^2 – 3\cdot1^2 = 1$
- $9^2 – 5\cdot4^2 = 1$
- $5^2 – 6\cdot2^2 = 1$
- $8^2 – 7\cdot3^2 = 1$
Hence, by considering minimal solutions in $x$ for $D \leq 7$, the largest $x$ is obtained when $D=5$.
Find the value of $D \leq 1000$ in minimal solutions of $x$ for which the largest value of $x$ is obtained.
What is missing from the problem statement but implicit in “Diophantine equations” is that $x$ and $y$ are both integers. For rational numbers, there would always be a solution. The equation $x^2 -D y^2 = 1$ with $x, y \in \mathbb N$ will always have a solution if $y^2 = (x^2 - 1)/D$ is a perfect square. We can therefore just iterate through all $x$ and check whether the expression $(x^2 - 1)/D$ is a perfect square.
To check whether it is a square is rather easy to implement:
def is_square(number: int) -> bool: floor = int(math.sqrt(number)) return floor**2 == number
And then we can write a solution like this:
def solution() -> int: max_x = 0 max_d = 0 for d in range(1, 1001): if is_square(d): continue for x in tqdm(itertools.count(1), desc=f"d={d}"): if is_square((x**2 - 1) / d): if x > max_x: max_x = x max_d = d break return max_d
In the beginning that works pretty well:
d=2: 1it [00:00, 27413.75it/s] d=3: 0it [00:00, ?it/s] d=5: 7it [00:00, 140479.08it/s] d=6: 3it [00:00, 65536.00it/s] d=7: 6it [00:00, 148910.20it/s] d=8: 1it [00:00, 31300.78it/s] d=10: 17it [00:00, 383350.37it/s] d=11: 8it [00:00, 204600.20it/s] d=12: 5it [00:00, 148734.18it/s] d=13: 647it [00:00, 1836072.18it/s] d=14: 13it [00:00, 288497.10it/s] d=15: 2it [00:00, 66576.25it/s] d=17: 31it [00:00, 577881.88it/s]
However, there are a few values of $D$ for which this doesn't work really sensibly:
d=59: 528it [00:00, 1949465.24it/s] d=60: 29it [00:00, 614317.25it/s] d=61: 335159610it [03:13, 1728850.74it/s] d=62: 61it [00:00, 1224174.85it/s] d=63: 6it [00:00, 299593.14it/s]
See how for $D=61$ it took 3:13 minutes to find the solution because we had to go up to $x = 335\,159\,612$ to find it. And then there are more such numbers like that. We need something more clever.
Reading up on Wikipedia on the Diophantine equation points us to the special case of Pell's equation which is exactly the form that we have. In the article it gives a relation to the approximation of square roots:
These solutions may be used to accurately approximate the square root of n by rational numbers of the form x/y.
In the section about solutions to the equation it says this:
Let $h_{i}/k_{i}$ denote the sequence of convergents to the regular continued fraction for $\sqrt {n}$. This sequence is unique. Then the pair $(x_{1},y_{1})$ solving Pell's equation and minimizing x satisfies $x_1 = h_i$ and $y_1 = k_i$ for some $i$. This pair is called the fundamental solution. Thus, the fundamental solution may be found by performing the continued fraction expansion and testing each successive convergent until a solution to Pell's equation is found.
This means that instead of having to iterate over all $x$, we can iterate over the convergents that we have programmed in Solution 65: Convergents of e. We combine this with the continued fractions of square roots from Solution 64: Odd period square roots. It is nice to see how we can directly reuse bits from previous problems and it nicely shows the progression within the problems.
Using this insight we can take the root expansion from Solution 64 and create an iterator which just loops over the periodic part:
def square_root_fraction_expansion(number: int) -> Iterator[int]: prefix, period = expand_root(number) yield from prefix while True: yield from period
Then we use that with the convergent generator from Solution 65 and then iterate through all the convergents to find the minimal solution in $x$ for a given $D$:
def minimal_solution(d: int) -> int: for x, y in convergents_series(square_root_fraction_expansion(number)): if x == 1: continue if x**2 - d * y**2 == 1: return x
It is important to use the explicit equation with integers here. The function is_square from above will fail for numbers larger than about $10^{17}$ as the FP64 data type then lacks the precision. Using the arbitrarily large integers makes this check possible.
To make sure that this really works, we verify that with a test using the examples from the problem statement:
def test_minimal_solution() -> None: assert minimal_solution(2) == 3 assert minimal_solution(3) == 2 assert minimal_solution(5) == 9 assert minimal_solution(6) == 5 assert minimal_solution(7) == 8
This checks out, so we can build our solution on top of that:
def solution() -> int: max_x = 0 max_d = 0 for d in range(1, 1001): if is_square(d): continue x = minimal_solution(d) if x > max_x: max_x = x max_d = d return max_d
This finds the solution within 34 ms, so that quite a difference to the brute force solution from above. This was a very interesting problem, and I have learned a lot while solving it!