Solving $x^2 + y^2 = M$ over integers In the last few days, I have been wondering about the following equation in $x,y \in N$: $$x^2+y^2 = M$$ with $M \in N$. More precisley, I don't understand how can we compute the number of solutions that this equation have and especially how can we find it. 
I've tried to use the formula: $(a^2+b^2)(c^2+d^2)=(ac-bd)^2+(ad+bc)^2$, but it looks me like as another complication, especially because I've to do this calculations for large value of $M$. I've also read the post Integer solutions to $x^2+y^2=N$?, but it hasn't really helped me alot.  
Any ideas or suggestions? Are there also sites which I can consult on this topic? Thanks.
ANOTHER QUESTION (not already answered)
If using the above formula, we can find the solutions to $x^2+y^2=M$, what if M is a very large prime (such as $M=10000000103$): do I have to do only bruteforcing?
EDIT:
As suggested, I have read the page on Wolphram Alpha linked above in the comments, but I still not uderstand in what way we can find the solutions of $x^2+y^2=M$. For example take the case: $x^2+y^2=85$ has $r_2(85)=16$, buit how can we obtain the value of $(x,y)$?
 A: Given $x^2+y^2=85$, and the well-known factorization formula [which you already obviously know] $$(a^2+b^2)(c^2+d^2) = (ac\pm bd)^2+(ad \mp bc)^2.$$
Now just work backwards:


*

*Factor $85=5 \cdot 17.$

*Note that each factor (i.e., $5$ and $17$) have representations as the sum of two squares: $5=4+1=(\pm 2)^2+(\pm 1)^2$ and $17=16+1=(\pm 4)^2+(\pm 1)^2$.

*Determine $x$ and $y$ from $(a,b,c,d)$, going through all possible permutations (including multiple signs).

A: This is in general a very hard problem because any
method which finds solutions to
$\, n = a^2+b^2 = c^2+d^2\,$ allows us to
factor $\,n.\,$
However, factoring large numbers is a very
hard problem according to Wikipedia
Integer factorization.
Using Gaussian integers we have that
 $\,a^2 + b^2 = (a+bi)(a-bi)\,$ and if we
can find that
 $\,n = a^2 + b^2 = c^2 + d^2,\,$ then we
must have that
$\, n = (a+bi)(a-bi)=(c+di)(c-di).\,$
Find the $\,g := \gcd(a+bi,c+di) = e+fi\,$
using the Euclidean algorithm.
Now we have that $\,(e+fi)|(a+bi)\,$ and hence
 $\,(e^2+f^2)|(a^2+b^2).\,$ 
But now $\,e^2+f^2\,$
is a factor of $\,n.$
