Efficiently finding two squares which sum to a prime

In practice this comes to the same algorithm as Moron's, but I'll give my two penn'orth anyway. There are two stages: (i) find $z$ with $z^2\equiv -1$ (mod $p$), and (ii) use $z$ to find $x$ and $y$.

Stage (i). If $a$ is a quadratic nonresidue modulo $p$ then $a^{(p-1)/2} \equiv-1$ (mod $p$) so we can take $z\equiv a^{(p-1)/4}$ (mod $p$) (and that is easily computed by modular exponentiation). How does one find quadratic nonresidues? Well just over half of all numbers $a$ such that $1 < a < p-1$ are quadratic nonresidues, so taking $a$ at random will need at most two goes on average.

Stage (ii). Compute the greatest common divisor of $p$ and $z+i$ using the Euclidean algorithm for the Gaussian integers. The answer will be $x+yi$ where $x^2+y^2=p$.

If $p\equiv 1$ (mod $2^k$) where $k\ge3$ one can speed up stage (i) a bit on average. Compute $w=a^{(p-1)/2^k}$ modulo $p$. If $w\equiv\pm1$ we lose but otherwise keep squaring $w$. Then we eventually get to $-1$ and before that we have the required $z$. This sometimes wins even if $a$ is a quadratic residue.

In effect the Gaussian gcd above and Moron's Hermite-Serret algorithm amount to the same thing.

Below is some rough-and-ready Python code to compute $a$ and $b$ such that $p=a^2+b^2$. The function that does this is 2sq. It will behave erratically if fed p not a prime congruent to 1 modulo 4.

def mods(a, n):
    if n <= 0:
        return "negative modulus"
    a = a % n
    if (2 * a > n):
        a -= n
    return a

def powmods(a, r, n):
    out = 1
    while r > 0:
        if (r % 2) == 1:
            r -= 1
            out = mods(out * a, n)
        r //= 2
        a = mods(a * a, n)
    return out

def quos(a, n):
    if n <= 0:
        return "negative modulus"
    return (a - mods(a, n))//n

def grem(w, z):
    # remainder in Gaussian integers when dividing w by z
    (w0, w1) = w
    (z0, z1) = z
    n = z0 * z0 + z1 * z1
    if n == 0:
        return "division by zero"
    u0 = quos(w0 * z0 + w1 * z1, n)
    u1 = quos(w1 * z0 - w0 * z1, n)
    return(w0 - z0 * u0 + z1 * u1,
           w1 - z0 * u1 - z1 * u0)

def ggcd(w, z):
    while z != (0,0):
        w, z = z, grem(w, z)
    return w

def root4(p):
    # 4th root of 1 modulo p
    if p <= 1:
        return "too small"
    if (p % 4) != 1:
        return "not congruent to 1"
    k = p//4
    j = 2
    while True:
        a = powmods(j, k, p)
        b = mods(a * a, p)
        if b == -1:
            return a
        if b != 1:
            return "not prime"
        j += 1

def sq2(p):
    a = root4(p)
    return ggcd((p,0),(a,1))

Just for amusement, here is an exact formula for finding $a$ so that $a^2 + b^2 = p$:

$a = \dfrac{1}{2}\sum_{x = 0}^{p-1} \bigl( \dfrac{x^3 - x}{p} \bigr),$

where $\bigl( \dfrac{a}{p} \bigr)$ denotes the Legendre symbol of $a$ mod $p$.

Added: At Bill Dubuque's suggestion, I am noting here that this exact formula is not an efficient means to compute the value of $a$ (hence my prefatory "Just for amusement ...").

We can quickly compute a representation of a prime $ \,p\equiv 1\pmod{\!4}\,$ as a sum of two squares by using the Euclidean GCD algorithm in $ \mathbb Z[i]$ and an algorithm for computing square roots $\!\pmod{\!p}.$

Theorem $\ $ Let $ \ \, c = \sqrt{-1}\pmod{\!p}\, $ and $\,\gcd(p,\,i-c) = a+b\,i\ $.$\ $ Then $ \ p = a^2 + b^2\,$.

Proof $\ $ We shall show: $\ (1)\, $ Representing $ \, p\,$ as a sum of squares is equivalent to finding a proper splitting $ \, p = \alpha\beta\ $ in $ \, \mathbb Z[i].\,$ $(2)\ $ Since $\,\Bbb Z[i]\,$ has a Euclidean algorithm, a proper splitting of $ \, p\,$ can be computed by GCD with a suitable splitting of some multiple of $ \, p.\,$ $ (3)\ $ A suitable splitting of some multiple of $ \, p\,$ arises by factoring $\, x^2+1\pmod{\!p},\,$ i.e. by computing $ \ \ \sqrt{-1}\pmod{\!p}.$

$(1)\ $ If $ \,\alpha\mid p\ $ properly then, conjugating, $ \,\alpha'\mid p\ $ properly. Multiplying $\,\Rightarrow\,\alpha\alpha'\mid p^2$ properly in $\,\Bbb Z,\,$ where $\,p\,$ is the only proper factor$>0$ of $\,p^2,\,$ so $ \, p = \alpha\alpha' = (a\!+\!b\,i)(a\!-\!b\,i) = a^2\! + b^2.$

$(2)\, \ p\gamma = \alpha\beta,\,$ $p\nmid\alpha, \beta\,$ $\Rightarrow\, \gcd(\alpha,p)\mid p\,$ properly. Else $\,\gcd(\alpha,p) = p\,$ or $1.\,$ If the gcd is $\,p\,$ then $ \, p\mid\alpha\,$ contra hypothesis. Hence $ \ \gcd(\alpha,p)=1\, $ so by Euclid's Lemma $\, \alpha\mid p\,\gamma\, \Rightarrow\, \alpha\mid\gamma\,$ thus $\, \gamma/\alpha = \beta/p\, \Rightarrow\, p\mid \beta,\,$ again contra hypothesis. $ $ Note: $ $ generally, in rings, GCDs are unique only up to unit multiples. Here the units are $ \ i^n = \pm 1,\,\pm i\,$.

$(3)\,\ x^2 + 1 = (x\!-\!c)(x\!+\!c) + p\, f(x)\ $ in $\,\Bbb Z[x]\, \overset{\large x\,=\,i}\Longrightarrow -p\, f(i) = (i\!-\!c)(i\!+\!c)\, $ in $\,\Bbb Z[i].\,$ By $(2)$ this splitting is suitable to split $ \,p\,$ since $ \,p\nmid i\pm c\,$ in $\,\Bbb Z[i].$

Remark $ $ There are many variations on the Euclidean algorithm in the Gaussian integers $\,\Bbb Z[i],\,$ e.g. employing continued fractions, binary quadratic forms, etc. Also, there are also at least a few algorithms for computing sqrts $\!\pmod{\!p},\,$ e.g. by factoring polynomials, by elliptic curves (Schoof), and algorithms of Tonelli and Shanks. For much further information see Henri Cohen's book "A course in computational algebraic number theory". Below is an excerpt on Cornacchia's algorithm.

You can try using the Hermite-Serret algorithm, which is briefly as follows:

Find a $z$ such that $z^2 = -1 \mod p$. Apply the Euclidean algorithm to $p$ and $z$ stopping when the first pair of remainders $(x,y)$ go below $\sqrt{p}$. Hermite proved that $x^2 + y^2 = p$.

To find $z^2 = -1 \mod p$, you can probably guess a $z$ and start with $z^{p-1} = 1 \mod p$ and start taking square roots or as Robin noted (in a now deleted comment/his answer) guess an $a$ and see if $z = a^{\frac{p-1}{4}}$ is a candidate.

There have been improvements to this by Brillhart. Derek was kind enough to hunt down the reference to that here: http://jstor.org/pss/2005889 and a free link (thanks to J.M) http://www.ams.org/journals/mcom/1972-26-120/S0025-5718-1972-0314745-6/S0025-5718-1972-0314745-6.pdf

I must add Gauss's construction, which I learned about from p.64 of Stark's "An Introduction to Number Theory". If by "best", you mean computationally most efficient, then certainly the following is not the best way. But there are other ways to measure quality...

Quoting Stark:

"In 1825 Gauss gave the following construction for writing a prime congruent to $1 \pmod{4}$ as a sum of two squares: Let $p=4k+1$ be a prime number. Determine $x$ (this is uniquely possible...) so that

$$ x = \frac{(2k)!}{2(k!)^2} \pmod{p}, \quad |x| < \frac{p}{2}$$

Now determine $y$ so that

$$ y = x \cdot (2k)! \pmod{p}, \quad |y| < \frac{p}{2}$$

Gauss showed that $x^2+y^2=p$."