Detailed mathematical approach

Problem Summary

For each non-square integer \(D \le 1000\), Pell's equation

$x^2 - Dy^2 = 1$

has infinitely many positive integer solutions, and among them there is a unique fundamental solution with the smallest positive \(x\). The problem asks for the value of \(D\) whose fundamental solution has the largest \(x\). The search over \(D \le 1000\) ends at \(D = 661\).

A brute-force search in \(x\) or \(y\) is not realistic, because the minimal solution can already be enormous for moderate \(D\). The right viewpoint is to study the periodic continued fraction of \(\sqrt{D}\), because its convergents produce exactly the candidates that can solve Pell's equation.

Mathematical Approach

The implementations do not guess solutions numerically. They generate the continued fraction data for \(\sqrt{D}\), turn that data into convergents, and stop at the first convergent that satisfies the Pell identity.

Why continued fractions are the right object

If \(D\) is not a square, then \(\sqrt{D}\) is a quadratic irrational, so its simple continued fraction is periodic:

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

The convergents of this expansion are the best rational approximations to \(\sqrt{D}\) in a precise number-theoretic sense. For Pell's equation, that matters because if

$\frac{p}{q} \approx \sqrt{D},$

then

$p^2 - Dq^2$

is forced to be small. The remarkable theorem behind the method is that the fundamental solution of \(x^2 - Dy^2 = 1\) appears among these convergents, so scanning convergents is enough.

The recurring state for \(\sqrt{D}\)

The continued fraction is generated from the standard complete-quotient state. Start with

$a_0 = \lfloor \sqrt{D} \rfloor,\qquad m_0 = 0,\qquad d_0 = 1.$

At each stage define

$\alpha_n = \frac{\sqrt{D} + m_n}{d_n},\qquad a_n = \lfloor \alpha_n \rfloor.$

The next state is given by the recurrences used directly in the code:

$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.$

Two invariants make this work. First, \(d_n\) always divides \(D - m_n^2\), so the division is exact. Second, each complete quotient keeps the form \((\sqrt{D}+m_n)/d_n\). Because only finitely many such states can occur for fixed \(D\), the sequence eventually repeats, which is exactly the periodicity of the continued fraction of \(\sqrt{D}\).

Convergents and the first time the Pell identity becomes \(1\)

Once the coefficients \(a_n\) are available, the convergents are built from the usual two-step recurrence:

$p_{-1}=1,\quad p_0=a_0,\qquad q_{-1}=0,\quad q_0=1,$

$p_{n+1}=a_{n+1}p_n+p_{n-1},\qquad q_{n+1}=a_{n+1}q_n+q_{n-1}.$

Every convergent \(p_n/q_n\) is a better approximation to \(\sqrt{D}\), so the quantity

$p_n^2 - D q_n^2$

keeps landing on small integers such as \(-1\), \(1\), \(3\), or \(-4\). The first time it becomes exactly \(1\), the current numerator \(p_n\) is the fundamental \(x\). Classical theory says that this happens after one traversal of the repeating block when the period length is even and after two traversals when it is odd, but the implementations do not need to compute the period length explicitly. They simply keep iterating until the Pell test succeeds.

Worked example: \(D = 13\)

The problem statement already highlights \(D=13\), and it is a useful example because the first solution is not small. Here

$\sqrt{13} = [3; \overline{1,1,1,1,6}].$

The complete-quotient state begins

$ (m,d,a) = (0,1,3) \to (3,4,1) \to (1,3,1) \to (2,3,1) \to (1,4,1) \to (3,1,6) \to \cdots $

and the convergents are

$3,\ \frac{4}{1},\ \frac{7}{2},\ \frac{11}{3},\ \frac{18}{5},\ \frac{119}{33},\ \frac{137}{38},\ \frac{256}{71},\ \frac{393}{109},\ \frac{649}{180}.$

Evaluating the Pell expression shows the progression

$-4,\ 3,\ -3,\ 4,\ -1,\ 4,\ -3,\ 3,\ -4,\ 1.$

So the first hit of \(1\) occurs at

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

which means the fundamental solution is \((x,y)=(649,180)\). This is exactly the kind of loop the implementations perform for every non-square \(D\).

How the Code Works

The C++, Python, and Java implementations scan \(D=2,3,\dots,1000\). If \(D\) is a perfect square, that case is skipped immediately, because \(x^2 - Dy^2 = 1\) would then factor in a way that leaves only the trivial \(y=0\) solution and no meaningful Pell search.

For each non-square \(D\), the implementation stores only the current continued-fraction state \((m,d,a)\) and the previous two convergents. After each update it checks whether the new pair satisfies \(p^2 - Dq^2 = 1\). The first success is returned as the minimal \(x\) for that \(D\). This avoids storing the entire period and avoids any separate parity logic.

Finally, the implementation compares that \(x\) against the best one seen so far and keeps the corresponding \(D\). The built-in sanity checks verify two small prefixes of the search: for \(D \le 7\) the answer is \(5\), and for \(D \le 13\) the answer is \(13\).

Complexity Analysis

Let \(\ell_D\) be the number of continued-fraction coefficients processed for a given non-square \(D\) before the first solution is found. Then that case performs \(O(\ell_D)\) recurrence steps. The state variables \(m\), \(d\), and \(a\) stay small, but the convergents grow with the size of the fundamental solution, so the dominant work is big-integer arithmetic on numbers with \(O(\log x_D)\) bits.

Across all \(D \le 1000\), the total running time is the sum of those per-\(D\) costs and is easily practical for this range. Memory usage is \(O(\log x_{\max})\) bits, because the algorithm keeps only a constant number of big integers at any moment rather than the whole history of coefficients or convergents.

Footnotes and References

1. Project Euler Problem 66: https://projecteuler.net/problem=66
2. Pell's equation: Wikipedia - Pell's equation
3. Continued fraction: Wikipedia - Continued fraction
4. Periodic continued fraction: Wikipedia - Periodic continued fraction
5. Convergent of a continued fraction: Wikipedia - Convergent
6. Pell equation overview: MathWorld - Pell Equation