Diophantine Equation¶

SageMath can solve these Diophantine equations for us. For equations of the form $x^2 - Dy^2 = 1$ (which are called Pell equations), it will give parameterizations in $t$ for both $x$ and $y$. Here's $D=2$ as an example:

In [1]:
var('x,y')
solve_diophantine(x^2 - 2*y^2 == 1)
Out[1]:
[(sqrt(2)*(2*sqrt(2) + 3)^t - sqrt(2)*(-2*sqrt(2) + 3)^t + 3/2*(2*sqrt(2) + 3)^t + 3/2*(-2*sqrt(2) + 3)^t,
  -3/4*sqrt(2)*(2*sqrt(2) + 3)^t + 3/4*sqrt(2)*(-2*sqrt(2) + 3)^t - (2*sqrt(2) + 3)^t - (-2*sqrt(2) + 3)^t),
 (-sqrt(2)*(2*sqrt(2) + 3)^t + sqrt(2)*(-2*sqrt(2) + 3)^t - 3/2*(2*sqrt(2) + 3)^t - 3/2*(-2*sqrt(2) + 3)^t,
  3/4*sqrt(2)*(2*sqrt(2) + 3)^t - 3/4*sqrt(2)*(-2*sqrt(2) + 3)^t + (2*sqrt(2) + 3)^t + (-2*sqrt(2) + 3)^t)]

$t=0$ seems to correspond to the minimal solution we are after.

In [2]:
def minimal_solution(d):
    sols = solve_diophantine(x^2 - d*y^2 == 1)
    u, v = sols[0]
    return (abs(u(t=0)), abs(v(t=0)))

Now we can just iterate (although it is pretty slow).

In [3]:
max((d for d in range(1, 1001) if not is_square(d)), key=lambda k: minimal_solution(k)[0])
Out[3]:
661

Let's dig a little deeper so we can optimize.

Solving Pell equations¶

Lagrange proved that if $(x_0, y_0)$ is a solution to $$x^2 - dy^2 = 1$$ then $\frac{x_0}{y_0}$ is a convergent of the continued fraction of $\sqrt{d}$. Specifically, if $p$ is the period of the continued fraction, then the first solution will be the $(p-1)$th convergent if $p$ is even, and the $(2p-1)$th convergent if $p$ is odd.

This is great for us, since there are algorithms to compute these convergents. We'll use SageMath here; see problem 64 for how to compute the partial denominators of the continued fraction of a square root, and see problem 65 for an algorithm that uses partial denominators to compute convergents of continued fractions (SageMath's constructions make this implementation a little slow, but it makes the code easier to read - and it's still considerably faster than using solve_diophantine).

In [4]:
def pell_fundamental_solution(d):
    K.<s> = QuadraticField(d)
    f = continued_fraction(s)
    p = f.period_length()
    if p % 2 == 0:
        n = p - 1
    else:
        n = 2 * p - 1
        
    x, y = f.convergent(n).as_integer_ratio()
    return x, y

Now we'll just iterate with this function instead.

In [5]:
max((d for d in range(1, 1001) if not is_square(d)), key=lambda x: pell_fundamental_solution(x)[0])
Out[5]:
661

And in case you want to know the minimal $x$ for $x^2 - 661y^2 = 1$:

In [6]:
pell_fundamental_solution(661)
Out[6]:
(16421658242965910275055840472270471049, 638728478116949861246791167518480580)

Relevant sequences¶

  • Minimal values of $x$ for solutions to the Pell equation: A002350

Copyright (C) 2025 filifa¶

This work is licensed under the Creative Commons Attribution-ShareAlike 4.0 International license and the BSD Zero Clause license.