Computational reproducibility

MotivationH

A result someone else can check

Amara in Accra, Minh in Hanoi and Sophie in Lyon run the same notebook from a published paper. Amara's and Minh's tables agree with the paper to sixteen digits; Sophie's differs in the ninth digit of one entry, and her Monte Carlo figure looks slightly different. None of them has made a mistake. Different library versions summed the same numbers in a different order, and the random generator was never seeded. Is the paper's claim wrong? Without more information, nobody can say — and that is the problem.

Mathematics has met this problem at the highest level. Appel and Haken's 1976 proof of the four colour theorem rested on a computer check of 1,936 configurations that no referee could repeat by hand; doubts persisted until Robertson, Sanders, Seymour and Thomas (1997) gave a new proof with new, independently written programs, and Gonthier (2005) checked a complete proof in the Coq proof assistant. The referees of Hales's 1998 proof of the Kepler conjecture reported that they were "99% certain" but could not verify all of the computations; the Flyspeck project then verified the whole proof formally, finishing in 2014.

A computation counts as mathematical evidence only if others can repeat it and know what kind of evidence it is. This lesson develops the mathematics and the practice behind that: floating-point error bounds that tell you which digits to trust, interval arithmetic and certificates that turn computations into theorems, and the engineering — pinned environments, seeds, tests, archives — that lets someone else get your answer.

IntuitionH

What can change between two runs

A computation is a function of more inputs than the ones we think about:

Code — fixed by version control, if the exact commit is recorded.
Data — fixed by keeping raw data read-only and recording checksums.
Environment — the interpreter, libraries, compiler, linear-algebra back end, operating system and hardware. A library update may change an algorithm; a compiler may fuse a multiplication and an addition into one correctly rounded operation; different processors may use different instruction sets.
Hidden nondeterminism — the order in which parallel threads add up partial sums, unseeded random numbers, iteration order over hash tables (Python randomises string hashing between runs, so the iteration order of a set of strings can change), and notebook cells run out of order.

Each of these is either pinned (fixed and recorded) or tolerated (its effect bounded and stated). Exact computations — integers, rationals, symbolic algebra — can be pinned completely and reproduced bit for bit. Floating-point results generally cannot be guaranteed bit for bit across environments; for them the honest statement is a result with a tolerance, and the tolerance should come from an error analysis, not from hope.

There are levels of checking: re-running the same code on the same machine; reproducing it from the same code and data elsewhere; replicating with an independent implementation (or new data); and, strongest of all for proofs, formal verification.

The lesson builds on floating-point arithmetic and error analysis from numerical methods, and on the habit, from mathematical writing, of stating exactly what has and has not been established.

DefinitionsH

Definitions

Definition— reproducible

A computational result is reproducible if others can obtain the same result — exactly, or within a stated tolerance — from the same code, data and environment. It is replicable if an independent study (new code or new data) reaches consistent conclusions. (Terminology varies between communities; these are the definitions of the US National Academies' 2019 report.)

Definition— pinned environment

An environment is pinned when the exact versions of all software dependencies are fixed and recorded — typically in a lock file listing every package, version and content hash, or as a container image or a functional package specification (Nix, Guix) that rebuilds the same software.

Definition— floating-point numbers

A binary floating-point system with precision $p$ consists of $0$ and the numbers $\pm m \cdot 2^{e - p + 1}$ with integer significand $2^{p-1} \le m \lt 2^p$ (normal numbers) and exponent $e_{\min} \le e \le e_{\max}$ , together with subnormal numbers below $2^{e_{\min}}$ . IEEE 754 double precision has $p = 53$ . Rounding to nearest maps a real number to a nearest floating-point number. The unit roundoff is $u = 2^{-p}$ , so $u = 2^{-53} \approx 1.1 \times 10^{-16}$ in double precision.

Definition— evidence types

A computational claim is

exact if it is computed with exact arithmetic (integers, rationals, exact symbolic algebra) or by exhaustive enumeration; it is then as reliable as the code and its trusted base;
numerical if it is an approximation of a real quantity; it is meaningful only with an error bound, which may be heuristic, proved a priori, or rigorous a posteriori (as with interval arithmetic);
statistical if it rests on random sampling; it holds with a stated confidence and is reproducible only given the generator and seed;
experimental if it reports a pattern observed in finitely many cases, with no error model at all.

A formal proof checked by a proof assistant is stronger than all of these for a mathematical statement.

Definition— certificate

A certificate for a claim is extra output from which the claim can be verified by a check that is simpler, cheaper or independently trusted — for example Bézout coefficients for a gcd, a factor for compositeness, or a primitive root for primality.

Definition— persistent identifier

A persistent identifier names a digital object permanently: a DOI (Digital Object Identifier) for a deposited release of code or data, or an intrinsic identifier such as a Software Heritage identifier (SWHID), computed from the content itself by a cryptographic hash.

Worked examplesH

Examples

Example— Classifying claims

" $2^{127} - 1$ is prime" (Lucas, 1876; checkable today in a millisecond with integer arithmetic) — exact.
" $\zeta(3) = 1.2020569031595942\ldots$ " printed by a double-precision program — numerical; it needs a bound on the error of truncating the series and of rounding.
"The giant component of the random graph $G(1000, 0.002)$ has about 800 vertices, averaged over 100 samples" — statistical; it needs the spread of the samples and the seed.
"The pattern continues for $n \le 10^9$ " — experimental, however large the bound; the bound must be stated.

Example— Catastrophic cancellation

In double precision, $\cos(10^{-8})$ rounds to exactly $1$ , so the program (1 - cos(x))/x**2 returns $0$ at $x = 10^{-8}$ , while the true value is very close to $\tfrac12$ . The subtraction itself is exact (see Sterbenz's lemma below); it merely exposes the rounding error of $\cos x$ , which is about $u$ in absolute size and therefore enormous relative to $1 - \cos x \approx x^2/2$ . The algebraically equal expression $2\sin^2(x/2)/x^2$ involves no cancellation and is accurate to a few units of $u$ . Rearranging formulas is an error-analysis decision, not a cosmetic one.

Example— The dependency problem in interval arithmetic

Interval arithmetic computes with intervals guaranteed to contain the true values. For $X = [0, 1]$ , evaluating $x(1 - x)$ as $X \cdot (1 - X) = [0,1] \cdot [0,1] = [0, 1]$ is a valid enclosure of the true range $[0, \tfrac14]$ , but four times too wide, because the two occurrences of $x$ are treated as independent. Splitting $X$ into small subintervals gives sharper enclosures, and rewriting the expression as $\tfrac14 - (x - \tfrac12)^2$ , with a single occurrence of $x$ and an interval squaring operation (which knows that a square is non-negative), gives the exact range.

Example— A reproducible project

A layout that someone else can use without asking you anything:

project/
  README.md         what is computed, how long it takes, the one command
  CITATION.cff      how to cite this software, with its DOI
  LICENSE           terms of reuse
  environment.lock  exact versions and hashes of every dependency
  data/raw/         read-only inputs, with SHA-256 checksums
  src/              the code, under version control (git)
  tests/            unit and regression tests, with stated tolerances
  Makefile          "make reproduce" regenerates every table and figure

Every output records the commit hash, the environment and the seeds that produced it; a released version is archived with a DOI (for example on Zenodo) and is also preserved by Software Heritage.

ProofH

Error bounds, enclosures and certificates

Theorem— Theorem 1 (the standard model of rounding)

If $z$ is a real number in the normal range and $\hat z$ is a nearest floating-point number, then $\hat z = z(1 + \delta)$ with $|\delta| \le u = 2^{-p}$ . Consequently IEEE 754 arithmetic, which rounds every basic operation correctly, satisfies $\mathrm{fl}(x \circ y) = (x \circ y)(1 + \delta)$ , $|\delta| \le u$ , for $\circ \in \{+, -, \times, \div\}$ , barring overflow and underflow.

Proof

Let $2^e \le |z| \lt 2^{e+1}$ . In that range consecutive floating-point numbers are $2^{e - p + 1}$ apart, so a nearest one satisfies $|\hat z - z| \le \tfrac12 \cdot 2^{e - p + 1} = 2^{-p} \cdot 2^{e} \le u\,|z|$ . Put $\delta = (\hat z - z)/z$ . IEEE 754 requires each basic operation to return the correctly rounded exact result, so the second statement follows. $\blacksquare$

∎

Theorem— Theorem 2 (recursive summation)

Let $\hat s_0 = x_0$ and $\hat s_k = (\hat s_{k-1} + x_k)(1 + \delta_k)$ with $|\delta_k| \le u$ , and let $s_n = \sum_{i=0}^{n} x_i$ . Then

|\hat s_n - s_n| \;\le\; \big((1+u)^n - 1\big) \sum_{i=0}^{n} |x_i| \;\le\; \gamma_n \sum_{i=0}^n |x_i|, \qquad \gamma_n = \frac{nu}{1 - nu} \ \ (nu \lt 1).

Proof

Write $e_k = \hat s_k - s_k$ and $A_k = \sum_{i \le k} |x_i|$ . For $n = 0$ both sides are $0$ . For the step,

e_{k} = (\hat s_{k-1} + x_k)(1 + \delta_k) - (s_{k-1} + x_k) = e_{k-1}(1 + \delta_k) + s_k\, \delta_k .

Using the induction hypothesis $|e_{k-1}| \le ((1+u)^{k-1} - 1) A_{k-1}$ , together with $A_{k-1} \le A_k$ and $|s_k| \le A_k$ ,

|e_k| \le \big((1+u)^{k-1} - 1\big)(1 + u) A_k + u A_k = \big((1+u)^{k} - 1\big) A_k .

Finally $(1 + u)^n \le e^{nu} \le 1/(1 - nu)$ for $0 \le nu \lt 1$ (since $e^{-t} \ge 1 - t$ ), so $(1+u)^n - 1 \le nu/(1 - nu)$ . $\blacksquare$

∎

Three lessons follow. The bound involves $\sum |x_i|$ , not $|\sum x_i|$ : when there is cancellation the relative error of the sum can be huge. The bound holds for every order of summation (apply it to the reordered list), but the actual results of two orders can differ by up to twice the bound; so a parallel reduction whose order depends on thread timing is not bitwise reproducible, and a test comparing results must use a tolerance — this theorem tells you how large. And the bound grows like $nu$ , while pairwise summation achieves $O(u \log n)$ and Kahan's compensated summation $2u + O(nu^2)$ (Higham 2002, chapter 4).

Theorem— Theorem 3 (Sterbenz's lemma, 1974)

If $x$ and $y$ are floating-point numbers with $y/2 \le x \le 2y$ , then $x - y$ is a floating-point number: the subtraction is exact.

Proof

The hypothesis forces $x, y \ge 0$ , and the case $y = 0$ is trivial; if $x \lt y$ we compute $y - x$ instead (then $x \le y \le 2x$ ), so assume $0 \lt y \le x \le 2y$ . Write $y = m_y \cdot 2^{e_y - p + 1}$ , where $2^{e_y - p + 1}$ is the spacing of floating-point numbers at $y$ and $m_y \lt 2^p$ is an integer. Every floating-point number $\ge y$ has an exponent at least $e_y$ , so it is an integer multiple of $2^{e_y - p + 1}$ ; in particular $x - y = N \cdot 2^{e_y - p + 1}$ for an integer $N \ge 0$ . From $x - y \le y$ we get $N \le m_y \lt 2^p$ , and every number $N \cdot 2^{e_y - p + 1}$ with $0 \le N \lt 2^p$ and $e_y \ge e_{\min}$ is representable (as a normal or subnormal number). $\blacksquare$

∎

So cancellation never creates error; it reveals the errors already present in its operands.

Theorem— Theorem 4 (interval arithmetic is sound)

For intervals $[a, b]$ , $[c, d]$ define $[a,b] + [c,d] = [a + c, b + d]$ and $[a,b] \cdot [c,d] = [\min P, \max P]$ with $P = \{ac, ad, bc, bd\}$ . If $x \in [a,b]$ and $y \in [c, d]$ then $x + y \in [a,b] + [c,d]$ and $xy \in [a,b]\cdot[c,d]$ . Consequently, if an expression $E(x_1, \dots, x_k)$ built from $+, -, \times$ is evaluated on intervals $X_i \ni x_i$ , the resulting interval contains $E(x_1, \dots, x_k)$ .

Proof

Addition is immediate. For multiplication, suppose $y \ge 0$ ; then $a y \le x y \le b y$ . If $a \ge 0$ then $ay \ge ac$ , and if $a \le 0$ then $ay \ge ad$ ; either way $xy \ge \min P$ . If $b \ge 0$ then $by \le bd$ , and if $b \le 0$ then $by \le bc$ ; so $xy \le \max P$ . The case $y \le 0$ is symmetric, with $by \le xy \le ay$ . The statement for expressions follows by induction on the expression, using these two cases at each operation (and $X - Y = X + (-1)\cdot Y$ ). $\blacksquare$

∎

On a computer the endpoints are themselves floating-point numbers, so each computed lower endpoint is rounded down and each upper endpoint up (IEEE 754 provides these directed rounding modes). With that outward rounding Theorem 4 survives, and a computation can establish strict inequalities, enclosures of integrals and zeros of functions, and existence results via fixed-point theorems: Tucker's computer-assisted proof (2002) that the Lorenz equations have a strange attractor, answering Smale's 14th problem, is built on rigorous interval computations.

Theorem— Theorem 5 (a Bézout identity certifies a gcd)

If $d \ge 0$ , $d \mid a$ , $d \mid b$ and $d = ax + by$ for some integers $x, y$ , then $d = \gcd(a, b)$ .

Proof

Since $d$ divides $a$ and $b$ , it divides $\gcd(a, b)$ . Since $\gcd(a, b)$ divides $a$ and $b$ , it divides $ax + by = d$ . Two non-negative integers that divide each other are equal. $\blacksquare$

∎

Theorem— Theorem 6 (Lucas's primality certificate) — proof sketched

Let $p \ge 2$ . If some $a$ satisfies $a^{p-1} \equiv 1 \pmod p$ and $a^{(p-1)/q} \not\equiv 1 \pmod p$ for every prime $q \mid p - 1$ , then $p$ is prime.

Proof

The conditions say that $a$ is a unit modulo $p$ of multiplicative order exactly $p - 1$ . So the unit group of $\Z/p\Z$ has at least $p - 1$ elements, i.e. every non-zero residue is a unit, which happens only when $p$ is prime. $\blacksquare$

∎

Pratt (1975) observed that the prime factors $q$ can be certified recursively in the same way, giving certificates of polynomial size: primality is in NP (and, by Agrawal, Kayal and Saxena, 2002, in P). The point for reproducibility is that checking a certificate is far easier than trusting the program that found it.

For statistical claims, the analogue of an error bound is a probability bound. If $X_1, \dots, X_n$ are independent with values in $[\alpha, \beta]$ and mean $\mu$ , Hoeffding's inequality (1963) gives $\Pr\big(|\bar X - \mu| \ge \varepsilon\big) \le 2 \exp\!\big(-2 n \varepsilon^2 / (\beta - \alpha)^2\big)$ for every $\varepsilon \gt 0$ , a non-asymptotic guarantee for Monte Carlo estimates; the familiar interval $\bar X \pm 1.96\, \hat\sigma/\sqrt n$ is the asymptotic version from the central limit theorem. Either is a statement about the method: about 5% of honest 95% intervals miss the truth.

ExplanationH

Packaging and reporting a computation

Version control. Keep code (and small data) in a version-control system such as git, and record in every output the commit that produced it. A result that cannot be traced to a commit cannot be reproduced.

Pinned environment and one command. Record exact versions with hashes in a lock file, or build a container or a Nix/Guix environment. Provide a single command that regenerates every number, table and figure from the raw data. Run it on a clean machine before publishing.

Tests. Unit tests check the pieces on cases with known answers; regression tests compare the outputs with recorded values — exactly for exact computations, and within a tolerance justified by an error analysis such as Theorem 2 for floating point. A test with tolerance 1e-6 chosen because "it passed" is not evidence.

Notebooks. A notebook that "ran once" may depend on cells run out of order, on variables that no longer exist in the code, and on an unrecorded environment. Restart the kernel and run all cells from the top, and pin its environment like any other code.

Randomness. Seed every generator and record the seeds, the generator and its version. Python guarantees that random.random() gives the same stream for the same seed across versions, but other methods may change; NumPy's Generator with explicitly spawned independent streams is the safe pattern for parallel work. Report the spread over seeds, not the luckiest run.

Data provenance. Keep raw data immutable, record where it came from and checksums of the files, and script every transformation.

Archiving and citation. Deposit each released version in an archive that issues a DOI (for example Zenodo); Software Heritage additionally preserves public source code with intrinsic identifiers. Cite software like a paper: authors, title, version, year and identifier — a CITATION.cff file makes this easy — and cite the libraries your result depends on. zbMATH's swMATH database indexes mathematical software and the papers that use it.

Reporting. State the evidence type of each computational claim, its bounds (searched ranges, error bounds, confidence levels) and what it does not establish. This site follows the same rule: every experiment declares its evidence types and ends with an "Epistemic status" paragraph, and its "formally verified" and "computationally tested" labels are derived from build reports, never written by hand.

Formal mathematicsF

Error bounds and certificates in Lean

The Lean module proves Theorem 2 in the standard model. Recursive summation with the rounding errors as data:

fsumLean 4Formally verified

/-- Recursive (left-to-right) summation in floating point, in the *standard model*: each
addition is exact up to a relative error, `fl(a + b) = (a + b)(1 + δ)`. The sequence `δ`
records the rounding error committed at each step; for IEEE double precision `|δ k| ≤ 2⁻⁵³`. -/
def fsum (x δ : ℕ → ℝ) : ℕ → ℝ
  | 0 => x 0
  | n + 1 => (fsum x δ n + x (n + 1)) * (1 + δ (n + 1))

Axioms used: propext, Classical.choice, Quot.soundformal/Lumenary/ComputationalReproducibility.leanTry it in the Lean web editor

The bound, by induction on $n$ with exactly the steps of the proof above:

fsum_errorLean 4Proved here from basic factsFormally verified

/-- **Error bound for recursive summation.** If every rounding error satisfies `|δ k| ≤ u`, the
computed sum of `x 0, …, x n` differs from the exact sum by at most
`((1 + u)ⁿ - 1) · Σ |xᵢ|`. (For `n u < 1` this is at most `γₙ Σ |xᵢ|` with `γₙ = nu/(1 - nu)`,
the form in Higham's book.) Induction on `n`: the new error is the old one, amplified by
`1 + δ`, plus the rounding error `(sₙ + xₙ₊₁) δ` of the new addition. -/
theorem fsum_error (x δ : ℕ → ℝ) (u : ℝ) (hu : 0 ≤ u) (hδ : ∀ k, |δ k| ≤ u) (n : ℕ) :
    |fsum x δ n - ∑ i ∈ range (n + 1), x i| ≤
      ((1 + u) ^ n - 1) * ∑ i ∈ range (n + 1), |x i| := by
  induction n with
  | zero => simp [fsum]
  | succ n ih =>
    rw [sum_range_succ (fun i => x i) (n + 1), sum_range_succ (fun i => |x i|) (n + 1)]
    set s := ∑ i ∈ range (n + 1), x i with hs
    set A := ∑ i ∈ range (n + 1), |x i| with hA
    set e := fsum x δ n - s with he
    have hA0 : 0 ≤ A := sum_nonneg fun i _ => abs_nonneg (x i)
    have hsA : |s| ≤ A := abs_sum_le_sum_abs _ _
    have hpow : 1 ≤ (1 + u) ^ n := one_le_pow₀ (by linarith)
    have hd : |δ (n + 1)| ≤ u := hδ (n + 1)
    have h1d : |1 + δ (n + 1)| ≤ 1 + u := by
      calc |1 + δ (n + 1)| ≤ |(1 : ℝ)| + |δ (n + 1)| := abs_add_le _ _
        _ ≤ 1 + u := by rw [abs_one]; linarith
    have hsx : |s + x (n + 1)| ≤ A + |x (n + 1)| := by
      calc |s + x (n + 1)| ≤ |s| + |x (n + 1)| := abs_add_le _ _
        _ ≤ A + |x (n + 1)| := by linarith
    have key : fsum x δ (n + 1) - (s + x (n + 1)) =
        e * (1 + δ (n + 1)) + (s + x (n + 1)) * δ (n + 1) := by
      simp only [fsum, he]; ring
    rw [key]
    calc |e * (1 + δ (n + 1)) + (s + x (n + 1)) * δ (n + 1)|
        ≤ |e| * |1 + δ (n + 1)| + |s + x (n + 1)| * |δ (n + 1)| := by
          rw [← abs_mul, ← abs_mul]; exact abs_add_le _ _
      _ ≤ (((1 + u) ^ n - 1) * A) * (1 + u) + (A + |x (n + 1)|) * u :=
          add_le_add
            (mul_le_mul ih h1d (abs_nonneg _) (mul_nonneg (sub_nonneg.2 hpow) hA0))
            (mul_le_mul hsx hd (abs_nonneg _) (by positivity))
      _ ≤ ((1 + u) ^ (n + 1) - 1) * (A + |x (n + 1)|) := by
          have hx := abs_nonneg (x (n + 1))
          have hp : (1 + u) ^ (n + 1) = (1 + u) ^ n * (1 + u) := pow_succ _ _
          rw [hp]
          nlinarith [mul_nonneg hx (sub_nonneg.2 hpow), mul_nonneg hu hx,
            mul_nonneg (mul_nonneg hx (sub_nonneg.2 hpow)) hu]

Axioms used: propext, Classical.choice, Quot.soundformal/Lumenary/ComputationalReproducibility.leanTry it in the Lean web editor

Interval arithmetic with rational endpoints (exact, so no outward rounding is needed) and Theorem 4 for $+$ and $\times$ :

IvalLean 4Formally verified

/-- A closed interval with rational endpoints (in practice: floating-point numbers, which are
rationals). A real number `x` lies in `I` when `I.lo ≤ x ≤ I.hi`. -/
structure Ival where
  /-- The lower endpoint. -/
  lo : ℚ
  /-- The upper endpoint. -/
  hi : ℚ

Axioms used: No axiomsformal/Lumenary/ComputationalReproducibility.leanTry it in the Lean web editor

MemLean 4Formally verified

/-- Membership of a real number in an interval. -/
def Ival.Mem (I : Ival) (x : ℝ) : Prop := (I.lo : ℝ) ≤ x ∧ x ≤ (I.hi : ℝ)

Axioms used: propext, Classical.choice, Quot.soundformal/Lumenary/ComputationalReproducibility.leanTry it in the Lean web editor

addLean 4Formally verified

/-- Interval addition: add the endpoints. -/
def Ival.add (I J : Ival) : Ival := ⟨I.lo + J.lo, I.hi + J.hi⟩

Axioms used: propext, Classical.choice, Quot.soundformal/Lumenary/ComputationalReproducibility.leanTry it in the Lean web editor

mulLean 4Formally verified

/-- Interval multiplication: the smallest and largest of the four endpoint products. -/
def Ival.mul (I J : Ival) : Ival :=
  ⟨min (min (I.lo * J.lo) (I.lo * J.hi)) (min (I.hi * J.lo) (I.hi * J.hi)),
   max (max (I.lo * J.lo) (I.lo * J.hi)) (max (I.hi * J.lo) (I.hi * J.hi))⟩

Axioms used: propext, Classical.choice, Quot.soundformal/Lumenary/ComputationalReproducibility.leanTry it in the Lean web editor

add_memLean 4Proved here from basic factsFormally verified

/-- **Soundness of interval addition**: if `x ∈ I` and `y ∈ J` then `x + y ∈ I + J`. -/
theorem Ival.add_mem {I J : Ival} {x y : ℝ} (hx : I.Mem x) (hy : J.Mem y) :
    (I.add J).Mem (x + y) := by
  obtain ⟨hx1, hx2⟩ := hx
  obtain ⟨hy1, hy2⟩ := hy
  simp only [Ival.Mem, Ival.add, Rat.cast_add]
  constructor <;> linarith

Axioms used: propext, Classical.choice, Quot.soundformal/Lumenary/ComputationalReproducibility.leanTry it in the Lean web editor

mul_memLean 4Proved here from basic factsFormally verified

/-- **Soundness of interval multiplication** (the "fundamental theorem of interval arithmetic"
for `×`): if `x ∈ [a, b]` and `y ∈ [c, d]` then `x y` lies between the least and the greatest
of `ac, ad, bc, bd`. For the lower bound: if `y ≥ 0` then `xy ≥ ay`, and `ay ≥ ac` or `ay ≥ ad`
according to the sign of `a`; if `y ≤ 0` argue with `b` instead. The upper bound is symmetric. -/
theorem Ival.mul_mem {I J : Ival} {x y : ℝ} (hx : I.Mem x) (hy : J.Mem y) :
    (I.mul J).Mem (x * y) := by
  obtain ⟨hx1, hx2⟩ := hx
  obtain ⟨hy1, hy2⟩ := hy
  simp only [Ival.Mem, Ival.mul, Rat.cast_min, Rat.cast_max, Rat.cast_mul]
  set a : ℝ := (I.lo : ℝ)
  set b : ℝ := (I.hi : ℝ)
  set c : ℝ := (J.lo : ℝ)
  set d : ℝ := (J.hi : ℝ)
  constructor
  · rcases le_total 0 y with hy0 | hy0
    · rcases le_total 0 a with ha | ha
      · exact min_le_of_left_le (min_le_of_left_le (by nlinarith))
      · exact min_le_of_left_le (min_le_of_right_le (by nlinarith))
    · rcases le_total 0 b with hb | hb
      · exact min_le_of_right_le (min_le_of_left_le (by nlinarith))
      · exact min_le_of_right_le (min_le_of_right_le (by nlinarith))
  · rcases le_total 0 y with hy0 | hy0
    · rcases le_total 0 b with hb | hb
      · exact le_max_of_le_right (le_max_of_le_right (by nlinarith))
      · exact le_max_of_le_right (le_max_of_le_left (by nlinarith))
    · rcases le_total 0 a with ha | ha
      · exact le_max_of_le_left (le_max_of_le_right (by nlinarith))
      · exact le_max_of_le_left (le_max_of_le_left (by nlinarith))

Axioms used: propext, Classical.choice, Quot.soundformal/Lumenary/ComputationalReproducibility.leanTry it in the Lean web editor

A numerical value as a theorem — the same kind of statement as Zu Chongzhi's bounds on $\pi$ (see the history below):

sqrt_two_boundsLean 4Proved here from basic factsFormally verified

/-- The numerical claim "`√2 = 1.41421356…`" as a theorem with explicit bounds. The proof is
exact rational arithmetic: `1.41421356² < 2 < 1.41421357²`. -/
theorem sqrt_two_bounds : (1.41421356 : ℝ) < √2 ∧ √2 < 1.41421357 := by
  constructor
  · rw [Real.lt_sqrt (by norm_num)]; norm_num
  · rw [Real.sqrt_lt' (by norm_num)]; norm_num

Axioms used: propext, Classical.choice, Quot.soundformal/Lumenary/ComputationalReproducibility.leanTry it in the Lean web editor

Certificates. Theorem 5 is proved here; Theorem 6 is Mathlib's lucas_primality, and is then used to certify that $65537$ is prime from the single witness $a = 3$ :

gcd_certificateLean 4Proved here from basic factsFormally verified

/-- **A Bézout identity certifies a gcd.** If `d ≥ 0` divides `a` and `b`, and `d = a x + b y`
for some integers `x, y`, then `d = gcd(a, b)`. A program that outputs `(d, x, y)` can thus be
checked with two divisions and two multiplications, without trusting the program. -/
theorem gcd_certificate {a b d x y : ℤ} (hd : 0 ≤ d) (hda : d ∣ a) (hdb : d ∣ b)
    (hbez : a * x + b * y = d) : d = Int.gcd a b := by
  have h1 : d ∣ (Int.gcd a b : ℤ) := Int.dvd_coe_gcd hda hdb
  have h2 : (Int.gcd a b : ℤ) ∣ d := by
    rw [← hbez]
    exact dvd_add (dvd_mul_of_dvd_left (Int.gcd_dvd_left a b) _)
      (dvd_mul_of_dvd_left (Int.gcd_dvd_right a b) _)
  exact Int.dvd_antisymm hd (Int.natCast_nonneg _) h1 h2

Axioms used: propext, Quot.soundformal/Lumenary/ComputationalReproducibility.leanTry it in the Lean web editor

lucas_certificateLean 4Applies a library theoremFormally verified

/-- **Lucas's primality certificate.** If `a^(p-1) = 1` in `ZMod p` but `a^((p-1)/q) ≠ 1` for
every prime `q ∣ p - 1`, then `p` is prime. Checking the certificate `(a, factorisation of p - 1)`
needs only modular exponentiations. This applies Mathlib's `lucas_primality`. -/
theorem lucas_certificate (p : ℕ) (a : ZMod p) (ha : a ^ (p - 1) = 1)
    (hd : ∀ q : ℕ, q.Prime → q ∣ p - 1 → a ^ ((p - 1) / q) ≠ 1) : p.Prime :=
  lucas_primality p a ha hd

Axioms used: propext, Classical.choice, Quot.soundformal/Lumenary/ComputationalReproducibility.leanTry it in the Lean web editor

prime_65537Lean 4Proved here from basic factsFormally verified

/-- `65537 = 2¹⁶ + 1` is prime, certified by `a = 3`: since `p - 1 = 2¹⁶`, the only prime
`q` is `2`, and the certificate consists of `3^(2¹⁶) ≡ 1` and `3^(2¹⁵) ≡ -1 (mod 65537)`.
The two exponentiations are evaluated exactly by the kernel's arithmetic on natural numbers. -/
theorem prime_65537 : Nat.Prime 65537 := by
  have hcast : ∀ n : ℕ, (3 : ZMod 65537) ^ n = ((3 ^ n % 65537 : ℕ) : ZMod 65537) := by
    intro n; rw [ZMod.natCast_mod]; push_cast; rfl
  apply lucas_certificate 65537 3
  · have h : 3 ^ 65536 % 65537 = 1 := by decide +kernel
    rw [hcast, h]; rfl
  · intro q hq hdvd
    have h2 : q = 2 := by
      have : q ∣ 2 ^ 16 := by simpa using hdvd
      exact (Nat.prime_dvd_prime_iff_eq hq Nat.prime_two).1 (hq.dvd_of_dvd_pow this)
    subst h2
    have h : 3 ^ 32768 % 65537 = 65536 := by decide +kernel
    rw [hcast, h]; decide

Axioms used: propext, Classical.choice, Quot.soundformal/Lumenary/ComputationalReproducibility.leanTry it in the Lean web editor

The last proof asks Lean's kernel to evaluate $3^{65536} \bmod 65537$ (decide +kernel); the kernel does this with built-in arbitrary-precision arithmetic on natural numbers, which is therefore part of the trusted base of that proof — a small, well-tested part, but a part.

What is formalised: Theorems 2, 4 (for $+$ and $\times$ ), 5 and the certified examples; Theorem 6 is from Mathlib. Not formalised: IEEE 754 itself (Theorem 1 and Sterbenz's lemma need a model of floating-point formats; the Flocq library for Coq provides one), outward rounding, and Hoeffding's inequality.

ComputationC

Experiments on reproducibility

The first experiment sums the same 100,000 numbers in several orders. Because every double is a rational number, the exact sum can be computed with integers, so every error is measured exactly. Watch how many different "answers" appear, how small they are compared with the worst-case bound of Theorem 2, and how pairwise and compensated summation do.

Summing the same numbers in different ordersPythonExact computationNumerical approximationComputationally tested

"""Floating-point summation: why the order of operations changes the answer.

Every double-precision float is a rational number with a power-of-two
denominator, so the exact sum of a list of floats can be computed with
integers. Against that exact value the experiment measures the error of
summing the same 100,000 numbers in different orders (as a parallel reduction
might), by pairwise summation, by Kahan's compensated summation and by
math.fsum, and compares the errors with the proven bound
((1 + u)^(n-1) - 1) * sum |x_i|, u = 2^-53, from the lesson. The exact sums are
exact; the floating-point results are what this machine's IEEE 754 arithmetic
produced, which other hardware or libraries need not reproduce bit for bit.
"""

import math
import random
from fractions import Fraction

U = 2.0 ** -53          # unit roundoff of IEEE double precision

def exact_sum(xs):
    """The exact sum of a list of floats, as a Fraction (via integers)."""
    ratios = [x.as_integer_ratio() for x in xs]
    big = max((d for _, d in ratios), default=1)   # all are powers of two
    return Fraction(sum(n * (big // d) for n, d in ratios), big)

def recursive_sum(xs):
    """Left-to-right summation, one rounding per addition."""
    s = 0.0
    for x in xs:
        s += x
    return s

def pairwise_sum(xs):
    """Pairwise (tree) summation: error grows like log n instead of n."""
    if len(xs) <= 8:
        return recursive_sum(xs)
    mid = len(xs) // 2
    return pairwise_sum(xs[:mid]) + pairwise_sum(xs[mid:])

def kahan_sum(xs):
    """Kahan's compensated summation: carries the lost low-order bits."""
    s, c = 0.0, 0.0
    for x in xs:
        y = x - c
        t = s + y
        c = (t - s) - y
        s = t
    return s

def recursive_bound(xs):
    """The lesson's bound ((1 + u)^(n-1) - 1) * sum |x_i|, rounded up a little."""
    n = len(xs)
    factor = math.expm1((n - 1) * math.log1p(U))
    return factor * math.fsum(abs(x) for x in xs) * (1 + 1e-12)

def make_data(n, seed=2024):
    """n seeded numbers of mixed sign and magnitude (1e-6 to 1e6)."""
    rng = random.Random(seed)
    return [rng.uniform(-1.0, 1.0) * 10.0 ** rng.randint(-6, 6)
            for _ in range(n)]

def main():
    a, b, c = 0.1, 0.2, 0.3
    print(f"(0.1 + 0.2) + 0.3 = {(a + b) + c!r}")
    print(f"0.1 + (0.2 + 0.3) = {a + (b + c)!r}")
    print(f"the float 0.1 is exactly {Fraction(a)}")
    print()
    xs = make_data(100_000)
    exact = exact_sum(xs)
    rng = random.Random(7)
    shuffled = xs[:]
    rng.shuffle(shuffled)
    methods = [
        ("given order", recursive_sum(xs)),
        ("shuffled order", recursive_sum(shuffled)),
        ("increasing |x|", recursive_sum(sorted(xs, key=abs))),
        ("decreasing |x|", recursive_sum(sorted(xs, key=abs, reverse=True))),
        ("pairwise", pairwise_sum(xs)),
        ("Kahan", kahan_sum(xs)),
        ("math.fsum", math.fsum(xs)),
    ]
    print(f"n = {len(xs)}; exact sum = {float(exact)!r} (rounded to print)")
    print(f"{'method':<16}{'computed sum':>26}{'|error|':>12}")
    for name, value in methods:
        err = abs(Fraction(value) - exact)
        print(f"{name:<16}{value!r:>26}{float(err):>12.2e}")
    distinct = len({value for _, value in methods})
    print(f"{distinct} different results from {len(methods)} methods.")
    print(f"Proven bound for any recursive order: {recursive_bound(xs):.2e}")
    print()
    print("Epistemic status: the exact sum and every error are exact rational")
    print("arithmetic; the computed sums are what this machine produced. The")
    print("observed errors lie far below the worst-case bound, which is a theorem")
    print("(proved in Lean) under the standard model of rounding; a different")
    print("order, compiler or library may legitimately give different last digits.")

if __name__ == "__main__":
    main()

Runs entirely in your browser, in an isolated worker. Nothing is sent to a server.

The second experiment fixes a seed and checks, with a 64-bit fingerprint of the output, that two runs agree bit for bit; then it shows what different seeds do, and how often the nominal 95% confidence intervals actually contain the exact answer.

Seeds, fingerprints and confidence intervals in Monte CarloPythonExact computationStatistical evidenceComputationally tested

"""Seeds, pseudorandom streams and honest statistical claims.

A pseudorandom generator is a deterministic function of its seed, so a
Monte Carlo computation is reproducible bit for bit once the generator, its
version and the seed are fixed - which the experiment checks with a 64-bit
fingerprint (FNV-1a) of the output. Different seeds give different estimates;
the honest report is an estimate with a standard error. Over 200 seeds the
experiment counts how often the nominal 95% interval contains the true value
of the integral of exp(-x^2) over [0, 1], known exactly through erf. That
coverage count is statistical evidence about the method, not a proof.
"""

import math
import random
import statistics

TRUE_VALUE = math.sqrt(math.pi) / 2 * math.erf(1.0)

def fnv1a64(text):
    """64-bit FNV-1a hash of a string: a cheap fingerprint of an output."""
    h = 0xCBF29CE484222325
    for byte in text.encode("utf-8"):
        h ^= byte
        h = (h * 0x100000001B3) % 2 ** 64
    return h

def monte_carlo(seed, n):
    """Estimate the integral of exp(-x^2) on [0, 1]; return (mean, std. error)."""
    rng = random.Random(seed)
    values = [math.exp(-rng.random() ** 2) for _ in range(n)]
    return statistics.fmean(values), statistics.stdev(values) / math.sqrt(n)

def coverage(seeds, n, z=1.96):
    """How many of the intervals mean +- z*SE (one per seed) contain the truth."""
    hits = 0
    for seed in seeds:
        mean, se = monte_carlo(seed, n)
        hits += abs(mean - TRUE_VALUE) <= z * se
    return hits

def minstd(seed, k):
    """First k outputs of the Park-Miller 'minimal standard' generator."""
    out, x = [], seed
    for _ in range(k):
        x = 16807 * x % (2 ** 31 - 1)
        out.append(x)
    return out

def main():
    print(f"Target: integral of exp(-x^2) on [0, 1] = {TRUE_VALUE:.15f}")
    print("Park-Miller generator from seed 42: "
          + ", ".join(map(str, minstd(42, 3))) + ", ...")
    print()
    run_a, run_b = monte_carlo(2024, 5000), monte_carlo(2024, 5000)
    fa, fb = fnv1a64(repr(run_a)), fnv1a64(repr(run_b))
    print(f"seed 2024, run A: {run_a[0]!r}  fingerprint {fa:016x}")
    print(f"seed 2024, run B: {run_b[0]!r}  fingerprint {fb:016x}")
    print(f"bit-for-bit identical: {run_a == run_b}")
    covers = abs(run_a[0] - TRUE_VALUE) <= 1.96 * run_a[1]
    print(f"its 95% interval +- {1.96 * run_a[1]:.6f} covers the truth: {covers}")
    print()
    print("Same method, five other seeds (estimate +- 1.96 standard errors):")
    for seed in range(1, 6):
        mean, se = monte_carlo(seed, 5000)
        print(f"  seed {seed}: {mean:.6f} +- {1.96 * se:.6f}")
    hits = coverage(range(100, 300), 1000)
    print(f"Nominal 95% intervals (n = 1000) covering the truth: {hits}/200")
    print()
    print("Epistemic status: the identity of the two seeded runs is exact (same")
    print("generator, same seed, same Python); the estimates are statistical, and")
    print("the coverage count is statistical evidence that the intervals behave")
    print("as advertised for this integrand. Python guarantees the random()")
    print("stream for a given seed, but other methods may change between versions.")

if __name__ == "__main__":
    main()

Runs entirely in your browser, in an isolated worker. Nothing is sent to a server.

Explore

Explore: cancellation and seeds

Both curves below are the same function in exact arithmetic, $\frac{1 - \cos x}{x^2} = \frac{2\sin^2(x/2)}{x^2}$ , which tends to $\tfrac12$ as $x \to 0$ . The window is $|x| \le 10^{-7}$ . The first formula, evaluated in double precision, dissolves into steps and noise as $x$ shrinks; the second stays at $\tfrac12$ . Use the trace to read off values near $x = 2 \times 10^{-8}$ .

Catastrophic cancellation, live

Move the sliders and trace along the curves to read off values.

(1 − cos x)/x², as written
2 sin²(x/2)/x², rearranged

Trace along xx = 1×10⁻⁷

(1 − cos x)/x², as written: y = 0.4996
2 sin²(x/2)/x², rearranged: y = 0.5

The exact pattern of the noise may differ between browsers: JavaScript's Math.cos is not required to be correctly rounded. That is itself a reproducibility lesson.

The simulation below estimates $\pi$ with Buffon's needle. Replay reruns the same seed and gives exactly the same result; New run changes the seed. Compare several runs with the reported margin: the estimates differ, but they should scatter within it about 95% of the time.

Same seed, same result

Drop needles on a floor ruled with lines one needle-length apart. The chance of crossing a line is 2/π, so needles can estimate π.

crosses a line
misses the lines

Number of needles1,000

623 of 1,000 needles cross a line

Simulated0.623

Exact2/π ≈ 0.637

Estimate of π = 2 × needles ÷ crossings (± 95% margin)

Simulated3.21 ± 0.15

Exactπ ≈ 3.14159

Actual error this time: 0.069. The margin shrinks only like 1/√N, so one more correct digit needs about 100 times as many needles.

Evidence, not proofThis is a simulation of 1,000 needles (run 1) using a seeded random number generator. It is evidence, not proof: the exact values come from mathematics, and a new run gives slightly different numbers.

The needle length equals the line spacing. The simulation chooses needle directions without using π (it picks random points in a disc), so the estimate does not sneak π in.

Exercises

1. Computational reproducibility

Warm-upHuman

Which of these reports state their evidence honestly? (Choose all that apply.)

Choose all that apply

"With seed 7, NumPy 1.26.4 and

10^6

samples, our Monte Carlo estimate is

0.7466 \pm 0.0011

(95% confidence interval).""We checked the conjecture for all

n \le 10^9

; it is therefore true.""Interval arithmetic with outward rounding shows

1.41421356 \lt \sqrt 2 \lt 1.41421357

.""Our formula agrees with the computed value to 15 digits, which proves it.""Summed in double precision; the result is exact because floating point is deterministic."

2. Computational reproducibility

Warm-upComputational

In IEEE double precision ( $p = 53$ ), let x = 1e16. What does (x + 1) - x evaluate to?

3. Computational reproducibility

CoreHuman

In double precision, 1 - math.cos(1e-8) evaluates to 0.0, although $1 - \cos(10^{-8}) \approx 5 \times 10^{-17}$ . Which explanation is correct?

Choose one

math.cos(1e-8) is correctly rounded to 1.0, because doubles next to

1

are about

10^{-16}

apart; the subtraction is then exact and reveals that earlier rounding.Subtracting two nearly equal numbers is itself an inaccurate floating-point operation.The library's cosine function has a bug.No floating-point formula can compute this quantity accurately.

4. Computational reproducibility

CoreHuman

Choose one

Because of cancellation: the exact sum can be

0

while the partial sums, and hence the rounding errors, are large.Because absolute values make the proof easier, but the stronger bound is also true.Because the bound must also hold for complex numbers.Because the rounding errors

\delta_k

might all have the same sign.

5. Computational reproducibility

CoreHuman

Using interval multiplication, compute $[-1, 2] \cdot [-1, 2]$ . What is the lower endpoint of the result?

6. Computational reproducibility

CoreHuman

A program you do not trust outputs: $\gcd(1071, 462) = 21$ , with certificate $x = -3$ , $y = 7$ . What can you conclude?

Choose one

The claim is proved:

1071\cdot(-3) + 462\cdot 7 = 21

, and

21

divides both

1071 = 51 \cdot 21

and

462 = 22 \cdot 21

.Nothing, unless we rerun Euclid's algorithm ourselves.The certificate is invalid, because Bézout coefficients must be positive.The certificate only shows that

21

is a common divisor.

7. Computational reproducibility

CoreFormal

In the Lean module, fsum x δ n is the floating-point recursive sum of x 0, …, x n with rounding errors δ. Which statement correctly formalises Theorem 2?

Choose one

∀ u, 0 ≤ u → (∀ k, |δ k| ≤ u) → |fsum x δ n - ∑ i ∈ range (n+1), x i| ≤ ((1+u)^n - 1) * ∑ i ∈ range (n+1), |x i|

|fsum x δ n - ∑ i ∈ range (n+1), x i| ≤ ((1+u)^n - 1) * ∑ i ∈ range (n+1), |x i| with no hypothesis on δ

∀ u, 0 ≤ u → (∀ k, |δ k| ≤ u) → |fsum x δ n - ∑ i ∈ range (n+1), x i| ≤ u * |∑ i ∈ range (n+1), x i|

∃ δ, fsum x δ n = ∑ i ∈ range (n+1), x i

8. Computational reproducibility

Warm-upHuman

You want a colleague to be able to regenerate every figure of your paper in five years' time. Which is the best plan?

Choose one

A versioned repository with a lock file of exact dependency versions, recorded seeds, tests, one command that rebuilds everything, and a release archived with a DOI.A requirements.txt saying numpy>=1.20 and the notebook as it was last run.High-resolution images of all the figures in the supplementary material.A detailed description of the algorithms in the paper; the code is not needed.

9. Computational reproducibility

StretchComputational

Implement Kahan's compensated summation kahan_sum(xs). Keep a running sum s and a compensation c for the low-order bits lost so far: for each x, set y = x - c, t = s + y, c = (t - s) - y, s = t. The tests compare your result with the exact rational sum.

Runs entirely in your browser, in an isolated worker. Nothing is sent to a server.

Common misconceptionsH

Common misconceptions

Misconception

"The notebook ran, so it is reproducible." A notebook that ran once on one machine may depend on out-of-order cells, on unrecorded library versions and on unseeded randomness. Reproducibility is demonstrated by someone else, from a clean environment, with one command.

Misconception

"Agreement to many digits is a proof." The number $e^{\pi\sqrt{163}} = 262537412640768743.99999999999925\ldots$ agrees with an integer to twelve decimal places, and is not an integer. Agreement is evidence that guides the search for a proof — or a certified computation — and is reported as numerical evidence.

Misconception

"Floating point is random noise, so its results mean nothing." Floating-point arithmetic is deterministic and precisely specified, and its errors can be bounded by theorems such as Theorem 2. Results are meaningful when reported with a justified tolerance.

Misconception

"Exact arithmetic means the result is proved." Exact integer or symbolic computation removes rounding error, not programming error. An exact computation is as reliable as the code, the libraries and the hardware; a certificate, an independent reimplementation or a formal proof is needed to go further.

Misconception

"Setting a seed makes a statistical result correct." A seed makes the run repeatable. Whether the conclusion holds is a statistical question about the spread over seeds; reporting only a favourable seed is a form of selective reporting.

ViewpointsH

Many ways of seeing

Analysis— Backward error and conditioning

Wilkinson's backward error analysis asks not "how wrong is the answer?" but "for which nearby input is the computed answer exact?". Theorem 2's proof shows that recursive summation computes the exact sum of perturbed data $x_i(1 + \theta_i)$ with $|\theta_i| \le \gamma_n$ : it is backward stable. The forward error then obeys $\text{forward error} \lesssim \kappa \cdot \text{backward error}$ , where the condition number $\kappa$ measures the sensitivity of the problem itself — for summation $\kappa = \sum|x_i| / |\sum x_i|$ , which is large precisely when there is cancellation. Good algorithms cannot fix ill-conditioned problems; they can only avoid adding instability. See numerical analysis.

Geometry— Rounding as a projection; ill-conditioning as distance

Floating-point numbers form a grid on the real line whose spacing doubles at each power of two, and rounding to nearest is the nearest-point projection onto that grid: relative distances, not absolute ones, are preserved, which is why the standard model is multiplicative. Conditioning has a geometric meaning too. For an invertible matrix $A$ , the relative distance to the nearest singular matrix in the spectral norm is exactly $1/\kappa_2(A) = \sigma_{\min}(A)/\sigma_{\max}(A)$ (Eckart–Young), and Demmel (1987) showed that for many problems the condition number is the reciprocal of the distance to the set of ill-posed problems. A problem is hard to compute reproducibly when it lies close to a singular one.

Algebra— Floating-point numbers are not a field

Floating-point addition is commutative but not associative ( $(0.1 + 0.2) + 0.3 \ne 0.1 + (0.2 + 0.3)$ ), and multiplication does not distribute over addition, so floats form no ring — algebraic rewritings that are identities over $\R$ change the result, which is why compilers may not reorder floating-point arithmetic without permission. Exact computation lives in genuine algebraic structures ( $\Z$ , $\Q$ , finite fields, polynomial rings), and certificates are often algebraic identities: Bézout's identity certifies a gcd, a Nullstellensatz identity $\sum g_i f_i = 1$ certifies that polynomial equations have no common complex solution, and a sum-of-squares decomposition certifies that a polynomial is non-negative. Each can be checked by exact multiplication and addition. See rings.

Probability— Statistical reproducibility

A Monte Carlo estimate is a random variable; reproducing it bit for bit only reproduces one sample. The meaningful claim is about its distribution: a confidence interval from Hoeffding's inequality or the central limit theorem, whose coverage can itself be checked by repeating the experiment over many seeds (as in the second experiment). Pseudorandom generators are deterministic, and their statistical quality is tested empirically (for example with L'Ecuyer and Simard's TestU01 suite); for rigorous statements one uses either certified randomness assumptions or bounds that hold for every outcome. See statistical inference.

Logic— Proof checking as the strongest reproducibility

A formal proof is a finite object that a small program — the kernel of a proof assistant — checks deterministically; anyone can re-check it, and the trusted base is explicit. The de Bruijn criterion asks that this checker be small enough to be trusted by inspection. Large computations can be brought into this framework by certificates: SAT solvers emit DRAT proofs of unsatisfiability that verified checkers validate, which is how the 200-terabyte proof of the Boolean Pythagorean triples theorem (Heule, Kullmann and Marek, 2016) was checked. See formalising published mathematics.

History & global context

History and global context

Reporting precision. Zu Chongzhi in fifth-century China bounded $\pi$ between $3.1415926$ and $3.1415927$ — a numerical result stated as a rigorous enclosure. In Kerala, Mādhava (c. 1400) gave series for $\pi$ with correction terms that control the truncation error, recorded in Jyeṣṭhadeva's Yuktibhāṣā. In Samarkand, Jamshīd al-Kāshī's Treatise on the Circumference (1424) computed $2\pi$ to the equivalent of sixteen decimal places, choosing a polygon of $3 \cdot 2^{28}$ sides so that the error was guaranteed to be small enough.
Rounding error as mathematics. Turing (1948) and von Neumann and Goldstine (1947) analysed rounding errors in matrix computations; Wilkinson developed backward error analysis (Rounding Errors in Algebraic Processes, 1963). Kahan led the design of the IEEE 754 standard (1985, revised 2008 and 2019), which made floating-point results specified and portable, and received the 1989 Turing Award.
Interval arithmetic was proposed by Warmus in Poland (1956) and Sunaga in Japan (1958), and developed by Moore (Interval Analysis, 1966); validated numerics later underpinned computer-assisted proofs such as Tucker's (2002).
Reproducible research. The geophysicist Jon Claerbout made his Stanford group's papers reproducible from their code around 1990; Buckheit and Donoho (1995) summarised his view: an article about computational science is "merely advertising of the scholarship", the scholarship being the code and data that produced its figures. Knuth's literate programming (1984) and the frozen, versioned design of TeX are early examples of software built to give the same output for decades.
Infrastructure. CERN launched Zenodo in 2013; Software Heritage was started by Inria in 2016 with UNESCO support; the FAIR principles (findable, accessible, interoperable, reusable) were published in 2016; and the US National Academies' report Reproducibility and Replicability in Science appeared in 2019.

A result someone else can check

What can change between two runs

Definitions

Examples

Error bounds, enclosures and certificates

Packaging and reporting a computation

Error bounds and certificates in Lean

Experiments on reproducibility

Explore: cancellation and seeds

Exercises

Common misconceptions

Many ways of seeing

History and global context

Further reading