## Cornacchia’s Algorithm

### April 3, 2012

We have a fun little problem from number theory today. Take a minute and try to find *x* and *y* so that *x*^{2} + 4 *y*^{2} = 1733. If that’s too easy for you, try *x*^{2} + 58 *y*^{2} = 3031201.

Equations of the form *x*^{2} + *d* *y*^{2} = *p*, with *p* prime, can be solved by an algorithm developed by Giuseppe Cornacchia in 1908 (actually, a fellow named Smith developed the same algorithm in 1885, but his work seems to be forgotten). Here’s the algorithm, which assumes that *p* is prime and 0 < *d* < *p*:

1. If the jacobi symbol (−

d/p) is negative, there is no solution; stop.2. Compute

xsuch thatx^{2}≡ −d(modp). If there are two such values, choose the larger. Then seta←pandb←x.3. While

b^{2}>p, seta←bandb=amodb. The two assignments are done simultaneously, so the values on the right-hand sides of the two assignments are theoldvalues of the variables. (You may notice that this is Euclid’s algorithm.)4. After the loop of Step 3 is complete, if

ddoes not evenly dividep−b^{2}or ifc= (p−b^{2}) /dis not a perfect square, there is no solution; stop. Otherwise,x=bandy= √c.

We can use Cornacchia’s algorithm with *d* = 1 to verify Fermat’s *Christmas* Theorem (because it was announced in a letter to Marin Mersenne on December 25, 1640) that all primes of the form 4*k*+1 can be represented as the sum of two squares; as usual, Fermat didn’t give the proof, which was finally published by Leonhard Euler in a letter to Goldbach in 1749 after seven years of effort.

Your task is to implement Cornacchia’s algorithm and use it to demonstrate that all primes of the form 4*k*+1 and less than 1000 can be written as the sum of two squares. When you are finished, you are welcome to read or run a suggested solution, or to post your own solution or discuss the exercise in the comments below.