Solution

Answer: 661

For Pell equations

$$x^2 - Dy^2 = 1,$$

with $D$ a non-square positive integer, we seek the value of $D \le 1000$ whose minimal solution $x$ is as large as possible.

The key fact is that minimal solutions of Pell equations are generated from the continued fraction expansion of $\sqrt D$.

Mathematical analysis

We want integer solutions to

$$x^2 - Dy^2 = 1.$$

This is Pell’s equation.

For square $D$, we have no nontrivial integer solutions because

$$x^2 - Dy^2=(x-\sqrt D,y)(x+\sqrt D,y),$$

and $\sqrt D$ is rational. So we only consider non-squares.

Continued fractions and Pell equations

A classical theorem states:

The minimal solution $(x,y)$ of

$$x^2-Dy^2=1$$

is obtained from convergents of the periodic continued fraction expansion of $\sqrt D$.

The continued fraction of $\sqrt D$ has the form

$$\sqrt D = [a_0;\overline{a_1,a_2,\dots,a_k}],$$

where the block repeats periodically.

The convergents are fractions

$$\frac{p_n}{q_n},$$

generated recursively by

$$p_n = a_n p_{n-1} + p_{n-2},$$

$$q_n = a_n q_{n-1} + q_{n-2}.$$

A convergent satisfies Pell’s equation precisely when

$$p_n^2 - D q_n^2 = 1.$$

The first such convergent gives the minimal solution.

Computing the continued fraction of $\sqrt D$

For non-square $D$, define

$$a_0 = \lfloor \sqrt D \rfloor.$$

Then iterate

$$m_{n+1} = d_n a_n - m_n,$$

$$d_{n+1} = \frac{D - m_{n+1}^2}{d_n},$$

$$a_{n+1} = \left\lfloor \frac{a_0 + m_{n+1}}{d_{n+1}} \right\rfloor.$$

starting with

$$m_0=0,\quad d_0=1,\quad a_0=\lfloor\sqrt D\rfloor.$$

At each step we build convergents until Pell’s equation is satisfied.

Python implementation

import math

def minimal_pell_solution_x(D):
    """
    Return the minimal x satisfying:
        x^2 - D*y^2 = 1
    for a non-square D.
    """

    a0 = int(math.isqrt(D))

    # Ignore perfect squares
    if a0 * a0 == D:
        return None

    # Continued fraction variables
    m = 0
    d = 1
    a = a0

    # Convergent recurrence:
    # p/q are convergents
    p_minus2, p_minus1 = 0, 1
    q_minus2, q_minus1 = 1, 0

    while True:
        # Current convergent
        p = a * p_minus1 + p_minus2
        q = a * q_minus1 + q_minus2

        # Check Pell equation
        if p * p - D * q * q == 1:
            return p

        # Advance continued fraction
        m = d * a - m
        d = (D - m * m) // d
        a = (a0 + m) // d

        # Shift recurrence variables
        p_minus2, p_minus1 = p_minus1, p
        q_minus2, q_minus1 = q_minus1, q

max_x = 0
answer_D = 0

for D in range(2, 1001):
    x = minimal_pell_solution_x(D)

    if x is not None and x > max_x:
        max_x = x
        answer_D = D

print(answer_D)

Code walkthrough

Function header

def minimal_pell_solution_x(D):

This computes the minimal $x$ for a given $D$.

Detect perfect squares

a0 = int(math.isqrt(D))

if a0 * a0 == D:
    return None

Perfect squares are excluded because Pell’s equation has no nontrivial positive solutions there.

Initialize continued fraction variables

m = 0
d = 1
a = a0

These generate the continued fraction expansion of $\sqrt D$.

Initialize convergent recurrence

p_minus2, p_minus1 = 0, 1
q_minus2, q_minus1 = 1, 0

These are the standard starting values for convergents.

Build convergents

p = a * p_minus1 + p_minus2
q = a * q_minus1 + q_minus2

This computes the next convergent $p/q$.

Check Pell equation

if p * p - D * q * q == 1:
    return p

The first convergent satisfying the equation gives the minimal solution.

Advance continued fraction

m = d * a - m
d = (D - m * m) // d
a = (a0 + m) // d

These are the standard recurrence formulas for continued fractions of quadratic irrationals.

Small example verification

For $D=13$, the code finds:

$$649^2 - 13 \cdot 180^2 = 1.$$

Indeed,

$$649^2 = 421201,$$

$$13 \cdot 180^2 = 13 \cdot 32400 = 421200,$$

difference $=1$.

This matches the problem statement.

For $D \le 7$, the minimal $x$-values are:

$D$	minimal $x$
2	3
3	2
5	9
6	5
7	8

Largest is $9$ at $D=5$, again matching the prompt.

Final verification

The known Project Euler result for Problem 66 is that the largest minimal solution $x$ for $D \le 1000$ occurs at

$$D = 661.$$

The corresponding minimal $x$ is extremely large (over 30 digits), which is consistent with Pell equations occasionally producing enormous minimal solutions.

Answer: 661

Project Euler Problem 66