Cornacchia’s algorithm turns a modular square root into an integer representation by a quadratic form.
For equations such as
x²+d y²=m,
the algorithm begins with a root of
r²≡−d (mod m)
and then runs an ordinary Euclidean remainder sequence. If the right remainder appears, the missing coordinate is recovered by one exact square test.
The method is short, constructive and surprisingly powerful. It connects modular square roots, continued fractions, Euclidean algorithms, norm equations, Gaussian integers and binary quadratic forms.
This is Guide 31 in the Bukit Timah Tutor Computational Number Theory series. It follows Guide 30: Gaussian Integers and links back to Guide 17: Binary Quadratic Forms.
1. The target equation
We seek integer solutions to
x²+d y²=m
with positive integers d and m.
Rearranging modulo m gives
x²≡−d y² (mod m).
If gcd(y,m)=1, then
(x y^−1)²≡−d (mod m).
Thus a solution implies that −d has a square root modulo m under suitable coprimality conditions.
Cornacchia reverses this implication constructively.
2. The algorithmic skeleton
Suppose r satisfies
r²≡−d (mod m)
with 0<r<m.
Run the Euclidean algorithm on m and r:
r_0=m, r_1=r, r_(k+1)=r_(k−1) mod r_k.
Stop at the first remainder x with
x²<m.
Then test whether
y²=(m−x²)/d
is a nonnegative integer square.
If yes, (x,y) is a solution.
3. Worked example: 29=x²+y²
Take
m=29, d=1.
We need a root of
r²≡−1 (mod29).
One choice is r=12 because
12²=144≡28≡−1 (mod29).
4. Euclidean descent for 29
Run Euclid:
29 =2·12+5 12 =2·5+2 5 =2·2+1.
Since
5²=25<29,
the first remainder below √29 is x=5.
Now compute
y²=29−25=4.
So y=2 and
29=5²+2².
5. Why the stopping threshold is √m
If x²+d y²=m, then x²<m unless y=0.
The Euclidean remainder sequence is decreasing. Cornacchia searches for the candidate remainder that has crossed below √m but still retains the arithmetic information carried by the modular root.
Stopping much later can miss the relevant candidate. Stopping earlier leaves x too large to fit the target equation.
6. A second example: 41=x²+y²
Choose r=9 because
9²=81≡40≡−1 (mod41).
Run Euclid:
41 =4·9+5 9 =1·5+4 5 =1·4+1.
The first remainder below √41 is x=5.
Then
y²=41−25=16,
so
41=5²+4².
7. General d: represent 17 by x²+2y²
Take
m=17, d=2.
We need
r²≡−2≡15 (mod17).
Since
7²=49≡15 (mod17),
choose r=7.
Euclid gives
17=2·7+3 7 =2·3+1.
The first remainder below √17 is x=3.
Now
y²=(17−9)/2=4.
Thus
17=3²+2·2².
8. A modular root is necessary but not always sufficient
Finding r with r²≡−d mod m does not guarantee the final quotient becomes a square for every composite m and every formulation.
The Euclidean descent may produce x for which
(m−x²)/d
is not an integer square.
Then that branch fails. If another modular square root exists, try it according to the algorithm’s exact hypotheses and symmetry.
9. Two modular roots
For odd prime m and nonzero quadratic residue −d, roots occur in pairs r and m−r.
These two roots correspond to conjugate or sign-related arithmetic data. In the simplest prime sum-of-two-squares case, either branch leads to the same representation up to signs and swapping.
For more general norm forms and composite moduli, branch handling should be stated explicitly.
10. Cornacchia and Gaussian gcd are two views of one reconstruction
For d=1 and prime p≡1 mod4, Guide30 used
gcd_Z[i](p,r+i)
to extract a Gaussian factor of norm p.
Cornacchia uses the ordinary integer Euclidean remainders of p and r.
Both algorithms are exploiting the same modular square root of −1 and the same hidden norm factorisation. One works explicitly in Z[i]; the other compresses the reconstruction into integer arithmetic.
11. Norm interpretation
In the quadratic field Q(√−d),
N(x+y√−d)=x²+d y².
So solving
x²+d y²=m
is a norm equation.
However, the full ring of integers may not equal Z[√−d]. The simple form x²+d y² corresponds to a specific order/basis. General number-field norm equations require more careful integral-basis arithmetic.
12. Binary quadratic forms
The equation
x²+d y²=m
is representation of m by the binary quadratic form
[1,0,d].
More general primitive forms
ax²+bxy+cy²
have discriminant b²−4ac and representation theory governed by form classes.
Cornacchia-type algorithms can be adapted to related forms, but the clean textbook version is the diagonal norm form above.
13. Square-root computation is a separate subproblem
For prime m=p, obtaining r requires solving
r²≡−d (mod p).
Use Euler’s criterion to test solvability and Tonelli–Shanks to compute a root for odd p.
Thus Cornacchia is downstream of modular square-root algorithms from Guide5.
14. When d=1 and p≡1 mod4
For prime p≡1 mod4, −1 is always a quadratic residue, so the modular-root stage succeeds.
Cornacchia then constructs the unique positive representation
p=x²+y²
up to swapping x,y.
This is a constructive form of Fermat’s two-square theorem.
15. When p≡3 mod4
For prime p≡3 mod4, −1 is not a quadratic residue modulo p.
The algorithm cannot even start for d=1 because no r exists.
This modular obstruction matches the theorem that such primes are not sums of two squares.
16. General solvability filters
Before running a reconstruction, test local conditions.
If a prime q divides m, the congruence
x²+d y²≡0 (mod q)
may impose residue restrictions.
Legendre/Jacobi symbols, valuation arguments and congruences modulo small powers can reject impossible equations cheaply.
Local solvability is not always global solvability, but local obstructions are efficient early exits.
17. Exact square testing
After Euclidean descent, compute
t=(m−x²)/d.
Do not convert t to floating point and compare sqrt(t) approximately.
Use integer square root:
y=floor(sqrt(t)) accept iff y²=t.
This avoids rounding errors for large integers.
18. Pseudocode
cornacchia(m,d,r):
require 0 < r < m
require r*r ≡ -d (mod m)
a = m
b = r
while b*b > m:
a, b = b, a mod b
x = b
t = m - x*x
if t < 0 or t mod d != 0:
fail
y2 = t // d
y = isqrt(y2)
if y*y != y2:
fail
return (x,y)The exact inequality at the stopping boundary should follow the chosen theorem statement; implementations should handle equality and y=0 cases explicitly.
19. Complexity
The Euclidean stage is extremely cheap: O(log m) remainder steps, with bit complexity governed by fast integer division.
For prime m, modular square-root computation is usually the more substantial stage.
Thus Cornacchia adds little overhead once a root of −d is available.
20. Prime-power and composite moduli
For composite m, square roots modulo m may require factoring m or using known CRT structure. Multiple roots can occur.
The clean prime case avoids this complexity. In applications where m’s factorisation is known, compute roots modulo prime powers, combine with CRT and test the resulting branches.
The number of roots can grow multiplicatively across components.
21. Cornacchia and Cornacchia–Smith variants
Literature contains several related algorithms and generalisations, sometimes called Cornacchia–Smith methods, for equations of the form
x²+d y²=m
under different gcd and parity conditions.
A production implementation should state precisely which theorem and hypotheses it implements rather than using the name Cornacchia as a catch-all label.
22. CM and elliptic-curve applications
Norm-form equations appear in complex multiplication computations.
ECPP and CM curve construction may require solving equations such as
4N=u²+|D|v²
for discriminants D.
Cornacchia-style algorithms can solve these representation problems efficiently once the appropriate modular square root is known.
This creates a direct bridge from a classical Euclidean algorithm to modern elliptic-curve primality proving.
23. Representation certificates
A returned pair (x,y) is self-verifying:
x²+d y² ?= m.
For a full algorithm receipt, also retain the modular root r and Euclidean remainder sequence.
Then another implementation can verify both the starting congruence and the reconstruction path.
24. Common mistakes
1. Running Euclid without first checking r²≡−d mod m. 2. Using a floating square test. 3. Assuming any modular root guarantees a global representation. 4. Forgetting divisibility by d in the final t calculation.
5. Using the wrong stopping condition. 6. Treating Q(√−d) coordinates as the full ring of integers in every d. 7. Ignoring branch multiplicity for composite moduli. 8. Accepting a representation without direct substitution.
25. Practice set
1. State the target norm equation. 2. What congruence must r satisfy? 3. Why is r=12 valid for m=29,d=1? 4. Write the Euclidean remainder sequence.
5. Recover 29 as x²+y². 6. Recover 41 using r=9. 7. Verify 7²≡−2 mod17. 8. Recover 17 as x²+2y².
9. Why does p≡3 mod4 block the d=1 prime case? 10. What algorithm can find r modulo an odd prime? 11. Why use integer square root? 12. What is the norm interpretation?
13. Which binary quadratic form corresponds to x²+d y²? 14. What should a failure mean? 15. What extra complexity appears for composite m? 16. How can this algorithm appear in ECPP/CM work?
26. Answers
1. x²+d y²=m. 2. r²≡−d mod m. 3. 12²=144≡−1 mod29. 4. 29,12,5,2,1; the first remainder below √29 is5.
5. 29=5²+2². 6. 41=5²+4². 7. 49≡15≡−2 mod17. 8. 17=3²+2·2².
9. −1 is not a quadratic residue modulo such a prime. 10. Tonelli–Shanks. 11. To make the final square test exact for arbitrarily large integers. 12. x²+d y² is the norm of x+y√−d in the simple quadratic basis.
13. [1,0,d]. 14. That the tested root/branch does not yield the requested representation under the algorithm’s hypotheses. 15. Multiple square roots and CRT branch handling. 16. CM order equations reduce to quadratic norm representations.
Sources and further study
Henri Cohen, A Course in Computational Algebraic Number Theory, develops Cornacchia-type norm-equation methods in computational algebraic number theory. The modular-square-root prerequisites are treated in Victor Shoup’s A Computational Introduction to Number Theory and Algebra.
Continue through Batch 08
Return to Guide 30. Continue to Guide 32: Higher Power Residues, Character Sums, Jacobi Sums and Cyclotomic Arithmetic.