Pell Equation

Prerequisites: continued fractions, quadratic fields

Introduction

This article discusses solving the (generalized) Pell equation. A generalized Pell equation is an indeterminate equation in $x$ and $y$ of the form

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

where $D$ is a positive integer that is not a perfect square¹, and $N$ is a nonzero integer. The narrow sense of the Pell equation refers specifically to the cases $N=1$ or $N=\pm 1$, and sometimes also includes $N=\pm 4$. Generalized Pell equations are closely related to finding quadratic integers of norm $N$ in the ring of real quadratic integers; the cases commonly called the (narrow) Pell equation can be viewed as finding units in that ring.

In this article, Pell equation specifically means the case $N=1$. Correspondingly, the case $N=-1$ is called the negative Pell equation² (negative Pell's equation).

Structure of the solutions

The integer solutions $(x,y)$ of a generalized Pell equation are closely connected with the quadratic integer $x+y\sqrt{D}$, so solutions of Pell equations are often written in the form $x+y\sqrt{D}$ in the literature. The norm of a quadratic integer is

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

so a generalized Pell equation is roughly the problem of finding quadratic integers of norm $N$. There is, however, a subtle difference. When $x$ and $y$ are both integers, $x+y\sqrt{D}$ is certainly a quadratic integer; conversely, a quadratic integer does not necessarily require $x$ and $y$ to be integers—in the case $D\equiv 1\pmod 4$, $x$ and $y$ may also both be half-integers³.

This distinction is especially important when finding fundamental units. A unit in a quadratic integer ring is a quadratic integer of norm $\pm 1$. For $D\equiv 2,3\pmod 4$, finding such units only requires solving the generalized Pell equation for $N=\pm 1$; for $D\equiv 1\pmod 4$, the cases $N=\pm 4$ must also be considered. The method for finding units will be discussed below.

To understand the structure of generalized Pell equation solutions, begin with the Brahmagupta identity:

$$(x_1^2-Dy_1^2)(x_2^2-Dy_2^2)=(x_1x_2+Dy_1y_2)^2-D(x_1y_2+x_2y_1)^2.$$

This is equivalent to the multiplicativity of the norm of a quadratic integer:

$$\begin{aligned} N\left(x_1+y_1\sqrt{D}\right)N\left(x_2+y_2\sqrt{D}\right) &= N\left((x_1+y_1\sqrt{D})(x_2+y_2\sqrt{D})\right) \\ &= N\left((x_1x_2+Dy_1y_2)+(x_1y_2+x_2y_1)\sqrt{D}\right). \end{aligned}$$

Using this identity, integer solutions of $x^2-Dy^2=N_1$ and $x^2-Dy^2=N_2$ can be combined to obtain an integer solution of $x^2-Dy^2=N_1N_2$. From the viewpoint of quadratic integers, combining solutions is simply multiplication of quadratic integers, which illustrates why writing Pell equation solutions in quadratic-integer form is convenient. In particular, taking $N_1=N$ and $N_2=1$ shows that if one solution of $x^2-Dy^2=N$ and all solutions of the corresponding Pell equation $x^2-Dy^2=1$ are known, more solutions of $x^2-Dy^2=N$ can be obtained. This method need not generate all solutions, but it at least shows that understanding the structure of Pell equation solutions is important for understanding the structure of generalized Pell equation solutions.

Pell equation

Geometrically, the equation $x^2-Dy^2=1$ is a hyperbola whose real axis is the $x$-axis and imaginary axis is the $y$-axis. Every point on the hyperbola corresponds uniquely to a nonzero value of $x+y\sqrt{D}$: the left branch corresponds to negative values and the right branch to positive values. Moreover, on each branch, the values of $x+y\sqrt{D}$ increase strictly from bottom to top. The values of quadratic integers therefore give the solutions of the Pell equation a natural order in $x+y\sqrt{D}$.

The hyperbola is symmetric about both the $x$-axis and the $y$-axis, so it is enough to consider the part in the first quadrant; the other solutions follow by symmetry. Equivalently, consider only solutions with $x+y\sqrt{D}>1$. If the equation has a nontrivial solution besides $(\pm 1,0)$, then there is a solution $(x_1,y_1)$ in the first quadrant for which $x+y\sqrt{D}$ is minimal. It is also the lattice point in the first quadrant, excluding the coordinate axes, with both coordinates minimal, and is called the fundamental solution of the Pell equation⁴. By the discussion above, every integer pair $(x_k,y_k)$ satisfying $x_k+y_k\sqrt{D}=(x_1+y_1\sqrt{D})^k$ is a Pell equation solution in the first quadrant. Conversely, these are exactly all Pell equation solutions in the first quadrant. Symmetry then gives the following conclusion:

Theorem

Let $(x_1,y_1)$ be the fundamental solution of $x^2-Dy^2=1$. Then all its solutions are

$$\{(x,y):x+y\sqrt{D}=\pm(x_1+y_1\sqrt{D})^k,k\in\mathbf Z\}.$$

Proof

First prove that there are no other solutions in the first quadrant. Suppose there is another solution $x+y\sqrt{D}$ and, for some $k\ge 0$,

$$x_k+y_k\sqrt{D}< x+y\sqrt{D}< x_{k+1}+y_{k+1}\sqrt{D}.$$

Geometrically, this means that the lattice point $(x,y)$ lies strictly between $(x_k,y_k)$ and $(x_{k+1},y_{k+1})$ on the hyperbola. Multiplying the inequality by $x_k-y_k\sqrt{D}=(x_k+y_k\sqrt{D})^{-1}$ gives

$$1< (x+y\sqrt{D})(x_k-y_k\sqrt{D})=(xx_k-Dyy_k)+(x_ky-xy_k)\sqrt{D} < x_1+y_1\sqrt{D}.$$

By the monotonicity mentioned above, this says that $(xx_k-Dyy_k,x_ky-xy_k)$ is an integer point between $(1,0)$ and $(x_1,y_1)$, contradicting the choice of $(x_1,y_1)$.

When extending the first-quadrant solutions to the whole plane, replacing $k$ by its opposite (taking the entire expression's inverse) gives reflection about the $x$-axis, while negating the entire expression gives reflection about the origin. Adding the trivial solution for $k=0$ gives all Pell equation solutions.

The preceding discussion only assumed that a fundamental solution exists. We now show that the Pell equation always has a nontrivial solution.

Theorem

The Pell equation $x^2-Dy^2=1$ always has an integer solution other than $(\pm 1,0)$.

Proof

First, Dirichlet's theorem states that there are infinitely many pairs of positive integers $(x,y)$ satisfying

$$\left|\dfrac{x}{y}-\sqrt{D}\right| \le \dfrac{1}{y^2}$$

and all of them satisfy

$$|x^2-Dy^2|=y^2\left|\dfrac{x}{y}-\sqrt{D}\right|\left|\dfrac{x}{y}+\sqrt{D}\right| \le \dfrac{1}{y^2}+2\sqrt{D}<1+2\sqrt{D}.$$

Hence there is an integer $m\in(-1-2\sqrt{D},1+2\sqrt{D})$ such that infinitely many pairs $(x,y)$ satisfy $x^2-Dy^2 = m$. Classify these $(x,y)$ by their residues modulo $m$. For some pair of integers $(x_0,y_0)$, infinitely many $(x,y)$ must satisfy $x\equiv x_0\pmod m$ and $y\equiv y_0\pmod m$. Choose two distinct pairs $(x_1,y_1)$ and $(x_2,y_2)$ satisfying these conditions. Then

$$\dfrac{x_1+y_1\sqrt{D}}{x_2+y_2\sqrt{D}}=\dfrac{x_1x_2-Dy_1y_2}{m}+\dfrac{x_2y_1-x_1y_2}{m}\sqrt{D}.$$

By the congruences,

$$\begin{aligned} x_1x_2-Dy_1y_2 &\equiv x_0^2-Dy_0^2 = m \equiv 0 \pmod{|m|},\\ x_2y_1-x_1y_2 &\equiv x_0y_0-x_0y_0 = 0 \pmod{|m|}, \end{aligned}$$

so the right-hand side is an integer quadratic integer. Since $(x_1,y_1)\neq(x_2,y_2)$, it is nontrivial. Thus the Pell equation indeed has a nontrivial solution.

The proof in this section is nonconstructive. In the discussion of solving Pell equations below, we directly use convergents of continued fractions to construct a Pell equation solution, providing another proof that a nontrivial solution exists. Although the structure of the Pell equation solutions obtained here agrees with the structure of units in the ring of real quadratic integers, for $D\equiv 1\pmod 4$ this section does not completely resolve the structure of units in the corresponding ring; that will be discussed below.

Generalized Pell equation

The graph of the generalized Pell equation $x^2-Dy^2=N$ is also a hyperbola in the plane, symmetric about both the $x$-axis and the $y$-axis. As noted above, some solutions of the equation $x^2-Dy^2=N$ may differ only by a factor that is a Pell equation solution, so the solutions of $x^2-Dy^2=N$ can be divided into equivalence classes. For two solutions $(x_1,y_1)$ and $(x_2,y_2)$ of $x^2-Dy^2=N$, if there is a Pell equation solution $(u,v)$ such that $x_2+y_2\sqrt{D}=(x_1+y_1\sqrt{D})(u+v\sqrt{D})$, then the solutions $(x_1,y_1)$ and $(x_2,y_2)$ are called equivalent. The necessary and sufficient conditions for two solutions to be equivalent are

$$N\mid (x_1x_2-Dy_1y_2),\ N\mid (x_2y_1-x_1y_2).$$

Since Pell equation solutions are relatively easy to find, a natural idea is to find one solution in each equivalence class. Once these solutions are known, all solutions of the generalized Pell equation can be obtained using the corresponding Pell equation solutions. Within the equivalence classes of generalized Pell equation solutions, symmetry guarantees a solution with nonnegative and as small as possible $y$: if it is unique, it is called the fundamental solution of that class; otherwise, the class has two solutions with nonnegative minimal $y$, symmetric about the $y$-axis, and the one with $x>0$ is chosen as the fundamental solution. Thus solving $x^2-Dy^2=N$ amounts to finding its set $U$ of fundamental solutions. Let $(r,s)$ be the fundamental solution of its corresponding Pell equation. Then all solutions of the generalized Pell equation are

$$\{(x,y):x+y\sqrt{D}=\pm(r+s\sqrt{D})^k(u+v\sqrt{D}),k\in\mathbf Z,u+v\sqrt{D}\in U\}.$$

There are only finitely many fundamental solutions of a generalized Pell equation. Indeed, the general solution above shows that $|u+v\sqrt{D}|$ must lie between $r-s\sqrt{D}$ and $r+s\sqrt{D}$. The references at the end give sharper estimates for the ranges of the coordinates of fundamental solutions. Unlike the Pell equation, the generalized Pell equation may have no solution.

Besides combining solutions, all solutions in the same equivalence class can be obtained from a generalized Pell solution $(u,v)$ and the fundamental Pell solution $(r,s)$ by the recurrence

$$x_{k} = 2rx_{k-1} - x_{k-2},\ y_{k} = 2ry_{k-1} - y_{k-2},$$

where $x_k+y_k\sqrt{D}=(r+s\sqrt{D})^k(u+v\sqrt{D})$. This is because $x_n$ and $y_n$ can each be written, for some real pair $(A,B)$, in the form $A(r+s\sqrt{D})^k+B(r-s\sqrt{D})^k$. By Vieta's theorem, $r\pm s\sqrt{D}$ are the two real roots of $x^2-2rx+1=0$, so $x_n$ and $y_n$ both satisfy the second-order constant-coefficient recurrence above. Compared with combining solutions, this recurrence uses fewer multiplications.

Methods of solution

Both Pell equations and generalized Pell equations can be solved using continued fractions.

The PQa algorithm

The algorithms discussed here are all based on the PQa algorithm, which can find the continued-fraction expansion of a particular quadratic irrational.

Let integers $P_0,Q_0,D$ satisfy $Q_0\neq 0$, $D>0$ and not a perfect square, and $P_0^2\equiv D\pmod{Q_0}$. Then the continued-fraction expansion of the quadratic irrational

$$\omega=\dfrac{P_0+\sqrt{D}}{Q_0}$$

can be found by the following recurrence, with notation $[a_0,a_1,\cdots]$:

$$a_k = \left\lfloor\dfrac{P_k+\sqrt{D}}{Q_k}\right\rfloor,\ P_{k+1} = a_kQ_k - P_k,\ Q_{k+1} = \dfrac{D-P_{k+1}^2}{Q_k}.$$

The numerator and denominator $A_k$ and $B_k$ of the $k$-th convergent of $\omega$ are then given by the following recurrence:

$$A_k = a_kA_{k-1} + A_{k-2},\ B_k = a_kB_{k-1} + B_{k-2}$$

with $A_{-1} = 1$, $A_{-2}=0$, $B_{-1}=0$, and $B_{-2}=1$.

The correctness of these formulas was proved in the article on continued fractions. It also explains that because a quadratic irrational has a periodic continued fraction, the triple $(P_k,Q_k,a_k)$ eventually enters a cycle, so the algorithm always terminates in finitely many steps. Let the minimum period length be $\ell$, and let the earliest starting position of the cycle be $k_0$. Then the continued-fraction expansion is

$$\omega=[a_0,\cdots,a_{k_0-1},\overline{a_{k_0},\cdots,a_{k_0+\ell-1}}].$$

To solve Pell equations with the PQa algorithm, we need the following result:

Theorem

Continue with the notation above. Let $G_k=Q_0A_k-P_0B_k$. Then the integer pair $(G_{k-1},B_{k-1})$ satisfies

$$G_{k-1}^2-DB_{k-1}^2=(-1)^{k}Q_0Q_{k},$$

and their greatest common divisor $\gcd(G_{k-1},B_{k-1})$ divides $Q_{k}$.

Proof

Let the $k$-th remainder (complete quotient) in the continued-fraction expansion of $\omega$ be $\omega_{k}$, so

$$\omega = [a_0,a_1,\cdots,a_{k-1},\omega_k] = \dfrac{\omega_k A_{k-1}+A_{k-2}}{\omega_k B_{k-1}+B_{k-2}}.$$

Substitute $\omega=(P_0+\sqrt{D})/Q_0$ and $\omega_k=(P_k+\sqrt{D})/Q_k$ into the equation above to obtain

$$\dfrac{P_0+\sqrt{D}}{Q_0} = \dfrac{(P_k+\sqrt{D})A_{k-1}+Q_kA_{k-2}}{(P_k+\sqrt{D})B_{k-1}+Q_kB_{k-2}}.$$

Cancel the denominators on both sides, compare the coefficients of the rational and irrational parts, and substitute the expression for $G_k$. This gives

$$\begin{aligned} G_{k-1} &= P_kB_{k-1} + Q_kB_{k-2},\\ DB_{k-1} &= P_kG_{k-1} + Q_kG_{k-2}. \end{aligned}$$

Multiply the first equation by $G_{k-1}$ and subtract the second multiplied by $B_{k-1}$:

$$\begin{aligned} G_{k-1}^2-DB_{k-1}^2 &= (B_{k-2}G_{k-1}-B_{k-1}G_{k-2})Q_k \\ &= (A_{k-1}B_{k-2}-B_{k-1}A_{k-2})Q_0Q_k \\ &= (-1)^kQ_0Q_k. \end{aligned}$$

The last step uses the difference formula for convergents. This proves the first conclusion.

To prove the second conclusion, substitute the expression for $G_k$ into the first conclusion:

$$(Q_0A_{k-1}-P_0B_{k-1})^2 - DB_{k-1}^2 = (-1)^kQ_0Q_k.$$

Since $Q_0\mid(P_0^2-D)$, we have

$$Q_0A_{k-1}^2 +\left(\dfrac{P_0^2-D}{Q_0}B_{k-1}- 2P_0A_{k-1}\right)B_{k-1} = (-1)^kQ_k.$$

Therefore, $\gcd(G_{k-1},B_{k-1}) = \gcd(Q_0A_{k-1},B_{k-1})$ divides $Q_k$.

This result gives a way to find solutions of $x^2-Dy^2=N$. Choose $Q_0>0$ appropriately and choose $P_0$ as a solution of the congruence $P_0^2\equiv D\pmod{Q_0}$. Run the PQa algorithm on $(P_0+\sqrt{D})/Q_0$ until $(-1)^kQ_0Q_{k}=N$ is found. Then $(G_{k-1},B_{k-1})$ is a solution of the original equation. Moreover, if $Q_k=\pm 1$, this solution is primitive, meaning that $G_{k-1}$ and $B_{k-1}$ are coprime.

This idea is the core of solving Pell equations and generalized Pell equations. With it understood, we can handle some algorithmic details and prove that all solutions can be obtained this way.

Pell equation

To solve $x^2-Dy^2=1$, run the PQa algorithm from $(P_0,Q_0,D)=(0,1,D)$ until $(-1)^kQ_k=1$. Then $(A_{k-1},B_{k-1})$ is a Pell equation solution, because at this point $G_{k-1}$ is exactly $A_{k-1}$. For the Pell equation, this process can be described more precisely.

First, the solution must occur at the end of a period. The process above is the continued-fraction expansion of $\sqrt{D}$. We have the following result:

$$\sqrt{D} = [\lfloor\sqrt{D}\rfloor,\overline{a_1,\cdots,a_{\ell-1},2\lfloor\sqrt{D}\rfloor}].$$

The period length is $\ell$, and it starts at term $1$ (with indices starting at $0$). The remainder at term $\ell$ is $\lfloor\sqrt{D}\rfloor+\sqrt{D}$, which shows that $Q_{\ell}=1$. Therefore, if $\ell$ is even, $(A_{\ell-1},B_{\ell-1})$ is a nontrivial Pell equation solution; if $\ell$ is odd, $(A_{2\ell-1},B_{2\ell-1})$ is a nontrivial Pell equation solution.

We next show that this solution is the fundamental solution. There are two reasons. First, the fractions $x/y$ corresponding to all positive integer Pell solutions $(x,y)$ occur among the convergents of $\sqrt{D}$, so $(x,y)$ must be some $(A_k,B_k)$ during the PQa process. Second, $Q_k=1$ does not occur at any other position besides the end of a period, because the recurrences for $A_k$ and $B_k$ make them increase with the index. Thus the smallest positive integer solution, the fundamental solution, must occur at the position just specified. These two facts follow from the following two theorems:

Theorem

Suppose $x^2-Dy^2=N$ has a positive integer solution $(x,y)$. If $|N|<\sqrt{D}$, then $\dfrac{x}{y}$ is a convergent of $\sqrt{D}$.

Proof

When $N>0$, $x^2-Dy^2>0$, so $x>y\sqrt{D}$. Hence

$$\left|\dfrac{x}{y}-\sqrt{D}\right| = \dfrac{N}{y(x+y\sqrt{D})}<\dfrac{N}{2y^2\sqrt{D}}<\dfrac{1}{2y^2}.$$

By Legendre's criterion, $\dfrac{x}{y}$ is a convergent of $\sqrt{D}$.

When $N<0$, $x>y\sqrt{D}$ no longer holds. Instead, consider a solution of $y^2-\dfrac{1}{D}x^2=-\dfrac{N}{D}$. Since $\dfrac{|N|}{D}<\sqrt{\dfrac{1}{D}}$, the preceding argument still applies. Thus $\dfrac{y}{x}$ is a convergent of $\dfrac{1}{\sqrt{D}}$. By the reciprocal theorem, $\dfrac{x}{y}$ is also a convergent of $\sqrt{D}$.

Theorem

During the PQa algorithm above starting from $(P_0,Q_0,D)=(0,1,D)$, $Q_k=1$ implies $\ell\mid k$.

Proof

In the continued-fraction expansion of $\sqrt{D}$, every remainder other than the $0$th is a purely periodic continued fraction. Let $Q_k=1$. By Galois's result, the remainder must satisfy $\omega_k=P_k+\sqrt{D}>1$ and its conjugate $-1<P_k-\sqrt{D}<0$, which means $P_k=\lfloor\sqrt{D}\rfloor$. Therefore $\omega_k$ equals $\omega_\ell$. Repetition of a remainder means that the continued fraction has entered its cycle. If $k$ were not a multiple of $\ell$, this would contradict the minimality of $\ell$. Hence $\ell\mid k$.

Combining the discussion in this section, expanding $\sqrt{D}$ as a continued fraction—that is, running the PQa algorithm from $(P_0,Q_0,D)=(0,1,D)$—reaches the end of the first period when $Q_\ell=1$ is first obtained. If $\ell$ is even, $(A_{\ell-1},B_{\ell-1})$ is the fundamental solution; otherwise $(A_{2\ell-1},B_{2\ell-1})$ is the fundamental solution. When $\ell$ is odd, it is unnecessary to continue the PQa algorithm to twice the period: we will soon see that $A_{2\ell-1}+B_{2\ell-1}\sqrt{D}=(A_{\ell-1}+B_{\ell-1}\sqrt{D})^2$, so the fundamental solution can be computed directly from $(A_{\ell-1},B_{\ell-1})$. All other Pell equation solutions can be computed from the fundamental solution.

Example

1. Solve $x^2-14y^2=1$.

   Running the PQa algorithm from $(P_0,Q_0,D)=(0,1,14)$ gives the following results (the red entries form the first period):

   $k$	$P$	$Q$	$a$	$A$	$B$	$G$	$G^2-DB^2$
   $0$	$0$	$1$	$3$	$3$	$1$	$3$	$-5$
   $1$	$3$	$5$	$\color{red}{1}$	$4$	$1$	$4$	$2$
   $2$	$2$	$2$	$\color{red}{2}$	$11$	$3$	$11$	$-5$
   $3$	$2$	$5$	$\color{red}{1}$	$15$	$4$	$15$	$1$
   $4$	$3$	$1$	$\color{red}{6}$	$101$	$27$	$101$	$-5$
   $5$	$3$	$5$	$1$	$116$	$31$	$116$	$2$

   The period length $\ell=4$ is even. The smallest positive integer solution is $(G_3,B_3)=(15,4)$.

2. Solve $x^2-41y^2=1$.

   Running the PQa algorithm from $(P_0,Q_0,D)=(0,1,41)$ gives the following results (the red entries form the first period):

   $k$	$P$	$Q$	$a$	$A$	$B$	$G$	$G^2-DB^2$
   $0$	$0$	$1$	$6$	$6$	$1$	$6$	$-5$
   $1$	$6$	$5$	$\color{red}{2}$	$13$	$2$	$13$	$5$
   $2$	$4$	$5$	$\color{red}{2}$	$32$	$5$	$32$	$-1$
   $3$	$6$	$1$	$\color{red}{12}$	$397$	$62$	$397$	$5$
   $4$	$6$	$5$	$2$	$826$	$129$	$826$	$-5$
   $5$	$4$	$5$	$2$	$2049$	$320$	$2049$	$1$
   $6$	$6$	$1$	$12$	$25414$	$3969$	$25414$	$-5$
   $7$	$6$	$5$	$2$	$52877$	$8258$	$52877$	$5$

   The period length $\ell=3$ is odd. The smallest positive integer solution is $(G_5,B_5)=(2049,320)$. It can also be computed from $(G_2,B_2)=(32,5)$:

   $$(32+5\sqrt{41})^2=2049+320\sqrt{41}.$$

Negative Pell equation

By the discussion in the preceding section, solutions of the negative Pell equation also correspond to convergents of $\sqrt{D}$, and can occur only when $(-1)^kQ_k=-1$. This can happen only at the end of a period. Therefore, the negative Pell equation has a solution if and only if the period length $\ell$ is odd. When a solution exists, $(A_{\ell-1},B_{\ell-1})$ is its fundamental solution, and it is solved in the same way as in the preceding section.

Using an argument similar to the proof of the structure of Pell equation solutions, we can prove:

Theorem

Suppose $x^2-Dy^2=-1$ has a solution with fundamental solution $(x_1,y_1)$. Then all integer solutions of $x^2-Dy^2=\pm 1$ belong to

$$\{(x,y):x+y\sqrt{D}=\pm(x_1+y_1\sqrt{D})^k,k\in\mathbf Z\}.$$

In particular, the integer solution $(x_2,y_2)$ satisfying $x_2+y_2\sqrt{D}=(x_1+y_1\sqrt{D})^2$ is the fundamental solution of $x^2-Dy^2=1$.

Proof

By symmetry, it is enough to consider positive integer solutions $(x,y)$, that is, $x+y\sqrt{D}>1$. However, since $x^2-Dy^2=\pm 1$ consists of two hyperbolas, $x+y\sqrt{D}$ does not correspond one-to-one with $(x,y)$. To handle this difficulty, first prove that the above $(x_2,y_2)$ is the fundamental solution of $x^2-Dy^2=1$.

Clearly, $(x_2,y_2)$ solves $x^2-Dy^2=1$. If $(z,w)$ is the fundamental solution of $x^2-Dy^2=1$, then $1<z+w\sqrt{D}\le x_2+y_2\sqrt{D}$. If the right inequality is strict, divide the inequality by $x_1+y_1\sqrt{D}$ to obtain $-x_1+y_1\sqrt{D}<(z+w\sqrt{D})(-x_1+y_1\sqrt{D})<x_1+y_1\sqrt{D}$. Expanding the middle expression gives $x'+y'\sqrt{D}$, whose norm is $-1$, and $(x',y')$ is also an integer solution. Taking the reciprocal of this inequality shows that $-x'+y'\sqrt{D}$ also lies between $-x_1+y_1\sqrt{D}$ and $x_1+y_1\sqrt{D}$. The quadratic integers $\pm x'+y'\sqrt{D}$ are reciprocals of each other, so one of them must be greater than $1$. But no other quadratic integer of norm $-1$ should occur between $1$ and $x_1+y_1\sqrt{D}$, contradicting the minimality of $x_1+y_1\sqrt{D}$. Thus $x_2+y_2\sqrt{D}=z+w\sqrt{D}$, and $(x_2,y_2)$ is the fundamental solution of $x^2-Dy^2=1$.

Based on this, if a solution of $x^2-Dy^2=\pm 1$ were not of the form $(x_1+y_1\sqrt{D})^k$, there would be a $k$ such that $(x_1+y_1\sqrt{D})^{2k}<x+y\sqrt{D}<(x_1+y_1\sqrt{D})^{2k+2}$. Canceling the factor $(x_1+y_1)^{2k+1}$ would produce a quadratic integer $x'+y'\sqrt{D}\neq 1$ of norm $\pm 1$ between $-x_1+y_1\sqrt{D}$ and $x_1+y_1\sqrt{D}$. Repeating the reciprocal argument from the preceding paragraph gives a contradiction to the minimality of $x_1+y_1\sqrt{D}$. The proposition follows.

Since $(A_{\ell-1},B_{\ell-1})$ is the smallest positive integer solution of the negative Pell equation, all positive integer solutions of $x^2-Dy^2=\pm 1$ lie in

$$\{(x,y):x+y\sqrt{D}=(A_{\ell-1}+B_{\ell-1}\sqrt{D})^k,k\in\mathbf N_+\}$$

and these positive integer solutions necessarily correspond to the convergents of $\sqrt{D}$ immediately before the end of a period. Since the numerators and denominators of convergents increase strictly, for every $k\in\mathbf N_+$ we have

$$(A_{\ell-1}+B_{\ell-1}\sqrt{D})^k = A_{k\ell-1}+B_{k\ell-1}\sqrt{D}.$$

Among these positive integer solutions, odd $k$ gives negative Pell solutions and even $k$ gives Pell solutions; they alternate.

Determining whether the negative Pell equation has a solution requires computing the period length of the continued-fraction expansion of $\sqrt{D}$, which is not easy. It is therefore desirable to find a simpler test. At present, however, there is no simple and easily computable criterion⁵. We give only one simple result here.

Theorem

If $x^2-Dy^2=-1$ has a solution, then $4$ does not divide $D$ and $D$ has no prime factor of the form $4k+3$. Conversely, if $D=2$ or $D$ is a prime of the form $4k+1$, then the equation necessarily has a solution.

Proof

First, a solution of the negative Pell equation means that $-1$ is a quadratic residue modulo $D$, and therefore $-1$ is also a quadratic residue modulo every factor $d$ of $D$. Thus $d\neq 4$ and $d$ is not a prime of the form $4k+3$. Conversely, $x^2-2y^2=-1$ has the nontrivial solution $(1,1)$. It remains to handle the case where $D$ is a prime of the form $4k+1$.

Let $D$ be a prime of the form $4k+1$. We prove that $x^2-Dy^2=-1$ has a solution by starting with the fundamental solution pair $(u,v)$ of $x^2-Dy^2=1$ and constructing a solution $(\alpha,\beta)$ of $x^2-Dy^2=-1$. If $u$ were even, reducing $u^2-Dv^2=1$ modulo $4$ would give $v^2\equiv -1\pmod 4$, but $-1$ is not a quadratic residue modulo $4$. Hence $u$ is odd. Consider $Dv^2=u^2-1=(u+1)(u-1)$. Since $u$ is odd, $\gcd(u+1,u-1)=\gcd(u+1,2)=2$. Therefore, distributing the factors of $Dv^2$ between $u+1$ and $u-1$, one must be $2\alpha^2$ and the other $2D\beta^2$, where $\alpha$ and $\beta$ are coprime positive integers and $v=2\alpha\beta$. Substituting $u=\alpha^2+D\beta^2$ and $v=2\alpha\beta$ into $u^2-Dv^2=1$ gives $\alpha^2-D\beta^2=\pm 1$. Since $(u,v)$ is the fundamental Pell solution and $(\alpha,\beta)$ is a smaller positive integer pair than $(u,v)$, the right-hand side cannot be $+1$, so it must be $-1$. Thus $x^2-Dy^2=-1$ has the solution $(\alpha,\beta)$.

If $D$ is composite, having no prime factor of the form $4k+3$ and having no square factor still does not guarantee a solution to $x^2-Dy^2=-1$; for example, $x^2-34y^2=-1$ has no solution.

Example

From the calculations in the preceding examples, $x^2-14y^2=-1$ has no solution, while the smallest positive integer solution of $x^2-41y^2=-1$ is $(G_2,B_2)=(32,5)$.

The case of norm ±4

We now discuss solutions of $x^2-Dy^2=\pm 4$. Their behavior depends on $D\bmod 4$.

Some cases are easy. If $D\equiv 0\pmod 4$, then $x$ is even, so $(x/2,y)$ solves $u^2-(D/4)v^2=\pm 1$. In the remaining cases, $x,y$ must be both odd or both even. If both are odd, reducing both sides modulo $4$ gives $D\equiv 1\pmod 4$. Thus, if $D\equiv 2,3\pmod 4$, $x,y$ must both be even, and $(x/2,y/2)$ solves $u^2-Dv^2=\pm 1$. Therefore, except when $D\equiv 1\pmod 4$, solutions of $x^2-Dy^2=\pm 4$ can be obtained from the corresponding (negative) Pell equation solutions for $x,y$.

Now consider $D\equiv 1\pmod 4$, which cannot be reduced directly to a case already solved. To find the fundamental solution, apply the PQa algorithm to $(P_0,Q_0,D)=(1,2,D)$. When $Q_\ell=2$ is first obtained, the end of the first period has been reached. If the period length $l$ is even, $(G_{\ell-1},B_{\ell-1})$ is the fundamental solution of $x^2-Dy^2=4$; otherwise the pair $(G_{\ell-1},B_{\ell-1})$ is the fundamental solution of $x^2-Dy^2=-4$. Starting from the pair $(G_{\ell-1},B_{\ell-1})$, all solutions of $x^2-Dy^2=\pm 4$ are given by the following expression:

$$\left\{(x,y):\dfrac{x+y\sqrt{D}}{2}=\pm\left(\dfrac{G_{\ell-1}+B_{\ell-1}\sqrt{D}}{2}\right)^k,k\in\mathbf Z\right\}.$$

If the period length $\ell$ is even, all of these solve $x^2-Dy^2=4$; otherwise, when $k$ is odd, $(x,y)$ solves $x^2-Dy^2=-4$, and when $k$ is even, $(x,y)$ solves $x^2-Dy^2=4$.

The correctness of this algorithm depends on the following fact:

Theorem

Suppose $x^2-Dy^2=\pm 4$ has a positive integer solution $(x,y)$. If $D\equiv 1\pmod 4$, then $\dfrac{(x+y)/2}{y}$ is a convergent of $\dfrac{1+\sqrt{D}}{2}$.

Proof

First note that $x,y$ must have the same parity, so $(x+y)/2$ is an integer. If $(x,y)$ solves $x^2-Dy^2=4$, then $x>y\sqrt{D}>2y$, and hence

$$\left|\dfrac{(x+y)/2}{y}-\dfrac{1+\sqrt{D}}{2}\right| = \dfrac{2}{y(x+y\sqrt{D})}<\dfrac{1}{2y^2}.$$

By Legendre's criterion, $\dfrac{(x+y)/2}{y}$ is a convergent of $\dfrac{1+\sqrt{D}}{2}$.

If $(x,y)$ solves $x^2-Dy^2=-4$, it is enough to prove $4y<x+y\sqrt{D}$ to establish the same inequality. This holds for all cases except $D=5,13$. For $D=5,13$, substituting $x=\sqrt{Dy^2-4}$ shows that the inequality is equivalent to $2(\sqrt{D}-2)y^2>1$. Apart from $(D,y)=(5,1)$, this holds for all $D=5,13$ and positive integers $y$. In the remaining case $(D,y)=(5,1)$, the solution of $x^2-5y^2=-4$ is $(x,y)=(1,1)$, and we only need to verify that $\dfrac{1}{1}$ is a convergent of $\dfrac{1+\sqrt{5}}{2}=[\overline{1}]$, which is immediate.

Theorem

Let $D$ be a positive nonsquare integer. The continued-fraction expansion of the quadratic irrational $\omega=\dfrac{1+\sqrt{D}}{2}$ has the form

$$\omega = [\lfloor\omega\rfloor,\overline{a_1,\cdots,a_{\ell-1},2\lfloor\omega\rfloor-1}],$$

where $\ell$ is the period length and $a_k=a_{\ell-k}$ for every $1<k<\ell$.

Proof

Since $\lfloor\omega\rfloor-1+\omega>1$ and its conjugate is $\lfloor\omega\rfloor - \omega$, which lies between $-1$ and $0$, Galois's result shows that $\lfloor\omega\rfloor-1+\omega$ is a purely periodic continued fraction:

$$\lfloor\omega\rfloor-1+\omega = [\overline{2\lfloor\omega\rfloor-1,a_1,\cdots,a_{\ell-1}}].$$

Galois's result about the reciprocal of the negative conjugate gives

$$\dfrac{1}{\omega-\lfloor\omega\rfloor} = [\overline{a_{\ell-1},\cdots,a_1,2\lfloor\omega\rfloor-1}].$$

Therefore, by the definition of continued fractions,

$$\lfloor\omega\rfloor-1+\omega = 2\lfloor\omega\rfloor-1 + \dfrac{1}{\dfrac{1}{\omega-\lfloor\omega\rfloor}} = [2\lfloor\omega\rfloor-1,\overline{a_{\ell-1},\cdots,a_1,2\lfloor\omega\rfloor-1}].$$

The uniqueness of continued-fraction expansions gives $a_k=a_{\ell-k}$ for all $1<k<\ell$, proving the claimed expansion.

Theorem

Let $D\equiv 1\pmod 4$. During the PQa algorithm above starting from $(P_0,Q_0,D)=(1,2,D)$, $Q_k=2$ implies $\ell\mid k$.

Proof

In the continued-fraction expansion of $\dfrac{1+\sqrt{D}}{2}$, every remainder other than the $0$th is a purely periodic continued fraction. Let $Q_k=2$. By Galois's result, the conjugate of $\omega_k=\dfrac{P_k+\sqrt{D}}{2}$ satisfies $-1<\dfrac{P_k-\sqrt{D}}{2}<0$, or equivalently $\sqrt{D}-2<P_k<\sqrt{D}$. Since $Q_k\mid P_k^2-D$ always holds in the PQa algorithm (see the proof of algorithm correctness), $P_k$ must be odd. Thus $P_k$ is uniquely determined, namely $P_k=P_0+2(\lfloor\omega\rfloor-1)$, so $\omega_k=\omega\ell$. Repetition of a remainder means that the continued fraction has entered its cycle. If $k$ were not a multiple of $\ell$, this would contradict the minimality of $\ell$. Hence $\ell\mid k$.

Theorem

Let $(x_1,y_1)$ be the smallest positive integer solution of $x^2-Dy^2=\pm 4$. Then all its solutions are

$$\left\{(x,y):\dfrac{x+y\sqrt{D}}{2}=\pm\left(\dfrac{x_1+y_1\sqrt{D}}{2}\right)^k,k\in\mathbf Z\right\}.$$

Proof

By symmetry, it is enough to consider positive integer solutions $(x,y)$. The only thing to prove is that every pair $(x,y)$ in the set is indeed an integer solution of $x^2-Dy^2=\pm 4$. The rest repeats the proof of the solution structure for $x^2-Dy^2=\pm 1$.

Specifically, for any two integer solutions $(x_1,y_1)$ and $(x_2,y_2)$ of $x^2-Dy^2=\pm 4$, define a positive real pair $(x_3,y_3)$ by

$$\dfrac{x_3+y_3\sqrt{D}}{2} = \dfrac{x_1+y_1\sqrt{D}}{2}\dfrac{x_2+y_2\sqrt{D}}{2}.$$

Expanding the right side and comparing the rational and irrational coefficients gives

$$x_3=\dfrac{x_1x_2+Dy_1y_2}{2},\ y_3=\dfrac{x_1y_2+x_2y_1}{2}.$$

For $i=1,2$, $x_i\equiv x_i^2\equiv Dy_i^2\equiv Dy_i\pmod 2$, so

$$\begin{aligned} 2x_3 &= x_1x_2+Dy_1y_2 \equiv D^2y_1y_2+Dy_1y_2=D(D+1)y_1y_2 \equiv 0 \pmod 2,\\ 2y_3 &= x_1y_2+x_2y_1 \equiv Dy_1y_2+Dy_2y_1 = 2Dy_1y_2 \equiv 0 \pmod 2. \end{aligned}$$

Thus $x_3$ and $y_3$ are integers. Multiplicativity of the norm then shows that $(x_3,y_3)$ solves $x^2-Dy^2=\pm 4$.

Combining these facts and repeating the arguments of the preceding sections proves the correctness of the algorithm for $x^2-Dy^2=\pm 4$. These results show that $x^2-Dy^2=\pm 4$ has a simple solution structure similar to that of $x^2-Dy^2=\pm 1$: all solutions can be represented by its smallest positive integer solution, without solving any other equation.

In fact, all solutions of $x^2-Dy^2=\pm 1$ can be found among the solutions of $x^2-Dy^2=\pm 4$, so in this sense $x^2-Dy^2=\pm 4$ is more fundamental. Clearly, $(x,y)$ solves $x^2-Dy^2=\pm 1$ if and only if $(2x,2y)$ solves $x^2-Dy^2=\pm 4$. The preceding analysis shows that when $D\equiv 2,3\pmod 4$, every solution of $x^2-Dy^2=\pm 4$ is even in both coordinates, and therefore corresponds to a solution of $x^2-Dy^2=\pm 1$.

When $D\equiv 0\pmod 4$, a solution $(x,y)$ of $x^2-Dy^2=\pm 4$ always has even $x$, but $y$ may be odd. If $y_1$ is even in the smallest positive solution $(x_1,y_1)$, then $y$ is even for every solution, and these integer solutions correspond one-to-one with the integer solutions of $x^2-Dy^2=\pm 1$. If $y_1$ is odd, the parity of $y_k$ agrees with $k$ and alternates, so only even $k$ correspond to solutions of $x^2-Dy^2=\pm 1$. If $y_1$ is odd in the smallest positive solution $(x_1,y_1)$ of $x^2-Dy^2=\pm 4$ and the norm of $x_1+y_1\sqrt{D}$ is $-4$, then for such $D$, the equation $x^2-Dy^2=\pm 4$ has the solution $x^2-Dy^2=-4$, but $x^2-Dy^2=-1$ does not.

When $D\equiv 1\pmod 4$, a solution $(x,y)$ of $x^2-Dy^2=\pm 4$ may be both odd or both even. If the smallest positive solution $(x_1,y_1)$ is already even in both coordinates, then all its integer solutions are also even in both coordinates, and therefore correspond to integer solutions of $x^2-Dy^2=\pm 1$. If the smallest positive solution $(x_1,y_1)$ is odd in both coordinates, we have the following result:

Theorem

Let $(x_1,y_1)$ be the smallest positive integer solution of $x^2-Dy^2=\pm 4$. If $x_1$ and $y_1$ are both odd, then $D\equiv 5\pmod 8$, and an integer solution $(x,y)$ of this equation is even in both coordinates if and only if

$$\dfrac{x+y\sqrt{D}}{2} = \pm\left(\dfrac{x_1+y_1\sqrt{D}}{2}\right)^{3k},k\in\mathbf Z.$$

Proof

Reducing both sides of $x_1^2-Dy_1^2=\pm 4$ modulo $8$ gives $D\equiv 5\pmod 8$. To prove the second conclusion, first show that $(x_3,y_3)$ are both even. Since

$$\dfrac{x_3+y_3\sqrt{D}}{2} = \left(\dfrac{x_1+y_1\sqrt{D}}{2}\right)^{3} = \dfrac{x_1^3+3Dxy_1^2}{8}+\dfrac{3x_1^2y_1+Dy_1^3}{8}\sqrt{D},$$

it is enough to prove that the right-hand side is integral. The square of an odd number is congruent to $1$ modulo $8$, so

$$\begin{aligned} &x_1^3+3Dxy_1 = x_1(x_1^2+3Dy_1) \equiv x_1(1+3\times 5\times 1) = 16x_1 = 0 \pmod 8,\\ &3x_1^2y_1+Dy_1^3 = y_1(3x_1^2+Dy_1^2) \equiv y_1(3\times 1+5\times 1) = 8y_1 = 0 \pmod 8. \end{aligned}$$

Hence $x_3,y_3$ are both even. Therefore, for every $k\in\mathbf Z$,

$$\dfrac{x+y\sqrt{D}}{2} = \pm\left(\dfrac{x_1+y_1\sqrt{D}}{2}\right)^{3k} = \pm\left(\dfrac{x_3+y_3\sqrt{D}}{2}\right)^k \in \mathbf Z,$$

so $(x,y)$ is even in both coordinates. Conversely, for $r=1,2$,

$$\pm\left(\dfrac{x_1+y_1\sqrt{D}}{2}\right)^{3k+r} = \pm\left(\dfrac{x_3+y_3\sqrt{D}}{2}\right)^k\left(\dfrac{x_r+y_r\sqrt{D}}{2}\right).$$

To prove that the corresponding $(x,y)$ is not even, it is enough to prove that the expression is not integral, that is, that the second factor in the product is not integral. For $r=1$, this is the given condition. For $r=2$,

$$\dfrac{x_2+y_2\sqrt{D}}{2} = \left(\dfrac{x_1+y_1\sqrt{D}}{2}\right)^2 = \dfrac{x_1^2+Dy_1^2}{4} + \dfrac{x_1y_1}{2}\sqrt{D},$$

and $x_1^2+Dy_1^2\equiv 1+1\times 1=2\pmod 4$, while $x_1y_1\equiv 1\pmod 2$, so this expression is also not integral. Thus only powers divisible by $3$ give solutions that are even in both coordinates.

In other words, one out of every three solutions of $x^2-Dy^2=\pm 4$ is even in both coordinates, and corresponds to an integer solution of $x^2-Dy^2=\pm 1$. This also shows that, for $D\equiv 1\pmod 4$, $x^2-Dy^2=-4$ has a solution if and only if $x^2-Dy^2=-1$ has a solution.

The discussion so far is enough to compute the fundamental units of a real quadratic integer ring. Let $D$ be a positive integer without a square factor. For $D\equiv 2,3\pmod 4$, it is enough to find the smallest positive integer solution of $x^2-Dy^2=\pm 1$; for $D\equiv 1\pmod 4$, it is enough to find the smallest positive integer solution of $x^2-Dy^2=\pm 4$. Given the smallest positive solution $(x,y)$, the fundamental units are $\pm x\pm y\sqrt{D}$ for $D\equiv 2,3\pmod 4$, and $\dfrac{\pm x\pm y\sqrt{D}}{2}$ for $D\equiv 1\pmod 4$.

Example

1. Solve $x^2-14y^2=\pm 4$.

   The preceding example shows that the smallest positive integer solution of $x^2-14y^2=4$ is $(30,8)$, while $x^2-14y^2=-4$ has no solution.

2. Solve $x^2-41y^2=\pm 4$.

   Running the PQa algorithm from $(P_0,Q_0,D)=(1,2,41)$ gives the following results (the red entries form the first period):

   $k$	$P$	$Q$	$a$	$A$	$B$	$G$	$G^2-DB^2$
   $0$	$1$	$2$	$3$	$3$	$1$	$5$	$-16$
   $1$	$5$	$8$	$\color{red}{1}$	$4$	$1$	$7$	$8$
   $2$	$3$	$4$	$\color{red}{2}$	$11$	$3$	$19$	$-8$
   $3$	$5$	$4$	$\color{red}{2}$	$26$	$7$	$45$	$16$
   $4$	$3$	$8$	$\color{red}{1}$	$37$	$10$	$64$	$-4$
   $5$	$5$	$2$	$\color{red}{5}$	$211$	$57$	$365$	$16$
   $6$	$5$	$8$	$1$	$248$	$67$	$429$	$-8$
   $7$	$3$	$4$	$2$	$707$	$191$	$1223$	$8$
   $8$	$5$	$4$	$2$	$1662$	$449$	$2875$	$-16$
   $9$	$3$	$8$	$1$	$2369$	$640$	$4098$	$4$
   $10$	$5$	$2$	$5$	$13507$	$3649$	$23365$	$-16$
   $11$	$5$	$8$	$1$	$15876$	$4289$	$27463$	$8$

   The period length $\ell=5$ is odd. The minimal positive integer solution of $x^2-41y^2=-4$ is $(G_4,B_4)=(64,10)$, and the minimal positive integer solution of $x^2-41y^2=4$ is $(G_9,B_9)=(4098,640)$. They satisfy:

   $$\dfrac{4098+640\sqrt{41}}{2} = \left(\dfrac{64+10\sqrt{41}}{2}\right)^2.$$

   Of course, because $D\equiv 1\pmod 8$, the conclusions above show that the minimal positive integer solutions of $x^2-41y^2=\pm 4$ in this case are always even, and are always the minimal positive integer solutions of $x^2-41y^2=\pm 1$ multiplied by $2$. Thus, the example above can also be used directly.

3. Solve the equation $x^2-13y^2=\pm 4$.

   Running the PQa algorithm with $(P_0,Q_0,D)=(1,2,13)$ gives the following result (the red portion is the first period):

   $k$	$P$	$Q$	$a$	$A$	$B$	$G$	$G^2-DB^2$
   $0$	$1$	$2$	$2$	$2$	$1$	$3$	$-4$
   $1$	$3$	$2$	$\color{red}{3}$	$7$	$3$	$11$	$4$
   $2$	$3$	$2$	$3$	$23$	$10$	$36$	$-4$
   $3$	$3$	$2$	$3$	$74$	$33$	$119$	$4$

   The period length $\ell=1$ is odd. The minimal positive integer solution of $x^2-13y^2=-4$ is $(G_0,B_0)=(3,1)$, and the minimal positive integer solution of $x^2-13y^2=4$ is $(G_1,B_1)=(11,3)$.

   Because the minimal positive integer solutions of this equation are all odd, the pair $(G_2,B_2)=(36,10)$ at the end of the third period can be used to obtain the corresponding (negative) Pell equation $x^2-13y^2=\pm 1$ with minimal positive integer solution $(18,5)$. It can also be obtained directly:

   $$\dfrac{36+10\sqrt{13}}{2}=\left(\dfrac{3+\sqrt{13}}{2}\right)^3.$$

   Moreover, this is a solution of the negative Pell equation. The minimal positive integer solution of the corresponding Pell equation is $(649,180)$.

4. Solve the equation $x^2-52y^2=\pm 4$.

   Because the minimal positive integer solutions of $x^2-13y^2=\pm 1$ are $(18,5)$ and $(649,180)$, the minimal positive integer solutions of $x^2-52y^2=\pm 4$ are $(36,10)$ and $(1298,360)$.

General case

Finally, we discuss how to solve the generalized Pell equation.

For the case $|N|<\sqrt{D}$, there is a simple solution. The conclusions above show that the equation $x^2-Dy^2=N$ has a solution $(x,y)$, which must satisfy that $\dfrac{x}{y}$ is a convergent of $\sqrt{D}$. Moreover, from the structure of solutions discussed above, every fundamental solution $(x,y)$ satisfies $x+y\sqrt{D}$ no greater than the corresponding fundamental solution of the Pell equation $x^2-Dy^2=1$, namely $x_1+y_1\sqrt{D}$. By the monotonicity of the denominator sequence $B_k$ in the PQa algorithm, these fundamental solutions of the generalized Pell equation must occur before the fundamental solution of the corresponding Pell equation. Therefore, run the PQa algorithm with $(P_0,Q_0,D)=(0,1,D)$ until $Q_{\ell'}=1$ and $\ell'$ is even. During this process, for every $(A_k,B_k)$ that appears, check whether there is an integer $f$ such that

$$A_k^2-DB_k^2 = (-1)^{k+1}Q_{k+1} = N/f^2$$

If this holds, record $(fA_{k},fB_{k})$ as a minimal positive integer solution. All $(fA_k,fB_k)$ recorded in this process are exactly the minimal positive integer solutions of $x^2-Dy^2=N$. Using $(A_{\ell'-1},B_{\ell'-1})$, the fundamental solution of the corresponding Pell equation, we can generate all solutions of the generalized Pell equation from these minimal positive integer solutions. Note that, depending on whether the period length $\ell$ is even or odd, the $\ell'$ above may be either $\ell$ or $2\ell$.

For a more general $N$, the above method no longer applies. First enumerate all square factors of $N$, writing each as $f^2$; let $m=N/f^2$, and enumerate all solutions to the congruence $z^2\equiv D\pmod{|m|}$ satisfying $-|m|/2<z \le |m|/2$, denoting each solution by $z$. Then run the PQa algorithm with $(P_0,Q_0,D)=(z,|m|,D)$ until $Q_k=\pm 1$ or one period has ended. In the second case, the equation associated with this pair $(f,z)$ has no solution. In the first case, further determine whether $(-1)^kQ_k=N/|N|$. If the signs agree, $(fG_{k-1},fB_{k-1})$ is a solution of $x^2-Dy^2=N$. Otherwise, it is a solution of $x^2-Dy^2=-N$, and it can be composed with the fundamental solution of the corresponding negative Pell equation to obtain a solution of $x^2-Dy^2=N$ if and only if that negative Pell equation has a solution. After traversing all pairs $(f,z)$, we obtain exactly one solution from each equivalence class of $x^2-Dy^2=N$, and that solution is the fundamental solution or minimal positive integer solution in its class. These solutions and the fundamental solution of the corresponding Pell equation can be used to generate all integer solutions of the equation. This algorithm is called the Lagrange–Matthews–Mollin algorithm.

The correctness of this algorithm follows from the theorem below:

Theorem

Suppose the equation $x^2-Dy^2=N$ has an integer solution $(x,y)$ with $x\ge 0, y>0,\gcd(x,y)=1$. Let $Q_0=|N|$, so $\gcd(Q_0,y)=1$. Let $P_0$ be a solution of the congruence $x\equiv -P_0y\pmod{Q_0}$ with $-Q_0/2<P_0\le Q_0/2$, and let the integer $X$ satisfy $x=Q_0X-P_0y$. Then, $P_0^2\equiv D\pmod{Q_0}$, and $\dfrac{X}{y}$ is a convergent of $\omega=\dfrac{P_0+\sqrt{D}}{Q_0}$, namely $\dfrac{A_{k-1}}{B_{k-1}}$, with $Q_k=(-1)^k\dfrac{N}{|N|}$.

Proof

Using $x\equiv -P_0y\pmod{Q_0}$ and $x^2-Dy^2=N\equiv 0\pmod{Q_0}$; it is clear that $P_0^2\equiv D\pmod{Q_0}$. Therefore,

$$P_0x+Dy\equiv -P_0^2y+Dy = (D-P_0^2)y\equiv 0\pmod{Q_0}.$$

Thus, consider the integer matrix

$$\begin{pmatrix}P & R \\ Q & S\end{pmatrix} = \begin{pmatrix}X & \dfrac{P_0x+Dy}{Q_0} \\ y & x\end{pmatrix}.$$

Its determinant is

$$PS-QR = \dfrac{x(x+P_0y)-y(P_0x+Dy)}{Q_0} = \dfrac{x^2-Dy^2}{Q_0} = \pm 1.$$

Moreover, let $\zeta =\sqrt{D} > 1$; then

$$\dfrac{P\zeta+R}{Q\zeta+S} = \dfrac{(x+P_0y)\sqrt{D}+(P_0x+Dy)}{(x+y\sqrt{D})Q_0} = \dfrac{P_0+\sqrt{D}}{Q_0} = \omega.$$

We now prove that $\dfrac{P}{Q}$ is a convergent of $\omega$. Without loss of generality, suppose $\dfrac{P}{Q}$ has continued fraction expansion

$$\dfrac{P}{Q} = [a_0,a_1,\cdots,a_k]$$

and $PS-QR = (-1)^{k-1}$. Suppose $\dfrac{p_k}{q_k}$ is its $k$-th convergent. Then $(p_k,q_k)=(P,Q)$, and the difference formula for convergents shows that $p_kq_{k-1}-q_kp_{k-1}=(-1)^{k-1}$. This shows that

$$p_k(S-q_{k-1}) = q_k(R-p_{k-1}).$$

Consider the following cases:

• If $S=0$, it is easy to verify $Q=R=1$, so $\omega=P+\zeta^{-1}=[P,\zeta]$, and therefore $\dfrac{P}{Q}=P$ is a convergent of $\omega$ with index $0$;
• If $Q=S>0$, then $Q=S=1$ and $P-R=\pm 1$. In this case,
  • If $P=R+1$, then $\omega=R+\dfrac{1}{1+\zeta^{-1}}=[R,1,\zeta]$, so $\dfrac{P}{Q}=\dfrac{R+1}{1}=[R,1]$ is a convergent of $\omega$ with index $1$;
  • If $P=R-1$, then $\omega=R-1+\dfrac{1}{1+\zeta}=[R-1,\zeta-1]$, so $\dfrac{P}{Q}=R-1$ is a convergent of $\omega$ with index $0$;
• If $Q\neq S>0$, then, since $Q=q_k\mid(S-q_{k-1})$, there is always an integer $\kappa$ such that $S=\kappa q_k+q_{k-1}$ and $R=\kappa p_k+p_{k-1}$. Since $q_k\ge q_{k-1}$ and $S>0$, we have $\kappa\ge 0$. Therefore, $\omega=\dfrac{(\kappa+\zeta)p_k+p_{k-1}}{(\kappa+\zeta)q_k+q_{k-1}}=[a_0,a_1,\cdots,a_k,\kappa+\zeta]$, so $\dfrac{P}{Q}$ is its $k$-th convergent.

In summary, $\dfrac{X}{y}$ is a convergent of $\omega=\dfrac{P_0+\sqrt{D}}{Q_0}$, denoted by $\dfrac{A_{k-1}}{B_{k-1}}$. Since $A_{k-1}^2-DB_{k-1}^2=(-1)^kQ_0Q_k$, so $Q_k=(-1)^k\dfrac{N}{|N|}$.

The theorem guarantees that all positive solutions of the equation occur among the convergents of the corresponding quadratic irrational. When convergents are computed with the PQa algorithm, once the period is reached, the convergents are guaranteed to be positive. Therefore, by enumerating all quadratic irrationals allowed by the conditions of the theorem and computing their convergents through one period, we can find a solution. Since two solutions appearing among the convergents of the same quadratic irrational must be equivalent, we can stop as soon as we obtain the first solution satisfying $(-1)^kQ_k=N/|N|$. Unlike all the preceding algorithms, here the qualifying $k$ may occur before the period has been entered.

Example

1. Solve the equation $x^2-157y^2=12$.

   Since $12^2<157$, running the PQa algorithm with $(P_0,Q_0,D)=(0,1,157)$ gives the following result (the red portion is the first period):

   $k$	$P$	$Q$	$a$	$A$	$B$	$G$	$G^2-DB^2$
   $0$	$0$	$1$	$12$	$12$	$1$	$12$	$-13$
   $1$	$12$	$13$	$\color{red}1$	$13$	$1$	$13$	$12$
   $2$	$1$	$12$	$\color{red}1$	$25$	$2$	$25$	$-3$
   $3$	$11$	$3$	$\color{red}7$	$188$	$15$	$188$	$19$
   $4$	$10$	$19$	$\color{red}1$	$213$	$17$	$213$	$-4$
   $5$	$9$	$4$	$\color{red}5$	$1253$	$100$	$1253$	$9$
   $6$	$11$	$9$	$\color{red}2$	$2719$	$217$	$2719$	$-12$
   $7$	$7$	$12$	$\color{red}1$	$3972$	$317$	$3972$	$11$
   $8$	$5$	$11$	$\color{red}1$	$6691$	$534$	$6691$	$-11$
   $9$	$6$	$11$	$\color{red}1$	$10663$	$851$	$10663$	$12$
   $10$	$5$	$12$	$\color{red}1$	$17354$	$1385$	$17354$	$-9$
   $11$	$7$	$9$	$\color{red}2$	$45371$	$3621$	$45371$	$4$
   $12$	$11$	$4$	$\color{red}5$	$244209$	$19490$	$244209$	$-19$
   $13$	$9$	$19$	$\color{red}1$	$289580$	$23111$	$289580$	$3$
   $14$	$10$	$3$	$\color{red}7$	$2271269$	$181267$	$2271269$	$-12$
   $15$	$11$	$12$	$\color{red}1$	$2560849$	$204378$	$2560849$	$13$
   $16$	$1$	$13$	$\color{red}1$	$4832118$	$385645$	$4832118$	$-1$
   $17$	$12$	$1$	$\color{red}{24}$	$118531681$	$9459858$	$118531681$	$13$
   $18$	$12$	$13$	$1$	$123363799$	$9845503$	$123363799$	$-12$
   $19$	$1$	$12$	$1$	$241895480$	$19305361$	$241895480$	$3$
   $20$	$11$	$3$	$7$	$1816632159$	$144983030$	$1816632159$	$-19$
   $21$	$10$	$19$	$1$	$2058527639$	$164288391$	$2058527639$	$4$
   $22$	$9$	$4$	$5$	$12109270354$	$966424985$	$12109270354$	$-9$
   $23$	$11$	$9$	$2$	$26277068347$	$2097138361$	$26277068347$	$12$
   $24$	$7$	$12$	$1$	$38386338701$	$3063563346$	$38386338701$	$-11$
   $25$	$5$	$11$	$1$	$64663407048$	$5160701707$	$64663407048$	$11$
   $26$	$6$	$11$	$1$	$103049745749$	$8224265053$	$103049745749$	$-12$
   $27$	$5$	$12$	$1$	$167713152797$	$13384966760$	$167713152797$	$9$
   $28$	$7$	$9$	$2$	$438476051343$	$34994198573$	$438476051343$	$-4$
   $29$	$11$	$4$	$5$	$2360093409512$	$188355959625$	$2360093409512$	$19$
   $30$	$9$	$19$	$1$	$2798569460855$	$223350158198$	$2798569460855$	$-3$
   $31$	$10$	$3$	$7$	$21950079635497$	$1751807067011$	$21950079635497$	$12$
   $32$	$11$	$12$	$1$	$24748649096352$	$1975157225209$	$24748649096352$	$-13$
   $33$	$1$	$13$	$1$	$46698728731849$	$3726964292220$	$46698728731849$	$1$
   $34$	$12$	$1$	$24$	$1145518138660728$	$91422300238489$	$1145518138660728$	$-13$
   $35$	$12$	$13$	$1$	$1192216867392577$	$95149264530709$	$1192216867392577$	$12$

   The period length $\ell=17$ is odd, so within two periods we need to consider $G_{k-1}^2-157B_{k-1}^2$ differs from $12$ by a square factor, namely $k=1,9,13,19,23,31$; the corresponding solutions are the values of $(fG,fB)$:

   $k$	$f$	$fG_{k-1}$	$fB_{k-1}$	$x$	$y$
   $1$	$1$	$13$	$1$	$13$	$1$
   $9$	$1$	$10663$	$851$	$10663$	$851$
   $13$	$2$	$579160$	$46222$	$579160$	$46222$
   $19$	$2$	$483790960$	$38610722$	$-579160$	$46222$
   $23$	$1$	$26277068347$	$2097138361$	$-10663$	$851$
   $31$	$1$	$21950079635497$	$1751807067011$	$-13$	$1$

   All $(fG,fB)$ are the minimal positive integer solutions in all equivalence classes of $x^2-157y^2=12$. To obtain all solutions from these, we can use the fundamental solution of the corresponding Pell equation $(46698728731849,3726964292220)$. For example, they can be transformed into the fundamental solutions $(x,y)$, and the corresponding solutions are also listed in the table above.

2. Solve the equation $x^2-157y^2=12$.

   This time, use the Lagrange–Matthews–Mollin algorithm to solve it. First, enumerate $N=12$ square factors:

   • $f^2=1^2$ gives $m=12$, and the congruence $P^2\equiv 157\pmod{12}$ has solutions $z=\pm 1,\pm 5$;
   • $f^2=2^2$ gives $m=3$, and the congruence $P^2\equiv 157\pmod{3}$ has solutions $z=\pm 1$.

   For every possible $(f,z)$ combination, run the PQa algorithm with initial parameters $(P_0,Q_0,D)=(z,|m|,D)$ and find the first occurrence of $(-1)^kQ_k=1$; the corresponding $(fG_{k-1},fB_{k-1})$ is a solution. The results are shown in the following table:

   $f$	$z$	$m$	$k$	$fG_{k-1}$	$fB_{k-1}$
   $1$	$1$	$12$	$32$	$21950079635497$	$1751807067011$
   $1$	$-1$	$12$	$2$	$13$	$1$
   $1$	$5$	$12$	$24$	$26277068347$	$2097138361$
   $1$	$-5$	$12$	$10$	$10663$	$851$
   $2$	$1$	$3$	$20$	$483790960$	$38610722$
   $2$	$-1$	$3$	$14$	$579160$	$46222$

   These are the minimal positive integer solutions in all equivalence classes listed above; the fundamental solution of the Pell equation can be used to transform them into fundamental solutions.

3. Solve the equation $x^2-79y^2=\pm 101$.

   Again, use the Lagrange–Matthews–Mollin algorithm to solve it. Since $N=101$ is prime, so necessarily $f=1$. Then, $m=101$, the corresponding congruence $P^2\equiv 79\pmod{101}$ has solutions $P=\pm 33$.

   Running the PQa algorithm with $(P_0,Q_0,D)=(33,101,79)$ gives the following result (the red portion is the first period):

   $k$	$P$	$Q$	$a$	$A$	$B$	$G$	$G^2-DB^2$
   $0$	$33$	$101$	$0$	$0$	$1$	$-33$	$1010$
   $1$	$-33$	$-10$	$2$	$1$	$2$	$35$	$909$
   $2$	$13$	$9$	$2$	$2$	$5$	$37$	$-606$
   $3$	$5$	$6$	$\color{red}2$	$5$	$12$	$109$	$505$
   $4$	$7$	$5$	$\color{red}3$	$17$	$41$	$364$	$-303$
   $5$	$8$	$3$	$\color{red}5$	$90$	$217$	$1929$	$1010$
   $6$	$7$	$10$	$\color{red}1$	$107$	$258$	$2293$	$-707$
   $7$	$3$	$7$	$\color{red}1$	$197$	$475$	$4222$	$909$
   $8$	$4$	$9$	$\color{red}1$	$304$	$733$	$6515$	$-606$
   $9$	$5$	$6$	$2$	$805$	$1941$	$17252$	$505$

   The period length $\ell=6$ is even. There is no $Q_k=\pm 1$, so this case has no solution. Similarly, running the PQa algorithm with $(P_0,Q_0,D)=(-33,101,79)$ also shows a similar situation. Therefore, this equation has no solution.

Exercises

• LOJ 6687. Project Euler 66 — Solve the equation
• SPOJ EQU2 - Yet Another Equation
• SPOJ PELL2 - Pell (Mid pelling)
• UVa 12909. Numeric Center
• UVa 10241. Semi-triangular and also Square

References and notes

• Pell's equation - Wikipedia
• John P. Robertson - Solving the generalized Pell equation $x^2-Dy^2=N$
• Keith Matthews - The Diophantine Equation $x^2-Dy^2=N$,$D>0$
• Existence of Solution to Pell’s Equation - Suryateja Gavva's Blog
• Calculating the simple continued fraction of a quadratic irrational - Number Theory Web(PQa algorithm)
• Solving the diophantine equation x2–Dy2 = N, D > 0 and not a perfect square, N ≠ 0 - Number Theory Web(Lagrange–Matthews–Mollin algorithm)

Footnotes

1. When $D$ is a perfect square, direct factorization gives $(x+y\sqrt{D})(x-y\sqrt{D})=N$, so all solutions can be found by enumerating the factors of $N$. In particular, when $N=1$, the equation has only the solutions $(\pm 1,0)$; when $N=-1$ and $D\neq 1$, the equation has no solution. ↩

2. Some Chinese sources also call it the negative Pell equation of the second type. ↩

3. That is, a rational number of the form $n+\dfrac12$ with $n\in\mathbf Z$. ↩

4. Note that the definition of a fundamental solution of a Pell equation is not the same as the definition of a fundamental unit in a real quadratic integer ring. First, in some real quadratic integer rings, a fundamental unit $x+y\sqrt{D}$ has half-integer $x,y$, so it is not a solution of a Pell equation. Second, a real quadratic integer ring has four fundamental units, but there is only one fundamental solution, because a fundamental solution requires $x,y$ to both be positive. ↩

5. A relatively practical test and tool can be found here and in its references. The equation $x^2-Dy^2=-1$ has a solution for exactly the positive integers $D$ listed in OEIS A031396. ↩