Problem 66: Diophantine Equation

Mathematical Analysis

Pell’s Equation

The equation $x^2 - Dy^2 = 1$ is known as Pell’s equation. For non-square $D > 0$, this equation always has infinitely many solutions in positive integers.

Connection to Continued Fractions

The fundamental (minimal) solution $(x_1, y_1)$ of Pell’s equation is obtained from the convergents of the continued fraction expansion of $\sqrt{D}$.

Theorem: The continued fraction expansion of $\sqrt{D}$ is periodic:

where $r$ is the period length.

Theorem: Let $p_k / q_k$ denote the $k$-th convergent of the continued fraction of $\sqrt{D}$. Then the fundamental solution is:

• If $r$ is even: $(x_1, y_1) = (p_{r-1}, q_{r-1})$
• If $r$ is odd: $(x_1, y_1) = (p_{2r-1}, q_{2r-1})$

Convergent Computation

The convergents satisfy the recurrence:

with $p_{-1} = 1, p_{-2} = 0$ and $q_{-1} = 0, q_{-2} = 1$.

Continued Fraction Algorithm for $\sqrt{D}$

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

Recurrence:

The period ends when $a_n = 2a_0$.

Proof of Correctness

Claim: The fundamental solution of $x^2 - Dy^2 = 1$ is given by $(p_{r-1}, q_{r-1})$ where $r$ is the period (or $2r$ if $r$ is odd).

Proof sketch:

1. By Lagrange’s theorem, $\sqrt{D}$ has a periodic continued fraction.
2. For any convergent $p_k/q_k$, we have $|p_k/q_k - \sqrt{D}| < 1/(q_k q_{k+1})$.
3. If $p_k^2 - D q_k^2 = 1$, then $(p_k + q_k\sqrt{D})(p_k - q_k\sqrt{D}) = 1$.
4. The theory of continued fractions guarantees that the first such solution occurs at the end of the first period (or second period if the period is odd).

Editorial

We scan all non-square values of $D$ up to $1000$. For each one, we compute the period of the continued fraction of $\sqrt{D}$ and use its parity to determine which convergent yields the fundamental solution of Pell’s equation. We then generate convergents up to that index, extract the minimal numerator $x$, and keep the value of $D$ for which this minimal $x$ is largest.

Pseudocode

Initialize the best pair as
    current best D = 0
    current best minimal x = 0

For each integer D from 2 through 1000:
    if D is a perfect square, move on to the next value

    compute the period length of the continued fraction of sqrt(D)
    from that parity, determine which convergent carries the fundamental solution

    regenerate the continued fraction up to the required convergent
    read off its numerator x

    if x is larger than the current best minimal x:
        replace the recorded best pair by (D, x)

Return the recorded value of D.

Complexity

• Time: $O(N \cdot P_{\max})$ where $N = 1000$ and $P_{\max}$ is the maximum period length of any continued fraction for $D \le 1000$.
  • The period length is $O(\sqrt{D} \log D)$ in the worst case.
• Space: $O(1)$ per $D$ (only need to track current and previous convergents).

Answer

#include <bits/stdc++.h>
using namespace std;

// Big integer class (base 10^9) for Pell's equation
// Needed since minimal x can be ~10^37
struct BigInt {
    vector<int> d; // digits in base 10^9, least significant first

    BigInt() {}
    BigInt(long long v) {
        if (v == 0) d.push_back(0);
        while (v > 0) { d.push_back(v % 1000000000); v /= 1000000000; }
    }

    bool operator>(const BigInt& o) const {
        if (d.size() != o.d.size()) return d.size() > o.d.size();
        for (int i = (int)d.size()-1; i >= 0; i--)
            if (d[i] != o.d[i]) return d[i] > o.d[i];
        return false;
    }

    BigInt operator+(const BigInt& o) const {
        BigInt res;
        int carry = 0;
        for (int i = 0; i < (int)max(d.size(), o.d.size()) || carry; i++) {
            long long s = carry;
            if (i < (int)d.size()) s += d[i];
            if (i < (int)o.d.size()) s += o.d[i];
            res.d.push_back(s % 1000000000);
            carry = s / 1000000000;
        }
        return res;
    }

    BigInt operator*(long long v) const {
        BigInt res;
        long long carry = 0;
        for (int i = 0; i < (int)d.size() || carry; i++) {
            long long cur = carry;
            if (i < (int)d.size()) cur += (long long)d[i] * v;
            res.d.push_back(cur % 1000000000);
            carry = cur / 1000000000;
        }
        while (res.d.size() > 1 && res.d.back() == 0) res.d.pop_back();
        return res;
    }
};

// Solve Pell's equation x^2 - D*y^2 = 1 using continued fractions
// Returns the fundamental (minimal positive) solution x
BigInt solve_pell(int D) {
    int a0 = (int)sqrt((double)D);
    if (a0 * a0 == D) return BigInt(0); // perfect square

    // Find period length of continued fraction of sqrt(D)
    int m = 0, d = 1, a = a0;
    int period = 0;
    do {
        m = d * a - m;
        d = (D - m * m) / d;
        a = (a0 + m) / d;
        period++;
    } while (a != 2 * a0);

    // If period is even, solution is convergent p_{period-1}
    // If period is odd, solution is convergent p_{2*period-1}
    int target = (period % 2 == 0) ? period - 1 : 2 * period - 1;

    // Compute convergent p_{target}
    m = 0; d = 1; a = a0;
    BigInt pm2(1), pm1(a0); // p_{-1}=1, p_0=a0

    if (target == 0) return pm1;

    for (int i = 1; i <= target; i++) {
        m = d * a - m;
        d = (D - m * m) / d;
        a = (a0 + m) / d;

        BigInt pn = pm1 * a + pm2;
        pm2 = pm1;
        pm1 = pn;
    }

    return pm1;
}

int main() {
    ios_base::sync_with_stdio(false);

    BigInt max_x(0);
    int best_D = 0;

    for (int D = 2; D <= 1000; D++) {
        int s = (int)sqrt((double)D);
        if (s * s == D) continue;

        BigInt x = solve_pell(D);
        if (x > max_x) {
            max_x = x;
            best_D = D;
        }
    }

    cout << best_D << endl;

    return 0;
}

"""
Project Euler Problem 66: Diophantine Equation (Pell's Equation)

Find the value of D <= 1000 for which the minimal solution x of
x^2 - D*y^2 = 1 is maximized.

Uses continued fraction expansion of sqrt(D) to find fundamental solutions.
"""

from math import isqrt

def solve_pell(D):
    """
    Solve x^2 - D*y^2 = 1 (Pell's equation) for the fundamental solution.
    Returns the minimal positive x using continued fraction expansion of sqrt(D).
    """
    a0 = isqrt(D)
    if a0 * a0 == D:
        return None  # D is a perfect square, no solution

    # Find the period length of the continued fraction of sqrt(D)
    m, d, a = 0, 1, a0
    period = 0
    while True:
        m = d * a - m
        d = (D - m * m) // d
        a = (a0 + m) // d
        period += 1
        if a == 2 * a0:
            break

    # If period is even, solution is at convergent index (period - 1)
    # If period is odd, solution is at convergent index (2*period - 1)
    target = period - 1 if period % 2 == 0 else 2 * period - 1

    # Compute convergent p_{target}
    m, d, a = 0, 1, a0
    p_prev, p_curr = 1, a0  # p_{-1} = 1, p_0 = a0

    if target == 0:
        return p_curr

    for i in range(1, target + 1):
        m = d * a - m
        d = (D - m * m) // d
        a = (a0 + m) // d

        p_prev, p_curr = p_curr, a * p_curr + p_prev

    return p_curr

def main():
    max_x = 0
    best_D = 0

    for D in range(2, 1001):
        x = solve_pell(D)
        if x is not None and x > max_x:
            max_x = x
            best_D = D

    print(best_D)

if __name__ == "__main__":
    main()