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