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

Question

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)$.

Answer 1

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.

Answer 2

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

Answer 3

(This elaborates on Gerry's answer.)

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.

Answer 4

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.

Answer 5

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!