Cornacchia's algorithm

The algorithm

First, find any solution to r 0 2 ≡ − d ( mod m ) (perhaps by using an algorithm listed here); if no such r 0 exist, there can be no primitive solution to the original equation. Without loss of generality, we can assume that r₀ ≤ m/2 (if not, then replace r₀ with m - r₀, which will still be a root of -d). Then use the Euclidean algorithm to find r 1 ≡ m ( mod r 0 ) , r 2 ≡ r 0 ( mod r 1 ) and so on; stop when r k < m . If s = m − r k 2 d is an integer, then the solution is x = r k , y = s ; otherwise there is no primitive solution.

To find non-primitive solutions (x, y) where gcd(x, y) = g ≠ 1, note that the existence of such a solution implies that g² divides m (and equivalently, that if m is square-free, then all solutions are primitive). Thus the above algorithm can be used to search for a primitive solution (u, v) to u² + dv² = m/g². If such a solution is found, then (gu, gv) will be a solution to the original equation.

Example

Solve the equation x 2 + 6 y 2 = 103 . A square root of −6 (mod 103) is 32, and 103 ≡ 7 (mod 32); since 7 2 < 103 and 103 − 7 2 6 = 3 , there is a solution x = 7, y = 3.