Chapter 1
$n=x^2+y^2$ and Fermat's Theorem
First, we shall prove that every prime of the form $p = 4k+1$ can be written as a sum of two squares. The proof itself (attributed to Richard Dedekind) is decently short; however, there are a fair number of prerequisites to cover. That being said, the concepts used in this proof are important for solving other Diophantine equations!
Feel free to skip this section if you are already familiar with the basics of modular arithmetic.
Throughout the text, I will be using the symbol $a\ |\ b$ to denote that $a$ divides $b$, or equivalently that $b=ka$ for some number $k$ of the same type of $a$ and $b$ (i.e., if $a$ and $b$ are integers, then we require $k$ to be an integer).
Given any pair of integers $a$ and $b$, the largest positive integer dividing both $a$ and $b$ is called the greatest common divisor of $a$ and $b$ and is denoted by $\gcd(a, b)$. We call two integers coprime if $\gcd(a, b)=1$.
Before proceeding, we must define what is meant by an equivalence relation.
An equivalence relation $\sim$ over a set $S$ is a relationA set of tuples $(x,y)$ where $x,y\in S$ indicating two objects are related in some way.
If $S$ is the set of humans, then we could define a relation $x\sim y$ to indicate that $x$ as a parent of $y$. that is reflexiveFor all $x\in S$, $x\sim x$., symmetricIf $x\sim y$, then $y\sim x$., and transitiveIf $x\sim y$ and $y\sim z$, then $x\sim z$..
These properties make sense if we think of $\sim$ as meaning “equals”. We can think of equivalence relations as formalizing the idea of two mathematical objects being “similar”.
We use the symbol $a \equiv b\pmod{n}$ to denote the equivalence relation for modular congruence. It represents the idea that two integers $a$ and $b$ are “the same” if they have equal remainders after division by $n$. More specifically:
For instance, each element in $S=\{1, 5, 9, 13, \cdots\}$ is congruent to $1$ since they leave a remainder of $1$ if we divide by $4$. Since all the elements of $S$ are “the same” as $1$, we can let $1$ be a representative for the whole set. We call this representative a residue modulo $n$ and we call the set of all residues the integers modulo $n$, denoted by $\mathbb{Z}/n\mathbb{Z}$.
Since the residue is (effectively) an integer, we can do arithmetic with it. Please note that I’m skipping over a lot of formalism and rigor necessary to prove this, but, in short, we can add, subtract, and multiply residues so long as we take the residue of the result modulo $n$.
For instance, $4+5\equiv 9\equiv 1\cdot 7+2\equiv 2\pmod{7}$ and $-4\cdot 5\equiv -20\equiv -3\cdot 7+1\equiv 1\pmod{7}$.
Due to a result known as Bézout’s identity, if $x$ and $n$ are coprime, then there exist integers $y$ and $m$ such that
I am omitting the proof of this for brevity, but one proof involves running the Euclidean algorithm on $a$ and $b$ and constructing the identity by running all the steps in reverse.
If we apply modulo $n$ to both sides of Bézout’s identity, it follows that $xy\equiv 1\pmod{n}$. Hence, any $x\in \mathbb{Z}/n\mathbb{Z}$ has a multiplicative inverse (or reciprocal) if and only if $x$ and $n$ are coprime. This will be the primary way to divide residues, since we can “divide” by $x$ by instead multiplying by $y$.
Much like how we can add, subtract, multiply, and divide over $\mathbb{Z}/n\mathbb{Z}$, it is intuitive to ask if we can do other arithmetic operations, such as taking square roots. Understanding when an integer has a square root modulo $n$ will be crucial in showing when a prime can be written as a sum of squares.
An integer $a$ has a square root modulo $n$, also known as a quadratic redsidue modulo $n$, if there exists some $x\in \mathbb{Z}/n\mathbb{Z}$ such that
For now, we will only be considering prime-valued moduli, and henceforth I will be using $p$ in place of $n$.
Not every integer has a square root in the land of residues! For instance, there is no solution to $x^2\equiv 5\pmod{13}$ but there are two solutions to $x^2\equiv 10\pmod{13}$. We use the Legendre symbol to denote if an integer has a quadratic residue or not.
Let $a$ and $p$ be coprime integers with $p$ prime. $$\left(\dfrac{a}{p}\right) = \begin{cases} 1 & \text{if } x^2\equiv a\pmod{p} \text{ for some } x\in \mathbb{Z}\\ -1 & \text{otherwise} \end{cases}$$
In fact, we can explicitly compute the Legendre symbol for arbitrary $a$ and $p$ thanks to a formula attributed to Euler.
Let $a$ and $p$ be coprime integers with $p$ prime.
Below are a few examples of how we would use this criterion:
- $\left(\dfrac{2}{17}\right) \equiv 2^8 \equiv 256 \equiv 1\pmod{17}$ and indeed $6^2\equiv 2\pmod{17}$.
- $\left(\dfrac{3}{17}\right) \equiv 3^8 \equiv 6561 \equiv -1\pmod{17}$ and so there is no quadratic residue of $3$ modulo $17$.
- $\left(\dfrac{7}{53}\right) \equiv 7^{26} \equiv 1\pmod{53}$ and indeed $22^2\equiv 7\pmod{53}$.
Raising a number to a potentially large power like this may seem laborious; however, there is an efficient method that works for modular exponentiation. Unfortunately, while this criterion tells us if a solution exists, it doesn’t tell us the solutions. There are algorithms that can give us a solution (such as Tonelli-Shanks and Cipolla’s algorithm), but they are beyond the scope of this text.
Proving Euler’s criterion relies on another theorem attributed to Fermat, often referred to as Fermat’s Little Theorem (to avoid confusion with the similarly named Last Theorem).
I will offer two proofs of Fermat's Little Theorem: one algebraic, and the other combinatorial.
-
Recall the Binomial Theorem:
$\displaystyle (x+y)^n=\sum_{k=0}^{n}\binom{n}{k}x^ky^{n-k}$where $\displaystyle \binom{n}{k}$ indicates the binomial coefficientThe number of ways to choose $k$ items from a set of $n$; also denoted by $_n\text{C}_k$ or $C(n,k)$
$\displaystyle \binom{n}{k}=\dfrac{n!}{k!(n-k)!}$.. If we consider $(x+y)^p$ over the integers modulo $p$, we arrive at an interesting result, often referred to as the Freshman's Dream.$\displaystyle (x+y)^p\equiv x^p+y^p\pmod{p}$This identity follows from the fact that $\displaystyle \binom{p}{k}\equiv 0\pmod{p}$ except when $k=0$ or $k=p$.
$\displaystyle \binom{p}{k}\equiv \dfrac{p!}{k!(p-k)!}\equiv p\left(\dfrac{(p-1)!}{k!(p-k)!}\right)\equiv 0\pmod{p}$Since no factor of $k!$ or $(p-k)!$ divides the prime $p$ (unless $k=0$ or $k=p$), we know that the quantity inside the parentheses is an integer since $\displaystyle \binom{p}{k}$ is an integer. So, $p$ divides $\displaystyle \binom{p}{k}$. $\Box$
So, all terms in the binomial expansion are eliminated except the ones containing just $x$'s and just $y$'s.
Then $a^p\equiv ((a-1)+1)^p \equiv (a-1)^p+1$. Likewise, $(a-1)^p+1\equiv ((a-2)+1)^p+1\equiv (a-2)^p+2$ and it can be proved by induction that $a^p\equiv (a-k)^p+k$. If we let $k=a$, it follows that $a^p\equiv a\pmod{p}$. $\Box$
-
Suppose we have $a$ different types of colored beads ($\text{\textcolor{red}{red}, \textcolor{green}{green}, \textcolor{blue}{blue}}$, etc.) and that there are infinitely many beads of each color. Suppose we want to pick $p$ of these beads to put on a necklace. How many different necklaces can we create?
We can represent a necklace by listing each bead in clockwise order. One such necklace of length 6 with 3 different colors is $\text{\textcolor{red}{o}-\textcolor{red}{o}-\textcolor{blue}{o}-\textcolor{green}{o}-\textcolor{blue}{o}-\textcolor{red}{o}}$. We can rotate any necklace clockwise by shifting each element right (in our example, this yields $\text{\textcolor{red}{o}-\textcolor{red}{o}-\textcolor{red}{o}-\textcolor{blue}{o}-\textcolor{green}{o}-\textcolor{blue}{o}}$). I've prepared a Desmos graph to help you visualize these necklaces.
Notice, however, that this doesn't actually give us a distinct necklace. That is, necklaces are cyclic permutationsPermutations whose structure does not change if you shift its elements in a cycle..
So, how many different necklaces of length $p$ can we create? There are $a^p$ possible strings of beads we can create, and $a$ of them consist of only a single color of bead. If a necklace has only one color of bead, no matter how we shift the necklace, the permutation is unchanged. However, for any necklace with more than one color, it can be rotated counterclockwise exactly $p$ times to yield a new arrangement. So, we must avoid overcounting necklaces with more than one color.
Removing all duplicate necklaces, we find that there are $\dfrac{a^p-a}{p}$ distinct necklaces. Since this quantity must be an integer, we conclude that $a^p\equiv a\pmod{p}$. $\Box$
If $a\neq 0$, we may also write Fermat’s Little Theorem as $a^{p-1}\equiv 1\pmod{p}$, which is often a more useful formulation. Fermat’s Little Theorem is a crucial result in number theory and is routinely used to simplify large powers modulo $p$.
For instance, we may instantly verify that $7^{23843902143243}\equiv 2\pmod{11}$ since
Importantly for us, Fermat’s Little Theorem allows us to derive Euler’s Criterion!
If $a\not\equiv 0\pmod{p}$, then $a^{p-1}\equiv 1\pmod{p}$.
The only two solutions of $x^2\equiv 1\pmod{p}$ are $x\equiv \pm 1\pmod{p}$. This is because $(x-1)(x+1)\equiv 0\pmod{p}$, implying that one of $x-1$ or $x+1$ must be congruent to $0$ (or both, in the case of $p=2$). Hence, $x\equiv \pm 1\pmod{p}$. Notice that this implies that $a^{(p-1)/2}\equiv \pm 1\pmod{p}$.
If there exists an $x$ such that $x^2\equiv a\pmod{p}$, then $a^{(p-1)/2}\equiv (x^2)^{(p-1)/2}\equiv x^{p-1}\equiv 1\pmod{p}$. The contrapositive implies that if $a^{(p-1)/2}\equiv -1\pmod{p}$, then there is no quadratic residue of $a$.
Conversely, suppose that $a^{(p-1)/2}\equiv 1\pmod{p}$. To proceed, we need an important lemma: every integer $a$ may be written as $g^k\pmod{p}$ for some integers $g, k$. You'll have to take my word for this, as the proof relies on some group theory.
Since $a\equiv g^k\pmod{p}$, then $1\equiv a^{(p-1)/2}\equiv g^{k(p-1)/2}\pmod{p}$ implying that $(p-1)\ |\ k(p-1)/2$ by Fermat's Little Theorem. This means that $k$ must be even, hence $k=2k'$ and so $a\equiv (g^{k'})^2\pmod{p}$. Thus, $a$ is a square of some number modulo $p$. The contrapositive implies that if there does not exist an $x$ such that $x^2\equiv a\pmod{p}$, then $a^{(p-1)/2}\equiv -1\pmod{p}$.
Thus, all parts of Euler's Criterion have been proved! $\Box$.
What’s the point of this whole aside? I thought the point was to figure out when we can write a number as a sum of squares? Well, to answer this question, we must make use of an important lemma.
The condition $p\ |\ m^2+1$ is equivalent to saying that there exists some $m$ such that $m^2\equiv -1\pmod{p}$. Proving this is a direct application of Euler's criterion.
Euler's Criterion tells us that a solution $m$ exists only if $(-1)^{(p-1)/2}\equiv 1\pmod{p}$. This can only occur if
One of the fundamental tools that number theorists use to solve Diophantine equations is factoring. For instance, consider $p=x^2-y^2$, which is a similar Diophantine equation to the sum of squares equation we wish to solve. We can factor the right-hand side using the identity $x^2-y^2=(x-y)(x+y)$. Since the only divisors of the prime $p$ are either $1$ or $p$, we know that $x-y = 1$ and $x+y=p$ (since of course $x-y<x+y$).
Thus, we have the following linear system of equations:
\[\begin{aligned} x-y&=1\\ x+y&=p\\ \end{aligned}\]Then, the solution to this system is given by $x=\dfrac{p+1}{2}$ and $y=\dfrac{p-1}{2}$. Notice that this tells us that $p=x^2-y^2$ only has a solution for odd primes.
This is a straightforward example, but it illustrates an important point. Some Diophantine equations can be easily solved if we are able to completely factor one side of the equation.
Now consider our equation $p=x^2+y^2$. If we stick with just the integers, it’s not immediately clear how we could factor the right-hand side. But if we extend the integers to include a number $i$ satisfying $i^2=-1$, we can factor the right-hand side as
Hence, if we knew how $p$ factors over this new number system, we could use a similar approach as earlier to find solutions for $x$ and $y$.
This extension of the integers is known as the Gaussian integersNumbers of the form $x+iy$ where $x,y\in \mathbb{Z}$ and $i^2=-1$. and are denoted by $\mathbb{Z}[i]$. As I will only be explaining these numbers at a basic level, I recommend Keith Conrad’s notes on the subject if you wish to learn more.
Working with the Gaussian integers is a lot like working with the complex numbers. The rules for arithmetic are the same:
Division, however, is a bit more complicated since we cannot guarantee that the result will have integral real and imaginary parts. For instance, consider dividing $1-2i$ by $3+4i$:
One would expect the Gaussian integers not to always divide evenly, as the same is true for the integers.
We can easily visualize the Gaussian integers using the complex plane, letting the real component correspond to the $x$-axis and the imaginary component correspond to the $y$-axis. The Gaussian integers therefore represent the lattice points — the intersection of the grid lines. You can explore this representation using the Desmos graph below.
With this visualization, adding Gaussian integers corresponds to translation, multiplying Gaussian integers corresponds to dilation and/or rotation, and the solutions to the Diophantine equation $x^2+y^2=n^2$ correspond to the lattice points on a circle with radius $n$ centered at the origin. I will not go into further detail on proving these properties, but it is worth mentioning as they are invaluable interpretations for visualizing Gaussian integers and for doing analytic geometry.
Unlike with integers, if we have a pair of Gaussian integers $z$ and $w$, we cannot easily say which is “larger”. Is $10+2i > 1 - 12i$ merely because both its real and imaginary components are larger? To give a meaningful answer to this question, we can use a Gaussian integer’s normInformally, a integer-valued function which in some sense measures length.. We can think about the norm of $z=a+bi$ as a measure of $z$’s distance from the origin.
The norm of a Gaussian integer $z$ is given by $\mathcal{N}(z) = z\overline{z} = a^2+b^2$ and is guaranteed to be a non-negative integer. To ensure the norm’s output is an integer, we do not take the square root, as one would expect when finding the length of a vector. A crucial property of this norm is that it is multiplicative. That is, $\mathcal{N}(zw) = \mathcal{N}(z)\mathcal{N}(w)$ for any two $z,w\in \mathbb{Z}[i]$.
Let $z=a+bi$ and $w=c+di$. Then $zw = (ac-bd)+(ad+bc)i$. To prove $\mathcal{N}(zw)=\mathcal{N}(z)\mathcal{N}(w)$, we must show that $(ac-bd)^2+(ad+bc)^2=(a^2+b^2)(c^2+d^2)$. $$ \begin{align*} (ac-bd)^2+(ad+bc)^2 &= (a^2c^2-2abcd+b^2d^2)+(a^2d^2+2abcd+b^2c^2)\\ &= a^2c^2+b^2d^2+a^2d^2+b^2c^2\\ &= a^2(c^2+d^2)+b^2(c^2+d^2)\\ &= (a^2+b^2)(c^2+d^2) \end{align*} $$ Thus, $\mathcal{N}(zw)=\mathcal{N}(z)\mathcal{N}(w)$. $\Box$
It is this multiplicative property of norms that will allow us to factor Gaussian integers completely. If $z\in \mathbb{Z}[i]$ with $\mathcal{N}(z)=n=ab$ for positive integers $a, b, n$, then provided we can find $x, y\in \mathbb{Z}[i]$ such that $\mathcal{N}(x)=a$ and $\mathcal{N}(y)=b$, the hope is that $z=xy$, since \(\mathcal{N}(z)=\mathcal{N}(x)\mathcal{N}(y)=ab=n\) The challenge is therefore figuring out how to choose the correct $x, y$ for a given $a,b$. To answer that question, we must understand how divisibility works over the Gaussian integers. Much like how we can decompose any integers into a product of prime numbers, we can decompose any Gaussian integer into a product of smaller Gaussian “primes” — atoms in $\mathbb{Z}[i]$ that cannot be reduced further. But what are the Gaussian primes? To answer this, we’ll first need to introduce some terminology.
We call any $u\in \mathbb{Z}[i]$ with the property that $\mathcal{N}(u)=1$ a unit in $\mathbb{Z}[i]$. Which $z\in \mathbb{Z}[i]$ have unit norm? Well, we can show that $\mathcal{N}(z)=1$ if and only if $z$ is one of $1,-1,i,$ or $-i$. Moreover, these are the only invertible Gaussian integers — they are the only integers $z$ such that $1/z\in \mathbb{Z}[i]$.
Let $z\in \mathbb{Z}[i]$ where $z=a+bi$ for $a,b\in \mathbb{Z}$. For $\mathcal{N}(z)=1$, it follows that $a^2+b^2=1$ or $a^2=1-b^2$. This forces $b^2\leq 1$, hence $b$ is either $0$ or $\pm 1$. If $b=0$, then $a=\pm 1$. If $b=\pm 1$, then $a=0$. Thus, the only valid Gaussian integers with unit norm are $z=\pm 1, \pm i$. One can check that each of $U=\{\pm 1, \pm i\}$ has a multiplicative inverse that also belongs to $U$.
Suppose $z$ is invertible, meaning that there exists some $w\in \mathbb{Z}[i]$ such that $zw=1$. By applying the norm to both sides, we obtain:
More generally, algebraists define a unit of a ringA set where addition and multiplication can be adequately defined. to be any invertible element. Defining units in terms of invertibility helps us generalize the concept to rings where a norm cannot be defined (such as $\mathbb{Z}/n\mathbb{Z}$ where $n$ is composite).
We say that $z\in \mathbb{Z}[i], z\neq 0$ is an irreducible if and only if $z=xy$ implies that either $x$ or $y$ is unit (that is, $\mathcal{N}(x)=1$ or $\mathcal{N}(y)=1$). Equivalently, we can say that $z$ is irreducible if it cannot be factored into two Gaussian integers, each with smaller norm. Irreducibility is a generalization of the idea that a prime integer may only be written as the product of $1$ and itself.
There is another desirable property of prime integers worth considering: if $p$ is a prime integer and $p\ |\ ab$ for some integers $a, b$, then either $p\ |\ a$ or $p\ |\ b$.
Recall Bézout's identity. If $\gcd(a, p)=1$, then there exist integers $x$ and $y$ such that $ax+py=1$. If we multiply both sides of this equation by $b$, we find that $abx+bpy=b$. Since $ab$ is divisible by $p$, the entire left side of this equation is divisible by $p$. Thus, $p\ |\ b$.
If $\gcd(a, p)\neq 1$, then there exists some integer $d>1$ which divides $p$. But the only divisors of $p$ are $1$ and $p$ itself, meaning that $d=p$ and so $p\ |\ a$. Thus, $p\ |\ ab$ implies that either $p\ |\ a$ or $p\ |\ b$. $\Box$
We can easily extend this property of primality to Gaussian integers. We say that $z\in \mathbb{Z}[i],z\neq 0$ is a prime if and only if $z\ | ab$ implies that $z\ |\ a$ or $z\ |\ b$.
Why make a distinction between irreducibles and primes? Well, there are plenty of valid number systems where not every irreducible is a prime number. For instance, if we consider the set $\mathbb{Z}[-\sqrt{5}]=\{a+b\sqrt{-5}\ |\ a,b\in \mathbb{Z}\}$ which extends the integers by the irrational complex number $\sqrt{-5}$. It can be shown that $1+\sqrt{-5}$ divides $6=2\cdot 3$ but $1+\sqrt{-5}$ divides neither $2$ nor $3$. $\mathbb{Z}[-\sqrt{5}]$ is not a unique number system — we can prove similar results for $\mathbb{Z}[-\sqrt{7}]$, $\mathbb{Z}[-\sqrt{11}]$, etc.
Luckily for us, a Gaussian integer is irreducible if and only if it is prime. This means that, for our purposes, we can treat primes and irreducibles interchangeably. The same is true for the integers (as well as any Unique Factorization Domain).
Suppose that $z$ is prime and that $z=xy$ for some $x,y\in \mathbb{Z}[i]$. Then $z$ trivially divides $xy$ and so, since $z$ is prime, $z$ divides $x$ or $z$ divides $y$. Without loss of generality, we may assume that $z$ divides $x$. Then, there exists some $w\in \mathbb{Z}[i]$ such that $zw=x$. Substituting this into our definition for $x$ yields $z=zwy$ or $z(1-wy)=0$. Since $z$ is non-zero, it follows that $wy=1$, hence $y$ is unit and $z$ is irreducible.
Fully proving that irreducibles are prime is a much longer process. Instead, I will provide a proof sketch. The idea involves showing that every Gaussian integer can be written as a product of a unit and powers of irreducibles. This would imply that irreducibles are primes, since any $ab$ could be split into its constituent irreducibles. If $z\ |\ ab$ with $z$ irreducible, then $z$ must divide one of the irreducibles belonging to $a$ or belonging to $b$.
Suppose a Gaussian integer $z$ has a prime norm, that is, $\mathcal{N}(z)=p$. Then we can guarantee that $z$ itself is irreducible in $\mathbb{Z}[i]$, since if $z=xy$ for some $x, y\in \mathbb{Z}[i]$, then $\mathcal{N}(x)\mathcal{N}(y)=p$. This forces one of $x$ or $y$ to be unit and thus forces $z$ to be irreducible.
Using this fact, we can factor every Gaussian integer $z$ as a product of a unit and irreducibles via induction on the norm of $z$.
- For our base cases, if $\mathcal{N}(z)=1$, then $z$ is simply a unit. If $\mathcal{N}(z)=2$, then it can be shown that $z=u(1+i)$ where $u$ is a unit.
- Let the inductive hypothesis be that there exists some $z\in \mathbb{Z}[i]$ with $\mathcal{N}(z)< n$ which can be factored into irreducibles. Now consider a $z$ such that $\mathcal{N}(z)=n$. If $n$ is a prime number, then we know that $z$ itself must be irreducible. If instead $n$ is a composite number, then $n=ab$ and, since $a<n$ and $b<n$, we know by the inductive hypothesis that there exist $x, y\in \mathbb{Z}[i]$ such that $\mathcal{N}(x)=a$ and $\mathcal{N}(y)=b$ and $x,y$ may be written as products of irreducibles.
Thus, by induction, every Gaussian integer permits a factorization into irreducibles. Further work can be done to show that this factorization is unique (up to unit). $\Box$
Thus, $z$ is irreducible if and only if it is prime. $\Box$
There is much, much more that can be said about this topic. However, for the sake of brevity, we must carry on. I highly encourage you to consult the references if you wish to know more.
After all that work, we finally have all the tools we need to prove Fermat’s Theorem!
Recall our lemma from earlier: if $p\equiv 1\pmod{4}$, then $p\ |\ m^2+1$. Earlier, we considered this statement only over the integers, but now we are free to explore its ramifications over the Gaussian integers.
Suppose that $p$ is a prime over $\mathbb{Z}[i]$ and note that $m^2+1=(m+i)(m-i)$. Then, since $p\ |\ (m+i)(m-i)$, it follows that either $p\ |\ (m+i)$ or $p\ |\ (m-i)$. But $p>1$, and can neither divide the imaginary parts of $m+i$ nor $m-i$. So, $p\not|\ (m-i)$ and $p\not|\ (m+i)$. We have reached a contradiction; hence $p$ cannot be prime.
Since $p$ is not prime (and thus not irreducible) it can be written as a product of two Gaussian integers, neither of which are unit. Let $p=(x+yi)(w+zi)$ where $x,y,w,z\in \mathbb{Z}$. Since $p$ is real, it follows that $w+zi$ must be the conjugate of $x+yi$, hence $p=(x+yi)(x-yi)$. Expanding the right-hand side tells us that $p=x^2+y^2$. So, if $p\equiv 1\pmod{4}$, then $p$ may be written as the sum of two squares.
But what about the converse? Suppose that $p=x^2+y^2$. Then $p\equiv x^2+y^2\pmod{4}$. Since the possible values of $x^2\pmod{4}$ are $0$ and $1$, the possible values of $p\mod{4}$ are $0,1,2\pmod{4}$. Since $p$ is an odd prime, this guarantees that $p\equiv 1\pmod{4}$.
Thus, $p=x^2+y^2$ for integers $x$ and $y$ if and only if $p\equiv 1\pmod{4}$.
Notice what this says about factoring over the Gaussian integers: If $p\equiv 3\pmod{4}$, then $p\neq (x+iy)(x-iy)$ for any $x, y\in \mathbb{Z}$. In the case of $p=2$, we may write it as $2=1^2+1^2=(1+i)(1-i)$. Thus, $1+i$ (and its unit multiples), the irreducible factors of $p\equiv 1\pmod{4}$ (and their unit multiples), and all primes $p\equiv 3\pmod{4}$ (and their unit multiples) comprise the entire list of irreducibles over the Gaussian integers.
This gives us an algorithm for factoring any Gaussian integer. Suppose we have the Gaussian integer $z=165+111i$. Its norm is given by $\mathcal{N}(165+111i)=39546=2\cdot3^2\cdot 13^3$, so $z$ must factor into irreducibles with norms dividing these prime powers (i.e., we must find an irreducible $w$ whose norm is $2$, divides $3^2$, or divides $13^3$). We can choose $1+i$, $3$, and $2+3i$ as our set of irreducibles. We then need to choose to multiply by $\pm 1$ or $\pm i$ in order to match the signs on each component. After some calculation, we arrive at $165+111i=-3(1+i)(2+3i)^3$.
Now suppose we wanted to know if any given integer $n$ can be written as a sum of squares. Furthermore, what if we wanted to know how many solutions there were?
Let $n=\displaystyle 2^\alpha\prod_{i=1}^mp_i^{\beta_i}\prod_{i=1}^kq_i^{\gamma_i}$ where $\alpha,\beta_i,\gamma_i \geq 0$, $p_i\equiv 1\pmod{4}$, and $q_i\equiv 3\pmod{4}$.
$n$ may be written as a sum of squares if and only if $\gamma_i$ is even for all $i$. That is, $n$ can be written as a sum of squares if and only if every prime $q_i$ dividing $n$ has even order.
Furthermore, if a solution exists, the number of solutions to $n=x^2+y^2$ is given by $$C(n)=\displaystyle 4\prod_{i=1}^m(\beta_i+1)=4(\beta_1+1)(\beta_2+1)\cdots (\beta_m+1)$$
Let $\mathcal{N}(z)=n$ where $z=x+yi$. The number of solutions $(x, y)$ to $n=x^2+y^2$ can be determined by factoring $z$ into irreducibles over the Gaussian integers. But to do that, we'll first need to factor $n$.
We can use the factorization $n=\displaystyle 2^\alpha\prod_{i=1}^mp_i^{\beta_i}\prod_{i=1}^kq_i^{\gamma_i}$ where $\alpha,\beta_i,\gamma_i \geq 0$, $p_i\equiv 1\pmod{4}$ and $q_i\equiv 3\pmod{4}$, as this gives a unique representation for any $n$ which allows us to directly apply Fermat's Theorem on the Sum of Squares.
First, we may simply write $2=(1+i)(1-i)$. Then, for each $p_i$, since $p_i\equiv 1\pmod{4}$, we can find a unique factorization of $p_i$ as $(p_i^{(x)}+ip_i^{(y)})(p_i^{(x)}-ip_i^{(y)})$. For each $q_i$, though we cannot write $q_i$ itself as a sum of squares, note that $q_i^2$ and any even power of $q_i$ may be trivially written as a sum of squares (e.g., $q_i^2=q_i^2+0^2$). Furthermore, we can conclude via induction that there cannot exist any Gaussian integer $w$ such that $\mathcal{N}(w)=q_i^{2k+1}$ for some integer $k>0$. Any $n$ which permits solutions must have even orders for its primes $q_i$.
To form a solution $(x, y)$ to $n=x^2+y^2$, we may write write $n=\mathcal{N}(z)=z\overline{z}$. Our goal is to choose factors of $\mathcal{N}(z)$ to form $z$ such that the product of the remaining factors is $\overline{z}$. For instance, if we wanted to solve $$x^2+y^2=2880=3^2(1+i)^6(1-i)^6(2+i)(2-i)$$ then $z=3(1+i)^6(2+i)$ and $\overline{z}=3(1-i)^6(2-i)$, which corresponds to the solution $(x, y) = (24,-48)$.
Let the product of all the Gaussian integers we didn't choose for $z$ be $w$. The goal is to find the number of ways to choose divisors for $z$ so that $w=\overline{z}$. Let's look at how we should select each type of prime.
- If $\alpha>0$ (so $2\ |\ n$), then $z$ must contain all of the powers of either $(1+i)$ or $(1-i)$. This is because any Gaussian integer formed from $z$ by increasing/decreasing the power of $1+i$ (and correspondingly decreasing/increasing the power of $1-i$) must take the form: $$z'=\dfrac{1\pm i}{1\mp i}z=\dfrac{(1\pm i)^2}{2}z=\pm i z$$ Since $z=x+yi$ and corresponds to the solution $(x, y)$, this means that $z'=\pm(-y+xi)$ and corresponds to the solutions $(-y, x)$ and $(y, -x)$. We will revisit this property of multiplying by units to yield alternate solutions later.
-
For each $p_i$, we can write $p_i=(p_i^{(x)}+ip_i^{(y)})(p_i^{(x)}-ip_i^{(y)})$ and choose $(p_i^{(x)}+ip_i^{(y)})^k(p_i^{(x)}-ip_i^{(y)})^{\beta_i-k}$ to belong to $z$ where $k\leq \beta_i$. This is the only way to ensure that the terms that go to $w$ will be the conjugate of the terms that go to $z$. This is because, if we sent $k$ powers of $p_i^{(x)}+ip_i^{(y)}$ and $l$ powers of $p_i^{(x)}-ip_i^{(y)}$ to $z$, $w$ gets $\beta_i-k$ of the former and $\beta_i-l$ of the latter. Then, in order for $w=\overline{z}$, $k=\beta_i-l$ and $l=\beta_i-k$.
So, $(p_i^{(x)}+ip_i^{(y)})^k(p_i^{(x)}-ip_i^{(y)})^{\beta_i-k}$ is a factor of $z$ and the number of possible factors of this type is $k=\beta_i+1$.
- Each $q_i$'s order is even, hence we must give $\gamma_i/2$ powers of $q_i$ to $z$. There is no other choice we can make, so each $q_i$ has no effect on the total number of solutions.
-
The final thing to consider is alternate factorizations of $2$, $p_i$, and $q_i$. For our purposes, $p$ factors uniquely since, if $p_i=((p_i')^{(x)}+i(p_i')^{(y)})((p_i')^{(x)}-i(p_i')^{(y)})$, then $(p_i')^{(x)}+i(p_i')^{(y)}$ must be a unit multiple of $p_i^{(x)}+ip_i^{(y)}$. The same is true for $2$, since $2=(-1+i)(-1-i)$ and $-1+i=i(1+i)$. For $q_i$, since its order is even, we can choose to use $\pm q_i$.
So, using unit multiples of each $2, p_i,$ and $q_i$ has the net effect of multiplying $z$ by $\pm 1$ or $\pm i$.
Combining everything together, we get $$C(n)=\displaystyle 4\prod_{i=1}^m(\alpha_i+1)=4(\alpha_1+1)(\alpha_2+1)\cdots (\alpha_m+1)$$ corresponding to four choices of what to multiply $z$ by (either $\pm 1$ or $\pm i$) and $\alpha_i+1$ distributions of factors of $p_i$ for each $i$. $\Box$
That was quite a lot of work! But what’s the takeaway? Well, given any integer $n$, we now have an algorithm for finding the number of ways to write $n$ as a sum of squares.
Suppose $n=50!^2$ — a ludicrously large number with 129 digits. Brute-forcing solutions with a computer is infeasible, but if we use our algorithm, it’s pretty straightforward. Using Legendre’s formula, we can factor $50!$ as
\[50!^2=2^{94}\cdot3^{44}\cdot5^{24}\cdot7^{16}\cdot11^{8}\cdot13^{6}\cdot17^{4}\cdot19^{4}\cdot23^{4}\cdot29^2\cdot31^2\cdot37^2\cdot41^2\cdot43^2\cdot47^2\]Note that every prime dividing $n$ which is congruent to $3\pmod{4}$ has even order, so solutions to $x^2+y^2=50!^2$ exist. The primes dividing $n$ which are congruent to $1\pmod{4}$ are $5,13,17,29,37$ and $41$. Thus,
\[C(50!^2)=4(24+1)(6+1)(4+1)(2+1)(2+1)=31500\]As one last observation, for any integer $n$, notice that $C(n^2)$ counts the number of lattice points on a circle centered at the origin with radius $n$. This fact was used by YouTuber 3Blue1Brown to give an alternative proof of a famous formula for $\pi$. I only remembered this video existed after writing this whole explainer, so consider this text to be a somewhat more formal take on the same ideas presented in 3B1B’s video.
That’s it for $n=x^2+y^2$! But why stop there? How might we solve equations of the form $n=x^2+my^2$ for arbitrary $m$? This question and many others will be covered in the next chapter.