The well known ["Sum of Squares Function"](http://mathworld.wolfram.com/SumofSquaresFunction.html) 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)$.



Enumerating ways to decompose an integer into the sum of two squares
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. 


