Pell's Equation

$x^2-D{y}^2=1,D>0\text{ not a perfect square.}$

Always has infinitely many integer solutions, generated from a single fundamental solution.

Continued fractions — the key tool

Every quadratic irrational $\sqrt{D}$ has a periodic continued fraction expansion:

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

and the fundamental solution appears among its convergents $p_k/q_k$.

// continued fraction expansion of sqrt(D)
vector<long long> cfSqrt(long long D) {
    long long a0 = (long long)sqrtl((long double)D);
    if (a0 * a0 == D) return {};                          // perfect square: no solution
    vector<long long> cf{a0};
    long long m = 0, d = 1, a = a0;
    while (a != 2 * a0) {
        m = d * a - m;
        d = (D - m * m) / d;
        a = (a0 + m) / d;
        cf.push_back(a);
    }
    return cf;                                            // period ends when a == 2*a0
}

The recurrence $(m,d,a)$ keeps everything in integers — no floating point beyond the initial $\lfloor \sqrt{D}\rfloor$.

The fundamental solution

Generate convergents until $p^2-D{q}^2=1$:

pair<long long,long long> pellFundamental(long long D) {
    vector<long long> cf = cfSqrt(D);
    long long p0 = 1, q0 = 0, p1 = cf[0], q1 = 1;
    for (int i = 1; ; i++) {
        if (p1 * p1 - D * q1 * q1 == 1) return {p1, q1};
        long long a = cf[i % (cf.size() - 1) + (i >= (int)cf.size() ? 1 : 0)];
        // (index into the periodic part; simplest is to extend cf lazily)
        long long p2 = a * p1 + p0, q2 = a * q1 + q0;
        p0 = p1; q0 = q1; p1 = p2; q1 = q2;
    }
}

The numbers explode

The fundamental solution can be enormous. For $D=61$ it is $(1766319049,226153980)$; for $D=991$, $x$ has 30 digits. Use __int128 or big integers, and be suspicious of any problem that asks for $D$ up to $10^4$ without a modulus.

Generating all solutions

If $(x_1,y_1)$ is fundamental, then

$x_k+y_k\sqrt{D}=(x_1+y_1\sqrt{D}{)}^k,$

which unrolls to the recurrence

$x_{k+1}=x_1{x}_k+D{y}_1{y}_k,y_{k+1}=x_1{y}_k+y_1{x}_k.$

For the $k$-th solution modulo something, this is a $2\times 2$ matrix power — $O(\log k)$:

$\left[\begin{array}{c}{x}_k\\ y_k\end{array}\right]={\left[\begin{array}{cc}{x}_1& D{y}_1\\ y_1& x_1\end{array}\right]}^{k-1}\left[\begin{array}{c}{x}_1\\ y_1\end{array}\right].$

This is how “find the $10^{18}$-th solution mod $p$” problems are solved.

The negative Pell equation

$x^2-D{y}^2=-1$

is solvable iff the continued fraction period of $\sqrt{D}$ has odd length. In that case the fundamental solution is the convergent at the end of the first period, and squaring it gives the fundamental solution of the positive equation.

Necessary (not sufficient) condition: $D$ has no prime factor $\equiv 3(\mathrm{mod}4)$.

Generalised Pell: $x^2-D{y}^2=N$

Harder. The solutions fall into finitely many classes, each generated by multiplying a fundamental class solution by powers of the Pell unit. The LMM (Lagrange-Matthews-Mollin) algorithm enumerates the classes. Rarely needed in contests.

Where it appears

Problem	Connection
”Smallest $x$ with $x^2-D{y}^2=1$“	directly
Square triangular numbers	$x^2-8y^2=1$
Almost-equilateral Heronian triangles	Pell-type
${\sum }_{i=1}^{n}i={\sum }_{i=n+1}^{m}i$	reduces to Pell
Best rational approximations to $\sqrt{D}$	the convergents themselves
Nearest integer to $(a+b\sqrt{D}{)}^n$	conjugate trick: the sum with its conjugate is an integer
Fibonacci-like recurrences from $\sqrt{5}$	the same algebraic-number machinery

The conjugate trick

Because $|x_1-y_1\sqrt{D}|<1$, the quantity $(x_1+y_1\sqrt{D}{)}^k+(x_1-y_1\sqrt{D}{)}^k$ is an integer and the second term is tiny. So

$\lfloor (x_1+y_1\sqrt{D}{)}^k\rfloor =\left(\text{that integer}\right)-1$

for appropriate parity. The integer satisfies a linear recurrence, so it is computable by matrix exponentiation — the standard way to handle “floor of an irrational power” problems.

See also: Euclidean Algorithm · Matrix Exponentiation · Linear Diophantine Equations