In this post, we’re interested in finding all integer solutions of a binary quadratic equation:

$$\begin{equation} ax^2 + by^2 = n, \label{eq:ellipse} \end{equation}$$

for $a, b, n > 0$.

Normalization

(Pairwise coprime) Equation \eqref{eq

} is reduced if $a, b, n$ are pairwise coprime. We can reduce it as follows:

• If $p | a$ and $p | b$, then $p$ must divide $n$ (otherwise there is no solution), so we can reduce it to $(a / p, b / p, n / p)$.
• If $p | a$ and $p | n$ but $p \nmid b$, then $p | y$, so we can reduce it to $(a / p, yp, n / p)$.
• Similarly for $p | b$ and $p | n$ but $p \nmid a$.

Every time we make a reduction, $|abn|$ decreases so it must terminate with a reduced equation.

(Primitive solutions) A solution $(x, y)$ is primitive if and only if $\gcd(x, y) = 1$. It’s trivial to use an algorithm which only finds all primitive solutions to find all solutions: If $g := \gcd(x, y) > 1$, then $(x / g, y / g)$ is a primitive solution to the equation $ax^2 + by^2 = n/g^2$.

From this point forward, we assume that Equation \eqref{eq

} is reduced unless stated otherwise, and we focus on the primitive solutions:

Method 1: Lattice

(Partition of solution space) We take the idea from a method on solving Legendre equations. Let $\mathcal{S}$ be a set of integer points in $\ZZ^2$:

$$S = \{(x, y): a x^2 + by^2 \equiv 0 \pmod n, \quad \gcd(x, y) = 1\}.$$

Obviously, $S$ contains all solutions to Equation \eqref{eq

}. We may partition $S$ into subsets: Let $\Lambda_{g} = \{\lambda: \lambda^2 \equiv -ab \pmod{n}\}$, then

$$S = \bigcup_{\lambda \in \Lambda} S_{\lambda}, \quad S_{\lambda} := \{(x, y): a x \equiv \lambda y \pmod{n}, \quad \gcd(x, y) = 1\}.$$

Next, we need to determine the solutions in each individual $S_{\lambda}$.

(Lattice) To find the solutions in $S_{\lambda}$, we do it in a superset of $S_{\lambda}$, denoted as $L_{\lambda}$, and check if the solution is in $L_{\lambda}$. The superset $L_{\lambda}$ is defined as:

$$L_{\lambda} := \{(x, y): a x \equiv \lambda y \pmod{n}\}.$$

It turns out that $L_{\lambda}$ is a lattice with basis $\mathbf{u} = (n, 0)$, $\mathbf{v} = (\lambda a^{-1} \bmod{n}, 1)$. It’s also easy to check that every point $(x, y)$ in $L_{\lambda}$ satisfies $ax^2 + by^2 \equiv 0 \pmod n$, thus, checking all minimizers of $ax^2 + by^2$ in $L_{\lambda}$ suffices to find all solutions to Equation \eqref{eq

}.

(Shortest vector) As in solving Legendre equation, we can use the same lemma for 3D lattice to efficiently find the shortest vector after reducing the basis by either Lenstra–Lenstra–Lovász algorithm or Lagrange-Gauss algorithm.

The following is the pseudocode for the complete algorithm. We need an efficient algorithm to solve $\lambda^2 \equiv -ab \pmod n$, which was covered before. It’s also interesting to explore whether it is possible to avoid repeatly solving $\lambda$‘s but that is left for future work. (Well, probably there won’t be any future work because the other method is much easier…)

def solve_primitive(a, b, n):
    """find all *primitive* solutions of (x, y) such that a x^2 + b y^2 = n"""
    assert gcd(a, b) == 1 and gcd(b, n) == 1 and gcd(a, n) == 1
    assert 1 <= a and 1 <= b and 1 <= n

    dot = lambda u, v: sum(u * v * [a, b])
    for λ in sqrtmod_all(-a * b, n):
        B = np.array([
            [n, 0],
            [λ * invmod(a, n) % n, 1],
        ])
        B = LLL(B, dot=dot)
        for w in product([-1, 0, 1], [-1, 0, 1]):
            if w != (0, 0) and gcd(x := w @ B) == 1 and dot(x, x) == n:
                 yield x

Method 2: The Hardy-Muskat-Williams algorithm

The most relevant paper is [HMW90], which introduces the Hardy-Muskat-Williams (HMW) algorithm. For pairwise copime integers $a, b, n > 0$ such that $n > a + b$, it finds all primitive solutions $(x, y)$ to $ax^2 + by^2 = n$ with $x \geq 1, y \geq 1$.

(Parition of solutions) Similarly, the HMW algorithm partitions the solutions using a linear congruence: For any $0 < \lambda < n / 2$ such that $a \lambda^2 + b \equiv 0 \pmod n$, define

$$S_\lambda = \{(x, y) \in \NN^{+2}: a x^2 + b y^2 = n, \gcd(x, y) = 1, x \equiv \pm \lambda y \pmod n\},$$

where pairs are unordered when $ab = 1$, that is, $(x, y)$ and $(y, x)$ is considered the same solution when $ab = 1$. It can be shown that $|S_\lambda| \leq 1$.

_(Structure of $S_\lambda$)_ The HMW algorithm is based on the following elegant theorem:

The implementation is surprisingly easy:

def solve_primitive(a, b, n):
    """find all positive primitive (unordered if a = b = 1) solutions such that ax^2 + by^2 = n"""
    for t in sqrtmod_all(-b * invmod(a, n) % n):
        if t >= n // 2: continue
        p, x = n, t
        while a * x * x >= n:
            p, x = x, p % x
        y = int(sqrt((n - a * x * x) // b))
        if a * x * x + b * y * y == n:
            yield x, y

A special case of the HMW algorithm is Cornacchia’s algorithm where $a = 1$. [Basilla04] proved the correctness of Cornacchia’s algorithm using same lattice as ours.

A failed attempt

The Brahmagupta–Fibonacci identity states that

$$(x_1 + d y_1^2)(x_2 + d y_2^2) = (x_1 x_2 - n y_1 y_2)^2 + n (x_1 y_2 + x_2 y_1)^2,$$

so the attempt is to find all prime factors of $n$ and solve subproblems for each prime factors. However, not all solutions can be generated this way. A simple counterexample is $1^2 + 5 \times 1^2 = 6$, while $x^2 + 5 y^2 = 2$ and $x^2 + 5 y^2 = 3$ have no solutions.

It’s also related to the fact that $\ZZ[\sqrt{ -D }]$ is not a UFD for $D \geq 3$.

Reference

Due to an incomplete literature review, I initially overlooked a significant portion of previous work before arriving with the lattice algorithm. Later I came across Cornacchia’s algorithm algorithm, which can be further generalized to the HMW algorithm.

• Rathnayake13: Rathnayake T., 2013. Improving the solving process using Cornacchia [online]
• Basilla04: Basilla, J.M., 2004. On the solution of $x^2+dy^2=m$.
• Lattice Reduction: https://en.wikipedia.org/wiki/Lattice_reduction
• Deng16: Deng, X., 2016. An introduction to lenstra-lenstra-lovasz lattice basis reduction algorithm. Massachusetts Institute of Technology (MIT).
• Nguyen09: Nguyen, P.Q., 2009. Hermite’s constant and lattice algorithms. In The LLL Algorithm: Survey and Applications (pp. 19-69). Berlin, Heidelberg: Springer Berlin Heidelberg.
• HMW90: Hardy, K., Muskat, J.B. and Williams, K.S., 1990. A deterministic algorithm for solving $n = fu^2 + gv^2$ in coprime integers 𝑢 and 𝑣. Mathematics of Computation, 55(191), pp.327-343.

---

Update @2025-03-13: Added the Hardy-Muskat-Williams algorithm.