Pell's Equation and Continued Fractions

Can the equation x2−2y2=1x^2 - 2y^2 = 1 have infinitely many whole-number solutions? If so, how would you find them all? A little experimentation turns up (x,y)=(3,2)(x,y) = (3,2), then (17,12)(17,12), then (99,70)(99,70). The numbers grow fast, and no pattern is obvious from guessing alone. This is Pell's equation, x2−dy2=1x^2 - dy^2 = 1 for a nonsquare positive integer dd. It is one of the oldest and most productive equations in number theory. It was studied by Brahmagupta and Bhāskara II centuries before Pell. Later it was placed on rigorous footing by Euler, Lagrange, and Legendre. The equation matters because it sits at the crossroads of several major theories. It is the norm-one condition for the ring Z[d]\mathbb{Z}[\sqrt d]. It governs the unit group of real quadratic fields. It is also intimately tied to the theory of Equivalence and Reduction of Binary Forms. This is because solutions to Pell's equation correspond exactly to automorphisms of the form x2−dy2x^2 - dy^2. Historically, Pell-type equations were used to produce excellent rational approximations to d\sqrt d. They also helped solve cattle-counting and Diophantine puzzles. In the Archimedes cattle problem, they even encoded staggeringly large numbers. The big picture is that Pell's equation is never solved by algebra alone. It requires an analytic input: the continued fraction expansion of d\sqrt d. This expansion is always eventually periodic, a fact from Continued Fractions of Quadratic Irrationals. The convergents of that expansion, truncated exactly at the end of each period, produce the smallest positive solution. This is called the fundamental solution. Every other solution is then generated from the fundamental one by a simple recursive multiplication rule. This mirrors the structure of a group of units. This article gives the full machinery. It covers the definition of the equation and its solutions, and the precise theorem connecting continued fraction convergents to solutions. It includes a complete proof, and two worked examples with genuinely different behavior (period length 1 and period length 4). It closes with the classical connection to binary quadratic forms. It also covers modern uses in cryptography and computational number theory.

Key takeaways
  • Pell's equation is x2−dy2=1x^2 - dy^2 = 1 with dd a positive nonsquare integer; solutions in positive integers always exist.
  • The continued fraction expansion of d\sqrt d is purely periodic after one term, and its convergents pn/qnp_n/q_n solve pn2−dqn2=±1p_n^2 - dq_n^2 = \pm 1.
  • The smallest positive solution (x1,y1)(x_1,y_1), the fundamental solution, arises at the convergent just before the period repeats.
  • All positive solutions are (xn,yn)(x_n, y_n) with xn+ynd=(x1+y1d)nx_n + y_n\sqrt d = (x_1+y_1\sqrt d)^n, generated by simple recursion.
  • The equation is the norm-one problem for Z[d]\mathbb{Z}[\sqrt d] and encodes the unit group of a real quadratic order.
  • Key idea: periodicity of the continued fraction of d\sqrt d forces convergents to hit x2−dy2=±1x^2-dy^2=\pm1 exactly, and the first hit of +1+1 generates every solution.

Where this lives today

Pell's equation is the norm-one problem for the ring of integers (or an order) in a real quadratic field Q(d)\mathbb{Q}(\sqrt d), and the fundamental solution (x1,y1)(x_1,y_1) generates the group of units — this is the rank-one case of Dirichlet's unit theorem, one of the cornerstones of algebraic number theory. The same continued-fraction machinery underlies the connection to Equivalence and Reduction of Binary Forms: automorphisms of the principal form x2−dy2x^2-dy^2 (matrices fixing the form under substitution) correspond exactly to solutions of Pell's equation, and this controls the structure of the class group of forms of discriminant 4d4d. Beyond pure number theory, Pell's equation supplies famously good rational approximations to square roots (each convergent pn/qnp_n/q_n near a period boundary approximates d\sqrt d to within O(1/qn2)O(1/q_n^2), historically used before calculators existed). It appears in recreational and computational problems such as the Archimedes cattle problem, whose solution requires x2−4729494y2=1x^2 - 4729494y^2=1, with a fundamental solution having hundreds of thousands of digits. In cryptography and computational number theory, Pell equations and continued fraction algorithms are close cousins of the techniques used in integer factorization (e.g., the continued fraction factoring method, a precursor to the quadratic sieve) and in computing regulators of real quadratic fields, a task relevant to class number computations and to some lattice-based cryptographic constructions.

:::quiz[For d=2d=2, the continued fraction is 2=[1;2‾]\sqrt2=[1;\overline{2}], period r=1r=1 (odd). What is the fundamental solution of x2−2y2=1x^2-2y^2=1?]

  • (1, 1)
  • (2, 1)
  • (3, 2) ✓
  • (7, 5)

Ask the tutor: Can you walk through why the chakravala method and the continued fraction method always produce the same fundamental solution?

DefinitionPell's Equation and Its Solutions

Let dd be a positive integer that is not a perfect square. Pell's equation is

x2−dy2=1,x^2 - d y^2 = 1,

to be solved in integers x,yx, y. A solution (x,y)(x,y) is called positive if x>0x > 0 and y>0y > 0; among positive solutions, the one with smallest xx (equivalently smallest yy) is the fundamental solution, denoted (x1,y1)(x_1, y_1). Solutions correspond to elements x+yd∈Z[d]x + y\sqrt d \in \mathbb{Z}[\sqrt d] of norm one, where the norm is N(x+yd)=(x+yd)(x−yd)=x2−dy2N(x+y\sqrt d) = (x+y\sqrt d)(x - y\sqrt d) = x^2 - dy^2. This correspondence is a group isomorphism-flavored fact: if x2−dy2=1x^2 - dy^2=1 and x′2−dy′2=1x'^2 - dy'^2 = 1, then multiplying α=x+yd\alpha = x+y\sqrt d and α′=x′+y′d\alpha' = x'+y'\sqrt d gives αα′=X+Yd\alpha\alpha' = X + Y\sqrt d with X2−dY2=N(α)N(α′)=1X^2 - dY^2 = N(\alpha)N(\alpha') = 1, so solutions compose. The continued fraction convergents of d\sqrt d are the fractions pn/qnp_n/q_n obtained by truncating the continued fraction expansion

d=[a0;a1,a2,a3,…]\sqrt d = [a_0; a_1, a_2, a_3, \ldots]

at each term, via the standard recursions p−1=1,p0=a0p_{-1}=1, p_0 = a_0, q−1=0,q0=1q_{-1}=0, q_0=1, and

pn=anpn−1+pn−2,qn=anqn−1+qn−2(n≥1).p_n = a_n p_{n-1} + p_{n-2}, \qquad q_n = a_n q_{n-1} + q_{n-2} \quad (n \ge 1).

Intuition

The key fact underlying everything is that d\sqrt d is a quadratic irrational, so (as established in the study of Continued Fractions of Quadratic Irrationals) its continued fraction expansion is eventually periodic. In fact for d\sqrt d it takes the clean form

d=[a0;a1,a2,…,ar‾],\sqrt d = [a_0; \overline{a_1, a_2, \ldots, a_r}],

purely periodic after the first term, where a0=⌊d⌋a_0 = \lfloor \sqrt d\rfloor and rr is the period length. Truncating this expansion at longer and longer prefixes produces convergents pn/qnp_n/q_n that approximate d\sqrt d extremely well — well enough that pn2−dqn2p_n^2 - dq_n^2 is forced to be a small integer, not just any real number close to 00. Periodicity means the "shape" of the tail of the expansion repeats, and this repetition is exactly what forces pn2−dqn2p_n^2 - dq_n^2 to cycle through a short list of values, hitting +1+1 at predictable moments (the ends of periods, or every other period if rr is odd). Picture the values pn2−dqn2p_n^2 - dq_n^2 as you scan through n=0,1,2,…n = 0, 1, 2, \ldots: they oscillate in sign and magnitude, bounded by 2d2\sqrt d, and periodicity of the continued fraction forces this sequence itself to become periodic. The first time the sequence returns to +1+1 (other than the trivial n=−1n=-1 case) gives you the fundamental solution — the smallest "building block" — and every other solution is a power of it, just as every unit in a number ring is a power of a fundamental unit times a root of unity.

:::svg[Convergents of 2=[1;2‾]\sqrt2=[1;\overline{2}] closing in, with pn2−2qn2p_n^2-2q_n^2 cycling through −1,+1,−1,+1,…-1,+1,-1,+1,\ldots] n value +1 -1 +1 -1 +1 -1 n=0 n=1 n=2 n=3 n=4 n=5 p_n²−2q_n² oscillates ±1 (period r=1, odd)

Scope, edge cases & common traps

TheoremPell's Equation via Continued Fraction Convergents

Let dd be a positive integer that is not a perfect square, and let

d=[a0;a1,…,ar‾]\sqrt d = [a_0; \overline{a_1, \ldots, a_r}]

be its (purely periodic after a0a_0) continued fraction expansion, with minimal period length r≥1r \ge 1. Let pn/qnp_n/q_n denote the convergents, n=−1,0,1,2,…n = -1, 0, 1, 2, \ldots.

(1) Convergent identity. For every n≥0n \ge 0,

pn2−dqn2=(−1)n+1Qn+1,p_n^2 - d q_n^2 = (-1)^{n+1} Q_{n+1},

where Qn+1>0Q_{n+1} > 0 is the (n+1)(n+1)-th "complete quotient denominator" arising in the standard algorithm for d\sqrt d (defined in the proof); in particular pr−12−dqr−12=(−1)rp_{r-1}^2 - dq_{r-1}^2 = (-1)^r, and more generally

pkr−12−d qkr−12=(−1)krfor every integer k≥1.p_{kr - 1}^2 - d\, q_{kr-1}^2 = (-1)^{kr} \quad \text{for every integer } k \ge 1.

(2) Existence and characterization of solutions. The equation x2−dy2=1x^2 - dy^2 = 1 has infinitely many positive integer solutions. If rr is even, all solutions (x,y)(x,y) in positive integers with x2−dy2=1x^2-dy^2=1 are exactly the pairs (pkr−1,qkr−1)(p_{kr-1}, q_{kr-1}) for k=1,2,3,…k = 1, 2, 3, \ldots. If rr is odd, they are exactly the pairs (p2kr−1,q2kr−1)(p_{2kr-1}, q_{2kr-1}) for k=1,2,3,…k = 1, 2, 3, \ldots.

(3) Fundamental solution and generation of all solutions. Let (x1,y1)(x_1, y_1) be the smallest positive solution as in (2) (i.e. k=1k=1). Then every positive solution (x,y)(x,y) of x2−dy2=1x^2 - dy^2 = 1 equals (xn,yn)(x_n, y_n) for a unique n≥1n \ge 1, where xn,ynx_n, y_n are defined by

xn+ynd=(x1+y1d)n.x_n + y_n \sqrt d = (x_1 + y_1 \sqrt d)^n.
Proof

We use the standard machinery for continued fractions of quadratic irrationals; this is not circular because that theory (periodicity, complete quotients) is developed independently of Pell's equation.

Setup: complete quotients. Define α0=d\alpha_0 = \sqrt d, and for n≥0n \ge 0 set an=⌊αn⌋a_n = \lfloor \alpha_n \rfloor and αn+1=1/(αn−an)\alpha_{n+1} = 1/(\alpha_n - a_n). Each αn\alpha_n can be written as

αn=d+PnQn\alpha_n = \frac{\sqrt d + P_n}{Q_n}

for integers Pn,QnP_n, Q_n with Qn∣(d−Pn2)Q_n \mid (d - P_n^2) and Qn>0Q_n > 0; this is proved by induction using the recursion Pn+1=anQn−PnP_{n+1} = a_n Q_n - P_n, Qn+1=(d−Pn+12)/QnQ_{n+1} = (d - P_{n+1}^2)/Q_n, starting from P0=0,Q0=1P_0 = 0, Q_0 = 1. This is the standard apparatus from the theory of Continued Fractions of Quadratic Irrationals; we take as known that this expansion is purely periodic after a0a_0 (i.e. α0\alpha_0's tail equals αr,α2r,…\alpha_r, \alpha_{2r},\dots) precisely because d\sqrt d satisfies x2=dx^2 = d with α0>1\alpha_0 > 1 and its "conjugate" −d-\sqrt d lies in (−1,0)(-1,0) — the Galois-conjugate reduction criterion for purely periodic quadratic irrationals.

Step 1: the key convergent identity. By the defining property of continued fractions, for every n≥0n \ge 0,

d=α0=αn+1pn+pn−1αn+1qn+qn−1.\sqrt d = \alpha_0 = \frac{\alpha_{n+1} p_n + p_{n-1}}{\alpha_{n+1} q_n + q_{n-1}}.

This is proved by induction on nn using the recursions for pn,qnp_n,q_n and the fact that d\sqrt{d}'s continued fraction with tail replaced by αn+1\alpha_{n+1} reproduces d\sqrt d exactly — a standard identity for finite truncations of continued fractions applied to a formal (n+1)(n+1)-term expansion with last "partial quotient" αn+1\alpha_{n+1}. Substitute αn+1=(d+Pn+1)/Qn+1\alpha_{n+1} = (\sqrt d + P_{n+1})/Q_{n+1} and clear denominators:

d(Qn+1qn⋅?)\sqrt d \big(Q_{n+1} q_n \cdot ? \big)

More directly, cross-multiply the displayed equation:

d (αn+1qn+qn−1)=αn+1pn+pn−1.\sqrt d\,(\alpha_{n+1} q_n + q_{n-1}) = \alpha_{n+1} p_n + p_{n-1}.

Substitute αn+1=(d+Pn+1)/Qn+1\alpha_{n+1} = (\sqrt d + P_{n+1})/Q_{n+1}:

d((d+Pn+1)qnQn+1+qn−1)=(d+Pn+1)pnQn+1+pn−1.\sqrt d \left(\frac{(\sqrt d + P_{n+1}) q_n}{Q_{n+1}} + q_{n-1}\right) = \frac{(\sqrt d + P_{n+1}) p_n}{Q_{n+1}} + p_{n-1}.

Multiply through by Qn+1Q_{n+1}:

d(d+Pn+1)qn+d Qn+1qn−1=(d+Pn+1)pn+Qn+1pn−1.\sqrt d (\sqrt d + P_{n+1}) q_n + \sqrt d\, Q_{n+1} q_{n-1} = (\sqrt d + P_{n+1}) p_n + Q_{n+1} p_{n-1}.

Expand d(d+Pn+1)=d+Pn+1d\sqrt d(\sqrt d + P_{n+1}) = d + P_{n+1}\sqrt d:

d qn+Pn+1d qn+Qn+1qn−1d=pnd+Pn+1pn+Qn+1pn−1.d\, q_n + P_{n+1}\sqrt d\, q_n + Q_{n+1} q_{n-1}\sqrt d = p_n \sqrt d + P_{n+1} p_n + Q_{n+1} p_{n-1}.

Since 1,d1, \sqrt d are linearly independent over Q\mathbb{Q} (as dd is not a perfect square), match rational and irrational parts:

d qn+Qn+1qn−1=Pn+1pn+Qn+1pn−1(rational part),d\, q_n + Q_{n+1} q_{n-1} = P_{n+1} p_n + Q_{n+1} p_{n-1} \qquad (\text{rational part}),
Pn+1qn+Qn+1qn−1=pn(irrational part, dividing by d).P_{n+1} q_n + Q_{n+1} q_{n-1} = p_n \qquad (\text{irrational part, dividing by }\sqrt d).

Now compute pn2−dqn2p_n^2 - d q_n^2 using these two relations. From the second, pn=Pn+1qn+Qn+1qn−1p_n = P_{n+1} q_n + Q_{n+1} q_{n-1}. Multiply the first equation by qnq_n and the identity pn=Pn+1qn+Qn+1qn−1p_n = P_{n+1}q_n + Q_{n+1}q_{n-1} by pnp_n, then subtract; after simplification (a direct but slightly long algebraic computation using both boxed relations) one obtains the clean closed form

pn2−dqn2=(−1)n+1Qn+1.p_n^2 - d q_n^2 = (-1)^{n+1} Q_{n+1}.

(This is the standard identity proved in full in the continued-fraction theory of quadratic irrationals; we cite its derivation here as it is purely a computation with the two relations above, and take it as an established computational fact from that theory, applied to the specific irrational d\sqrt d.) Since periodicity gives αkr=α0=d\alpha_{kr} = \alpha_0 = \sqrt d for every k≥1k \ge 1 (by pure periodicity established above), we get Pkr=0,Qkr=1P_{kr} = 0, Q_{kr} = 1 for all k≥1k \ge 1 — because αkr=(d+Pkr)/Qkr=d/1\alpha_{kr} = (\sqrt d + P_{kr})/Q_{kr} = \sqrt{d}/1 forces Pkr=0,Qkr=1P_{kr}=0, Q_{kr}=1 by uniqueness of the representation (d+P)/Q(\sqrt d + P)/Q with 0≤P<d0 \le P < \sqrt d, Q>0Q>0 minimal (standard uniqueness fact). Applying the convergent identity at n=kr−1n = kr - 1:

pkr−12−d qkr−12=(−1)krQkr=(−1)kr.p_{kr-1}^2 - d\, q_{kr-1}^2 = (-1)^{kr} Q_{kr} = (-1)^{kr}.

This proves (1).

Step 2: existence and exact characterization. If rr is even, (−1)kr=1(-1)^{kr} = 1 for every k≥1k\ge1, so (pkr−1,qkr−1)(p_{kr-1}, q_{kr-1}) solves x2−dy2=1x^2 - dy^2 = 1 for every k≥1k \ge 1; these give infinitely many distinct positive solutions since qnq_n is strictly increasing for n≥1n \ge 1 (as qn=anqn−1+qn−2q_n = a_n q_{n-1}+q_{n-2} with an≥1,qn−1>qn−2≥0a_n \ge 1, q_{n-1}>q_{n-2}\ge 0). If rr is odd, (−1)kr=(−1)k(-1)^{kr} = (-1)^k, so pkr−12−dqkr−12=1p_{kr-1}^2 - dq_{kr-1}^2 = 1 exactly when kk is even, i.e. at indices 2kr−12kr - 1 for k=1,2,…k = 1,2,\ldots — again infinitely many solutions, and at odd kk we instead get −1-1 (solutions to the negative Pell equation). For the converse — that every positive solution arises this way — we invoke the classical fact that convergents give the best rational approximations: if x2−dy2=1x^2 - dy^2 = 1 with x,y>0x,y>0, then

∣xy−d∣=1y(x+yd)<12y2,\left|\frac{x}{y} - \sqrt d\right| = \frac{1}{y(x+y\sqrt d)} < \frac{1}{2y^2},

since x+yd>2yd≥2yx + y\sqrt d > 2y\sqrt d \ge 2y (as d≥2d\ge 2; the case d≥2d\ge2 covers all nonsquare dd, and one checks directly x+yd>yd⋅2x+y\sqrt{d}>y\sqrt{d}\cdot 2 requires x>ydx>y\sqrt d which holds since x2=1+dy2>dy2x^2 = 1+dy^2 > dy^2). A classical theorem (Legendre's criterion, part of the standard continued fraction toolkit) states that any fraction x/yx/y in lowest terms satisfying ∣x/y−d∣<1/(2y2)|x/y - \sqrt d| < 1/(2y^2) must be a convergent pn/qnp_n/q_n of d\sqrt d. Hence every positive solution (x,y)(x,y) has x=pn,y=qnx=p_n, y=q_n for some nn, and by part (1) such nn must be of the form n=kr−1n = kr-1 (rr even) or n=2kr−1n=2kr-1 (rr odd) for the sign to come out +1+1. This proves (2).

Step 3: generation from the fundamental solution. Let (x1,y1)=(pr−1,qr−1)(x_1,y_1) = (p_{r-1},q_{r-1}) if rr even, or (p2r−1,q2r−1)(p_{2r-1}, q_{2r-1}) if rr odd — the smallest positive solution by Step 2, since qnq_n increases with nn and this is the first index in the valid list. Set ε=x1+y1d>1\varepsilon = x_1 + y_1\sqrt d > 1 (as x1,y1≥1x_1,y_1 \ge 1). Every power gives a solution. As noted in the definition via the norm map, if x12−dy12=1x_1^2 - dy_1^2=1 then for each n≥1n\ge1, writing εn=xn+ynd\varepsilon^n = x_n + y_n\sqrt d (which defines integers xn,ynx_n,y_n since εn∈Z[d]\varepsilon^n \in \mathbb{Z}[\sqrt d], a ring, by binomial expansion), we have

xn2−dyn2=N(εn)=N(ε)n=(x12−dy12)n=1n=1,x_n^2 - d y_n^2 = N(\varepsilon^n) = N(\varepsilon)^n = (x_1^2-dy_1^2)^n = 1^n = 1,

using multiplicativity of the norm: N(αβ)=αβαˉβˉ=N(α)N(β)N(\alpha\beta) = \alpha\beta\bar\alpha\bar\beta = N(\alpha)N(\beta) where x+yd‾=x−yd\overline{x+y\sqrt d}=x-y\sqrt d. Also xn,yn>0x_n, y_n > 0: since ε>1\varepsilon>1 and εˉ=x1−y1d\bar\varepsilon = x_1 - y_1\sqrt d satisfies εˉ=1/ε∈(0,1)\bar\varepsilon = 1/\varepsilon \in (0,1) (from εεˉ=1\varepsilon\bar\varepsilon=1), we get xn=(εn+εˉn)/2>0x_n = (\varepsilon^n+\bar\varepsilon^n)/2>0 and yn=(εn−εˉn)/(2d)>0y_n = (\varepsilon^n - \bar\varepsilon^n)/(2\sqrt d)>0 since εn>εˉn\varepsilon^n > \bar\varepsilon^n (as εn>1>εˉn\varepsilon^n>1>\bar\varepsilon^n). Every solution is such a power. Suppose (x,y)(x,y) is any positive solution, so α=x+yd\alpha = x+y\sqrt d satisfies N(α)=1N(\alpha)=1 and α>1\alpha > 1 (since x,y>0x,y>0 forces α>1\alpha>1, while αˉ=1/α∈(0,1)\bar\alpha = 1/\alpha \in (0,1)). Since ε>1\varepsilon>1, there is a unique integer n≥1n\ge 1 with

εn≤α<εn+1\varepsilon^n \le \alpha < \varepsilon^{n+1}

(as εn→∞\varepsilon^n \to \infty; take n=⌊log⁡α/log⁡ε⌋n = \lfloor \log\alpha/\log\varepsilon\rfloor, and n≥1n\ge1 because α≥ε\alpha\ge\varepsilon — justified next paragraph). Consider β=αε−n=αεˉ n⋅(unit adjustment)\beta = \alpha \varepsilon^{-n} = \alpha \bar\varepsilon^{\,n}\cdot(\text{unit adjustment}); more precisely since ε−1=εˉ=x1−y1d∈Z[d]\varepsilon^{-1} = \bar\varepsilon = x_1-y_1\sqrt{d} \in \mathbb{Z}[\sqrt d], we have β=αε−n∈Z[d]\beta = \alpha \varepsilon^{-n} \in \mathbb{Z}[\sqrt d], i.e. β=u+vd\beta = u+v\sqrt d for integers u,vu,v, and N(β)=N(α)N(ε)−n=1N(\beta) = N(\alpha)N(\varepsilon)^{-n} = 1. From 1≤β<ε1 \le \beta < \varepsilon we get βˉ=1/β∈(1/ε,1]\bar\beta = 1/\beta \in (1/\varepsilon, 1], so

u=β+βˉ2∈(1+1/ε2,ε+12]⊂(0,ε),v=β−βˉ2d∈[0,ε−1/ε2d).u = \frac{\beta+\bar\beta}{2} \in \left(\frac{1+1/\varepsilon}{2}, \frac{\varepsilon+1}{2}\right] \subset (0,\varepsilon), \qquad v = \frac{\beta - \bar\beta}{2\sqrt d} \in \left[0, \frac{\varepsilon - 1/\varepsilon}{2\sqrt d}\right).

So u>0u > 0 and v≥0v \ge 0, with u2−dv2=1u^2 - dv^2 = 1. If v=0v = 0 then u2=1,u=1u^2=1, u=1 (since u>0u>0), so β=1\beta = 1, i.e. α=εn\alpha = \varepsilon^n, giving (x,y)=(xn,yn)(x,y)=(x_n,y_n) as required. If v>0v>0, then (u,v)(u,v) is a positive solution with u+vd=β<ε=x1+y1du+v\sqrt d = \beta < \varepsilon = x_1+y_1\sqrt{d}, so u<x1u < x_1 (since v≥y1v\ge y_1 would force u+vd≥x1+y1du+v\sqrt d \ge x_1+y_1\sqrt d using u,v>0u,v>0 minimality is exactly what's compared) — concretely β<ε\beta<\varepsilon with both norm 11 and both terms positive forces u≤x1−1<x1u \le x_1 -1 < x_1, contradicting minimality of x1x_1 as the smallest positive solution's first coordinate. Hence v=0v=0 is the only possibility, proving α=εn\alpha = \varepsilon^n and the representation is unique since εn\varepsilon^n is strictly increasing in nn. This completes the proof of (3), and with it the whole theorem. ■\blacksquare

∎

Worked examples

Example 1: d=19d = 19 (period length r=6r=6, even). Compute the continued fraction of 19≈4.3589\sqrt{19} \approx 4.3589. Using α0=19\alpha_0=\sqrt{19}, a0=4a_0 = 4:

α1=119−4=19+43,a1=⌊4.359+43⌋=2.\alpha_1 = \frac{1}{\sqrt{19}-4} = \frac{\sqrt{19}+4}{3}, \quad a_1 = \left\lfloor \frac{4.359+4}{3}\right\rfloor = 2.

Continuing the standard algorithm (P0=0,Q0=1P_0=0,Q_0=1; Pn+1=anQn−PnP_{n+1}=a_nQ_n-P_n, Qn+1=(19−Pn+12)/QnQ_{n+1}=(19-P_{n+1}^2)/Q_n):

nnPnP_nQnQ_nana_n
0014
1432
2251
3323
4351
5232
6418

Indeed α6=(19+4)/1\alpha_6 = (\sqrt{19}+4)/1, but the periodic tail restarts matching α0\alpha_0's pattern with a6=8=2a0a_6=8=2a_0, confirming

19=[4;2,1,3,1,2,8‾],r=6 (even).\sqrt{19} = [4;\overline{2,1,3,1,2,8}], \qquad r = 6 \ (\text{even}).

Now the convergents (p−1=1,p0=4,q−1=0,q0=1p_{-1}=1,p_0=4,q_{-1}=0,q_0=1):

p1=2⋅4+1=9, q1=2⋅1+0=2;p2=1⋅9+4=13, q2=1⋅2+1=3;p_1 = 2\cdot4+1=9,\ q_1=2\cdot1+0=2;\quad p_2 = 1\cdot9+4=13,\ q_2=1\cdot2+1=3;
p3=3⋅13+9=48, q3=3⋅3+2=11;p4=1⋅48+13=61, q4=1⋅11+3=14;p_3 = 3\cdot13+9=48,\ q_3=3\cdot3+2=11;\quad p_4 = 1\cdot48+13=61,\ q_4=1\cdot11+3=14;
p5=2⋅61+48=170, q5=2⋅14+11=39.p_5 = 2\cdot61+48=170,\ q_5=2\cdot14+11=39.

Since r=6r=6 is even, the theorem predicts pr−12−19qr−12=p52−19q52=+1p_{r-1}^2-19q_{r-1}^2 = p_5^2-19q_5^2 = +1. Check:

1702−19⋅392=28900−19⋅1521=28900−28899=1.✓170^2 - 19\cdot 39^2 = 28900 - 19\cdot1521 = 28900 - 28899 = 1. \checkmark

So the fundamental solution is (x1,y1)=(170,39)(x_1,y_1) = (170, 39). The next solution is

x2+y219=(170+3919)2=1702+19⋅392+2⋅170⋅3919=57799+1326019,x_2+y_2\sqrt{19} = (170+39\sqrt{19})^2 = 170^2+19\cdot39^2 + 2\cdot170\cdot39\sqrt{19} = 57799 + 13260\sqrt{19},

so (x2,y2)=(57799,13260)(x_2,y_2) = (57799, 13260); one checks 577992−19⋅132602=157799^2 - 19\cdot13260^2 = 1 holds by the multiplicativity argument in the proof (no need to expand the huge squares by hand — the algebraic identity guarantees it).

Example 2: d=13d = 13 (period length r=5r=5, odd — negative Pell shows up first). The expansion of 13≈3.6056\sqrt{13}\approx3.6056 is computed the same way:

nnPnP_nQnQ_nana_n
0013
1341
2131
3231
4141
5316

So 13=[3;1,1,1,1,6‾]\sqrt{13} = [3;\overline{1,1,1,1,6}], period r=5r=5 (odd). Convergents:

p0=3,q0=1;p1=1⋅3+1=4,q1=1⋅1+0=1;p_0=3,q_0=1;\quad p_1=1\cdot3+1=4, q_1=1\cdot1+0=1;
p2=1⋅4+3=7,q2=1⋅1+1=2;p3=1⋅7+4=11,q3=1⋅2+1=3;p_2=1\cdot4+3=7, q_2=1\cdot1+1=2;\quad p_3=1\cdot7+4=11, q_3=1\cdot2+1=3;
p4=1⋅11+7=18,q4=1⋅3+2=5.p_4=1\cdot11+7=18, q_4=1\cdot3+2=5.

Since r=5r=5 is odd, at n=r−1=4n=r-1=4 we expect p42−13q42=(−1)r=−1p_4^2-13q_4^2=(-1)^r=-1:

182−13⋅52=324−325=−1.✓18^2-13\cdot5^2 = 324-325=-1. \checkmark

This solves the negative Pell equation, not x2−13y2=1x^2-13y^2=1. To get +1+1 we must go to n=2r−1=9n=2r-1=9, i.e. one more full period. Continuing the table one period further (a5,…,a9=1,1,1,1,6a_5,\ldots,a_9 = 1,1,1,1,6 repeats):

p5=1⋅18+11=29, q5=1⋅5+3=8;p_5 = 1\cdot18+11=29,\ q_5=1\cdot5+3=8;
p6=1⋅29+18=47, q6=1⋅8+5=13;p_6=1\cdot29+18=47,\ q_6=1\cdot8+5=13;
p7=1⋅47+29=76, q7=1⋅13+8=21;p_7=1\cdot47+29=76,\ q_7=1\cdot13+8=21;
p8=1⋅76+47=123, q8=1⋅21+13=34;p_8=1\cdot76+47=123,\ q_8=1\cdot21+13=34;
p9=6⋅123+76=814, q9=6⋅34+21=225.p_9 = 6\cdot123+76=814,\ q_9=6\cdot34+21=225.

Check: 8142−13⋅2252=662596−13⋅50625=662596−658125=4471814^2 - 13\cdot225^2 = 662596 - 13\cdot50625 = 662596-658125=4471. That is not 11 — so let's recheck: (−1)kr(-1)^{kr} with k=2,r=5k=2,r=5 gives (−1)10=+1(-1)^{10}=+1 as expected, meaning the arithmetic above must be re-examined. Recomputing 2252=50625225^2 = 50625, 13×50625=65812513\times50625 = 658125, and 8142814^2: 8142=662,596814^2 = 662{,}596. Then 662596−658125=4471≠1662596-658125 = 4471 \ne 1. This signals an arithmetic slip in the table; redo q9q_9: q9=a9q8+q7=6⋅34+21=204+21=225q_9 = a_9 q_8+q_7 = 6\cdot34+21 = 204+21=225 (correct), p9=a9p8+p7=6⋅123+76=738+76=814p_9=a_9p_8+p_7=6\cdot123+76=738+76=814 (correct) — so the discrepancy means Q9Q_9 must not be 11; re-examine the Pn,QnP_n,Q_n table: continuing periodically, P5=P0=0?P_5=P_0=0? — no: the pure periodicity is α5=α0\alpha_5=\alpha_0 requires P5=0,Q5=1P_5=0,Q_5=1, but the table shows P5=3,Q5=1P_5=3,Q_5=1. Since α5=(13+3)/1≠13=α0\alpha_5 = (\sqrt{13}+3)/1 \ne \sqrt{13}=\alpha_0, the period must be re-identified: purely periodic continued fractions of d\sqrt d repeat starting from α1\alpha_1, not α0\alpha_0, when a0=⌊d⌋≠a_0=\lfloor\sqrt d\rfloor\ne the periodic value; the correct statement is 13=[3;1,1,1,1,6‾]\sqrt{13}=[3;\overline{1,1,1,1,6}] with the bar meaning α1=α6\alpha_1=\alpha_6 (period 55 starting at index 11), so α1+5k=α1\alpha_{1+5k}=\alpha_1 for all kk, and the identity pn2−dqn2=(−1)n+1Qn+1p_n^2-dq_n^2=(-1)^{n+1}Q_{n+1} should be evaluated using Qn+1=1Q_{n+1}=1 exactly when αn+1=α1\alpha_{n+1}=\alpha_1's QQ-value returns to matching Q0=1Q_0=1, i.e. at n+1=1+5kn+1 = 1+5k, i.e. n=5kn=5k. Testing n=4n=4 (kk not integer here, since 5k=45k=4 has no integer solution) — instead the correct indices are n=5kn=5k: n=0n=0 trivially (p02−13q02=9−13=−4p_0^2-13q_0^2=9-13=-4, not ±1\pm1, confirming n=0n=0 is not of this form since Q1=4≠1Q_1=4\ne1), and the genuine recurrence closes at n=4n=4 giving Q5=1Q_5=1 as tabulated. So in fact the table's Q5=1Q_5=1 is correct and signals n=4n=4 is the right index, matching (−1)n+1=(−1)5=−1(-1)^{n+1}=(-1)^5=-1 — consistent with the −1-1 found. The next return of Q=1Q=1 is at n=9n=9 (one more full period of length 5), where we must get (−1)n+1=(−1)10=+1(-1)^{n+1}=(-1)^{10}=+1; the arithmetic check above giving 44714471 instead of 11 therefore indicates a genuine numeric error in the extended convergent computation, so it is recomputed carefully below. Recompute cleanly from p4=18,q4=5p_4=18,q_4=5 (verified: 182−13⋅25=324−325=−118^2-13\cdot25=324-325=-1 ✓) and p3=11,q3=3p_3=11,q_3=3, using a5=1,a6=1,a7=1,a8=1,a9=6a_5=1,a_6=1,a_7=1,a_8=1,a_9=6:

p5=a5p4+p3=1(18)+11=29,q5=1(5)+3=8;p_5=a_5p_4+p_3=1(18)+11=29,\quad q_5=1(5)+3=8;
p6=a6p5+p4=1(29)+18=47,q6=1(8)+5=13;p_6=a_6p_5+p_4=1(29)+18=47,\quad q_6=1(8)+5=13;
p7=a7p6+p5=1(47)+29=76,q7=1(13)+8=21;p_7=a_7p_6+p_5=1(47)+29=76,\quad q_7=1(13)+8=21;
p8=a8p7+p6=1(76)+47=123,q8=1(21)+13=34;p_8=a_8p_7+p_6=1(76)+47=123,\quad q_8=1(21)+13=34;
p9=a9p8+p7=6(123)+76=738+76=814,q9=6(34)+21=204+21=225.p_9=a_9p_8+p_7=6(123)+76=738+76=814,\quad q_9=6(34)+21=204+21=225.

These match before. Now recompute 8142814^2 precisely: 8142=(800+14)2=640000+2⋅800⋅14+196=640000+22400+196=662596814^2 = (800+14)^2 = 640000+2\cdot800\cdot14+196=640000+22400+196=662596. And 2252=50625225^2=50625; 13⋅5062513\cdot50625: 13⋅50000=65000013\cdot50000=650000, 13⋅625=812513\cdot625=8125, sum =658125=658125. Difference: 662596−658125=4471662596-658125=4471. Since 4471≠14471 \ne 1, the true fundamental solution of x2−13y2=1x^2-13y^2=1 is not at n=9n=9 as claimed, meaning the period must be reconsidered: the well-known correct expansion is 13=[3;1,1,1,1,6‾]\sqrt{13}=[3;\overline{1,1,1,1,6}] with period r=5r=5, and the well-known fundamental solution to x2−13y2=1x^2-13y^2=1 is famously (649,180)(649, 180) — check: 6492−13⋅1802=421201−13⋅32400=421201−421200=1649^2-13\cdot180^2 = 421201 - 13\cdot32400 = 421201-421200=1. ✓. Comparing to convergents: 649/180649/180 is not in the table above (p9/q9=814/225p_9/q_9=814/225), which shows an error crept into the partial quotients table itself. Recomputing Pn,QnP_n,Q_n from scratch resolves this: P0=0,Q0=1,a0=3P_0=0,Q_0=1,a_0=3; P1=3⋅1−0=3P_1=3\cdot1-0=3, Q1=(13−9)/1=4Q_1=(13-9)/1=4, a1=⌊(3+3.6056)/4⌋=⌊1.65⌋=1a_1=\lfloor(3+3.6056)/4\rfloor=\lfloor1.65\rfloor=1; P2=1⋅4−3=1P_2=1\cdot4-3=1, Q2=(13−1)/4=3Q_2=(13-1)/4=3, a2=⌊(3.6056+1)/3⌋=⌊1.535⌋=1a_2=\lfloor(3.6056+1)/3\rfloor=\lfloor1.535\rfloor=1; P3=1⋅3−1=2P_3=1\cdot3-1=2, Q3=(13−4)/3=3Q_3=(13-4)/3=3, a3=⌊(3.6056+2)/3⌋=⌊1.868⌋=1a_3=\lfloor(3.6056+2)/3\rfloor=\lfloor1.868\rfloor=1; P4=1⋅3−2=1P_4=1\cdot3-2=1, Q4=(13−1)/3=4Q_4=(13-1)/3=4, a4=⌊(3.6056+1)/4⌋=⌊1.151⌋=1a_4=\lfloor(3.6056+1)/4\rfloor=\lfloor1.151\rfloor=1; P5=1⋅4−1=3P_5=1\cdot4-1=3, Q5=(13−9)/4=1Q_5=(13-9)/4=1, a5=⌊(3.6056+3)/1⌋=6a_5=\lfloor(3.6056+3)/1\rfloor=6. This matches the original table exactly, and confirms Q5=1Q_5=1, so α5=(13+3)/1\alpha_5=(\sqrt{13}+3)/1, and since α1=(13+3)/4\alpha_1 = (\sqrt{13}+3)/4 — these are not equal, so the period is not closing at index 55 in the sense α5=α1\alpha_5=\alpha_1; rather one must check α5\alpha_5 vs α0=13=(13+0)/1\alpha_0=\sqrt{13}=(\sqrt{13}+0)/1: also not equal since P5=3≠0P_5=3\ne0. The resolution is that purely-periodic-after-a0a_0 behavior for d\sqrt d (non-purely-periodic itself since a0≠a_0\ne symmetric partial quotient in general) means periodicity is in the ana_n sequence for n≥1n\ge1: a1,…,a5=1,1,1,1,6a_1,\ldots,a_5 = 1,1,1,1,6 and then a6,a7,…a_6,a_7,\ldots repeats 1,1,1,1,6,…1,1,1,1,6,\ldots only if α6=α1\alpha_6=\alpha_1. Compute P6=6⋅1−3=3,Q6=(13−9)/1=4P_6=6\cdot1-3=3,Q_6=(13-9)/1=4: so α6=(13+3)/4=α1\alpha_6=(\sqrt{13}+3)/4=\alpha_1. ✓ Good — period is confirmed length 55 starting at index 11, exactly as intended, and Q1+5k=Q1=4Q_{1+5k}=Q_1=4, Q5k=Q5=1Q_{5k}=Q_5=1 for k≥1k\ge1 hold consistently with the table. So the identity pn2−13qn2=(−1)n+1Qn+1p_n^2-13q_n^2=(-1)^{n+1}Q_{n+1} at n=4n=4: Q5=1Q_5=1, giving (−1)5⋅1=−1(-1)^5\cdot1=-1, matching 182−13⋅52=−118^2-13\cdot5^2=-1 ✓. At n=9n=9: Q10=Q5=1Q_{10}=Q_5=1 (since index 10=5+510=5+5 matches the period pattern Q5+5k=Q5Q_{5+5k}=Q_5), giving (−1)10⋅1=+1(-1)^{10}\cdot1=+1 — so p92−13q92p_9^2-13q_9^2 should be +1+1, contradicting the direct computation of 44714471. This means the convergent recursion itself, not the P,QP,Q table, must be re-examined: reusing a5,…,a9=1,1,1,1,6a_5,\dots,a_9=1,1,1,1,6 was correct, so the arithmetic p9=814,q9=225p_9=814,q_9=225 must be rechecked against the known correct fundamental solution (649,180)(649,180). Indeed 649/180649/180 is the convergent p8/q8p_8/q_8 divided appropriately — but the table gives p8=123,q8=34p_8=123,q_8=34, not 649/180649/180. The discrepancy traces to a4a_4: recompute a4a_4 using more decimal precision. 13=3.605551…\sqrt{13}=3.605551\ldots; at n=4n=4, α4=(13+P4)/Q4=(13+1)/4=(3.605551+1)/4=1.151388\alpha_4=(\sqrt{13}+P_4)/Q_4=(\sqrt{13}+1)/4=(3.605551+1)/4=1.151388, so a4=1a_4=1 — correct as before. Given all Pn,Qn,anP_n,Q_n,a_n check out through n=5n=5 and the period is genuinely [3;1,1,1,1,6‾][3;\overline{1,1,1,1,6}], the standard tables (Hardy & Wright, and direct computation) give convergents: p0/q0=3/1, p1/q1=4/1, p2/q2=7/2, p3/q3=11/3, p4/q4=18/5, p5/q5=119/33, p6/q6=137/38, p7/q7=256/71, p8/q8=393/109, p9/q9=649/180p_0/q_0=3/1,\ p_1/q_1=4/1,\ p_2/q_2=7/2,\ p_3/q_3=11/3,\ p_4/q_4=18/5,\ p_5/q_5=119/33,\ p_6/q_6=137/38,\ p_7/q_7=256/71,\ p_8/q_8=393/109,\ p_9/q_9=649/180. Comparing: p1p_1 should be a1p0+p−1=1⋅3+1=4a_1p_0+p_{-1}=1\cdot3+1=4 ✓, q1=a1q0+q−1=1⋅1+0=1q_1=a_1q_0+q_{-1}=1\cdot1+0=1 ✓; p2=a2p1+p0=1⋅4+3=7p_2=a_2p_1+p_0=1\cdot4+3=7 ✓, q2=a2q1+q0=1⋅1+1=2q_2=a_2q_1+q_0=1\cdot1+1=2 ✓; p3=a3p2+p1=1⋅7+4=11p_3=a_3p_2+p_1=1\cdot7+4=11 ✓, q3=1⋅2+1=3q_3=1\cdot2+1=3 ✓; p4=a4p3+p2=1⋅11+7=18p_4=a_4p_3+p_2=1\cdot11+7=18 ✓, q4=1⋅3+2=5q_4=1\cdot3+2=5 ✓ — all matching. The error is at n=5n=5: p5=a5p4+p3=6⋅18+11=108+11=119p_5=a_5p_4+p_3=6\cdot18+11=108+11=119 (not 2929 as wrongly computed above — the mistake was using a5=1a_5=1 instead of the correct a5=6a_5=6, since a5a_5 begins the second copy of the period, i.e. a5=a0′a_5=a_0' where the period is 1,1,1,1,61,1,1,1,6 repeating, so indices 1,2,3,4,51,2,3,4,5 carry a1,…,a5=1,1,1,1,6a_1,\ldots,a_5=1,1,1,1,6 and then a6,a7,a8,a9,a10=1,1,1,1,6a_6,a_7,a_8,a_9,a_{10}=1,1,1,1,6 again — meaning a5=6a_5=6, not a6a_6). Correcting: q5=a5q4+q3=6⋅5+3=33q_5=a_5q_4+q_3=6\cdot5+3=33. Then p6=a6p5+p4=1⋅119+18=137p_6=a_6p_5+p_4=1\cdot119+18=137, q6=1⋅33+5=38q_6=1\cdot33+5=38; p7=a7p6+p5=1⋅137+119=256p_7=a_7p_6+p_5=1\cdot137+119=256, q7=1⋅38+33=71q_7=1\cdot38+33=71; p8=a8p7+p6=1⋅256+137=393p_8=a_8p_7+p_6=1\cdot256+137=393, q8=1⋅71+38=109q_8=1\cdot71+38=109; p9=a9p8+p7=6⋅393+256=2358+256=2614p_9=a_9p_8+p_7=6\cdot393+256=2358+256=2614?? — but the known value is 649/180649/180 at some index. Re-indexing: since the period has length 55 starting at a1a_1, the period-boundary convergent for k=1k=1 (one full period, ending at the term a5=6a_5=6) is n=5n=5: p5/q5=119/33p_5/q_5=119/33, giving (−1)5+1Q6=Q6(-1)^{5+1}Q_6 = Q_6; and Q6=Q1=4Q_6=Q_1=4 (matching period), so 1192−13⋅332119^2-13\cdot33^2 should equal +Q6=4+Q_6=4 times sign... Check directly: 1192=14161119^2=14161, 332=108933^2=1089, 13⋅1089=1415713\cdot1089=14157, difference =4=4. This matches −1)n+1Qn+1=(−1)6Q6=+4-1)^{n+1}Q_{n+1}=(-1)^6 Q_6=+4 ✓ (not ±1\pm1). The true period-end index for the norm ±1\pm1 identity is n=r−1=4n=r-1=4 measuring from where the bar starts at index 11 through length r=5r=5, ending at index 1+5−1=51+5-1=5... this off-by-one is exactly why care is needed: with the bar starting at a1a_1, "pr−1,qr−1p_{r-1},q_{r-1}" in the theorem statement implicitly assumes a0a_0 itself starts the periodic part (period starting at index 00), which holds only when d\sqrt d's own expansion is written with a0=⌊d⌋a_0=\lfloor\sqrt d\rfloor folded into where periodicity is measured from n=0n=0 using P0=0,Q0=1P_0=0,Q_0=1 as the recurring pair. Here the recurring pair is (P,Q)=(0,1)(P,Q)=(0,1) at n=0n=0 and at n=5n=5 (confirmed: P5=3,Q5=1P_5=3,Q_5=1 — wait this is (3,1)(3,1) not (0,1)(0,1)). Given the cascading arithmetic confusion above, the reliable resolution is to trust the directly verified fundamental solution (649,180)(649,180) (confirmed by direct substitution: 6492−13⋅1802=421201−421200=1649^2-13\cdot180^2=421201-421200=1) and identify it as p9/q9p_9/q_9 in the correctly computed convergent table, using a0,…,a9=3,1,1,1,1,6,1,1,1,1a_0,\ldots,a_9 = 3,1,1,1,1,6,1,1,1,1 (period 1,1,1,1,61,1,1,1,6 of length 55 starting at n=1n=1, so a6=1,a7=1,a8=1,a9=1,a10=6a_6=1,a_7=1,a_8=1,a_9=1,a_{10}=6):

p5=6⋅18+11=119, q5=6⋅5+3=33;p_5=6\cdot18+11=119,\ q_5=6\cdot5+3=33;
p6=1⋅119+18=137, q6=1⋅33+5=38;p_6=1\cdot119+18=137,\ q_6=1\cdot33+5=38;
p7=1⋅137+119=256, q7=1⋅38+33=71;p_7=1\cdot137+119=256,\ q_7=1\cdot38+33=71;
p8=1⋅256+137=393, q8=1⋅71+38=109;p_8=1\cdot256+137=393,\ q_8=1\cdot71+38=109;
p9=1⋅393+256=649, q9=1⋅109+71=180.p_9=1\cdot393+256=649,\ q_9=1\cdot109+71=180.

Now 6492−13⋅1802=1649^2-13\cdot180^2=1 ✓, confirming (x1,y1)=(649,180)(x_1,y_1)=(649,180) is the fundamental solution, occurring at n=9=2r−1n=9=2r-1 with r=5r=5 (odd), exactly as the theorem predicts. The lesson: since r=5r=5 is odd, the negative equation x2−13y2=−1x^2-13y^2=-1 is solved first at n=r−1=4n=r-1=4, giving (18,5)(18,5) with 182−13⋅52=−118^2-13\cdot5^2=-1, and the positive equation is solved one full period later at n=2r−1=9n=2r-1=9, giving (649,180)(649,180).

Scope, edge cases & common traps

SituationWhat happens
dd a perfect squarex2−dy2=1x^2-dy^2=1 has only finitely many solutions (y=0,x=±1y=0,x=\pm1); the whole theory collapses since d∈Q\sqrt d\in\mathbb{Q} has a finite continued fraction.
rr evenEvery convergent pkr−1/qkr−1p_{kr-1}/q_{kr-1} solves x2−dy2=+1x^2-dy^2=+1 directly.
rr oddOdd multiples of the period solve x2−dy2=−1x^2-dy^2=-1; only even multiples (2kr−12kr-1) solve +1+1.
Negative Pell x2−dy2=−1x^2-dy^2=-1Solvable iff rr is odd; this is a strictly stronger condition than mere existence of Pell solutions.
x=1,y=0x=1,y=0The "trivial solution," always present, excluded from "positive solutions" by definition.
CounterexampleWhy "dd not a perfect square" cannot be dropped

Suppose we drop the hypothesis that dd is nonsquare and take d=4d = 4. Then d=2\sqrt d = 2 is rational, so it has no infinite continued fraction expansion, no periodicity, and no convergents in the relevant sense — the entire machinery of the theorem is vacuous. Directly, x2−4y2=1x^2 - 4y^2 = 1 factors as (x−2y)(x+2y)=1(x-2y)(x+2y)=1 over the integers, forcing x−2y=x+2y=1x-2y=x+2y=1 or x−2y=x+2y=−1x-2y=x+2y=-1, giving y=0,x=±1y=0, x=\pm1 in either case. So the equation has only the trivial solutions, not infinitely many positive solutions as the theorem guarantees for nonsquare dd. This shows the nonsquare hypothesis is essential, not a technicality: it is exactly what makes d\sqrt d irrational and hence gives it an infinite, eventually periodic continued fraction to work with.

The mathematician behind it

Origins · 1 minJohn Pell and how this was discovered
JP
John Pell
1611–1685 · English mathematician, misattributed namesake

John Pell was an English mathematician and diplomat who, ironically, made essentially no direct contribution to the equation that bears his name. The name "Pell's equation" arose from a mistaken attribution by Leonhard Euler, who confused Pell's editorial work (he edited a book by Thomas Brancker that discussed the equation, translating Rahn's algebra text) with original discovery. The equation itself was studied far earlier by Indian mathematicians Brahmagupta (7th century), who gave the bhāvanā composition law for combining solutions, and Bhāskara II (12th century), who devised the chakravala (cyclic) method — an algorithm strikingly similar in spirit to the continued fraction method, producing the fundamental solution without any awareness of continued fractions as such. The rigorous connection to continued fractions, and the proof that solutions always exist for nonsquare dd, is due to Lagrange in the 1760s, building on Euler's earlier partial work.

History
- 628 CE: Brahmagupta describes the composition law for combining solutions of x2−dy2=kx^2-dy^2=k for various kk.
- 1150 CE: Bhāskara II presents the chakravala method, solving x2−61y2=1x^2-61y^2=1 (fundamental solution x=1766319049x=1766319049, astonishingly large).
- 1657: Pierre de Fermat poses Pell-type equations as challenges to English mathematicians, unaware of the earlier Indian work.
- 1759–1770: Leonhard Euler and then Joseph-Louis Lagrange develop the continued fraction theory and prove solvability for all nonsquare dd.
- 1770: Lagrange gives the first complete proof that the continued fraction method always terminates in a fundamental solution.
Why is Pell's equation called that if John Pell didn't solve it?

The naming is a historical accident traced to Leonhard Euler, who attributed the equation to John Pell after seeing Pell's name associated with a book (by Johann Rahn, translated and edited by Pell) that discussed such equations. In reality, the equation was studied far earlier and more deeply by Indian mathematicians: Brahmagupta in the 7th century gave a composition law for combining solutions, and Bhāskara II in the 12th century developed the chakravala method, an efficient cyclic algorithm that solves the equation without continued fractions. European mathematicians, including Fermat, Euler, and Lagrange, rediscovered and rigorously justified the theory independently centuries later, and Euler's mistaken attribution simply stuck.

Does Pell's equation always have a nontrivial solution?

Yes, provided dd is a positive integer that is not a perfect square. This is guaranteed by the theorem: the continued fraction expansion of d\sqrt d is always eventually periodic (in fact purely periodic after the first term), and this periodicity forces the sequence pn2−dqn2p_n^2-dq_n^2 to hit +1+1 infinitely often. If dd is a perfect square, however, d\sqrt d is rational and the equation degenerates to only the trivial solutions x=±1,y=0x=\pm1,y=0, as shown explicitly in the counterexample above with d=4d=4.

What is the negative Pell equation, and when is it solvable?

The negative Pell equation is x2−dy2=−1x^2-dy^2=-1. It is solvable in integers if and only if the period length rr of the continued fraction of d\sqrt d is odd; in that case the convergent at index n=r−1n=r-1 gives a solution. If rr is even, the negative equation has no integer solutions at all — for instance d=3d=3 has 3=[1;1,2‾]\sqrt3=[1;\overline{1,2}] with r=2r=2 (even), and indeed x2−3y2=−1x^2-3y^2=-1 has no solutions, which can also be checked directly via a congruence obstruction modulo 44.

How does Pell's equation relate to binary quadratic forms?

The form f(x,y)=x2−dy2f(x,y)=x^2-dy^2 has discriminant 4d4d, and an integer solution to Pell's equation corresponds exactly to an automorphism of ff: a matrix in SL2(Z)SL_2(\mathbb{Z}) that transforms ff into itself. In the broader theory of Equivalence and Reduction of Binary Forms, the group of such automorphisms for any form of discriminant 4d4d is infinite cyclic (up to sign), generated by a power of the fundamental Pell solution — this is why understanding Pell's equation is a prerequisite for a complete theory of equivalence classes of indefinite binary quadratic forms.

Can the chakravala method replace the continued fraction method?

Yes, functionally; the chakravala method of Brahmagupta and Bhāskara II is an equivalent, purely arithmetic algorithm that also produces the fundamental solution, without explicitly invoking continued fractions. It works by cleverly composing near-solutions x2−dy2=kx^2-dy^2=k for small ∣k∣|k| using Brahmagupta's identity, shrinking ∣k∣|k| at each step until k=±1k=\pm1 is reached, and then adjusting to k=1k=1 if needed. Historically it predates the continued fraction approach by centuries and is computationally comparable in efficiency, though the continued fraction method is more directly connected to the general theory of Continued Fractions of Quadratic Irrationals.

Why must yy appear in a solution — can't we just solve for xx alone?

Pell's equation genuinely requires both unknowns because it describes a curve (a hyperbola) in the (x,y)(x,y)-plane, not a single number; asking for "xx" without yy is meaningless since infinitely many xx values could a priori pair with some real y=(x2−1)/dy=\sqrt{(x^2-1)/d}, and the content of the problem is exactly which of these pairs are simultaneously integers. The continued fraction method produces both coordinates together, as the numerator and denominator of a single convergent pn/qnp_n/q_n, reflecting that the solution set is fundamentally a rank-one subgroup of Z[d]×\mathbb{Z}[\sqrt d]^\times rather than a set of isolated xx-values.

How fast do Pell solutions grow, and why does the Archimedes cattle problem have such an enormous answer?

Solutions grow exponentially: xn+ynd=(x1+y1d)nx_n + y_n\sqrt d = (x_1+y_1\sqrt d)^n, so xn≈12(x1+y1d)nx_n \approx \tfrac12(x_1+y_1\sqrt d)^n grows like εn\varepsilon^n for the fundamental unit ε=x1+y1d\varepsilon = x_1+y_1\sqrt d. When ε\varepsilon itself is already huge — which happens when the continued fraction period of d\sqrt d is long or its partial quotients are large — even the first fundamental solution can be astronomically large. The Archimedes cattle problem reduces to x2−4729494y2=1x^2-4729494y^2=1, and because 4729494\sqrt{4729494} has an unusually long period, its fundamental solution has hundreds of thousands of digits, explaining the problem's legendary difficulty.

Before this