Skip to content

Pell's Equation

Prerequisite knowledge: Continued fractions, Quadratic fields

Introduction

This article discusses the solution of (generalized) Pell's equation. The generalized Pell's equation refers to the Diophantine equation in \(x\) and \(y\):

\[ x^2 - D y^2 = N, \]

where \(D\) is a positive integer and not a perfect square1, and \(N\) is a nonzero integer. The narrow sense Pell's equation specifically refers to the special cases \(N = 1\) or \(N = \pm 1\), and sometimes also includes \(N = \pm 4\). The generalized Pell's equation is closely related to finding quadratic integers with norm \(N\) in the ring of integers of a real quadratic field, while the cases often called (narrow) Pell's equation can be viewed as finding units in the ring of integers of a real quadratic field.

When this article mentions Pell's equation, it specifically refers to the case \(N = 1\). Correspondingly, the case \(N = -1\) is called the negative Pell's equation2.

Structure of Solutions

The integer solutions \((x, y)\) of the generalized Pell's equation are closely related to the quadratic integer \(x + y\sqrt{D}\). Therefore, solutions to Pell's equation in the literature are often written in the form \(x + y\sqrt{D}\). Because the norm of a quadratic integer is

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

the generalized Pell's equation roughly corresponds to solving for quadratic integers with norm \(N\). However, there is a subtle difference. When \(x\) and \(y\) are both integers, \(x + y\sqrt{D}\) must be a quadratic integer; conversely, a quadratic integer does not necessarily require \(x\) and \(y\) to be integers—when \(D \equiv 1 \pmod 4\), \(x\) and \(y\) can also be half-integers3.

This distinction is particularly important when solving for fundamental units. Because units in the ring of quadratic integers are those with norm \(\pm 1\). For \(D \equiv 2, 3 \pmod 4\), to find such units, it suffices to solve the generalized Pell's equation with \(N = \pm 1\); but for \(D \equiv 1 \pmod 4\), we also need to consider the case \(N = \pm 4\). Later we will discuss methods for finding units.

To understand the structure of solutions to the generalized Pell's equation, we start from Brahmagupta's identity:

\[ (x_1^2 - D y_1^2)(x_2^2 - D y_2^2) = (x_1 x_2 + D y_1 y_2)^2 - D (x_1 y_2 + x_2 y_1)^2. \]

This is equivalent to the norm preserving multiplication for quadratic integers, i.e.,

\[ \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_1 x_2 + D y_1 y_2) + (x_1 y_2 + x_2 y_1)\sqrt{D}\right). \end{aligned} \]

Using this identity, we can combine integer solutions of \(x^2 - D y^2 = N_1\) and \(x^2 - D y^2 = N_2\) to obtain integer solutions of \(x^2 - D y^2 = N_1 N_2\). Of course, from the perspective of quadratic integers, combining solutions is just multiplication of quadratic integers, which illustrates the convenience of writing Pell's equation solutions as quadratic integers. In particular, taking \(N_1 = N\) and \(N_2 = 1\) shows that if we know one solution of \(x^2 - D y^2 = N\) and all solutions of the corresponding Pell's equation \(x^2 - D y^2 = 1\), we can obtain more solutions of \(x^2 - D y^2 = N\). However, this method may not generate all solutions. But it at least shows that understanding the structure of Pell's equation solutions is important for understanding the structure of generalized Pell's equation solutions.

Pell's Equation

The geometric meaning of the equation \(x^2 - D y^2 = 1\) is a hyperbola with the real axis on the \(x\)-axis and the imaginary axis on the \(y\)-axis. Each point on the hyperbola uniquely corresponds to a nonzero value of \(x + y\sqrt{D}\): the left branch corresponds to negative values of \(x + y\sqrt{D}\), and the right branch corresponds to positive values. Moreover, on each branch, the value of \(x + y\sqrt{D}\) corresponding to the hyperbola from bottom to top is strictly increasing. The values of quadratic integers give a natural ordering to the solutions of Pell's equation.

The hyperbola is symmetric about both the \(x\)-axis and the \(y\)-axis, so discussing solutions of Pell's equation only requires considering the segment in the first quadrant; the remaining solutions can be obtained by symmetry. This is equivalent to only considering solutions with \(x + y\sqrt{D} > 1\). If the equation has nontrivial solutions besides \((\pm 1, 0)\), then there must exist a solution \((x_1, y_1)\) in the first quadrant with the smallest value of \(x + y\sqrt{D}\). This is also the integer point in the first quadrant (excluding the axes) with the smallest coordinates, and it is called the fundamental solution4 of Pell's equation. According to the previous discussion, the integer pairs \((x_k, y_k)\) satisfying \(x_k + y_k\sqrt{D} = (x_1 + y_1\sqrt{D})^k\) are all solutions of Pell's equation and lie in the first quadrant. Conversely, these are indeed all solutions in the first quadrant. Using symmetry, we obtain the following conclusion:

Theorem

Let the fundamental solution of Pell's equation \(x^2 - D y^2 = 1\) be \((x_1, y_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 exists another solution \(x + y\sqrt{D}\) such that 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 the integer point \((x, y)\) lies on the hyperbola between \((x_k, y_k)\) and \((x_{k+1}, y_{k+1})\) (excluding endpoints). 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}) = (x x_k - D y y_k) + (x_k y - x y_k)\sqrt{D} < x_1 + y_1\sqrt{D}. \]

By the monotonicity mentioned earlier, this inequality implies that \((x x_k - D y y_k, x_k y - x y_k)\) is an integer point between \((1, 0)\) and \((x_1, y_1)\). This contradicts the choice of \((x_1, y_1)\).

Extending the solutions from the first quadrant to the entire plane, taking the exponent \(k\) to its opposite (i.e., taking the reciprocal) corresponds to symmetry about the \(x\)-axis, and taking the negative corresponds to symmetry about the origin. Adding the trivial solution when \(k = 0\) yields all solutions of Pell's equation.

The previous discussion only assumed the existence of a fundamental solution. Now we show that Pell's equation always has nontrivial solutions.

Theorem

Pell's equation \(x^2 - D y^2 = 1\) always has integer solutions other than \((\pm 1, 0)\).

Proof

First, Dirichlet's theorem shows that there are infinitely many pairs of positive integers \((x, y)\) such that

\[ \left|\frac{x}{y} - \sqrt{D}\right| \le \frac{1}{y^2} \]

holds. They all satisfy the inequality

\[ |x^2 - D y^2| = y^2 \left|\frac{x}{y} - \sqrt{D}\right| \left|\frac{x}{y} + \sqrt{D}\right| \le \frac{1}{y^2} + 2\sqrt{D} < 1 + 2\sqrt{D}. \]

Thus, there must exist an integer \(m \in (-1 - 2\sqrt{D}, 1 + 2\sqrt{D})\) such that infinitely many pairs \((x, y)\) satisfy \(x^2 - D y^2 = m\). Classifying these \((x, y)\) by their residues modulo \(m\), we know that for some pair of integers \((x_0, y_0)\), there exist infinitely many pairs \((x, y)\) such that \(x \equiv x_0 \pmod m\) and \(y \equiv y_0 \pmod m\). Take two distinct pairs \((x_1, y_1)\) and \((x_2, y_2)\) satisfying these conditions. Then

\[ \frac{x_1 + y_1\sqrt{D}}{x_2 + y_2\sqrt{D}} = \frac{x_1 x_2 - D y_1 y_2}{m} + \frac{x_2 y_1 - x_1 y_2}{m}\sqrt{D}. \]

Because by the congruence relations,

\[ \begin{aligned} x_1 x_2 - D y_1 y_2 &\equiv x_0^2 - D y_0^2 = m \equiv 0 \pmod{|m|},\\ x_2 y_1 - x_1 y_2 &\equiv x_0 y_0 - x_0 y_0 = 0 \pmod{|m|}, \end{aligned} \]

this shows that the right-hand side yields an integer solution. Moreover, since \((x_1, y_1) \neq (x_2, y_2)\), it is nontrivial. Thus, Pell's equation indeed has nontrivial solutions.

Of course, this section provides a non-constructive proof. In the following discussion of methods for solving Pell's equation, we will directly use convergents of continued fractions to construct solutions, providing another proof of the existence of nontrivial solutions. Additionally, although the structure of Pell's equation solutions obtained in this section is consistent with the structure of units in the ring of integers of a real quadratic field, for the case \(D \equiv 1 \pmod 4\), this section does not fully resolve the structure of units in the corresponding ring of integers. This will be further discussed below.

Generalized Pell's Equation

The graph of the generalized Pell's equation \(x^2 - D y^2 = N\) is also a hyperbola in the plane, symmetric about the \(x\)-axis and \(y\)-axis. As noted earlier, some solutions of \(x^2 - D y^2 = N\) may differ only by a factor that is a solution of Pell's equation, meaning we can partition the solutions of \(x^2 - D y^2 = N\) into equivalence classes. For two solutions \((x_1, y_1)\) and \((x_2, y_2)\) of \(x^2 - D y^2 = N\), if there exists a solution \((u, v)\) of Pell's equation 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 said to be equivalent. A necessary and sufficient condition for two solutions to be equivalent is

\[ N \mid (x_1 x_2 - D y_1 y_2),\quad N \mid (x_2 y_1 - x_1 y_2). \]

Since solutions of Pell's equation are relatively easy to find, a natural idea is to find one solution in each equivalence class. Knowing these solutions, we can use the corresponding Pell's equation solutions to obtain all solutions of the generalized Pell's equation. Due to symmetry, in each equivalence class of solutions of the generalized Pell's equation, there exists a solution with nonnegative \(y\) but as small as possible: if such a solution is unique, it is called the fundamental solution of that equivalence class; otherwise, the equivalence class must have two solutions with nonnegative and minimal \(y\), symmetric about the \(y\)-axis, and we choose the one with \(x > 0\) as the fundamental solution. Thus, solving the generalized Pell's equation \(x^2 - D y^2 = N\) is equivalent to finding its set of fundamental solutions \(U\). Let the fundamental solution of the corresponding Pell's equation be \((r, s)\). Then the set of all solutions of the generalized Pell's equation is

\[ \{(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\}. \]

The fundamental solutions of the generalized Pell's equation must be finite. This is because from the general solution expression above, the absolute value \(|u + v\sqrt{D}|\) must lie between \(r - s\sqrt{D}\) and \(r + s\sqrt{D}\). The references at the end provide stricter bounds on the coordinates of fundamental solutions. Of course, unlike Pell's equation, the generalized Pell's equation may have no solutions.

Besides using solution composition, another method to obtain all solutions in the same equivalence class from one solution \((u, v)\) of the generalized Pell's equation and the fundamental solution \((r, s)\) of Pell's equation is to use the recurrence relation:

\[ x_{k} = 2r x_{k-1} - x_{k-2},\quad y_{k} = 2r y_{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 be written in the form \(A (r + s\sqrt{D})^k + B (r - s\sqrt{D})^k\) for some pair of real numbers \((A, B)\). By Vieta's theorem, \(r \pm s\sqrt{D}\) are the two real roots of \(x^2 - 2r x + 1 = 0\), so \(x_n\) and \(y_n\) satisfy the above second-order linear recurrence relation. Compared to solution composition, this recurrence formula involves fewer multiplications.

Solution Methods

Both Pell's equation and the generalized Pell's equation can be solved based on continued fractions.

PQa Algorithm

The algorithms discussed in this article are based on the PQa algorithm, which can be used to find the continued fraction expansion of a specific 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 \([a_0, a_1, \cdots]\) of the quadratic irrational

\[ \omega = \frac{P_0 + \sqrt{D}}{Q_0} \]

can be obtained by the following recurrence formulas:

\[ a_k = \left\lfloor\frac{P_k + \sqrt{D}}{Q_k}\right\rfloor,\quad P_{k+1} = a_k Q_k - P_k,\quad Q_{k+1} = \frac{D - P_{k+1}^2}{Q_k}. \]

Furthermore, the numerator and denominator \(A_k\) and \(B_k\) of the \(k\)-th convergent of \(\omega\) are given by the following recurrence formulas:

\[ A_k = a_k A_{k-1} + A_{k-2},\quad B_k = a_k B_{k-1} + B_{k-2} \]

with \(A_{-1} = 1\), \(A_{-2} = 0\), \(B_{-1} = 0\), \(B_{-2} = 1\).

The correctness of these formulas has been proven in the continued fractions article. It also shows that because quadratic irrationals are periodic continued fractions, the triple \((P_k, Q_k, a_k)\) will eventually enter a cycle, and the algorithm always terminates in finite steps. Let the minimal period length be \(\ell\), and the earliest starting position of the cycle be \(k_0\). Then the continued fraction expansion of the quadratic irrational can be written as

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

To use the PQa algorithm to solve Pell's equation, we need the following result:

Theorem

Continue with the above notation. Let \(G_k = Q_0 A_k - P_0 B_k\). Then the integer pair \((G_{k-1}, B_{k-1})\) satisfies

\[ G_{k-1}^2 - D B_{k-1}^2 = (-1)^k Q_0 Q_k, \]

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

Proof

Let the \(k\)-th complete quotient in the continued fraction expansion of \(\omega\) be \(\omega_k\), i.e.,

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

Substituting \(\omega = (P_0 + \sqrt{D})/Q_0\) and \(\omega_k = (P_k + \sqrt{D})/Q_k\) into the above gives

\[ \frac{P_0 + \sqrt{D}}{Q_0} = \frac{(P_k + \sqrt{D}) A_{k-1} + Q_k A_{k-2}}{(P_k + \sqrt{D}) B_{k-1} + Q_k B_{k-2}}. \]

Clearing denominators on both sides and comparing rational and irrational parts, then substituting the expression for \(G_k\), we obtain the following equations:

\[ \begin{aligned} G_{k-1} &= P_k B_{k-1} + Q_k B_{k-2},\\ D B_{k-1} &= P_k G_{k-1} + Q_k G_{k-2}. \end{aligned} \]

Therefore, multiplying the first equation by \(G_{k-1}\) and subtracting the second equation multiplied by \(B_{k-1}\) yields

\[ \begin{aligned} G_{k-1}^2 - D B_{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_0 Q_k \\ &= (-1)^k Q_0 Q_k. \end{aligned} \]

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

To prove the second claim, substitute the expression for \(G_k\) into the first claim:

\[ (Q_0 A_{k-1} - P_0 B_{k-1})^2 - D B_{k-1}^2 = (-1)^k Q_0 Q_k. \]

So, using \(Q_0 \mid (P_0^2 - D)\),

\[ Q_0 A_{k-1}^2 + \left(\frac{P_0^2 - D}{Q_0} B_{k-1} - 2 P_0 A_{k-1}\right) B_{k-1} = (-1)^k Q_k. \]

Hence, \(\gcd(G_{k-1}, B_{k-1}) = \gcd(Q_0 A_{k-1}, B_{k-1})\) divides \(Q_k\).

This result provides a method to find solutions of \(x^2 - D y^2 = N\). If we choose \(Q_0 > 0\) appropriately and choose \(P_0\) as a solution to the congruence \(P_0^2 \equiv D \pmod{Q_0}\), then by running the PQa algorithm on \((P_0 + \sqrt{D})/Q_0\) until we find \((-1)^k Q_0 Q_k = N\), the pair \((G_{k-1}, B_{k-1})\) becomes a solution of the original equation. Moreover, if \(Q_k = \pm 1\), then the solution obtained is primitive, meaning \(G_{k-1}\) and \(B_{k-1}\) are coprime.

This idea is central to solving Pell's equation and the generalized Pell's equation. Understanding this, we now address some details of the algorithm and prove that all solutions can be obtained this way.

Pell's Equation

To solve Pell's equation \(x^2 - D y^2 = 1\), we run the PQa algorithm with \((P_0, Q_0, D) = (0, 1, D)\) until \((-1)^k Q_k = 1\). Then \((A_{k-1}, B_{k-1})\) is a solution of Pell's equation (since \(G_{k-1}\) is just \(A_{k-1}\) in this case). For Pell's equation, this process can be described more precisely.

First, the solution always appears at the end of a period. The above process is equivalent to the continued fraction expansion of \(\sqrt{D}\). For this, we already have the result:

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

Here, the period length is \(\ell\), and it starts at the first term (index starting from 0). Moreover, its \(\ell\)-th complete quotient equals \(\lfloor\sqrt{D}\rfloor + \sqrt{D}\), which implies \(Q_\ell = 1\). Therefore, if \(\ell\) is even, then \((A_{\ell-1}, B_{\ell-1})\) is a nontrivial solution of Pell's equation; if \(\ell\) is odd, then \((A_{2\ell-1}, B_{2\ell-1})\) is a nontrivial solution of Pell's equation.

Next, we show that the solution just obtained is indeed the fundamental solution. This is based on two reasons: First, all positive integer solutions \((x, y)\) of Pell's equation correspond to fractions \(x/y\) that appear as convergents of \(\sqrt{D}\), ensuring that \((x, y)\) must be some \((A_k, B_k)\) during the PQa algorithm; Second, aside from the end of the period, no other position will have \(Q_k = 1\), because the recurrence for \(A_k\) and \(B_k\) ensures they increase with index, so the smallest positive integer solution (the fundamental solution) must appear at the position indicated. These two reasons follow from the following two theorems:

Theorem

Suppose the equation \(x^2 - D y^2 = N\) has a positive integer solution \((x, y)\). If \(|N| < \sqrt{D}\), then \(\frac{x}{y}\) must be a convergent of \(\sqrt{D}\).

Proof

When \(N > 0\), since \(x^2 - D y^2 > 0\), we have \(x > y\sqrt{D}\). Hence,

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

By Legendre's criterion, \(\frac{x}{y}\) is a convergent of \(\sqrt{D}\).

When \(N < 0\), \(x > y\sqrt{D}\) no longer holds. Instead, consider the equation \(y^2 - \frac{1}{D} x^2 = -\frac{N}{D}\). Since \(\frac{|N|}{D} < \sqrt{\frac{1}{D}}\), the above argument still holds. This shows that \(\frac{y}{x}\) is a convergent of \(\frac{1}{\sqrt{D}}\). By the reciprocal theorem, \(\frac{x}{y}\) is also a convergent of \(\sqrt{D}\).

Theorem

In the PQa algorithm with \((P_0, Q_0, D) = (0, 1, D)\), if \(Q_k = 1\), then necessarily \(\ell \mid k\).

Proof

In the continued fraction expansion of \(\sqrt{D}\), all complete quotients except the 0-th are purely periodic. Suppose \(Q_k = 1\). By Galois' result, there must be a complete quotient \(\omega_k = P_k + \sqrt{D} > 1\) with its conjugate \(-1 < P_k - \sqrt{D} < 0\), implying \(P_k = \lfloor\sqrt{D}\rfloor\). Therefore, the complete quotient \(\omega_k\) equals \(\omega_\ell\). But repetition of complete quotients implies periodicity. If \(k\) is not a multiple of \(\ell\), it contradicts that \(\ell\) is the minimal positive period. Hence, \(\ell \mid k\).

Combining the discussions in this section, to solve Pell's equation, we perform the continued fraction expansion of \(\sqrt{D}\), i.e., run the PQa algorithm starting with \((P_0, Q_0, D) = (0, 1, D)\). When we first obtain \(Q_\ell = 1\), we have reached the end of the first period. If \(\ell\) is even, then \((A_{\ell-1}, B_{\ell-1})\) is the fundamental solution of Pell's equation; otherwise, \((A_{2\ell-1}, B_{2\ell-1})\) is the fundamental solution. For the case where the period length \(\ell\) is odd, we do not need to continue the PQa algorithm to twice the period; it will soon be shown that \(A_{2\ell-1} + B_{2\ell-1}\sqrt{D} = (A_{\ell-1} + B_{\ell-1}\sqrt{D})^2\), so we can directly compute the fundamental solution of Pell's equation from \((A_{\ell-1}, B_{\ell-1})\). All other solutions of Pell's equation can be computed from the fundamental solution.

Examples
  1. Solve \(x^2 - 14 y^2 = 1\).

    Running PQa with \((P_0, Q_0, D) = (0, 1, 14)\) yields: (red part indicates 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\)

    Period length \(\ell = 4\) is even. The minimal positive solution is \((G_3, B_3) = (15, 4)\).

  2. Solve \(x^2 - 41 y^2 = 1\).

    Running PQa with \((P_0, Q_0, D) = (0, 1, 41)\) yields: (red part indicates 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\)

    Period length \(\ell = 3\) is odd. The minimal positive 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's Equation

From the previous section, solutions of the negative Pell's equation must also correspond to convergents of \(\sqrt{D}\), and can only appear when \((-1)^k Q_k = -1\). This can only happen at the end of a period. Therefore, the negative Pell's 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 the fundamental solution. The solution method is the same as in the previous section.

Using a similar argument as in the proof of the structure of Pell's equation solutions, we can prove:

Theorem

Suppose the equation \(x^2 - D y^2 = -1\) has a solution, and let its fundamental solution be \((x_1, y_1)\). Then all integer solutions of \(x^2 - D y^2 = \pm 1\) belong to the set

\[ \{(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 exactly the fundamental solution of \(x^2 - D y^2 = 1\).

Proof

By symmetry, we only need to consider positive integer solutions, i.e., those with \(x + y\sqrt{D} > 1\). But because \(x^2 - D y^2 = \pm 1\) corresponds to two hyperbola branches, \(x + y\sqrt{D}\) cannot be put in one-to-one correspondence with \((x, y)\). To handle this difficulty, first prove that the above \((x_2, y_2)\) is the fundamental solution of \(x^2 - D y^2 = 1\).

Clearly, \((x_2, y_2)\) is a solution of \(x^2 - D y^2 = 1\). Let \((z, w)\) be the fundamental solution of \(x^2 - D y^2 = 1\). Then necessarily \(1 < z + w\sqrt{D} \le x_2 + y_2\sqrt{D}\). If the inequality on the right is strict, dividing by \(x_1 + y_1\sqrt{D}\) gives \(-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 term yields a form \(x' + y'\sqrt{D}\) with norm \(-1\) and \((x', y')\) also an integer solution. Taking the reciprocal of this inequality, we find 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, so one of them is greater than 1. But there should be no other quadratic integer with norm \(-1\) between 1 and \(x_1 + y_1\sqrt{D}\), contradicting the minimality of \(x_1 + y_1\sqrt{D}\). Hence, we must have \(x_2 + y_2\sqrt{D} = z + w\sqrt{D}\), i.e., \((x_2, y_2)\) is the fundamental solution of \(x^2 - D y^2 = 1\).

Based on this, if there exists a solution \((x, y)\) of \(x^2 - D y^2 = \pm 1\) not corresponding to some \((x_1 + y_1\sqrt{D})^k\), then there exists \(k\) such that \((x_1 + y_1\sqrt{D})^{2k} < x + y\sqrt{D} < (x_1 + y_1\sqrt{D})^{2k+2}\). Removing the factor \((x_1 + y_1\sqrt{D})^{2k+1}\) shows there exists a quadratic integer \(x' + y'\sqrt{D} \neq 1\) with norm \(\pm 1\) lying between \(-x_1 + y_1\sqrt{D}\) and \(x_1 + y_1\sqrt{D}\). Repeating the reciprocal argument as above leads to a contradiction with the minimality of \(x_1 + y_1\sqrt{D}\). Thus, the proposition is proved.

Since \((A_{\ell-1}, B_{\ell-1})\) is the minimal positive solution of the negative Pell's equation, and all positive integer solutions of \(x^2 - D y^2 = \pm 1\) appear in the set

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

and because these positive integer solutions must correspond to convergents of \(\sqrt{D}\) at the end of the period (previous position), and the numerators and denominators of convergents are strictly increasing, for all \(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 all these positive integer solutions, when \(k\) is odd, they are solutions of the negative Pell's equation; when \(k\) is even, they are solutions of Pell's equation, alternating.

Determining whether the negative Pell's equation has a solution requires computing the period length of the continued fraction expansion of \(\sqrt{D}\), which is not easy. Therefore, we hope to find simpler criteria. However, there is currently no simple, easily computable criterion5. Here we only provide a simple result.

Theorem

If the equation \(x^2 - D y^2 = -1\) has a solution, then \(4 \nmid 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, if the negative Pell's equation has a solution, then \(-1\) is a quadratic residue modulo \(D\), so \(-1\) is also a quadratic residue modulo any factor \(d\) of \(D\). Hence, \(d \neq 4\) and \(d\) is not a prime of the form \(4k+3\). Conversely, the equation \(x^2 - 2 y^2 = -1\) has the nontrivial solution \((1, 1)\). The remaining case is when \(D\) is a prime of the form \(4k+1\).

Let \(D\) be a prime of the form \(4k+1\). We want to prove that \(x^2 - D y^2 = -1\) has a solution. The idea is to start from the fundamental solution \((u, v)\) of Pell's equation \(x^2 - D y^2 = 1\) and construct a solution \((\alpha, \beta)\) of \(x^2 - D y^2 = -1\). If \(u\) is even, reducing \(u^2 - D v^2 = 1\) modulo 4 gives \(v^2 \equiv -1 \pmod 4\), but \(-1\) is not a quadratic residue modulo 4. This contradiction shows \(u\) is odd. Consider the equality \(D v^2 = u^2 - 1 = (u+1)(u-1)\). Since \(u\) is odd, \(\gcd(u+1, u-1) = \gcd(u+1, 2) = 2\). Using this fact, distributing the factors of \(D v^2\) to \(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 - D v^2 = 1\) yields \(\alpha^2 - D\beta^2 = \pm 1\). Since \((u, v)\) is the fundamental solution of Pell's equation and \((\alpha, \beta)\) is a smaller positive integer pair, the right-hand side cannot be \(+1\), so it must be \(-1\). This proves that \(x^2 - D y^2 = -1\) has a solution \((\alpha, \beta)\).

If \(D\) is composite, not having factors of the form \(4k+3\) and being squarefree does not guarantee that \(x^2 - D y^2 = -1\) has a solution. For example, \(x^2 - 34 y^2 = -1\) has no solution.

Examples

Using the calculations from the previous examples, the equation \(x^2 - 14 y^2 = -1\) has no solution, and the minimal positive solution of \(x^2 - 41 y^2 = -1\) is \((G_2, B_2) = (32, 5)\).

Norm ±4 Case

Next, we discuss solutions of \(x^2 - D y^2 = \pm 4\). The behavior of solutions depends on \(D \bmod 4\).

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

Now consider \(D \equiv 1 \pmod 4\), which cannot be simply reduced to already solved cases. To find the fundamental solution, we apply the PQa algorithm with \((P_0, Q_0, D) = (1, 2, D)\). When we first obtain \(Q_\ell = 2\), we have reached the end of the first period. If the period length \(\ell\) is even, then \((G_{\ell-1}, B_{\ell-1})\) is the fundamental solution of \(x^2 - D y^2 = 4\); otherwise, \((G_{\ell-1}, B_{\ell-1})\) is the fundamental solution of \(x^2 - D y^2 = -4\). From \((G_{\ell-1}, B_{\ell-1})\), we can obtain all solutions of \(x^2 - D y^2 = \pm 4\):

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

If the period length \(\ell\) is even, all these are solutions of \(x^2 - D y^2 = 4\); otherwise, when \(k\) is odd, \((x, y)\) is a solution of \(x^2 - D y^2 = -4\), and when \(k\) is even, \((x, y)\) is a solution of \(x^2 - D y^2 = 4\).

The correctness of this algorithm relies on the following facts:

Theorem

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

Proof

First note that in this case \(x\) and \(y\) have the same parity, so \((x+y)/2\) is an integer. If \((x, y)\) is a solution of \(x^2 - D y^2 = 4\), then \(x > y\sqrt{D} > 2y\), so

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

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

If \((x, y)\) is a solution of \(x^2 - D y^2 = -4\), to establish the above inequality, we need to prove \(4y < x + y\sqrt{D}\). This holds for all cases except possibly \(D = 5, 13\). For \(D = 5, 13\), substituting \(x = \sqrt{D y^2 - 4}\) into the inequality shows it is equivalent to \(2(\sqrt{D} - 2) y^2 > 1\). Except for \((D, y) = (5, 1)\), this holds for all \(D = 5, 13\) and positive integers \(y\). The remaining case is \((D, y) = (5, 1)\), where the solution of \(x^2 - 5 y^2 = -4\) is \((x, y) = (1, 1)\). We need to verify that \(\frac{1}{1}\) is a convergent of \(\frac{1+\sqrt{5}}{2} = [\overline{1}]\), which is obvious.

Theorem

Let \(D\) be a positive integer not a perfect square. The continued fraction expansion of the quadratic irrational \(\omega = \frac{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 all \(1 < k < \ell\).

Proof

Since \(\lfloor\omega\rfloor - 1 + \omega > 1\) and its conjugate \(\lfloor\omega\rfloor - \omega\) lies between \(-1\) and \(0\), by Galois' result, \(\lfloor\omega\rfloor - 1 + \omega\) is purely periodic and can be written as

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

Galois' result about the negative conjugate of the reciprocal gives

\[ \frac{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 + \frac{1}{\frac{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 implies \(a_k = a_{\ell-k}\) for all \(1 < k < \ell\), and the desired expansion follows.

Theorem

Let \(D \equiv 1 \pmod 4\). In the PQa algorithm with \((P_0, Q_0, D) = (1, 2, D)\), if \(Q_k = 2\), then necessarily \(\ell \mid k\).

Proof

In the continued fraction expansion of \(\frac{1+\sqrt{D}}{2}\), all complete quotients except the 0-th are purely periodic. Suppose \(Q_k = 2\). By Galois' result, there must be a complete quotient \(\omega_k = \frac{P_k + \sqrt{D}}{2}\) with conjugate \(-1 < \frac{P_k - \sqrt{D}}{2} < 0\), i.e., \(\sqrt{D} - 2 < P_k < \sqrt{D}\). Since in the PQa algorithm we always have \(Q_k \mid P_k^2 - D\) (see proof of algorithm correctness), \(P_k\) must be odd, which implies \(P_k\) is uniquely determined: \(P_k = P_0 + 2(\lfloor\omega\rfloor - 1)\). That is, the complete quotient \(\omega_k = \omega_\ell\). But repetition of complete quotients implies periodicity. If \(k\) is not a multiple of \(\ell\), it contradicts that \(\ell\) is the minimal positive period. Hence, \(\ell \mid k\).

Theorem

Let the minimal positive solution of \(x^2 - D y^2 = \pm 4\) be \((x_1, y_1)\). Then all its solutions are

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

By symmetry, we only need to consider positive integer solutions \((x, y)\). Here we only need to prove that the pairs \((x, y)\) in the set are indeed integer solutions of \(x^2 - D y^2 = \pm 4\). The rest follows by repeating the proof for the structure of solutions of \(x^2 - D y^2 = \pm 1\).

Actually, we need to prove that for any integer solutions \((x_1, y_1)\) and \((x_2, y_2)\) of \(x^2 - D y^2 = \pm 4\), the positive real pair \((x_3, y_3)\) defined as follows is also an integer solution:

\[ \frac{x_3 + y_3\sqrt{D}}{2} = \frac{x_1 + y_1\sqrt{D}}{2} \cdot \frac{x_2 + y_2\sqrt{D}}{2}. \]

Expanding the right-hand side and comparing rational and irrational parts gives

\[ x_3 = \frac{x_1 x_2 + D y_1 y_2}{2},\quad y_3 = \frac{x_1 y_2 + x_2 y_1}{2}. \]

Since for \(i = 1, 2\), we have \(x_i \equiv x_i^2 \equiv D y_i^2 \equiv D y_i \pmod 2\),

\[ \begin{aligned} 2x_3 &= x_1 x_2 + D y_1 y_2 \equiv D^2 y_1 y_2 + D y_1 y_2 = D(D+1) y_1 y_2 \equiv 0 \pmod 2,\\ 2y_3 &= x_1 y_2 + x_2 y_1 \equiv D y_1 y_2 + D y_2 y_1 = 2D y_1 y_2 \equiv 0 \pmod 2. \end{aligned} \]

This shows \(x_3\) and \(y_3\) are integers. Using the norm-preserving property of multiplication, \((x_3, y_3)\) is a solution of \(x^2 - D y^2 = \pm 4\).

Combining these facts and repeating the arguments from previous sections, we can justify the correctness of the above algorithm for solving \(x^2 - D y^2 = \pm 4\). These results show that \(x^2 - D y^2 = \pm 4\) has a simple solution structure similar to \(x^2 - D y^2 = \pm 1\): all its solutions can be expressed via its minimal positive solution without solving other equations.

In fact, all solutions of \(x^2 - D y^2 = \pm 1\) can be found among solutions of \(x^2 - D y^2 = \pm 4\), so from this perspective, \(x^2 - D y^2 = \pm 4\) is more fundamental. Clearly, \((x, y)\) is a solution of \(x^2 - D y^2 = \pm 1\) if and only if \((2x, 2y)\) is a solution of \(x^2 - D y^2 = \pm 4\). Previous analysis shows that when \(D \equiv 2, 3 \pmod 4\), all solutions of \(x^2 - D y^2 = \pm 4\) are both even, corresponding to solutions of \(x^2 - D y^2 = \pm 1\).

When \(D \equiv 0 \pmod 4\), in a solution \((x, y)\) of \(x^2 - D y^2 = \pm 4\), \(x\) is always even, but \(y\) may be odd. If in the minimal positive solution \((x_1, y_1)\) of \(x^2 - D y^2 = \pm 4\), \(y_1\) is even, then in all solutions \(y\) is also even, and these integer solutions correspond one-to-one with solutions of \(x^2 - D y^2 = \pm 1\). However, if in the minimal solution \((x_1, y_1)\), \(y_1\) is odd, then the parity of \(y_k\) alternates with \(k\), so only when \(k\) is even do they correspond to solutions of \(x^2 - D y^2 = \pm 1\). If the minimal positive solution of \(x^2 - D y^2 = \pm 4\) has odd \(y_1\) and the norm of \(x_1 + y_1\sqrt{D}\) is \(-4\), then for such \(D\), \(x^2 - D y^2 = -4\) has a solution, but \(x^2 - D y^2 = -1\) has no solution.

When \(D \equiv 1 \pmod 4\), solutions \((x, y)\) of \(x^2 - D y^2 = \pm 4\) may be both odd or both even. If the minimal positive solution \((x_1, y_1)\) is both even, then all its integer solutions are also both even, always corresponding to solutions of \(x^2 - D y^2 = \pm 1\). If the minimal positive solution \((x_1, y_1)\) is both odd, we have the following result:

Theorem

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

\[ \frac{x + y\sqrt{D}}{2} = \pm \left(\frac{x_1 + y_1\sqrt{D}}{2}\right)^{3k}, \quad k \in \mathbf Z. \]
Proof

Reducing \(x_1^2 - D y_1^2 = \pm 4\) modulo 8 gives \(D \equiv 5 \pmod 8\). To prove the second claim, first show that \((x_3, y_3)\) are both even, because

\[ \frac{x_3 + y_3\sqrt{D}}{2} = \left(\frac{x_1 + y_1\sqrt{D}}{2}\right)^3 = \frac{x_1^3 + 3D x_1 y_1^2}{8} + \frac{3x_1^2 y_1 + D y_1^3}{8}\sqrt{D}, \]

so we only need to prove the right-hand side is an integer. Since the square of an odd number is congruent to 1 modulo 8,

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

This shows \(x_3, y_3\) are both even. Then for all \(k \in \mathbf Z\),

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

so in this case \((x, y)\) are both even. Conversely, for \(r = 1, 2\), we always have

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

To prove the corresponding \((x, y)\) are not integers, it suffices to show the expression is not an integer, i.e., the second factor on the right is not an integer. For \(r = 1\), this is given; for \(r = 2\), because

\[ \frac{x_2 + y_2\sqrt{D}}{2} = \left(\frac{x_1 + y_1\sqrt{D}}{2}\right)^2 = \frac{x_1^2 + D y_1^2}{4} + \frac{x_1 y_1}{2}\sqrt{D}, \]

and \(x_1^2 + D y_1^2 \equiv 1 + 1 \times 1 = 2 \pmod 4\), \(x_1 y_1 \equiv 1 \pmod 2\), this expression is also not an integer. Hence, only when the exponent is a multiple of 3 are the corresponding solutions both even.

That is, every three solutions of \(x^2 - D y^2 = \pm 4\), one is both even, corresponding to an integer solution of \(x^2 - D y^2 = \pm 1\). This also shows that for \(D \equiv 1 \pmod 4\), the equation \(x^2 - D y^2 = -4\) has a solution if and only if \(x^2 - D y^2 = -1\) has a solution.

The discussion so far is sufficient to compute the fundamental unit of the ring of integers of a real quadratic field. Let \(D\) be a positive integer without square factors. For \(D \equiv 2, 3 \pmod 4\), we only need to find the minimal positive solution of \(x^2 - D y^2 = \pm 1\); for \(D \equiv 1 \pmod 4\), we only need to find the minimal positive solution of \(x^2 - D y^2 = \pm 4\). When we obtain the minimal positive solution \((x, y)\), for \(D \equiv 2, 3 \pmod 4\), the fundamental unit is \(\pm x \pm y\sqrt{D}\); for \(D \equiv 1 \pmod 4\), the fundamental unit is \(\frac{\pm x \pm y\sqrt{D}}{2}\).

Examples
  1. Solve \(x^2 - 14 y^2 = \pm 4\).

    Using previous calculations, the minimal positive solution of \(x^2 - 14 y^2 = 4\) is \((30, 8)\), and \(x^2 - 14 y^2 = -4\) has no solution.

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

    Running PQa with \((P_0, Q_0, D) = (1, 2, 41)\) yields: (red part indicates 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\)

    Period length \(\ell = 5\) is odd. The minimal positive solution of \(x^2 - 41 y^2 = -4\) is \((G_4, B_4) = (64, 10)\), and that of \(x^2 - 41 y^2 = 4\) is \((G_9, B_9) = (4098, 640)\). They are related by:

    \[ \frac{4098 + 640\sqrt{41}}{2} = \left(\frac{64 + 10\sqrt{41}}{2}\right)^2. \]

    Of course, since \(D \equiv 1 \pmod 8\), by the previous result, the minimal positive solutions of \(x^2 - 41 y^2 = \pm 4\) are both even and are exactly twice the minimal positive solutions of \(x^2 - 41 y^2 = \pm 1\), so we could also directly use the previous example.

  3. Solve \(x^2 - 13 y^2 = \pm 4\).

    Running PQa with \((P_0, Q_0, D) = (1, 2, 13)\) yields: (red part indicates 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\)

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

    Since the minimal positive solutions are both odd, we can use the pair at the end of the third period, \((G_2, B_2) = (36, 10)\), to obtain the minimal positive solutions of the (negative) Pell's equation \(x^2 - 13 y^2 = \pm 1\) as \((18, 5)\). This can also be computed directly:

    \[ \frac{36 + 10\sqrt{13}}{2} = \left(\frac{3 + \sqrt{13}}{2}\right)^3. \]

    And this is a solution of the negative Pell's equation. The corresponding minimal positive solution of Pell's equation is \((649, 180)\).

  4. Solve \(x^2 - 52 y^2 = \pm 4\).

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

General Case

Finally, we discuss solving the generalized Pell's equation.

For \(|N| < \sqrt{D}\), there is a simple method. The previous result shows that a solution \((x, y)\) of \(x^2 - D y^2 = N\) must satisfy \(\frac{x}{y}\) being a convergent of \(\sqrt{D}\). Moreover, from the solution structure discussed earlier, each fundamental solution \((x, y)\) satisfies \(x + y\sqrt{D} \le x_1 + y_1\sqrt{D}\), where \(x_1 + y_1\sqrt{D}\) is the fundamental solution of the corresponding Pell's equation \(x^2 - D y^2 = 1\). Using the monotonicity of the denominator sequence \(B_k\) in the PQa algorithm, these fundamental solutions of the generalized Pell's equation must appear before the fundamental solution of the corresponding Pell's equation appears. Thus, we run the PQa algorithm with \((P_0, Q_0, D) = (0, 1, D)\) until \(Q_{\ell'} = 1\) and \(\ell'\) is even. During the process, for each \((A_k, B_k)\) we check whether there exists an integer \(f\) such that

\[ A_k^2 - D B_k^2 = (-1)^{k+1} Q_{k+1} = N / f^2 \]

holds. If so, record \((f A_k, f B_k)\) as a minimal positive solution. All recorded \((f A_k, f B_k)\) are all minimal positive solutions of \(x^2 - D y^2 = N\). Using \((A_{\ell'-1}, B_{\ell'-1})\), the fundamental solution of the corresponding Pell's equation, we can generate all solutions of the generalized Pell's equation from these minimal positive solutions. Note that depending on whether the period length \(\ell\) is even or odd, \(\ell'\) may be \(\ell\) or \(2\ell\).

For more general \(N\), the above method no longer applies. First, enumerate all square factors \(f^2\) of \(N\), let \(m = N / f^2\), and enumerate all solutions \(z\) of the congruence \(z^2 \equiv D \pmod{|m|}\) satisfying \(-|m|/2 < z \le |m|/2\). Then run the PQa algorithm with \((P_0, Q_0, D) = (z, |m|, D)\) until \(Q_k = \pm 1\) or a full period has been completed. In the second case, no solution related to that \((f, z)\) exists. In the first case, we further check whether \((-1)^k Q_k = N / |N|\). If the signs match, then \((f G_{k-1}, f B_{k-1})\) is a solution of \(x^2 - D y^2 = N\). Otherwise, it is a solution of \(x^2 - D y^2 = -N\), and only if the corresponding negative Pell's equation has a solution can we combine it with the fundamental solution of the negative Pell's equation to obtain a solution of \(x^2 - D y^2 = N\). After iterating over all pairs \((f, z)\), we obtain exactly one solution from each equivalence class of solutions of \(x^2 - D y^2 = N\), and this solution is the fundamental or minimal positive solution of that class. Using them and the fundamental solution of the corresponding Pell's equation, we can generate all integer solutions of the equation. This algorithm is called the Lagrange–Matthews–Mollin algorithm.

The correctness of this algorithm is guaranteed by the following theorem:

Theorem

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

Proof

Using \(x \equiv -P_0 y \pmod{Q_0}\) and \(x^2 - D y^2 = N \equiv 0 \pmod{Q_0}\), clearly \(P_0^2 \equiv D \pmod{Q_0}\). Thus,

\[ P_0 x + D y \equiv -P_0^2 y + D y = (D - P_0^2) y \equiv 0 \pmod{Q_0}. \]

Hence, consider the integer matrix

\[ \begin{pmatrix}P & R \\ Q & S\end{pmatrix} = \begin{pmatrix}X & \dfrac{P_0 x + D y}{Q_0} \\ y & x\end{pmatrix}. \]

Its determinant

\[ PS - QR = \frac{x(x + P_0 y) - y(P_0 x + D y)}{Q_0} = \frac{x^2 - D y^2}{Q_0} = \pm 1. \]

Moreover, let \(\zeta = \sqrt{D} > 1\). Then

\[ \frac{P\zeta + R}{Q\zeta + S} = \frac{(x + P_0 y)\sqrt{D} + (P_0 x + D y)}{(x + y\sqrt{D}) Q_0} = \frac{P_0 + \sqrt{D}}{Q_0} = \omega. \]

We now prove that \(\frac{P}{Q}\) is a convergent of \(\omega\). Suppose \(\frac{P}{Q}\) has the simple continued fraction expansion

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

and \(PS - QR = (-1)^{k-1}\). If \(\frac{p_k}{q_k}\) is its \(k\)-th convergent, then \((p_k, q_k) = (P, Q)\), and by the difference formula for convergents, \(p_k q_{k-1} - q_k p_{k-1} = (-1)^{k-1}\). This implies

\[ p_k (S - q_{k-1}) = q_k (R - p_{k-1}). \]

Consider cases:

  • If \(S = 0\), then easily \(Q = R = 1\), so \(\omega = P + \zeta^{-1} = [P, \zeta]\), and \(\frac{P}{Q} = P\) is the 0-th convergent of \(\omega\).
  • If \(Q = S > 0\), then \(Q = S = 1\) and \(P - R = \pm 1\).
    • If \(P = R + 1\), then \(\omega = R + \frac{1}{1 + \zeta^{-1}} = [R, 1, \zeta]\), so \(\frac{P}{Q} = \frac{R+1}{1} = [R, 1]\) is the 1st convergent of \(\omega\).
    • If \(P = R - 1\), then \(\omega = R - 1 + \frac{1}{1 + \zeta} = [R-1, \zeta-1]\), so \(\frac{P}{Q} = R-1\) is the 0-th convergent of \(\omega\).
  • If \(Q \neq S > 0\), then since \(Q = q_k \mid (S - q_{k-1})\), there exists an integer \(\kappa\) such that \(S = \kappa q_k + q_{k-1}\) and \(R = \kappa p_k + p_{k-1}\). Because \(q_k \ge q_{k-1}\) and \(S > 0\), \(\kappa \ge 0\). Thus, \(\omega = \frac{(\kappa + \zeta) p_k + p_{k-1}}{(\kappa + \zeta) q_k + q_{k-1}} = [a_0, a_1, \cdots, a_k, \kappa + \zeta]\), so \(\frac{P}{Q}\) is its \(k\)-th convergent.

In summary, \(\frac{X}{y}\) is always a convergent of \(\omega = \frac{P_0 + \sqrt{D}}{Q_0}\), denoted as \(\frac{A_{k-1}}{B_{k-1}}\) in PQa notation. Since \(A_{k-1}^2 - D B_{k-1}^2 = (-1)^k Q_0 Q_k\), we have \(Q_k = (-1)^k \frac{N}{|N|}\).

This theorem ensures that all positive solutions of the equation exist among the convergents of the corresponding quadratic irrational. Since when computing convergents via the PQa algorithm, once entering the period, the convergents are always positive. So by enumerating all quadratic irrationals allowed by the theorem conditions and computing their convergents up to one period, we can find a solution. Because two solutions appearing in the convergents of the same quadratic irrational are necessarily equivalent, we can stop further computation once we find the first solution satisfying \((-1)^k Q_k = N/|N|\). Unlike all previous algorithms, here the satisfying \(k\) may appear before entering the period.

Examples
  1. Solve \(x^2 - 157 y^2 = 12\).

    Since \(12^2 < 157\), run PQa with \((P_0, Q_0, D) = (0, 1, 157)\) yields: (red part indicates 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\)

    Period length \(\ell = 17\) is odd, so we need to examine two periods for cases where \(G_{k-1}^2 - 157 B_{k-1}^2\) differs from 12 by a square factor, i.e., \(k = 1, 9, 13, 19, 23, 31\). The corresponding solutions are \((fG, fB)\) in the table:

    \(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 solutions in each equivalence class of \(x^2-157y^2=12\). To get all solutions, use the fundamental solution of Pell's equation \((46698728731849,3726964292220)\). For example, transform them into the fundamental solution \((x,y)\) of that class, also listed in the table.

  2. Solve \(x^2-157y^2=12\) using Lagrange–Matthews–Mollin.

    First, enumerate square factors of \(N=12\):

    • \(f^2=1^2\): \(m=12\), congruence \(P^2\equiv 157\pmod{12}\) has solutions \(z=\pm 1,\pm 5\);
    • \(f^2=2^2\): \(m=3\), congruence \(P^2\equiv 157\pmod{3}\) has solutions \(z=\pm 1\).

    For each \((f,z)\) run PQa with \((P_0,Q_0,D)=(z,|m|,D)\) and find first \((-1)^kQ_k=1\). The corresponding \((fG_{k-1},fB_{k-1})\) is a solution. Results:

    \(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 all minimal positive solutions in each equivalence class, which can be transformed to fundamental solutions using Pell's fundamental solution.

  3. Solve \(x^2-79y^2=\pm 101\).

    Use Lagrange–Matthews–Mollin. Since \(N=101\) is prime, \(f=1\). Then \(m=101\), congruence \(P^2\equiv 79\pmod{101}\) has solutions \(P=\pm 33\).

    Run PQa with \((P_0,Q_0,D)=(33,101,79)\) yields: (red part indicates 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\)

    Period length \(\ell=6\) even. Up to one period, no \(Q_k=\pm 1\), so no solution. Similarly for \((-33,101,79)\). Hence, no solution.

Exercises

References and Notes


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

  2. Some Chinese literature also calls it the second-type Pell's equation. ↩

  3. i.e., rational numbers of the form \(n+\frac12\) with \(n\in\mathbf Z\). ↩

  4. Note that the definition of fundamental solution in Pell's equation is not consistent with the definition of fundamental unit in the ring of integers of a real quadratic field. First, in some rings, the fundamental unit \(x+y\sqrt{D}\) has \(x,y\) half-integers, so it is not a solution of Pell's equation. Second, there are four fundamental units in the same ring, but only one fundamental solution, because the fundamental solution requires \(x,y\) both positive. ↩

  5. A more practical criterion and tool can be found here and its references. The list of positive integers \(D\) for which \(x^2-Dy^2=-1\) has a solution is OEIS A031396. ↩