Copying from http://mathworld.wolfram.com/SumofSquaresFunction.html: 

>To find in how many ways a positive integer $n>1$ can be expressed as a sum of $k=2$ squares ignoring order and signs, factor it as
$$n=2^{a_0}p_1^{2a_1}...p_r^{2a_r}q_1^{b_1}...q_s^{b_s}$$ 	
where the $p_i$ are primes of the form $4k+3$ and the $q_i$ are primes of the form $4k+1$. If $n$ does not have such a representation with integer $a_i$ because one or more of the powers of $p_i$ is odd, then there are no representations. Otherwise, define
 $$B=(b_1+1)(b_2+1)...(b_r+1)$$ 	
The number of representations of $n$ as the sum of two squares ignoring order and signs is then given by:  

>[![enter image description here][1]][1]

You might also be interested in Theorem 4.4, p. 10 from [these notes][2] and as for some code the function [PowersRepresentations[n, 2, 2]][3] might be also useful for your purposes. 


  [1]: https://i.sstatic.net/XyU8x.png
  [2]: https://www.math.uchicago.edu/~may/VIGRE/VIGRE2008/REUPapers/Bhaskar.pdf
  [3]: https://reference.wolfram.com/language/ref/PowersRepresentations.html