Integer and Rational Solutions of a Binary Quadratic Equation (4): Pell equation and generalized Pell equation

Posted on In , Disqus: Word count in article: 1.8k Reading time ≈ 6 mins.

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} \tag{Pell Equation} \end{equation} \] and generalized Pell equation \[ \begin{equation} x^2 - d y^2 = n, \label{eq:generalized-pell} \tag{Generalized Pell Equation} \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.