How can I solve $m^2+n^2=5077$ $\forall (m,n)\in\mathbb{Z}$, I'm looking for an efficient way to solve this equation 
$$ m^2+n^2=5077 $$
 A: There are algorithms for writing $rs$ as a sum of two squares if you already have $r$ and $s$ as sums of two squares.  By using the identity on complex numbers
$$(a+bi)(c+di) = (ac+bd) + (bc-ad)i$$
and taking absolute values, you get
$$(a^2+b^2) (c^2+d^2) = (ac+bd)^2 + (bc-ad)^2.$$
Unfortunately 5077 is prime, so this doesn't help.  But as Vadim has pointed out there aren't too many combinations to try.  Also useful here is Fermat's theorem on sums of two squares: an odd prime can be written as a sum of two squares if and only if it's of the form $4k+1$.  Since 5077 is of that form you're guaranteed to find a solution.
A: Since $p = 5077$ is a prime of the form $4k+1$, we know it can be factorized as
$$5077 = (m + ni)(m - ni),\quad m, n \in \mathbb{Z}$$
over the Gaussian integers $\mathbb{Z}[i]$.
Furthermore, the two factors $m \pm ni$ are primes over $\mathbb{Z}[i]$.
To find $m + ni$, one first find an integer $x : 1 \le x < p$ such that
$$x^2 + 1 = 0 \pmod p \quad \iff x^2 \equiv -1 \pmod p\tag{*1}$$
We known $(\mathbb{Z}/p\mathbb{Z})^{\times}$, the non-zero elements of the field $\mathbb{Z}/p\mathbb{Z}$, forms a cyclic group with respect to multiplication. If we randomly picking an integer $y : 2 \le y < p$ and by chance picked a generator of $(\mathbb{Z}/p\mathbb{Z})^{\times}$, then $y^k$ will be one possible choice of $x$.
In fact, if we pick $y = 2$, we find
$$x = \mod(2^{\frac{p-1}{4}}, p ) = \mod( 2^{1269}, 5077 ) = 4219$$
satisfies $(*1)$. Now 
$$\begin{align}
& 4219^2 + 1 \equiv 0 \pmod{5077} \\
\implies & 5077\;|\; 4219^2 + 1\\
\implies & \gcd{}_{\mathbb{Z}[i]}(5077,4219+i) \ne 1\\
\implies & \gcd{}_{\mathbb{Z}[i]}(5077,4219+i) = m + n i \text{ or } m - n i
\quad( \text{ up to units over } \mathbb{Z}[i] )
\end{align}$$
The Gaussian integers $\mathbb{Z}[i]$ is known to be an Euclidean domain and we can compute the GCD between
any pair of its elements using some form of Euclidean algorithm. 
Let's take $5077$ and $4219 + i$ as an example.


*

*Initialize 
$$\begin{cases}
a &\leftarrow 4219 + i\\
b &\leftarrow 5077\end{cases}$$

*The closest Gaussian integer to $\frac{a}{b} = \frac{4219 + i}{5077}$ is $1$.
Update the values of $(a,b)$ as
$$\begin{cases}
a &\leftarrow b = 5077\\
b &\leftarrow a - 1 \times b = -858 + i\\
\end{cases}$$

*The closest Gaussian integer to $\frac{a}{b} = \frac{5077}{-858+i}$ is $-6$.
Update the values of $(a,b)$ as
$$\begin{cases}
a & \leftarrow b = -858+i\\
b & \leftarrow a - (-6) \times b = -71 + 6i\\
\end{cases}
$$

*Since $\frac{a}{b} = \frac{-858+i}{-71+6i} = 12 + i$ is an Gaussian integer,
the algorithm terminates and the desired GCD is $-71 + 6i$.


As a consequence, $$5077 = 71^2 + 6^2$$
A: I have found the solution  $6^2+71^2=5077$
I have found it by trial and error
A: This is called Quadratic Diophantine Equation.
There is online solver for such equations - please see here. Also, there is a lot of info about methods of solving this kind of equations on this page.
