Putting everything together, here is how you count number of ways to decompose $n$ into sum of two squares.
Divide the number by highest power of $4$ in it. If the number is a power of $4$, return $0$.
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.
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$.
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$.
This will count decompositions of the form $0^2+a^2$ too, if you want to avoid counting trivial decompositions, just take floor instead of ceiling.
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!