Programming Praxis

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

Equations of the form x² + d y² = 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 x such that x² ≡ −d (mod p). If there are two such values, choose the larger. Then set a ← p and b ← x.

3. While b² > p, set a ← b and b = a mod b. The two assignments are done simultaneously, so the values on the right-hand sides of the two assignments are the old values of the variables. (You may notice that this is Euclid’s algorithm.)

4. After the loop of Step 3 is complete, if d does not evenly divide p − b² or if c = (p − b²) / d is not a perfect square, there is no solution; stop. Otherwise, x = b and y = √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 4k+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 4k+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.

Pages: 1 2