In this post, we are interested in finding all integer solutions of Pell equation

$$\begin{equation} x^2 - d y^2 = 1, \label{eq:pell} \end{equation}$$

and generalized Pell equation

$$\begin{equation} x^2 - d y^2 = n, \label{eq:generalized-pell} \end{equation}$$

for a non-square integer $d > 0$ and $n \neq 0$.

Connection to $\ZZ[\sqrt{ d }]$

For any element $s = x + y \sqrt{ d }$ in $\ZZ[\sqrt{ d }]$, we define $\bar{s} := x - y \sqrt{ d }$. It is easy to show that $\bar{s} \bar{t} = \overline{st}$ for any $s, t \in \ZZ[\sqrt{d}]$.

Note that for $s = x + y \sqrt{ d }$, $x^2 - dy^2 = s \bar{s}.$ Now we define

$$S_{n} := \{ s \in \ZZ[\sqrt{ d }]: s \bar{s} = n \}.$$

so solving Pell equations is equivalent to calculating $S_1$ and solving generalized Pell equations is equivalent to calculating $S_n$.

Most references I found online use the $(x, y)$ notation, but in certain cases, it might be cleaner to present the results using operations in $\ZZ[\sqrt{ d }]$, especially when presenting the structure of solutions.

Primitive solutions

(Primitive solutions) A solution $(x, y)$ is primitive if $\gcd(x, y) = 1$. If $(x, y)$ is primitive, $\gcd(x, n) = \gcd(y, n) = 1$.

Similar to solving $ax^2 + by^2 = n$, once we have an algorithm to find all primitive solutions, we can derive all solutions. Let $S'_n$ be the set of all primitive solutions of $x^2 - d y^2 = n$:

$$S'_{n}:= \left\{ s\in \ZZ[\sqrt{ d }]: s \bar{s} = n, \quad s = x + y \sqrt{ d }, \quad \gcd(x, y) = 1 \right\},$$

then it is straightforward to see that

$$S_{n} = \left\{ gs: g \in \NN, \quad g^2 | n, \quad s \in S'_{n / g^2} \right\} .$$

Clearly, $S_1 = S'_1$.

Pell equation $x^2 - d y^2 = 1$

The Pell equation has been studied extensively over time. See [Lenstra02] for a comprehensive overview. There are so many tutorials (e.g. [Conrad22A] and [Conrad22B]) on solving Pell equations, so here I only list a few key conclusions and the algorithm.

(Structure of solutions) It follows from Dirichlet’s unit theorem that $S_{1}$ has the structure

$$S_{1} = \{ \pm s^{*k}: k \in \ZZ \},$$

for a unique fundamental solution $s^* := x^* + y^* \sqrt{d}$ where $x^* > 0, y^* > 0$. Note $n$ can be negative, as $\bar{s} = s^{*-1}$. See [Conrad22A, Theorem 5.3] for proof.

(Fundamental solution). It is also shown that $x^*/y^*$ must be a convergent of $\sqrt{d}$. We can use the PQa algorithm below to find all convergents. See [Conrad22B, Theorem 5.1] for proof.

(PQa algorithm) Given $p_0, q_0$ and $d$ such that $p_0^2 \equiv d \pmod {q_{0}}$, the PQa algorithm computes the continued fraction $[ a_0, a_1, a_2, \dots]$ of $(p_0 + \sqrt{d}) / q_0$. It works as follows:

$$\begin{aligned} a_i & = \floor{(p_i + \sqrt{d}) / q_i}, \\ p_{i+1} & = a_i q_i - p_i, \\ q_{i+1} & = (d -p_{i+1}^2)/q_i. \\ \end{aligned}$$

Note that the division for calculating $q_{i+1}$ is always exact.

Below is an implementation of the PQa algorithm optimized for solving the Pell equation. See [Robertson04] for more details.

def pell(d):
    """find the fundamental solution of x^2 - dy^2 = 1 and -1"""
    p, q = 0, 1

    sqrt_d = int(d ** 0.5)
    x0, x1 = 0, 1
    y0, y1 = 1, 0

    l = 0
    while l == 0 or q != 1:
        a = (p + sqrt_d) // q
        p = a * q - p
        q = (d - p ** 2) // q
        l += 1

        x0, x1 = x1, a * x1 + x0
        y0, y1 = y1, a * y1 + y0

    if l % 2 == 1:
        pos = x1 * x1 + y1 * y1 * d, x1 * y1 * 2
        neg = x1, y1
    else:
        pos = x1, y1
        neg = None

    return pos, neg

Generalized Pell equation $x^2 -d y^2 = n$

(The idea in this section is taken from [Matthews00].)

(Nagell equivalence) For any $s \in S'_{n}$ and $t \in S'_1$, we observe

$$(s t) \overline{(st)} = (s \overline{s})(t \overline{t}) = n, \quad (st) \bar{t} = s \in S'_{n},$$

implying that $s t \in S'_{n}$. This motivates us to define equivalence relation, the Nagell equivalence, for elements in $S'_{n}$: Two solutions $s \sim s'$ if and only if there exists an element $t \in S'_{1}$, such that $st = s'$.

(Partition of solution space) Let $s = x + y \sqrt{ d }$ and $s' = x' + y' \sqrt{ d }$ be two primitive solutions, and define $\lambda(s) := -x y^{-1} \bmod |n|$. Since $\lambda(s)^2 \equiv d \pmod {|n|}$, it follows that

$$\begin{aligned} s \sim s' \quad & \Longleftrightarrow \quad s s'^{-1} \in \ZZ[\sqrt{ d }] \\ & \Longleftrightarrow \quad x x' \equiv yy' d \pmod{|n|}, \quad xy' \equiv x' y \pmod {|n|} \\ & \Longleftrightarrow \quad \tau(s) = \tau(s'). \end{aligned}$$

So an alternative definition of Nagall equivalence for primitive solutions is that $s \sim s'$ if and only if $\lambda(s) = \lambda(s')$.

(Structure of solutions) The alternative definition of Nagall equivalence directly implies the number of equivalent classes in $S'_n$ is finite, that is, there exists a finite set $F_{n} \subset \ZZ[\sqrt{ d }]$ such that

$$S'_{n} = \left\{f t: f \in F_{n}, t \in S'_{1} \right\} = \left\{ \pm f s^{*n}: F \in F_{n, }n \in \ZZ \right\} ,$$

where $s^*$ is the fundamental solution of $x^2 - dy^2 = 1$.

(Finding a minimum positive solution) Now we focus on one equivalence class, represented by $\lambda \equiv -x y^{-1} \pmod{|n|}$ (with abuse of notation) where $-|n| / 2 < \lambda \leq |n| / 2$. If we write $x = x' |n| - \lambda y$, it is proved in [Matthews00, Theorem 1] that $x' / y$ is a convergent of $(\lambda + \sqrt{ d }) / |n|$. Thus, we can find a minimum positive solution $s = x + y \sqrt{ d }$ with smallest positive $x$ and $y$ using the PQa algorithm.

(Finding a fundamental solution) It might be more convenient to find a fundamental solution for an equivalence class, with smallest positive $y$. If multiple solutions have the same smallest positive $y$, choose the one with positive $x$. A fundamental solution can be found as follows:

• If $\lambda \in \{ 0, n / 2\}$, the minimum positive solution is the fundamental solution;
• Otherwise, let $(x, y)$ and $(x', y')$ be the minimum positive solutions of $\lambda$ and $-\lambda$. If $y < y'$, the fundamental solutions are $(x, y)$ and $(-x, y)$; otherwise they’re $(-x', y')$ and $(x', y')$.

The algorithm is almost the same as [Mollin01, Algorithm 4.1]. The Lagrange-Matthews-Mollin (LMM) algorithm, proposed in [Matthews00, Section 5], improves the process by considering properties of the continued fraction of $(\lambda + \sqrt{ d }) / |n|$.

Here is the pseudocode:

def solve_generalized_pell_primitive(d, n):
    """find all fundamental solutions of x^2 - d y^2 = n"""

    for λ in sqrtmod_all(d, n):
        if (s := PQa(d, λ, n)) is None:
            continue

		if λ == 0 or λ * 2 == n:
	        yield s
        elif λ * 2 < n:
            p1, q1 = s
            p2, q2 = PQa(d, -λ, n)

            if q1 < q2:
                yield from (p1, q1), (-p1, q1)
            else:
                yield from (p2, q2), (-p2, q2)

A special generalized Pell equation $x^2 - dy^2 = 4$

We also discuss special generalized Pell equation $x^2 - dy^2 = 4$, as it also shares many properties as the Pell equation $x^2 - dy^2 = 1$.

As we have mentioned above, the set $S_1$ has a special structure $S_1 = \{s_{1}^{*n}: n \in \ZZ\}$ for a fundamental solution $s_{1}^*$. Similarly, the solutions for $x^2 - dy^2 = 4$ have a similar structure

$$T_{4} = \left\{ \pm t_{4}^{*k}: k \in \ZZ \right\}, \quad T_{4} := \left\{ \frac{x+y\sqrt{ d }}{2} : x^2 - dy^2 = 4 \right\}.$$

for a fundamental solution $t_4^* = \frac{x^* + y^* \sqrt{ d }}{2}$ where $(x^*, y^*)$ is minimum positive solution to $x^2 - dy^2 = 4$ [Stolt57, p.4]. There is actually a connection between $s_{1}^{*}$ and $t_4^{*}$: $s_1^{*} = t_{4}^{*k}$ for $k \in \left\{ 1,2,3 \right\}$ ([Stolt57, p.4], [Matthews23, Theorem 1.1]).

(Equivalence relation) Observing that Negall equivalence is the associated equivalence relation of the group action of group $\{\pm s_1^{*k}: k \in \ZZ\}$ on the set $S'_n$. As $T_4$ is also a group, we can define a similar equivalence relation, the Stolt equivalence. In some sense, the Stolt equivalence is more fundamental than the Negall equivalence, in the sense that $s_1^*$ can be generated by $t_4^*$ but $t_4^*$ might not be generated by $s_1^*$. As a consequence, the Stolt equivalence might have fewer equivalence classes than the Negall equivalence.

To find the fundamental solution $t_4^*$, check out [Johnson04] for an efficient algorithm implemented below:

def pell4(d):
    """find the mininum positive solution to x^2 - d y^2 = 4 or -4"""
    if d % 4 == 1:
        p, q = 1, 2

        sqrt_d = int(d ** 0.5)
        x0, x1 = -1, 2
        y0, y1 = 1, 0

        l = 0
        while l == 0 or q != 2:
            a = (p + sqrt_d) // q
            p = a * q - p
            q = (d - p ** 2) // q
            l += 1

            x0, x1 = x1, a * x1 + x0
            y0, y1 = y1, a * y1 + y0

        if l % 2 == 1:
            pos = (x1 * x1 + y1 * y1 * d) // 2, x1 * y1
            neg = x1, y1
        else:
            pos = x1, y1
            neg = None

    elif d % 4 == 0:
        pos, neg = pell(d // 4)
        pos = (pos[0] * 2, pos[1])
        if neg:
            neg = (neg[0] * 2, neg[1])

	else:
        pos, neg = pell(d)
        pos = (pos[0] * 2, pos[1] * 2)
        if neg:
            neg = (neg[0] * 2, neg[1] * 2)

    return pos, neg

References

• Rathnayake13: Rathnayake T., 2013. Quadratic Diophantine equation – I [online]

• Lenstra02: Lenstra Jr, H.W., 2002. Solving the Pell equation. Notices of the AMS, 49(2), pp.182-192.

• Conrad22A: Conrad, K., 2022. Pell’s equation, I [online]

• Conrad22B: Conrad, K., 2022. Pell’s equation, II [online]

• Matthews00: Matthews, K., 2000. The Diophantine Equation $x^ 2-Dy^ 2= N$, $D> 0$. Expositiones Mathematicae, 18(4), pp.323-332.

  • slides.

• Tamang21: Tamang, B.B., 2021. On the study of quadratic Diophantine equations (Doctoral dissertation, Tribhuvan University Kathmandu).

  • Discussed this topic in more details and shows more examples.

• Mollin01: Mollin, R.A., 2001. Simple continued fraction solutions for Diophantine equations. Expositiones Mathematicae, 19(1), pp.55-73.

• Robertson04: Robertson, J.P., 2004. Solving the generalized Pell equation $x^2− Dy^2= N$. Unpublished manuscript.

  • It gives a complete description of the LMM algorithm.

• Stolt57: Stolt, B., 1957. On a Diophantine equation of the second degree. Arkiv för Matematik, 3(4), pp.381-390.

---

Update @2025-03-20: Added the discussion of $x^2 - dy^2 = 4$ and the optimized algorithm for the Pell equation.