Given a large prime, how can I find two perfect squares that add up to that prime? I have a large prime number, $106202791239577$.  How can I find two perfect squares that when added together,  equal this number?  Note that $106202791239577$  is of the form $4n{+}1$ (it is congruent to $1 \bmod 4$ ), so it should be possible to find two squares that  add up to this prime.
 A: I'll illustrate one of the techniques at the question lulu linked although due to limitations on the tools I have available I'll use smaller numbers:
Say we want to find the squares that add to the prime $9874577$. Then we can use an exponentiation by squaring (towards $p^{9874576}$) on some small primes to observe that $2^{2468644} \equiv 9874576 \equiv -1 \bmod 9874577$ and that the step before, $2^{1234322} \equiv 1698670$ which is thus a square root of $-1 \bmod 9874577$.
Then we can use the Euclidean GCD algorithm for Gaussian integers on $9874577$ and $1698670+i$ to find that $\gcd(9874577,1698670+i) = 2924+1151i$ and so get $2924^2+1151^2 = 9874577$.
A: $$106202791239577=9996044^2+2506371^2$$
The target number is, as I expect you know, a Pythagorean prime https://en.wikipedia.org/wiki/Pythagorean_prime
As such, it has only one way to be expressed as the sum of two squares,  
This document http://eulerarchive.maa.org/docs/translations/E228en.pdf in section 44 starting on page 17, as far as I can understand, gives a manual method to test a number of the form $4n+1$ to be prime, by finding then number of ways it can be expressed as the sum of two squares. I’ve not tried this method.
Consider the equation $n=x^2+y^2$, where $n$ does not have to be prime.
We could take the brute force methodology, and look at $y=(n-x^2)^{0.5}$ for all values of $x$ from $1$ to $n^{0.5}$, but this is not only wasteful, but gives two solutions for each “real” solution (in this case $9996044^2+2506371^2$ and $2506371^2+9996044^2$).
If we define $x>=y$, it’s sufficient to constrain $x$ to the range 
$$((0.5n)^{0.5}),(n^{0.5})$$
In this case it’s $x=(7287069,10305473)$.
Perhaps surprisingly, I used Excel 2016 to find the solution, copying thousands of rows, each calculating the new $x$ from the line above, then the value of $y$. Next I searched for $y$ values containing $.000000$, and copied those values.
Clearly, there’s too much data for one sheet, so I copied the bottom $x$ value to the top.
This sounds a lot of trouble, but took less than ten minutes, easily found the solution and two near misses. Please, do let me know if you need more details.
There is a product which I’ve used, years ago, to solve this type of problem: Excel Solver Add In.
