\documentclass[runningheads]{llncs}
\usepackage[T1]{fontenc}
\usepackage{lmodern}
\usepackage{amsmath,amssymb}
\usepackage{booktabs}
\usepackage{listings}
\usepackage{xcolor}
\usepackage{tikz}
\usepackage{microtype}
\usepackage{hyperref}
\usepackage{xurl}
\hypersetup{pdftitle={Formal Verification of a Computer-Assisted Proof in Celestial Mechanics},pdfauthor={Tejasvi Singh Tomar}}
\emergencystretch=1.5em
\renewcommand{\topfraction}{0.9}
\renewcommand{\floatpagefraction}{0.8}
\renewcommand{\textfraction}{0.1}
\renewcommand\UrlFont{\color{blue}\rmfamily}
\urlstyle{rm}

\spnewtheorem*{thmA}{Theorem A}{\bfseries}{\itshape}
\spnewtheorem*{thmB}{Theorem B}{\bfseries}{\itshape}
\providecommand{\doi}[1]{\url{https://doi.org/#1}}

\newcommand{\code}[1]{\texttt{#1}}
\newcommand{\tr}{\operatorname{tr}}
\newcommand{\R}{\mathbb R}
\newcommand{\Odir}{\Omega_{\mathrm{dir}}}
\newcommand{\Och}{\Omega_{\mathrm{ch}}}

\lstdefinelanguage{lean}{
  keywords={def,theorem,structure,where,class,extends,macro,let,fun,match,
    with,by,abbrev,instance,noncomputable,if,then,else,deriving,variable,
    set_option,in},
  sensitive=true,
  comment=[l]{--},
  morecomment=[s]{/-}{-/},
}
\lstset{
  language=lean,
  basicstyle=\ttfamily\scriptsize,
  keywordstyle=\bfseries,
  commentstyle=\itshape\color{black!65},
  columns=fullflexible,
  keepspaces=true,
  captionpos=t,
  frame=tb,
  framesep=3pt,
  aboveskip=6pt,
  belowskip=4pt,
  escapeinside={(*}{*)},
  upquote=true,
}

\begin{document}

\title{Formal Verification of a Computer-Assisted\\Proof in Celestial Mechanics}
\titlerunning{Formal Verification of a Proof in Celestial Mechanics}
\author{Tejasvi Singh Tomar\orcidID{0000-0003-4668-3639}}
\authorrunning{T.\,S. Tomar}
\institute{}

\maketitle

\begin{abstract}
Albouy, Fu and Sun conjectured in 2008 that four positive masses admit, for each cyclic order of the bodies, exactly one convex central configuration up to similarity; the question is Problem~10 of a 2012 list of open problems on the $n$-body problem. We recently proved the conjecture by reducing it to an inequality for a trace at the real zeros of an equation in four variables; a branch and bound that starts from a grid of 136{,}448 boxes establishes the inequality. In this paper we describe a formalization of the whole proof in Lean~4 and Mathlib in which Lean's kernel checks the computation, using \code{decide +kernel}, so that no compiler has to be trusted. Our checker stores an interval as two natural numbers in a biased fixed-point format. The kernel can then hand all the arithmetic to GMP. Functions are enclosed over boxes by first-order slopes. An unverified search records the decisions for each box in one natural number, and the checker replays them. We write each formula once over a type class; evaluating it at a type that pairs a real number with an interval then proves the enclosure correct. Checking the 3418 certificate theorems takes 7.6 hours of processor time, and the independent type checker nanoda accepts the proof. A version rebuilt on one eighth of the boxes passed three type checkers in 2.9 hours on a public registry of formal proofs.

\keywords{Lean 4 \and Interval arithmetic \and Proof by reflection \and Central configurations \and Computer-assisted proof}
\end{abstract}

\section{Introduction}

A central configuration (CC) of $n$ point masses is a configuration in which the acceleration of each body under Newtonian gravity points towards the centre of mass and is proportional to the distance from it. CCs are the shapes of the solutions of the $n$-body problem that rotate rigidly, and they are a basic object of celestial mechanics~\cite{Moe15}. For four bodies, MacMillan and Bartky~\cite{MB32} showed that for any positive masses and any cyclic order of the bodies there is a convex CC with the bodies in that order; Xia~\cite{Xia04} gave a shorter proof. Albouy, Fu and Sun~\cite{AFS08} conjectured that this convex CC is unique up to similarity, and the conjecture is Problem~10 in the list of Albouy, Cabral and Santos~\cite{ACS12}. It had been proved for four equal masses~\cite{Alb96}, for equal masses at the ends of a diagonal~\cite{LS02,PCS07,AFS08,Lea03}, and, with interval arithmetic, when each mass is at least $0.6$ times the largest~\cite{SXY23}.

In~\cite{Tom26} we proved the conjecture for all masses. We bounded the Hessian of the potential at a convex CC from below by a quarter of a form that vanishes only on rotations and translations (Theorem~B in Sect.~\ref{sec:math}); a covering argument then turns this nondegeneracy into uniqueness (Theorem~A). An exact identity for the Hessian and the Cauchy--Schwarz inequality reduce Theorem~B to an inequality, $\tr S(y)<3/4$ for some $y\in\R^2$, that has to hold at each real zero of an explicit function $P$ of four variables in a bounded region. We proved this inequality by a branch and bound in Python that uses Arb ball arithmetic~\cite{Joh17}. The search starts from a grid of 136{,}448 boxes and takes about four minutes to run on 32 hardware threads.

A proof of this kind asks the reader to trust the program and the libraries below it, and to accept that the formulas in the code are those in the paper. To remove these assumptions, we formalized the whole proof in Lean~4~\cite{dMU21} and Mathlib~\cite{mC20}. The mathematical part is a conventional Mathlib development. For the computation we follow the reflection approach~\cite{Bou97,BB02}: we define a checker as a function in the logic, prove it sound once, and let the kernel run it on each certificate. Lean offers two ways to run such a checker. The tactic \code{native\_decide} compiles it and executes machine code, but the proof then depends on the axiom \code{Lean.ofReduceBool}, through which the compiler and the run-time system have to be trusted. The tactic \code{decide +kernel} instead has the kernel evaluate the Boolean term, and it adds no axiom. We use only \code{decide +kernel}. The kernel evaluates terms by rewriting. However, it hands fourteen operations on natural-number literals to GMP, among them addition, multiplication, division, comparisons, bitwise operations and shifts, and we wrote the checker so that all of its arithmetic is done by these operations.

In this paper we report what it took to make a nonlinear computer-assisted proof of this size checkable by Lean's kernel, and what the checks cost. Our checker stores an interval as two natural numbers that are both upper bounds, in a biased fixed-point format, so that each step of the interval arithmetic is an accelerated operation. We write the formulas of the proof once over a type class and evaluate them at the reals in the proofs and at intervals and slope jets in the checker. A carrier type that pairs a real number with its enclosure turns one evaluation into the soundness proof of a formula. The search is done by an unverified program that writes its decisions into a natural number; the checker replays them. We also describe the problems we met in keeping the proofs about the checker cheap for the kernel and for a second type checker, and how we cut the certificates so that the memory needed to check them stays bounded. On 32 hardware threads the certificate modules build in 19 minutes, or 7.6 hours of processor time, and the independent type checker nanoda also accepts the proof.

Section~\ref{sec:math} summarizes the mathematics and Sect.~\ref{sec:statement} the formal statement. Section~\ref{sec:checker} describes the checker and Sect.~\ref{sec:results} the checks, including a version that fits the limits of a public registry. Section~\ref{sec:related} discusses related work, and Sect.~\ref{sec:discussion} the limits of the design.

\section{The Theorem and Its Computer-Assisted Proof}\label{sec:math}

\subsection{Central Configurations}

Let $m_1,\dots,m_4>0$ be masses at distinct points $q_1,\dots,q_4\in\R^2$. We write $r_{ij}=|q_i-q_j|$, $U=\sum_{i<j}m_im_j/r_{ij}$, $M=\sum_im_i$, $c=M^{-1}\sum_im_iq_i$ and $I=\sum_im_i|q_i-c|^2$. The configuration $q$ is a CC if there is $\lambda\in\R$ with
\begin{equation}\label{eq:cc}
\nabla_{q_i}U+\lambda m_i(q_i-c)=0\qquad(i=1,\dots,4);
\end{equation}
then $\lambda=U/I>0$. The CCs are the critical points of $UI^{1/2}$, a function invariant under similarities. At a CC, we fix $\lambda$ and define two quadratic forms on $v\in(\R^2)^4$,
\[
Q(v)=D^2\bigl(U+\tfrac\lambda2 I\bigr)(q)[v,v],\qquad
K(v)=\sum_{i<j}\frac{3m_im_j}{r_{ij}^3}\,\dot r_{ij}^2,
\]
where $\dot r_{ij}=\langle q_i-q_j,v_i-v_j\rangle/r_{ij}$ is the derivative of $r_{ij}$ in the direction~$v$. At a convex configuration, $K(v)=0$ exactly when $v$ is an infinitesimal rotation or translation.

\begin{thmB}[\cite{Tom26}]
At every convex CC of positive masses, $Q\ge K/4$.
\end{thmB}

\begin{thmA}[\cite{Tom26}]
For all positive masses $m_1,\dots,m_4$ there is, up to similarity, exactly one convex CC with the bodies in the cyclic order $1234$. It is a nondegenerate local minimum of $UI^{1/2}$ on the space of shapes, and it depends analytically on the masses.
\end{thmA}

Theorem~A follows from Theorem~B by a covering argument. We fix the similarity by placing $q_1$ at $(0,0)$ and $q_2$ at $(1,0)$, fix the orientation by the signs of the triangle areas, and let $Z$ be the set of pairs $(m,q)$ of positive masses and a convex CC $q$ for $m$ in this slice, with the cyclic order $1234$. By Theorem~B and the implicit function theorem, the projection $Z\to(0,\infty)^4$ is a local homeomorphism. It is also proper: a limit of such CCs with masses in a compact set has no collision, by a lemma of Shub~\cite{Shu71} applied with varying masses, and it is not collinear, by Theorem~B again. Since $(0,\infty)^4$ is connected, all fibres then have the same number of points, and the fibre over four equal masses is the square alone~\cite{Alb96},~\cite[Theorem~1]{AFS08}. The argument needs only $Q\ge\varepsilon K$ for some $\varepsilon>0$; Sect.~\ref{sec:discussion} returns to this.

\subsection{From the Hessian to a Bounded Computation}

At a convex CC, Dziobek's equations~\cite{Dzi00} give $m_im_j(r_{ij}^{-3}-\lambda')=\sigma A_iA_j$ for all $i<j$, where $\lambda'=\lambda/M$, $A_l$ is the signed area of the triangle formed by the three bodies other than $l$, and $\sigma<0$. From them we obtain the identity
\begin{equation}\label{eq:identity}
Q=K-|\sigma|\,|L|^2,\qquad L=\sum_lA_lv_l\in\R^2 .
\end{equation}
The six rates $\dot r_e$, $e=ij$, satisfy one linear relation $\sum_e\tau_e\dot r_e=0$, with $\tau_{ij}=A_iA_jr_{ij}$, and hence $L=\sum_e(\beta_e+\tau_ey)\dot r_e$ for all $y\in\R^2$, with vectors $\beta_e\in\R^2$ that depend only on the shape. The Cauchy--Schwarz inequality with the weights of $K$ gives
\begin{equation}\label{eq:majorant}
|\sigma|\,|L|^2\le\lambda_{\max}(S(y))\,K\le\tr S(y)\,K,
\end{equation}
where $S(y)=\sum_eD_e(\beta_e+\tau_ey)(\beta_e+\tau_ey)^{\mathsf T}$ and $D_{ij}=|\sigma|r_{ij}^3/(3m_im_j)=-(r_{ij}^{-3}-\lambda')\,r_{ij}^3/(3A_iA_j)$. On convex CCs, $\lambda'=(s_{12}s_{34}-s_{13}s_{24})/(s_{12}+s_{34}-s_{13}-s_{24})$ with $s_e=r_e^{-3}$. Thus $S(y)$ does not depend on the masses. By~\eqref{eq:identity} and~\eqref{eq:majorant}, it is therefore enough to show that each convex CC has a $y$ with $\tr S(y)<3/4$. The minimum of $\tr S(y)$ over $y$ has supremum about $0.593$ over convex CCs, and the constant $3/4$ leaves a margin of about $0.157$.

\subsection{Coordinates and Regions}

\begin{figure}[t]
\centering
\begin{tikzpicture}[scale=2.5, every node/.style={font=\small}]
  \coordinate (O) at (0,0);
  \coordinate (q1) at (1,0);
  \coordinate (q2) at (0.18,0.8818);
  \coordinate (q3) at (-0.5,0);
  \coordinate (q4) at (-0.08,-0.3919);
  \draw[thick] (q1)--(q2)--(q3)--(q4)--cycle;
  \draw[dashed] (q1)--(q3) (q2)--(q4);
  \draw (0.22,0) arc[start angle=0,end angle=78.46,radius=0.22];
  \node at (0.15,0.1) {$\theta$};
  \fill (q1) circle (0.025) node[right] {$q_1$};
  \fill (q2) circle (0.025) node[above] {$q_2$};
  \fill (q3) circle (0.025) node[left] {$q_3$};
  \fill (q4) circle (0.025) node[below] {$q_4$};
  \fill (O) circle (0.015);
  \node[below] at (0.6,0) {$1$};
  \node[below] at (-0.25,0) {$b$};
  \node[left] at (0.11,0.5) {$a$};
  \node[right] at (-0.04,-0.21) {$c$};
\end{tikzpicture}
\caption{The coordinates of~\cite{CCR19}, drawn at $a=0.9$, $b=0.5$, $c=0.4$ and $x=\cos\theta=0.2$, a point where all six $g_i$ are positive. The diagonals meet at the origin.}\label{fig:coords}
\end{figure}

We use the coordinates of Corbera, Cors and Roberts~\cite{CCR19}. We label the bodies so that $q_1q_2$ is a longest side. At a convex CC the opposite side $q_3q_4$ is then a shortest side, and both diagonals are longer than all four sides. A similarity brings the configuration to the form
\[
q_1=(1,0),\quad q_2=a(x,s),\quad q_3=(-b,0),\quad q_4=-c(x,s),\qquad s=\sqrt{1-x^2},
\]
with the diagonals meeting at the origin (Fig.~\ref{fig:coords}). The squared distances $R_e=r_e^2$ are polynomials, such as $R_{12}=1+a^2-2ax$ and $R_{24}=(a+c)^2$, and the order of the six distances becomes six inequalities $g_i\ge0$, each affine in $x$; for instance $g_1=R_{13}-R_{12}=b^2+2b-a^2+2ax$. If we eliminate $\lambda'$ from Dziobek's equations, one equation remains,
\begin{equation}\label{eq:P}
P=\eta_{14}\eta_{23}(\eta_{13}+\eta_{24}+\eta_{34})+\eta_{13}\eta_{24}(\eta_{14}+\eta_{23}-\eta_{34})=0,
\end{equation}
where $\eta_e=|s_e-s_{12}|$. Our checker computes $\eta_e$ as $R_e-R_{12}$ or $R_{12}-R_e$, whichever is nonnegative, written as a combination of the $g_i$ with nonnegative coefficients, times $\psi(R_e,R_{12})$ with $\psi(X,Y)=(X+\sqrt{XY}+Y)/((\sqrt X+\sqrt Y)(XY)^{3/2})$, a form without cancellation. The points with $g\ge0$ and $P=0$ are exactly the convex CCs, up to relabelling and similarity~\cite[Lemma~3.1]{CCR19}.

A priori bounds~\cite[Lemma~3.2]{CCR19} put every convex CC with $b\ge1/8$ in the box $\Odir=[\frac12,\frac74]\times[\frac18,1]\times[0,\frac74]$ of the coordinates $(a,b,c)$. As $b\to0$, the bodies $3$ and $4$ collide at the origin, which forms an equilateral triangle with $q_1$ and $q_2$, and $m_3$ and $m_4$ tend to $0$ like $b^3$. For $b\le1/8$ we use the chart $a=1+b\alpha$, $c=b\gamma$, $x=\frac12+b\xi$, in which the bounds give $(b,\alpha,\gamma)\in\Och=[0,\frac18]\times[-2,1]\times[0,3]$, and in which $\hat P=bP$ extends to $b=0$ with $\partial\hat P/\partial\xi<0$ there. The computation has to establish two claims:
\begin{itemize}
\item[(C1)] if $(a,b,c)\in\Odir$, $-1\le x\le1$, $g_i(a,b,c,x)\ge0$ for all $i$ and $P(a,b,c,x)=0$, then $\tr S(y)<3/4$ for some $y\in\R^2$;
\item[(C2)] the same holds in the chart, for $(b,\alpha,\gamma)\in\Och$, with $\hat P$ in place of $P$.
\end{itemize}

\subsection{Why the Computation Is Awkward to Formalize}\label{sec:awkward}

The claims are not inequalities over boxes. They quantify over the zeros of $P$ in $x$ and ask for a witness $y$ at each zero. The checker therefore has to localize the zeros in each box, and we need a $y$ that works on a whole piece; we find it in floating point and only have to verify it in the checker. The formulas also divide and take square roots. Division and square root in Mathlib are total, with $x/0=0$ and $\sqrt x=0$ for $x<0$, so we have to relate the checker's enclosures to total functions. Finally, the volume is large for a kernel that evaluates by rewriting: the Python search evaluated $P$ 17.8~million times.

\section{The Formal Statement}\label{sec:statement}

\begin{lstlisting}[float=tp,caption={The definitions behind the main theorems and their statements (from \code{C4/\allowbreak Defs.lean}, \code{ShapeHess.lean}, \code{TheoremB.lean} and \code{TheoremA.lean}). The variables \mbox{\code{m : Masses}} and \mbox{\code{q : Conf}} are section variables.},label=lst:defs]
def IsCC : Prop :=
  CollisionFree q (*$\wedge$*)
    (*$\exists$*) lam : (*$\mathbb{R}$*), (*$\forall$*) i, ((*$\sum$*) j, (m i * m j * ss q i j) (*$\bullet$*) (q j - q i)) +
      (lam * m i) (*$\bullet$*) (q i - cm m q) = 0
noncomputable def hessK (v : Conf) : (*$\mathbb{R}$*) :=
  esum fun i j => 3 * m i * m j * ss q i j * dr q v i j ^ 2
noncomputable def hessQ (v : Conf) : (*$\mathbb{R}$*) :=
  hessK m q v -
    esum fun i j => m i * m j * wgeo m q i j * dot (v i - v j) (v i - v j)
def Order1234 (q : Conf) : Prop :=
  0 < area q 0 * area q 2 (*$\wedge$*) 0 < area q 1 * area q 3 (*$\wedge$*) area q 0 * area q 1 < 0
def Similar (q q' : Conf) : Prop :=
  (*$\exists$*) a b : (*$\mathbb{R}$*), (*$\exists$*) t : V2, a ^ 2 + b ^ 2 (*$\ne$*) 0 (*$\wedge$*)
    (q' = simc a b t q (*$\vee$*) q' = simc a b t (mirror q))
theorem hessian_Phi {m : Masses} (hm : (*$\forall$*) i, 0 < m i) {q : Conf}
    (hq : CollisionFree q) (v : Conf) :
    fderiv (*$\mathbb{R}$*) (fderiv (*$\mathbb{R}$*) (fun p : Conf => Upot m p + lamC m q / 2 * Iner m p))
        q v v = hessQ m q v
theorem theoremB (m : Masses) (hm : (*$\forall$*) i, 0 < m i) (q : Conf) (hcc : IsCC m q)
    (hconv : IsConvex q) (v : Conf) : hessK m q v / 4 (*$\le$*) hessQ m q v
theorem theoremA (m : Masses) (hm : (*$\forall$*) i, 0 < m i) :
    (*$\exists$*) q, IsCC m q (*$\wedge$*) Order1234 q (*$\wedge$*) (*$\forall$*) q', IsCC m q' (*$\to$*) Order1234 q' (*$\to$*) Similar q q'
\end{lstlisting}

In Listing~\ref{lst:defs} we show the definitions that a reader has to check and the statements of Theorems~A and~B. Positions are in \code{V2}${}=\R\times\R$ and configurations are maps from \code{Fin 4} to \code{V2}, so the bodies $1,\dots,4$ of Sect.~\ref{sec:math} are \code{0},~\dots,~\code{3} here. In \code{IsCC}, \code{ss q i j} is $r_{ij}^{-3}$ and \code{cm m q} is the centre of mass. \code{esum} sums over the six pairs, \code{dr} is $\dot r_{ij}$, and \code{wgeo} is $r_{ij}^{-3}-\lambda/M$ with $\lambda=U/I$. The form \code{hessQ} is an explicit formula, and \code{hessian\_Phi} proves that it is the second Fr\'echet derivative of $U+\frac\lambda2I$, so the formula does not have to be trusted. \code{Order1234} is stated with signed areas, and \code{order1234\_iff} proves it equivalent to the crossing of the diagonals $q_1q_3$ and $q_2q_4$. \code{Similar} allows reflections. We also prove uniqueness for every cyclic order, the analytic dependence on the masses, and the nondegeneracy as a local minimum on the shape space, which we construct as a compact real-analytic manifold of dimension~4.

For the main theorems, \code{\#print axioms} lists \code{propext}, \code{Classical.choice} and \code{Quot.sound}; for the two theorems that collect the certificates (Sect.~\ref{sec:certs}) it lists \code{propext} alone. No source file uses \code{sorry}, \code{admit}, \code{axiom} declarations, \code{native\_decide}, \code{implemented\_by}, \code{extern}, \code{unsafe} or \code{partial}.

The Mathlib part has 39 files with 12{,}790 lines, and 1{,}900 of these lines prove the checker sound. The checker \code{C4Check} has 941 lines in five files: \code{Num} (162 lines) for numbers, \code{Jet} (158) for jets and the class \code{Ops}, \code{Formulas} (300), \code{Search} (212) and \code{Grid} (109). We kept it free of Mathlib: it uses core Lean only, and the certificate modules import \code{C4Check} alone, so checking them never loads Mathlib.

\section{The Checker}\label{sec:checker}

\subsection{Evaluation by the Kernel}\label{sec:eval}

The tactic \code{decide +kernel} proves a proposition from its \code{Decidable} instance by having the kernel reduce the instance to \code{isTrue}. Each of our certificates states \code{f h = true} for a Boolean function \code{f} and a literal \code{h}. To check it, the kernel evaluates \code{f h} by weak-head reduction: it unfolds definitions, reduces matches and recursors, and calls GMP when one of \code{Nat.add}, \code{sub}, \code{mul}, \code{div}, \code{mod}, \code{gcd}, \code{beq}, \code{ble}, \code{land}, \code{lor}, \code{xor}, \code{shiftLeft}, \code{shiftRight} and \code{pow} is applied to two literals. It has no other arithmetic. In particular, \code{Int} has no accelerated operations of its own, and a \code{Nat} operation whose arguments do not reduce to literals, as happens in a proof with free variables, unfolds to the unary definition of \code{Nat}. Both facts shaped the design that we describe in the rest of this section.

\subsection{Numbers}\label{sec:num}

\begin{lstlisting}[float=tp,caption={Fixed-point intervals (from \code{C4Check/\allowbreak Num.lean}, abridged).},label=lst:num]
def Bb : Nat := 2 ^ 256
def C1 : Nat := 2 ^ 512 + 2 ^ 352      -- Bb^2 + Bb 2^96
def RND : Nat := 2 ^ 96 - 1
macro "sel " b:term:max t:term:max e:term:max : term =>
  `(Bool.casesOn (motive := fun _ => _) $b $e $t)
structure I where
  L : Nat
  H : Nat
/-- `(X - Bb) (Y - Bb) + Bb 2^96`, truncated at `0` -/
def pp (X Y : Nat) : Nat :=
  Nat.sub (Nat.add (Nat.mul X Y) C1) (Nat.shiftLeft (Nat.add X Y) 256)
def up (P : Nat) : Nat := Nat.shiftRight (Nat.add P RND) 96
def add (A B : I) : I :=
  (*$\langle$*)Nat.sub (Nat.add A.L B.L) Bb, Nat.sub (Nat.add A.H B.H) Bb(*$\rangle$*)
def neg (A : I) : I := (*$\langle$*)A.H, A.L(*$\rangle$*)
/-- the product, by the signs of the factors (`L (*$\le$*) Bb` means that every element
is `(*$\ge$*) 0`, `H (*$\le$*) Bb` that every element is `(*$\le$*) 0`) -/
def mul (A B : I) : I :=
  sel (Nat.ble A.L Bb)
    (sel (Nat.ble B.L Bb) (*$\langle$*)up (pm A.L B.L), up (pp A.H B.H)(*$\rangle$*)
      (sel (Nat.ble B.H Bb) (*$\langle$*)up (pp A.H B.L), up (pm A.L B.H)(*$\rangle$*)
        (*$\langle$*)up (pp A.H B.L), up (pp A.H B.H)(*$\rangle$*)))
    (sel (Nat.ble A.H Bb) (*\textrm{\dots}*) (*\textrm{\dots}*))     -- six more cases
def isqD (n : Nat) : Nat := let r := isq n; sel (Nat.ble (Nat.mul r r) n) r 0
def isqU (n : Nat) : Nat :=
  let s := Nat.succ (isq n); sel (Nat.ble n (Nat.mul s s)) s n
def sqrt (A : I) : Option I :=
  sel (Nat.ble A.L Bb)
    (some (*$\langle$*)Nat.sub Bb (isqD (Nat.shiftLeft (Nat.sub Bb A.L) 96)),
      Nat.add Bb (isqU (Nat.shiftLeft (Nat.sub A.H Bb) 96))(*$\rangle$*))
    none
\end{lstlisting}

Since \code{Int} is not accelerated, we represent a real number by a natural number $X$ with value $(X-2^{256})/2^{96}$. Values thus start at $-2^{160}$ and have absolute precision $2^{-96}$, and the bias keeps all quantities natural numbers. An interval $\langle L,H\rangle$ contains $v$ if $-v\le\mathrm{val}(L)$ and $v\le\mathrm{val}(H)$ (Listing~\ref{lst:num}). Since both components are upper bounds, of $-v$ and of $v$, one rounding direction is enough: we always round up. The truncated subtraction of \code{Nat} can only increase a result, so it is in fact safe wherever it occurs. We do not need to maintain an invariant such as $\mathrm{val}(L)+\mathrm{val}(H)\ge0$, because an interval whose bounds cross contains no number, and statements about its elements then hold vacuously.

We form products directly on the codes. For $X=2^{256}+x\,2^{96}$ and $Y=2^{256}+y\,2^{96}$,
\[
XY+2^{512}+2^{352}-2^{256}(X+Y)=2^{96}\bigl(2^{256}+xy\,2^{96}\bigr),
\]
and therefore \code{pp X Y} followed by \code{up}, a shift by 96 bits that rounds up, gives the code of $xy$. The function \code{mul} chooses among nine cases by the signs of the factors~\cite{MKC09}, and the tests are written with \code{sel}, a macro for a bare \code{Bool.casesOn}. To compute a reciprocal, we divide $2^{192}$ by the distance of a code from $2^{256}$ with \code{Nat.div}, rounding down or up as the bound requires. Integer square roots are computed by Newton's iteration from above, starting at a power of two that is found by bisection on the exponent. The result is checked by squaring it once. If this check fails, the checker falls back to a bound that is valid but useless ($0$ below, the radicand above). Division and the square root return an \code{Option I}, which is \code{none} when the divisor may vanish or the radicand may be negative.

\subsection{Formulas Written Once}\label{sec:carrier}

\begin{lstlisting}[float=tp,caption={The class \code{Ops}, a generic formula, and the carrier \code{EI} (from \code{C4Check/\allowbreak Jet.lean}, \code{C4Check/\allowbreak Formulas.lean}, \code{C4/\allowbreak Carrier.lean} and \code{C4/\allowbreak ModeSound.lean}, abridged). In \code{Formulas.lean}, $\alpha$ has an instance \code{[Ops $\alpha$]}.},label=lst:ops]
class Ops ((*$\alpha$*) : Type) extends Add (*$\alpha$*), Sub (*$\alpha$*), Mul (*$\alpha$*), Div (*$\alpha$*), Neg (*$\alpha$*) where
  sqrt : (*$\alpha$*) (*$\to$*) (*$\alpha$*)
  sq : (*$\alpha$*) (*$\to$*) (*$\alpha$*)
  ofInt : Int (*$\to$*) (*$\alpha$*)
instance : Ops (Option I) where
  div x y := match x, y with | some a, some b => a.div b | _, _ => none
  (*\textrm{\dots}*)
/-- `hd X Y = (Y^(-3/2) - X^(-3/2)) / (X - Y)`, in a form without cancellation -/
def hd (X Y : (*$\alpha$*)) : (*$\alpha$*) :=
  let sx := Ops.sqrt X
  let sy := Ops.sqrt Y
  (X + sx * sy + Y) / ((sx + sy) * X * sx * Y * sy)
def dirP (a b c x : (*$\alpha$*)) : (*$\alpha$*) :=
  let g := dirG a b c x
  let R12 := 1 + a * a - 2 * a * x
  (*\textrm{\dots}*)
  E14 * E23 * (Dl13 + E34 + Dl24) + Dl13 * Dl24 * (E14 + E23 - E34)
def OEnc (v : (*$\mathbb{R}$*)) (o : Option I) : Prop := (*$\forall$*) A, o = some A (*$\to$*) I.mem v A
structure EI where
  r : (*$\mathbb{R}$*)
  i : Option I
  ok : OEnc r i
@[instance_reducible] def eiOps : Ops EI where
  mul x y := (*$\langle$*)x.r * y.r, x.i * y.i, by
    rcases x with (*$\langle$*)rx, ix, hx(*$\rangle$*); rcases y with (*$\langle$*)ry, iy, hy(*$\rangle$*)
    intro A h
    cases ix <;> cases iy <;> cases h
    exact I.mem_mul (hx _ rfl) (hy _ rfl)(*$\rangle$*)
  (*\textrm{\dots}*)
theorem dirP_enc :
    OJEnc S zc (fun z => dirPR (z 0) (z 1) (z 2) (z 3))
      (@dirP _ (jetOps R) o0 o1 o2 o3) :=
  (@dirP (EJ S zc) (ejOps R hR)
    (*$\langle$*)_, o0, h0(*$\rangle$*) (*$\langle$*)_, o1, h1(*$\rangle$*) (*$\langle$*)_, o2, h2(*$\rangle$*) (*$\langle$*)_, o3, h3(*$\rangle$*)).ok
\end{lstlisting}

We write the formulas of Sect.~\ref{sec:math}, that is, the $g_i$, $P$, the coefficients $D_e$, $\tau_e$, $\beta_e$, and $\tr S(y)$, only once, in \code{Formulas.lean}, for an arbitrary type with an instance of the class \code{Ops} (Listing~\ref{lst:ops}). Integer literals are turned into \code{Ops.ofInt} by a low-priority \code{OfNat} instance. The same definitions can then be evaluated at four types: at $\R$ with the instance \code{opsReal} in the proofs, at \code{Option I} and at \code{Option J}, the jets of Sect.~\ref{sec:jets}, in the checker, and at \code{Float} in numerical tests. In \code{Option I}, an undefined operation is recorded as \code{none}. We chose not to make \code{opsReal} a global instance, because with it in scope ordinary real arithmetic would be elaborated through \code{Ops}. The Mathlib modules \code{Chart} and \code{NormalForm} prove that the real evaluations equal the geometric quantities of Sect.~\ref{sec:math}. As a result, no formula is transcribed by hand between the mathematics and the checker.

We obtain the soundness of the checker's evaluations from carrier types. An element of \code{EI} consists of a real number, an \code{Option I}, and a proof that the option, when it is \code{some A}, contains the number. The instance \code{eiOps} performs each operation on both components and proves the enclosure with a lemma of \code{NumSound}; Listing~\ref{lst:ops} shows multiplication. If we evaluate a generic formula at \code{EI}, the first component of the result is the real formula and the second is the checker's computation, both by definition. The third component is then the enclosure theorem for that formula. The carrier \code{EJ} does the same for jets. This is why the enclosure of $P$ over a box has a proof of one line (\code{dirP\_enc}), and this proof does not change when the formula does. The real side needs no side conditions: division and square root on $\R$ are total, and the enclosure is claimed only when the checker returns \code{some}. For a quotient, this happens only if the divisor's interval excludes~$0$.

The usual alternative, used for instance in CoqInterval~\cite{Mel08,MDM16}, is to reify expressions into a syntax tree and to prove an evaluator of the tree correct by induction. Carriers need neither a syntax nor reification. The price is that every formula has to be written against \code{Ops}, also on the Mathlib side, where we have to supply \code{opsReal} explicitly.

\subsection{Slope Jets}\label{sec:jets}

\begin{lstlisting}[float=tp,caption={Slope jets (from \code{C4Check/\allowbreak Jet.lean}, abridged).},label=lst:jet]
structure J where
  U : I
  C : I
  s0 : I
  s1 : I
  s2 : I
  s3 : I
/-- centred form `C + (*$\Sigma$*)_i s_i [-r_i, r_i]` -/
def cf (R : Rad) (C s0 s1 s2 s3 : I) : I :=
  let E := Nat.add (Nat.add (I.eterm s0 R.r0) (I.eterm s1 R.r1))
    (Nat.add (I.eterm s2 R.r2) (I.eterm s3 R.r3))
  (*$\langle$*)Nat.add C.L E, Nat.add C.H E(*$\rangle$*)
def range (R : Rad) (u : J) : I := u.U.inter (cf R u.C u.s0 u.s1 u.s2 u.s3)
/-- `u v - (*$\bar u\bar v$*) = u (v - (*$\bar v$*)) + (*$\bar v$*) (u - (*$\bar u$*))` -/
def mul (R : Rad) (u v : J) : J :=
  let s0 := (u.U.mul v.s0).add (v.C.mul u.s0)
  (*\textrm{\dots}*)
  let s3 := (u.U.mul v.s3).add (v.C.mul u.s3)
  let C := u.C.mul v.C
  (*$\langle$*)(u.U.mul v.U).inter (cf R C s0 s1 s2 s3), C, s0, s1, s2, s3(*$\rangle$*)
\end{lstlisting}

Plain interval evaluation of $P$ over a box of side $1/32$ overestimates its range badly, because each variable occurs many times in $P$. To reduce the overestimation, we use first-order centred forms with slopes~\cite{KN85}. A jet over a box $Z\subset\R^4$ with centre $\bar z$ consists of intervals $U$, $C$ and $S_0,\dots,S_3$; the coordinates of $Z$ are the three coordinates of the cell and $x$. The jet encloses a function $u$ if $u(z)\in U$ for $z\in Z$, $u(\bar z)\in C$, and for each $z\in Z$ there are $s_i\in S_i$ with $u(z)-u(\bar z)=\sum_is_i(z_i-\bar z_i)$. The rules for the operations are algebraic identities, such as $uv-\bar u\bar v=u(v-\bar v)+\bar v(u-\bar u)$ (Listing~\ref{lst:jet}) and $\sqrt u-\sqrt{\bar u}=(u-\bar u)/(\sqrt u+\sqrt{\bar u})$. As a consequence, their soundness proofs in \code{JetSound} need no mean value theorem. Each operation intersects the plain interval result with the centred form $C+\sum_iS_i[-r_i,r_i]$, where $r_i$ bounds $|z_i-\bar z_i|$. Slopes also give us the interval Newton step in $x$: at a zero, $0=P(\bar z)+\sum_{i<3}s_i(z_i-\bar z_i)+s_3(x-\bar x)$, and this equation can be used to confine $x$ when $S_3$ excludes~$0$.

\subsection{Hints}\label{sec:hints}

\begin{lstlisting}[float=tp,caption={Replaying a hint (from \code{C4Check/\allowbreak Search.lean} and \code{C4/\allowbreak SearchSound.lean}, abridged).},label=lst:search]
def locF (M : Mode) (B : Box)
    (rec : Nat (*$\to$*) List I (*$\to$*) List I (*$\to$*) Option (List I (*$\times$*) Nat)) (h : Nat)
    (X : I) (work kept : List I) : Option (List I (*$\times$*) Nat) :=
  sel (Nat.ble (Nat.land h 3) 1)
    (sel (Nat.beq (Nat.land h 3) 0)
      ((pstep M B X).bind fun X' =>
        sel (I.isEmpty X') (rec (Nat.shiftRight h 2) work kept) none)
      (rec (Nat.shiftRight h 2) work (X :: kept)))
    (sel (Nat.beq (Nat.land h 3) 2)
      ((pstep M B X).bind fun X' => rec (Nat.shiftRight h 2) (X' :: work) kept)
      (rec (Nat.shiftRight h 2)
        ((I.halves X).1 :: (I.halves X).2 :: work) kept))
def locH (M : Mode) (B : Box) :
    Nat (*$\to$*) Nat (*$\to$*) List I (*$\to$*) List I (*$\to$*) Option (List I (*$\times$*) Nat)
  | _, h, [], kept => some (kept, h)
  | 0, _, _ :: _, _ => none
  | fuel + 1, h, X :: work, kept => locF M B (locH M B fuel) h X work kept
def leafH (M : Mode) (B : Box) (h : Nat) : Option Nat :=
  (xrange M B).bind fun X =>
    sel (I.isEmpty X) (some h)
      ((locH M B 4096 h [X] []).bind fun r => certAllH M B r.1 r.2)
def checkBoxH (M : Mode) : Nat (*$\to$*) Box (*$\to$*) Nat (*$\to$*) Option Nat
  | 0, _, _ => none
  | fuel + 1, B, h =>
    sel (Nat.beq (Nat.land h 3) 0) (leafH M B (Nat.shiftRight h 2))
      ((checkBoxH M fuel (splitBox B (Nat.land h 3)).1 (Nat.shiftRight h 2)).bind
        fun h' => checkBoxH M fuel (splitBox B (Nat.land h 3)).2 h')
theorem sel_pos {(*$\alpha$*) : Sort _} {b : Bool} (t e : (*$\alpha$*)) (h : b = true) :
    Bool.casesOn (motive := fun _ => (*$\alpha$*)) b e t = t := by
  subst h; rfl
def AllGood (M : Mode) (Sm : Sem) (B : Box) : Prop :=
  (*$\forall$*) z0 z1 z2 x, B.mem z0 z1 z2 (*$\to$*) (M.clamp = true (*$\to$*) -1 (*$\le$*) x (*$\wedge$*) x (*$\le$*) 1) (*$\to$*)
    ((*$\forall$*) i, 0 (*$\le$*) (Sm.G z0 z1 z2 x).get i) (*$\to$*) Sm.P z0 z1 z2 x = 0 (*$\to$*)
    (*$\exists$*) y0 y1, trSR (Sm.C z0 z1 z2 x) y0 y1 < 3 / 4
theorem checkBoxH_sound (hM : ModeSound M Sm) : (*$\forall$*) (fuel : Nat) (B : Box)
    (h h' : Nat), checkBoxH M fuel B h = some h' (*$\to$*) AllGood M Sm B
\end{lstlisting}

The search is not part of the checker. An unverified program, \code{Gen.lean}, decides where to split a box, how to localize the zeros of $P$ in $x$, and which $y$ to use. We run \code{Gen} beforehand. It searches in the arithmetic of the checker and records the decisions in a \emph{hint}: a natural number that the checker reads two bits at a time, starting from the least significant end (Listing~\ref{lst:search}). At a box, the symbol $0$ makes the box a leaf, and the symbols $1$, $2$, $3$ split it in half along coordinate $0$, $1$, $2$. At a leaf, the checker first computes the interval of $x$ allowed by the conditions $g_i\ge0$ (they are affine in $x$). Unless this interval is empty, the hint continues with one symbol per step of the localization ($0$ drop, $1$ keep, $2$ Newton step, $3$ bisect). Then, for each kept interval, it gives two 32-bit fields that place $y$ on a grid of step $2^{-24}$ in $[-128,128)^2$. The function \code{checkBoxH} replays the hint and checks each decision. An interval can be dropped only if one step of the localization empties it, i.e., if the enclosure of some $g_i$ is negative, the enclosure of $P$ excludes~$0$, or the Newton step leaves nothing. A kept interval has to satisfy $\tr S(y)<3/4$ in jet arithmetic. A wrong hint can of course make the check fail, but it cannot make it succeed wrongly: the soundness theorem \code{checkBoxH\_sound} holds for all hints. The recursions carry explicit fuel (depth 100 for boxes, 4096 localization steps per leaf). This makes them structural, so that they can be unfolded by the kernel.

\code{Gen} replays each hint before writing it. Its search differs from the Python one in several respects, and on the direct region it splits 3{,}800 times, against 318{,}688 times for Python (Table~\ref{tab:search}). We do not know which of the differences actually accounts for this.

\subsection{Writing for the Kernel}\label{sec:kernel}

The same kernel also checks the soundness proofs, and in these proofs the checker's functions are applied to free variables. This caused us several problems. When the kernel meets a \code{Bool.casesOn} or a \code{match} at the head of a term, it reduces the condition or the discriminant to weak head normal form. If the condition compares an interval computed from variables, the arithmetic unfolds to unary numbers, and that is infeasible for numbers of this size. We dealt with this and with related problems in four ways.

\paragraph{Tests take their interval as an argument.} We put each test on a computed interval into a small definition that takes the interval as a parameter (\code{negH}, \code{ifExcl}, \code{newtonStep}, \code{trOK}, \code{xtight}), and every other \code{match} is on an argument or goes through \code{Option.bind}. In \code{SearchSound} we decide each test by rewriting with \code{sel\_pos} or \code{sel\_neg} (Listing~\ref{lst:search}) from a hypothesis on the condition, so that the kernel never has to evaluate a condition.

\paragraph{Equation lemmas.} With the step of the localization loop inlined into \code{locH}, Lean could not prove the equation lemmas of \code{locH}. We therefore moved the step into a separate definition, \code{locF}, and \code{locH} now recurses structurally on its fuel and calls \code{locF}.

\paragraph{Large constants in the elaborator.} Unlike the kernel, the elaborator refuses to evaluate powers whose exponent exceeds the option \code{exponentiation.threshold} (256 by default). When it elaborates the enclosure theorems for the chart formulas, it meets the constants $2^{352}$ and $2^{512}$ of \code{I.mul} through the literal $1/2$. We raised the threshold to 512 in two theorems of \code{ModeSound}, and to 1024 in \code{NumSound}.

\paragraph{A second type checker.} nanoda (Sect.~\ref{sec:indep}) at first ran out of memory on our proofs, after using more than 80\,GB. To locate the cause, we checked the export in windows of consecutive declarations, with the theorems before each window turned into axioms. Two test theorems that differed only in the name of one bound variable then showed the problem: nanoda accepted one and ran out of memory on the other. It does not identify terms that differ only in the names of bound variables until it has unfolded them, and on the tests of intervals computed from free variables this unfolding again reaches unary arithmetic. The lemmas \code{negH\_of\_ble}, \code{ifExcl\_pos} and \code{ifExcl\_neg} now state these tests once for a variable interval, and we changed the proofs to use them.

\subsection{Certificates}\label{sec:certs}

\begin{lstlisting}[float=tp,caption={Cells, a certificate, and the use of the certificates (from \code{C4Check/\allowbreak Grid.lean}, \code{C4Cert/\allowbreak Dir000.lean}, \code{C4Cert/\allowbreak Cover.lean} and \code{C4/\allowbreak Cap.lean}; hints shortened).},label=lst:cert]
def iv32 (n : Nat) : I :=
  (*$\langle$*)Nat.sub Bb (Nat.mul n (2 ^ 91)), Nat.add Bb (Nat.mul (Nat.succ n) (2 ^ 91))(*$\rangle$*)
def dirBox (i j k : Nat) : Box := (*$\langle$*)iv32 (16 + i), iv32 (4 + j), iv32 k(*$\rangle$*)
def dirCell (n h : Nat) : Bool :=
  (checkBoxH dirMode depth (dirBox (n / 1568) (n / 56 % 28) (n % 56)) h).isSome
def allCells (P : Nat (*$\to$*) Nat (*$\to$*) Bool) : Nat (*$\to$*) Nat (*$\to$*) List Nat (*$\to$*) Bool
  | n, e, [] => Nat.ble e n
  | n, e, h :: hs => sel (P n h) (allCells P (Nat.succ n) e hs) false
def Cover (P : Nat (*$\to$*) Nat (*$\to$*) Bool) (a b : Nat) : Prop :=
  (*$\forall$*) n, a (*$\le$*) n (*$\to$*) n < b (*$\to$*) (*$\exists$*) h, P n h = true
theorem c17 : allCells dirCell 4889 4944 [
    10287459482542943953378407389815845107842950544629(*\ldots*)407766484658,
    18361862748823562884835754746870110948299090702878(*\ldots*)800100330722,
    38864081836796090953348552457253283813803933986586(*\ldots*)810783856434,
    147494596195641123668, 0, 0, 0, (*\ldots*)] = true := by
  decide +kernel
theorem dirCover : Cover dirCell 0 62720 := Cert.Dir000.cover.trans <| (*\ldots*)
theorem dir_cap_of (hC : Cover dirCell 0 62720) {a b c x : (*$\mathbb{R}$*)}
    (ha1 : 1 / 2 (*$\le$*) a) (ha2 : a (*$\le$*) 7 / 4) (hb1 : 1 / 8 (*$\le$*) b) (hb2 : b (*$\le$*) 1)
    (hc1 : 0 (*$\le$*) c) (hc2 : c (*$\le$*) 7 / 4) (hx1 : -1 (*$\le$*) x) (hx2 : x (*$\le$*) 1)
    (hG : (*$\forall$*) i, 0 (*$\le$*) (dirGR a b c x).get i) (hP : dirPR a b c x = 0) :
    (*$\exists$*) y0 y1, trSR (dirCertR a b c x) y0 y1 < 3 / 4
\end{lstlisting}

We cut the regions into grids of cells: $40\times28\times56$ cells of side $1/32$ in $\Odir$ and $8\times96\times96$ cells of side $1/64\times1/32\times1/32$ in $\Och$, 136{,}448 in all. \code{Grid.lean} numbers the cells and defines \code{dirCell n h} to run \code{checkBoxH} on cell \code{n} with hint \code{h} (Listing~\ref{lst:cert}). A certificate theorem states that \code{allCells dirCell n e} $[h_n,\dots,h_{e-1}]$ is \code{true}, i.e., that each cell from $n$ to $e-1$ passes with the hint listed for it. Such a theorem is proved by \code{decide +kernel}. Most hints are $0$, for a cell with an empty range of $x$ or one that can be dropped in one step; the longest hint has 11{,}114 bits. The lemma \code{allCells\_cover} turns a certificate into a statement \code{Cover dirCell n e}. The modules chain these statements with \code{Cover.trans}, and \code{Cover.lean} proves \code{dirCover : Cover dirCell 0 62720} and \code{chCover : Cover chCell 0 73728}. From these two theorems, \code{checkBoxH\_sound} and the carriers, \code{Cap.lean} then proves (C1) and (C2) for the real functions (\code{dir\_cap\_of} in Listing~\ref{lst:cert}).

We can predict the kernel time of a certificate from the replay. Our script \code{mkcert.py} estimates the time of a cell as $1.15\,(0.07\,n_P+0.21\,n_T)+0.005$ seconds, where $n_P$ and $n_T$ count the jet evaluations of $P$ and of $\tr S(y)$ in the replay of the cell, that is, about 80\,ms and 240\,ms per evaluation. Since the memory of a kernel check grows with the length of the theorem, the cells are cut into theorems of about 6\,s of predicted time. The theorems are then grouped into 211 modules: \code{Dir000} to \code{Dir158} with 2573 theorems and \code{Ch000} to \code{Ch051} with 845. Even so, a module can need up to about 5\,GB while the kernel checks its longest theorem, and the build script therefore checks the 52 modules with the longest theorems first, in batches. The model is optimistic: it predicted 20{,}209\,s of kernel time for all modules, and they took 27{,}483\,s, 36\% more.

\section{Results}\label{sec:results}

\subsection{Search and Kernel Check}

\begin{table}[t]
\caption{The Python search of~\cite{Tom26} and the generator \code{Gen}. Python counts are evaluations made during the search, on 32 processes. \code{Gen} counts are the jet evaluations in the replay of its final hints, which is the work the kernel repeats.}\label{tab:search}
\centering
\begin{tabular}{@{}lrrrr@{}}
\toprule
& \multicolumn{2}{c}{Python} & \multicolumn{2}{c}{\code{Gen}}\\
\cmidrule(lr){2-3}\cmidrule(l){4-5}
& direct & chart & direct & chart\\
\midrule
Initial cells & 62{,}720 & 73{,}728 & 62{,}720 & 73{,}728\\
Splits & 318{,}688 & 34 & 3{,}800 & 0\\
Leaves & 381{,}408 & 73{,}762 & 66{,}520 & 73{,}728\\
Evaluations of $P$ & 17{,}422{,}672 & 379{,}455 & 36{,}673 & 1{,}749\\
Evaluations of $\tr S(y)$ & 645{,}551 & 18{,}608 & 49{,}527 & 18{,}523\\
Wall-clock time & 245\,s & 10\,s & 10\,min & 6\,min\\
\bottomrule
\end{tabular}
\end{table}

In Table~\ref{tab:search} we compare the Python search with \code{Gen}. The kernel in fact replays only the final decisions of \code{Gen}. Over both regions, they amount to 38{,}422 jet evaluations of $P$ and 68{,}050 of $\tr S(y)$. The Python search, in contrast, evaluated $P$ 17.8~million times. A jet evaluation of $P$ costs about 80\,ms in the kernel.

Our machine has 32 hardware threads and 128\,GB of memory. On it, the 211 certificate modules took 19 minutes of wall-clock time and 7.6 hours of processor time. The slowest module took 173\,s. The rest of the development built in about two minutes.

\subsection{Independent Checks}\label{sec:indep}

The certificates rely on the kernel's GMP arithmetic far more than ordinary proofs do, so we checked the proof again with nanoda~0.4.17\footnote{\url{https://github.com/ammkrn/nanoda_lib}}, a type checker for Lean written in Rust that comes with the toolchain. Its arithmetic on natural numbers is done by the Rust library num-bigint, and it reads the exports of lean4export. The check has three parts. Part~A is the closure of the main theorems, with the 3418 certificates replaced by axioms of the same names and types. Part~B is the closure of each certificate module with the standard axioms; these closures are checked 16 modules at a time. We then use a script to check that each certificate assumed in part~A is proved exactly once in part~B with the same type, that the declarations shared by the two parts are identical, and that part~B uses no axiom except \code{propext}; declarations are compared by SHA-256 hashes of a canonical form. lean4export panics on a missing root name but still exits with status~0, so the script also searches the error output for a panic. The three parts took 29 minutes. As a control, a script adds 1 to a hint in the statement of a certificate in the export of its module, and nanoda rejects the result.

\subsection{A Version for a Public Registry}\label{sec:registry}

\begin{table}[t]
\caption{lean4body v1.0 and the version in the Palomar registry. Times are wall-clock unless marked. The local times of v1.0 are on a machine with 32 hardware threads, and those of the registry version on a machine with 384 hardware threads that was also running other jobs; the registry version was built with at most 16 modules at a time.}\label{tab:versions}
\centering
\begin{tabular}{@{}lrr@{}}
\toprule
& v1.0 & registry version\\
\midrule
Cells & 136{,}448 & 17{,}056\\
Certificate theorems / modules & 3418 / 211 & 7307 / 185\\
Build & 21\,min & 12\,min\\
Processor time & 7.6\,h (certificates) & 2.8\,h (all)\\
Memory of one module & about 5\,GB & about 3.4\,GB\\
nanoda, parts A, B and glue & 29\,min & 9\,min\\
Continuous integration, GitHub & -- & 28\,min\\
Registry verification job & -- & 2\,h 55\,min\\
\bottomrule
\end{tabular}
\end{table}

We also submitted a version of the proof to the Palomar registry\footnote{\url{https://palomar-registry.org}} of Lean formalizations. The registry requires the module system, at most 10{,}000 lines per file and the standard axioms, and its verification job must finish within 5.5 hours on 16 processors with 32\,GB of memory. The job runs \code{lake comparator}. This tool checks that the statements in a \code{Challenge} file, which imports only Mathlib and copies the definitions, are those proved in a \code{Solution} file; it also checks the axioms and replays the solution in Lean's kernel and in nanoda. The registry replays the solution in con-ron as well, another type checker written in Rust that has its own arithmetic. We first ported lean4body to the module system as it was, but unfortunately the port did not fit: \code{leanchecker} alone took 8 hours, and con-ron needed more than 32\,GB.

The registered version proves the same statements, Theorem~B with its constant $1/4$ included, and differs in the certificates. Cells twice as large in each direction turned out to be enough: $20\times14\times28$ and $4\times48\times48$ cells, 17{,}056 in all, one eighth of the original number. We wrote a new, greedy generator for the hints. It splits a box while a split along some dimension makes the box cheaper to check and localizes each leaf as cheaply as possible, and where the $y$ chosen at the centre of a box does not certify it, it searches the 32-bit fields of $y$ for one that does; this search lowered the predicted checking time by about 6\%. A further program, \code{Cut.lean}, replays the hints and predicts the kernel time of each box in each cell's tree of splits, in units of about 0.75\,s. A cell over 4 units is cut into the topmost boxes of the tree that take at most 4 units, and each of these boxes is certified by a separate theorem. We chose this fine cut because of con-ron, whose memory grows with the time it spends on a single declaration. The predicate \code{Certified} in \code{Grid.lean} accepts a box if \code{checkBoxH} succeeds on it or if both halves of one of its splits are certified, and \code{Cap.lean} proves it sound. The result has 7307 theorems in 185 modules: 1309 for runs of cells and 5998 for pieces. Of the leaf boxes without a split, 44 remain above 4 units, the largest at 4.90. Outside the certificates the changes are small: the grids change in the checker, \code{Cap.lean} assembles the new certificates, every file became a module, and four of the definitions that the \code{Challenge} file copies moved to \code{Defs.lean}. Table~\ref{tab:versions} compares the two versions.

The registry's run took 10{,}507\,s, 53\% of its limit, with 508\,s for the build, 17\,s for the export and 9{,}859\,s for the comparator. On a machine with 384 hardware threads, \code{leanchecker} replayed the registered solution in 17{,}752\,s, against 9{,}855\,s for the solution of version~1, whose certificates, made by an earlier version of the generator, prove only $\tr S(y)<127/128$; the two ran side by side on the two hardware threads of one core. Their ratio, 1.80, is close to that of the predicted units, 1.83 (20{,}472 against 11{,}193). The control also had to be changed. Since a hint is replayed rather than trusted, a changed hint can still certify the cell: in version~1, adding 1 to the hint of the control changed only the coordinate of the first split, Lean's kernel proved the changed statement true, and nanoda, correctly, accepted it. The control therefore adds the smallest number for which Lean's kernel proves the changed statement false.

The repository of the registered version also runs a continuous-integration check on GitHub-hosted runners with four processors and 16\,GB of memory. It runs three checkers: Lean's kernel (through \code{leanchecker}), nanoda and con-ron. Con-ron declines any axiom other than the standard three, so it checks the certificates, in part~B, and not part~A. The run took 28 minutes: 20 shards of 10 to 17 minutes each, then a final job of 11 minutes. Summed over the shards, the builds took 1.1\,h, \code{leanchecker} 3.0\,h, nanoda 2.8\,h and con-ron 7.1\,h, and no single process in a shard needed more than 0.9\,GB of memory. All three checkers rejected the control.

\section{Related Work}\label{sec:related}

Proof by reflection~\cite{Bou97} and the use of a proof assistant's evaluator for large computations~\cite{BB02} are old ideas. Coq gained a compiled implementation of strong reduction~\cite{GL02} and later machine integers and arrays in the kernel~\cite{AGST10}, which made checks such as those of SAT certificates feasible. Lean~4~\cite{dMU21} accelerates arithmetic on \code{Nat} literals in the kernel, and we rely on this for all of our certificates. The tactic \code{native\_decide} offers compiled evaluation at the price of an axiom.

In Coq, Melquiond's methods for proving bounds by computation~\cite{Mel08} grew into CoqInterval~\cite{MDM16}; in HOL Light, Solovyev and Hales verified nonlinear inequalities with Taylor intervals~\cite{SH13}. Their tool served the formal proof of the Kepler conjecture~\cite{HAB+17}, whose more than 23{,}000 inequalities took approximately 5000 processor-hours on Microsoft Azure. Our computation is much smaller, but the claims it has to establish are not inequalities over boxes (Sect.~\ref{sec:awkward}). This is why our hints record localizations and witnesses as well as splits. Tucker's computer-assisted proof that the Lorenz attractor exists~\cite{Tuc02} was checked again by Immler in Isabelle/HOL with a verified ODE solver~\cite{Imm15,Imm18}, and Gonthier's proof of the four colour theorem~\cite{Gon08} runs its case analysis inside Coq.

The closest relatives in method are the SAT-based proofs in combinatorics, of the Boolean Pythagorean triples theorem~\cite{HKM16}, of Keller's conjecture~\cite{BHMN20} and of the empty hexagon theorem~\cite{HS24}. Their proofs were checked by proof checkers, including the verified checker cake\_lpr~\cite{THM21}, and the encoding of the empty hexagon theorem was verified in Lean~\cite{SNG+24}. There, as here, an untrusted search produces a certificate that a small trusted checker validates. In our case the checker is a Lean function, and the kernel is the only trusted executor. Lean4Lean~\cite{Car24} is a step toward a verified type checker for Lean, which would shrink the trusted base of results like ours further.

\section{Discussion}\label{sec:discussion}

The constant $1/4$ in Theorem~B can be traded for checking time. The covering argument needs only some $\varepsilon>0$, and on the same grids, with the same generator, certificates of $\tr S(y)<127/128$ have about half the predicted checking time of those of $\tr S(y)<3/4$ (11{,}193 units against 22{,}025). Version~1 of the same registry entry uses this threshold and proves only $Q\ge K/128$.

Parts of the checker can be reused. \code{Num.lean} and \code{Jet.lean} do not depend on the problem, and the carrier construction applies to any formula written against \code{Ops}. \code{Search.lean}, however, is specific to claims of the form (C1), with one variable localized and a witness chosen for each localized piece.

The design also has limits. The fixed-point format has absolute precision $2^{-96}$ and no values below $-2^{160}$. This is enough here because all quantities are bounded on the two regions, but the format would not be suitable for quantities of widely different sizes. Since the jets are first order, the boxes have to be small; higher-order Taylor models~\cite{SH13} might reduce the number of cells. Finally, Lean's kernel trusts GMP for the arithmetic of the certificates. This is one reason why we also ran nanoda and con-ron: nanoda does its arithmetic with num-bigint, and con-ron with its own code.

\begin{credits}
\subsubsection{Use of Generative AI.}
Generative AI tools were used to find and prove the results of~\cite{Tom26}, to write the programs and the Lean formalization described in this paper, and to draft the text.

\subsubsection{\discintname}
The author has no competing interests to declare that are relevant to the content of this article.
\end{credits}

\section*{Data-Availability Statement}

The formalization, the Python programs, and the logs of the build and the nanoda check are in the repository \url{https://github.com/tejasvi/lean4body}, release v1.0 (commit \code{58ac0fbf4d66}), archived on Zenodo~\cite{Tom26z}. The registry version is \url{https://github.com/tejasvi/lean4body-uniqueness} at commit \code{353653eceec4}, version~2 of entry PALOMAR-2026-10-02-000011 of the Palomar registry (\url{https://palomar-registry.org/entry.html?id=PALOMAR-2026-10-02-000011&version=2}), verified in run 37187361609 of \url{https://github.com/PalomarRegistry/PalomarSubmission}; its continuous-integration run is run 37182678749 of \url{https://github.com/tejasvi/lean4body-uniqueness}. Version~1 of the entry, at commit \code{9da70faf3423}, uses the threshold $127/128$. The toolchain is Lean v4.35.0-rc3 with Mathlib at the matching tag; the Python programs use Python 3.14.4 and python-flint 0.9.0.

\begin{thebibliography}{10}

\bibitem{Alb96}
Albouy, A.: The symmetric central configurations of four equal masses. In: Hamiltonian Dynamics and Celestial Mechanics. Contemp. Math., vol.~198, pp. 131--135. AMS (1996)

\bibitem{ACS12}
Albouy, A., Cabral, H.E., Santos, A.A.: Some problems on the classical $n$-body problem. Celest. Mech. Dyn. Astron. \textbf{113}, 369--375 (2012)

\bibitem{AFS08}
Albouy, A., Fu, Y., Sun, S.: Symmetry of planar four-body convex central configurations. Proc. R. Soc. A \textbf{464}, 1355--1365 (2008)

\bibitem{AGST10}
Armand, M., Gr\'egoire, B., Spiwack, A., Th\'ery, L.: Extending Coq with imperative features and its application to SAT verification. In: Interactive Theorem Proving (ITP 2010). LNCS, vol.~6172, pp. 83--98. Springer (2010)

\bibitem{BB02}
Barendregt, H., Barendsen, E.: Autarkic computations in formal proofs. J. Autom. Reasoning \textbf{28}(3), 321--336 (2002)

\bibitem{Bou97}
Boutin, S.: Using reflection to build efficient and certified decision procedures. In: Theoretical Aspects of Computer Software (TACS 1997). LNCS, vol.~1281, pp. 515--529. Springer (1997)

\bibitem{BHMN20}
Brakensiek, J., Heule, M., Mackey, J., Narv\'aez, D.: The resolution of Keller's conjecture. In: Automated Reasoning (IJCAR 2020). LNCS, vol.~12166, pp. 48--65. Springer (2020)

\bibitem{Car24}
Carneiro, M.: Lean4Lean: Towards a verified typechecker for Lean, in Lean. arXiv:2403.14064 (2024)

\bibitem{CCR19}
Corbera, M., Cors, J.M., Roberts, G.E.: Classifying four-body convex central configurations. Celest. Mech. Dyn. Astron. \textbf{131}, Paper~34 (2019)

\bibitem{dMU21}
de~Moura, L., Ullrich, S.: The Lean 4 theorem prover and programming language. In: Automated Deduction (CADE 28). LNCS, vol.~12699, pp. 625--635. Springer (2021)

\bibitem{Dzi00}
Dziobek, O.: \"Uber einen merkw\"urdigen Fall des Vielk\"orperproblems. Astron. Nachr. \textbf{152}, 33--46 (1900)

\bibitem{Gon08}
Gonthier, G.: Formal proof---the four-color theorem. Notices Amer. Math. Soc. \textbf{55}(11), 1382--1393 (2008)

\bibitem{GL02}
Gr\'egoire, B., Leroy, X.: A compiled implementation of strong reduction. In: International Conference on Functional Programming (ICFP 2002), pp. 235--246. ACM (2002)

\bibitem{HAB+17}
Hales, T., Adams, M., Bauer, G., Dang, T.D., Harrison, J., Hoang, L.T., Kaliszyk, C., Magron, V., McLaughlin, S., Nguyen, T.T., Nguyen, Q.T., Nipkow, T., Obua, S., Pleso, J., Rute, J., Solovyev, A., Ta, T.H.A., Tran, N.T., Trieu, T.D., Urban, J., Vu, K., Zumkeller, R.: A formal proof of the Kepler conjecture. Forum Math. Pi \textbf{5}, e2 (2017)

\bibitem{HKM16}
Heule, M.J.H., Kullmann, O., Marek, V.W.: Solving and verifying the Boolean Pythagorean triples problem via cube-and-conquer. In: Theory and Applications of Satisfiability Testing (SAT 2016). LNCS, vol.~9710, pp. 228--245. Springer (2016)

\bibitem{HS24}
Heule, M.J.H., Scheucher, M.: Happy ending: An empty hexagon in every set of 30 points. In: Tools and Algorithms for the Construction and Analysis of Systems (TACAS 2024). LNCS, vol.~14570, pp. 61--80. Springer (2024)

\bibitem{Imm15}
Immler, F.: Verified reachability analysis of continuous systems. In: Tools and Algorithms for the Construction and Analysis of Systems (TACAS 2015). LNCS, vol.~9035, pp. 37--51. Springer (2015)

\bibitem{Imm18}
Immler, F.: A verified ODE solver and the Lorenz attractor. J. Autom. Reasoning \textbf{61}, 73--111 (2018)

\bibitem{Joh17}
Johansson, F.: Arb: efficient arbitrary-precision midpoint-radius interval arithmetic. IEEE Trans. Comput. \textbf{66}(8), 1281--1292 (2017)

\bibitem{KN85}
Krawczyk, R., Neumaier, A.: Interval slopes for rational functions and associated centered forms. SIAM J. Numer. Anal. \textbf{22}(3), 604--616 (1985)

\bibitem{Lea03}
Leandro, E.S.G.: Finiteness and bifurcations of some symmetrical classes of central configurations. Arch. Ration. Mech. Anal. \textbf{167}, 147--177 (2003)

\bibitem{LS02}
Long, Y., Sun, S.: Four-body central configurations with some equal masses. Arch. Ration. Mech. Anal. \textbf{162}, 25--44 (2002)

\bibitem{MB32}
MacMillan, W.D., Bartky, W.: Permanent configurations in the problem of four bodies. Trans. Amer. Math. Soc. \textbf{34}, 838--875 (1932)

\bibitem{MDM16}
Martin-Dorel, \'E., Melquiond, G.: Proving tight bounds on univariate expressions with elementary functions in Coq. J. Autom. Reasoning \textbf{57}(3), 187--217 (2016)

\bibitem{Mel08}
Melquiond, G.: Proving bounds on real-valued functions with computations. In: Automated Reasoning (IJCAR 2008). LNCS, vol.~5195, pp. 2--17. Springer (2008)

\bibitem{Moe15}
Moeckel, R.: Central configurations. In: Llibre, J., Moeckel, R., Sim\'o, C.: Central Configurations, Periodic Orbits, and Hamiltonian Systems, pp. 105--167. Advanced Courses in Mathematics CRM Barcelona, Birkh\"auser (2015)

\bibitem{MKC09}
Moore, R.E., Kearfott, R.B., Cloud, M.J.: Introduction to Interval Analysis. SIAM (2009)

\bibitem{PCS07}
P\'erez-Chavela, E., Santoprete, M.: Convex four-body central configurations with some equal masses. Arch. Ration. Mech. Anal. \textbf{185}, 481--494 (2007)

\bibitem{Shu71}
Shub, M.: Appendix to Smale's paper: Diagonals and relative equilibria. In: Manifolds -- Amsterdam 1970. Lecture Notes in Math., vol.~197, pp. 199--201. Springer (1971)

\bibitem{SH13}
Solovyev, A., Hales, T.C.: Formal verification of nonlinear inequalities with Taylor interval approximations. In: NASA Formal Methods (NFM 2013). LNCS, vol.~7871, pp. 383--397. Springer (2013)

\bibitem{SNG+24}
Subercaseaux, B., Nawrocki, W., Gallicchio, J., Codel, C., Carneiro, M., Heule, M.J.H.: Formal verification of the empty hexagon number. In: Interactive Theorem Proving (ITP 2024). LIPIcs, vol.~309, pp. 35:1--35:19 (2024)

\bibitem{SXY23}
Sun, S., Xie, Z., You, P.: On the uniqueness of convex central configurations in the planar 4-body problem. Regul. Chaotic Dyn. \textbf{28}, 512--532 (2023)

\bibitem{THM21}
Tan, Y.K., Heule, M.J.H., Myreen, M.O.: cake\_lpr: Verified propagation redundancy checking in CakeML. In: Tools and Algorithms for the Construction and Analysis of Systems (TACAS 2021). LNCS, vol.~12652, pp. 223--241. Springer (2021)

\bibitem{mC20}
The mathlib Community: The Lean mathematical library. In: Certified Programs and Proofs (CPP 2020), pp. 367--381. ACM (2020)

\bibitem{Tom26}
Tomar, T.S.: Uniqueness of four-body convex central configurations for all masses. arXiv:2609.35632 (2026)

\bibitem{Tom26z}
Tomar, T.S.: lean4body: programs, certificates and Lean formalization for the uniqueness of convex four-body central configurations, version 1.0. Zenodo (2026). \doi{10.5281/zenodo.23004297}

\bibitem{Tuc02}
Tucker, W.: A rigorous ODE solver and Smale's 14th problem. Found. Comput. Math. \textbf{2}(1), 53--117 (2002)

\bibitem{Xia04}
Xia, Z.: Convex central configurations for the $n$-body problem. J. Differential Equations \textbf{200}, 185--190 (2004)

\end{thebibliography}

\end{document}