# 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). Jun 26, 2010 at 22:07
• ("Subfactors" refers to a completely different mathematical concept, so I have removed the tag.) 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. Jun 27, 2010 at 8:01
• Note also $65^2+65^2$ Feb 17, 2018 at 8:07

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? 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? 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. 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]$. Jun 27, 2010 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)$. Dec 7, 2014 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),

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.

Another point of view, which might be easier from an algorithmic point of view (in my opinion), is to look at the decomposition $$n = a^2+b^2$$ using complex numbers $$n = (a+bi)(a-bi)$$. Suppose you know how to find the solutions $$a,b$$ when $$n$$ is prime. Then note that the identity $$(a^2+b^2)(c^2+d^2)=(ac+bd)^2+(ad-bc)^2$$ can be seen as the equality of the product of modules of two complex numbers.

In order to find all possible decompositions of a number $$n$$ into a sum of squares just write the factorization of $$n$$ in $$\Bbb{Z}[i]$$: $$n = \prod (a_j+ib_j)$$ and note that in this product every factor comes with its conjugate. First, ignore the powers of $$2$$ and the primes of the form $$4k+3$$. Now, in order to find all factorizations, just split all factors into two columns with conjugate pairs being on different columns. Doing this, when taking the product on each column we'll get a pair of conjugate numbers whose product equals to $$n$$, so we find a solution of the representation $$n=a^2+b^2$$. The number of ways to split the factors in the two columns with conjugate pairs on different sides will be equal the number of divisors of $$n$$ in $$\Bbb{Z}[i]$$ (divided by 4, since you can multiply factors by $$1,-1,i,-i$$) so all factors will be generated.

If $$n$$ contains primes of the forms $$4k+3$$ or powers of $$2$$ look at Will Jagy's answer. Note that if you reduce all powers of $$4$$ and you still have a $$2$$ left in $$n$$ you can write $$2 = (1+i)(1-i)$$ and split this on different sides of the two columns.

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 May 12, 2015 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? May 12, 2015 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. May 14, 2015 at 2:00
• @KalEl please post the python code. I would like to try this .
– john
Sep 2, 2015 at 3:32
• @GHfromMO Indeed it is, I just wanted to highlight the fact in a comment, because there is a solution missing from the original post - the correct formula counts it. The question was linked from Math Stackexchange and I didn't want any misconceptions from people following the link. Feb 17, 2018 at 18:23