Integer and Rational Solutions of a Binary Quadratic Equation (3): $ax^2+by^2=n$

Posted on 2025-03-11 In CS/Math , Number Theory Disqus: Word count in article: 913 Reading time ≈ 3 mins.

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

Equation $\eqref{eq:ellipse}$ 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.

From this point forward, we assume that Equation $\eqref{eq:ellipse}$ is reduced unless stated otherwise.

Solve a reduced equation

(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$: \[ \mathcal{S} = \{(x, y): a x^2 + by^2 \equiv 0 \pmod n\}. \] Obviously, $\mathcal{S}$ contains all solutions to Equation $\eqref{eq:ellipse}$. We may partition $\mathcal{S}$ into subsets: For any $g > 0$ such that $g^2 | n$, let $\Lambda_{g} = \{\lambda: \lambda^2 \equiv -ab \pmod{n/g^2}\}$, then \[ \mathcal{S} = \bigcup_{\substack{g^2 | n \\ \lambda \in \Lambda}} \mathcal{S}_{g, \lambda}, \quad \mathcal{S}_{g,\lambda} := \{(g x, g y): \gcd(x, y) = 1, a x \equiv \lambda y \pmod{n/g^2}\}. \] Next, we need to determine the solutions in each individual $\mathcal{L}_{g,\lambda}$.

(Lattice) To find the solutions in $\mathcal{S}_{g, \lambda}$, we do it in a superset of $\mathcal{S}_{g, \lambda}$, denoted as $\mathcal{L}_{g, \lambda}$, and check if the solution is in $\mathcal{L}_{g, \lambda}$. The superset $\mathcal{L}_{g, \lambda}$ is defined as: \[ \mathcal{L}_{g,\lambda} := \{(g x, g y): a x \equiv \lambda y \pmod{n/g^2}\}. \] It turns out that $\mathcal{L}_{g, \lambda}$ is a lattice with basis $\mathbf{u} = (n/g, 0)$, $\mathbf{v} = (g (\lambda a^{-1} \bmod{n / g^2}), g)$. It's also easy to check that every point $(x, y)$ in $\mathcal{L}_{g, \lambda}$ satisfies $ax^2 + by^2 \equiv 0 \pmod n$, thus, checking all minimizers of $ax^2 + by^2$ in $\mathcal{L}_{g, \lambda}$ suffices to fine all solutions to Equation $\eqref{eq:ellipse}$.

(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.

Let $\mathcal{L} = \{\mathbf{w} \mathbf{B}: \mathbf{w} \in \ZZ^2\}$ be a 2D lattice where the (row-based) basis $\mathbf{B}$ is LLL-reduced. Then all shortest vectors in $\mathcal{L}$ have form of $\mathbf{w} \mathbf{B}$ where $\mathbf{w} \in \{-1, 0, 1\}^2$.

The following is the pseudocode for the complete algorithm. We need an efficient algorithm to solve $\lambda^2 \equiv -ab \pmod n$, which is covered before. It's also interesting to see if it's possible to avoid solving $\lambda$'s again and again but that's left for future work.

def solve(a, b, n):
    """find all 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 g in divisors(n):
        if n % (g**2) != 0: continue
        n_2 = n // g**2
        for λ in sqrtmod_all(-a * b, n_2):
            B = np.array([
                [n // g, 0],
                [g * (λ * invmod(a, n_2) % n_2), g],
            ])
            B = LLL(B, dot=dot)
            for w in product([-1, 0, 1], [-1, 0, 1]): 
                if w != (0, 0) and gcd(x := w @ B) == g and dot(x, x) == n: 
                    yield x

A special case: $x^2 + d y^2 = n$

Due to incomplete literature review, I missed a large part of works before coming up with the algorithm above. The most related one is Cornacchia's algorithm on solving $x^2 + dy^2 = n$ for $1 \leq d < n$ and $\gcd(d, n) = 1$. [Basilla04] proved the correctness of Cornacchia's algorithm using lattice constructed in the similar way to ours (i.e., the special case of $a = 1$ and $g = 1$).

A failed attempt

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. One 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 solution.

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

Reference

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.

Normalization

Solve a reduced equation

A special case: \(x^2 + d y^2 = n\)

A failed attempt

Reference