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

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 4

One may quickly compute a representation of a prime $\rm\:p\equiv 1\ (mod\ 4)\:$ as a sum of two squares by employing the Euclidean GCD algorithm in $\rm\mathbb Z[\:i\:]\:$ and an algorithm for computing square roots $\rm\:(mod\ p)\:$.

THEOREM $\ $ Let $\rm\ \: c = \sqrt{-1}\ \:(mod\ p)\ $ and $\rm\ gcd(p,\:i-c)\: = \: a+b\:i\ $.$\ $ Then $\rm\ p = a^2 + b^2\:$.

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

(1) $\:$ If $\rm\:\alpha|p\ $ properly then, conjugating, $\rm\:\alpha'|p\ $ properly. Multiplying them yields $\rm\:\alpha\:\alpha'|p^2\:$ properly in $\:\mathbb Z\:$. But in $\rm\mathbb{Z}\:$ the only proper factor$>0$ of $\rm\:p^2\:$ is $\rm\: p\:$, thus $\rm\: p = \alpha\:\alpha' = (a+b\:i)(a-b\:i) = a^2 + b^2\:$.

(2) $\ $ If $\rm\ p\ \gamma \:=\: \alpha\:\beta\ $ and $\rm\:p\:$ does not divide $\rm\:\alpha\:$ nor $\rm\:\beta\:$,$\ $ then $\rm\: gcd(\alpha,p)\:$ is a $\:$ proper$\ $ factor of $\rm\:p\ $.
Else $\rm\ gcd(\alpha,p) = p\:$ or $1$. If the gcd is $\rm p$ then $\rm\: p\:|\:\alpha\:$ contra hypothesis. Otherwise, if $\rm\ gcd(\alpha,p)=1\ $ then by Euclid's Lemma, $\rm\ \alpha\:|\:p\:\gamma \ \Rightarrow\ \alpha\:|\:\gamma\:$,$\ $ so $\rm\ \gamma/\alpha = \beta/p\ \Rightarrow\ p\:|\:\beta\:$,$\ $ again contra hypothesis. Note: generally, in rings, GCDs are unique only up to unit multiples. Here the units are $\rm\ i^n = \pm 1,\:\pm i\:$.

(3) $\rm\ \ x^2 + 1 \ =\ (x-c)\:(x+c) + p\ f(x)\ $ in $\rm\mathbb\ Z[x]\ \: \Rightarrow\: -p\ f(i)\ =\ (i-c)\:(i+c)\ $ in $\rm\mathbb\ Z[\:i\:]\ \ $ by evaluation at $\rm\: x = i\:$. This splitting is suitable to split $\rm\:p\:$ by (2) since $\rm\:p\:$ doesn't divide $\rm\: i\pm c\ \:$ in $\rm\ \: Z[\:i\:]\:$.

REMARK $\ $ There are many variations on the Euclidean algorithm in $\rm\mathbb \ Z[\:i\:]\ $ that are employed in practice, e.g. employing continued fractions, binary quadratic forms, etc. Also, there are also at least a few algorithms for computing sqrts $\rm\:(mod\ p)\:$, e.g. by factoring polynomials $\rm\:(mod\ p)\:$, by elliptic curves (Schoof), and and algorithm of Tonelli and Shanks. For much further information see Henri Cohen's book "A course in computational algebraic number theory". Here is an excerpt on the beautiful algorithm of Cornacchia:

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 $\rm\ p = a^2 + b^2\ $ for an integer $\rm\ p = 1\ \: (mod\ 4)\:.\ \ $ The algorithm is coded in $\ $ $\rm Macsyma$ $\ $ but it can easily be translated into almost any another language. The only Macsyma specific subroutines that you might need to implement in another language are exponentiation $\rm (mod\ p)$ (by repeated squaring) and the floor function in $\rm \mathbb Q[\:i\:]\:$. While this algorithm is reasonably efficient, the point is not to optimize efficiency but, rather, to emphasize the simplicity of the solution employing the Euclidean algorithm in $\rm\ \mathbb Z[\:i\:]\:.\ $ For the underlying theory see my companion post. 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 = 1\ \: (mod\ 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\quad\quad\quad\quad\quad\quad\quad\quad\quad\ \ \text{do}\ \: \text{ j : rat(a)^e ),}$

$\rm\quad \text{x : p,}\ \ \: \text{y : j - %i,} \quad\quad\quad\quad\quad\quad\quad $ /* 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 )) }$