# 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 ?

• do you know the prime and you are wondering which two squares add up to it? Oct 18, 2016 at 19:11
• @user25406: Yes that was the idea; you're about six years late to this question though. Oct 18, 2016 at 19:37
• here's the solution then. math.stackexchange.com/questions/1972771/… Oct 18, 2016 at 19:37
• better late than never but 6 years ago I did not know the answer but I do now. Oct 18, 2016 at 19:39
• One can use Gauss lattice reduction (equivalently, LLL algorithm in dimension 2) to find a short vector in the 2-dimensional lattice L generated by (p, 0) and (alpha, 0) where alpha^2 = -1 mod p. Jun 10 at 5:27

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. Oct 2, 2010 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 Oct 2, 2010 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. Oct 8, 2010 at 23:56
• @user477343: The above comments make more sense if you know that the user who made the first comment (currently using the name "Aryabhata") used to use the name "Moron". Jan 20, 2019 at 17:25
• @timday I had always been curious about the Moron... until this day... Jan 22, 2019 at 19:29

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$. :-) Oct 2, 2010 at 16:25
• where does this formula come from?
– anon
Oct 2, 2010 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. Oct 2, 2010 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). Oct 4, 2010 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.) Oct 4, 2010 at 5:31

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.   • Surely this is as fast as any algorithm could be for this task? Aug 25, 2012 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. Oct 2, 2010 at 18:42
• 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, 2010 at 12:07
• @J.M: Added to the answer. Thanks. Oct 5, 2010 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. Oct 5, 2010 at 14:13
• @Derek: FWIW, AMS has been very generous in giving free access to a substantial portion of the Mathematics of Computation archive. Enjoy! Oct 5, 2010 at 14:34

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