Efficiently finding two squares which sum to a prime

Question

The web is littered with any number of pages (example) giving an existence and uniqueness proof that a pair of squares can be found summing to primes congruent to 1 mod 4 (and also that there are no such pairs for primes congruent to 3 mod 4).

However, none of the stuff I've read on the topic offers any help with actually efficiently finding (ie other than a straight search up to sqrt(p)) the concrete values of such squares.

What's the best way to actually find them ?

Answer 1

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

Answer 2

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 ...").

Answer 3

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 a $\!\bmod p\,$ square root algorithm).

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 $\,c = \sqrt{-1}\pmod{\!p}.$

$ [\![1]\!]\ $ If $ \,\alpha\mid p\ $ properly, conjugating $\:\!\Rightarrow \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), Tonelli, 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.

Answer 4

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

Answer 5

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

Answer 6

Below is an algorithm to compute a representation $\, p = a^2 + b^2\, $ for an integer $\, p = 1\pmod{\!4}.\, $ The algorithm is coded in $\rm Macsyma$ but it can easily be translated into most any another language. The only Macsyma specific subroutines that you might need to implement in another language are exponentiation $\!\bmod\ p\,$ (e.g. by repeated squaring) and the floor function in $\,\mathbb Q[i].\,$ While this algorithm is reasonably efficient, the point is not to optimize efficiency but, rather, to emphasize the simplicity of a solution employing the Euclidean algorithm in $ \mathbb Z[i].\,$ For the underlying theory see my companion post here. For a faster algorithm see e.g. Collins, A fast Euclidean algorithm for Gaussian integers, and see also Agarwal and Frandsen, Binary GCD Like Algorithms for Some Complex Quadratic Rings.

/* $\:$ For $\rm\ p\equiv 1\pmod{\!4}\,$ compute $\rm\ \alpha = a + b\:i\ $ with $\rm \ p = \alpha\:\!\alpha' = a^2 + b^2\ \ $*/

$\rm\text {psplit(p) := block([}\ \:\text{j : 0,}\ \ \ \text{e : floor(p/4),}\ \ \ \text{algebraic : true},\ \ x,\ y\, ], $

$\rm\quad\text{block([ modulus : p ],}\quad\quad\quad\quad\ $ /* search for $\rm\ j = \sqrt{-1}\ \ (mod\ p)\ $ */

$\rm\quad\quad\ \text{for}\ \: \text{ a : 2}\ \ \text{while}\ \ \text{j^2 # -1}$

$\rm\qquad\qquad\qquad\quad\ \ \text{do}\ \: \text{ j : rat(a)^e ),}$

$\rm\quad \text{x : p,}\ \ \: \text{y : j - %i,} \qquad\qquad\qquad\ \ $ /* find $\rm\ gcd(p,\:j - i)\ \in\ \mathbb Z[\:i\:]\ $ by Euclidean algorithm */

$\rm\quad \text{do ( if}\ \text{ y = 0}\ \ \text{then}\ \text{ return(x), }$

$\rm\quad\quad\quad \text{j : rat(x - y * floor(rat(x/y))), }$

$\rm\quad\quad\quad \text{x : y,}\ \text{ y : j )) }$

Answer 7

I found an alternative way to obtain $ \gamma, \alpha, \beta $ such that $ p\gamma = \alpha \beta $, $ p \not \mid \alpha, \beta $ as in (2) in the comment above.

Consider numbers of the form $ (k^2 + 1)p $, where $ k \geq 0 $. I claim that there exists some $ k $ such that $ (k^2 + 1)p = A^2 + B^2 $ and $ A^2 $ is the largest square $ < (k^2 + 1)p $, i.e. $ A = \lfloor \sqrt{(k^2 + 1)p} \rfloor $.

Proof of claim: Suppose that $ p = a^2 + b^2 $, then $ (k^2 + 1)p = (ak + b)^2 - (a - bk)^2 $. The required condition is then $ \sqrt{p(k^2 + 1)} - 1 < ak + b $ or equivalently \begin{align} b^2 k^2 - 2a(b + 1)k + a^2 - 2b - 1 < 0. \end{align} One can check that the last inequality holds for some $ k \geq 0 $, so we are done.

However, I am not sure how efficient is this method compared to finding quadratic non residue.