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:  

>$$
r'_2(n)=
\left\{
  \begin{array}{lll}
   0 & \textrm{if any }  a_i \textrm{ is a half integer},  \\ 
\frac{1}{2}B & \textrm{if all }  a_i \textrm{ are integers and } B \textrm{ is even}, \\ 
\frac{1}{2}\big(B-(-1)^{a_0}\big) & \textrm{if all }  a_i \textrm{ are integers and } B \textrm{ is odd}
  \end{array}
\right.
$$

which implies that the missing rule in your algorithm (which btw is just the [Sum of two squares theorem][2]), is the following:  
In order to have $r'_2(n)=1$, there can only be three cases: 

 - Either: $n$ has a single $4k+1$ prime factor with multiplicity $b=1$, hence  $B=b+1=1+1=2$ and $r'_2(n)=\frac{1}{2}B=1$  
($n=5$ provides an example for that case),
 - or: $n$ has a single $4k+1$ prime factor with multiplicity $b=2$ and the multiplicity (in the prime factorization of $n$) of $2$ is $a_0$:even, hence  $B=b+1=2+1=3$ and $r'_2(n)=\frac{1}{2}(B-(-1)^{a_0})=\frac{1}{2}(3-1)=1$  
 ($n=100=2^2\cdot 5^2$ provides an example for that case), 
 - or: there are no $4k+1$ prime factors in the prime factorization of $n$ and the multiplicity (in the prime factorization of $n$) of $2$ is $a_0$:odd, hence in that case, all $b_j$ are zero which gives $B=1$ thus $r'_2(n)=\frac{1}{2}(B-(-1)^{a_0})=\frac{1}{2}(1+1)=1$.  
However, this last case should be excluded from your considerations since you are asking for numbers $n$ which can be written in a unique way as a sum of two distinct squares while this last case leads to numbers which are twice a square ($n=18=2\cdot 3^2=3^2+3^2$ for example). 

You might also be interested in Theorem 4.4, p. 10 from [these notes][3] and as for some code the function [PowersRepresentations[n, 2, 2]][4] could also be of some use to your purposes (however note that this function also allows zeros). 


  [2]: https://en.wikipedia.org/wiki/Sum_of_two_squares_theorem
  [3]: https://www.math.uchicago.edu/~may/VIGRE/VIGRE2008/REUPapers/Bhaskar.pdf
  [4]: https://reference.wolfram.com/language/ref/PowersRepresentations.html