# Efficiently finding two squares which sum to a prime

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 ?

-

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

-
+1: If you would care to expand on Hermite-Serret too, I will delete my answer. –  Aryabhata Oct 2 '10 at 16:19
Don't do that, Moron. It's rather tedious to prove it in detail but if you do some examples you see that the same numbers come out as remainders in H-S as the real parts of remainders in the Gaussian gcd calculation –  Robin Chapman Oct 2 '10 at 16:20
Thanks especially for the python. I have to admit I actually went on to solve the problem which originally prompted this question by simply creating a big array to hold pairs up to the max value m of interest, iterating p<-(1 to sqrt(m)),q<-(1 to p) and storing (p,q) in the appropriate entry for p^2+q^2, and then just looking up the pairs for values as needed. Of course that approach only scales up so far... it's nice to see the general solution. –  timday Oct 8 '10 at 23:56

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

-
Yes, just counting points on the CM elliptic curve $y^2=x^3-x$ mod $p$. :-) –  Robin Chapman Oct 2 '10 at 16:25
where does this formula come from? –  anon Oct 2 '10 at 17:01
@muad: Dear muad, My memory is that this formula dates back to the 1700's, but I forget to whom it is originally due. But from a modern viewpoint, as Robin notes, it is a special case of the general theory of counting points on CM elliptic curves. A comment box is not very conducive to explaining the details, but if you asked a question, I (and I'm sure others) could try to explain more. –  Matt E Oct 2 '10 at 19:58
@Matt E: I think it would be wise to add a note to your answer that explicitly emphasizes that this formula is worthless for "efficiently finding" the representation - which is the OP's question. Indeed, it's worst than a brute-force search. I fear some readers may not realize that at quick glance (it's not necessarily implied by your "for amusement" preamble). –  Bill Dubuque Oct 4 '10 at 5:07
@Bill: Dear Bill, I though that "Just for amusement" conveyed this fairly clearly, and with your comment (and my agreement here!) there is no room left for ambiguity. (Still, I will add a disclaimer in the answer itself, just as a fail-safe.) –  Matt E Oct 4 '10 at 5:31

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:

-
Surely this is as fast as any algorithm could be for this task? –  paul garrett Aug 25 '12 at 22:03

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

-
Here's the Brillhart reference: jstor.org/pss/2005889 Also, here's a very relevant an interesting article by Alf van der Poorten on this topic maths.mq.edu.au/~alf/SomeRecentPapers/155.pdf that I'm part way through thanks to this discussion. –  Derek Jennings Oct 2 '10 at 18:42
@Derek: Thanks! –  Aryabhata Oct 2 '10 at 20:41
Here's a freely accessible link to the Brillhart paper, from AMS: ams.org/journals/mcom/1972-26-120/S0025-5718-1972-0314745-6/… –  Ｊ. Ｍ. Oct 5 '10 at 12:07
@J.M: Added to the answer. Thanks. –  Aryabhata Oct 5 '10 at 13:54
@J.M: Thanks for that. I was interested to see what it said but had not yet got hold of a copy. –  Derek Jennings Oct 5 '10 at 14:13

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

-

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

-