# 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 '10 at 22:07
("Subfactors" refers to a completely different mathematical concept, so I have removed the tag.) –  Qiaochu Yuan Jun 26 '10 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 '10 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 '10 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 '10 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 '10 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 '10 at 14:32
@pts, why not experiment a bit, and see for yourself? $5=2^1+1^1$, $13^2=12^2+5^2=13^2+0^2$, mix 'n' match, see what happens? Or, $(2\pm i)(3\pm2i)(3+2i)$. –  Gerry Myerson Dec 7 '14 at 5:43

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),

-

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 '12 at 22:25
Hi... it is 13 because you have 3 primes with powers of 2 each, and $\lfloor(2+1)(2+1)(2+1)/2\rfloor=13$. Likewise for $N=5*5*5*13*13$, the count is $\lfloor(3+1)(2+1)/2\rfloor=6$ and so on. Hope that helps. See my detailed answer below. –  KalEl May 12 at 15:24

This article describes how to solve the $p=x^2+y^2$ equation quickly if $p\equiv 1$ mod 4 and $p$ is a prime.

John Brillhart: Note on representing a prime as a sum of two squares

It also explains how $x^2\equiv-1$ mod $p$ can be solved.

-

Putting everything together, here is how you count number of ways to decompose $n$ into sum of two squares.

For this we are gonna count including trivial representation $0^2+a^2$ for square numbers. (To get rid of it, you may subtract $1$ if the number is a perfect square.)

1. Divide the number by highest power of $4$ in it. If the number is a power of $4$, return $1$.

2. Decompose what remains into prime factors.

a. If there is a prime factor of the form $4n+3$ with odd power, return $0$.

b. Discard all prime factors form $4n+3$ with even power.

3. Now you have all prime factors of the form $4n+1$, and possibly a $2$ hanging around in the decomposition. Let's say you have $2^{n_0}\prod_{k=1}^m p_k^{n_k}$ with $p_k\equiv1\mod 4$, and $n_0$ being either $0$ or $1$.

4. Then number of ways $n$ can be decomposed in sum of square of pairs is $\left\lceil\frac{\prod_{k=1}^m (n_k+1)}{2}\right\rceil$.

If you want to actually enumerate instead of count, you will need two things, 1) To be keep track of powers you discarded and 2) To be able to extract root of $-1$ modulo $p$, and use it to factorize $4n+1$ into a gaussian integer and its conjugate. It's just a bit more of work but isn't difficult - I wrote the code based on the discussion here, and some papers referred here, it works pretty well!

-
This is classical and can be found in many textbooks. See e.g. mathworld.wolfram.com/SumofSquaresFunction.html –  GH from MO May 12 at 15:54
As you are counting decompositions of the form $0^2+a^2$, shouldn't a power of 4 return 1, not 0? –  Gerry Myerson May 12 at 23:17
Gerry - you are right! Corrected it. Anyone looking for the enumeration algorithm, it will be clear at once if you see how the counting formula works. If you need it I can share my Python code. –  KalEl May 14 at 2:00