## Enumerating ways to decompose an integer into the sum of two squares

The well known "Sum of Squares Function" tells you the number of ways you can represent an integer as the sum of two squares. See the link for details, but it is based on counting the factors of the number N into powers of 2, powers of primes = 1 mod 4 and powers of primes = 3 mod 4.

Given such a factorization, it's easy to find the number of ways to decompose N into two squares. But how do you efficiently enumerate the decompositions?

So for example, given N=2*5*5*13*13=8450 , I'd like to generate the four pairs:

13*13+91*91=8450

23*23+89*89=8450

35*35+85*85=8450

47*47+79*79=8450

The obvious algorithm (I used for the above example) is to simply take i=1,2,3,...,$\sqrt{N/2}$ and test if (N-i*i) is a square. But that can be expensive for large N. Is there a way to generate the pairs more efficiently? I already have the factorization of N, which may be useful.

(You can instead iterate between $i=\sqrt{N/2}$ and $\sqrt{N}$ but that's just a constant savings, it's still $O(\sqrt N)$.

-
 The prime factorization of N tells you its prime factorization over the Gaussian integers (en.wikipedia.org/wiki/Gaussian_integer), and then you're just counting all the ways to split N into the product of two Gaussian integers (up to units). – Qiaochu Yuan Jun 26 2010 at 22:07 ("Subfactors" refers to a completely different mathematical concept, so I have removed the tag.) – Qiaochu Yuan Jun 26 2010 at 22:08 If one can obtain two essentially distinct representations: $n=a^2+b^2=c^2+d^2$, then one can factor $n$ nontrivially. Just take the gcd of $a+bi$ and $c+di$ in the Gaussian integers, and take the norm. The moral: it cannot be much harder to factor $n$ first and build up from representations of primes as sums of two squares as suggested by Gerry. – Robin Chapman Jun 27 2010 at 8:01

The factorization of $N$ is useful, since $$(a^2+b^2)(c^2+d^2)=(ac+bd)^2+(ad-bc)^2$$ There are good algorithms for expressing a prime as a sum of two squares or, what amounts to the same thing, finding a square root of minus one modulo $p$. See, e.g., http://www.emis.de/journals/AMEN/2005/030308-1.pdf

Edit: Perhaps I should add a word about solving $x^2\equiv-1\pmod p$. If $a$ is a quadratic non-residue (mod $p$) then we can take $x\equiv a^{(p-1)/4}\pmod p$. In practice, you can find a quadratic non-residue pretty quickly by just trying small numbers in turn, or trying (pseudo-)random numbers.

-
 So what would the algorithm itself be? It sounds like I should enumerate all possible $xy=N$ factorings (both prime and composite). Then for each, decompose $x$ into each possible $a^2+b^2$ and each y into each possible $c^2+d^2$, and use the above formula to find one answer to the top level N decomposition. Finally after iterating over all such factors, and over the two inner loops of all decompositions of those factors, I should take all the answers and sort them and eliminate duplicates. Is this the right algorithm or is it doing unnecessary work? – MathMonkey Jun 27 2010 at 2:50 And finally, is it guaranteed that the above algorithm will actually find ALL of the top level N decompositions? The formula just tells us that given one factoring we get one sum of two squares decompositions, but does that mean that all factorings will give us all decompositions? – MathMonkey Jun 27 2010 at 2:53 I think it will be easier for you to learn by doing than by me explaining (since I'm not so hot at explaining). Take an example, say, $N=5\times13\times17$, and use the expressions $5=2^2+1^2$, $13=3^2+2^2$, $17=4^2+1^2$, and see what you have to do to get all 4 distinct representations. Qiaochu Yuan's remark may help guide you; in effect, we're finding all 4 values of $a+bi=(2+i)(3\pm2i)(4\pm i)$ where $i$ is the square root of minus one, and our representations are $N=a^2+b^2$. Yes, this is guaranteed to get everything, and without duplicates if you set it up right. Try it and see. – Gerry Myerson Jun 27 2010 at 6:12 Yes, all expressoins as a sum of squares occur in this way. This is an immediate consequence of unique factorization in $\mathbb{Z}[i]$. – David Speyer Jun 27 2010 at 14:32 Hi, Gerry. I remembered the Stan Wagon piece with this material, jstor.org/pss/2323912 – Will Jagy Jun 28 2010 at 1:21

I wouldn't mind an elaboration on Gerry's hints. For example, N=5*5*13*13*17*17 is going to have 13 representations. What variants of x +/- yi are we multiplying together to come up with those? In Gerry's example, how come we don't look at (2-i) as a term and end up with 8 representations?

-
 The reason the conjugate wasn't used on one term was because there is a symmetry on the whole set by conjugating everything. If it's confusing, write out the whole set without worrying about the conjugation symmetry first, then note that you have solutions which only differ by conjugation or equivalently, switching $a^2$ and $b^2$. There are many choices for how to eliminate the duplicates. For $5*5*x$, you have two choices of sign in $(2\pm i)(2 \pm i)$ and what matters is the count of plus signs. Afterwards, worry about the duplicates from the conjugation symmetry. – Douglas Zare Jul 20 at 22:25

This is the simplest case of the Hardy-Muskat-Williams algorithm. Anyway, here is a link to a 1995 paper by Kenneth S. Williams, http://www.mathstat.carleton.ca/~williams/papers/pdf/202.pdf and to the original HMW paper http://www.ams.org/journals/mcom/1990-55-191/S0025-5718-1990-1023762-3/S0025-5718-1990-1023762-3.pdf .

As I'm not sure you are aware of these details, let me point out that if $$4^k \;| \; \; x^2 + y^2$$ then $2^k \; | \; x$ and $2^k \; | \; y.$ That is, you might as well divide your target by powers of 4 before doing anything difficult. Then after you are finished multiply $x,y$ by the appropriate power of $2.$

This is very similar. If there is a prime $$q \equiv 3 \pmod 4$$ and $q | n,$ then keep dividing the target by powers of $q^2$ until it is no longer divisible by $q^2.$ If the remaining number is divisible by $q$ there is actually no representation at all. But if $$q^{2k} \;\parallel \; \; x^2 + y^2$$ then $q^k \; | x$ and $q^k \; | y.$ The notation $q^{2k} \;\parallel \; \; x^2 + y^2$ means $q^{2k} \; | \; \; x^2 + y^2$ but it is not true that $q^{2k +1} \; | \; \; x^2 + y^2$

Well, that is enough caution. What you really need to know is expressing primes $$p \equiv 1 \pmod 4$$ and indeed $p^m,$ which is not much more difficult. Once you can do that, combine my notes with all possible ways of applying Gerry's multiplication formula (by changing $\pm$ signs and order),

-