Cornacchia's Algorithm

Purpose: Solve the Diophantine equation
$x^2+d{y}^2=m$
for a given $d>0$ and $m$, in $O({\log }^{2}m)$ — essentially one modular square root plus one Euclidean descent.

Algorithm

1. Compute $r_0$, a square root of $-d$ modulo $m$ (via Tonelli-Shanks or Cipolla). If none exists, there is no solution. Take $m/2<r_0<m$.
2. Run the Euclidean algorithm on $(m,r_0)$, stopping at the first remainder $r$ with $r<\sqrt{m}$.
3. Let $s^2=(m-r^2)/d$. If this is an exact integer square, the solution is $(x,y)=(r,s)$; otherwise no solution exists.

Code

// solve x^2 + d*y^2 = m ; returns {-1,-1} if no solution
pair<long long,long long> cornacchia(long long d, long long m) {
    long long r0 = tonelli((m - d % m + m) % m, m);   // sqrt(-d) mod m
    if (r0 < 0) return {-1, -1};
    if (2 * r0 < m) r0 = m - r0;                      // want m/2 < r0 < m
 
    long long a = m, b = r0, lim = (long long)sqrtl((long double)m);
    while (b > lim) { long long t = a % b; a = b; b = t; }
 
    long long rest = m - b * b;
    if (rest % d) return {-1, -1};
    long long s = (long long)sqrtl((long double)(rest / d));
    while (s * s < rest / d) s++;
    while (s * s > rest / d) s--;
    if (s * s * d + b * b != m) return {-1, -1};
    return {b, s};
}

Complexity

• Time: $O({\log }^{2}m)$ — dominated by the modular square root
• Space: $O(1)$

Why It Works

The Euclidean remainders of $(m,r_0)$ correspond to a reduction of the binary quadratic form $x^2+d{y}^2$ of discriminant $-4d$. Gauss’s theory of form reduction guarantees that if $m$ is representable by this form at all, the representation appears at exactly the first remainder dropping below $\sqrt{m}$. Everything else is verification.

Classic special cases

• $d=1$: sum of two squares. $p=x^2+y^2$ has a solution iff $p=2$ or $p\equiv 1(\mathrm{mod}4)$ (Fermat’s theorem on sums of two squares). Cornacchia constructs the representation, which the classical proof does not.
• $d=3$: Loeschian numbers. $p=x^2+3y^2$ iff $p\equiv 1(\mathrm{mod}3)$ — the numbers representable as $a^2+ab+b^2$.
• $d=2$: $p=x^2+2y^2$ iff $p\equiv 1,3(\mathrm{mod}8)$.

Variants / Use Cases

• Sum of two squares decomposition — the most common contest appearance; combine with a prime factorisation and the Gaussian-integer multiplication identity $(a^2+b^2)(c^2+d^2)=(ac-bd{)}^2+(ad+bc{)}^2$ to represent composites
• Sum of four squares — every positive integer is a sum of four squares (Lagrange); a randomised reduction plus Cornacchia gives an efficient construction
• Counting lattice points on circles — $r_2(n)$, the number of representations, follows from the factorisation of $n$ in $\mathbb{Z}[i]$
• Complex multiplication / primality proving (ECPP) — Cornacchia finds the CM discriminant representation needed to build curves of known order
• Quadratic residues — the theory the first step depends on