\section{Lithium integral certificates}\label{sec:11-lithium-integrals}

This appendix defines the trial functions and sources used in
Section~\ref{sec:06-lithium-binding} (Definition~\ref{def:li-trials}). It
then proves the rational enclosures of their moments
(Lemma~\ref{lem:li-physical-moments}).

\begin{definition}[Fixed lithium trials and sources]\label{def:li-trials}
For $x_i\in\mathbb R^3$, $1\le i\le3$, put $r_i=|x_i|$,
$d=(r_1,r_2,r_3,|x_1-x_2|,|x_1-x_3|,|x_2-x_3|)$,
and let $A_{13}=1-P_{13}$, where $P_{13}$ exchanges positions
$1$ and $3$. For $a\in(0,\infty)^3$ and $\nu\in\mathbb Z_{\ge0}^6$
put $F_{a,\nu}=e^{-\sum_i a_i r_i}\prod_{k=1}^6d_k^{\nu_k}$.
Let $(\nu,c_\nu)$ be the 100 exponent/integer pairs in
Table~\ref{tab:li-s-trial}, and let $(b_t)_{t=1}^{140}$ be the ordered
integers in Table~\ref{tab:li-p-trial}. Write $x_{i,3}$ for the
third Cartesian coordinate of electron $i$. Define the ordered
mode triples $(a_t,\nu_t,i_t)$, $1\le t\le140$, as follows.
The first 112 modes have $i_t=3$: for each rate triple
$(14/5,14/5,13/20)$, $(14/5,14/5,6/5)$,
$(31/10,12/5,4/5)$, $(12/5,31/10,4/5)$, in that order,
use the factors $1,d_1,\ldots,d_6,d_jd_k$
($1\le j\le k\le6$, lexicographically).
The last 28 modes use rates $(14/5,14/5,3/4)$ and
$(14/5,14/5,5/4)$, in order; for each use $i_t=1$ then $2$,
each with factors $1,d_1,\ldots,d_6$.
The multi-index $\nu_t$ is the exponent vector of the indicated factor.
Define
\[
 u_S=10^{-14}\sum_\nu c_\nu A_{13}F_{(14/5,14/5,17/20),\nu},
 \qquad
 u_P=10^{-18}\sum_{t=1}^{140}b_t A_{13}
                    [x_{i_t,3}F_{a_t,\nu_t}].
\]
The first sum ranges over the 100 pairs of Table~\ref{tab:li-s-trial},
read left to right.

Set $Q=8/3$. For $2\le n\le6$ put
\[
 s_n=A_{13}\left[
 e^{-Qr_1-Qr_2-Qr_3/n}
 \sum_{j=0}^{n-1}\frac{(-2Q/n)^j}{j!}
                 \binom{n}{n-1-j}r_3^j\right],
 \quad
 C_{S,n}=\frac{Q^3}{\pi}\frac{(2Q/n)^{3/2}}{4n\sqrt\pi}.
\]
Put
\[
 t_2=A_{13}[x_{3,3}e^{-Qr_1-Qr_2-Qr_3/2}],\qquad
 t_3=A_{13}[x_{3,3}(1-Qr_3/6)e^{-Qr_1-Qr_2-Qr_3/3}],
\]
\[
 C_{P,2}=\frac{Q^3}{\pi}\sqrt{\frac{(Q/2)^5}{2\pi}},\qquad
 C_{P,3}=\frac{Q^3}{\pi}\sqrt{\frac{4(Q/3)^5}{3\pi}},
\]
Define $p_{S,n}=C_{S,n}s_n$ for $2\le n\le6$ and
$p_{P,n}=C_{P,n}t_n$ for $n\in\{2,3\}$.
All square roots are positive. Spatial norms use ordinary Lebesgue
measure. Each trial and source $f$, which satisfies $P_{13}f=-f$,
is represented isometrically
as a physical spinor with $\uparrow\downarrow\uparrow$ component
$f/\sqrt3$, extended by signed simultaneous permutations to the
other two strings with two up spins, and zero on the remaining
spin strings. Use the same symbols for these spinors.
The labels $S$ and $P$ here distinguish the two trial families: their
orbital angular momenta are $0$ and $1$, respectively, their orbital
magnetic number is zero, their parities are even and odd, respectively,
and their spin magnetic number is $1/2$. No total-spin restriction is imposed.
\end{definition}

\begingroup
\scriptsize
\setlength{\tabcolsep}{4pt}
\begin{longtable}{@{}cr@{\quad}cr@{\quad}cr@{\quad}cr@{}}
\caption{Lithium $S$ trial: ordered pairs $(\nu,c_\nu)$, read left to right. The six digits of $\nu$ give the exponents of $d_1,\ldots,d_6$. The coefficient of $A_{13}F_{(14/5,14/5,17/20),\nu}$ is $10^{-14}c_\nu$.}\label{tab:li-s-trial}\\
\hline
$\nu$ & $c_\nu$ & $\nu$ & $c_\nu$ & $\nu$ & $c_\nu$ & $\nu$ & $c_\nu$\\
\hline
\endfirsthead
\caption[]{Lithium $S$ trial: ordered pairs $(\nu,c_\nu)$, read left to right. The six digits of $\nu$ give the exponents of $d_1,\ldots,d_6$. The coefficient of $A_{13}F_{(14/5,14/5,17/20),\nu}$ is $10^{-14}c_\nu$. (continued)}\\
\hline
$\nu$ & $c_\nu$ & $\nu$ & $c_\nu$ & $\nu$ & $c_\nu$ & $\nu$ & $c_\nu$\\
\hline
\endhead
\hline
\endfoot
\hline
\endlastfoot
\texttt{000000} & 192578507642049 & \texttt{100000} & -18637105584193 & \texttt{010000} & -4893191508198 & \texttt{001000} & -377772895588427\\
\texttt{000100} & 21687913877446 & \texttt{101000} & 55883362330885 & \texttt{011000} & 64991563194886 & \texttt{002000} & -38457760107399\\
\texttt{001100} & -97869737076541 & \texttt{000200} & -2817165576809 & \texttt{000110} & -6574049558047 & \texttt{000011} & -21175160412110\\
\texttt{210000} & 128411633891504 & \texttt{201000} & -111268970809246 & \texttt{200100} & 85958547690281 & \texttt{200001} & -11135646743619\\
\texttt{120000} & -33435030557499 & \texttt{111000} & 99213941867312 & \texttt{102000} & 22771619995637 & \texttt{101100} & -40915170411464\\
\texttt{101010} & -33137323291035 & \texttt{100110} & -15706633630960 & \texttt{100020} & 13271380434227 & \texttt{100011} & 8708426264134\\
\texttt{021000} & -33637679816705 & \texttt{020001} & -96379393879640 & \texttt{010101} & -35130689066215 & \texttt{010002} & 5752002624058\\
\texttt{002100} & -14573043963146 & \texttt{001200} & 86372395390789 & \texttt{001101} & -39908942457146 & \texttt{000300} & -41645540686493\\
\texttt{000201} & 36940233903379 & \texttt{000102} & 7502802960424 & \texttt{000021} & 3160700248761 & \texttt{300100} & -63250958657558\\
\texttt{220000} & -41890829949681 & \texttt{210100} & -114611305957702 & \texttt{201100} & 33697170918180 & \texttt{201010} & -17594454979840\\
\texttt{200110} & 96496437822204 & \texttt{200020} & -6023803910370 & \texttt{130000} & -49526101674321 & \texttt{121000} & 109082871822618\\
\texttt{120100} & -97003483595728 & \texttt{112000} & 6662575813563 & \texttt{111100} & -42994396428613 & \texttt{111010} & 5261593757741\\
\texttt{110200} & 102382385187551 & \texttt{100210} & -61857093841605 & \texttt{100003} & -421408036738 & \texttt{030001} & 16481601753808\\
\texttt{022000} & -23563953921328 & \texttt{020101} & 121205576172158 & \texttt{011002} & 1806540114768 & \texttt{010201} & -65251992063050\\
\texttt{010003} & -1563301644880 & \texttt{002200} & 14483571075089 & \texttt{001300} & -49994233247290 & \texttt{001102} & 4084328942332\\
\texttt{000400} & 34836644219690 & \texttt{000310} & 8704218112017 & \texttt{000301} & -7540626026929 & \texttt{000220} & 1416699560466\\
\texttt{000040} & -144518102488 & \texttt{401000} & 5004468341455 & \texttt{320000} & -26938458762230 & \texttt{311000} & -16259637468816\\
\texttt{310100} & 76817730215707 & \texttt{301100} & -16853901957893 & \texttt{230000} & 34638676982410 & \texttt{211010} & 12384132460967\\
\texttt{210110} & -35174471648651 & \texttt{202100} & 5279457589922 & \texttt{200300} & 4306233837560 & \texttt{200210} & -20559937525297\\
\texttt{200030} & 684896568552 & \texttt{131000} & -26028992998374 & \texttt{130100} & 53347598075651 & \texttt{121100} & -24773719544820\\
\texttt{120200} & 10934434016173 & \texttt{111200} & 26102826113440 & \texttt{111110} & -11143103209542 & \texttt{110300} & -52094422510428\\
\texttt{110210} & 24973988134754 & \texttt{100310} & 21086139523661 & \texttt{050000} & 5869589803234 & \texttt{030110} & -2399864140321\\
\texttt{030101} & -24678742221128 & \texttt{023000} & 718307824117 & \texttt{022100} & 7839439226000 & \texttt{020201} & -17103988441973\\
\texttt{013100} & -541959493170 & \texttt{010301} & 22558696338943 & \texttt{002300} & -3794756813227 & \texttt{002201} & -388502500303\\
\texttt{001400} & 9048105445461 & \texttt{000500} & -6805192162237 & \texttt{000410} & -4526294280747 & \texttt{000104} & -170281908585\\
\end{longtable}
\endgroup

\begingroup
\scriptsize
\setlength{\tabcolsep}{4pt}
\begin{longtable}{@{}rrrrr@{}}
\caption{Lithium $P$ trial: successive integers $b_t$, read left to right in the mode order of Definition~\ref{def:li-trials}. The coefficient of mode $t$ is $10^{-18}b_t$.}\label{tab:li-p-trial}\\
\hline
$t$ & $b_t$ & $b_{t+1}$ & $b_{t+2}$ & $b_{t+3}$\\
\hline
\endfirsthead
\caption[]{Lithium $P$ trial: successive integers $b_t$, read left to right in the mode order of Definition~\ref{def:li-trials}. The coefficient of mode $t$ is $10^{-18}b_t$. (continued)}\\
\hline
$t$ & $b_t$ & $b_{t+1}$ & $b_{t+2}$ & $b_{t+3}$\\
\hline
\endhead
\hline
\endfoot
\hline
\endlastfoot
1 & -120015870113720000 & -5527719335490280000 & -5662662512097800000 & 186960680244002000\\
5 & 8993657481192360000 & 98711850003572900 & 194083685119928000 & 232359157147092000\\
9 & 887314913863710000 & 467761626184719000 & 196480110615253000 & -172039858373383000\\
13 & -12662948266792200 & -179329958693696000 & 169456497556216000 & 158124641883264000\\
17 & 106771804491281000 & 16061905386670600 & -98175363322029100 & -545263189686892000\\
21 & 13237502766589500 & 51642919598783700 & -627437759755436000 & 4512159831252850\\
25 & 37625439684789800 & -32487525562657900 & 49351572655156400 & -49796552667561200\\
29 & -229036377859776000 & -1439576674902220000 & -8135482863011680000 & 4503775487779260000\\
33 & 7517381434037780000 & -14823664189105800 & -621137013893765000 & 310741670530606000\\
37 & 2255510321923590000 & -825346181192189000 & -365230774580283000 & 708810332821348000\\
41 & -652864814227281000 & -24564618439033400 & -3057666884818770000 & 990392680619719000\\
45 & -492573803101796000 & 956933229999546000 & 1929900881087360000 & 2796120011114770000\\
49 & -249434812347310000 & -226786886787875000 & -1491807709059100000 & 206995224193631000\\
53 & 189556961807521000 & -351229476439412000 & 316334819589271000 & -453976751964670000\\
57 & 3867279158789010000 & 1370756521039610000 & 9641774677603950000 & -3418116940985900000\\
61 & -8762790980746690000 & 355658931719882000 & -112416908162737000 & -187364998886040000\\
65 & 307814858046278000 & 100137319601905000 & -652904699417170000 & -38550485008248200\\
69 & -77087997812530800 & -1920733570648110000 & 109637964067611000 & 1236605291154880000\\
73 & 74764536721586000 & 94990533335603500 & 524887971724604000 & 115363373071532000\\
77 & -224987633941586000 & -64168674989026800 & 508559138614380000 & -39733921285032700\\
81 & 7788471490130710 & 88814579642896900 & -84895271159553000 & 48270722631235900\\
85 & 124474002247613000 & 6232093768174070000 & 2375559818539930000 & -2510325407821860000\\
89 & -6783982643852250000 & 371120307392657000 & 1354444402963760000 & -1048106464222570000\\
93 & -477595009973174000 & 29630064983615600 & 808180795456249000 & -40943101887638200\\
97 & -125834307341012000 & -47601453984732800 & -120020287293432000 & -439909122232083000\\
101 & 88065815728008700 & 59447311705798700 & 411993905925623000 & 17219661067493700\\
105 & 29577249255805000 & -339613999279857000 & 509378795942874000 & 10691020494274200\\
109 & -88322402063989500 & 37903998405299800 & -62188273102382800 & 127034235858460000\\
113 & -351015977798309000 & -85738453397099300 & -83769320762075600 & -103928148905266000\\
117 & 79327390110805700 & -554173931860421 & 156401268414680000 & 42381359401224900\\
121 & 13903846915434600 & 28602294905942200 & 173299294555551000 & 48493875456655500\\
125 & -62897580204762500 & -113554182280085000 & -140903867796034000 & 313730933322561000\\
129 & 470068723557144000 & 1090711016080320000 & -416801262744941000 & 53182997120479400\\
133 & -836467615918667000 & 216214282604895000 & -38009874696524800 & -143874674567923000\\
137 & -709080785460829000 & -216177827252217000 & 195116461020833000 & 381035505002165000\\
\end{longtable}
\endgroup

\subsection{Physical expansion and arithmetic conventions}

The complete integer arrays, rational group inputs, returned integral
intervals and independent readers used below are supplied in the lithium
part of the accompanying code/data release. We specify their mathematical
contents and the predicates checked by the readers here. In particular,
the trial coefficients in Tables~\ref{tab:li-s-trial}
and~\ref{tab:li-p-trial} are interpreted as rationals.

\begin{setup}[Integral executions]\label{setup:lithium-integrals-execution}
We specialize Set-up~\ref{setup:execution} as follows. The exact computations
use unbounded integers, reduced rational arithmetic, integer factorials,
binomial coefficients, floor division and integer square roots. Parsing,
ordered containers, finite iteration, compression, file operations and
preservation of the recorded inputs, sources, outputs and executions are
assumed faithful; assertions are enabled. Caches retain results for the same
arguments and arithmetic environment.

The native radial executions use the specified C++ operation sequence,
signed 32-bit C and interface integers, binary64 doubles, matching argument
and return interfaces, and writable result buffers. Compilation, loading
and dispatch are faithful. Elementary binary64 operations and rational
conversions round to nearest with even ties, with gradual underflow and
flush-to-zero and denormals-are-zero disabled. There is no reassociation,
fused contraction, disruptive excess precision or concurrent change of
the floating environment. Adjacent-number operations, binary decomposition
and scaling by powers of two have their specified exact meanings.

The double-logarithmic calculation uses the same binary64 semantics for
each NumPy component operation, with faithful shapes, indices and
broadcasting. A reduction of $m$ nonnegative entries performs at most
$m-1$ additions, in any order. Any participating threads preserve these
operations and data ownership. The retained environment uses CPython
3.12.3 and GCC 13.3.0; native compilation uses
\texttt{-O2 -fno-fast-math -ffp-contract=off}.
The supplement records the remaining environment and source versions.
These are execution assumptions. Integral accuracy, admissibility of
arguments, coverage and normalization are proved below.
\end{setup}

For intervals we use endpointwise addition and the four-corner product
\[
 [a,b][c,d]=[\min\{ac,ad,bc,bd\},\max\{ac,ad,bc,bd\}].
\]
This convention retains the signs of all coefficients. A rational string
denotes its exact rational value. A numeric token in a saved binary64
interval denotes the exact dyadic value obtained by nearest-even parsing.
It need not equal the rational decimal written in that token.

\begin{definition}[Scalar integral arrays]\label{def:lithium-integrals-arrays}
For positive integers $R_i$, integers $n_i,p_{ij}$ and
$k=(R_1,R_2,R_3,n_1,n_2,n_3,p_{12},p_{13},p_{23})$, put
\begin{equation}\label{eq:lithium-integrals-master}
 \mathcal I(k)=(4\pi)^{-3}\int_{(\mathbb R^3)^3}
 e^{-\sum_iR_ir_i/180}\prod_i r_i^{n_i}
 \prod_{i<j}d_{ij}^{p_{ij}}\,dX.
\end{equation}
A coefficient dictionary is a finite rational sum of the scalar
primitives $F_{a,\nu}$ in Definition~\ref{def:li-trials}, or of
$x_{i,3}F_{a,\nu}$. After formal differentiation, we allow integer
distance exponents $\nu\in\mathbb Z^6$, interpreting the same primitive
formula away from collisions. The denominator of a dictionary is the least
common multiple of
its reduced coefficient denominators; the zero dictionary has denominator
one. Let $U_T$ be the dictionary of $u_T$ and $V_T$ the dictionary
obtained by applying the formal differential expression
$H_3(3)+15/2$ to $u_T$. Let $S_n,T_n$ be those of $s_n,t_n$.
The ordered base pairs are
\[
\begin{split}
 S:&\ (U_S,V_S),(V_S,V_S),(S_2,U_S),(S_2,V_S),
                 (S_3,V_S),\ldots,(S_6,V_S),\\
 P:&\ (U_P,V_P),(V_P,V_P),(T_2,U_P),(T_2,V_P),(T_3,V_P).
\end{split}
\]
The extra pairs are $(U_S,U_S),(U_S,S_3),\ldots,(U_S,S_6)$ and
$(U_P,U_P),(U_P,T_3)$.

For a pair $(f,g)$, multiply its dictionaries, adding rates and powers.
In a vector product replace $x_{i,3}x_{j,3}$ by
$(r_i^2+r_j^2-d_{ij}^2)/6$, with $d_{ii}=0$.
Set $D_c=\operatorname{den}(f)\operatorname{den}(g)$ in $S$ and
$D_c=6\operatorname{den}(f)\operatorname{den}(g)$ in $P$.
Multiply the coefficients by $D_c$ and the rates by $180$.
Replace each key by the lexicographically least of its six simultaneous
electron permutations, and add its signed integer weights. Delete a key
only when its entire column vector is zero. Denote the resulting base
key sets by $\mathcal L_T$, their weights by $W_c(k)$, and the extra
sets by $\mathcal L_T^+$. Missing column entries have weight zero.
The source factor $C_c$ is one for a norm or a trial--residual product,
and is the positive $C_{T,n}$ of Definition~\ref{def:li-trials} for a
column involving that source.
\end{definition}

\subsection{Integral formulas and error bounds}

Write $P_\ell$ for the Legendre polynomial normalized by
$P_\ell(1)=1$. For $p\ge-2$ and $0<t<1$ define
\[
 A_\ell^{(p)}(t)=\frac12\int_{-1}^1
       (1+t^2-2tx)^{p/2}P_\ell(x)\,dx,
 \qquad B_\ell(t)=t^{-\ell}A_\ell^{(-2)}(t).
\]
For a radial order $\sigma=(h,m,b)$ the pair ratios are
$t_{hm}=s,t_{mb}=v,t_{hb}=sv$ under
$r_h=R,r_m=Rs,r_b=Rsv$. Put
\[
 \kappa=9+\sum_i n_i+\sum_{i<j}p_{ij},\qquad
 C_\sigma=n_b+2,\qquad D_\sigma=n_m+n_b+p_{mb}+5.
\]
These letters denote scalar degrees, not source constants or column
denominators. All ordered integrals below have positive rational rates.
The conditions on the clusters that we use are
\begin{equation}\label{eq:lithium-integrals-clusters}
 9+\sum_i n_i+\sum_{i<j}p_{ij}>0,\qquad
 n_i+n_j+p_{ij}>-6\ (i<j),\qquad \sum_{i<j}p_{ij}>-6.
\end{equation}
For $j\ge0$ write $\mathsf h_j=\sum_{i=1}^j1/i$, with
$\mathsf h_0=0$.

\begin{lemma}[Angular reduction]\label{lem:lithium-integrals-angular}
Suppose $R_i>0$ are integers, $n_i,p_{ij}\ge-2$ are integers and
\eqref{eq:lithium-integrals-clusters} holds. Suppose further that there
are at most two pair powers $-2$, and that if there are two then the
remaining pair power is nonnegative and even. The master integral is
finite.
Its normalized angular average at distinct radii is
\begin{equation}\label{eq:lithium-integrals-angular-average}
 \prod_{i<j}r_{\mathrm{large}(ij)}^{p_{ij}}
 \sum_{\ell\ge0}(2\ell+1)\prod_{i<j}A_\ell^{(p_{ij})}(t_{ij}).
\end{equation}
If $p\ge-1$, $A_\ell^{(p)}$ is a rational polynomial with no
powers below $\ell$, and it vanishes for $2\ell>p$ when $p$ is even.
For $p=2m-1$, $m\ge0$, it equals
$\sum_{j=0}^m c_jt^{\ell+2j}$, where
\begin{equation}\label{eq:lithium-integrals-odd-coefficients}
 c_0=\frac{(-1)^m(2m-1)!!}{\prod_{j=0}^m(2\ell+1-2j)},\qquad
 c_{j+1}=-c_j\frac{(m-j)(2\ell+1-2m+2j)}{(j+1)(2\ell+3+2j)}.
\end{equation}
Here $(-1)!!=1$. For the inverse-square coefficient,
\begin{equation}\label{eq:lithium-integrals-inverse-square}
 B_\ell(t)=\sum_{j\ge0}
 \frac{\prod_{i=1}^{\ell}(2j+2i)}{\prod_{i=0}^{\ell}(2j+2i+1)}t^{2j},
 \qquad 0<B_\ell(t)\le1-\log(1-t).
\end{equation}
An even pair power makes \eqref{eq:lithium-integrals-angular-average}
terminate at the smallest half of an even pair power.
\end{lemma}

\begin{proof}
Rodrigues' formula gives
\[
 \int_{S^2}P_\ell(\omega\cdot u)P_j(\omega\cdot v)
             \frac{d\omega}{4\pi}
 =\frac{\mathbf1_{\ell=j}}{2\ell+1}P_\ell(u\cdot v).
\]
For example, the rotation-invariant reproducing kernel of degree $\ell$
spherical harmonics is $(2\ell+1)P_\ell$. Its trace determines this
coefficient. Applying the identity twice gives
\eqref{eq:lithium-integrals-angular-average}.
To compute a coefficient, expand $P_\ell(x)=\sum_{j=0}^\ell b_jx^j$
and substitute $y=1+t^2-2tx$. For $p\ge-1$ the result is
\begin{equation}\label{eq:lithium-integrals-angular-polynomial}
 \sum_{j=0}^\ell b_j\sum_{i=0}^j
 \frac{(-1)^i\binom ji(1+t^2)^{j-i}
 [(1+t)^{p+2i+2}-(1-t)^{p+2i+2}]}
 {2^{j+1}(p+2i+2)t^{j+1}}.
\end{equation}
Analyticity at zero cancels negative powers, and Legendre orthogonality
cancels powers below $\ell$. Orthogonality also proves termination for
even $p$. The leading coefficient for $p=2m-1$ is
$(-1)^\ell4^\ell\ell!\prod_{j=0}^{\ell-1}(m-\frac12-j)/(2\ell+1)!$.
It equals $c_0$ in \eqref{eq:lithium-integrals-odd-coefficients},
including when $\ell<m$. Projecting
$\Delta|x-e_3|^{2m-1}=2m(2m-1)|x-e_3|^{2m-3}$ onto degree $\ell$
shows that the remaining coefficients obey the displayed recurrence:
the radial operator sends $t^{\ell+2j}$ to
$2j(2\ell+2j+1)t^{\ell+2j-2}$.

Direct integration gives
$A_0^{(-2)}=\log((1+t)/(1-t))/(2t)$ and
$A_1^{(-2)}=(1+t^2)A_0^{(-2)}/(2t)-1/(2t)$.
The Legendre recurrence then proves the positive series
\eqref{eq:lithium-integrals-inverse-square} by induction. Its coefficient
is at most $1/(2j+1)$, which gives its logarithmic majorant.
All angular series converge absolutely at distinct radii.

For $\ell\ge m$, expansion of a beta integral gives
\[
 A_\ell^{(2m-1)}(t)=c_0t^\ell
 \frac{\int_0^1z^{\ell-m-1/2}(1-z)^m(1-t^2z)^m\,dz}
      {\int_0^1z^{\ell-m-1/2}(1-z)^m\,dz}.
\]
The quotient lies in $[0,1]$. This proves the sign $(-1)^m$ and the
bound $|A_\ell^{(2m-1)}(t)|\le|c_0|t^\ell$ used in the tails below.
The conditions on the clusters imply $C_\sigma,D_\sigma\ge0$.
After the scale integral, a finite angular term is bounded by a constant
times $s^{D_\sigma}v^{C_\sigma}$ times at most two logarithmic
majorants. For two distinct ratios in $\{s,v,sv\}$ their product is
bounded by $(1-\log(1-s))(1-\log(1-v))$. This is integrable.
An infinite angular sum has either three odd powers or one inverse-square
power and two odd powers. In the first case, the beta bounds and the
factor $2\ell+1$ give a summable $O(\ell^{-2})$ envelope after any
fixed initial segment. In the second, the two odd coefficients give
$O((2\ell+1)^{-2})$, and the product of all three ratios to power
$\ell$ is $(sv)^{2\ell}$. The remaining sum is bounded by a constant
times $1-\log(1-sv)$. Multiplication by the inverse-square majorant
again gives a bound by a constant times
$(1-\log(1-s))(1-\log(1-v))$, whichever edge has power $-2$.
The scale integral contributes a bounded factor, since its denominator
is at least the largest-radius rate. These bounds prove absolute
integrability and justify termwise integration.
\end{proof}

\begin{lemma}[Ordered radial integrals]\label{lem:lithium-integrals-radial}
Let $a,b,c>0$ be rational, and let $\alpha,\beta,\gamma$ be integers with
$\gamma\ge0$, $d=\beta+\gamma+1\ge0$ and
$s=\alpha+\beta+\gamma+3\ge1$. Define
\[
 J=\int_{0<z<y<x}x^\alpha y^\beta z^\gamma e^{-ax-by-cz}\,dz\,dy\,dx,
 \quad \lambda=a+b+c,\quad \xi=\frac{b+c}{\lambda},\quad
 \eta=\frac c\lambda.
\]
Then, with $(s)_j=s(s+1)\cdots(s+j-1)$ and $(s)_0=1$,
\begin{align}
 t_j&=\frac{(s)_j\gamma!}{(d+j+1)!}
       \sum_{i=0}^j\frac{(d+i)!}{(\gamma+i+1)!}\xi^{j-i}\eta^i,
 &f_j&=\frac{(s)_j\xi^j}{j!(d+1)(\gamma+1)},\notag\\
 J&\in\frac{(s-1)!}{\lambda^s}
 \left[\sum_{j=0}^N t_j,
       \sum_{j=0}^N t_j+\frac{f_{N+1}}{1-\rho_N}\right],
 &\rho_N&=\frac{\xi(s+N)}{N+1}<1.\label{eq:lithium-integrals-radial-series}
\end{align}
If $\alpha,\beta,\gamma\ge0$, the exact rational alternative is
\begin{equation}\label{eq:lithium-integrals-rational-chamber}
 J=\sum_{i=0}^{\alpha}\sum_{j=0}^{\beta+i}
 \frac{\alpha!(\beta+i)!(\gamma+j)!}
 {i!j!a^{\alpha+1-i}(a+b)^{\beta+i+1-j}\lambda^{\gamma+j+1}}.
\end{equation}
For a master with $n_i\ge-2$, $p_{ij}\ge-1$ all odd, and
\eqref{eq:lithium-integrals-clusters}, put
$m_{ij}=(p_{ij}+1)/2$, $M=\max m_{ij}$, $h=2+\sum m_{ij}$ and
\[
 b_L=\frac{\prod_{i<j}(2m_{ij}-1)!!}
 {2(h-1)(2L+1-2M)^{h-1}},\qquad L\in\mathbb Z,\quad L\ge M.
\]
The absolute error after angular index $L$ is at most $b_L$ times
the sum over the six radial orders of $J$ with baseline degrees
$(n_h+2+p_{hm}+p_{hb},n_m+2+p_{mb},n_b+2)$ and the shift
$(-2L-2,0,2L+2)$, with the correspondingly ordered rates $R_i/180$.
\end{lemma}

\begin{proof}
Set $x=R,y=Ru,z=Ruv$ and integrate $R$. The remaining denominator is
\[
 \lambda^s[1-\xi(1-u)-\eta u(1-v)]^s.
\]
Its binomial expansion has nonnegative terms. Integrating the term of
degree $j$ by beta integrals gives $t_j$, and bounding its argument by
$\xi$ gives $t_j\le f_j$. The ratios
$f_{j+1}/f_j=\xi(s+j)/(j+1)$ decrease to $\xi<1$. Hence the
geometric remainder is \eqref{eq:lithium-integrals-radial-series}.
Repeated integration from $x=y$ and then $y=z$ gives
\eqref{eq:lithium-integrals-rational-chamber}.

For the angular tail multiply the three beta-integral bounds in the
preceding proof. Cancel one factor $2\ell+1$ and bound the remaining
$h$ factors below by $2\ell+1-2M$. Since
$t_{12}t_{13}t_{23}=(r_b/r_h)^2$, summing the decreasing coefficient
by its integral on $[L,\infty)$ gives $b_L(r_b/r_h)^{2L+2}$.
Every resulting radial integral has nonnegative lowest degree and
middle-plus-lowest degree plus one, with total shape $\kappa$.
The same is true for the finite angular terms: their polynomial exponents
are nonnegative. This proves integrability and the asserted remainder.
\end{proof}

\begin{lemma}[Native radial arithmetic]\label{lem:lithium-integrals-native}
Assume the native execution model of
Set-up~\ref{setup:lithium-integrals-execution}. In
Lemma~\ref{lem:lithium-integrals-radial}, suppose
$(a,b,c)=(R_h,R_m,R_b)/180$ with $233\le R_i\le1116$,
$0\le\gamma,d\le256$ and $1\le s\le32$.
The positive interval recurrence described in the proof returns an
enclosure of $J$ within $4096$ iterations. It has no integer overflow,
invalid denominator or nonfinite numerical operand. For masters with
$-2\le n_i\le10$, $-1\le p_{ij}\le10$, $1\le\kappa\le32$ and
these rate bounds, the finite angular wrapper and the all-odd wrapper
at $M+1\le L\le24$ return enclosures.
Their sequential signed contraction remains enclosing and finite for
at most $2^{24}$ terms with rational coefficients of magnitude at most
$2^{128}$.
\end{lemma}

\begin{proof}
\emph{Outward recurrence.}
Write $\operatorname{RN}$ for nearest-even binary64 rounding and
$\operatorname{pred},\operatorname{succ}$ for the adjacent numbers.
For nonnegative expressions use
$\max\{0,\operatorname{pred}(\operatorname{RN}(x))\}$ and
$\operatorname{succ}(\operatorname{RN}(x))$ as lower and upper bounds.
Rounding cells lie between their adjacent representable numbers, including
at zero and in the subnormal range. Addition, multiplication and division
of positive intervals therefore enclose their exact operations.

Set $\chi=c/(b+c)$. The native loop maintains intervals for $t_j,f_j$
of \eqref{eq:lithium-integrals-radial-series}, its last summand
$\upsilon_j=(s)_j\gamma!\eta^j/[(d+j+1)(\gamma+j+1)!]$, and
$\sum_{i<j}t_i$. Initially $t_0=\upsilon_0=f_0=1/[(d+1)(\gamma+1)]$.
At iteration $j$ it first adds $t_j$, updates
$f_{j+1}=\rho_j f_j$, and, if the upper endpoint of $\rho_j$ is
less than one, bounds the remainder by
\[
 \operatorname{succ}\!\left(\operatorname{RN}
 \frac{(f_{j+1})_+}
 {\operatorname{pred}(\operatorname{RN}(1-(\rho_j)_+))}\right).
\]
The denominator is greater than $2^{-54}$. The loop returns when this
bound is less than the rounded product of the current lower sum and
the binary64 literal $10^{-15}$. It then multiplies the enclosing sum
by the interval for $(s-1)!/\lambda^s$, built by successive factors
$i/\lambda$. The comparison only chooses a return time. Enclosure
already follows from the remainder formula. Otherwise the updates are
\[
 a_j=\frac{\xi(s+j)}{d+j+2},\qquad
 b_j=\frac{\chi(d+j+1)}{\gamma+j+2},\qquad
 \upsilon_{j+1}=a_jb_j\upsilon_j,\quad t_{j+1}=a_jt_j+\upsilon_{j+1}.
\]
Cancellation of factorials in the exact series proves these identities,
and induction proves the interval invariant. A failure return is never
used as an interval.

\emph{Range and termination.}
All integer expressions, including those formed before the guard, have
magnitude at most $4353$ on this box. With $\mu=2^{-1074}$ an outward
scalar operation has absolute error at most $2^{-50}|x|+2\mu$.
The constructed upper ratio satisfies $\xi_+<29/32$ since
$\xi\le2232/2465$ and $29\cdot2465-32\cdot2232=61$.
An induction with
$K_j=4(1+2^{-20})^j\binom{s+j-1}{j}$ gives
\[
 (\upsilon_j)_+,(f_j)_+\le K_j,\qquad (t_j)_+\le(j+1)K_j,\qquad
 K_j<2^{406}\quad(0\le j\le4096).
\]
Indeed each recurrence factor incurs a relative multiplier below
$1+2^{-40}$. The larger multiplier $1+2^{-20}$ absorbs those errors
and the additive $2\mu$ terms. The accumulated sum is below $2^{432}$,
the scale below $2^{156}$, and the returned value below $2^{640}$.
For termination the sharper envelope bound is
\[
 (f_N)_+<4(29/32)^N\binom{s+N-1}{N}
       +4\mu N\binom{s+N-1}{N}.
\]
At $N=4096$ use $29^8<2^{39}$ and
$\binom{s+N-1}{N}<2^{403}$: the envelope is below $2^{-106}$.
At $j=4095$ the ratio is below $15/16$, the remainder below
$2^{-100}$, and the lower sum remains above $2^{-18}$.
The stopping threshold exceeds $2^{-69}$, so the return occurs.

\emph{Angular and signed arithmetic.}
The exact angular lists are
\eqref{eq:lithium-integrals-angular-polynomial} and
\eqref{eq:lithium-integrals-odd-coefficients}. Equal radial powers
are collected before conversion. Signed rational conversion uses both
neighbors of its rounded value. Signed multiplication takes the four
rounded corner products and then their exterior neighbors. Addition
widens each rounded endpoint once. These rules enclose through
cancellation. The odd wrapper adds the symmetric tail interval of
Lemma~\ref{lem:lithium-integrals-radial}, with its upper endpoint
rounded outwards.

Here $m_{ij}\le5$. The largest angular exponent is $L+2m_{ij}\le34$,
so each native call has lowest degree at most $80$ and combined degree
at most $103$. A tail call has bounds $62$ and $85$. The finite
wrapper has smaller bounds. The coefficient absolute sum of an odd
edge is at most $9!!\,4^5<2^{20}$; for even edges the Legendre
recurrence and binomial expansion give the same bound. There are fewer
than $2^{16}$ terms in a wrapper and every collected coefficient has
magnitude below $2^{80}$. The preceding outward-error estimate bounds
all master endpoints by $2^{740}$ and all stated signed-contraction
intermediates by $2^{900}$. This also proves finiteness of the outward
neighbors used throughout.
\end{proof}

\begin{lemma}[Finite logarithmic integrals]\label{lem:lithium-integrals-log}
Let $\alpha,\beta,\gamma$ be nonnegative integers and $a,b,c>0$ be
rational. On $x>y>z>0$ let $(r'_0,r'_1,r'_2)=(x,y,z)$ and choose
$0\le h<k\le2$. The integral
\[
 \int x^\alpha y^\beta z^\gamma e^{-ax-by-cz}
       \log\frac{r'_h+r'_k}{r'_h-r'_k}\,dz\,dy\,dx
\]
has the exact finite gamma-mixture evaluation and rational error bounds
given in the proof. They apply also when auxiliary gamma rates coincide.
\end{lemma}

\begin{proof}
\emph{Finite reduction.}
Put $g_0=x-y,g_1=y-z,g_2=z$. Expand
$(g_0+g_1+g_2)^\alpha(g_1+g_2)^\beta g_2^\gamma$.
Each monomial is a positive rational multiple of independent gamma
densities of integer shapes and rates $(a,a+b,a+b+c)$.
The logarithm is the difference of the logarithms of
$\sum_{i\ge h}g_i+\sum_{i\ge k}g_i$ and
$\sum_{h\le i<k}g_i$. Thus it suffices to integrate
$\log Y$ for $Y=\sum_i\ell_i X_i$, with integer gamma shapes and
nonnegative rational $\ell_i$, not all zero.

Merge equal positive rates $\rho_i/\ell_i$ into rates $v$ with
total shapes $K_v$. The Laplace transform is
$\prod_v(v/(v+t))^{K_v}$. At the pole $-v$ compute the first
$K_v$ Taylor coefficients $c_{v,j}$ of
$\prod_{w\ne v}(w/(w-v+z))^{K_w}$. A factor's coefficients obey
\[
 s_0=\left(\frac w{w-v}\right)^{K_w},\qquad
 s_j=-s_{j-1}\frac{K_w+j-1}{j(w-v)}.
\]
Finite convolution gives all $c_{v,j}$ exactly. Principal-part subtraction
gives the partial fractions with weights $w_{v,j}=v^jc_{v,j}$,
and $\sum_{v,j}w_{v,j}=1$ by evaluation at $t=0$.
Consequently,
\begin{equation}\label{eq:lithium-integrals-gamma-log}
 \mathbb E\log Y=\mathfrak c+
 \sum_{v,j}w_{v,j}\bigl(\mathsf h_{K_v-j-1}-\log v\bigr),
 \quad \mathfrak c=\int_0^\infty e^{-t}\log t\,dt.
\end{equation}
This follows by integrating the partial-fraction gamma densities.
Integration by parts gives their log moments. The identity of densities
can equivalently be obtained by inverse finite convolution of exponential
polynomials. Their absolute logarithmic moments are finite. The common
constant $\mathfrak c$ cancels in the required difference, so its
numerical value is unnecessary. Equal rates were merged before division.
Signed partial-fraction weights are retained.

\emph{Error and rounding.}
For $1\le v\le2$ put $z=(v-1)/(v+1)$. Then
\begin{equation}\label{eq:lithium-integrals-log-series}
 0\le\log v-2\sum_{j=0}^{127}\frac{z^{2j+1}}{2j+1}
 \le\frac{2z^{257}}{257(1-z^2)}.
\end{equation}
Integrating the geometric series proves this bound. Reduce any positive
rational rate to $2^kv$ with $1\le v<2$ and use the same formula for
$\log2$. Reverse endpoints if $k<0$. Let $\delta(v)$ be the resulting
width for that rate. Formula~\eqref{eq:lithium-integrals-gamma-log}
then has width at most
$\sum_v|\sum_jw_{v,j}|\delta(v)+2^{-119}$ after outward rounding
to multiples of $2^{-120}$. Multiply the two log-mean intervals by
the exact positive masses of the gap monomials, subtract with reversed
endpoints, and add. This proves the claimed ordered-log enclosure with
an explicit finite error, without a floating logarithm.
\end{proof}

\begin{lemma}[Compact inverse-square integrals]\label{lem:lithium-integrals-compact}
Let $\ell,D,C$ be nonnegative integers and let
$\epsilon\in\{(1,0),(0,1),(1,1)\}$. Define
\[
 \mathcal J_{\ell,\epsilon}(D,C)=
 \int_0^1\int_0^1s^Dv^C B_\ell(s^{\epsilon_1}v^{\epsilon_2})\,dv\,ds.
\]
Then $0<\mathcal J_{\ell,\epsilon}(D,C)\le2$, and it is a rational
linear combination of $1,\log2,6\zeta(2),\zeta(3)$.
For positive integer rate numerators $(R_h,R_m,R_b)$, integer shape
$m\ge1$ and integer precision $d\ge0$, the exact compact algorithm
specified in the proof encloses
\begin{equation}\label{eq:lithium-integrals-compact-ordered}
 (m-1)!\int_0^1\int_0^1
 \frac{s^Dv^CB_\ell(s^{\epsilon_1}v^{\epsilon_2})}
 {(a+bs+csv)^m}\,dv\,ds,\qquad
 (a,b,c)=(R_h,R_m,R_b)/180.
\end{equation}
Its arithmetic uses integers and rationals, with denominator scale
$U=2^{640}$ and the constant intervals below.
\end{lemma}

\begin{proof}
We evaluate the compact moments by partial fractions and approximate only
the radial denominator, with its integrated remainder included explicitly.

\emph{Compact moments and constants.}
The coefficient in \eqref{eq:lithium-integrals-inverse-square} has the
partial fractions
\[
 \frac{\prod_{i=1}^\ell(x+2i)}{\prod_{i=0}^\ell(x+2i+1)}
 =\sum_{h=0}^\ell\frac{a_{\ell h}}{x+2h+1},\qquad
 a_{\ell h}=
 \frac{\prod_{i=1}^\ell(2i-2h-1)}
      {\prod_{0\le i\le\ell,\ i\ne h}(2i-2h)}>0.
\]
Each active variable contributes a factor $(2j+D+1)^{-1}$ or
$(2j+C+1)^{-1}$ to the integrated positive series. An inactive one
contributes $1/(D+1)$ or $1/(C+1)$. Thus each term reduces to
\[
 \sum_{j\ge0}\prod_{v\in\mathcal P}(2j+v)^{-1},
\]
where $\mathcal P$ is a list of two or three positive integers, with
repetitions allowed. Partial fractions at a repeated pole $-v$ are
computed by the Taylor coefficients of
$\prod_{w\ne v}(w-v+z)^{-m_w}$. Their recurrence is the one in
the preceding proof without its rate numerators. The sum of the
simple-pole coefficients is zero, since the original rational function
is $O(x^{-2})$. Subtract $1/(2j+1)$ from each simple-pole series before
summing. Odd and even poles then give respectively finite harmonic
tails and $-\log2$ minus finite even harmonic tails. The square-pole
tails are $6\zeta(2)/8$ and $6\zeta(2)/24$ minus their finite
initial sums. The cube-pole tails are $7\zeta(3)/8$ and
$\zeta(3)/8$ minus their initial sums. This proves the linear-combination
formula, including triple poles. No divergent series is subtracted
after evaluation.

For two distinct active degrees use
\[
 \mathcal J_{\ell,(1,1)}(D,C)=
 \frac{\mathcal J_{\ell,(1,0)}(D,0)
       -\mathcal J_{\ell,(1,0)}(C,0)}{C-D}.
\]
For equal degrees use the repeated-pole formula. Outward integer division
implements this identity and division by the inactive degree plus one.
This also proves the correctness of caching these two reference types.
The bound $\mathcal J\le2$ follows from
$B_\ell(t)\le1-\log(1-t)$: the active ratio is at most one of
$s,v$, and $\int_0^1(1-\log(1-t))\,dt=2$.

Here are exact definitions of the constant inputs. Put
$\mathsf h_j^{(2)}=\sum_{i=1}^j1/i^2$, and set
\[
\begin{aligned}
 E_1&=\sum_{j=1}^{769}\frac{2^{-j}}j,&
 E_2&=\sum_{j=1}^{769}\frac{2^{-j}\mathsf h_j}j,&
 E_3&=\sum_{j=1}^{769}
       \frac{2^{-j}(\mathsf h_j^2+\mathsf h_j^{(2)})}{2j}.
\end{aligned}
\]
For $(\gamma_1,\gamma_2,\gamma_3)=(1,12,4/3)$ use the integer pairs
\begin{equation}\label{eq:lithium-integrals-constant-guard}
 a_i=\lfloor U\gamma_iE_i\rfloor,\qquad
 b_i=\lceil U\gamma_i(E_i+2^{-769})\rceil.
\end{equation}
They enclose $U(\log2,6\zeta(2),\zeta(3))$.
Indeed
\[
 b_{s,k}=\frac1{(s-1)!}\int_0^1(1-t)^k(-\log t)^{s-1}\,dt
 \in[0,1]
\]
and $\sum_{k\ge0}2^{-k-1}b_{s,k}$ is $\log2$ for $s=1$ and
$(1-2^{1-s})\zeta(s)$ for $s=2,3$. The omitted tail is at most
$2^{-769}$. Differentiating
$\int_0^1t^a(1-t)^k\,dt=k!/\prod_{j=1}^{k+1}(a+j)$ at $a=0$
gives the three displayed finite sums. Signed substitution of
\eqref{eq:lithium-integrals-constant-guard} and outward rounding
therefore returns integers enclosing $U\mathcal J$.

\emph{Denominator and its error.}
Put $c_0=a+(b+c)/2$, $h=(b+c)/(2c_0)<1$ and
$S_0=(m-1)!/c_0^m$. Expand
\[
 \frac{(m-1)!}{(a+bs+csv)^m}
 =S_0\sum_{k\ge0}\binom{m+k-1}{k}
       \left(h-\frac{bs+csv}{c_0}\right)^k.
\]
The argument has absolute value at most $h$. Write
$u_k=S_0\binom{m+k-1}{k}h^k$ and
$\rho_k=h(m+k-1)/k$. Choose the first $n\ge1$ with
$\rho_n<1$ and $\tau=u_n/(1-\rho_n)<10^{-d}$. Such $n$ exists.
The polynomial through degree $n-1$ has uniform error at most $\tau$.
Expand it as $\sum g_{ij}s^iv^j$ and set $w_{ij}=\lfloor Ug_{ij}\rfloor$.
If there are $N_\triangle$ coefficients, the integrated denominator error
is at most
\begin{equation}\label{eq:lithium-integrals-compact-error}
 \frac eU,\qquad e=\lceil2U\tau\rceil+2N_\triangle.
\end{equation}
This follows by multiplying the uniform Taylor and coefficient-floor
errors by a nonnegative weight whose integral is at most two.

For correspondence with the implemented polynomial shift, let
$P_N(x)=\sum_{j=0}^N\binom{m+j-1}{j}x^j$, $N=n-1$, and
$P_N(h+y)=\sum t_i y^i$. Horner evaluation gives $t_0$ and
\[
 t_{i+1}=\frac{(m+i)t_i-d_i}{(1-h)(i+1)},\qquad
 d_i=(m+N)\binom{m+N-1}{N}\binom Ni h^{N-i}.
\]
These identities follow by comparing coefficients in
$(1-x)P_N'-mP_N=-(m+N)\binom{m+N-1}{N}x^N$.
They give
$g_{ij}=S_0t_i(-1/(180c_0))^i\binom ijR_m^{i-j}R_b^j$.
The consecutive integer binomial updates are exact divisions.

If $[l_{ij},u_{ij}]$ encloses $U\mathcal J(D+i,C+j)$, form
\[
 v_- =\sum_{i,j}\min(w_{ij}l_{ij},w_{ij}u_{ij}),\qquad
 v_+ =\sum_{i,j}\max(w_{ij}l_{ij},w_{ij}u_{ij}).
\]
Then $[(v_--eU)/U^2,(v_++eU)/U^2]$ encloses
\eqref{eq:lithium-integrals-compact-ordered}. This is precisely the
integer contraction and error allowance used in the compact evaluations.
\end{proof}

\begin{lemma}[Signed odd-logarithmic tail]\label{lem:lithium-integrals-oddlog-tail}
Let a master have one pair power $-2$ and the other two powers
$2m_1-1,2m_2-1$, with nonnegative integers $m_i$ and
$\mu=m_1+m_2>0$.
Assume $n_i\ge-2$, positive integer rate numerators and the cluster
conditions \eqref{eq:lithium-integrals-clusters}. Put
$M=\max(m_1,m_2)$ and $\lambda=\sum_iR_i/180$.
Choose an integer $L\ge M$ such that for every radial order
\[
 d_\sigma=D_\sigma+2L+2-\kappa\ge0,\qquad
 c_\sigma=C_\sigma+2L+2-\kappa\ge0.
\]
Let $\mathcal I_L$ be the angular expression through $\ell=L$,
using the compact integrals of
Lemma~\ref{lem:lithium-integrals-compact}. Set
\[
 v_\sigma=
 \begin{cases}c_\sigma,&\text{the inverse-square edge is }mb,\\
 d_\sigma,&\text{otherwise},\end{cases}
 \qquad
 T_\sigma=
 \frac{(2m_1-1)!!(2m_2-1)!!(\kappa-1)!
              (1+\mathsf h_{v_\sigma+1})}
 {2\mu(2L+1-2M)^\mu\lambda^\kappa
              (d_\sigma+1)(c_\sigma+1)}.
\]
Then
\begin{equation}\label{eq:lithium-integrals-oddlog-tail}
 0\le(-1)^\mu(\mathcal I-\mathcal I_L)\le\sum_\sigma T_\sigma.
\end{equation}
\end{lemma}

\begin{proof}
For $\ell>L$, the two odd angular coefficients have signs
$(-1)^{m_1},(-1)^{m_2}$ by the beta-integral formula, and the
inverse-square coefficient is positive. Cancel one factor $2\ell+1$
from their coefficient bound. The remaining sum is bounded by
\[
 \sum_{\ell>L}(2\ell+1-2M)^{-\mu-1}
 \le[2\mu(2L+1-2M)^\mu]^{-1}.
\]
The radius-ratio factor is $(sv)^{2L+2}$. Also
$a+bs+csv\ge\lambda sv$, and
$B_\ell(t)\le1-\log(1-t)$, with $t=v$ on edge $mb$ and
$t\le s$ otherwise. The remaining two monomial integrals are
$1/(j+1)$ and
$(1+\mathsf h_{j+1})/(j+1)$. They give $T_\sigma$.
The bound is integrable, justifying the angular integration and retaining
its common sign. The finite part uses
$A_\ell^{(-2)}=t^\ell B_\ell(t)$ and the two rational polynomials.
Their exponents ensure $D,C\ge0$ in every compact call.
\end{proof}

\begin{lemma}[Double-logarithmic quadrature]\label{lem:lithium-integrals-doublelog}
Let $1\le a,b,c\le7$ be positive rational rates, let $1\le m\le16$,
$0\le X,Y\le16$ and $0\le\ell\le3$ be integers with
$X+Y+4\ell\le36$, and let $t_1,t_2$ be distinct members of
$\{x,y,xy\}$. The order-ten tensor Gauss algorithm specified in the
proof, with endpoint cutoff $2^{-34}$, encloses
\begin{equation}\label{eq:lithium-integrals-doublelog-integrand}
 \int_0^1\int_0^1
 \frac{x^Xy^YA_\ell^{(-2)}(t_1)A_\ell^{(-2)}(t_2)}
 {(a+bx+cxy)^m}\,dx\,dy
\end{equation}
under the NumPy execution model of
Set-up~\ref{setup:lithium-integrals-execution}.
Consequently its finite signed angular assembly encloses every master
with pair powers a permutation of $(-2,-2,2k)$,
$k\in\{0,1,2,3\}$, whose
ordered terms satisfy these bounds and whose nonzero angular
coefficients have magnitude below $2^{45}$.
\end{lemma}

\begin{proof}
We bound the Gauss remainder and omitted strips, then verify the
floating allowances for the angular functions and the positive sums.

\emph{Rule and analytic error.}
The ten nodes are enclosed by rational sign brackets for $P_{10}$,
starting with width $2\cdot10^{-12}$ and bisecting $48$ times.
The brackets must be disjoint, and their derivative lower bounds
positive. Proposals for the brackets have no assumed accuracy.
The polynomial is
\[
 P_{10}(t)=
 \frac{46189t^{10}-109395t^8+90090t^6-30030t^4+3465t^2-63}{256}.
\]
Its ten sign changes certify all ten roots. Since
$|P_{10}''|\le1485<10^4$ on $[-1,1]$, a bracket $[l,u]$ gives
$|P_{10}'(r)|\in[|P_{10}'((l+u)/2)|-10^4(u-l)/2,
|P_{10}'((l+u)/2)|+10^4(u-l)/2]$.
Use this interval and the endpoint bounds for $r^2$ in
$w(r)=2/[(1-r^2)P_{10}'(r)^2]$ to enclose each positive weight.
The exact sign/weight certificate in the supplement gives
\[
 |r|<.974,\qquad .066671<w(r)<.295525.
\]
Its rational comparisons use the displayed polynomial and derivative
formula. Any successful final source bracket has width
$1/(140737488355328000000000000)<10^{-26}$, and hence lies in
$(-.975,.975)$. The weight bounds change by less than $3\cdot10^{-14}$.
Thus all source brackets lie in $(-.99,.99)$ and their weight
endpoints in $(1/16,1)$. Orthogonality and division by $P_{10}$ prove
the rule exact through degree $19$. The weight formula follows by
integrating its Lagrange cardinal polynomials.

The mesh consists of $0,1/32,\ldots,1/2$, followed by the midpoint
and right endpoint of each
$[1-2^{-j},1-2^{-j-1}]$, $1\le j\le33$.
There are $82$ panels and $820$ nodes in each direction, hence
$N=672400$ tensor nodes. Let $\delta=2^{-34}$.
On a panel of midpoint $z_0$ and half-width $h$, use the complex disk
$|z-z_0|\le4h$. The positive upper bound on $|z|$ is less than
$1-2^{-36}$, the lower real-part bound is at least $-3/64$, and
the real part of $a+bx+cxy$ is at least $11/32$ when one variable
is in its disk and the other in its real panel.

Write $\Lambda(t)=1-\log(1-t)$. The positive series gives
$|A_\ell^{(-2)}(z)|\le\Lambda(|z|)$ for $|z|<1$.
For two distinct ratios their product is bounded by
$\Lambda(x_+)\Lambda(y_+)$, where $x_+,y_+$ are the disk/real
upper moduli. If $d_-$ is the denominator lower bound on that disk,
an analytic majorant is
$x_+^Xy_+^Y\Lambda(x_+)\Lambda(y_+)/d_-^m$.
Let $\mathcal M$ be the sum of rectangle areas times these majorants
for both choices of the complex variable, replacing
$\Lambda(t)$ by the integer upper bound $3-e$ when
$\operatorname{frexp}(1-t)$ has exponent $e$.
Cauchy's coefficient bound gives a Taylor remainder
$(4/3)4^{-20}$ times the majorant. The integral and positive rule
each have mass equal to the panel length. Applying
$(I_x-Q_x)I_y+Q_x(I_y-Q_y)$ on each rectangle therefore gives error
$(8/3)4^{-20}\mathcal M$.

Define $M_Z(A,B)=\max_{0\le t\le1}t^Z/(A+Bt)^m$.
For $Z=0$ it is $A^{-m}$. Otherwise the maximizer is $1$ if
$m\le Z$ or $B=0$, and
$\min\{1,ZA/((m-Z)B)\}$ in the remaining case.
Differentiation proves this formula, so $M_Z$ is evaluated rationally.
The two omitted strips have total integral at most
\begin{equation}\label{eq:lithium-integrals-doublelog-strip}
 2\delta(2+34)
 \left[M_Y(a+b(1-\delta),c(1-\delta))
       +M_X(a,b+c(1-\delta))\right].
\end{equation}
Here $\int_0^1\Lambda=2$ and
$\int_{1-\delta}^1\Lambda=\delta(2-\log\delta)\le36\delta$.
An overlap of the two strips only enlarges the bound.

\emph{Angular evaluation.}
Put $u=2^{-53}$ and $\gamma_j=ju/(1-ju)$.
All node-ratio endpoints lie in $(2^{-27},1-2^{-35})$.
For $t<5/8$, evaluate the positive series for
$A_\ell^{(-2)}(t)$ through $j=64$ by Horner's rule. The lower
polynomial has terminal coefficient $c_{\ell,64}$ and the upper
one has $c_{\ell,64}+2$, since the omitted tail is less than
$2t^\ell(t^2)^{64}$. The operation and coefficient-conversion error
is at most $\gamma_{200}<3\cdot10^{-14}$ relatively. Multiplication
by $1\mp2\cdot10^{-12}$ covers it in the respective directions.

For $t\ge5/8$, put $z=(1+t^2)/(2t)$ and evaluate
\[
 Q_0=\log(1+t)-\log(1-t),\quad Q_1=zQ_0-1,\quad
 Q_{j+1}=\frac{(2j+1)zQ_j-jQ_{j-1}}{j+1},
 \quad A_j^{(-2)}=Q_j/(2t).
\]
Each logarithm is computed by exact binary argument reduction and
$2\sum_{i=0}^{24}w^{2i+1}/(2i+1)$, with
$w=(v-1)/(v+1)$ and $1\le v<2$. Its omitted tail is bounded by
$2/[51\,3^{51}(1-1/9)]$. The binary64 constant for $\log2$ is
$6243314768165359/2^{53}$. The same finite series at $w=1/3$
places $\log2$ strictly between this value and this value plus $u$.
Counting the argument, Horner and final-addition errors gives less
than $400u$ absolute error per logarithm and $900u$ for $Q_0$.
Since $1\le z\le89/80$ and $0<Q_j<25$, the first recurrence
has local error below $200u$ and the next two below $1000u$ each.
Start the error allowance at $2\cdot10^{-12}$, update it first as
$1.2E+2\cdot10^{-12}$ and then as
$2.3E+E_{\rm previous}+2\cdot10^{-12}$.
Induction bounds the actual error by half this allowance, including
rounding of the allowance. Subtracting/adding it, dividing by $2t$
and multiplying by $1\mp10^{-14}$ gives enclosing endpoints.
Both branches have lower endpoint greater than $t^\ell/8$ and
upper endpoint less than $32$. Clipping at zero changes neither.

\emph{Positive reductions and directed assembly.}
Use lower node, weight, monomial and angular endpoints for the lower
tensor sum and the upper denominator endpoint, and reverse those
choices for the upper sum. Before reduction the relative error in
each positive summand is below $\gamma_{400}$. Including the
nonnegative reduction and coefficient multiplication gives
$\gamma_{N+512}<2u(N+512)$. The source multiplier
\begin{equation}\label{eq:lithium-integrals-doublelog-rounding}
 1\mp32u(N+1024)
\end{equation}
covers this error and the final rounding in the indicated direction.
For the disk denominator, possible cancellation has absolute error
at most $21\gamma_5$ and relative error below $\gamma_{1024}$,
using its lower bound $11/32$. Powering and inversion cost less
than $\gamma_{33000}$. The two disk reductions and all remaining
operations give a combined bound $\gamma_{60000}<7\cdot10^{-12}$
for the error expression. Its source factor $1+10^{-9}$ therefore
makes the sum of \eqref{eq:lithium-integrals-doublelog-strip} and
$(8/3)4^{-20}\mathcal M$ an upper bound after floating evaluation.

These relative estimates are valid throughout: node coordinates exceed
$2^{-13}$, ratio coordinates exceed $2^{-27}$, weights exceed
$2^{-41}$, and the positive denominator lies in $[1,21]$.
Using $X+Y+4\ell\le36$ bounds every positive tensor intermediate
between $2^{-700}$ and $2^{100}$. Disk and strip intermediates lie
between $2^{-400}$ and $2^{110}$. Positive angular intermediates
also stay normal, and exact zeros in argument reduction stay zero.
Thus no harmful underflow or overflow was assumed in the relative
estimates. Finally subtract/add the error using the adjacent exterior
binary64 number, reverse endpoints for a negative angular coefficient,
and add with the same directed rule. At most $60$ angular terms occur
in a master. These last signed operations remain enclosing even at
zero or subnormal cancellation. This proves the assertion.
\end{proof}

\begin{lemma}\label{lem:lithium-integrals-expansion}
For the trials and sources of Definition~\ref{def:li-trials}, all functions
belong to $D(H_3(3))$, the sources have norm one, and the spin representation
preserves their inner products and Hamiltonian inner products. Every moment
specified by a pair in Definition~\ref{def:lithium-integrals-arrays} equals
\begin{equation}\label{eq:lithium-integrals-physical-sum}
 (4\pi)^3 C_c\sum_k\frac{W_c(k)}{D_c}\mathcal I(k).
\end{equation}
Every integral in this sum converges absolutely. All its distance powers
are at least $-2$, and an inverse-square pair power occurs only in a
base residual-norm column.
\end{lemma}

\begin{proof}
We identify the Hamiltonian dictionaries, justify their domains, and then
check the spin and scalar normalizations.

\emph{Hamiltonian dictionaries.}
Use the six oriented edges $(1,0),(2,0),(3,0),(1,2),(1,3),(2,3)$,
where $x_0=0$. For $e=(j,k)$ write $y_e=x_j-x_k$,
$\epsilon_{ei}=\mathbf1_{j=i}-\mathbf1_{k=i}$,
$w_e=\sum_i\epsilon_{ei}^2$, and
$g_{ef}=\sum_i\epsilon_{ei}\epsilon_{fi}$.
Put $a_e=a_i$ on $(i,0)$ and zero on the pair edges, and put
$v_e=-3$ on nuclear edges and $1$ on pair edges. Direct differentiation
of a scalar primitive $F=F_{a,\nu}$ gives
\begin{align}
 (H_3(3)+15/2)F=F\bigg[&\frac{15}{2}
 -\frac12\sum_e w_ea_e^2
 +\sum_e\frac{v_e+w_ea_e(\nu_e+1)}{d_e}
 -\frac12\sum_e\frac{w_e\nu_e(\nu_e+1)}{d_e^2}\notag\\
 &-\sum_{e<f}g_{ef}
  \left(\frac{\nu_e}{d_e^2}-\frac{a_e}{d_e}\right)
  \left(\frac{\nu_f}{d_f^2}-\frac{a_f}{d_f}\right)y_e\cdot y_f
 \bigg].\label{eq:lithium-integrals-hamiltonian}
\end{align}
For a Cartesian primitive one adds the product-rule term
\begin{equation}\label{eq:lithium-integrals-cartesian}
 (H_3(3)+15/2)(x_{j,3}F)
 =x_{j,3}(H_3(3)+15/2)F
 -F\sum_e\epsilon_{ej}
       \left(\frac{\nu_e}{d_e^2}-\frac{a_e}{d_e}\right)y_{e,3}.
\end{equation}
The identity
$2y_{(j,k)}\cdot y_{(s,t)}=d_{jt}^2+d_{ks}^2-d_{js}^2-d_{kt}^2$
turns these into finite rational dictionaries. A power lowered by two
has a nonzero coefficient only if its original power is positive.
Thus every output monomial has powers at least $-1$ and at most two
negative powers. In a product, all powers are at least $-2$, at most two
can equal $-2$, and the sum of all six powers is at least $-4$.
Cartesian averaging increases one power by two. In particular,
\eqref{eq:lithium-integrals-clusters} holds term by term.
A product with an undifferentiated trial or source
cannot have a pair power $-2$.

\emph{Domains and convergence.}
The primitives and their first derivatives are bounded by an exponentially
decaying polynomial. Each second derivative is bounded on compact sets by
a constant times $1+\sum_e d_e^{-1}$. These bounds are square integrable:
each collision has three transverse coordinates. Cut off tubes of radius
$\varepsilon$ around the six collisions. Their volumes on compact sets
are $O(\varepsilon^3)$, whereas cutoff gradients are
$O(\varepsilon^{-1})$. The integration-by-parts boundary errors against
bounded first derivatives are $O(\varepsilon^2)$. Thus the displayed
derivatives are also the distributional derivatives. The functions are
in $H^2$, and their Coulomb potential products are in $L^2$.
The distributional identity therefore gives the defining form identity
for $D(H_3(3))$, first against smooth test functions and then by
$H^1$ density. Indeed integration of
$\operatorname{div}(x/|x|^2)=|x|^{-2}$ and Cauchy--Schwarz give
$\||x|^{-1}\phi\|_2\le2\|\nabla\phi\|_2$. Applying this also
in relative collision coordinates controls each Coulomb form term.

For the expanded products, Lemma~\ref{lem:lithium-integrals-angular}
gives absolute convergence. To check its remaining hypothesis, two
pair powers $-2$ require two Hamiltonian monomials each having the
same two pair powers $-1$. Such a monomial comes from differentiating
two pair factors of power one and retaining the square of the third
distance in the polarization identity. The $S$ table has no term
with all three pair powers positive, so the third power in the
product is four. A $P$ scalar feature has degree at most two, so
the third power is four or, after Cartesian averaging, six.
It is therefore even. This proves termwise integrability without
cancellation of singular monomials.

\emph{Spin, sources and measure.}
Two signed permutations producing the same spin string differ by the
exchange of the two up spins. The identity $P_{13}f=-f$ makes their
prescriptions agree. The three nonzero components have squared norms
$\|f\|^2/3$, so their sum is $\|f\|^2$. Polarization and commutation
of $H_3(3)$ with permutations prove the asserted isometry also for graph
products. No assertion about total spin is used.

The normalized $1s$ orbital is $Q^{3/2}e^{-Qr}/\sqrt\pi$.
The $nS$ orbital is
\[
 \frac{(2Q/n)^{3/2}}{\sqrt2 n\sqrt{4\pi}}
 e^{-Qr/n}L_{n-1}^1(2Qr/n),
\]
where
$L_{n-1}^1(y)=\sum_{j=0}^{n-1}(-1)^j\binom n{j+1}y^j/j!$.
Integration by parts in its Rodrigues formula gives
$\int_0^\infty e^{-y}y^2(L_{n-1}^1(y))^2\,dy=2n^2$.
The overlap with $1s$ vanishes for $n\ge2$ by the same finite
polynomial integral. For $2P$ and $3P$, integrate $x_3^2$
as $r^2/3$ and use $\int_0^\infty r^j e^{-ar}\,dr=j!/a^{j+1}$.
Their overlap with $1s$ is zero by odd parity. Multiplying the three
normalized orbitals and applying $A_{13}/\sqrt2$ gives exactly the
positive constants $C_{T,n}$ in Definition~\ref{def:li-trials}.
It proves unit source norms and preserves the stated source phases.

Rotational invariance of the scalar factors gives
$\int x_{i,3}x_{j,3}G=\frac13\int x_i\cdot x_jG$,
which proves the factor $1/6$ in the array definition. A simultaneous
electron permutation leaves a scalar integral unchanged. It introduces
neither a sign nor an orbit-size factor. All antisymmetrizer signs are
already in the coefficients. Finally, the three normalized angular
measures in \eqref{eq:lithium-integrals-master} require exactly the
restoration factor $(4\pi)^3$. This proves
\eqref{eq:lithium-integrals-physical-sum}.
\end{proof}

\begin{lemma}[Physical factors]\label{lem:lithium-integrals-factors}
Put $A=(4\pi)^3$, $F_{T,n}=AC_{T,n}$ and $B_G=[k_G/2^{128},(k_G+1)/2^{128}]$
for each factor in Table~\ref{tab:lithium-integrals-factors}.
Then $G\in B_G$ for all eight factors.
\end{lemma}

\begin{table}[htbp]
\centering\small
\caption{Integer numerators for the physical factor intervals.}
\label{tab:lithium-integrals-factors}
\begin{tabular}{c r}
\hline
$G$&$k_G$\\\hline
$A$&675256909963386956951116526350654582672596\\
$F_{S,2}$&1251738287056864148917692673584208527429128\\
$F_{S,3}$&454240014043670239193980669330780944972364\\
$F_{S,4}$&221278157762185462735926447436975349160194\\
$F_{S,5}$&126667008689195376140461016289088989981668\\
$F_{S,6}$&80299048554137951591401791333167863507093\\
$F_{P,2}$&3337968765484971063780513796224556073144343\\
$F_{P,3}$&1978055564731834704462526694058996191492944\\\hline
\end{tabular}
\end{table}

\begin{proof}
Squaring the positive constants gives
\[
 F_{S,n}^2=\frac{2^{38}\pi^3}{3^9n^5},\qquad
 F_{P,2}^2=\frac{2^{39}\pi^3}{3^{11}},\qquad
 F_{P,3}^2=\frac{2^{47}\pi^3}{3^{17}}.
\]
The rational arithmetic certificate in the supplement bounds $\pi$
between $P_*/2^{256}$ and $(P_*+1)/2^{256}$, where
\[
 P_*=363771576891766324280234942777729862653393377328392429958772151117938894466185.
\]
Its construction uses
\[
 \pi=16\arctan(1/5)-4\arctan(1/239),\qquad
 0<\arctan(1/v)-\sum_{j=0}^{63}\frac{(-1)^j}{(2j+1)v^{2j+1}}
   <\frac1{129v^{129}}.
\]
The first identity follows by multiplying
$(5+i)^4(239-i)=114244(1+i)$ and fixing the angle in $(0,\pi/2)$.
The second follows by integrating a finite geometric series.
For completeness, the integer certificate lists the quotients
$f_{v,j}=\lfloor2^{288}/((2j+1)v^{2j+1})\rfloor$, their positive
remainders and complementary remainders, and their signed prefixes.
If $S_v=\sum_{j=0}^{63}(-1)^jf_{v,j}$, it verifies
\[
 2^{32}P_*<16(S_5-32)-4(S_{239}+33),\qquad
 16(S_5+33)-4(S_{239}-32)<2^{32}(P_*+1).
\]
These are exact integer inequalities and prove the stated $\pi$ interval.
Cubing its endpoints and testing the squares of $k_G/2^{128}$ and
$(k_G+1)/2^{128}$ proves the seven radical bounds. The bound for $A$
uses $64\pi^3$ directly. The complete positive cross-product
remainders are supplied with the factor table. Positive square roots
preserve both inequalities, so no source phase changes in this operation.
\end{proof}

\subsection{Complete finite certificates}

We now define the family partition used by the computations. With no
pair power $-2$, the $P$ families are \emph{Pol} (all three pair
powers even), \emph{Fin} (some but not all even), and \emph{Odd}
(all odd). Their union is \emph{Splain} for $S$.
With two pair powers $-2$ the family is \emph{Doublelog}.
With exactly one, use \emph{Oddlog} when the other two powers are
odd. In the remaining case use the following exact predicate.

Write $\Lambda_0(t)=\log((1+t)/(1-t))$ and expand
$A_\ell^{(-2)}=U_\ell\Lambda_0+V_\ell$ as Laurent polynomials:
\[
 U_0=(2t)^{-1},\quad V_0=0,\quad
 U_1=(1+t^2)/(4t^2),\quad V_1=-(2t)^{-1}.
\]
Both sequences obey
$(\ell+1)Z_{\ell+1}=(2\ell+1)(1+t^2)Z_\ell/(2t)-\ell Z_{\ell-1}$.
For every radius order, multiply these lists and the other angular
polynomials through the finite even cutoff, include the radial
Jacobian, collect equal monomials and logarithm labels, and remove
zero coefficients. Call the tuple \emph{supported} if each of the
six lists is nonempty and every surviving term has at most one log
label, degrees $(\alpha,\beta,\gamma)$ satisfying
$\gamma\ge0$, $\beta+\gamma+1\ge0$, $\alpha+\beta+\gamma+3>0$,
and all three degrees nonnegative when a log label is present.
The supported part is \emph{Log}. Its complement is \emph{Single}.
This predicate is a finite rational calculation, with no integral
evaluation. It tests the collected Laurent expression, not the
number of inverse-square factors alone.

\begin{proposition}[Evaluated family certificates]\label{prop:lithium-integrals-certificates}
Assume Set-up~\ref{setup:lithium-integrals-execution} for the retained
executions and readers. For the exact arrays of
Definition~\ref{def:lithium-integrals-arrays}, the family partition is
exhaustive and has the sizes in Table~\ref{tab:lithium-integrals-coverage}.
The complete raw interval tables in the accompanying code/data release,
denoted $J_{T,c,E}$ for family $E$, contain
\[
 X_{T,c,E}=\sum_{k\in\mathcal L_T\cap E}\frac{W_c(k)}{D_c}\mathcal I(k).
\]
For Splain the supplied intervals $R_c$ instead contain the physical
quantities $(4\pi)^3C_cX_{S,c,\mathrm{Splain}}$.
The seven extra moments have the physical interval certificates described
in the proof. All statements concern the fixed trials simultaneously.
\end{proposition}

\begin{table}[htbp]
\centering\small
\caption{Complete base-family coverage and retained evaluation parameters.}
\label{tab:lithium-integrals-coverage}
\begin{tabular}{l r r l}
\hline
Family&$S$ keys&$P$ keys&Evaluation\\\hline
Splain&61344&---&finite or odd, $L=14$\\
Pol&---&17277&exact rational\\
Fin&---&108710&$17$ decimal digits, $96$-bit endpoints\\
Odd&---&15165&$L=24$, $120$-bit endpoints\\
Log&6983&21318&$17$ digits for plain terms, $96$-bit endpoints\\
Single&3128&7749&compact, $18$ digits; one $15$-digit group\\
Oddlog&2248&3376&$L=8$, compact $13$ digits, $200$-bit endpoints\\
Doublelog&50&72&Gauss order $10$, endpoint power $34$\\\hline
Total&73753&173667&\\\hline
\end{tabular}
\end{table}

\begin{proof}
The proof checks the exact input arrays, applies the integral error bounds
to all recorded calls, and reconstructs every signed contraction.
The supplement contains the full arrays and group/row records used in
these comparisons, together with the independent readers. The counts
in Table~\ref{tab:lithium-integrals-coverage} summarize those full
comparisons.

\emph{Common compiler and input identity.}
The compiler implements
\eqref{eq:lithium-integrals-hamiltonian}--\eqref{eq:lithium-integrals-cartesian}
by rational dictionary addition. Its exchange reverses nuclear slots
$1,3$, interchanges pair slots $12,23$, and also exchanges the
Cartesian label in $P$. In a square, it either enumerates both ordered
products or doubles an off-diagonal group once. Cartesian products
include both copies of $r_i^2$ when $i=j$, giving $r_i^2/3$.
Its integer packing is exact tuple addition: the scalar power digits
after offset $3$ lie in $[1,25]$ in base $32$. The vector digits after
offset $4$ and the square increment lie in $[2,24]$. Rate sums are
at most $1116<2048$, the rate-packing base. Thus there is no carry
or borrow into another field. Signed collision addition and removal
only of whole zero vectors implement the array definition exactly.

The independent input reader compared every rational primitive and
every mode, both entire generated and retained base arrays, and each
denominator with the least-common-multiple definition. Headers have
all nine key fields and all eight or five signed weights. Widths are
exactly $17$ and $14$. It checks strictly increasing canonical keys,
compares both streams to exhaustion, and rejects a missing or extra
row. For $P$ it also accumulated all $450585$ raw rows by signed
canonical addition and obtained the same $173667$-key array.
For every family, the reader removes a key from a copy of the full
map only after comparing its complete vector and its exact selector.
The remaining map is empty at the end. A repeated key cannot be
removed twice. This proves coverage and disjointness, including zero
entries within active vectors.

The corresponding extra arrays have union sizes $6830$ and $9878$.
Their nonzero column sizes are $(4256,450,579,708,837)$ and
$(9055,823)$, with denominators
\[
 (10^{28},81\cdot10^{14},81\cdot10^{14},151875\cdot10^{14},885735\cdot10^{14}),
 \qquad (6\cdot10^{36},54\cdot10^{18}).
\]
A separate uncompressed product expansion compares every extra signed
coefficient. The reader also compares the original first-insertion
dictionary order used for floating summation, not merely its sorted
map. Source comparisons retain every executable expression. Only
locations, unused docstrings and one diagnostic character-encoding
difference are irrelevant. Completed compiler stages are linked to
their actual outputs by their synchronous return and full later
comparisons. One scalar driver failed after completing its compiler
stage; that later failure is not represented as a successful whole run.

The structural power bounds in
Lemma~\ref{lem:lithium-integrals-expansion} exclude three inverse-square
pairs and put all inverse-square families in residual norms only.
Two such pair powers in one Hamiltonian monomial arise from differentiating
two pair factors of power one and retaining the square of the third
distance. The $S$ table has no monomial with three positive pair
powers and has total degree at most five. The $P$ scalar features
have degree at most two. Hence the Doublelog powers are $(-2,-2,4)$
with nuclear degree at most six in $S$, and $(-2,-2,4)$ with nuclear
degree two or $(-2,-2,6)$ with nuclear degree zero in $P$.
These restrictions also exclude Doublelog from the supported selector:
when the shared endpoint of the inverse-square edges has middle
radius, the coefficient of the two log labels, at $s=v=1/2$, is
$243/2048$ for the remaining power four and $3645/65536$ for power
six. Both are nonzero. Thus the concrete source selectors give exactly
the partition stated above.

\emph{Exact polynomial, finite and Log families.}
For Pol expand each even distance power by the multinomial theorem.
The angular moment of $u^j$ against $P_\ell$ is
$j!/[(j-\ell)!!(j+\ell+1)!!]$ when $j\ge\ell$ and $j-\ell$ is
even, and zero otherwise. Rodrigues integration by parts proves this
formula. The angular triple contraction then leaves independent
radial gamma integrals. Exact evaluation of all $17277$ selected
keys and all five signed columns gives the supplied Pol singletons.
The independent reader checks their full input and rational outputs.

For Fin, the even cutoff and
Lemma~\ref{lem:lithium-integrals-radial} give an exact finite angular
sum. Nonpolynomial radial calls stop only when their explicit tail
is at most $10^{-17}$, before the $10000$-iteration cap. Every
recorded call returned. Grouping uses a positive common angular
denominator. The reader reconstructs the same group independently
as a rational vector sum and deletes each observed nonzero radial
key only after comparing all five coefficients. All $59$ rate groups
and $50613$ nonzero radial calls were covered. The producer completed
these evaluations, output writes and all endpoint comparisons before
a later environment-metadata error. The separate reader exhaustively
checked those completed records. No failed integral return is used.

For Log, the supported predicate gives exactly the hypotheses of
Lemma~\ref{lem:lithium-integrals-log} on log terms and of
Lemma~\ref{lem:lithium-integrals-radial} on plain terms.
The plain calls use $17$ digits. Logarithmic calls use the fixed
$128$-term and $120$-bit bounds in their lemma. The $12$ scalar
rate groups give $10025$ nonzero radials ($8273$ logarithmic and
$1752$ plain). The $59$ vector groups give $48797$ ($42060$ and
$6737$). Empty groups remain in the manifest. The complete independent
regrouping and recorded-call comparisons cover every term.
In both Fin and Log each rational radial interval is rounded outward
to multiples of $2^{-96}$, multiplied with endpoint reversal by
its signed group coefficient, and divided once by the angular and
column denominators. Reconstructing all zero-start group sums gives
exactly the unchanged supplied raw endpoint pairs.

\emph{Native Odd, Splain and three extra moments.}
On the entire native base support, the rates lie in $[233,1116]$,
the nuclear powers in $[-2,10]$, the pair powers in $[-1,10]$,
and the shapes in $[5,19]$. Each column has fewer than $2^{23}$
terms and coefficient absolute sum less than $2^{52}$.
The extra native inputs have nonnegative powers at most ten,
shape in $[9,19]$, fewer than $2^{18}$ terms and coefficient
absolute sum less than $2^{26}$. These bounds follow by taking
products of the primitive dictionaries, and are also checked on
the complete consumed arrays. Thus all calls are in
Lemma~\ref{lem:lithium-integrals-native}. The cutoffs are $24$ for
$P$ Odd, $14$ for the odd part of Splain, and $12$ for the odd
parts of the extra norms. The remaining calls use the finite branch.

The observers call the unchanged native/wrapper function, record its
actual returned hexadecimal and rational endpoints, and return those
same values. Their before/after records identify the loaded library,
interface widths, rounding state and disabled flush modes. Such
observations instantiate, rather than prove, the execution assumptions.
The independent readers use integer arithmetic to decode every endpoint
and reconstruct every signed prefix from zero. For $P$ Odd they
check $A_j=\lfloor2^{120}l_j\rfloor$ and
$B_j=\lceil2^{120}u_j\rceil$, then add
$\min(W_cA_j,W_cB_j)$ and $\max(W_cA_j,W_cB_j)$.
They verify all $15165$ rows, all $949$ retained prefixes, and
equality of the five final pairs with the supplied raw intervals.

For Splain they check all $61344$ distinct selected masters and
$121$ prefixes, reconstructing both the exact rational contraction
of returned intervals and the original floating accumulation.
Nearest-even decoding is implemented by determining
$2^e\le|x|<2^{e+1}$ and rounding $|x|$ to the grid
$2^{\max(e-52,-1074)}$, with an even quotient at a midpoint.
The sign and adjacent bit patterns then reconstruct each outward
operation, including zero and subnormals. This proves the meaning
of the independent rounding reconstruction.
Multiplication by the factor intervals of
Lemma~\ref{lem:lithium-integrals-factors} puts all eight reconstructed
Splain intervals inside their saved physical $R_c$. Both signed
cross-product margins are positive in every column. In particular
the negative Splain residual-norm contribution is retained.
No interval-library normalization claim is needed.

The same method covers the three extra native columns in their
original source order: $4256$ scalar-norm calls ($4027$ finite and
$229$ odd), $9055$ vector-norm calls ($8890$ finite and $165$ odd),
and $823$ finite $P$ source-3 overlap calls. All $14134$ calls and
$115$ prefixes are checked. The exact physical scaling gives
intervals inside both the saved dyadics and the final printed
targets for these three moments.

\emph{Single.}
Keep the inverse-square angular factor as $t^\ell B_\ell(t)$
and expand the two other factors through their finite even cutoff.
For every chamber term the compact degrees satisfy $D,C\ge0$.
The shape is at least five in $S$ and seven in $P$.
Lemma~\ref{lem:lithium-integrals-compact} applies to every grouped
vertex. The six constant integers used in these executions were
compared exactly with \eqref{eq:lithium-integrals-constant-guard}.
Each upper integer equals the lower integer plus one. Thus their
guard is an evaluated identity, not an unevaluated hypothesis.

The full $12$ scalar and $59$ vector rate groups contain $22825$
and $63216$ nonzero compact vertices. There are respectively two
and $37$ nonempty groups. The other ten and $22$ have the exact
interval $[0,0]$. Each reader comparison includes the ordered
vertex arguments, signed coefficient, rational return, $96$-bit
outward group rounding and column denominator. All groups use
denominator precision $18$ except the first vector group, identified
as $P$-000 in the supplementary manifest, whose selected complete
evaluation uses precision $15$. It has exactly the same input
vertices and constants. Only that replacement is included. The
original wider precision-$18$ interval is retained for provenance
and contributes nothing. The compact error proof applies to either
precision. Reconstructing these selected outputs gives intervals
strictly inside both unchanged Single targets. This proves their
actual containments.

\emph{Oddlog.}
On every generated Oddlog key, $\mu\ge2$ in
Lemma~\ref{lem:lithium-integrals-oddlog-tail}. The complete domain
checks give $5\le\kappa\le17$, $M\le7$ and
$d_\sigma,c_\sigma\ge1$ in $S$; in $P$ they give
$7\le\kappa\le13$, $M\le5$ and $d_\sigma,c_\sigma\ge5$.
These are the conditions at the original cutoff $L=8$.
Every selected key contributes its six exact rational tails, with
sign $(-1)^\mu W_c(k)$, before any finite-part cancellation.
The reader independently recomputes all six terms and every tail prefix.

The finite parts have $54$ nonzero scalar groups with $89467$
compact rows and $390$ vector groups with $251107$ rows. Ordered
returns use compact precision $13$ and the same exact constant
tuple as Single. A group with denominator $\Delta>0$ and signed
row coefficient $w$ contributes the integer endpoints
\[
 \left\lfloor2^{200}\frac{\min(wl,wu)}\Delta\right\rfloor,\qquad
 \left\lceil2^{200}\frac{\max(wl,wu)}\Delta\right\rceil
\]
from a returned interval $[l,u]$.
The retained reference compact moments and complete denominator
triangles reconstruct every ordered return. Where a group spans
several invocations, each carried row is compared with its preceding
actual return, the prefix is recomputed from zero, and fresh calls
start at the next row. Induction over those parts proves the
completed sum uses each row once. The complete group manifest,
not a list of whichever outputs exist, fixes final coverage.
The finite intervals plus the exact signed tails lie inside precisely
the unchanged cutoff-$8$ Oddlog targets in both sectors. This
establishes both actual family containments.

\emph{Doublelog.}
The complete selected lists have $50$ scalar and $72$ vector keys.
The structural forms proved above give the ordered quadrature bounds
\[
\begin{array}{c|c|c|c|c}
 &m&X&Y&X+Y+4\ell\\\hline
 S&9\text{--}15&0\text{--}15&0\text{--}12&\le33\\
 P&11&0\text{--}11&0\text{--}8&\le28
\end{array}
\]
with $\ell\le2$ and $\ell\le3$, respectively. For example, if
the even edge is $mb$, its polynomial power $j$ gives
$X=n_m+n_b+5+2k$, $Y=n_b+2+j$. The other two placements give
$X=n_m+n_b+3+j$ and $Y=n_b+2$ or $n_b+2+j$.
Expanding $(1+t^2-2tx)^k$ for $k=2,3$ verifies these bounds
for every angular term. The rates lie in $[234,1116]/180$.
Thus Lemma~\ref{lem:lithium-integrals-doublelog} applies to every
recorded order-ten, endpoint-power-$34$ call. All $122$ complete
input/return records are compared independently. Their exact signed
rational contractions equal the unchanged two Doublelog endpoint
pairs. All other columns are zero.

\emph{Four extra scalar overlaps.}
Every row of Table~\ref{tab:li-s-trial} has a zero pair power.
Multiplication by $s_n$, $3\le n\le6$, preserves this property,
so its angular cutoff is zero. For a pair with radii $R>r$, its
average is
\[
 \sum_{j=0}^{\lfloor(p+1)/2\rfloor}
 \frac{\binom{p+2}{2j+1}}{p+2}R^{p-2j}r^{2j}.
\]
Each edge lowers its larger-radius exponent by at most one. The
radial Jacobian starts it at at least two. All ordered degrees are
therefore nonnegative. An independent gap-coordinate expansion
$x=g_0+g_1+g_2$, $y=g_1+g_2$, $z=g_2$ evaluates these integrals
as finite gamma products, independently of the original evaluator.
Writing $u_S=A_{13}f$ and $s_n=A_{13}g_n$ gives
\[
 \langle p_{S,n},u_S\rangle
 =2C_{S,n}\bigl(\langle g_n,f\rangle
                  -\langle P_{13}g_n,f\rangle\bigr).
\]
All $4\cdot2\cdot100=800$ signed primitive contributions, their
$80$ block sums and prefixes, and the direct-minus-exchange totals
were evaluated and compared as exact rationals. Applying
Lemma~\ref{lem:lithium-integrals-factors}, with reversed endpoints
for the negative totals, gives the four printed physical targets.
This argument establishes those targets directly. It does not
claim containment in narrower saved scalar-overlap dyadics.

The common array identity identifies the physical meaning of every
sum. The kernel lemmas prove the enclosures of every used return,
and the exhaustive independent signed contractions prove the
stated family and extra-moment containments.
\end{proof}

\begin{lemma}[Physical lithium moments]\label{lem:li-physical-moments}
Assume Set-up~\ref{setup:execution} for the specified integral
and endpoint computations. In the fixed-nucleus spin-$1/2$
Hartree model let $H=H_3(3)$ and use the trials and sources of
Definition~\ref{def:li-trials}. Put $h_T=(H+15/2)u_T$
for $T\in\{S,P\}$. Then the trials and sources belong to $D(H)$,
the sources have norm one, and the following twenty intervals
contain the indicated physical moments simultaneously.
Every displayed decimal endpoint is an exact rational.
\begingroup\small
\[
\begin{array}{c|c|c}
\text{moment for trial }T&S&P\\ \hline
\langle u_T,h_T\rangle&
[.633872331183,.633872821461]&[1.058162598846,1.058162654330]\\
\|h_T\|^2&
[.100526504813,.107640835047]&[.381799258126,.383836614894]\\
\|u_T\|^2&
[25.861096622130,25.861096625746]&[8.460567526663,8.460567530500]\\
\langle p_{T,2},u_T\rangle&
[4.369279500354,4.369279500364]&[2.498720734985,2.498720734986]\\
\langle p_{T,2},h_T\rangle&
[.162740889161,.162740939941]&[.467043822051,.467043822207]\\
\langle p_{T,3},u_T\rangle&
[-2.519811045442,-2.519811045409]&[-1.451977551312,-1.451977551218]\\
\langle p_{T,3},h_T\rangle&
[.035012429476,.035012489149]&[.074296072339,.074296072543]\\
\langle p_{S,4},u_S\rangle&
[-.351885973731,-.351885973646]&--\\
\langle p_{S,4},h_S\rangle&
[-.018396812622,-.018396752955]&--\\
\langle p_{S,5},u_S\rangle&
[-.194319038379,-.194319038203]&--\\
\langle p_{S,5},h_S\rangle&
[-.009880462815,-.009880405542]&--\\
\langle p_{S,6},u_S\rangle&
[-.129253886225,-.129253885908]&--\\
\langle p_{S,6},h_S\rangle&
[-.006454330094,-.006454274838]&--
\end{array}
\]
\endgroup
\end{lemma}

\begin{proof}
Lemma~\ref{lem:lithium-integrals-expansion} proves the domain,
source-norm and spin assertions and identifies the physical moments
with the signed scalar sums. We assemble the thirteen base columns,
then use the seven extra-moment certificates.

\emph{Base columns.}
Use the raw family intervals $J_{T,c,E}$ and physical Splain
intervals $R_c$ of Proposition~\ref{prop:lithium-integrals-certificates}.
Let $G_{T,c}=A$ for the first two base columns and
$G_{T,c}=F_{T,n}$ for a column involving source $n$.
The physical base moments are enclosed by
\begin{align*}
 Y_{P,c}&=B_{G_{P,c}}
   (J_{P,c,\mathrm{Pol}}+J_{P,c,\mathrm{Fin}}
    +J_{P,c,\mathrm{Odd}}+J_{P,c,\mathrm{Log}}
    +J_{P,c,\mathrm{Single}}+J_{P,c,\mathrm{Oddlog}}
    +J_{P,c,\mathrm{Doublelog}}),\\
 Y_{S,c}&=R_c+B_{G_{S,c}}
   (J_{S,c,\mathrm{Log}}+J_{S,c,\mathrm{Single}}
    +J_{S,c,\mathrm{Oddlog}}+J_{S,c,\mathrm{Doublelog}}).
\end{align*}
Here $B_G$ is the factor interval of
Lemma~\ref{lem:lithium-integrals-factors}, and products use all four
corners. Each raw family already includes its column denominator.
The factor $(4\pi)^3C_c$ enters once. In particular, $R_c$ is added
in physical units, with its sign retained.

The independent endpoint reader consumes every column of these
complete family tables. It checks that every inverse-square family
is zero outside the residual-norm column, that the Oddlog cutoff is
$8$, and that the headers and family sizes agree with the complete
arrays. It decodes numeric endpoints as exact nearest-even dyadics:
for a positive rational token $x$, choose the integer $e$ with
$2^e\le x<2^{e+1}$ and round $x/2^{\max(e-52,-1074)}$ to the
nearest integer, taking the even integer at a tie. Negative tokens
are decoded by symmetry. Rational strings remain rational.
Thus no decimal-to-binary uncertainty is omitted.

The retained exact calculation forms all thirteen $Y_{T,c}$ and
compares their endpoints by integer cross multiplication. For
$a/b\ge c/d$, with $b,d>0$, the tested integer is $ad-bc$.
It also reproduces exactly all thirteen final endpoint pairs of
the original collector, using that collector's rational factor
intervals. The tighter factor intervals used here give the first
thirteen outer boxes in Table~\ref{tab:lithium-integrals-assembly}.
Each box is formed as
$[\lfloor10^{14}\inf Y_{T,c}\rfloor,
\lceil10^{14}\sup Y_{T,c}\rceil]/10^{14}$.
The complete exact endpoint and cross-product data are in the
accompanying supplement. Every displayed outer box is contained
in its claimed interval.

\begin{table}[htbp]
\centering\small
\caption{Exact outer numerators for the final endpoint comparisons;
both columns are divided by $10^{14}$.}
\label{tab:lithium-integrals-assembly}
\begin{tabular}{l r r}
\hline
Moment&Lower numerator&Upper numerator\\\hline
$\langle u_S,h_S\rangle$&63387233118325&63387282146003\\
$\|h_S\|^2$&10052650481383&10764083504630\\
$\langle p_{S,2},u_S\rangle$&436927950035475&436927950036339\\
$\langle p_{S,2},h_S\rangle$&16274088916118&16274093994007\\
$\langle p_{S,3},h_S\rangle$&3501242947645&3501248914850\\
$\langle p_{S,4},h_S\rangle$&-1839681262135&-1839675295536\\
$\langle p_{S,5},h_S\rangle$&-988046281448&-988040554218\\
$\langle p_{S,6},h_S\rangle$&-645433009391&-645427483813\\
$\langle u_P,h_P\rangle$&105816259884616&105816265432999\\
$\|h_P\|^2$&38179925812656&38383661489309\\
$\langle p_{P,2},u_P\rangle$&249872073498570&249872073498571\\
$\langle p_{P,2},h_P\rangle$&46704382205165&46704382220679\\
$\langle p_{P,3},h_P\rangle$&7429607233961&7429607254257\\\hline
$\|u_S\|^2$&2586109662213062&2586109662574588\\
$\|u_P\|^2$&846056752666313&846056753049941\\
$\langle p_{P,3},u_P\rangle$&-145197755131125&-145197755121823\\\hline
\end{tabular}
\end{table}

\emph{Extra columns.}
The three native extra moments are contained in their saved physical
intervals by Proposition~\ref{prop:lithium-integrals-certificates}.
The endpoint reader encloses those dyadics by the last three rows of
Table~\ref{tab:lithium-integrals-assembly} and checks their containment
in the printed targets by the same cross products. The four remaining
scalar overlaps, with $3\le n\le6$, are enclosed directly in their
printed targets by the exact gap-coordinate calculation in that
proposition. This uses no stronger claim about the saved overlap
dyadics. The thirteen base enclosures, the three native extra
enclosures and the four exact overlap enclosures therefore give all
twenty physical moments of the same fixed trials and sources.
\end{proof}