Lagrange’s theorem of 1770 — every natural number is a sum of four integer squares — is one of the cleanest examples in number theory of the gap between beauty and machinery. The statement is one line. The classical proof bolts together two tools that look completely different from each other: Euler’s algebraic identity, and Fermat’s method of infinite descent. This note traces how those two tools survive when the proof is moved into Lean 4 and Mathlib.

1. The classical proof

The statement is

Requiring $a,b,c,d \in \mathbb{N}$ is equivalent, since squares eat negative signs.

The first pillar of the proof is Euler’s four-square identity:

This identity says that the property “sum of four squares” is closed under multiplication. It reduces the whole problem to the prime case: if every prime is a sum of four squares, then so is every positive integer, by prime factorisation and the identity above.

It remains to prove the prime case. For $p=2$ we have $2 = 1^2 + 1^2 + 0^2 + 0^2$. For an odd prime $p$, we first prove an “initial multiple” lemma: there exist integers $a, b$ and $m$ with

The construction is the usual pigeonhole over $\mathbb{F}_p$. Consider

Each set has $(p+1)/2$ elements; their disjoint union would have $p+1$ elements, which cannot fit in $\mathbb{F}_p$. So $A$ and $B$ intersect, giving $a, b$ with $a^2 \equiv -1 - b^2 \pmod p$, i.e. $p \mid a^2 + b^2 + 1$. Choosing $a, b$ as minimal-absolute-value representatives gives $0 < m < p$.

The second step is Fermat’s descent. Let $m$ be the smallest positive integer such that $mp = a^2 + b^2 + c^2 + d^2$ has a solution. The previous step shows that such an $m$ exists and is less than $p$. If $m = 1$, we have the prime case. Assume $m > 1$; we derive a contradiction by producing a smaller witness with $r < m$.

If $m$ is even, the argument is short. From $a^2+b^2+c^2+d^2 = 2M$ we can pair $a,b,c,d$ by parity; for same-parity $x, y$,

so $M$ is again a sum of four squares. Applied to $mp = 2 \cdot (m/2) \cdot p$, this gives a representation of $(m/2)\, p$, contradicting the minimality of $m$.

If $m$ is odd, replace $a, b, c, d$ by their minimal-absolute-value representatives modulo $m$, say $A, B, C, D$, with $|A|, |B|, |C|, |D| < m/2$. Since the original sum is divisible by $m$,

for some non-negative integer $r$. The absolute-value bound gives $0 \le r < m$. We rule out $r = 0$: if $A = B = C = D = 0$, then $m$ divides all of $a, b, c, d$, so $m^2 \mid mp$, hence $m \mid p$, contradicting $1 < m < p$. Therefore $0 < r < m$.

Now Euler’s identity returns. Multiplying $A^2+B^2+C^2+D^2 = rm$ and the original $a^2+b^2+c^2+d^2 = mp$ produces a sum-of-four-squares representation of $m^2 (rp)$, with the four linear combinations chosen so that each is divisible by $m$. Dividing through by $m^2$ gives

with $0 < r < m$. This contradicts the minimality of $m$. Hence $m = 1$, and the odd prime $p$ is a sum of four squares.

2. The formalization in Lean 4 / Mathlib

The corresponding file in Mathlib is Mathlib.NumberTheory.SumFourSquares. The main theorem is

theorem Nat.sum_four_squares (n : ℕ) :
    ∃ a b c d : ℕ, a^2 + b^2 + c^2 + d^2 = n

and the prime version is

protected theorem Nat.Prime.sum_four_squares {p : ℕ} (hp : p.Prime) :
    ∃ a b c d : ℕ, a^2 + b^2 + c^2 + d^2 = p

Euler’s identity lives as a separate tool: euler_four_squares holds over an arbitrary commutative ring, and the natural-number version Nat.euler_four_squares uses Int.natAbs to pull the possibly-signed intermediate expressions back into $\mathbb{N}$.

A typical use looks like this:

import Mathlib.NumberTheory.SumFourSquares

open scoped BigOperators

def FourSq (n : ℕ) : Prop :=
  ∃ a b c d : ℕ, a^2 + b^2 + c^2 + d^2 = n

example (n : ℕ) : FourSq n :=
  Nat.sum_four_squares n

example {p : ℕ} (hp : Nat.Prime p) : FourSq p :=
  hp.sum_four_squares

-- Closure under multiplication via Euler's identity.
theorem FourSq.mul {m n : ℕ} (hm : FourSq m) (hn : FourSq n) :
    FourSq (m * n) := by
  rcases hm with ⟨a, b, c, d, hm⟩
  rcases hn with ⟨x, y, z, w, hn⟩
  refine ⟨
    ((a : ℤ) * x - b * y - c * z - d * w).natAbs,
    ((a : ℤ) * y + b * x + c * w - d * z).natAbs,
    ((a : ℤ) * z - b * w + c * x + d * y).natAbs,
    ((a : ℤ) * w + b * z - c * y + d * x).natAbs,
    ?_⟩
  calc
    _ = (a^2 + b^2 + c^2 + d^2) * (x^2 + y^2 + z^2 + w^2) := by
        simpa using Nat.euler_four_squares a b c d x y z w
    _ = m * n := by rw [hm, hn]

Mathlib’s main proof does not iterate over the prime factorisation explicitly. It uses Nat.recOnMul: if you prove $P(0)$, $P(1)$, $P(p)$ for every prime $p$, and the multiplicative step $P(m) \wedge P(n) \Rightarrow P(mn)$, then $P$ holds for every natural number. The source code’s final branches are exactly the zero / one / prime / mul cases, with the multiplicative branch calling euler_four_squares. The textbook phrase “by prime factorisation and closure” becomes an induction principle.

3. Three real obstacles in the formalization

Granularity of the existence step. The textbook says “two sets of size $(p+1)/2$ in $\mathbb{F}_p$ must intersect” and moves on. Writing that out in Lean from first principles requires Fintype.card, Finset.image, injectivity of squaring on representatives, and pigeonhole, all together. Mathlib avoids the rewrite by calling existing finite-field facts: ZMod.sq_add_sq gives a sum-of-two-squares representation in $\mathbb{Z}/p\mathbb{Z}$, and Nat.sq_add_sq_zmodEq lifts it to a natural-number statement with size estimates on the representatives.

Encoding descent. Lean refuses to accept the phrase “by infinite descent, contradiction”. The proof has to be reorganised around “take the minimal counterexample”. Mathlib uses Nat.findX to pin down the minimal $m$, then constructs a strictly smaller $r < m$ witness, feeds it to the minimality hypothesis, and lands the contradiction. The structure mirrors the textbook proof exactly, but every “there exists a minimum” and “strictly less than” has to be made explicit.

Switching between $\mathbb{N}$ and $\mathbb{Z}$. Euler’s identity’s four linear combinations naturally live in $\mathbb{Z}$ because of the minus signs, while the target theorem lives in $\mathbb{N}$. Mathlib’s resolution is to use .natAbs^2, with proofs full of push_cast, norm_cast, sq_abs, and mod_cast shuffling values back and forth. The even case adds a further wrinkle: dividing by $2$ requires Int.sq_add_sq_of_two_mul_sq_add_sq, which turns $2m = x^2 + y^2$ into $m = ((x-y)/2)^2 + ((x+y)/2)^2$ — and a parity lemma is needed first to guarantee the division is integral.

What the experience feels like in practice: each one-line step of the paper proof — “an existence in a finite field”, “by minimality”, “by parity”, “by squaring” — expands in Lean into a named lemma you have to point to, but every named lemma lands. The two tools — Euler’s identity managing multiplication, Fermat’s descent managing the prime case — are still visible inside the formalisation. Lagrange’s original 1770 structure, two and a half centuries later, is still the canonical path along which this theorem is taught, proved, and machine-verified.

积土成山，风雨兴焉；积水成渊，蛟龙生焉。——荀子《劝学》

Pile up earth into a mountain, and wind and rain rise from it; pool water into an abyss, and dragons breed there.