I need to emulate this sequence for a program: <http://oeis.org/A025302>


Stuff that I've taken into account:


* After finding the prime divisors of a number. I take any divisor as **p** and apply the following rule:


**p is an odd prime divisor of n, then either p appears to an even power in n, or p≡1(mod4)**


* Also I only consider sums of 2 distinct nonzero squares


The sequence in the described link only goes as far as 229. Until here my simulator using these 2 rules work, but I'm not so sure after that since I get these numbers ( ***which I'm not sure that they belong in my sequence*** ) when I loop until *n <= 1000*..



```
[0] 325 int
[1] 425 int
[2] 650 int
[3] 725 int
[4] 845 int
[5] 850 int
[6] 925 int

```

So my question is:


Do these 7 digits belong originally to the <http://oeis.org/A025302> sequence, or have I gotten wrong my algorithm, if so, can you point the missing rule ?


Any question, feel free to comment



Numbers that are the sum of 2 distinct nonzero squares in exactly 1 way
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](https://en.wikipedia.org/wiki/Sum_of_two_squares_theorem)), 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](https://www.math.uchicago.edu/~may/VIGRE/VIGRE2008/REUPapers/Bhaskar.pdf) and as for some code the function [PowersRepresentations[n, 2, 2]](https://reference.wolfram.com/language/ref/PowersRepresentations.html) could also be of some use to your purposes (however note that this function also allows zeros). 


