Solving $x^2 + y^2 = M$ over integers

Question

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)$?

Answer 1

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:

1. Factor $85=5 \cdot 17.$
2. 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$.
3. Determine $x$ and $y$ from $(a,b,c,d)$, going through all possible permutations (including multiple signs).

Answer 2

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.$