Reading the mathematical literature

In 1990 Don Zagier published in the American Mathematical Monthly a note whose entire proof is a single sentence. It proves Fermat's theorem that every prime $p \equiv 1 \pmod 4$ is a sum of two squares:

The involution on the finite set $S = \{(x, y, z) \in \N^3 : x^2 + 4yz = p\}$ defined by $(x, y, z) \mapsto (x + 2z,\ z,\ y - x - z)$ if $x \lt y - z$; $(2y - x,\ y,\ x - y + z)$ if $y - z \lt x \lt 2y$; $(x - 2y,\ x - y + z,\ y)$ if $x \gt 2y$ has exactly one fixed point, so $|S|$ is odd and the involution defined by $(x, y, z) \mapsto (x, z, y)$ also has a fixed point.

Every word is correct, and almost nothing is explained. Why is $S$ finite? Why does the map send $S$ to $S$, and why is it an involution? What about the cases $x = y - z$ and $x = 2y$, which the sentence skips? Which fixed point, and why only one? Why does a fixed point of the second involution give two squares? Where is it used that $p$ is prime, and that $p \equiv 1 \pmod 4$?

Research mathematics is written like this — not always so tersely, but always for readers who will supply routine steps. Reading it is a skill: finding the right sources among millions of documents, reading them in the right order, reconstructing what was left out, and knowing at every point whether a statement is proved, cited, conjectured or hoped. This lesson develops that skill, with Zagier's sentence as a fully worked case study (checked by computer and in Lean).

Nobody reads a research paper once, from the first line to the last. Experienced readers make several passes:

1. Survey. Title, abstract, introduction, the statements of the main theorems, the section headings, the bibliography. What is claimed? Is it new? Is this the paper I need?
2. Map. Read the statements of all lemmas and propositions, and draw the dependency graph: which results are used to prove which, which are proved here and which are cited. Identify the one or two key ideas.
3. Read the proofs that matter, with pencil and paper: check examples and special cases, reconstruct omitted steps, and mark every place a hypothesis is used.
4. Place it. Follow the references backwards (where do the tools come from?) and the citations forwards (who used, corrected or improved this?).

Getting stuck is normal, not a sign of inadequacy: "it is easy to see" may hide an hour of work, and occasionally an error. The aim of a first reading is understanding the structure; full verification comes later, and only for what you need. Reading is the other side of mathematical writing: the conventions a good writer follows are the ones a good reader relies on, and checks.

Start from expository sources. For a field you do not know, read a survey, a good textbook, the relevant article of the Princeton Companion to Mathematics, or a talk in the proceedings of the International Congress of Mathematicians. Surveys in the Bulletin of the AMS and the Séminaire Bourbaki are written to be read by non-specialists.

Databases. The arXiv (since 1991) holds preprints in all areas of mathematics, freely available; it moderates but does not referee. zbMATH Open — the former Zentralblatt für Mathematik, freely accessible since 2021 — indexes the published literature with reviews and summaries, reaching back to the nineteenth century through the Jahrbuch über die Fortschritte der Mathematik; MathSciNet (the American Mathematical Society's Mathematical Reviews, since 1940) does the same by subscription. Both classify papers by the Mathematics Subject Classification (MSC2020), so that, for example, 11A41 collects papers on primes; both link authors' publication lists, references and citations.

Following citations. Backward citations (the bibliography) lead to the tools and the history; forward citations ("cited by", in zbMATH, MathSciNet or general search engines) lead to later corrections, simplifications and generalisations — essential before relying on an old result. Errata and retractions are attached to papers in the databases.

Specialised resources. The OEIS for integer sequences; the NIST Digital Library of Mathematical Functions; the LMFDB for $L$-functions and modular forms; the Stacks Project, whose results have permanent "tags" that can be cited; Mathlib's documentation for formalised mathematics; MathOverflow for research-level questions.

Languages. Important work exists in Russian, French, German, Chinese, Japanese and other languages. Many Russian journals have cover-to-cover English translations (for example Uspekhi Matematicheskikh Nauk as Russian Mathematical Surveys), and reviews in zbMATH and MathSciNet are in English; machine translation now makes the original accessible too. Mathematical notation travels well; terminology does not always — see the conventions below.

Here is the sentence, reconstructed as a sequence of lemmas. Throughout, $p$ is prime with $p \equiv 1 \pmod 4$, and (as in Zagier's note) $\N$ means the positive integers; the Lean formalisation uses $\N = \{0, 1, 2, \dots\}$, which gives the same set $S$ by Lemma 1.

So the three cases of the map cover $S$ completely — the first place where primality is used, silently.

The reconstruction also shows what the hypotheses buy. For $n \equiv 3 \pmod 4$ the set $S$ is empty, since $x^2 + 4yz \equiv 1 \pmod 4$ for odd $x$. For composite $n \equiv 1 \pmod 4$ there is one fixed point of $f$ for each factorisation $n = x(x + 4z)$, so $|S|$ can be even (for $n = 21$, which is not a sum of two squares) or odd (for $n = 45 = 3^2 + 6^2$). For squares the boundary cases occur and $f$ is not an involution of $S$. The first computation below checks all of this. Zagier's proof is itself a distillation of an argument of Heath-Brown (1984), inspired by Liouville — something a reader learns by following his references.

Annotate the claims. In the introduction of a typical paper you will meet sentences of every kind: "We prove that … (Theorem 1.1)" — proved here; "By a theorem of Siegel [12] …" — cited (check that the cited result really has the hypotheses needed); "Assuming the generalised Riemann hypothesis, we show … (Theorem 1.3)" — a proved implication, not a proved conclusion; "We conjecture that the constant $2$ is optimal" — conjectured; "A heuristic argument suggests …" — heuristic; "Numerical computations up to $10^{12}$ support …" — numerical, with its bound; "It is well known that …" — folklore, if no source is given. An abstract or a title may claim more than the theorems prove; the theorems are what count.

Read the quantifiers and constants. "For all sufficiently large $N$" may hide an ineffective constant, as in Vinogradov's 1937 theorem that every sufficiently large odd number is a sum of three primes — closing that gap for every odd number greater than $5$ took Helfgott's work of 2013, combining improved analysis with large computations. Implied constants in $O(\cdot)$ and $\ll$ may or may not be computable; the paper should say which.

Know the status of preprints. A preprint's claims are claims. Perelman's 2002–2003 arXiv preprints on the Poincaré and geometrisation conjectures were never submitted to a journal; they became accepted mathematics after several teams (Kleiner and Lott; Morgan and Tian; Cao and Zhu) worked through every step. By contrast, Mochizuki's 2012 claimed proof of the abc conjecture, published in 2021 in a journal edited at his own institute, is not accepted by most experts: Scholze and Stix (2018) argued that a key step does not hold. Status is decided by the community's checking, not by the place of publication.

Know that published results can be wrong. Kempe's 1879 "proof" of the four colour theorem was accepted for eleven years until Heawood found the error. Grunwald's 1933 theorem in algebraic number theory was false; Shianghao Wang found the counterexample in 1948, and the corrected statement is now called the Grunwald–Wang theorem. Dulac's 1923 proof that a polynomial vector field has finitely many limit cycles contained a gap noticed by Ilyashenko around 1981; complete proofs were given by Écalle and by Ilyashenko in the early 1990s. Before relying on a result, look for later papers that cite it.

Know what evidence can do. Goldbach's conjecture has been verified for every even number up to $4 \times 10^{18}$ (Oliveira e Silva, Herzog and Pardi, 2014) and is unproved. The Mertens conjecture, $|M(n)| \lt \sqrt n$ for all $n \gt 1$, holds for every $n$ ever checked — and is false (Odlyzko and te Riele, 1985), although no counterexample is known explicitly. Numerical support makes a conjecture worth taking seriously; it does not make it a theorem. See experimental mathematics.

Mind the conventions. Is $0 \in \N$? (Zagier's $\N$ starts at $1$; Lean's at $0$.) In French, positif means $\ge 0$ and strictement positif means $\gt 0$. Bourbaki's "compact" includes the Hausdorff condition, which English-language texts often call "compact Hausdorff"; Bourbaki's "quasi-compact" is the English "compact". "Ring" may or may not require a unit; "$\log$" may be natural or base $2$; "smooth" may mean $C^\infty$ or "sufficiently differentiable". In formal libraries $a/0 = 0$ and natural-number subtraction is truncated, so identities carry side conditions that informal texts leave implicit. A statement is only meaningful once its conventions are fixed.

The Lean module carries out the reconstruction above, step by step. The set $S$ and Zagier's map, with the boundary cases assigned to branches (they never occur):

/-- The finite set `S = {(x, y, z) : x² + 4yz = p}`, written as a filter of the box
`{0, …, p}³` so that its finiteness is evident. -/
def zagierSet (p : ℕ) : Finset (ℕ × ℕ × ℕ) :=
  (Finset.range (p + 1) ×ˢ Finset.range (p + 1) ×ˢ Finset.range (p + 1)).filter
    fun t => t.1 ^ 2 + 4 * t.2.1 * t.2.2 = p

/-- Zagier's map, in his three cases:
`(x, y, z) ↦ (x + 2z, z, y - x - z)` if `x < y - z`;
`(2y - x, y, x - y + z)` if `y - z < x < 2y`; and `(x - 2y, x - y + z, y)` if `x > 2y`.
(The boundary cases `x = y - z` and `x = 2y` are assigned to the second and third branch;
they never occur for prime `p`.) -/
def zagier (t : ℕ × ℕ × ℕ) : ℕ × ℕ × ℕ :=
  if t.1 + t.2.2 < t.2.1 then (t.1 + 2 * t.2.2, t.2.2, t.2.1 - t.1 - t.2.2)
  else if t.1 < 2 * t.2.1 then (2 * t.2.1 - t.1, t.2.1, t.1 + t.2.2 - t.2.1)
  else (t.1 - 2 * t.2.1, t.1 + t.2.2 - t.2.1, t.2.1)

Lemma 1 — the silent uses of primality — and the consequence that the box bounds are automatic:

/-- **The omitted case analysis.** If `x² + 4yz = p` with `p` prime, then `x, y, z ≥ 1`, and
neither boundary case `x + z = y` nor `x = 2y` can occur: the first would make `p` a square,
`(x + 2z)²`, and the second would make `4 ∣ p`. This is where Zagier's sentence silently uses
that `p` is prime. -/
theorem mem_facts {p x y z : ℕ} (hp : p.Prime) (h : x ^ 2 + 4 * y * z = p) :
    1 ≤ x ∧ 1 ≤ y ∧ 1 ≤ z ∧ x + z ≠ y ∧ x ≠ 2 * y := by
  have hp2 := hp.two_le
  -- `p` is not a perfect square
  have nsq : ∀ m : ℕ, m ^ 2 ≠ p := by
    intro m hm
    have hdvd : m ∣ p := ⟨m, by rw [← hm]; ring⟩
    rcases hp.eq_one_or_self_of_dvd m hdvd with h1 | h1
    · subst h1; omega
    · subst h1
      have : m * m = m * 1 := by rw [← sq, mul_one]; exact hm
      have hm0 : m ≠ 0 := by omega
      have := Nat.eq_of_mul_eq_mul_left (Nat.pos_of_ne_zero hm0) this
      omega
  -- `4` does not divide `p`
  have n4 : ¬ 4 ∣ p := by
    intro h4
    have h2 : 2 ∣ p := Nat.dvd_trans ⟨2, rfl⟩ h4
    rcases hp.eq_one_or_self_of_dvd 2 h2 with h1 | h1
    · omega
    · rw [← h1] at h4; omega
  refine ⟨?_, ?_, ?_, ?_, ?_⟩
  · by_contra hx
    have hx0 : x = 0 := by omega
    subst hx0
    exact n4 ⟨y * z, by rw [← h]; ring⟩
  · by_contra hy
    have hy0 : y = 0 := by omega
    subst hy0
    exact nsq x (by simpa using h)
  · by_contra hz
    have hz0 : z = 0 := by omega
    subst hz0
    exact nsq x (by simpa using h)
  · intro hxy
    subst hxy
    exact nsq (x + 2 * z) (by rw [← h]; ring)
  · intro hx
    subst hx
    exact n4 ⟨y * y + y * z, by rw [← h]; ring⟩

/-- For prime `p`, membership in `S` is just the equation: the box bounds are automatic. -/
theorem mem_zagierSet {p : ℕ} (hp : p.Prime) (t : ℕ × ℕ × ℕ) :
    t ∈ zagierSet p ↔ t.1 ^ 2 + 4 * t.2.1 * t.2.2 = p := by
  constructor
  · intro h
    exact (Finset.mem_filter.1 h).2
  · intro h
    obtain ⟨hx, hy, hz, -, -⟩ := mem_facts hp h
    simp only [zagierSet, Finset.mem_filter, Finset.mem_product, Finset.mem_range]
    refine ⟨⟨?_, ?_, ?_⟩, h⟩
    · nlinarith
    · nlinarith
    · nlinarith

Lemmas 2–4:

/-- **The map preserves `x² + 4yz`.** In each case this is a polynomial identity, e.g.
`(x + 2z)² + 4z(y - x - z) = x² + 4yz`; the subtractions are legitimate by the case
conditions. -/
theorem zagier_mem {p : ℕ} (hp : p.Prime) {t : ℕ × ℕ × ℕ} (ht : t ∈ zagierSet p) :
    zagier t ∈ zagierSet p := by
  obtain ⟨x, y, z⟩ := t
  rw [mem_zagierSet hp] at ht ⊢
  dsimp only at ht
  by_cases h1 : x + z < y
  · rw [zagier_of_lt h1]
    dsimp only
    obtain ⟨k, rfl⟩ : ∃ k, y = x + z + k := ⟨y - x - z, by omega⟩
    rw [show x + z + k - x - z = k by omega, ← ht]
    ring
  by_cases h2 : x < 2 * y
  · rw [zagier_of_mid h1 h2]
    dsimp only
    rw [← ht]
    zify [show x ≤ 2 * y by omega, show y ≤ x + z by omega]
    ring
  · rw [zagier_of_ge h1 h2]
    dsimp only
    rw [← ht]
    zify [show 2 * y ≤ x by omega, show y ≤ x + z by omega]
    ring

/-- **Zagier's map is an involution on `S`.** The first case is sent to the third, the third
to the first, and the second to itself; once the boundary cases are excluded (`mem_facts`),
each verification is linear arithmetic. -/
theorem zagier_zagier {p : ℕ} (hp : p.Prime) {t : ℕ × ℕ × ℕ} (ht : t ∈ zagierSet p) :
    zagier (zagier t) = t := by
  obtain ⟨x, y, z⟩ := t
  rw [mem_zagierSet hp] at ht
  dsimp only at ht
  obtain ⟨hx, hy, hz, hb1, hb2⟩ := mem_facts hp ht
  by_cases h1 : x + z < y
  · rw [zagier_of_lt h1, zagier_of_ge]
    · simp only [Prod.mk.injEq, and_true]; omega
    all_goals omega
  by_cases h2 : x < 2 * y
  · rw [zagier_of_mid h1 h2, zagier_of_mid]
    · simp only [Prod.mk.injEq, true_and]; omega
    all_goals omega
  · rw [zagier_of_ge h1 h2, zagier_of_lt]
    · simp only [Prod.mk.injEq, true_and]; omega
    omega

/-- **The fixed points.** A point of `S` is fixed by Zagier's map iff it is `(1, 1, k)` where
`p = 4k + 1`. (In the first and third cases a fixed point would need `z = 0` or `y = 0`; in the
second it needs `x = y`, so `p = x(x + 4z)`, and primality forces `x = 1`.) -/
theorem zagier_fixed_iff {p : ℕ} (hp : p.Prime) {t : ℕ × ℕ × ℕ} (ht : t ∈ zagierSet p) :
    zagier t = t ↔ t = (1, 1, p / 4) := by
  obtain ⟨x, y, z⟩ := t
  rw [mem_zagierSet hp] at ht
  dsimp only at ht
  obtain ⟨hx, hy, hz, hb1, hb2⟩ := mem_facts hp ht
  constructor
  · intro hfix
    by_cases h1 : x + z < y
    · rw [zagier_of_lt h1] at hfix
      simp only [Prod.mk.injEq] at hfix; omega
    by_cases h2 : x < 2 * y
    · rw [zagier_of_mid h1 h2] at hfix
      simp only [Prod.mk.injEq] at hfix
      have hxy : x = y := by omega
      subst hxy
      have hdvd : x ∣ p := ⟨x + 4 * z, by rw [← ht]; ring⟩
      rcases hp.eq_one_or_self_of_dvd x hdvd with h1' | h1'
      · subst h1'
        simp only [Prod.mk.injEq, true_and]
        omega
      · -- `x = p` is impossible: then `p² ≤ x² + 4xz = p`
        have hle : p * p ≤ p := by
          calc p * p = x ^ 2 := by rw [h1']; ring
            _ ≤ x ^ 2 + 4 * x * z := Nat.le_add_right _ _
            _ = p := ht
        have := hp.two_le
        nlinarith
    · rw [zagier_of_ge h1 h2] at hfix
      simp only [Prod.mk.injEq] at hfix; omega
  · intro h
    simp only [Prod.mk.injEq] at h
    obtain ⟨rfl, rfl, rfl⟩ := h
    rw [zagier_of_mid]
    · simp only [Prod.mk.injEq, true_and]; omega
    all_goals omega

The two parity arguments, using Mathlib's Equiv.Perm.card_fixedPoints_modEq for Lemma 5 (stated there for any map whose $p^n$-th power is the identity), and the theorem:

/-- **`|S|` is odd** when `p ≡ 1 (mod 4)`: Zagier's map is an involution of the finite set `S`
with exactly one fixed point, and the number of fixed points of an involution of a finite set
has the same parity as the set (`Equiv.Perm.card_fixedPoints_modEq` with `2¹`). -/
theorem card_odd {p : ℕ} (hp : p.Prime) (h4 : p % 4 = 1) :
    Fintype.card (zagierSet p) % 2 = 1 := by
  let f : Function.End (zagierSet p) := fun t => ⟨zagier t.1, zagier_mem hp t.2⟩
  have hf : f ^ 2 ^ 1 = 1 := by
    funext t
    exact Subtype.ext (zagier_zagier hp t.2)
  have hmod := Equiv.Perm.card_fixedPoints_modEq (p := 2) (n := 1) hf
  have t0mem : ((1, 1, p / 4) : ℕ × ℕ × ℕ) ∈ zagierSet p := by
    rw [mem_zagierSet hp]
    show 1 ^ 2 + 4 * 1 * (p / 4) = p
    norm_num
    omega
  have hone : Fintype.card (Function.fixedPoints f) = 1 := by
    rw [Fintype.card_eq_one_iff]
    refine ⟨⟨⟨_, t0mem⟩, ?_⟩, ?_⟩
    · exact Subtype.ext ((zagier_fixed_iff hp t0mem).2 rfl)
    · rintro ⟨⟨t, htm⟩, hfix⟩
      have h' : zagier t = t := congrArg Subtype.val hfix
      exact Subtype.ext (Subtype.ext ((zagier_fixed_iff hp htm).1 h'))
  rw [hone] at hmod
  exact hmod

/-- **Fermat's two-squares theorem, by Zagier's proof, with every step reconstructed.** If `p` is
prime and `p ≡ 1 (mod 4)`, then `p = a² + b²`. Since `|S|` is odd, the involution
`(x, y, z) ↦ (x, z, y)` of `S` also has a fixed point `(x, y, y)`, and then
`p = x² + 4y² = x² + (2y)²`. -/
theorem zagier_two_squares {p : ℕ} (hp : p.Prime) (h4 : p % 4 = 1) :
    ∃ a b : ℕ, a ^ 2 + b ^ 2 = p := by
  have swap_mem : ∀ t ∈ zagierSet p, ((t.1, t.2.2, t.2.1) : ℕ × ℕ × ℕ) ∈ zagierSet p := by
    intro t ht
    rw [mem_zagierSet hp] at ht ⊢
    simp only
    rw [← ht]; ring
  let g : Function.End (zagierSet p) := fun t => ⟨(t.1.1, t.1.2.2, t.1.2.1), swap_mem _ t.2⟩
  have hg : g ^ 2 ^ 1 = 1 := by
    funext t
    exact Subtype.ext rfl
  have hmod := Equiv.Perm.card_fixedPoints_modEq (p := 2) (n := 1) hg
  have hodd := card_odd hp h4
  have hpos : 0 < Fintype.card (Function.fixedPoints g) := by
    have : Fintype.card (Function.fixedPoints g) % 2 = 1 := by
      rw [← hodd]; exact hmod.symm
    omega
  obtain ⟨⟨⟨⟨x, y, z⟩, hmem⟩, hfix⟩⟩ := Fintype.card_pos_iff.1 hpos
  have hyz : z = y := by
    have := congrArg (fun t => t.1.2.1) hfix
    simpa [g] using this
  rw [mem_zagierSet hp] at hmem
  simp only at hmem
  subst hyz
  exact ⟨x, 2 * z, by rw [← hmem]; ring⟩

The same theorem in Mathlib, proved through the arithmetic of the Gaussian integers — a different literature — and stated with the hypothesis $p \bmod 4 \ne 3$, which also covers $p = 2$:

/-- The same theorem as it appears in Mathlib (`Nat.Prime.sq_add_sq`), proved there through the
Gaussian integers — a different proof from a different literature. Note its hypothesis
`p % 4 ≠ 3`, which also covers `p = 2 = 1² + 1²`: the same theorem can be stated in slightly
different, non-equivalent ways, and a careful reader checks which one a source proves. -/
theorem sq_add_sq_library {p : ℕ} [Fact p.Prime] (hp : p % 4 ≠ 3) :
    ∃ a b : ℕ, a ^ 2 + b ^ 2 = p :=
  Nat.Prime.sq_add_sq hp

Two conventions that change the meaning of familiar formulas:

/-- In Lean (as in much of the formal literature) `a / 0 = 0`, so the "identity"
`(a / b) · b = a` is false without the hypothesis `b ≠ 0`, which informal texts often leave
implicit. Counterexample: `a = 1, b = 0`. -/
theorem div_mul_cancel_needs_hypothesis : ¬ ∀ a b : ℚ, a / b * b = a := by
  intro h
  have := h 1 0
  norm_num at this

/-- Subtraction of natural numbers is truncated (`3 - 5 = 0`), so `(a - b) + b = a` needs
`b ≤ a`. A text that works in `ℕ` and one that works in `ℤ` may write the same formula and mean
different things. -/
theorem nat_sub_add_needs_hypothesis : ¬ ∀ a b : ℕ, a - b + b = a := by
  intro h
  have := h 3 5
  omega

What is formalised: the whole of Zagier's argument, with every omitted step; only the parity lemma (Lemma 5) and the library version of the theorem are cited from Mathlib. Formalising a proof is the most thorough way of reading it — every case the author skipped must be supplied — and a formal proof of this length took longer to write than the sentence did to read. Not formalised: anything about the literature itself — the claims about databases, history and status in this lesson are ordinary scholarly claims, with references.

The first experiment tests our reading of Zagier's proof: it builds $S$ and both involutions and checks every claim of the sentence for all primes $p \equiv 1 \pmod 4$ below $3000$; then it runs the same code where a hypothesis fails, to see exactly what breaks.

"""Reading a terse proof with a computer: Zagier's one-sentence proof.

Zagier (1990) proved that every prime p = 1 (mod 4) is a sum of two squares
in one sentence: an explicit involution of the finite set
S = {(x, y, z) in positive integers : x^2 + 4yz = p} has exactly one fixed
point, so |S| is odd, so the involution (x, y, z) -> (x, z, y) has a fixed
point too, giving p = x^2 + (2y)^2. The experiment checks every claim of the
sentence for all primes p = 1 (mod 4) below 3000 and then runs the same
machinery on numbers for which a hypothesis fails - p = 3 (mod 4), composite
numbers and squares - to see exactly where each hypothesis is used. This is
exhaustive-finite evidence that we have read the proof correctly; the proof
itself is the human argument (and the lesson's Lean formalisation).
"""

import math

def zagier_set(n):
    """All (x, y, z) of positive integers with x^2 + 4yz = n."""
    out = []
    for x in range(1, math.isqrt(n) + 1):
        rest = n - x * x
        if rest <= 0 or rest % 4:
            continue
        m = rest // 4
        out.extend((x, y, m // y) for y in range(1, m + 1) if m % y == 0)
    return out

def zagier_map(t):
    """Zagier's map; None on the boundary x = y - z or x = 2y."""
    x, y, z = t
    if x < y - z:
        return (x + 2 * z, z, y - x - z)
    if y - z < x < 2 * y:
        return (2 * y - x, y, x - y + z)
    if x > 2 * y:
        return (x - 2 * y, x - y + z, y)
    return None

def analyse(n):
    """Facts about S for the number n (a dict)."""
    s = zagier_set(n)
    members = set(s)
    images = {t: zagier_map(t) for t in s}
    boundary = [t for t, u in images.items() if u is None]
    involution = not boundary and all(
        images[t] in members and images[images[t]] == t for t in s)
    fixed = sorted(t for t, u in images.items() if u == t)
    swap_fixed = sorted(t for t in s if t[1] == t[2])
    squares = sorted({tuple(sorted((x, 2 * y))) for x, y, _ in swap_fixed})
    return {"size": len(s), "boundary": len(boundary),
            "involution": involution, "fixed": fixed,
            "swap_fixed": swap_fixed, "squares": squares}

def is_prime(n):
    """Trial division (n is small here)."""
    return n >= 2 and all(n % d for d in range(2, math.isqrt(n) + 1))

def check_primes(limit):
    """Check all claims of the proof for primes p = 1 (mod 4) below limit."""
    primes = [p for p in range(5, limit, 4) if is_prime(p)]
    bad = []
    for p in primes:
        a = analyse(p)
        ok = (a["involution"] and a["fixed"] == [(1, 1, (p - 1) // 4)]
              and a["size"] % 2 == 1 and a["squares"])
        if not ok:
            bad.append(p)
    return len(primes), bad

def main():
    count, bad = check_primes(3000)
    print(f"Primes p = 1 (mod 4) below 3000: {count}; claims failing for: "
          f"{bad if bad else 'none'}")
    print("(involution of S, unique fixed point (1, 1, (p-1)/4), |S| odd,")
    print(" and a fixed point of the swap giving p = a^2 + b^2)")
    print()
    print(f"{'n':>4} {'kind':<12}{'|S|':>4} {'bnd':>4}  {'Zagier fixed pts':<28}"
          f"{'n = a^2 + b^2'}")
    cases = [(13, "prime 1 mod 4"), (97, "prime 1 mod 4"),
             (19, "prime 3 mod 4"), (21, "3 * 7"), (45, "3^2 * 5"),
             (65, "5 * 13"), (9, "square"), (25, "square")]
    for n, kind in cases:
        a = analyse(n)
        fixed = ", ".join(f"({x},{y},{z})" for x, y, z in a["fixed"]) or "-"
        sq = ", ".join(f"{u}^2+{v}^2" for u, v in a["squares"]) or "-"
        print(f"{n:>4} {kind:<12}{a['size']:>4} {a['boundary']:>4}  "
              f"{fixed:<28}{sq}")
    print()
    print("bnd = points on Zagier's boundary, where his map is undefined.")
    print("p = 3 mod 4: S is empty (x^2 + 4yz is 1 mod 4). Composites: one fixed")
    print("point per factorisation n = x(x + 4z), so |S| may be even (21, 65)")
    print("or odd (45). Squares: boundary points appear and the map is not an")
    print("involution of S. Primality is what makes the fixed point unique.")
    print()
    print("Epistemic status: exhaustive-finite checks of each step of the proof")
    print("for all primes below 3000, plus counterexamples showing that the")
    print("hypotheses are used. This tests our reading of the proof; it does not")
    print("replace the proof, which is checked in full in the Lean module.")

if __name__ == "__main__":
    main()

The second experiment computes the kind of evidence a paper might quote for two famous conjectures — the Mertens conjecture and Goldbach's conjecture. The data look equally reassuring. One conjecture is false and the other is open.

"""Two conjectures with impressive numerical support: one false, one open.

The Mertens conjecture says |M(n)| < sqrt(n) for all n > 1, where
M(n) = mu(1) + ... + mu(n) sums the Moebius function. Goldbach's conjecture
says every even n > 2 is a sum of two primes. For both, the experiment
computes the evidence a paper might quote: the largest |M(n)|/sqrt(n) for
n <= 10^6 (it stays far below 1), and the smallest number of Goldbach
representations in each decade up to 2 * 10^5 (it grows). Odlyzko and te
Riele proved in 1985 that the Mertens conjecture is false; Goldbach's
conjecture is open. The computations are exact and exhaustive up to the
stated bounds - and they cannot tell the two cases apart.
"""

import numpy as np

def mobius_sieve(limit):
    """mu(n) for 0 <= n <= limit (mu(0) set to 0), by a sieve over primes."""
    mu = np.ones(limit + 1, dtype=np.int64)
    mu[0] = 0
    is_prime = np.ones(limit + 1, dtype=bool)
    is_prime[:2] = False
    for p in range(2, limit + 1):
        if is_prime[p]:
            is_prime[2 * p::p] = False
            mu[p::p] *= -1
            mu[p * p::p * p] = 0
    return mu

def mertens(limit):
    """M(n) for 0 <= n <= limit."""
    return np.cumsum(mobius_sieve(limit))

def worst_mertens_ratio(m, lo, hi):
    """(n, |M(n)|/sqrt(n)) maximising the ratio over lo <= n <= hi."""
    n = np.arange(lo, hi + 1)
    ratio = np.abs(m[lo:hi + 1]) / np.sqrt(n)
    k = int(np.argmax(ratio))
    return lo + k, float(ratio[k])

def goldbach_counts(limit):
    """r[n] = #{(p, q) primes, p <= q, p + q = n} for 0 <= n <= limit."""
    is_prime = np.ones(limit + 1, dtype=bool)
    is_prime[:2] = False
    for p in range(2, int(limit ** 0.5) + 1):
        if is_prime[p]:
            is_prime[p * p::p] = False
    f = is_prime.astype(float)
    size = 1 << (2 * limit + 1).bit_length()
    ordered = np.rint(np.fft.irfft(np.fft.rfft(f, size) ** 2, size))
    ordered = ordered[: limit + 1].astype(np.int64)    # ordered pairs
    squares = np.zeros(limit + 1, dtype=np.int64)       # pairs with p = q
    evens = 2 * np.nonzero(is_prime[: limit // 2 + 1])[0]
    squares[evens] = 1
    return (ordered + squares) // 2

def decade_minima(r, limit):
    """Smallest r(n) over even n in [lo, hi) for decades up to limit."""
    rows, lo = [], 4
    while lo <= limit:
        hi = min(10 if lo == 4 else lo * 10, limit + 1)
        evens = np.arange(lo + lo % 2, hi, 2)
        rows.append((lo, hi - 1, int(r[evens].min())))
        lo = hi
    return rows

def main():
    limit = 10 ** 6
    m = mertens(limit)
    print(f"Mertens function M(n), n <= {limit}:")
    for k in range(1, 7):
        print(f"  M(10^{k}) = {int(m[10 ** k]):>4}")
    for lo in (2, 10 ** 4):
        n_star, worst = worst_mertens_ratio(m, lo, limit)
        print(f"  largest |M(n)|/sqrt(n) for {lo} <= n <= 10^6: {worst:.4f} "
              f"(at n = {n_star})")
    print()
    glimit = 2 * 10 ** 5
    r = goldbach_counts(glimit)
    print(f"Goldbach representations r(n), even 4 <= n <= {glimit}:")
    for lo, hi, least in decade_minima(r, glimit):
        print(f"  min r(n) for even n in [{lo}, {hi}]: {least}")
    print(f"  r(n) = 0 for some even n in range: {bool((r[4::2] == 0).any())}")
    print()
    print("Same shape of evidence: a quantity that stays safely inside its")
    print("bound, or grows. Yet |M(n)| > sqrt(n) for some n (Odlyzko and")
    print("te Riele, 1985), no explicit such n is known, and Goldbach's")
    print("conjecture, verified to 4 * 10^18, is still unproved.")
    print()
    print("Epistemic status: exact integer computations, exhaustive for the")
    print("stated ranges (the FFT counts are rounded from floating point and")
    print("cross-checked in the tests). They are evidence about those ranges")
    print("only; they cannot distinguish a true conjecture from a false one.")

if __name__ == "__main__":
    main()

Papers state implications, and readers often silently reverse them. Let $p$ stand for "$n$ is a prime with $n \equiv 1 \pmod 4$" and $q$ for "$n$ is a sum of two squares". Fermat's theorem is $p \to q$. Its converse $q \to p$ is false — $25 = 3^2 + 4^2$ is not prime — and the table shows the row where the two differ. Then type not q -> not p into the comparison box: the contrapositive is equivalent to the theorem, which is why proofs "by contrapositive" are valid. Remember also: "$p$ only if $q$" means $p \to q$; "$q$ is necessary for $p$" means $p \to q$; "$q$ is sufficient for $p$" means $q \to p$.

• Commentaries. Reconstructing what an author left out is an ancient practice. Liu Hui's commentary (263 CE) on the Chinese Nine Chapters on the Mathematical Art supplied justifications for its procedures; Bhāskara I's commentary (629 CE) explained the terse verses of Āryabhaṭa's Āryabhaṭīya; scholars writing in Arabic and Persian, among them Naṣīr al-Dīn al-Ṭūsī in the thirteenth century, produced critical editions of Euclid that filled gaps and examined the parallel postulate.
• Journals and reviews. The first scientific journals appeared in 1665 (the Journal des sçavans in Paris and the Philosophical Transactions in London); Crelle's Journal für die reine und angewandte Mathematik (1826) was among the first devoted to mathematics. The Jahrbuch über die Fortschritte der Mathematik reviewed the literature from 1868. Otto Neugebauer founded the Zentralblatt für Mathematik in 1931 and, after leaving Nazi Germany, Mathematical Reviews in the United States in 1940.
• Self-taught readers. Srinivasa Ramanujan learned much of his mathematics in India from G. S. Carr's Synopsis of Elementary Results, a compendium of thousands of results with few proofs; reconstructing them shaped his style.
• Preprints. Paul Ginsparg started the arXiv at Los Alamos in 1991, and mathematicians began posting there in 1992; today most new research in mathematics appears there first, free to readers everywhere.
• The theorem of this lesson. Stated by Girard (1625) and Fermat (1640), proved by Euler (1747–1749), recast by Lagrange and Gauss, and proved in one sentence by Zagier (1990), after Heath-Brown (1984).

• D. Zagier, "A one-sentence proof that every prime $p \equiv 1 \pmod 4$ is a sum of two squares", American Mathematical Monthly 97 (1990), p. 144.
• Martin Aigner and Günter M. Ziegler, Proofs from THE BOOK (6th ed., 2018) — several proofs of the two-squares theorem, including Zagier's.
• Timothy Gowers (ed.), The Princeton Companion to Mathematics (2008) — an entry point to every field.
• Imre Lakatos, Proofs and Refutations (1976) — how proofs, definitions and counterexamples evolve together.
• Nicholas J. Higham, Handbook of Writing for the Mathematical Sciences (3rd ed., 2020) — includes chapters on reading, searching and refereeing.
• Paul R. Halmos, "How to write mathematics", L'Enseignement Mathématique (1970); Steven G. Krantz, A Primer of Mathematical Writing (2nd ed., 2017) — the writer's side of the same contract. See also mathematical writing.