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=2551313=8450 , I'd like to generate the four pairs:

1313+9191=8450

2323+8989=8450

3535+8585=8450

4747+7979=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 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=2551313=8450 , I'd like to generate the four pairs:

1313+9191=8450

2323+8989=8450

3535+8585=8450

4747+7979=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)$.