Given a large prime, how can I find two perfect squares that add up to that prime?

Question

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.

Answer 1

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

Answer 2

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