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