Expressing a prime as the sum of two squares

I saw where Elon Musk posted Grok's answer to the prompt What are the most beautiful theorems." I looked at the list, and there were no surprises, as you'd expect from a program that works by predicting the most likely sequence of words based on analyzing web pages.

There's only one theorem on the list that hasn't appeared on this blog, as far as I can recall, and that's Fermat's theorem that an odd primep can be written as the sum of two squares if and only ifp = 1 mod 4. The only if" direction is easy [1] but the if" direction takes more effort to prove.

Ifp is a prime andp = 1 mod 4, Fermat's theorem guarantees the existence ofx andy such that

Gauss' formula

Stan Wagon [2] gave an algorithm for finding a pair (x,y) to satisfy the equation above [2]. He also presents a beautiful formula due to Gauss" which does not seem to be of any value for computation." Gauss' formula says that ifp = 4k + 1, then a solution is

For x and y we choose the residues mod p with |x| and |y| less thanp/2.

Why would Wagon say Gauss' formula is computationally useless? The number of multiplications required is apparently on the order ofp and the size of the numbers involved grows likep!.

You can get around the problem of intermediate numbers getting too large by carrying out all calculations modp, but I don't see a way of implementing Gauss' formula with less thanO(p) modular multiplications [3].

Wagon's algorithm

If we want to express a large primep as a sum of two squares, an algorithm requiring O(p) multiplications is impractical. Wagon's algorithm is much more efficient.

You can find the details of Wagon's algorithm in [3], but the two key components are finding a quadratic non-residue mod p (a numberc such thatc x^2 mod p for any x) and the Euclidean algorithm. Since half the numbers between 1 andp - 1 are quadratic non-residues, you're very likely to find a non-residue after a few attempts.

[1] The square of an integer is either equal to 0 or 1 mod 4, so the sum of two squares cannot equal 3 mod 4.

[2] Stan Wagon. The Euclidean Algorithm Strikes Again. The American Mathematical Monthly, Vol. 97, No. 2 (Feb., 1990), pp. 125-129.

[3] Wilson's theorem gives a fast way to compute (n - 1)! modn. Maybe there's some analogous identity that could speed up the calculation of the necessary factorials mod p, but I don't know what it would be.

The post Expressing a prime as the sum of two squares first appeared on John D. Cook.