==== TABLE OF CONTENTS ====
A-sharp-Fourier-certificate-for-planar-circle-packing//W_minus.tex
A-sharp-Fourier-certificate-for-planar-circle-packing//W_plus.tex
A-sharp-Fourier-certificate-for-planar-circle-packing//certificates.tex
A-sharp-Fourier-certificate-for-planar-circle-packing//coefficients.tex
A-sharp-Fourier-certificate-for-planar-circle-packing//fourier.tex
A-sharp-Fourier-certificate-for-planar-circle-packing//half-gaps.tex
A-sharp-Fourier-certificate-for-planar-circle-packing//independent-certification.tex
A-sharp-Fourier-certificate-for-planar-circle-packing//interpolation.tex
A-sharp-Fourier-certificate-for-planar-circle-packing//introduction.tex
A-sharp-Fourier-certificate-for-planar-circle-packing//main.tex
A-sharp-Fourier-certificate-for-planar-circle-packing//normalization.tex
A-sharp-Fourier-certificate-for-planar-circle-packing//rational-check.tex
A-sharp-Fourier-certificate-for-planar-circle-packing//signs.tex
A-sharp-Fourier-certificate-for-planar-circle-packing//references.bib
===========================
==== BEGIN A-sharp-Fourier-certificate-for-planar-circle-packing//W_minus.tex ====
13265 & -1358 & -385 & -381 & 815 & -164 \\
616 & 10426 & 686 & -484 & 182 & 775 \\
-146 & 4 & 9842 & 215 & 267 & 179 \\
-2761 & 2483 & -958 & 5174 & -328 & 206 \\
-615 & -861 & 730 & -1388 & 5619 & -4370 \\
874 & 820 & 1148 & -431 & 905 & 11927 \\
==== END A-sharp-Fourier-certificate-for-planar-circle-packing//W_minus.tex ====
==== BEGIN A-sharp-Fourier-certificate-for-planar-circle-packing//W_plus.tex ====
4159 & 4452 & -2003 & -10467 & -2421 & -1683 \\
2714 & 6573 & 1689 & 8258 & 104 & -335 \\
-4675 & 4785 & 7003 & -13100 & -1355 & -1163 \\
36853 & -37142 & 27514 & 105973 & 1100 & 2012 \\
57838 & -56941 & 36296 & 158623 & 29097 & 17293 \\
-4705 & 3589 & -3311 & -11932 & -2562 & 6919 \\
==== END A-sharp-Fourier-certificate-for-planar-circle-packing//W_plus.tex ====
==== BEGIN A-sharp-Fourier-certificate-for-planar-circle-packing//certificates.tex ====
\section{Finite certificates}\label{sec:certificates}
This appendix completes the finite evidence used in
Lemma~\ref{interpolation:finite} and Lemma~\ref{signs:bernstein}.
We give the exact inputs, prove a uniform error bound for the finite
quadrature sums, and state the interval comparisons for the matrices,
residuals, and Bernstein coefficients. The final subsection records the
scalar substitutions used for the infinite tails and continuum signs.
All terminating decimals denote exact rational numbers.
The primary calculation uses interval evaluation of the $N=256$ quadrature
sums, enlarged by the analytic moment remainder proved below. It therefore
reduces the required numerical evidence to bounded finite sums, without
assuming convergence of a numerical integration or optimization procedure.
The infinite inverse and the signs between and beyond interpolation nodes
remain consequences of the analytic arguments in the main text.
Appendix~\ref{ind:section} gives the rational alternatives and their
separately stated scope.
\subsection{Exact input lists and matrices}\label{cert:inputs}
Table~\ref{tab:coefficients} gives the coefficients at every node of
$I_f=\mathcal N\cap[1,100]$. Its four coefficient columns are multiplied
by $10^{10}$: an integer $a$ in the table means the coefficient
$a\,10^{-10}$. At all nodes $n>100$ both approximate lists are zero.
The data at zero and the constant terms are, exactly,
\begin{equation}\label{cert:fixed-inputs}
(\widetilde c_{1,0},\widetilde d_{1,0},C_1)
=(1,.44,-.013),\qquad
(\widetilde c_{2,0},\widetilde d_{2,0},C_2)
=(0,-.368,.017).
\end{equation}
The zero-node data are included in all evaluations of the tabulated functions.
They are excluded when forming a column of the operator $S$.
Direct integer addition gives
\begin{equation}\label{cert:table-norms}
\begin{aligned}
\sum_{n\ge0}(|\widetilde c_{1,n}|+|\widetilde d_{1,n}|)
&=2.2516666247<2.29,\\
\sum_{n\ge0}(|\widetilde c_{2,n}|+|\widetilde d_{2,n}|)
&=.8927477314<2.29,\\
\sum_{n>40}(|\widetilde c_{1,n}|+|\widetilde d_{1,n}|)
&=.0000197471<.00002,\\
\sum_{n>40}(|\widetilde c_{2,n}|+|\widetilde d_{2,n}|)
&=.0000161623<.00002.
\end{aligned}
\end{equation}
Here and below an index $n\ge0$ in a coefficient sum ranges over
$\{0\}\cup\mathcal N$.
\begingroup
\small
\setlength{\tabcolsep}{7pt}
\begin{longtable}{r r r r r}
\caption{Rational coefficient approximations. Each of the four coefficient columns is
scaled by $10^{10}$.}\label{tab:coefficients}\\
\toprule
$n$ & $10^{10}\widetilde c_{1,n}$ & $10^{10}\widetilde d_{1,n}$
& $10^{10}\widetilde c_{2,n}$ & $10^{10}\widetilde d_{2,n}$\\
\midrule
\endfirsthead
\multicolumn{5}{c}{Table~\thetable\ (continued)}\\
\toprule
$n$ & $10^{10}\widetilde c_{1,n}$ & $10^{10}\widetilde d_{1,n}$
& $10^{10}\widetilde c_{2,n}$ & $10^{10}\widetilde d_{2,n}$\\
\midrule
\endhead
\midrule
\multicolumn{5}{r}{Continued on the next page}\\
\endfoot
\bottomrule
\endlastfoot
\input{data/coefficients.tex}
\end{longtable}
\endgroup
Write $J=\{1,3,4\}$ and $E=I_f\setminus J$. Matrix rows and columns on a
node set always list all $c$ entries first and all $d$ entries second, with
nodes increasing within each group. In particular the order on $J$ is
$(c_1,c_3,c_4,d_1,d_3,d_4)$. The two rational matrices used to certify the
small inverses are
\begin{equation}\label{cert:Wplus}
W_+=10^{-4}\begin{pmatrix}
\input{data/W_plus.tex}
\end{pmatrix},
\end{equation}
\begin{equation}\label{cert:Wminus}
W_-=10^{-4}\begin{pmatrix}
\input{data/W_minus.tex}
\end{pmatrix}.
\end{equation}
These matrices are exact inputs. The $W$-based inverse argument below
uses their products and column sums; it does not require a computed
matrix inverse. We describe the supplied program's additional inverse
checks in Section~\ref{cert:evaluation}.
\subsection{Moments and the integration rule}\label{cert:quadrature}
For a finite input $p=PR$ supported on nodes at most $100$, put
\begin{equation}\label{cert:moments}
U_{p,l}(m,y)=\frac{y^l p^{(l)}(m)}{l!},\qquad
X_{p,l}(m,y)=\frac{y^l K_p^{(l)}(m)}{l!}.
\end{equation}
The domain for the quadrature error estimate is stated explicitly in
Lemma~\ref{cert:moment-error}.
We use the first Fej\'er quadrature rule, whose cosine-sum weights are
described in \cite[Section~2, formula~(2.3)]{Waldvogel2006}. Our weights
are half the mass-two weights there, normalized to average over an interval;
for even $N=256$ the terminal cosine term in that formula vanishes.
We prove the exactness, positivity and error bound needed here.
Set $N=256$, and for $0\le v<N$ define
\begin{equation}\label{cert:rule}
\begin{aligned}
\alpha_v&=\frac{v+1/2}{N},\qquad
o_v=\frac1N\left(1-2\sum_{a=1}^{127}
\frac{\cos(2a\pi\alpha_v)}{4a^2-1}\right),\\
\tau_{jv}&=\frac{2j-1+\cos(\pi\alpha_v)}{12}
\qquad(1\le j\le6).
\end{aligned}
\end{equation}
Integration over $(t_{j-1},t_j)$ is replaced by evaluation at $\tau_{jv}$
with weights $o_v/6$. To check exactness, expand a polynomial in the
Chebyshev basis $T_k(\cos\theta)=\cos(k\theta)$. Discrete cosine
orthogonality at $\theta_v=\pi\alpha_v$ recovers its coefficients through
degree $255$. The integral averages of this basis on $[-1,1]$ are zero
at odd $k$ and $1/(1-k^2)$ at even $k$. These are exactly the coefficients
in~\eqref{cert:rule}, so the rule is exact through degree $255$.
Moreover,
\[
o_v\ge\frac1{256}\left(1-\sum_{a=1}^{127}\frac2{4a^2-1}\right)
=\frac1{256\cdot255}>0,\qquad \sum_v o_v=1.
\]
Thus the six pieces together have total weight one.
For clarity, the complete density on piece $j$ used in this appendix is
\begin{equation}\label{cert:density}
f_j(t)=2\pi\sum_n e^{-i\pi tn}
\left\{-\pi c_n\sum_{l=j}^6(t_l-t)P_l e^{i\pi t_ln}
+i d_n\sum_{l=j}^6 P_l e^{i\pi t_ln}\right\}.
\end{equation}
It includes node $0$ when that coefficient is part of the input. The
folded atoms have masses $e_jCP_j$ at $t_j$, where $e_0=1$ and
$e_j=2$ for $j>0$. The factors multiplying this density and these masses
for $U$ and $X$, respectively, are
\begin{equation}\label{cert:analytic-factors}
e^{i\pi tm}\frac{(i\pi ty)^l}{l!},\qquad
\lambda(t)e^{i\pi z(t)m}\frac{(i\pi z(t)y)^l}{l!}.
\end{equation}
The real part is taken after summation.
\begin{lemma}[Quadrature moment error]\label{cert:moment-error}
Let $p=PR$ have a real coefficient list supported in $\{0\}\cup I_f$,
with $\|(c,d)\|_1\le3$, and let $C\in\RR$ satisfy $|C|\le1$.
Let $m,y\in\RR$ and $l\in\{0,\ldots,30\}$ satisfy
$0\le m\le100$ and $|y|\le3/2$.
The quadrature rule~\eqref{cert:rule}, with exact evaluations at its nodes
and at the atoms, has error less than $10^{-22}$ in each moment
$U_{p,l}(m,y)$ and $X_{p,l}(m,y)$.
\end{lemma}
This domain contains all basis inputs and the two tabulated inputs used below.
\begin{proof}
Put $\rho=1/6$, $c_j=(2j-1)/12$, $a_j=c_j+ih$, and
$\Delta_j=|a_j|^2-\rho^2$. These last quantities increase from
$\Delta_1=167/1200$ to $\Delta_6=389/400$, so the disk
$|t-c_j|\le\rho$ does not contain the pole $-ih$. All factors in
\eqref{cert:density} and~\eqref{cert:analytic-factors} are analytic there.
The reciprocal-disk calculation is exact: for $\zeta=1/(t+ih)$,
\[
|1-a_j\zeta|^2-\rho^2|\zeta|^2
=\Delta_j\left|\zeta-\frac{\overline a_j}{\Delta_j}\right|^2
-\frac{\rho^2}{\Delta_j}.
\]
Thus inversion maps the disk for $t+ih$ onto the disk with center
$(c_j-ih)/\Delta_j$ and radius $\rho/\Delta_j$. In particular,
\[
|\zeta|\le\frac{|a_j|+\rho}{\Delta_j}
=\frac1{|a_j|-\rho}.
\]
The right side decreases with $c_j$, and $|a_1|=\sqrt{601}/60$.
Consequently the following are analytic bounds on every disk:
\[
\begin{aligned}
|t|&\le\frac{13}{12}<1.1,\\
|\lambda(t)|&\le\frac{120}{\sqrt3(\sqrt{601}-10)}<4.774<5,\\
|z(t)|&\le\frac{80}{\sqrt{601}-10}+\frac25<5.912<6.2.
\end{aligned}
\]
For the two finer comparisons, it suffices to square the rational lower
bounds $1.732^2=2.999824<3$ and $24.515^2=600.985225<601$, and then use
\[
\begin{gathered}
4.774(1.732)(24.515-10)=120.01826452>120,\\
(5.912-.4)(24.515-10)=80.00668>80.
\end{gathered}
\]
The image disk for $z=-B\zeta-ih$ has center
$-B(c_j-ih)/\Delta_j-ih$ and radius $B\rho/\Delta_j$. Its lowest
imaginary part is therefore $B(h-\rho)/\Delta_j-h$. Since $h>\rho$,
this decreases with $\Delta_j$, giving
\begin{equation}\label{cert:disk-imaginary}
\operatorname{Im}z(t)\ge
B\frac{h-1/6}{c_j^2+h^2-(1/6)^2}-h
\ge-\frac{1402}{17505}>-.081.
\end{equation}
Since $\sum_{l=1}^6|P_l|=10$ and $|t_l-t|\le13/12$ for $l\ge j$,
the analytic continuation of~\eqref{cert:density} on this disk is bounded
by
\[
3\max\{2\pi^2\cdot10\cdot13/12,\,2\pi\cdot10\}
e^{100\pi/6}<1000e^{100\pi/6}.
\]
Use $|a^l/l!|\le e^{|a|}$ for the divided powers. Including all factors,
the direct and transformed integrands are bounded, respectively, by
\[
\begin{aligned}
1000e^{200\pi/6+1.5\pi\cdot1.1}&<5.38\cdot10^{50},\\
5000e^{100\pi/6+100\pi\cdot.081+1.5\pi\cdot6.2}
&<1.51\cdot10^{50}.
\end{aligned}
\]
We may therefore use the common bound $M=10^{53}$.
The real interval has half-length $1/12$, half the disk radius.
Cauchy's coefficient estimate bounds the tail of the Taylor polynomial
of degree $255$ by $M\sum_{k=256}^{\infty}2^{-k}=M2^{-255}$.
Exactness, positivity and total mass one for the integral and quadrature
rule give discrepancy at most
\begin{equation}\label{cert:quadrature-error}
2^{-254}10^{53}<3.455\cdot10^{-24}<10^{-22}.
\end{equation}
Atoms are evaluated separately and incur no quadrature error.
\end{proof}
\subsection{Interval evaluation}\label{cert:evaluation}
For a finite input $p$ in the real finite-input domain of
Lemma~\ref{cert:moment-error} and a node row $m\in I_f$, form the two real numbers
\begin{equation}\label{cert:pair}
\mathcal S_m(p)=
\left(-\frac{X_{p,0}(m,1)}{Q_m},\;
\frac{D_mX_{p,0}(m,1)-X_{p,1}(m,1)}{Q_m}\right).
\end{equation}
With input $c_n=1$ or $d_n=1$ and all other data zero, these are the
entries of the corresponding column of $S$ at row node $m$.
Let $D=S_{J,J}$ and $U=S_{J,E}$.
For the Bernstein certificates, use exactly the triples $(i,m,y)$ and
orders $\nu$ specified in Lemma~\ref{signs:bernstein}; in particular $y$
is signed and no leftward interval at $m=1$ is used for index $1$.
Express its power coefficients in terms of the moments as
\begin{equation}\label{cert:power-from-moments}
T_r=\sigma_i y^{-\nu}
\left(U_{\widetilde p_i,r+\nu}(m,y)
+X_{\widetilde p_{3-i},r+\nu}(m,y)\right).
\end{equation}
The signed scaling is exact. Definition~\eqref{cert:moments} gives
$U_{p,l}(m,y)=y^lU_{p,l}(m,1)$ and the same identity for $X$.
Therefore, for positive or negative $y$,
\[
T_r=\sigma_i y^r\left(U_{\widetilde p_i,r+\nu}(m,1)
+X_{\widetilde p_{3-i},r+\nu}(m,1)\right).
\]
This is the identity used when the implementation caches unscaled moments
at $y=1$; it does not discard the orientation of a half-gap.
Equations~\eqref{cert:rule}, \eqref{cert:density},
\eqref{cert:analytic-factors}, \eqref{cert:pair}, and
\eqref{cert:power-from-moments} already specify all
sums in the finite certificates below. The following pseudocode records their
algebraic organization. A measure is a list of pairs $(t,g)$ in the order
used by the supplied interval implementation: first the $6\cdot256$
quadrature pairs $(\tau_{jv},o_v f_j(\tau_{jv})/6)$, with
$j=1,\ldots,6$ outermost and $v=0,\ldots,255$ innermost, and then the seven
atom pairs $(t_j,e_jCP_j)$ for $j=0,\ldots,6$. The atom masses are zero for
a basis input. The pseudocode is not a byte-for-byte execution trace.
The supplied source fixes the detailed arithmetic and summation order for
an actual interval run, including the factoring and caching described below.
\noindent
\begin{minipage}{\linewidth}
\begin{lstlisting}[basicstyle=\ttfamily\footnotesize]
Moments(measure, m, y, L):
U[0..L] = 0; X[0..L] = 0
for (t,g) in measure:
lam = i/(b*(t+i*h))
z = -B/(t+i*h)-i*h
u = g*exp(i*pi*t*m)
v = g*lam*exp(i*pi*z*m)
for r = 0..L:
U[r] += real(u); X[r] += real(v)
u *= (i*pi*t*y)/(r+1)
v *= (i*pi*z*y)/(r+1)
return U,X
Pair(measure,m):
U,X = Moments(measure,m,1,1)
return -X[0]/Q[m], (D[m]*X[0]-X[1])/Q[m]
for basis column c[n] or d[n], n in I_f:
measure = Measure(basis column; all other data zero)
column = stack Pair(measure,m), m in J
if n in J:
also form rows m in E intersect [7,16]
test column sums of W[eta] and the defect/product matrices
for i in {1,2}, m in I_f:
c,d = Pair(Measure(full table input 3-i),m)
test abs(table_c[i,m]-c)
test abs(table_d[i,m]+indicator(i==1 and m==1)-d)
for (i,m,y,nu) in the prescribed half-gap cases:
U,X = Moments for the two full table inputs, through L=30
T[r] = sigma[i]*(U[i,r+nu]+X[3-i,r+nu])/y**nu
for k = 0..28:
Bk = T[0]; a = 1
for r = 1..k:
a *= (k+1-r)/(29-r)
Bk += a*T[r]
test the appropriate lower bound for Bk
\end{lstlisting}
\end{minipage}
Here \texttt{stack} means all first entries followed by all second entries.
The recurrence for $a$ equals $\binom{k}{r}/\binom{28}{r}$.
\texttt{Measure} denotes this finite list and introduces no additional
approximation. The fixed data at zero are used only for the two full
table inputs. Reusing a measure means reusing its whole finite list.
For the accompanying certified evaluation, every real or complex quantity
in these sums is enclosed using interval arithmetic. The pair calculation
forms the moments at $y=1$. For the Bernstein calculation, the
implementation also forms and caches the unscaled moments at $y=1$, then
uses the signed $y^r$ identity following~\eqref{cert:power-from-moments}.
Each unscaled real moment enclosure is enlarged by
$[-10^{-22},10^{-22}]$, as justified by Lemma~\ref{cert:moment-error}.
Consequently the analytic-remainder contribution to the enclosure radius
of a Bernstein coefficient is bounded by
\[
\begin{aligned}
2\cdot10^{-22}\sum_{r=0}^k
\frac{\binom{k}{r}}{\binom{28}{r}}|y|^r
&\le2\cdot10^{-22}\sum_{r=0}^{28}(3/2)^r\\
&=4\cdot10^{-22}\bigl((3/2)^{29}-1\bigr)<5.2\cdot10^{-17}.
\end{aligned}
\]
The last inequality is an exact rational comparison. It is not correct
to regard $10^{-22}$ as an unamplified per-Bernstein error.
Arb rounding and dependency effects can further widen the balls; all of
them are propagated by the subsequent products, absolute values, column
sums, and Bernstein combinations. An upper bound is accepted only when
the complete enclosure lies strictly below the asserted rational bound,
and a lower bound only when it lies strictly above. This procedure uses
enclosures of the exact integrals, not a quadrature error inferred from
observed numerical agreement.
The ball-arithmetic inclusion methodology is described in
\cite{Johansson2017}. The supplied implementation uses Arb through
python-flint at $192$ bits of precision; the citation does not by itself
certify this local environment, its call domains, or a run.
The entry point \path{verification/verify_all.py} runs two verifiers in
the same directory, \path{interval_certificate.py} and
\path{interval_bounds.py}. Their output records, \path{finite.json} and
\path{scalar.json}, and the aggregate record \path{summary.json} are
placed in its \path{results/} subdirectory.
The first verifier checks the exact bounds in~\eqref{cert:matrix-table},
\eqref{cert:exterior-block}, \eqref{cert:residual-bounds}, and
\eqref{cert:bernstein-table}, as well as all 2436 individual Bernstein
inequalities in Lemma~\ref{signs:bernstein}. It also asserts
four finite brackets: the two defect brackets, the exterior bracket, and
the bracket for $(2,16,-1/2,28)$. All its comparisons and extrema use
complete enclosures.
The mathematical $W$ argument in the next subsection uses only products
and column sums. The supplied program also calls a certified inverse
for $I-\eta D$ and tests its inverse and cross-product norms as
additional checks; these calls are not needed for the $W$ deduction.
\subsection{Matrix and residual certificates}\label{cert:matrix-bounds}
The following are strict upper bounds
for the indicated column norms. The last three columns use only
columns whose node lies in the stated range. They refer to the exact
operator entries in~\eqref{cert:pair}.
\begin{equation}\label{cert:matrix-table}
\begin{array}{c|cc|ccc}
\eta&
\|W_\eta\|&
\|I-W_\eta(I-\eta D)\|&
\|W_\eta U_{7:9}\|&
\|W_\eta U_{12:16}\|&
\|W_\eta U_{19:100}\|\\ \hline
+1&30.84&.001&3.54&.44&.25\\
-1&1.83&.001&.73&.32&.17
\end{array}
\end{equation}
The unmultiplied exterior block also satisfies
\begin{equation}\label{cert:exterior-block}
\|S_{E\cap[7,16],J}\|<.088.
\end{equation}
The certified enclosure calculation described in
Section~\ref{cert:evaluation} proves these inequalities.
For example, the two defect norms are enclosed respectively in
\[
\begin{gathered}
[.00025422624910652544,\ .00025422624910652546],\\
[.00042288165102553355,\ .00042288165102553357],
\end{gathered}
\]
and the exterior block norm is enclosed in
$[.08688748182289641591,\ .08688748182289641592]$.
These intervals are stated with rational endpoints rounded outwards.
The three printed product ranges partition $E$, so $\|W_\eta U\|$ is the
maximum of the three grouped product norms. The exact bounds in
\eqref{cert:matrix-table} and~\eqref{cert:exterior-block} give the primary
matrix conclusion directly. Put $A_\eta=I-\eta D$.
Since $\|I-W_\eta A_\eta\|<.001$, the finite square
matrix $W_\eta A_\eta$ is invertible by a geometric series, hence so
is $A_\eta$, and
\[
A_\eta^{-1}=(W_\eta A_\eta)^{-1}W_\eta.
\]
It follows from~\eqref{cert:matrix-table} that
\[
\|A_\eta^{-1}\|<\frac{30.84}{1-.001}<32,\qquad
\|A_\eta^{-1}U\|<\frac{3.54}{1-.001}<3.7.
\]
Equation~\eqref{cert:exterior-block} gives the remaining bound $<.09$
in Lemma~\ref{interpolation:finite}.
For the two full tabulated inputs, including~\eqref{cert:fixed-inputs},
the residual certificate is
\begin{equation}\label{cert:residual-bounds}
\begin{aligned}
\left|\widetilde c_{i,m}+\frac{X_{\widetilde p_{3-i},0}(m,1)}{Q_m}\right|
&<3\cdot10^{-10},\\
\left|\widetilde d_{i,m}+\mathbf1_{\{i=1,\ m=1\}}
+\frac{X_{\widetilde p_{3-i},1}(m,1)
-D_mX_{\widetilde p_{3-i},0}(m,1)}{Q_m}\right|
&<3\cdot10^{-10}
\qquad(i=1,2,\ m\in I_f).
\end{aligned}
\end{equation}
The primary evaluation in Section~\ref{cert:evaluation} proves
\eqref{cert:residual-bounds} for these exact residuals. Since
$3\cdot10^{-10}<10^{-9}$, these inequalities, the preceding primary
matrix conclusion, and~\eqref{cert:table-norms} prove every finite
assertion of Lemma~\ref{interpolation:finite}.
\subsection{Bernstein coefficient certificates}\label{cert:bernstein-bounds}
\begin{samepage}
Let $B_{i,m,y,k}$ denote the sum $B_k$ in
\eqref{signs:bernstein-bound}. Each entry below is a strict lower bound
for every such sum in its row, for all $0\le k\le28$ and both half-gap
directions when applicable.
\begin{equation}\label{cert:bernstein-table}
\begin{array}{c|cc}
m& i=1&i=2\\ \hline
0&\text{---}&.762\\
1&.314&.191\\
3,4&.028&.024\\
7,9&.114&.112\\
12,13,15,16&.0105&.0094\\
19,21&.123&.113\\
24,25,27,28&.0111&.0113\\
31,33&.121&.126\\
36,37,39,40&.0107&.0115
\end{array}
\end{equation}
\end{samepage}
These bounds concern the exact derivatives of the tabulated functions,
not derivatives rounded to satisfy interpolation identities.
For example, the coefficient at $(i,m,y,k)=(2,16,-1/2,28)$ is enclosed in
\[
[.00951875955680,\ .00951875955684].
\]
The primary evaluation in Section~\ref{cert:evaluation} proves the grouped
bounds~\eqref{cert:bernstein-table} for the exact derivatives of the tabulated
functions. Each printed lower bound is stronger than its corresponding
threshold in~\eqref{signs:bernstein-bound}. Therefore these exact grouped
bounds imply Lemma~\ref{signs:bernstein}.
\subsection{Scalar envelope and sign bounds}\label{cert:scalar-bounds}
The remaining finite substitutions involve elementary real expressions.
For the endpoint and density envelopes in
\eqref{interpolation:piece-envelopes}, it suffices to test the finitely
many endpoints and residues specified there. The endpoint formulas
for $|\lambda|$, $|z|$, and $\operatorname{Im}z$ use rational numbers and
square roots. The density expressions additionally use sine and cosine
at rational multiples of $\pi$. The endpoint tests are certified by
the same elementary enclosures used above; their monotonicity and
convexity extensions to whole intervals were proved in
Section~\ref{sec:interpolation}.
For the row envelope in~\eqref{interpolation:row-envelope}, certified
outward bounds are
\begin{equation}\label{cert:row-substitutions}
\begin{array}{c|c}
\text{row cut }x&\text{upper bound for the finite envelope sum}\\ \hline
0&27.618723\\
4&.269094\\
16&.014122\\
100&4.617010\cdot10^{-10}
\end{array}
\end{equation}
The atom-column analogue at cut $100$ is less than
$3.131186\cdot10^{-8}$. These bounds imply all the row-tail constants
used in Section~\ref{sec:interpolation}.
For the sign proof, the following outward bounds record the relevant
substitutions. The expressions themselves, and the reason they bound
whole intervals or tails, are given at the cited equations.
\begin{center}
\begin{tabular}{p{.65\linewidth}r}
\toprule
Expression & Certified bound\\
\midrule
\eqref{signs:rational-lower-sum}, index $1$ & $>.01206834$\\
\eqref{signs:rational-lower-sum}, index $2$ & $>.01656468$\\
$\sum_{j=1}^6 M_j(\pi t_j)^2/6$ & $<76.043$\\
Right side of \eqref{signs:tail-second-derivative}
before its final comparison & $<.005303$\\
\eqref{signs:perturbation-envelope} divided by $\nu!$,
$(s_*,\nu)=(0,0)$ & $<.003189$\\
Same expression, $(s_*,\nu)=(1/2,1)$ & $<.003736$\\
Same expression, $(s_*,\nu)=(1/2,2)$ & $<.009146$\\
Same expression, $(s_*,\nu)=(2,2)$ & $<.001796$\\
$2.29\sum_{j=1}^6M_j/6+.017\sum_{j=0}^6m_j^*$ & $<72.669$\\
\bottomrule
\end{tabular}
\end{center}
In the fourth row the expression means the explicit sum with
$.00005\cdot77$, not the derivative supremum to its left.
For the five cases of~\eqref{signs:taylor-representatives}, the
corresponding bounds for $\mathcal E(m,Y,\nu)$ in
\eqref{signs:taylor-bound} are
\begin{equation}\label{cert:taylor-substitutions}
\begin{array}{c|c}
(m,Y,\nu)&\text{upper bound}\\ \hline
(0,1/2,0)&1.578\cdot10^{-10}\\
(1,1,1)&2.007\cdot10^{-5}\\
(1,1,2)&5.735\cdot10^{-6}\\
(3,1,2)&9.581\cdot10^{-11}\\
(4,3/2,2)&6.133\cdot10^{-7}
\end{array}
\end{equation}
The first two rational-part bounds can be obtained by exact fraction
addition alone. All the other displayed scalar bounds follow by
substitution and outward enclosure of the explicit exponential,
trigonometric and rational expressions. The scalar verifier records
22 named strict comparisons: the four row substitutions, the atom tail,
the nine sign-table rows, the five Taylor representatives, and the three
disk-integrand and quadrature comparisons. It also checks the endpoint
and residue envelopes, the coarser row-tail comparisons, the finite-inverse
substitution giving a bound below $90$,
and the distance substitution using $91$ and $2540$. The derivation of
these whole-space coefficients from the Schur estimates remains the
analytic argument in Section~\ref{sec:interpolation}.
Its disk assertions use the coarser bounds $5$ and $6.2$ for
$|\lambda|$ and $|z|$; the finer $4.774$ and $5.912$ bounds were proved
analytically above and are not literal scalar-program assertions.
At the six gap midpoints the program evaluates the product defining $P$
and checks that each quotient is greater than $.68$; it does not assert
the symbolic radical identities for those values. No sampling of an
unknown function between interpolation nodes is involved. The
logarithmic concavity argument and exact midpoint values in
Lemma~\ref{signs:P-lower-bound} supply the global lower bound for $P$.
==== END A-sharp-Fourier-certificate-for-planar-circle-packing//certificates.tex ====
==== BEGIN A-sharp-Fourier-certificate-for-planar-circle-packing//coefficients.tex ====
1 & -670689128 & -6398628799 & -250018302 & 4107732082 \\
3 & 93524468 & 564897270 & -186215974 & -65204793 \\
4 & -64146977 & 121133118 & 82882809 & -327822296 \\
7 & -1358193 & -16329465 & 1979670 & 21637038 \\
9 & -2932910 & -940875 & 3636920 & 2145806 \\
12 & 3646952 & -34846939 & -6005382 & 34369042 \\
13 & 3183411 & 63030995 & -574123 & -75237042 \\
15 & 11390977 & 27049982 & -11312377 & -34985096 \\
16 & -5593246 & 13405071 & 5915710 & -10980954 \\
19 & 54905 & -1515007 & -120136 & 1498351 \\
21 & -210208 & -634307 & 177992 & 747946 \\
24 & 976077 & -804949 & -1037705 & 105643 \\
25 & -824472 & 7003719 & 1011732 & -6510015 \\
27 & 611891 & 4536095 & -428528 & -4440927 \\
28 & -472882 & -107613 & 399655 & 403751 \\
31 & 30744 & -112653 & -31773 & 90515 \\
33 & -4974 & -105568 & 1240 & 99081 \\
36 & 128313 & 246018 & -114876 & -267388 \\
37 & -176685 & 505187 & 166591 & -396209 \\
39 & -21638 & 472037 & 32128 & -402283 \\
40 & -18461 & -165597 & 10836 & 164974 \\
43 & 4565 & -2296 & -4050 & 547 \\
45 & 1477 & -10734 & -1524 & 8978 \\
48 & 10748 & 50050 & -8611 & -45048 \\
49 & -20506 & 7495 & 17354 & -266 \\
51 & -10361 & 29035 & 9641 & -21534 \\
52 & 2121 & -26523 & -2269 & 23159 \\
55 & 425 & 858 & -344 & -839 \\
57 & 302 & -632 & -266 & 463 \\
60 & 382 & 5658 & -231 & -4688 \\
61 & -1481 & -4320 & 1140 & 4098 \\
63 & -1494 & -352 & 1265 & 647 \\
64 & 558 & -2609 & -491 & 2114 \\
67 & 21 & 169 & -14 & -145 \\
69 & 34 & 13 & -28 & -18 \\
72 & -55 & 394 & 53 & -301 \\
73 & -21 & -805 & 5 & 685 \\
75 & -136 & -379 & 108 & 338 \\
76 & 70 & -144 & -58 & 104 \\
79 & -1 & 19 & 1 & -15 \\
81 & 2 & 9 & -2 & -8 \\
84 & -13 & 3 & 11 & 0 \\
85 & 13 & -87 & -11 & 70 \\
87 & -6 & -60 & 4 & 49 \\
88 & 6 & 5 & -4 & -5 \\
91 & 0 & 1 & 0 & -1 \\
93 & 0 & 1 & 0 & -1 \\
96 & -2 & -4 & 1 & 3 \\
97 & 2 & -6 & -2 & 4 \\
99 & 0 & -6 & 0 & 5 \\
100 & 0 & 2 & 0 & -2 \\
==== END A-sharp-Fourier-certificate-for-planar-circle-packing//coefficients.tex ====
==== BEGIN A-sharp-Fourier-certificate-for-planar-circle-packing//fourier.tex ====
\section{Gaussian cardinal functions}\label{sec:fourier}
The construction uses a family of entire functions whose values and first
derivatives can be prescribed at a periodic set of nodes. Multiplication by
a fixed Gaussian turns these functions into radial Schwartz functions. Their
Fourier transforms have strictly stronger Gaussian decay, which will make the
coupled interpolation equations tractable.
Set
\begin{equation}\label{fourier:parameters}
b=\frac{\sqrt3}{2},\qquad h=\frac25,\qquad B=b^{-2}=\frac43,
\qquad s=b\abs{x}^2.
\end{equation}
The same coordinate, \(s=b\abs{\xi}^2\), is used on the Fourier side.
To see the lattice geometry behind this choice, let \(\Lambda\) have basis
\begin{equation}\label{fourier:lattice}
u=\frac{1}{\sqrt b}(1,0),\qquad
v=\frac{1}{\sqrt b}(1/2,b).
\end{equation}
This oriented basis has determinant one. If \(\mathcal R\) is
counterclockwise rotation through \(90\) degrees, its dual basis is
\(-\mathcal R v,\mathcal R u\). Hence
\(\Lambda^*=\mathcal R\Lambda\): the lattice and its dual have the same
shell radii. For \(a=ju+kv\), with \(j,k\in\ZZ\), the squared radius in
the common coordinate is
\begin{equation}\label{fourier:shell-coordinate}
s_a:=b\abs{a}^2=j^2+jk+k^2.
\end{equation}
We impose interpolation conditions at the positive nodes of a periodic set:
\begin{equation}\label{fourier:nodes}
A=\{0,1,3,4,7,9\},\qquad
\mathcal N=(12\ZZ+A)\cap(0,\infty).
\end{equation}
The set \(\mathcal N\) contains every positive value of \(j^2+jk+k^2\).
Indeed,
\(j^2+jk+k^2\equiv(j-k)^2\pmod3\), so it is never \(2\pmod3\).
Modulo \(4\), it is \(0\) when \(j,k\) are both even and is odd otherwise.
Combining these restrictions gives precisely the six allowed
residues in \(A\). The inclusion is strict: \(15\in\mathcal N\), but
\(j^2+jk+k^2=15\) would imply
\((2j+k)^2\equiv2k^2\pmod5\). Since \(2\) is not a square modulo \(5\),
this forces \(5\mid k\) and then \(5\mid j\), which would make \(25\)
divide \(15\). Thus we require zeros at some nodes that are not lattice
shell coordinates.
This deliberate enlargement uses a periodic residue pattern with an
elementary trigonometric vanishing function. Covering planar shell values by
periodic sets also occurs in Sardari's Fourier interpolation construction
\cite[Section~6.1]{Sardari2021}; the Gaussian transformation and the
coupled inversion needed here will be proved directly.
\subsection{The cardinal family}
The values and first derivatives will be encoded by double-pole and
simple-pole terms multiplied by a common vanishing factor. This is the
same cardinal pattern that occurs in Gaussian Beurling--Selberg
interpolation \cite[formulas~(1.2)--(1.3)]{CarneiroLittmannVaaler2013}.
The nodes and simultaneous Fourier constraints here are different.
For \(-6\le j\le6\), put \(t_j=j/6\), and define
\begin{equation}\label{fourier:polynomial}
P(s)=\prod_{a\in A}\left(2\sin\frac{\pi(s-a)}{12}\right)^2
=\sum_{j=-6}^{6}P_j e^{i\pi t_js}.
\end{equation}
Here \(P_{-j}=\overline{P_j}\), and
\begin{equation}\label{fourier:coefficients}
(P_0,P_1,\ldots,P_6)
=\left(5,-1+2ib,1+2ib,-2,-\frac12+ib,-1-2ib,1\right).
\end{equation}
For example, with \(w=e^{i\pi s/6}\) and \(\rho=e^{i\pi/6}\), each factor is
\[
\left(2\sin\frac{\pi(s-a)}{12}\right)^2
=-\rho^{-a}w^{-1}(w-\rho^a)^2.
\]
The six minus signs cancel, and \(\rho^{-\sum_{a\in A}a}=\rho^{-24}=1\).
Thus the product in \eqref{fourier:polynomial} equals
\[
w^{-6}\prod_{a\in A}(w-\rho^a)^2.
\]
Expanding this Laurent polynomial gives
\eqref{fourier:coefficients}.
Every zero \(n\in12\ZZ+A\) of \(P\) is double. Define \(Q_n,D_n\) by
\begin{equation}\label{fourier:node-expansion}
P(n+u)=Q_nu^2\bigl(1+D_nu+O(u^2)\bigr).
\end{equation}
If \(a\in A\) is the residue of \(n\) modulo \(12\), then
\begin{equation}\label{fourier:node-data}
\begin{split}
Q_n&=\left(\frac\pi6\right)^2
\prod_{\substack{a'\in A\\a'\ne a}}
\left(2\sin\frac{\pi(a-a')}{12}\right)^2,\\
D_n&=\frac\pi6\sum_{\substack{a'\in A\\a'\ne a}}
\cot\frac{\pi(a-a')}{12}.
\end{split}
\end{equation}
Thus \(Q_n,D_n\) are periodic in \(n\), and their exact values are
\[
\begin{array}{c|c|c}
a&Q_n&D_n\\\hline
0&\pi^2/3&-7\pi\sqrt3/18\\
1&(2\pi^2/3)(2-\sqrt3)&\pi(3+\sqrt3)/18\\
3&(2\pi^2/3)(2-\sqrt3)&-\pi(3+\sqrt3)/18\\
4&\pi^2/3&7\pi\sqrt3/18\\
7&(2\pi^2/3)(2+\sqrt3)&-\pi(3-\sqrt3)/18\\
9&(2\pi^2/3)(2+\sqrt3)&\pi(3-\sqrt3)/18
\end{array}
\]
The expansion follows by isolating the vanishing sine factor; logarithmic
differentiation of the remaining factors gives \(D_n\). The table then uses
the elementary trigonometric values at multiples of \(\pi/12\).
In particular,
\begin{equation}\label{fourier:node-bounds}
Q_n>\frac74=1.75,\qquad \abs{D_n}<\frac{53}{25}=2.12.
\end{equation}
To check these inequalities with rational bounds, use
\(157/50<\pi<22/7\) and \(\sqrt3<26/15\). The minimum of \(Q_n\) is
greater than \(197192/112500>7/4\), and the maximum of \(\abs{D_n}\) is
less than \(286/135<53/25\).
Let \(C\in\RR\), and let \((c_n,d_n)_{n\in\{0\}\cup\mathcal N}\) be a
real coefficient list with
\begin{equation}\label{fourier:list-norm}
\norm{(c,d)}_1
:=\sum_{n\in\{0\}\cup\mathcal N}(\abs{c_n}+\abs{d_n})<\infty.
\end{equation}
Away from the nodes, define
\begin{equation}\label{fourier:cardinal}
p(s)=P(s)R(s),\qquad
R(s)=C+\sum_{n\in\{0\}\cup\mathcal N}
\left(\frac{c_n}{(s-n)^2}+\frac{d_n}{s-n}\right).
\end{equation}
We next show that the apparent singularities are removable, with uniform
control of the resulting entire function.
\begin{lemma}[Finite-measure representation]\label{fourier:measure-lemma}
The function in \eqref{fourier:cardinal} extends to an entire function with
the representation
\begin{equation}\label{fourier:measure}
p(s)=\int_{-1}^{1}e^{i\pi ts}\,\dd\mu(t),\qquad s\in\CC,
\end{equation}
where \(\mu\) is a finite complex measure, supported in \([-1,1]\), with
\(\mu(-E)=\overline{\mu(E)}\) for Borel sets \(E\).
Its total variation satisfies
\begin{equation}\label{fourier:variation}
\norm{\mu}_{\mathrm{TV}}
\le25\abs C+\frac{65\pi^2}{18}\sum_n\abs{c_n}
+\frac{32\pi}{3}\sum_n\abs{d_n}.
\end{equation}
The atoms of \(\mu\) have mass \(CP_j\) at \(t_j\), for \(-6\le j\le6\).
For real-axis evaluations, the remaining measure can be folded onto
\([0,1]\): on \(t_{j-1}<t<t_j\), its folded density is
\(\sum_n(c_nu_n(t)+d_nv_n(t))\), where
\begin{equation}\label{fourier:folded}
\begin{split}
u_n(t)&=-2\pi^2e^{-i\pi tn}
\sum_{l=j}^{6}(t_l-t)P_l e^{i\pi t_ln},\\
v_n(t)&=2i\pi e^{-i\pi tn}
\sum_{l=j}^{6}P_l e^{i\pi t_ln}.
\end{split}
\end{equation}
The folded atoms have masses \(CP_0\) at \(0\) and \(2CP_j\) at \(t_j\),
\(1\le j\le6\). Integrate against this folded measure and take the real
part to obtain \(p(s)\) for real \(s\), or any of its real-axis derivatives.
\end{lemma}
\begin{proof}
Fix a node \(n\) and write \(u=s-n\). For real \(\tau\), the identities
\begin{equation}\label{fourier:division-identities}
\begin{split}
\frac{e^{i\pi\tau u}-1}{u}
&=i\pi\int_0^\tau e^{i\pi tu}\,\dd t,\\
\frac{e^{i\pi\tau u}-1-i\pi\tau u}{u^2}
&=-\pi^2\int_0^\tau(\tau-t)e^{i\pi tu}\,\dd t
\end{split}
\end{equation}
hold for complex \(u\), with the values at zero taken by continuity.
When \(\tau<0\), the integrals are oriented integrals. Inserting
\(\tau=t_l\) and multiplying by \(P_le^{i\pi t_ln}\), the subtracted
terms sum to zero because \(P(n)=P'(n)=0\). On the positive interval
\((t_{j-1},t_j)\), only frequencies with \(l\ge j\) contribute. The
resulting density for \(P(s)/(s-n)^2\) is \(u_n/2\), and that for
\(P(s)/(s-n)\) is \(v_n/2\). The
negative-frequency densities are their conjugate reflections, which proves
\eqref{fourier:folded} and accounts for its factors of two.
All phase factors involving \(n\) have modulus one on the real integration
interval. Hence the measures for the two individual basis functions satisfy
\[
\begin{split}
\int_0^1\abs{u_n(t)}\,\dd t
&\le2\pi^2\sum_{l=1}^{6}\abs{P_l}
\int_0^{t_l}(t_l-t)\,\dd t
=\pi^2\sum_{l=1}^{6}\abs{P_l}t_l^2
=\frac{65\pi^2}{18},\\
\int_0^1\abs{v_n(t)}\,\dd t
&\le2\pi\sum_{l=1}^{6}\abs{P_l}t_l
=\frac{32\pi}{3}.
\end{split}
\]
These are also bounds for the total variations of the full measures on
\([-1,1]\). The atomic measure for \(CP\) has total variation
\(\abs C\sum_{j=-6}^{6}\abs{P_j}=25\abs C\). The coefficient
summability therefore gives convergence of the sum of basis measures in
total variation and proves \eqref{fourier:variation}.
For any compact subset of \(\CC\), the integrand in
\eqref{fourier:measure}, and each of its derivatives with respect to \(s\),
is bounded uniformly in \(t\in[-1,1]\). The integral is consequently
entire, and the finite cardinal sums converge to it locally uniformly.
Away from the nodes, the series in \eqref{fourier:cardinal} converges
locally uniformly and has the same limit. Thus the integral supplies its
entire extension. Conjugate symmetry of the measure makes that extension
real on the real axis.
\end{proof}
In particular, for every integer \(k\ge0\),
\begin{equation}\label{fourier:derivative-bound}
\abs{p^{(k)}(s)}\le\pi^k\norm{\mu}_{\mathrm{TV}}
\qquad(s\in\RR).
\end{equation}
The finite cardinal sums and their derivatives converge locally uniformly:
the assertion for derivatives follows either directly from the measure
bound or from Cauchy's integral formula on a slightly larger compact set.
Passing their node values and derivatives to the limit in
\eqref{fourier:node-expansion} gives the cardinal identities
\begin{equation}\label{fourier:node-values}
p(n)=Q_nc_n,\qquad p'(n)=Q_n(d_n+D_nc_n)
\qquad(n\in\{0\}\cup\mathcal N).
\end{equation}
The other summands in \eqref{fourier:cardinal} vanish to second order at
\(n\), which explains why these identities involve only its own pair of
coefficients.
\subsection{Fourier transformation}
For a cardinal function \(p\), let
\[
q_p(x)=e^{-\pi hs}p(s),\qquad s=b\abs{x}^2.
\]
Define, for real \(-1\le t\le1\),
\begin{equation}\label{fourier:frequencies}
\lambda(t)=\frac{i}{b(t+ih)},\qquad
z(t)=-\frac{B}{t+ih}-ih,
\end{equation}
and, using the measure from Lemma~\ref{fourier:measure-lemma}, put
\begin{equation}\label{fourier:kernel}
K_p(s)=\int_{-1}^{1}\lambda(t)e^{i\pi z(t)s}\,\dd\mu(t).
\end{equation}
\begin{proposition}[Gaussian Fourier pairing]\label{fourier:transform}
The function \(q_p\) is real, radial, and Schwartz on \(\RR^2\), and
\[
\widehat{q_p}(\xi)=e^{-\pi hs}K_p(s),\qquad s=b\abs{\xi}^2.
\]
The function \(K_p\) is entire and real on the real axis. Moreover,
\begin{equation}\label{fourier:damping}
\operatorname{Im}z(t)
=h\left(\frac{B}{t^2+h^2}-1\right)
\ge\frac{26}{435}>0\qquad(-1\le t\le1).
\end{equation}
Consequently every derivative of \(K_p(s)\) decays exponentially as
\(s\to+\infty\).
\end{proposition}
\begin{proof}
By \eqref{fourier:measure},
\[
q_p(x)=\int_{-1}^{1}e^{-\pi b(h-it)\abs{x}^2}\,\dd\mu(t).
\]
Every spatial derivative of the integrand is a polynomial in \(x\), whose
coefficients are bounded uniformly for \(t\in[-1,1]\), times a Gaussian of
modulus \(e^{-\pi bh\abs{x}^2}\). The finite variation of \(\mu\)
justifies differentiation under the integral and bounds every Schwartz
seminorm. Reality follows from the real values of \(p\) on the real axis,
and radiality from the dependence on \(\abs{x}^2\).
For \(\operatorname{Re}a>0\), the two-dimensional Gaussian transform in
our Fourier convention is
\[
\widehat{e^{-\pi a\abs{x}^2}}(\xi)
=a^{-1}e^{-\pi\abs{\xi}^2/a}.
\]
This follows from the one-dimensional Gaussian integral by taking products
and analytically continuing in \(a\) to the right half-plane.
Apply this with \(a=b(h-it)\). Fubini's theorem applies because the
integrands have a common integrable Gaussian majorant. We have
\[
a^{-1}=\lambda(t),\qquad
-\frac{\pi\abs{\xi}^2}{a}+\pi hs=i\pi z(t)s,
\]
which gives \eqref{fourier:kernel} and the asserted Fourier transform.
The curve \(z([-1,1])\) is compact, so \(K_p\) is entire by differentiation
under its defining integral. Also,
\[
\lambda(-t)=\overline{\lambda(t)},\qquad
z(-t)=-\overline{z(t)}.
\]
Thus the same folded measure and real-part convention as in
Lemma~\ref{fourier:measure-lemma} apply to \(K_p\) and its real-axis
derivatives.
Finally, the imaginary part in \eqref{fourier:damping} is smallest at
\(\abs t=1\), where its value is \(26/435\). Differentiating
\eqref{fourier:kernel} gives, for \(s\ge0\),
\[
\abs{K_p^{(k)}(s)}
\le\norm{\mu}_{\mathrm{TV}}
\max_{\abs t\le1}\bigl(\abs{\lambda(t)}\,
(\pi\abs{z(t)})^k\bigr)
e^{-26\pi s/435}.
\]
This proves the decay assertion.
\end{proof}
No cardinal representation of \(K_p\) is assumed. Only its entireness, its values and
derivatives, and the estimates obtained from \eqref{fourier:kernel} will be
used in the interpolation equations. In particular, the measure in that
formula is supported on the real interval \([-1,1]\); later complex-disk
estimates concern specified finite integrands, not an analytic continuation
in the integration variable of an arbitrary \(\ell^1\) density.
\begin{remark}[A real-radius Hankel representation]\label{fourier:hankel}
For real \(s\ge0\), the same transform has the representation
\begin{equation}\label{fourier:hankel-formula}
K_p(s)=e^{\pi hs}\frac{\pi}{b}\int_0^\infty
e^{-\pi hu}p(u)
J_0\!\left(\frac{2\pi\sqrt{us}}{b}\right)\,\dd u,
\end{equation}
where
\(J_0(v)=(2\pi)^{-1}\int_0^{2\pi}e^{-iv\cos\theta}\,\dd\theta\)
for real \(v\). To prove this, use polar coordinates in the Fourier
integral for \(q_p\), perform the angular integral, and put \(u=b r^2\).
The resulting radial integral is
\(\widehat q_p(\xi)=(\pi/b)\int_0^\infty e^{-\pi hu}p(u)
J_0(2\pi\sqrt{us}/b)\,\dd u\), with \(s=b\abs{\xi}^2\).
Absolute convergence follows from \(\abs{p(u)}\le\norm{\mu}_{\mathrm{TV}}\)
and \(\abs{J_0(v)}\le1\). Multiplication by \(e^{\pi hs}\) gives
\eqref{fourier:hankel-formula}. This is an identity on the nonnegative real
axis; it is not used as a numerical sign certificate below.
\end{remark}
\subsection{The coupled interpolation conditions}
We use two cardinal functions \(p_1,p_2\), with coefficients
\(c_{i,n},d_{i,n},C_i\) for function index \(i\in\{1,2\}\). Their
coefficients at positive nodes will be determined later; fix
\begin{equation}\label{fourier:fixed}
(c_{1,0},d_{1,0},C_1)=(1,0.44,-0.013),\qquad
(c_{2,0},d_{2,0},C_2)=(0,-0.368,0.017).
\end{equation}
All terminating decimals in the construction denote exact rational numbers.
Write \(K_i=K_{p_i}\) and \(H_i=p_i+K_{3-i}\). The associated pair is
\begin{equation}\label{fourier:pair}
F(x)=e^{-\pi hs}H_1(s),\qquad
\widehat F(\xi)=e^{-\pi hs}H_2(s),
\quad s=b\abs{x}^2\ \text{or}\ b\abs{\xi}^2,
\end{equation}
respectively. Indeed, \(F=q_{p_1}+\widehat{q_{p_2}}\), and Fourier
inversion gives \(\widehat F=\widehat{q_{p_1}}+q_{p_2}\), since the
functions are even. Thus both functions in \eqref{fourier:pair} are real,
radial, and Schwartz for every pair of admissible coefficient lists.
The required interpolation conditions are
\begin{equation}\label{fourier:interpolation}
\begin{split}
H_1(n)&=H_2(n)=0,\\
H_2'(n)&=0,\qquad H_1'(n)=-Q_n\mathbf1_{\{n=1\}}
\qquad(n\in\mathcal N).
\end{split}
\end{equation}
Here primes denote derivatives with respect to \(s\).
They give zeros of order at least two for both functions at every positive
node except that \(H_1\) has a prescribed simple zero at \(1\).
\begin{lemma}[Coefficient form of the interpolation]\label{fourier:equivalence}
For two real \(\ell^1\) coefficient lists, conditions
\eqref{fourier:interpolation} are equivalent to
\begin{equation}\label{fourier:coefficient-equations}
\begin{split}
c_{i,n}&=-\frac{K_{3-i}(n)}{Q_n},\\
d_{i,n}&=\frac{D_nK_{3-i}(n)-K_{3-i}'(n)}{Q_n}
-\mathbf1_{\{i=1,\ n=1\}}
\qquad(i=1,2,\ n\in\mathcal N).
\end{split}
\end{equation}
\end{lemma}
\begin{proof}
By \eqref{fourier:node-values},
\[
H_i(n)=Q_nc_{i,n}+K_{3-i}(n),\qquad
H_i'(n)=Q_n(d_{i,n}+D_nc_{i,n})+K_{3-i}'(n).
\]
The first equation in \eqref{fourier:coefficient-equations} is therefore
equivalent to \(H_i(n)=0\). Substituting it into the second displayed
identity gives the derivative condition in
\eqref{fourier:interpolation} exactly when the second coefficient equation
holds. Division by \(Q_n\) is legitimate by
\eqref{fourier:node-bounds}.
\end{proof}
==== END A-sharp-Fourier-certificate-for-planar-circle-packing//fourier.tex ====
==== BEGIN A-sharp-Fourier-certificate-for-planar-circle-packing//half-gaps.tex ====
\begin{figure}[htbp]
\centering
\setlength{\unitlength}{0.67cm}
\begin{picture}(14,2.35)
\put(0,1.2){\vector(1,0){13.4}}
\put(13.6,1.05){$s$}
\put(0,1.2){\circle*{.13}}
\put(1,1.2){\circle*{.13}}
\put(3,1.2){\circle*{.13}}
\put(4,1.2){\circle*{.13}}
\put(7,1.2){\circle*{.13}}
\put(9,1.2){\circle*{.13}}
\put(12,1.2){\circle*{.13}}
\put(13,1.2){\circle*{.13}}
\put(0,.65){\makebox(0,0){$0$}}
\put(1,.65){\makebox(0,0){$1$}}
\put(3,.65){\makebox(0,0){$3$}}
\put(4,.65){\makebox(0,0){$4$}}
\put(7,.65){\makebox(0,0){$7$}}
\put(9,.65){\makebox(0,0){$9$}}
\put(12,.65){\makebox(0,0){$12$}}
\put(13,.65){\makebox(0,0){$13$}}
\color{blue!55!black}
\put(.5,1.05){\line(0,1){.3}}
\put(2,1.05){\line(0,1){.3}}
\put(3.5,1.05){\line(0,1){.3}}
\put(5.5,1.05){\line(0,1){.3}}
\put(8,1.05){\line(0,1){.3}}
\put(10.5,1.05){\line(0,1){.3}}
\put(12.5,1.05){\line(0,1){.3}}
\put(4,1.8){\line(1,0){1.5}}
\put(4,1.65){\line(0,1){.3}}
\put(5.5,1.65){\line(0,1){.3}}
\put(6,1.8){\makebox(0,0)[l]{\small one half-gap}}
\end{picture}
\caption{Interpolation nodes (dots) and adjacent gap midpoints (short
marks) in the squared-radius coordinate. Each node-to-midpoint segment
is a half-gap. The finite sign certificates control a Taylor quotient
on each whole half-gap; a separate argument covers the remaining unbounded region.}
\label{fig:half-gaps}
\end{figure}
==== END A-sharp-Fourier-certificate-for-planar-circle-packing//half-gaps.tex ====
==== BEGIN A-sharp-Fourier-certificate-for-planar-circle-packing//independent-certification.tex ====
\section{Alternative certification methods}\label{ind:section}
The primary certificate does not use the methods in this appendix.
We develop several ways to assess or reproduce its finite quantities.
The first gives a fixed rational evaluation of every $N=256$ moment and
proves a uniform error theorem. A separate exact-rational
program checks the first $6\times6$ block with an $N=64$ rule.
Rational-center Taylor integration gives a different, arbitrarily refinable
enclosure method, and the final proposition gives a local formula for the
direct cardinal derivatives.
All terminating decimals denote exact rational numbers. The full
$N=256$ rational recipe and the Taylor method specify computations and
their error bounds without asserting unreported finite values. The
$N=64$ program has the narrower scope stated with its proof.
\subsection{Rational exponential primitives}\label{cert:rational-primitives}
The fixed rational procedures below use rational arithmetic alone, with a
complex number represented by a pair of rationals. Put
\[
\begin{aligned}
p_*&=3.141592653589793238462643383279502884,\\
b_*&=.866025403784438646763723170752936183,
\end{aligned}
\]
and retain $B=4/3$ and $h=2/5$ exactly. Both $|p_*-\pi|$ and $|b_*-b|$
are less than $10^{-36}$. For example, the identity
$\pi/4=4\arctan(1/5)-\arctan(1/239)$ and alternating series through
power $57$ certify the first assertion; comparing $(b_*\pm10^{-36})^2$
with $3/4$ certifies the second.
Define
\[
[w]_{40}=10^{-40}\lfloor10^{40}w\rfloor
\]
for real $w$, and apply it componentwise to complex numbers. Define an
approximation $e_*(u)$ to $e^{i\pi u}$ by the following finite recursion:
\begin{equation}\label{cert:estar}
v_0=\left[\sum_{q=0}^{32}\frac{(ip_*u/1024)^q}{q!}\right]_{40},
\qquad v_{a+1}=[v_a^2]_{40}\ (0\le a<10),\qquad e_*(u)=v_{10}.
\end{equation}
Real and imaginary parts of $e_*(u)$ give the cosine and sine
replacements when $u$ is real.
\begin{lemma}[Rounded exponential]\label{cert:rounded-exponential}
For every $u\in\mathbb Q(i)$ satisfying
$\Re(ip_*u)\le0$ and $|p_*u|\le1000$, the recursion
\eqref{cert:estar} satisfies
\[
|e_*(u)-e^{ip_*u}|<3.542\cdot10^{-34}<10^{-32}.
\]
\end{lemma}
\begin{proof}
The reduced argument has modulus less than one, so the Taylor remainder
in the first step of~\eqref{cert:estar} is less than $3/33!$.
A componentwise floor contributes less than $2\cdot10^{-40}$ in complex
absolute value. Each exact exponential subsequently squared has modulus
at most one. Thus the recurrence
\[
e_0=3/33!+2\cdot10^{-40},\qquad
e_{a+1}=2e_a+e_a^2+2\cdot10^{-40}
\]
bounds the accumulated error. Exact rational evaluation gives
$e_{10}<3.542\cdot10^{-34}$, as claimed.
\end{proof}
\paragraph{Positive exponential upper bounds.}
For $x\ge0$ and an integer $n\ge0$ with $x<n+2$, define
\begin{equation}\label{ind:positive-exponential}
\mathcal T_n(x)=\frac{x^{n+1}}{(n+1)!}\frac1{1-x/(n+2)},\qquad
\mathcal E_n(x)=\sum_{k=0}^n\frac{x^k}{k!}+\mathcal T_n(x).
\end{equation}
After the first omitted term, the ratio of each term to its predecessor is
at most $x/(n+2)<1$, so $e^x\le\mathcal E_n(x)$. These are rational
upper bounds when $x$ is rational. They will be used both in the fixed
first-block check and in the arbitrary-refinement Taylor method.
\subsection{A fixed rational evaluation}\label{cert:rational}
We now apply the fixed primitives to the full $N=256$ rule, using rational
coefficients and rational $C,m,y$ in the finite domain of
Lemma~\ref{cert:moment-error}. After defining the rounded sums, we prove
their uniform approximation error. No matrix, residual, or Bernstein
values are asserted by that error theorem.
In~\eqref{cert:rule}, replace each cosine $\cos(\pi u)$ by $\Re e_*(u)$,
with $u=\alpha_v$ or $u=2a\alpha_v$ as appropriate, and then apply
$[\ ]_{40}$ to the final weight $o_v$ only. Denote the results by
$o_v^{\#},\tau_{jv}^{\#}$. In the displayed coefficient list for $P_l$,
replace $b$ by $b_*$; call the resulting numbers $P_l^{\#}$.
Evaluate~\eqref{cert:density} by replacing $\pi$ by $p_*$, $P_l$ by
$P_l^{\#}$, and every exponential by $e_*$, at the approximate point
$t=\tau_{jv}^{\#}$. Call the resulting density $f_j^{\#}(t)$.
For either these points or the exact atom points $t=t_j$, use
\begin{equation}\label{cert:lambda-z-star}
\lambda_*=[i/(b_*(t+ih))]_{40},\qquad
z_*=[-B/(t+ih)-ih]_{40},
\end{equation}
and set
\begin{equation}\label{cert:rational-factors}
E_l(t)=e_*(tm)\frac{(ip_*ty)^l}{l!},\qquad
G_l(t)=\lambda_*e_*(z_*m)\frac{(ip_*z_*y)^l}{l!}.
\end{equation}
The rational moment sums are therefore
\begin{equation}\label{cert:rational-moment-sums}
\begin{aligned}
U_{p,l}^{\#}
&=\Re\left\{\frac16\sum_{j=1}^6\sum_{v=0}^{255}
o_v^{\#}f_j^{\#}(\tau_{jv}^{\#})E_l(\tau_{jv}^{\#})
+C\sum_{j=0}^6e_jP_j^{\#}E_l(t_j)\right\},\\
X_{p,l}^{\#}
&=\Re\left\{\frac16\sum_{j=1}^6\sum_{v=0}^{255}
o_v^{\#}f_j^{\#}(\tau_{jv}^{\#})G_l(\tau_{jv}^{\#})
+C\sum_{j=0}^6e_jP_j^{\#}G_l(t_j)\right\}.
\end{aligned}
\end{equation}
No additional floors are used. For $Q_n^{\#},D_n^{\#}$ use the product
and cotangent formulas defining $Q_n,D_n$, with the same sine and cosine
replacements and residue $a=n\bmod12\in A$; cotangent means the ratio of
the replaced cosine to the replaced sine.
For the finite pseudocode in Section~\ref{cert:evaluation}, the rational
prescription uses the same measure order: the $6\cdot256$ quadrature pairs
followed by the seven atom pairs. Put $\#$ on the prescribed quantities
and make exactly the replacements above; no other rounding is performed.
\begin{lemma}[Rational moment and node errors]\label{cert:rational-moment-error}
Let $p=PR$ have a rational coefficient list supported in
$\{0\}\cup I_f$, with $\|(c,d)\|_1\le3$, and let $C\in\mathbb Q$
satisfy $|C|\le1$. Let $m,y\in\mathbb Q$ and
$l\in\{0,\ldots,30\}$ satisfy $0\le m\le100$ and $|y|\le3/2$.
The prescription in Section~\ref{cert:rational} gives rational numbers
$U_{p,l}^{\#}$ and $X_{p,l}^{\#}$ such that
\[
|U_{p,l}^{\#}-U_{p,l}|<10^{-14},\qquad
|X_{p,l}^{\#}-X_{p,l}|<10^{-14}.
\]
Separately, for every node $n\in12\ZZ+A$, using its residue $a\in A$
in~\eqref{fourier:node-data}, the prescribed rational values satisfy
\[
|Q_n^{\#}-Q_n|<10^{-25},\qquad |D_n^{\#}-D_n|<10^{-25}.
\]
\end{lemma}
\begin{proof}
We verify the call hypotheses of Lemma~\ref{cert:rounded-exponential}
before using its error bound, and then propagate the errors through the
finite sums. The cosine calls defining the nodes and weights have real arguments
$u=\alpha_v$ or $u=2a\alpha_v$, with $|u|<254$; hence
$|p_*u|<3.142\cdot254<1000$ and $\Re(ip_*u)=0$.
For real $u$, replacing $p_*$ by $\pi$ adds at most $|u|10^{-36}$.
The telescoping coefficient sum in the weights is less than one, so
\[
\begin{aligned}
|\tau_{jv}^{\#}-\tau_{jv}|
&<\frac{3.542\cdot10^{-34}+10^{-36}}{12}<9\cdot10^{-34},\\
|o_v^{\#}-o_v|
&<\frac{3.542\cdot10^{-34}+254\cdot10^{-36}}{256}
+10^{-40}<\frac{10^{-29}}{256}.
\end{aligned}
\]
The exact nodes have a uniform distance from their piece endpoints:
\[
\begin{aligned}
\min\{\tau_{jv}-t_{j-1},t_j-\tau_{jv}\}
&\ge\frac{1-\cos(\pi/512)}{12}
=\frac{\sin^2(\pi/1024)}6\\
&\ge\frac1{6\cdot512^2}=\frac1{1572864}>9\cdot10^{-34}.
\end{aligned}
\]
Here $\sin x\ge2x/\pi$ on $[0,\pi/2]$ follows from concavity.
It follows that every real rounded node $\tau_{jv}^{\#}$ remains strictly
inside $(t_{j-1},t_j)$. The density is therefore evaluated on its stated
piece, and all rounded nodes and exact atoms lie in $[0,1]$.
The remaining real exponential calls in the density and direct factors
now have $|u|\le100$; the residue arguments for $Q_n,D_n$ have $|u|<1$.
They too satisfy the call hypotheses.
For real $t\in[0,1]$, the exact identities
\[
\begin{gathered}
|\lambda(t)|\le\frac5{\sqrt3}<3,\\
|z(t)|^2=h^2+\frac{B^2-2Bh^2}{t^2+h^2}\le\left(\frac{44}{15}\right)^2,\\
\operatorname{Im}z(t)=h\left(\frac B{t^2+h^2}-1\right)
\ge\frac{26}{435}>.0597
\end{gathered}
\]
use $B^2-2Bh^2=304/225>0$. They apply in particular at every rounded
node before the floor in~\eqref{cert:lambda-z-star}. That floor lowers the
imaginary part by less than $10^{-40}$ and changes the complex value by
less than $2\cdot10^{-40}$. Thus
\[
\operatorname{Im}z_*\ge\frac{26}{435}-10^{-40}>0,
\qquad |z_*|<\frac{44}{15}+2\cdot10^{-40}<3.
\]
For the exact transformed arguments, $100\pi|z(t)|<922$; for the rounded
ones, $|p_*z_*m|<100(3.142)(3)<1000$ and
$\Re(ip_*z_*m)=-p_*m\operatorname{Im}z_*\le0$.
The exponential estimate already proved therefore applies to every
transformed call as well.
For a quadrature node write $t=\tau_{jv}$ and $t^{\#}=\tau_{jv}^{\#}$;
at an atom put $t=t^{\#}=t_j$. Along the intervening real segment,
$|t+ih|\ge h$ and $b,b_*>.86$. The derivative bounds for the reciprocal
factors and the change of $b$ give
\[
\begin{aligned}
|\lambda_* -\lambda(t)|
&<8(9\cdot10^{-34})+4\cdot10^{-36}+2\cdot10^{-40}<10^{-29},\\
|z_*-z(t)|
&<\frac{25}{3}(9\cdot10^{-34})+2\cdot10^{-40}<10^{-29}.
\end{aligned}
\]
Indeed $1/(b_*h^2)<8$, $B/h^2=25/3$, and the change of $1/b$ contributes
less than $10^{-36}/(.86^2h)<4\cdot10^{-36}$. At atoms the node-error
terms vanish. Both transformed exponential arguments lie in the closed
left half-plane after multiplication by $i$. The integral first-difference
formula there gives, for example,
\[
|e_*(z_*m)-e^{i\pi z(t)m}|
<10^{-32}+100(3.142\cdot10^{-29}+3\cdot10^{-36})<4\cdot10^{-27}.
\]
The direct exponential-factor error is smaller. The arguments of the
divided powers differ by less than
$1.5(3.142\cdot10^{-29}+3\cdot10^{-36})<5\cdot10^{-29}$.
The divided powers and their derivatives are bounded along the connecting
segment by $e^{1.5\cdot3.142\cdot3}<1.5\cdot10^6$; this also gives the
stated bound $e^{1.5\pi\cdot3}<1.5\cdot10^6$. One rational check of this
exponential bound uses $1.5(3.142)(3)/16<.884$ and
\[
\begin{gathered}
\sum_{q=0}^{8}\frac{.884^q}{q!}
+\frac{.884^9}{9!}\frac1{1-.884/10}<2.421,\\
2.421^{16}<1.5\cdot10^6.
\end{gathered}
\]
The first expression majorizes the exponential series because every
successive ratio in its tail is at most $.884/10$.
The divided-power errors are therefore less than
$7.5\cdot10^{-23}<1.5\cdot10^{-22}$.
Using $|\lambda|<3$, modulus at most one for the exact exponentials, and
modulus less than two for their approximations, the transformed-factor
error is at most
\[
2(1.5\cdot10^6)10^{-29}
+3(1.5\cdot10^6)(4\cdot10^{-27})
+3(1.5\cdot10^{-22})<3\cdot10^{-20}.
\]
The same bound applies to the direct factor.
For completeness, the density estimate can be obtained without relying on
cancellation. Its real exponential factors, including evaluation at a
rounded node, have error less than
$3.542\cdot10^{-34}+100(3.142(9\cdot10^{-34})+10^{-36})<3\cdot10^{-31}$.
Let $A_n,B_n$ denote the inner sums in~\eqref{cert:density} with and without
the factor $t_l-t$, respectively, and put $\#$ on their replacements.
On the common piece, $|t_l-t|,|t_l-t^{\#}|\le1$ for $l\ge j$.
The displayed list for $P_l$ gives
$\sum_{l=1}^6|P_l^{\#}-P_l|<7\cdot10^{-36}$ and
$\sum_{l=1}^6|P_l^{\#}|<11$. Hence $|A_n^{\#}|,|B_n^{\#}|<22$ and
\[
\begin{aligned}
|A_n^{\#}-A_n|,\ |B_n^{\#}-B_n|
&<22(9\cdot10^{-34})+2(7\cdot10^{-36})\\
&\quad+10(3\cdot10^{-31})<4\cdot10^{-30}.
\end{aligned}
\]
The bracket in~\eqref{cert:density} has approximate magnitude at most
$70(|c_n|+|d_n|)$ and error less than
$1.3\cdot10^{-29}(|c_n|+|d_n|)$. Expanding the outer product once and
using list norm at most three gives the pointwise bound
\[
\begin{aligned}
|f_j^{\#}(t^{\#})-f_j(t)|
&<2(10^{-36})(3)(2)(70)\\
&\quad+2(3.142)(3)\bigl(70(3\cdot10^{-31})+1.3\cdot10^{-29}\bigr)\\
&<7\cdot10^{-28}<10^{-22}.
\end{aligned}
\]
On the real segment $|f_j|<1000$, and each atom has magnitude at most six
and error less than $10^{-28}$. The latter error also follows directly
from the displayed $P_l$ list. The exact transformed factor has magnitude
below $3\cdot1.5\cdot10^6$, so every approximate factor has magnitude
below $4.6\cdot10^6$. The exact quadrature weights have total mass one,
and the sum of the absolute weight errors over all six pieces is less than
$10^{-29}$. Separating weight, density, factor, and atom errors, and then
adding the analytic quadrature error, bounds either moment error by
\[
\begin{split}
&10^{-29}(1001)(4.6\cdot10^6)+10^{-22}(4.6\cdot10^6)
+1000(3\cdot10^{-20})\\
&\quad+7\bigl(10^{-28}(4.6\cdot10^6)+6(3\cdot10^{-20})\bigr)
+10^{-22}<5\cdot10^{-16}<10^{-14}.
\end{split}
\]
Finally, for the residue arguments the sine and cosine replacement errors
are less than $3.542\cdot10^{-34}+10^{-36}<4\cdot10^{-34}$.
Every nonzero exact sine has magnitude at least $\sin(\pi/12)>1/4$, and
its replacement has magnitude greater than $1/5$. Each squared factor
$(2\sin)^2$ in $Q_n$ and its replacement is less than five in magnitude,
and their difference is less than $5\cdot10^{-33}$. Telescoping the five
factors and allowing for $(p_*/6)^2-(\pi/6)^2$ gives a $Q_n$ error less
than $2\cdot10^{-29}<10^{-25}$. For a cotangent, division by the two sine
lower bounds gives an error less than
$5(4\cdot10^{-34})+20(4\cdot10^{-34})=10^{-32}$.
Summing five terms and allowing for the prefactor gives a $D_n$ error less
than $6\cdot10^{-32}<10^{-25}$, as required.
\end{proof}
\subsection{Propagation to finite certificates}\label{cert:propagation}
We record stability estimates before applying any numerical margin.
The first controls the two operator coordinates under errors in the moments
and node constants. The second controls Bernstein sums from errors in
the signed-$y$ moments.
To propagate a moment radius, abbreviate $X_\ell=X_{p,\ell}(m,1)$ and
write $\mathcal S_{m,c},\mathcal S_{m,d}$ for the two components of
\eqref{cert:pair}, with hats for their rational approximations.
Suppose $|X_0|,|X_1|\le H$,
$|D_m|\le d$, and rational approximations have errors at most
$\epsilon,\epsilon,\eta_Q,\eta_D$ in $X_0,X_1,Q_m,D_m$, respectively.
If both exact and approximate $Q_m$ exceed $q>0$, expanding the numerator
and then the reciprocal in \eqref{cert:pair} gives
\begin{equation}\label{ind:pair-stability}
\max\{|\widehat{\mathcal S}_{m,c}-\mathcal S_{m,c}|,
|\widehat{\mathcal S}_{m,d}-\mathcal S_{m,d}|\}
\le \frac{(d+1+\eta_D)\epsilon+H\eta_D}{q}
+\frac{(d+1)H\eta_Q}{q^2}.
\end{equation}
For example, the perturbed second numerator differs by at most
$(d+\eta_D)\epsilon+H\eta_D+\epsilon$, and the exact numerator is bounded
by $(d+1)H$; this proves the inequality, which also dominates the first
entry's error. The Taylor budgets and the separated and common-kernel
first-block estimates below use $q=7/4$ and $d=2.12$. The separate-rounding
comparison uses $q=7/4$ and $d=4.3$.
The rational lower bound $Q_m>197192/112500$ proved after
\eqref{fourier:node-bounds} exceeds $7/4$ by $1268/450000$, so
$\eta_Q\le10^{-25}$ also ensures that the approximate $Q_m$ exceeds $q$.
\paragraph{Analytic pair-error allowance.}
For an input $p$ in the rational finite-input domain of
Lemma~\ref{cert:rational-moment-error} and a row $m\in I_f$, define the rational
version $\mathcal S_m^\#$ by using the quantities with superscript $\#$
throughout~\eqref{cert:pair}. We record the division estimate explicitly
for both basis and full-table inputs. By
\eqref{interpolation:piece-envelopes} and the atom masses there, the total
variation of the folded measure for any input in the real finite-input
domain of Lemma~\ref{cert:moment-error} is at most
\[
3\sum_{j=1}^6\frac{M_j}{6}+\sum_{j=0}^6m_j^*=119.642<133.
\]
On the real segment, $|\lambda|,|z|<3$ and
$|e^{i\pi z m}|\le1$ for $m\ge0$. Thus, using $\pi<22/7$,
\[
|X_{p,0}(m,1)|<133(3)<400,\qquad
|X_{p,1}(m,1)|<133(3)(3)(22/7)=3762<4000.
\]
For this calculation abbreviate $X_j=X_{p,j}(m,1)$, $Q=Q_m$, and $D=D_m$.
Write $\varepsilon=10^{-14}$ and $\delta=10^{-25}$. By
Lemma~\ref{cert:rational-moment-error}, the errors in $X_0,X_1$ are less than
$\varepsilon$ and those in $Q,D$ are less than $\delta$.
Equation~\eqref{fourier:node-bounds} gives $Q>1.75$, $|D|<2.12$, and hence
$Q^{\#}>1.7$. For $\mathcal A=DX_0-X_1$ and
$\mathcal A^{\#}=D^{\#}X_0^{\#}-X_1^{\#}$,
\[
|\mathcal A|<4848,\qquad
|\mathcal A^{\#}-\mathcal A|<2.12\varepsilon+(400+\varepsilon)\delta+\varepsilon
<3.121\cdot10^{-14}.
\]
The two quotient differences are consequently bounded by
\[
\begin{aligned}
\left|\frac{X_0^{\#}}{Q^{\#}}-\frac{X_0}{Q}\right|
&<\frac{10^{-14}}{1.7}+\frac{400\cdot10^{-25}}{1.75(1.7)}
<6\cdot10^{-15},\\
\left|\frac{\mathcal A^{\#}}{Q^{\#}}-\frac{\mathcal A}{Q}\right|
&<\frac{3.121\cdot10^{-14}}{1.7}
+\frac{4848\cdot10^{-25}}{1.75(1.7)}
<1.84\cdot10^{-14}<2\cdot10^{-14}.
\end{aligned}
\]
In particular, the uniform error in each entry of $\mathcal S_m^{\#}$ is
strictly less than $10^{-12}$; we retain this looser allowance below.
\subsubsection{Bernstein multipliers}\label{ind:bernstein}
For every signed half-gap in Lemma~\ref{signs:bernstein}, $y\ne0$ and
$|y|^{-\nu}\le4$. Directly from \eqref{cert:moments},
\[
y^{-\nu}\bigl(U_{p,r+\nu}(m,y)+X_{q,r+\nu}(m,y)\bigr)
=y^r\bigl(U_{p,r+\nu}(m,1)+X_{q,r+\nu}(m,1)\bigr).
\]
This identity also holds for negative $y$ and gives exactly the
power coefficients of \eqref{cert:power-from-moments}; it does not replace
their derivatives by interpolated jets. Two moment errors less than
$10^{-14}$ give a $T_r$ error less than $8\cdot10^{-14}$.
For fixed $0\le k\le28$, the rational recurrence $a_0=1$ and
$a_r=a_{r-1}(k+1-r)/(29-r)$ satisfies
\begin{equation}\label{ind:bernstein-multiplier}
a_r=\prod_{j=1}^r\frac{k+1-j}{29-j}
=\frac{k!}{(k-r)!}\frac{(28-r)!}{28!}
=\frac{\binom{k}{r}}{\binom{28}{r}}\qquad(0\le r\le k).
\end{equation}
Each factor is in $[0,1]$. There are
$\sum_{k=0}^{28}k=28\cdot29/2=406$ nonconstant multipliers, in addition
to $a_0=1$ for each $k$. Every Bernstein sum has at most $29$ terms,
so its error is strictly less than
$29(8\cdot10^{-14})=2.32\cdot10^{-12}$.
\input{sections/rational-check}
\subsection{Comparison of first-block majorants}\label{ind:small}
The preceding $N=64$ proof separates the value and derivative integrands.
Two coarser arguments instead use one majorant for both. They concern the
same basis inputs $n\in J$, row nodes $m\in J$, $y=1$, and orders zero and
one, with no atoms. The table compares their disk and node bounds with
the resulting entry errors:
\begingroup
\small
\[
\begin{array}{c|c|c|c}
\text{route}&\text{disk majorant}&|D_m|&\text{entry error}\\ \hline
\text{separated}&M_0<24028,\ M_1<468015&<2.12&<7\cdot10^{-14}\\
\text{common kernel}&M=10^6&<2.12&<4.1\cdot10^{-13}\\
\text{separate rounding}&M=750000&<4.3&<5.23\cdot10^{-13}
\end{array}
\]
\endgroup
The separated route uses the explicit density bound $D_0$ from
Section~\ref{cert:rational-small}. The two common routes use density and
exponential bounds $(1800,2.8)$ and $(2500,3)$, respectively.
\subsubsection{A common kernel majorant}\label{ind:small-kernel}
The rational comparisons
\begin{equation}\label{ind:small-exponentials}
\mathcal E_{32}(44/21)<42/5,\qquad
\mathcal E_{32}(891/875)<14/5
\end{equation}
follow by clearing positive denominators in~\eqref{ind:positive-exponential};
both arguments are nonnegative and less than $34$. Since
$4\pi/6<44/21$ and $4\pi(.081)<891/875$, the two basis densities on the
complex disks satisfy
\[
\begin{aligned}
|u_n(t)|&<2(22/7)^2(65/6)(42/5)=12584/7<1800,\\
|v_n(t)|&<2(22/7)10(42/5)=3696/7<1800.
\end{aligned}
\]
The bounds $|\lambda|<5$, $|e^{i\pi zm}|<2.8$, and
$\pi|z|<6.2\pi<20$ therefore bound both transformed integrands by
$1800(5)(2.8)(20)=504000<10^6$. The positive degree-63 rule and the
half-radius Cauchy estimate give, for each moment,
\begin{equation}\label{ind:small-kernel-quadrature}
|X_\ell-X_{\ell,\mathrm{quad}}|\le10^6 2^{-62}<2.17\cdot10^{-13}.
\end{equation}
The factor two here compares the integral and quadrature averages; a direct
Taylor average has only one tail, as in Lemma~\ref{ind:taylor-average}.
The rounding proof in Section~\ref{cert:rational-small} applies to these
same inputs and gives
finite moment-rounding error below $10^{-26}$, hence below $10^{-14}$.
It also gives $|X_0|,|X_1|<10^4$, $|\delta Q|<3\cdot10^{-35}$,
$|\delta D|<4\cdot10^{-36}$, and both denominators above $7/4$.
Thus the total moment allowance is $2.27\cdot10^{-13}$. In
\eqref{ind:pair-stability}, the additional $Q,D$ contribution is less than
\[
\frac{10^4(4\cdot10^{-36})+(4\cdot10^{-36})(2.27\cdot10^{-13})}{7/4}
+\frac{3.12\cdot10^4(3\cdot10^{-35})}{(7/4)^2}<4\cdot10^{-31}.
\]
Consequently
\begin{equation}\label{ind:small-kernel-entry}
\frac{3.12(2.27\cdot10^{-13})}{1.75}+4\cdot10^{-31}
<4.1\cdot10^{-13}<10^{-12}.
\end{equation}
\subsubsection{A separate rounding majorant}\label{ind:small-rounding}
The weaker bounds $e^{4\pi/6}<9$ and $e^{4\pi(.081)}<3$ also follow
from~\eqref{ind:small-exponentials}. They bound either basis density by
$2500$; for the larger density,
$2(22/7)^2(65/6)9<2500$. With the kernel bound $5$ and derivative
factor $6.2\pi<20$, this gives the distinct common majorant
$M=2500(5)(3)(20)=750000$. Its moment quadrature radius is
\begin{equation}\label{ind:small-rounding-quadrature}
\varepsilon_q=750000\,2^{-62}
<1.626303258728257\cdot10^{-13}.
\end{equation}
Using $|D_m|<4.3$ and $Q_m>7/4$, the pair quadrature error is less than
$(5.3/(7/4))\varepsilon_q<4.926\cdot10^{-13}$.
The same finite rounding allowance below $10^{-14}$ applies because the
schedule and inputs have not changed. Substitution of the same polynomial
$Q,D$ errors in~\eqref{ind:pair-stability}, now with $d=4.3$, bounds their
additional contribution by $6\cdot10^{-31}$. Therefore
\begin{equation}\label{ind:small-rounding-entry}
\frac{5.3}{7/4}\bigl(750000\,2^{-62}+10^{-14}\bigr)
+6\cdot10^{-31}<5.23\cdot10^{-13}<10^{-12}.
\end{equation}
All three routes support the program's $10^{-12}$ entry radius. They do
not extend its coverage: the public rational algorithm still checks only
the 36 first-block entries and the two inverse certificates, not the
residual or Bernstein tables.
\subsection{Rational-center Taylor integration}\label{ind:taylor}
For a finite input in the real moment domain of Lemma~\ref{cert:moment-error}, let
$F^U_{j,\ell}$ and $F^X_{j,\ell}$ be the density \eqref{cert:density}
multiplied by the respective factors in \eqref{cert:analytic-factors}.
The input is supported on nodes at most $100$, has coefficient-list norm at
most $3$ and $|C|\le1$, and the parameters satisfy
$0\le m\le100$, $|y|\le3/2$, and $\ell\in\{0,\ldots,30\}$. The coefficients,
$C,m,y$ are rational when a rational enclosure is requested. The six
centers and the two radii are
\begin{equation}\label{ind:centers}
c_j=\frac{2j-1}{12},\qquad R=\frac16,\qquad a=\frac1{12}
\quad(1\le j\le6).
\end{equation}
Each $F^U_{j,\ell}$ and $F^X_{j,\ell}$ is analytic on a neighborhood of
$|t-c_j|\le R$: the input list is finite and the only rational-kernel pole,
$-ih$, is outside that disk. This analytic assertion is made for finite
inputs. An arbitrary $\ell^1$ list need not give a complex-frequency
density analytic on these disks.
\begin{lemma}[Taylor averages with rational coefficient radii]
\label{ind:taylor-average}
Suppose $F(c_j+x)=\sum_{k\ge0}A_kx^k$ is analytic on a neighborhood of
$|x|\le R$ and $|F|\le M$ there. Let $D\ge0$ be an integer. For
$0\le k\le D$, suppose
$\widehat A_k\in\mathbb Q(i)$ and $e_k\in\mathbb Q_{\ge0}$ satisfy
$|A_k-\widehat A_k|\le e_k$. Then
\begin{align}\label{ind:taylor-average-error}
\left|6\int_{t_{j-1}}^{t_j}F(t)\,\dd t-
\sum_{2r\le D}\widehat A_{2r}\frac{a^{2r}}{2r+1}\right|
\le M2^{-D}+\sum_{2r\le D}e_{2r}\frac{a^{2r}}{2r+1}.
\end{align}
If the final sum of radii is at most $\delta$ for each of the six pieces,
their combined integral has error at most $M2^{-D}+\delta$. Discrete
atoms are not included in this estimate and are to be enclosed separately.
\end{lemma}
\begin{proof}
Since $2a=1/6$, the expression $6\int F$ is the interval average. The
average of $x^k$ on $[-a,a]$ is zero for odd $k$ and is $a^k/(k+1)$ for
even $k$. Thus the degree-$D$ Taylor polynomial has the exact average
\begin{equation}\label{ind:even-average}
\sum_{2r\le D} A_{2r}\frac{a^{2r}}{2r+1}.
\end{equation}
Cauchy's estimate gives $|A_k|\le M R^{-k}$. Bounding the average of the
absolute tail by its supremum on $[-a,a]$ gives the conservative bound
\[
M\sum_{k>D}(a/R)^k=M\sum_{k>D}2^{-k}=M2^{-D}.
\]
Replacing the even coefficients adds the displayed weighted sum of their
radii. Finally, each of the six integrals is $1/6$ times its average, so
the six factors $1/6$ sum to one. Taking a real part cannot increase the
error. Atoms have no integral tail of this kind.
\end{proof}
Here are explicit rational operations for obtaining the coefficient
enclosures in this lemma. Write
$\mathcal B(a,r)=\{w\in\mathbb C:|w-a|\le r\}$, where
$a\in\mathbb Q(i)$ and $r\in\mathbb Q_{\ge0}$, and put
$\langle a\rangle_1=|\Re a|+|\Im a|$. The following enclosure operations
use rational centers and radii:
\begin{align}\label{ind:ball-operations}
\mathcal B(a,r)+\mathcal B(b,s)&\subseteq\mathcal B(a+b,r+s),\nonumber\\
\mathcal B(a,r)\mathcal B(b,s)&\subseteq
\mathcal B\bigl(ab,\langle a\rangle_1s+\langle b\rangle_1r+rs\bigr),\nonumber\\
\mathcal B(a,r)^{-1}&\subseteq\mathcal B(a^{-1},r/\ell_0^2)
\quad\text{if }0<\ell_0\le |a|-r.
\end{align}
In the last line $a^{-1}\in\mathbb Q(i)$, and the condition on a rational
$\ell_0$ can be checked by squaring rational quantities. The bound follows
from $|w^{-1}-a^{-1}|=|w-a|/(|w||a|)$. These operations allow repeated
use of the same enclosure; possible dependence between its occurrences
does not invalidate inclusion.
Constant exponentials also admit arbitrarily small rational radii.
Use the rational bounds $\mathcal T_n,\mathcal E_n$ from
\eqref{ind:positive-exponential}. If
$w\in\mathcal B(w_0,r)$, choose rational $q\ge|w_0|$, for example
$q=\langle w_0\rangle_1$, and nonnegative integers $n,s$
with $q<n+2$ and $q+r<s+2$. Then
\begin{equation}\label{ind:constant-exponential}
e^w\in\mathcal B\left(\sum_{k=0}^n\frac{w_0^k}{k!},
\mathcal T_n(q)+r\mathcal E_s(q+r)\right).
\end{equation}
Indeed, the first radius bounds the complex Taylor tail at $w_0$.
The integral first-difference formula for $e^w$ bounds the change from
$w_0$ to $w$ by $r e^{q+r}$, giving the second radius. Both preconditions
in \eqref{ind:positive-exponential} are explicit; no estimate of a positive
exponential by a truncated sum alone is being used. For bounded $q$,
$\mathcal T_n(q)\to0$ as $n\to\infty$.
For clarity, the coefficient recurrences needed for the present integrands
are as follows. If $a(x)=\sum a_kx^k$ and $b(x)=\sum b_kx^k$, multiplication
is convolution,
\begin{equation}\label{ind:convolution}
(ab)_k=\sum_{v=0}^k a_v b_{k-v}.
\end{equation}
If $a_v\in\mathcal B(\widehat a_v,e_v)$ and
$b_v\in\mathcal B(\widehat b_v,f_v)$, a rational radius for the center
$\sum_{v=0}^k\widehat a_v\widehat b_{k-v}$ is explicitly
\[
\sum_{v=0}^k\left(\langle\widehat a_v\rangle_1f_{k-v}
+\langle\widehat b_{k-v}\rangle_1e_v+e_vf_{k-v}\right).
\]
For a reciprocal series $q(x)=a(x)^{-1}$ with $a_0\ne0$, and for an
exponential series $E(x)=e^{g(x)}$, the recurrences are
\begin{align}\label{ind:reciprocal-exponential-recurrences}
q_0&=a_0^{-1},&
q_k&=-q_0\sum_{v=1}^k a_vq_{k-v}\quad(k\ge1),\nonumber\\
E_0&=e^{g_0},&
E_k&=\frac1k\sum_{v=1}^k v g_v E_{k-v}\quad(k\ge1).
\end{align}
The second line follows by comparing coefficients in $E'=g'E$.
Use \eqref{ind:constant-exponential} for $E_0$ and
\eqref{ind:ball-operations} in these recurrences. Along with the displayed
convolution radius, this specifies a rational radius for every coefficient,
not merely a formal power series.
More particularly, $d_j=c_j+ih\in\mathbb Q(i)$ and
\begin{equation}\label{ind:kernel-recurrences}
r_0=d_j^{-1},\quad r_{k+1}=-r_k/d_j,\qquad
\lambda_k=ir_k/b,\quad z_0=-Br_0-ih,\quad z_k=-Br_k\ (k\ge1)
\end{equation}
are the coefficients of $(t+ih)^{-1}$, $\lambda(t)$, and $z(t)$ at
$t=c_j+x$. The $r_k,z_k$ are rational complex numbers; an enclosure for
$b^{-1}$ gives the radii for $\lambda_k$. On piece $j$, put
\[
L_{n,j}=\sum_{v=j}^6P_v e^{i\pi t_vn},\qquad
H_{n,j}=\sum_{v=j}^6t_vP_v e^{i\pi t_vn}.
\]
For $\Phi_n(x)=e^{-i\pi n(c_j+x)}$, use
$\phi_{n,0}=e^{-i\pi nc_j}$ and
$\phi_{n,k+1}=(-i\pi n)\phi_{n,k}/(k+1)$. The density is then the finite
series expression
\begin{equation}\label{ind:density-recurrence}
f_j(c_j+x)=2\pi\sum_n\Phi_n(x)
\left\{-\pi c_n\bigl(H_{n,j}-(c_j+x)L_{n,j}\bigr)+id_nL_{n,j}\right\}.
\end{equation}
The constant exponentials in $L_{n,j},H_{n,j},\phi_{n,0}$ are enclosed by
\eqref{ind:constant-exponential}. For a series $v(x)$ equal to $c_j+x$
or $z(c_j+x)$, the divided powers obey
\begin{equation}\label{ind:divided-power-recurrence}
V_0(x)=1,\qquad V_{\ell+1}(x)=
\frac{V_\ell(x)\,i\pi y v(x)}{\ell+1}.
\end{equation}
Finally apply the exponential recurrence to $g=i\pi m(c_j+x)$ or
$g=i\pi m z(c_j+x)$, and multiply the resulting series by the density,
the divided power, and, in the transformed case, $\lambda$. Truncating
these recurrences after degree $D$ gives every coefficient through $D$
exactly at the level of formal series; omitted higher coefficients cannot
enter a lower convolution coefficient.
The seven possible atoms are handled separately. At $t_j$, enclose their
masses $e_jCP_j$ and the appropriate factor in
\eqref{cert:analytic-factors} by the same rational operations and constant
exponential bound. If their respective centers and radii are
$(a_j,r_j)$ and $(g_j,s_j)$, their product has radius
\begin{equation}\label{ind:atom-radii}
\eta_j=\langle a_j\rangle_1s_j+\langle g_j\rangle_1r_j+r_js_j.
\end{equation}
The sum $\eta=\sum_{j=0}^6\eta_j$ is added to the six-piece error; there is
no Taylor integration tail for an atom. When $C=0$, all these terms vanish.
All of these radii can be refined effectively. Alternating rational
arctangent sums in Machin's identity enclose $\pi$ to arbitrarily small
width, and rational bisection of the positive square root of $3/4$ does the
same for $b$. The kernel denominators have $|c_j+ih|\ge h=2/5$, and
$b>4/5$. Taking real and imaginary parts of the constant exponential
enclosures also encloses cosine and sine with the same radius. The nonzero
sine denominators in \eqref{fourier:node-data} have modulus greater than
$1/4$. Consequently sufficiently refined enclosures certify fixed
positive lower bounds for the moduli of all denominators used above or in
computing $Q_m,D_m$, and separately a positive lower bound for $Q_m$ at
pair division. For each fixed $D$ there are finitely many operations in
\eqref{ind:convolution}--\eqref{ind:divided-power-recurrence}. Their
radius formulas tend to zero when their input radii tend to zero, while
\eqref{ind:constant-exponential} has radius tending to zero by refining
its input and increasing its two truncation orders. Thus one may refine
the rational inputs and exponential orders, recompute the displayed
rational radii, and stop when any prescribed finite collection of positive
coefficient and atom budgets holds. The convergence just proved ensures
that this rational stopping test terminates. This is a computability
statement for each specified finite set of inputs, not a claim that any
of its moment sums have been evaluated.
We next make sufficient budgets explicit. Put
$h_D=\lfloor D/2\rfloor+1$. For rational $\delta>0$, the conditions
\begin{equation}\label{ind:coefficient-budget}
e_{j,2r}\le \frac{\delta(2r+1)12^{2r}}{h_D}
\quad(0\le2r\le D)
\end{equation}
on each piece and moment imply that the weighted coefficient-radius sum in
\eqref{ind:taylor-average-error} is at most $\delta$.
We now combine these integration budgets with the pair estimate
\eqref{ind:pair-stability} and the Bernstein estimate from
Section~\ref{ind:bernstein}.
\begin{proposition}[Two conditional Taylor budgets]\label{ind:taylor-budgets}
The following choices give sufficient rational error enclosures for the
same finite quantities used in Appendix~\ref{sec:certificates}.
For a basis input at $n\in J=\{1,3,4\}$, a row $m\in J$, $y=1$, and
$\ell=0,1$, take $D=63$, $h_D=32$, $M=10^6$, and $\delta=10^{-14}$ in
\eqref{ind:coefficient-budget}. There are no atoms. If $Q_m,D_m$ have
rational errors at most $10^{-25}$, either entry of \eqref{cert:pair}
has error less than $2.13\cdot10^{-13}<10^{-12}$.
For the full rational finite-input domain at the start of this subsection,
take $D=255$, $h_D=128$, $M=10^{53}$, and $\delta=4\cdot10^{-15}$ in
\eqref{ind:coefficient-budget}, and require each of the seven atom-product
radii in \eqref{ind:atom-radii} to be at most $4\cdot10^{-15}/7$.
Each $U$ and $X$ moment then has error less than $10^{-14}$. For node
rows $m\in I_f$, with rational $Q_m,D_m$ errors at most $10^{-25}$, either
entry of \eqref{cert:pair} has error less than $1.8\cdot10^{-14}$.
\end{proposition}
\begin{proof}
The first-block majorant $10^6$ is proved in
Section~\ref{ind:small-kernel}. Lemma~\ref{ind:taylor-average} gives
\[
\epsilon\le10^6 2^{-63}+10^{-14}<1.19\cdot10^{-13}.
\]
On the real interval a basis density has modulus below $214$,
$|\lambda|<3$, and $\pi|z|<10$; thus the two exact moments have modulus
below $10^4$. Substitution of $H=10^4$, $q=7/4$, $d=2.12$, and
$\eta_Q=\eta_D=10^{-25}$ in \eqref{ind:pair-stability} gives
\[
\frac{(3.12+10^{-25})1.19\cdot10^{-13}+10^4 10^{-25}}{7/4}
+\frac{3.12\cdot10^4 10^{-25}}{(7/4)^2}<2.13\cdot10^{-13}.
\]
The pointwise majorant $10^{53}$ for the full domain was proved in
Lemma~\ref{cert:moment-error}. In the second case the total moment radius
is at most
\[
10^{53}2^{-255}+8\cdot10^{-15}
<1.728\cdot10^{-24}+8\cdot10^{-15}<10^{-14}.
\]
For the value and first-derivative moments in \eqref{cert:pair}, the real
density is bounded by $1000$, the sum of seven atom moduli by $42$, and
the two factors by $3$ and $30$. Hence both exact moments are bounded
by $H=4\cdot10^4$. The right side of \eqref{ind:pair-stability} is then
at most
\[
\frac{(3.12+10^{-25})10^{-14}+4\cdot10^4 10^{-25}}{7/4}
+\frac{3.12\cdot4\cdot10^4 10^{-25}}{(7/4)^2}<1.8\cdot10^{-14}.
\]
All inequalities here are comparisons of the displayed rational bounds.
The refinement argument above proves that their coefficient and atom
budgets can be met for any prescribed finite set of these moments.
\end{proof}
For a matrix with $r$ rows and entry error at most $\varepsilon$, the
induced column-norm error is at most $r\varepsilon$, and multiplication
by $W_\eta$ adds the factor $\|W_\eta\|_1$. In particular $r=6$ for the
$J$ row block and $r=12$ for the rows on $E\cap[7,16]$. In a tabulated
residual, the other summands in \eqref{cert:residual-bounds} are exact
rationals, so its error is just the corresponding pair error. Finally,
for the signed half-gaps of Lemma~\ref{signs:bernstein}, the two moment
errors below $10^{-14}$ give a power-coefficient error below
$8\cdot10^{-14}$ and a Bernstein-coefficient error below
$2.32\cdot10^{-12}$, as proved in Section~\ref{ind:bernstein}.
These bounds specify the margins required in a Taylor-based comparison;
they do not themselves assert the matrix, residual, or Bernstein values
printed in Appendix~\ref{sec:certificates}.
\subsection{Direct node-Taylor representation}\label{ind:node-taylor}
There is also a local Taylor representation for the direct cardinal
function that does not integrate its spectral density. It applies to
the full $\ell^1$ class of \eqref{fourier:cardinal}. Write
$\mathcal N_0=\{0\}\cup\mathcal N$ and, for $m\in\mathcal N_0$, put
\begin{equation}\label{ind:node-radius}
\rho_m=\min_{\substack{n\in\mathcal N_0\\n\ne m}}|m-n|.
\end{equation}
The consecutive gaps in one period of $\mathcal N_0$ are
$1,2,1,3,2,3$. Thus the minimum exists and is positive:
$\rho_m=1$ for $m=0$ and for residues $0,1,3,4$, while $\rho_m=2$
for residues $7,9$. This is a guaranteed radius to the nearest other
node; it may be smaller than the analytic radius for a particular list
whose corresponding coefficients vanish.
\begin{proposition}\label{ind:node-taylor-proposition}
Let $p=PR$ have a real $\ell^1$ coefficient list as in
\eqref{fourier:list-norm}. Define
\begin{equation}\label{ind:node-P-coefficients}
A_k(m)=\frac{P^{(k)}(m)}{k!}
=\frac1{k!}\sum_{v=-6}^6P_v(i\pi t_v)^k e^{i\pi t_vm}.
\end{equation}
Then $A_0(m)=A_1(m)=0$. The regular series has coefficients
\begin{equation}\label{ind:node-regular-coefficients}
r_j(m)=C\mathbf1_{\{j=0\}}+
(-1)^j\sum_{\substack{n\in\mathcal N_0\\n\ne m}}
\left(\frac{(j+1)c_n}{(m-n)^{j+2}}
+\frac{d_n}{(m-n)^{j+1}}\right),
\end{equation}
and converges for $|u|<\rho_m$. For $0<|u|<\rho_m$, the Laurent
expansion is
\[
R(m+u)=\frac{c_m}{u^2}+\frac{d_m}{u}
+\sum_{j=0}^{\infty}r_j(m)u^j.
\]
For every integer $\ell\ge0$ the direct derivative is therefore
\begin{equation}\label{ind:node-derivatives}
\frac{p^{(\ell)}(m)}{\ell!}
=c_mA_{\ell+2}(m)+d_mA_{\ell+1}(m)
+\sum_{k=2}^{\ell}A_k(m)r_{\ell-k}(m),
\end{equation}
where the last sum is empty for $\ell=0,1$.
\end{proposition}
\begin{proof}
For $n\ne m$, the geometric series for $(m-n+u)^{-1}$ and its derivative
give
\[
\frac1{m-n+u}=\sum_{j\ge0}\frac{(-1)^ju^j}{(m-n)^{j+1}},\qquad
\frac1{(m-n+u)^2}=\sum_{j\ge0}
\frac{(-1)^j(j+1)u^j}{(m-n)^{j+2}},
\]
whenever $|u|<|m-n|$. To justify summing these series over the infinite
list, fix $0\le r<\rho_m$ and set $q=r/\rho_m<1$. Since
$|m-n|\ge\rho_m$ for every other node, the sum of absolute values on
$|u|\le r$ is at most
\[
\frac{\sum_{n\ne m}|c_n|}{\rho_m^2}
\sum_{j\ge0}(j+1)q^j
+\frac{\sum_{n\ne m}|d_n|}{\rho_m}\sum_{j\ge0}q^j
=
\frac{\sum_{n\ne m}|c_n|}{\rho_m^2(1-q)^2}
+\frac{\sum_{n\ne m}|d_n|}{\rho_m(1-q)}<\infty.
\]
This proves absolute uniform convergence on every smaller closed disk,
permits exchanging the node and Taylor sums, and gives
\eqref{ind:node-regular-coefficients}. It also proves precisely why the
nearest-other-node distance is a valid common radius.
The finite Fourier expansion \eqref{fourier:polynomial} gives
\eqref{ind:node-P-coefficients}, and the double zero of $P$ gives
$A_0=A_1=0$. Multiply its entire Taylor series by the convergent Laurent
series for $R$. The two principal terms are canceled by the double zero;
collecting the coefficient of $u^\ell$ gives
\eqref{ind:node-derivatives}. In particular $A_2=Q_m$ and $A_3=Q_mD_m$,
so the cases $\ell=0,1$ reproduce \eqref{fourier:node-values}.
\end{proof}
For finite rational input, the constants in this direct formula can be
enclosed by the rational operations above. For an arbitrary $\ell^1$
input the formula is an exact convergence statement; a computable error
for truncating its node sum additionally requires a bound for that list's
tail. It represents the direct function $p$ only. The transformed
function $K_p$ still has the kernel representation
\eqref{fourier:kernel}, and this local identity neither places $K_p$ in
the cardinal class nor replaces the exact-derivative certificates for the
tabulated functions.
==== END A-sharp-Fourier-certificate-for-planar-circle-packing//independent-certification.tex ====
==== BEGIN A-sharp-Fourier-certificate-for-planar-circle-packing//interpolation.tex ====
\section{Solving the interpolation equations}\label{sec:interpolation}
The equations \eqref{fourier:coefficient-equations} couple two coefficient
lists through the transformed part of the construction. We solve them on
$\ell^1$, using the strict damping of the Fourier kernel to control every
row beyond a fixed finite set.
Let
\[
X=\ell^1(\mathcal N;\RR^2),\qquad
\norm{x}_1=\sum_{n\in\mathcal N}\bigl(\abs{c_n}+\abs{d_n}\bigr)
\quad\text{for }x=((c_n,d_n))_{n\in\mathcal N}.
\]
When writing matrices, we order all $c$ coordinates first and all $d$
coordinates second, with increasing nodes within each group. All operator
norms below are the induced $\ell^1$ norms, equivalently the maximum
absolute column sum for a finite matrix. A subscript specifying a set of
nodes includes both coordinates at each node; a star denotes all positive
nodes.
For $x\in X$, let $p_x$ be the function \eqref{fourier:cardinal} with these
positive-node coefficients and with $C=c_0=d_0=0$. Define the real linear
operator $S$ by
\begin{equation}\label{interpolation:S}
(Sx)_{c,n}=-\frac{K_{p_x}(n)}{Q_n},\qquad
(Sx)_{d,n}=\frac{D_nK_{p_x}(n)-K_{p_x}'(n)}{Q_n}.
\end{equation}
The bounds proved below show, in particular, that $S$ is a bounded operator
on $X$.
We use two successive decompositions of this operator. Set
\begin{equation}\label{interpolation:node-blocks}
J=\{1,3,4\},\qquad E=\mathcal N\cap[7,100],\qquad
I_f=J\cup E=\mathcal N\cap[1,100],\qquad D=S_{J,J}.
\end{equation}
A finite inverse handles the coordinates on $J$. A first Schur complement
then incorporates $E$, using separate bounds for its rows through $16$
and beyond $16$; a second incorporates the infinite tail beyond $100$.
The tail bounds apply to output rows uniformly over all input columns.
The first output rows require a separate, larger bound even for columns
indexed by large nodes.
\subsection{Uniform envelopes for the Fourier terms}
We use the following rational upper and lower bounds. As throughout the
construction, every finite decimal denotes its exact rational value.
\begin{table}[htbp]
\centering
\begin{tabular}{c|rrr|rr}
$j$ & $b_j^*$ & $a_j^*$ & $g_j$ & $M_j$ & $\beta_j$\\ \hline
0 & 2.887 & 2.934 & 2.933 & --- & ---\\
1 & 2.665 & 2.713 & 2.440 & 56 & 0\\
2 & 2.218 & 2.269 & 1.567 & 48 & 4.8\\
3 & 1.804 & 1.859 & .900 & 35 & 3.17\\
4 & 1.486 & 1.548 & .482 & 25 & 1.946\\
5 & 1.250 & 1.320 & .224 & 19 & 1.217\\
6 & 1.073 & 1.151 & .0597 & 6.284 & .792
\end{tabular}
\caption{Envelopes for the folded Fourier integrals.}
\label{interpolation:envelopes}
\end{table}
At each endpoint $t_j=j/6$,
\[
\abs{\lambda(t_j)}\le b_j^*,\qquad
\abs{z(t_j)}\le a_j^*,\qquad
\operatorname{Im}z(t_j)\ge g_j.
\]
On the interval $t_{j-1}\le t\le t_j$, the corresponding bounds are
\begin{equation}\label{interpolation:piece-envelopes}
\abs{\lambda(t)}\le b_{j-1}^*,\quad
\abs{z(t)}\le a_{j-1}^*,\quad
\operatorname{Im}z(t)\ge g_j+\beta_j(t_j-t),\quad
\abs{u_n(t)},\abs{v_n(t)}\le M_j.
\end{equation}
The density bounds hold uniformly for
$n\in\{0\}\cup\mathcal N$.
Here are direct ways to verify these estimates. The identities
\[
\abs{\lambda(t)}=\frac{1}{b\sqrt{t^2+h^2}},\qquad
\abs{z(t)}^2=h^2+\frac{B^2-2Bh^2}{t^2+h^2},\qquad
\operatorname{Im}z(t)=h\left(\frac{B}{t^2+h^2}-1\right)
\]
show that the first two quantities decrease on $[0,1]$. The nonnegative
quantity $2hBt/(t^2+h^2)^2$, which is minus the derivative of
$\operatorname{Im}z$, has its minimum on each of the six intervals at an
endpoint. Its only interior critical point on the positive axis is a
maximum. These observations give the first three bounds in
\eqref{interpolation:piece-envelopes} by endpoint substitution.
For the densities, $e^{i\pi t_l n}$ depends only on $n$ modulo $12$,
whereas $\abs{e^{-i\pi tn}}=1$ on the real interval. The norm of an
affine function of $t$ is convex. Thus an upper bound for both folded
densities on the $j$th interval is
\begin{equation}\label{interpolation:density-check}
2\max_{a\in A}\max\left\{
\pi\abs{\sum_{l=j}^6P_le^{i\pi t_la}},\quad
\pi^2\max_{u\in\{t_{j-1},t_j\}}
\abs{\sum_{l=j}^6(t_l-u)P_le^{i\pi t_la}}
\right\}.
\end{equation}
Substitution gives a value less than the displayed $M_j$ in every case.
For the constant input $C=1$, the absolute folded atom masses are
\[
(m_0^*,\ldots,m_6^*)=(5,4,4,4,2,4,2).
\]
For $s\ge0$, set
\begin{equation}\label{interpolation:L}
\begin{aligned}
L_j(s)&=b_{j-1}^*e^{-\pi g_js}
\min\left\{\frac16,\frac{1}{\pi\beta_js}\right\}
&&(1\le j\le6),\\
L_j^P(s)&=m_j^*b_j^*e^{-\pi g_js}
&&(0\le j\le6).
\end{aligned}
\end{equation}
When $\beta_js=0$, the minimum in the first line is understood to be
$1/6$. Integrating the exponential envelope in
\eqref{interpolation:piece-envelopes} over a piece of length $1/6$ gives
its factor $L_j(s)$. Differentiating the kernel contributes
$(i\pi z(t))^\ell$. Consequently, for every integer $\ell\ge0$ and every
input $p$ of the form \eqref{fourier:cardinal},
\begin{equation}\label{interpolation:moment-bound}
\abs{K_p^{(\ell)}(s)}
\le w\sum_{j=1}^6 M_jL_j(s)(\pi a_{j-1}^*)^\ell
+\abs C\sum_{j=0}^6L_j^P(s)(\pi a_j^*)^\ell,
\qquad s\ge0,
\end{equation}
where $w=\sum_{n\in\{0\}\cup\mathcal N}(\abs{c_n}+\abs{d_n})$.
The argument applies first to finite lists and then to all such lists by
the finite-measure representation and absolute summability.
Both $Q_n$ and $D_n$ are $12$-periodic on the nodes. Moreover,
\[
L_j(n+12k)\le e^{-12k\pi g_j}L_j(n)\qquad(k\ge0),
\]
because the minimum factor in \eqref{interpolation:L} is nonincreasing.
For every real cutoff $x\ge0$, the value and derivative terms in
\eqref{interpolation:S}, together with
\eqref{interpolation:moment-bound}, therefore imply
\begin{equation}\label{interpolation:row-envelope}
\norm{S_{>x,*}}
\le
\sum_{n\in\mathcal N\cap(x,x+12]}
\sum_{j=1}^6
\frac{1+\abs{D_n}+\pi a_{j-1}^*}{Q_n}\,
\frac{M_jL_j(n)}{1-e^{-12\pi g_j}}.
\end{equation}
Here $>x$ means a strict inequality for the row node. This is a bound
uniform over all input columns. It applies as well to a single input coefficient
$c_0$ or $d_0$.
For completeness, this estimate also justifies the operator definition on
the whole space. Apply the bound first to one positive-node basis column,
then sum with the absolute coefficients of a finitely supported list.
The case $x=0$ bounds the resulting output norm by a constant times the
input norm. Finite lists are dense in $X$, so this map has a unique bounded
extension to $X$. The basis measures converge in total variation by
Lemma~\ref{fourier:measure-lemma}; hence $K_{p_x}(n)$ and $K_{p_x}'(n)$
converge at each fixed node as the lists are truncated. The extended map
therefore has exactly the coordinates in \eqref{interpolation:S}, not merely
the coordinates of an auxiliary matrix limit.
For $x=0,4,16,100$, respectively, the right side of
\eqref{interpolation:row-envelope} is less than
\[
27.7,\qquad .272,\qquad .0145,\qquad 4.7\cdot10^{-10}.
\]
We shall use the slightly weaker bounds
\begin{equation}\label{interpolation:tails}
\norm{S_{>0,*}}<28,\qquad
\norm{S_{>4,*}}<.28,\qquad
\norm{S_{>16,*}}<.015,\qquad
\norm{S_{>100,*}}<5\cdot10^{-10}.
\end{equation}
For the input $C=1$, the analogue of
\eqref{interpolation:row-envelope} replaces $M_jL_j$ by $L_j^P$,
uses $a_j^*$ in place of $a_{j-1}^*$, and sums over $0\le j\le6$.
For rows greater than $100$, it gives the bound
\begin{equation}\label{interpolation:atom-tail}
\sum_{n\in\mathcal N,\ n>100}
\left(\frac{\abs{K_P(n)}}{Q_n}
+\frac{\abs{D_nK_P(n)-K_P'(n)}}{Q_n}\right)
<3.2\cdot10^{-8}.
\end{equation}
The scalar substitutions in these envelope estimates use only the explicit
constants and the node data. Section~\ref{cert:scalar-bounds} records the
finite scalar bounds and the precise coverage of their supplied verifier.
\begin{remark}[Compactness of the interpolation operator]
\label{interpolation:compactness}
Let $P_M$ keep the output coordinates with node at most $M$. Then $P_MS$
has finite rank. The right side of \eqref{interpolation:row-envelope}
tends to zero as $x\to\infty$: there are only six nodes in any period,
the periodic $Q_n,D_n$ are uniformly bounded with $Q_n$ bounded away from
zero, and every $L_j(n)$ contains the factor $e^{-\pi g_jn}$ with $g_j>0$.
Thus $\norm{S-P_MS}\to0$. The operator $S$ is a norm limit of finite-rank
operators and is compact on $X$. This property describes the tail
structure; the explicit inverse below is obtained independently from the
quantitative block estimates.
\end{remark}
\subsection{The finite certificates}
There are $51$ nodes in $I_f$ and hence $102$ scalar coordinates.
The approximate lists $\widetilde x_1,\widetilde x_2$ are given by
Table~\ref{tab:coefficients}, with the indicated factor $10^{-10}$,
and are zero at positive nodes beyond $100$. Append the fixed data
\eqref{fourier:fixed} to obtain $\widetilde p_i$, and write
$\widetilde K_i=K_{\widetilde p_i}$.
\begin{lemma}[Finite certificates]\label{interpolation:finite}
For each $\eta\in\{1,-1\}$, the matrix $I-\eta D$ is invertible, and
\begin{equation}\label{interpolation:finite-matrices}
\norm{(I-\eta D)^{-1}}<32,\qquad
\norm{(I-\eta D)^{-1}S_{J,E}}<3.7,\qquad
\norm{S_{E\cap[7,16],J}}<.09.
\end{equation}
For $i=1,2$ and $n\in I_f$, the table satisfies
\begin{equation}\label{interpolation:finite-residuals}
\begin{split}
\abs{\widetilde c_{i,n}+\frac{\widetilde K_{3-i}(n)}{Q_n}}
&<10^{-9},\\
\abs{\widetilde d_{i,n}+\mathbf1_{\{i=1,\ n=1\}}
+\frac{\widetilde K_{3-i}'(n)-D_n\widetilde K_{3-i}(n)}{Q_n}}
&<10^{-9}.
\end{split}
\end{equation}
For each $i$, the norm of the tabulated list, including node zero, is less
than $2.29$, and its norm on nodes greater than $40$ is less than $.00002$.
\end{lemma}
\begin{proof}
The primary interval calculation in Appendix~\ref{sec:certificates} proves
the exact bounds~\eqref{cert:matrix-table}, \eqref{cert:exterior-block}, and
\eqref{cert:residual-bounds}, using the quadrature assertion of
Lemma~\ref{cert:moment-error} for the analytic moment remainder.
The elementary argument in Section~\ref{cert:matrix-bounds} gives the
matrix conclusions, and~\eqref{cert:table-norms} gives the list norms.
\end{proof}
\subsection{Inversion on the whole coefficient space}
We first prove an inverse bound that keeps finite and tail residuals
separate.
\begin{lemma}\label{interpolation:inverse}
For each $\eta\in\{1,-1\}$, $I-\eta S$ has a bounded inverse on $X$.
If $r=r_f+r_t$, with $r_f$ supported on $I_f$ and $r_t$ on
$\mathcal N\cap(100,\infty)$, then
\begin{equation}\label{interpolation:inverse-bound}
\norm{(I-\eta S)^{-1}r}_1
\le91\norm{r_f}_1+2540\norm{r_t}_1.
\end{equation}
\end{lemma}
\begin{proof}
Fix $\eta$ and first restrict to the finite block $I_f$. Put
\[
V=(I-\eta D)^{-1},\qquad L=S_{E,J},\qquad U=S_{J,E}.
\]
The first block of tail estimates and the finite certificates give
\[
\norm L
\le\norm{S_{E\cap[7,16],J}}+\norm{S_{E\cap(16,100],J}}
<.09+.015<.11.
\]
Eliminating the $J$ coordinates from
$(I-\eta S_{I_f,I_f})x_f=r_f$ leaves on $E$ the operator
\[
I-\eta S_{E,E}-LVU.
\]
Its difference from the identity has norm less than
\[
\norm{S_{E,E}}+\norm L\norm{VU}
<.28+.11\cdot3.7=.687<1.
\]
It is therefore invertible by a Neumann series. Solving this block and
substituting into the $J$ equations yields an inverse $V_f$ for the whole
finite compression, with
\begin{equation}\label{interpolation:finite-inverse-bound}
\begin{split}
\norm{V_f}
&\le\max\left\{
32+\frac{(1+3.7)\cdot.11\cdot32}{1-.687},\quad
\frac{1+3.7}{1-.687}\right\}\\
&<90.
\end{split}
\end{equation}
Indeed, the two entries in the maximum bound the output norm arising
from a unit input on $J$ and on $E$, respectively.
Now write $X=X_f\oplus X_t$, where $X_t$ is the infinite tail space.
The full equation $(I-\eta S)x=r$ is equivalent, after elimination of
$x_f$, to
\begin{equation}\label{interpolation:tail-schur}
\bigl(I-\eta S_{t,t}-S_{t,f}V_fS_{f,t}\bigr)x_t
=r_t+\eta S_{t,f}V_fr_f.
\end{equation}
Put $\delta=5\cdot10^{-10}$. By
\eqref{interpolation:tails},
\[
\norm{S_{t,t}}<\delta,\qquad
\norm{S_{t,f}}<\delta,\qquad
\norm{S_{f,t}}<28.
\]
The operator subtracted from the identity in
\eqref{interpolation:tail-schur} thus has norm less than
\[
\delta(1+90\cdot28)=2521\delta<1.
\]
Its Neumann series converges in the bounded operators on $X_t$.
Together with the finite inverse this proves invertibility on the full
space $X$, without any passage to a limit of finite-section solutions.
Quantitatively,
\[
\norm{x_t}_1
\le\frac{\norm{r_t}_1+90\delta\norm{r_f}_1}{1-2521\delta},
\qquad
\norm{x_f}_1\le90\norm{r_f}_1+2520\norm{x_t}_1.
\]
Adding these inequalities gives coefficients less than $91$ and $2540$
for $\norm{r_f}_1$ and $\norm{r_t}_1$, respectively, which proves
\eqref{interpolation:inverse-bound}.
\end{proof}
\begin{proposition}[Exact interpolation and coefficient enclosure]
\label{interpolation:solution}
There are unique real lists $x_1,x_2\in X$ satisfying
\eqref{fourier:coefficient-equations} with the fixed data
\eqref{fourier:fixed}. For each $i=1,2$,
\begin{equation}\label{interpolation:coefficient-error}
\sum_{n\in\mathcal N}
\bigl(\abs{c_{i,n}-\widetilde c_{i,n}}
+\abs{d_{i,n}-\widetilde d_{i,n}}\bigr)<.00003.
\end{equation}
The norm of each exact list, including its node-zero coefficients, is
less than $2.3$. The corresponding functions satisfy all the
interpolation conditions \eqref{fourier:interpolation}.
\end{proposition}
The uniqueness here is within the specified coefficient space, with the
fixed node-zero data and constants. Triangular-shell values and first
derivatives do not determine an arbitrary radial Schwartz function and
its Fourier transform uniquely \cite[Theorem~1.7]{Sardari2021}.
\begin{proof}
Let $p_i^0$ denote \eqref{fourier:cardinal} with the fixed values
$(C_i,c_{i,0},d_{i,0})$ and all positive-node coefficients zero. Let
$e\in X$ have its $d_1$ coordinate equal to one and all other coordinates
zero. Define $g_i\in X$ by
\[
(g_i)_{c,n}=-\frac{K_{p_{3-i}^0}(n)}{Q_n},\qquad
(g_i)_{d,n}=
\frac{D_nK_{p_{3-i}^0}(n)-K_{p_{3-i}^0}'(n)}{Q_n}
-\mathbf1_{\{i=1\}}e_{d,n}.
\]
The estimates for the node-zero and constant columns show that these
are summable lists. The equations to solve are precisely
\[
x_1-Sx_2=g_1,\qquad x_2-Sx_1=g_2.
\]
Their sum is an equation for $I-S$, and their difference is an equation
for $I+S$. Lemma~\ref{interpolation:inverse} gives unique real solutions
and hence unique $x_1,x_2$.
It remains to locate these solutions relative to the table. Set
\[
\rho_i=\widetilde x_i-S\widetilde x_{3-i}-g_i.
\]
There are $102$ scalar coordinates in $I_f$, so
\eqref{interpolation:finite-residuals} gives
\[
\norm{(\rho_i)_f}_1<1.1\cdot10^{-7}=:R_f.
\]
On the tail, both the tabulated positive-node list and the vector $e$
vanish. The density-column estimate, now including node zero, and
\eqref{interpolation:atom-tail} imply
\[
\norm{(\rho_i)_t}_1
<2.29\cdot5\cdot10^{-10}+.017\cdot3.2\cdot10^{-8}=:R_t.
\]
Here $2.29$ bounds the whole tabulated coefficient list, and $.017$
bounds both fixed constants in absolute value.
The errors $x_i-\widetilde x_i$ are recovered by taking half the sum
and half the difference of the solutions with residuals
$-(\rho_1+\rho_2)$ and $-(\rho_1-\rho_2)$. Applying
\eqref{interpolation:inverse-bound} to both gives, for each list,
\[
\norm{x_i-\widetilde x_i}_1
\le2(91R_f+2540R_t)
=.00002860012<.00003.
\]
This proves \eqref{interpolation:coefficient-error}. Since the fixed
node-zero coefficients agree, the total exact coefficient norm is
less than $2.29+.00003<2.3$. Finally,
\eqref{fourier:coefficient-equations} is equivalent to
\eqref{fourier:interpolation} by Lemma~\ref{fourier:equivalence}.
\end{proof}
==== END A-sharp-Fourier-certificate-for-planar-circle-packing//interpolation.tex ====
==== BEGIN A-sharp-Fourier-certificate-for-planar-circle-packing//introduction.tex ====
\section{Introduction}
\label{sec:introduction}
The planar circle-packing problem asks what fraction of the plane can be
occupied by congruent disks with disjoint interiors. Scale their radius
to $1/2$, so their centers form a set $\mathcal C\subset\RR^2$ with
$|x-y|\ge1$ for distinct centers. If $B_R$ is the disk of radius $R$
centered at the origin, we use the upper covered-area density
\[
\overline\delta(\mathcal C)
=\limsup_{R\to\infty}
\frac{\vol\bigl(B_R\cap\bigcup_{c\in\mathcal C}B(c,1/2)\bigr)}
{\vol B_R}.
\]
The triangular lattice attains $\pi/(2\sqrt3)$, and Thue's theorem gives
the matching upper bound. Hales presents an elementary geometric proof
based on an idea of Rogers~\cite[pp.~442--443]{Hales2000}.
The Fourier-analytic approach of Cohn and Elkies~\cite{CohnElkies2003}
asks whether the same optimum can be certified by one function with
compatible signs before and after Fourier transformation. The density
theorem alone does not supply such a function. Our result is a sharp
Fourier certificate for the classical optimum.
Throughout, the Fourier transform is
\[
\widehat f(\xi)=\int_{\RR^2}f(x)e^{-2\pi i x\cdot\xi}\,\dd x.
\]
For real radial integrable functions it is real and radial. A Schwartz
function is a smooth function whose derivatives decay faster than every
inverse power of the radius; this regularity is stronger than the
admissibility conditions in the Cohn--Elkies theorem.
\begin{theorem}[Sharp Fourier certificate]\label{intro:main}
There is a real radial Schwartz function $f:\RR^2\to\RR$ such that
\begin{equation}\label{intro:conditions}
\widehat f(0)=1,\qquad f(0)=\frac{2}{\sqrt3},\qquad
\widehat f(\xi)\ge0\quad(\xi\in\RR^2),\qquad
f(x)\le0\quad(|x|\ge1).
\end{equation}
\end{theorem}
Theorem~\ref{intro:main} gives an affirmative answer to the
two-dimensional case of Cohn and Elkies's Conjecture~7.3
\cite{CohnElkies2003}, concerning sharp auxiliary functions for their
linear programming bound. Their Definition~2.2 and Theorem~3.1 apply to
this Schwartz function and to arbitrary packings, not only periodic ones.
Their opposite Fourier phase gives the same transform for these even
functions.
For disks of radius $1/2$ they give
\[
\vol B(0,1/2)\,\frac{f(0)}{\widehat f(0)}
=\frac{\pi}{4}\frac{2}{\sqrt3}
=\frac{\pi}{2\sqrt3}.
\]
This is the density of the lattice generated by $(1,0)$ and
$(1/2,\sqrt3/2)$. Thus the theorem identifies an extremizer for this
analytic bound. To put it in the equal-origin formulation of
Cohn and Elkies's Theorem~3.2, which their Conjecture~7.3 invokes, set
\[
g(x)=\frac{\sqrt3}{2}\,
f\bigl((\sqrt3/2)^{1/2}x\bigr).
\]
The two-dimensional scaling rule gives $g(0)=\widehat g(0)=1$,
$\widehat g\ge0$, and $g(x)\le0$ for
$|x|\ge(2/\sqrt3)^{1/2}$, exactly the requested planar witness.
The sign estimates also determine the zeros in their required domains.
For $f$, every zero with $|x|\ge1$ has integral squared radius. This
restriction recovers a classical equality consequence: among periodic
packings, equality in the density bound forces the triangular lattice
(Corollary~\ref{normalization:periodic-uniqueness}). The equal-origin
witness $g$ has zeros outside the shells of the triangular lattice of
covolume one, so it does not meet the additional zero-set requirement of
Cohn and Elkies's Conjecture~8.1. Section~\ref{sec:fourier} gives an explicit
extra node, and Section~\ref{sec:normalization} records the zero sets in
their respective sign domains and final radial coordinates.
\subsection{Fourier certificates and interpolation}
Cohn and Elkies formulated their sharpness conjecture in dimensions
$2$, $8$, and $24$~\cite[Conjecture~7.3]{CohnElkies2003}.
Their finite Laguerre--Gaussian approximations already combine prescribed
double roots, numerical searches for the remaining parameters, and rigorous
sign checks by rational arithmetic and Sturm root counts
\cite[Section~7, pp.~704--705]{CohnElkies2003}. Exact attainment of the
conjectured sharp bound remained unresolved in that approach.
Viazovska constructed the required function in dimension
$8$~\cite[Theorems~3 and~4]{Viazovska2017}, and Cohn, Kumar, Miller,
Radchenko, and Viazovska constructed its dimension-$24$ counterpart
\cite[Section~4]{CKMRV2017}.
Their constructions use modular forms to control a function and its
Fourier transform simultaneously.
Cohn and Elkies report that Gorbachev independently obtained essentially
the same linear programming bound \cite[p.~694]{CohnElkies2003}.
Interpolation gives a broader setting for these constructions.
Radchenko and Viazovska recovered even Schwartz functions on the real
line from their values and Fourier values at the square roots of
nonnegative integers~\cite[Theorem~1]{RadchenkoViazovska2019}.
Cohn, Kumar, Miller, Radchenko, and Viazovska developed interpolation
using values and first derivatives in dimensions $8$ and $24$ and used
it to prove universal optimality of the $E_8$ and Leech lattices
\cite[Theorems~1.4, 1.7, and~1.9]{CKMRV2022}.
The planar shell data have a different feature.
Cohn, Kumar, Miller, Radchenko, and Viazovska conjectured that values
and first derivatives of a radial Schwartz function and its Fourier
transform on the nonzero shells of the triangular lattice normalized to
covolume one do not determine
the function uniquely~\cite[Conjecture~7.5]{CKMRV2022}.
Sardari proved this nonuniqueness by constructing infinitely many
linearly independent radial Schwartz functions with double zeros,
before and after Fourier transformation, on those shells
\cite[Theorem~1.7]{Sardari2021}.
Thus prescribing the shell data alone does not single out a planar
certificate. We solve a specified interpolation problem within an
explicit family of Gaussian-damped cardinal functions and establish
the signs separately.
A 2026 bachelor's thesis abstract by Zhitniaia~\cite{Zhitniaia2026}
announces an integral representation of a magic function for the
hexagonal lattice using weakly holomorphic modular forms and theta
functions for $\Gamma_0(12)$. Our comparison is limited to the
institutional abstract, which does not specify the complete sign,
normalization, and regularity statement. The precise relation of that
claim to Theorem~\ref{intro:main} remains to be determined from the
full thesis.
\subsection{The construction}
For a sharp Fourier bound on a lattice, Poisson summation forces
the auxiliary function and its Fourier transform to vanish on the
nonzero shells of the lattice and its dual, respectively
\cite[Section~5]{CohnElkies2003}.
After normalizing the triangular lattice to covolume one, it is
isometric to its dual, so the two sets of shell radii coincide.
This normalization makes the coupled Fourier equations symmetric.
At the end we dilate the resulting function to put the first
physical-space shell at the exclusion radius one in
Theorem~\ref{intro:main}.
We cover their squared radii, in a fixed scaled coordinate, by a
periodic set of interpolation nodes.
Periodic supersets of the values of the triangular quadratic form
also occur in Sardari's construction
\cite[Section~6.1]{Sardari2021}; here a fixed periodic set permits
an explicit trigonometric cardinal family.
An interior zero of a differentiable function with a fixed sign must
also have zero first derivative. We therefore impose zeros of order
at least two on the Fourier side and at every physical-space shell
beyond the first. The first physical-space shell is the boundary of
the required sign region, so a simple zero is possible there; its
prescribed negative derivative will fix the normalization.
A trigonometric polynomial vanishes doubly on the chosen set.
Dividing it by linear and quadratic factors produces functions
that prescribe values and first derivatives at individual nodes.
This Hermite-cardinal pattern has a precedent in the Gaussian
Beurling--Selberg constructions of Carneiro, Littmann, and Vaaler
\cite[formulas~(1.2)--(1.3)]{CarneiroLittmannVaaler2013}; the Fourier
coupling and sign estimates below are proved for the present family.
Multiplication by a fixed Gaussian gives radial Schwartz functions.
Their Fourier transforms admit integral representations with
additional Gaussian decay. We use that decay to solve the coupled
interpolation equations: a finite matrix handles the first nodes,
and explicit estimates control the remaining infinite system. The inversion
acts on the full unweighted $\ell^1$ coefficient space. Its uniqueness
concerns the prescribed coefficient lists within this chosen family.
The resulting functions have exact zeros, but the coefficient
lists used for computation are finite approximations.
An explicit per-list $\ell^1$ error bound connects these approximations
to the exact solution. Near a zero, a raw error divided by its vanishing
power would be singular. We instead subtract the lower Taylor terms of
the error before division; positive Bernstein coefficients then certify
the signs on complete half-gaps. A separate sine-product barrier and
second-derivative estimate control the entire unbounded region.
Poisson summation first equates the two origin values.
Differentiating its dilation identity then determines their common
positive value from the prescribed first-shell derivative, so the
final normalization is exact.
Section~\ref{sec:fourier} develops the cardinal functions and their
Fourier transforms. Section~\ref{sec:interpolation} solves the
interpolation equations, and Section~\ref{sec:signs} proves the
global signs and identifies their zeros. Section~\ref{sec:normalization}
establishes the exact normalization, performs the final dilation, and
derives the packing consequences. Appendix~\ref{sec:certificates} gives
the finite data and the primary interval certificate, including the
analytic integration error that connects its finite sums to exact moments.
Appendix~\ref{ind:section} gathers the rational alternatives and
independent checks, with their precise mathematical and computational
scopes.
==== END A-sharp-Fourier-certificate-for-planar-circle-packing//introduction.tex ====
==== BEGIN A-sharp-Fourier-certificate-for-planar-circle-packing//main.tex ====
\documentclass[11pt]{article}
\ifdefined\pdfinfoomitdate\pdfinfoomitdate=1\fi
\ifdefined\pdftrailerid\pdftrailerid{}\fi
\ifdefined\pdfsuppressptexinfo\pdfsuppressptexinfo=15\fi
\ifdefined\pdfglyphtounicode
\input{glyphtounicode}
\pdfglyphtounicode{parenleftbig}{0028}
\pdfglyphtounicode{parenleftBig}{0028}
\pdfglyphtounicode{parenleftbigg}{0028}
\pdfglyphtounicode{parenleftBigg}{0028}
\pdfglyphtounicode{parenrightbig}{0029}
\pdfglyphtounicode{parenrightBig}{0029}
\pdfglyphtounicode{parenrightbigg}{0029}
\pdfglyphtounicode{parenrightBigg}{0029}
\pdfglyphtounicode{braceleftbig}{007B}
\pdfglyphtounicode{braceleftbigg}{007B}
\pdfglyphtounicode{bracerightbig}{007D}
\pdfglyphtounicode{bracerightbigg}{007D}
\pdfglyphtounicode{summationtext}{2211}
\pdfglyphtounicode{summationdisplay}{2211}
\pdfglyphtounicode{integraltext}{222B}
\pdfglyphtounicode{integraldisplay}{222B}
\pdfglyphtounicode{productdisplay}{220F}
\pdfglyphtounicode{uniontext}{22C3}
\pdfglyphtounicode{hatwide}{005E}
\pdfglyphtounicode{hatwider}{005E}
\pdfglyphtounicode{tildewide}{007E}
\pdfgentounicode=1
\fi
\usepackage[OT1]{fontenc}
\usepackage[utf8]{inputenc}
\usepackage[margin=1in]{geometry}
\usepackage{amsmath,amssymb,amsthm,mathtools}
\usepackage{booktabs,longtable,array}
\usepackage{listings}
\usepackage{microtype}
\usepackage{xcolor}
\usepackage[numbers,sort&compress]{natbib}
\usepackage[colorlinks=true,linkcolor=blue!45!black,citecolor=blue!45!black,urlcolor=blue!45!black]{hyperref}
\hypersetup{pdftitle={A sharp Fourier certificate for planar circle packing},pdfauthor={OpenAI}}
\newcommand{\RR}{\mathbb R}
\newcommand{\CC}{\mathbb C}
\newcommand{\ZZ}{\mathbb Z}
\newcommand{\NN}{\mathbb N}
\newcommand{\dd}{\mathop{}\!\mathrm d}
\newcommand{\norm}[1]{\left\lVert#1\right\rVert}
\newcommand{\abs}[1]{\left\lvert#1\right\rvert}
\DeclareMathOperator{\supp}{supp}
\DeclareMathOperator{\vol}{vol}
\DeclareMathOperator{\Rea}{Re}
\DeclareMathOperator{\Ima}{Im}
\newtheorem{theorem}{Theorem}[section]
\newtheorem{proposition}[theorem]{Proposition}
\newtheorem{lemma}[theorem]{Lemma}
\newtheorem{corollary}[theorem]{Corollary}
\theoremstyle{definition}
\newtheorem{definition}[theorem]{Definition}
\theoremstyle{remark}
\newtheorem{remark}[theorem]{Remark}
\numberwithin{equation}{section}
\setlength{\emergencystretch}{1.5em}
\lstset{basicstyle=\ttfamily\footnotesize,columns=fullflexible,keepspaces=true,
showstringspaces=false,breaklines=true,frame=single,framerule=.2pt}
\title{A sharp Fourier certificate for planar circle packing}
\author{OpenAI}
\date{September 23, 2026}
\begin{document}
\maketitle
\begin{abstract}
We resolve the planar Cohn--Elkies sharpness conjecture: the two-point
Fourier bound attains the optimal circle-packing density $\pi/(2\sqrt3)$.
We construct a radial Schwartz certificate and prove its global sign
conditions using rigorous interval arithmetic and analytic estimates.
The same certificate recovers the classical uniqueness of the triangular
packing among periodic equality cases.
\end{abstract}
\begingroup
\small
\tableofcontents
\endgroup
\clearpage
\input{sections/introduction}
\input{sections/fourier}
\input{sections/interpolation}
\input{sections/signs}
\input{sections/normalization}
\appendix
\input{sections/certificates}
\input{sections/independent-certification}
\bibliographystyle{plainnat}
\let\packingbibsection\bibsection
\renewcommand{\bibsection}{%
\packingbibsection
\addcontentsline{toc}{section}{\refname}%
}
\bibliography{references}
\end{document}
==== END A-sharp-Fourier-certificate-for-planar-circle-packing//main.tex ====
==== BEGIN A-sharp-Fourier-certificate-for-planar-circle-packing//normalization.tex ====
\section{Poisson normalization}
\label{sec:normalization}
The prescribed first-shell derivative determines the remaining scalar
normalization. Poisson summation first equates the values at the origin;
its differentiated dilation identity then proves that their common value is
positive, without estimating that value numerically.
Let the exact coefficient lists be those of
Proposition~\ref{interpolation:solution}, and put
\[
G_i(s)=e^{-\pi hs}H_i(s),\qquad F(x)=G_1(b|x|^2).
\]
Section~\ref{sec:fourier} proves that $F$ is real, radial, and Schwartz,
with $\widehat F(\xi)=G_2(b|\xi|^2)$. The interpolation conditions give
\begin{equation}\label{normalization:jets}
\begin{split}
G_1(n)&=G_2(n)=0,\\
G_2'(n)&=0,\qquad
G_1'(n)=-Q_1e^{-\pi h}\mathbf1_{\{n=1\}}
\qquad(n\in\mathcal N).
\end{split}
\end{equation}
Indeed, $G_i'(n)=e^{-\pi hn}H_i'(n)$ because $H_i(n)=0$.
\subsection{The dilation identity}
Use the lattice $\Lambda$ of covolume one and the shell coordinate $s_a$ from
\eqref{fourier:lattice}--\eqref{fourier:shell-coordinate}. Every positive
$s_a$ lies in $\mathcal N$. Exactly six lattice vectors have
$s_a=1$: from
$4(j^2+jk+k^2)=(2j+k)^2+3k^2=4$, one obtains
\[
(j,k)=(1,0),\ (-1,0),\ (0,1),\ (0,-1),\ (1,-1),\ (-1,1).
\]
As shown in Section~\ref{sec:fourier}, $\Lambda^*=\mathcal R\Lambda$.
Apply Poisson summation on $\Lambda$ to
$F_t(x)=F(\sqrt t\,x)$, whose Fourier transform is
$t^{-1}\widehat F(\xi/\sqrt t)$. Radiality makes the sums over
$\Lambda^*$ and $\Lambda$ equal, and gives, for every $t>0$,
\begin{equation}\label{normalization:poisson}
\sum_{a\in\Lambda}G_1(ts_a)
=\frac1t\sum_{a\in\Lambda}G_2(s_a/t).
\end{equation}
We will differentiate this identity twice. Here is the needed convergence
justification. On a compact subinterval of $(0,\infty)$, each of the first
two $t$ derivatives of $F(\sqrt t\,x)$ is a finite sum of spatial derivatives
of $F$ multiplied by polynomials in $x$ and bounded functions of $t$.
The same statement holds for $t^{-1}\widehat F(x/\sqrt t)$.
The Schwartz bounds therefore imply, for every $M$, a common bound
$C_M(1+|x|)^{-M}$ for these functions and their first two $t$ derivatives.
The number of points of a fixed planar lattice in a radius-$R$ disk is
$O(R^2)$, so $\sum_{a\in\Lambda}(1+|a|)^{-M}$ converges for $M>2$.
Uniform convergence on compact $t$ intervals now permits both termwise
differentiations of \eqref{normalization:poisson}.
\begin{lemma}[Exact origin value]\label{normalization:origin-lemma}
The Fourier pair constructed above satisfies
\begin{equation}\label{normalization:origin}
F(0)=\widehat F(0)=6Q_1e^{-\pi h}>0.
\end{equation}
\end{lemma}
\begin{proof}
At $t=1$, all nonzero summands in \eqref{normalization:poisson} vanish by
\eqref{normalization:jets}. Hence $G_1(0)=G_2(0)$.
The derivative of its left side at $1$ is
\[
\sum_{a\in\Lambda}s_aG_1'(s_a)=-6Q_1e^{-\pi h},
\]
because only the six vectors on the first shell contribute. On the right,
\[
\left.\frac{\dd}{\dd t}\left(t^{-1}G_2(s/t)\right)\right|_{t=1}
=-G_2(s)-sG_2'(s).
\]
For $s=s_a>0$, both terms vanish by the exact jets, while $s_a=0$ contributes
$-G_2(0)$. Equality of the derivatives gives
$G_2(0)=6Q_1e^{-\pi h}$. The positivity follows from $Q_1>0$ in
\eqref{fourier:node-bounds}. This calculation uses the Fourier pairing and
the exact interpolation conditions, but not the global sign estimates.
\end{proof}
\subsection{The final dilation}
\begin{proof}[Proof of Theorem~\ref{intro:main}]
Define
\begin{equation}\label{normalization:scaling}
f(x)=\frac{F(x/\sqrt b)}{bF(0)},\qquad
\widehat f(\xi)=\frac{\widehat F(\sqrt b\,\xi)}{F(0)}.
\end{equation}
The second identity is the two-dimensional Fourier scaling rule. Since
$F(0)>0$, the dilation and normalization preserve real-valuedness,
radiality, and Schwartz regularity. They give
$\widehat f(0)=1$ and $f(0)=1/b=2/\sqrt3$.
By Proposition~\ref{signs:global} and the positive factor $e^{-\pi hs}$,
$G_2(s)=e^{-\pi hs}H_2(s)\ge0$ for $s\ge0$ and
$G_1(s)=e^{-\pi hs}H_1(s)\le0$ for $s\ge1$.
Thus $\widehat f\ge0$ everywhere. If $|x|\ge1$,
then $b|x/\sqrt b|^2=|x|^2\ge1$, so $f(x)\le0$. These are exactly the
four assertions of \eqref{intro:conditions}.
\end{proof}
The scaling also records the prescribed contact data in the final radial
coordinates. Write $f_{\mathrm{rad}}(r)=f(r,0)$ and
$\widehat f_{\mathrm{rad}}(\rho)=\widehat f(\rho,0)$ for $r,\rho\ge0$. Then
\[
f_{\mathrm{rad}}(r)=\frac{G_1(r^2)}{bF(0)},\qquad
\widehat f_{\mathrm{rad}}(\rho)=\frac{G_2(b^2\rho^2)}{F(0)}.
\]
Proposition~\ref{signs:global} therefore gives the exact zero sets in the
respective sign domains:
\begin{equation}\label{normalization:zero-sets}
\begin{aligned}
\{r^2:r\ge1,\ f_{\mathrm{rad}}(r)=0\}&=\mathcal N,\\
\{\rho^2:\rho\ge0,\ \widehat f_{\mathrm{rad}}(\rho)=0\}
&=b^{-2}\mathcal N=\frac43\mathcal N.
\end{aligned}
\end{equation}
The physical zero at $r=1$ is simple, every other listed physical zero is
exactly double, and every listed Fourier zero is exactly double. Indeed,
the Gaussian and scalar factors do not vanish, while the changes of
variable $s=r^2$ and $s=b^2\rho^2$ have nonzero derivatives at these positive
radii, so they preserve the orders established for $H_1,H_2$.
At the exclusion radius, the radial derivative is
\begin{equation}\label{normalization:radial-slope}
f_{\mathrm{rad}}'(1)
=\frac{2G_1'(1)}{bF(0)}=-\frac{1}{3b}=-\frac{2}{3\sqrt3}.
\end{equation}
\subsection{The packing consequence}
The relation between the signs and packing density is particularly visible
for a periodic packing. Suppose the center set is
$\bigcup_{j=1}^N(L+t_j)$, where $L$ has covolume $V$, the $t_j$ are distinct
modulo $L$, and distinct centers have distance at least one. Absolute
convergence follows from the Schwartz bounds. Every summand below except
the $N$ diagonal terms at the origin is nonpositive, so
\[
\mathcal S:=\sum_{j,k=1}^N\sum_{v\in L}f(v+t_j-t_k)\le Nf(0).
\]
Poisson summation on each translated lattice gives
\[
\mathcal S=\frac1V\sum_{w\in L^*}\widehat f(w)
\left|\sum_{j=1}^N e^{2\pi i w\cdot t_j}\right|^2
\ge\frac{N^2}{V}\widehat f(0),
\]
because every Fourier summand is nonnegative. It follows that
$N/V\le f(0)/\widehat f(0)=2/\sqrt3$. The disks have area $\pi/4$, hence
their covered-area density is at most $\pi/(2\sqrt3)$.
For an arbitrary, not necessarily periodic, packing, the same implication
is Theorem~3.1 of Cohn and Elkies~\cite{CohnElkies2003}. Its admissibility
hypotheses hold because both $f$ and $\widehat f$ are Schwartz.
For disks of radius $1/2$, the Cohn--Elkies center-density normalization is
one quarter of the center-number density in their upper-density convention.
Multiplying it by the unit-disk area $\pi$ therefore gives the corresponding
covered-area density, equivalently $\pi/4$ times that center-number density.
Their upper density takes a supremum over observation centers, so it bounds
the fixed-origin limsup used in the introduction.
For the fixed-origin convention and separated centers, the
covered area inside $B_R$ lies between $(\pi/4)$ times the number of centers
in $B_{R-1/2}$ and $(\pi/4)$ times the number in $B_{R+1/2}$; division by
$\pi R^2$ and passage to the limsup matches the convention in the
introduction. The lattice generated by $(1,0)$ and $(1/2,b)$ has one center
per area $b$ and attains this value. This is the classical packing-density
theorem recovered from the new Fourier certificate, not a new density claim.
The restriction of the physical zero set to integer squared radii also
recovers the classical periodic equality case. We use the even integral
lattice argument of Cohn and Elkies~\cite[Section~8]{CohnElkies2003}.
\begin{corollary}[Periodic equality case]\label{normalization:periodic-uniqueness}
If a periodic packing of disks of radius $1/2$ has covered-area density
$\pi/(2\sqrt3)$, then its center set is the image of
$\ZZ(1,0)+\ZZ(1/2,b)$ under a Euclidean isometry.
\end{corollary}
\begin{proof}
Write the center set as
$\mathcal C=\bigcup_{j=1}^N(L+t_j)$, with the notation of the periodic
packing argument above. Equality of its density with $\pi/(2\sqrt3)$ gives
$N/V=f(0)/\widehat f(0)$. Thus both inequalities in the Poisson comparison
are equalities. Absolute convergence permits us to write the defect in
the physical inequality as
\[
0=Nf(0)-\mathcal S
=\sum_{\substack{1\le j,k\le N,\ v\in L\\
v\ne0\ \text{or}\ j\ne k}}
-f(v+t_j-t_k).
\]
Every summand is nonnegative, so every one is zero. Any difference of
distinct centers occurs in this sum, and \eqref{normalization:zero-sets}
therefore gives
\[
|x-y|^2\in\mathcal N\subset\ZZ
\qquad(x,y\in\mathcal C,\ x\ne y).
\]
Only the physical zero set is used here.
Translate one center to the origin and dilate by $\sqrt2$. Put
\[
L_0=\sqrt2\,L,\qquad
\mathcal C_0=\sqrt2\,(\mathcal C-t_1),\qquad
\Gamma=\operatorname{span}_{\ZZ}\mathcal C_0.
\]
Here $\Gamma$ consists of finite integer linear combinations of the
scaled centers, and $L_0\subset\mathcal C_0$ because the first coset becomes
$L_0$. All squared distances in $\mathcal C_0$ are even integers. Since
$0\in\mathcal C_0$, polarization gives
\[
\langle x,y\rangle
=\frac{|x|^2+|y|^2-|x-y|^2}{2}\in\ZZ
\qquad(x,y\in\mathcal C_0).
\]
Integer linear combinations preserve these integral pairings, and their
squared norms remain even because the cross terms carry a factor of two.
Thus every pairing in $\Gamma$ is integral and every squared norm in
$\Gamma$ is even.
In particular, pairing with a basis of $L_0$ shows
$L_0\subset\Gamma\subset L_0^*$. Both $L_0$ and $L_0^*$ are lattices of
rank two, so $L_0^*/L_0$ is finite. Hence $\Gamma$ is a finite union of
$L_0$ cosets and is itself a lattice of rank two. Put
$d:=[\Gamma:L_0]$. A basis of $\Gamma$ has a Gram matrix
\[
\begin{pmatrix}a&m\\m&c\end{pmatrix},
\qquad a,c\in2\ZZ_{>0},\quad m\in\ZZ.
\]
Its determinant $ac-m^2=\operatorname{covol}(\Gamma)^2$ is a positive
integer congruent to $0$ or $3$ modulo $4$. It is therefore at least $3$,
and $\operatorname{covol}(\Gamma)\ge\sqrt3$.
On the other hand, $\mathcal C_0$ consists of $N$ distinct $L_0$ cosets,
all among the $d$ cosets of $\Gamma$. Since
$\operatorname{covol}(L_0)=2V$, we obtain
\[
\sqrt3\le\operatorname{covol}(\Gamma)
=\frac{2V}{d}\le\frac{2V}{N}=\sqrt3.
\]
Equality holds throughout, so $d=N$ and the two sets of cosets coincide:
$\mathcal C_0=\Gamma$.
It remains to identify this lattice of Gram determinant $3$. Choose a
shortest nonzero vector $w_1\in\Gamma$ and extend it to a lattice basis
$w_1,w_2$; this is possible because a shortest vector is primitive.
Put $a=|w_1|^2$, and replace $w_2$ by an integer translate along $w_1$
so that $m=\langle w_1,w_2\rangle$ satisfies $|m|\le a/2$.
With $c=|w_2|^2\ge a$ and the even integer $a\ge2$,
\[
3=ac-m^2\ge\frac34a^2.
\]
Hence $a=2$. The equation $2c-m^2=3$, with $c$ even and $|m|\le1$,
forces $c=2$ and $m=\pm1$. This is the Gram matrix of the triangular
lattice with shortest length $\sqrt2$. Scaling back by $1/\sqrt2$ proves
the assertion.
\end{proof}
\subsection{A second-derivative consistency identity}
The same dilation identity contains an additional relation among the exact
curvatures. It is not needed for the construction or its sign proof.
For $n\in\mathcal N$, let
$r_\Lambda(n)=\#\{a\in\Lambda:s_a=n\}$, allowing the value zero for nodes
that are not lattice radii; in particular, $r_\Lambda(1)=6$.
\begin{proposition}[Second differentiated Poisson identity]
\label{normalization:second-poisson}
The exact pair satisfies the absolutely convergent identity
\begin{equation}\label{normalization:curvature-identity}
\sum_{n\in\mathcal N}r_\Lambda(n)n^2
\bigl(G_1''(n)-G_2''(n)\bigr)=2G_2(0).
\end{equation}
Using the global signs, it follows that
\begin{equation}\label{normalization:curvature-inequality}
H_1''(1)-H_2''(1)\ge2(1-\pi h)Q_1.
\end{equation}
\end{proposition}
\begin{proof}
The convergence argument above permits a second differentiation of
\eqref{normalization:poisson}. At $t=1$,
\[
\left.\frac{\dd^2}{\dd t^2}\left(t^{-1}G_2(s/t)\right)\right|_{t=1}
=2G_2(s)+4sG_2'(s)+s^2G_2''(s).
\]
On the left the corresponding term is $s^2G_1''(s)$. At every positive
node, the first two terms on the right vanish by
\eqref{normalization:jets}; at the origin they contribute $2G_2(0)$.
Grouping the remaining absolutely convergent lattice sums by their squared
radii proves \eqref{normalization:curvature-identity}.
For a positive node $n>1$, $G_1$ has a local maximum equal to zero and
$G_2$ has a local minimum equal to zero, by Proposition~\ref{signs:global}
and $e^{-\pi hs}>0$.
Therefore $G_1''(n)\le0\le G_2''(n)$. All terms in
\eqref{normalization:curvature-identity} beyond the first shell are
nonpositive. Since $r_\Lambda(1)=6$ and $G_2(0)=6Q_1e^{-\pi h}$,
\[
G_1''(1)-G_2''(1)\ge2Q_1e^{-\pi h}.
\]
Finally, differentiating $G_i=e^{-\pi hs}H_i$ twice and inserting the
exact value and slope at $1$ gives
\[
G_1''(1)-G_2''(1)
=e^{-\pi h}\bigl(H_1''(1)-H_2''(1)+2\pi hQ_1\bigr).
\]
This yields \eqref{normalization:curvature-inequality}. Unlike the origin
calculation, this last inequality is conditional on the global signs that
were proved in Section~\ref{sec:signs}.
\end{proof}
==== END A-sharp-Fourier-certificate-for-planar-circle-packing//normalization.tex ====
==== BEGIN A-sharp-Fourier-certificate-for-planar-circle-packing//rational-check.tex ====
\subsection{An independent rational check of the first block}
\label{cert:rational-small}
The optional program \path{verification/check_sjj_rational.py} checks
the $6\times6$ block on $J=\{1,3,4\}$ by exact rational arithmetic.
It uses the cosine rule \eqref{cert:rule} with $N=64$ and the sum over
$1\le a\le31$, and otherwise the rounded exponential, kernels and
basis densities specified above. There are no atom terms for these
basis inputs. The following error bound makes this smaller check
independent of the full Arb calculation.
\begin{lemma}\label{cert:rational-small-error}
For each basis input $c_n=1$ or $d_n=1$, $n\in J$, with all other data
zero, and each output node $m\in J$, the rational program computes
both entries of $\mathcal S_m(p)$ with absolute error less than
$10^{-12}$.
\end{lemma}
\begin{proof}
We separate quadrature error from rounding in the finite sum. Use the
same complex disks of radius $1/6$ as in
Lemma~\ref{cert:moment-error}. The basis density on any disk is bounded by
\[
D_0=2\pi^2\cdot10\cdot\frac{13}{12}e^{4\pi/6}.
\]
This also bounds a simple-pole basis density. The bounds
$|\lambda|<5$, $|z|<6.2$ and $\operatorname{Im}z>-.081$ proved there
show that the transformed value and first-derivative integrands are
bounded respectively by
\[
M_0=5D_0e^{4\pi\cdot.081}<24028,
\qquad M_1=6.2\pi M_0<468015.
\]
Use the rational upper bound $e^x\le\mathcal E_n(x)$ from
\eqref{ind:positive-exponential}. Put
$\overline\pi=p_*+10^{-36}>\pi$. Here
$x=743\overline\pi/750$ lies in $[0,4)$, so the formula applies with $n=40$.
Using $M_0=(325/3)\pi^2e^{743\pi/750}$, the exact rational comparisons
\[
\frac{325}{3}\overline\pi^2 \mathcal E_{40}(743\overline\pi/750)<24028,
\qquad
\frac{31}{5}\overline\pi\,
\frac{325}{3}\overline\pi^2 \mathcal E_{40}(743\overline\pi/750)<468015
\]
prove the displayed bounds for $M_0,M_1$. Each comparison is a fraction
inequality after clearing positive denominators; no rounded decimal
evaluation of an exponential is used.
Discrete cosine orthogonality now gives exactness through degree $63$.
The weights sum to one and each is at least $1/(64\cdot63)>0$ by
the same telescoping identity. Cauchy's estimate on the half-radius
real segments bounds the degree-$63$ Taylor tails by $M_i2^{-63}$.
The integral and quadrature have total mass one, after combining all
six pieces. Thus their discrepancies in the two moments are at most
\[
\epsilon_0=2^{-62}M_0<5.211\cdot10^{-15},\qquad
\epsilon_1=2^{-62}M_1<1.015\cdot10^{-13}.
\]
Using $Q_m>1.75$ and $|D_m|<2.12$, the quadrature error in either
entry of \eqref{cert:pair} is consequently less than
\begin{equation}\label{cert:rational-small-quadrature}
\frac{2.12\epsilon_0+\epsilon_1}{1.75}<6.432\cdot10^{-14}.
\end{equation}
For completeness, two changes in the rational implementation require
care in reusing the rounding estimates. First, for a real argument $u$
the program evaluates $e_*(r)$ with $r=u\bmod2\in[0,2)$.
Periodicity is used for the exact target:
$e^{i\pi r}=e^{i\pi u}$. The approximation $p_*$ is not treated as
giving an exactly periodic exponential. Since $r$ is real and
$|p_*r|<2(3.142)<1000$, Lemma~\ref{cert:rounded-exponential}, followed by
$|e^{ip_*r}-e^{i\pi r}|\le2\cdot10^{-36}$, gives an absolute error
less than $2\cdot10^{-32}$ for each such replacement. When $u$
contains an approximate quadrature node, first compare at that
approximate argument and then use the first-difference bound for the
exact exponential. This avoids any assumption of continuity of the
modulo operation.
Second, with $N=64$, the sum of the absolute cosine coefficients in
one weight is less than one. Including the final floor, the weight
error is therefore at most
\[
\frac{2\cdot10^{-32}}N+10^{-40}
<\frac{10^{-29}}N.
\]
Summing the stronger first expression, rather than the relaxed bound on
the right, gives the total weight error
\begin{equation}\label{cert:rational-small-weight-total}
N\left(\frac{2\cdot10^{-32}}N+10^{-40}\right)
=2\cdot10^{-32}+64\cdot10^{-40}<3\cdot10^{-32}.
\end{equation}
The node error is less than $2\cdot10^{-32}/12<2\cdot10^{-33}$.
The exact nodes are at distance at least $1/98304$ from the endpoints
of their own interval $[(j-1)/6,j/6]$, since
$1-\cos(\pi/128)\ge2/128^2$. Thus the approximate nodes remain in their
own intervals, in particular in $[0,1]$. Their transformed frequencies
have positive imaginary part even after the floor, because
$\operatorname{Im}z\ge26/435$ there. The real-segment estimates in the
proof of Lemma~\ref{cert:rational-moment-error} also give $|z_*|<3$.
Thus every transformed call has its argument in the closed left
half-plane and satisfies $|p_*z_*m|<3.142(3)(4)<1000$.
Here is a direct propagation bound for these smaller inputs. The
real-segment estimates are $|\lambda|<3$, $|z|<3$, and basis-density
modulus less than $214$. For one piece put
\[
A_n=\sum_{l=j}^6P_le^{i\pi t_ln},\qquad
B_n=\sum_{l=j}^6t_lP_le^{i\pi t_ln}.
\]
Their moduli are at most ten. Replacing $b$ and the real exponentials
changes each sum by less than $3\cdot10^{-31}$; including the node
movement changes $B_n-tA_n$ by less than $7\cdot10^{-31}$.
The phase $e^{-i\pi tn}$ changes by less than $5\cdot10^{-32}$,
since $n\le4$. Substitution in \eqref{fourier:folded} therefore gives
a density error less than $3\cdot10^{-29}$.
The derivatives of the rational kernels, using $b,b_*>4/5$ and
$|t+ih|\ge2/5$, give errors less than $2\cdot10^{-32}$ in both
$\lambda$ and $z$, including their floors. The first-difference
formula for exponentials in the closed left half-plane, with $m\le4$,
then bounds the errors in
\[
g_0(t)=\lambda(t)e^{i\pi mz(t)},\qquad
g_1(t)=g_0(t)i\pi z(t)
\]
by $2\cdot10^{-30}$ and $3\cdot10^{-29}$, respectively.
The exact factors have moduli below $3$ and $36$; their approximations
have moduli below $4$ and $37$. Product differences bound each
derivative-integrand error by
\[
3\cdot10^{-29}\cdot37+214\cdot3\cdot10^{-29}
<8\cdot10^{-27};
\]
the value-integrand bound is smaller. By
\eqref{cert:rational-small-weight-total}, the total absolute weight error
is below $3\cdot10^{-32}$, and the approximate integrands have modulus
below $10^4$. Hence each rounded moment differs from its exact finite
quadrature moment by less than $10^{-26}$.
Finally, for residues $1,3,4$, the exact polynomial forms are
\[
\begin{gathered}
Q_1=Q_3=\frac{4\pi^2}{3}(1-b),\qquad Q_4=\frac{\pi^2}{3},\\
D_1=\frac{\pi(3+2b)}{18},\qquad
D_3=-\frac{\pi(3+2b)}{18},\qquad D_4=\frac{7\pi b}{9}.
\end{gathered}
\]
The program replaces $\pi,b$ by $p_*,b_*$ in these expressions.
Using $|p_*^2-\pi^2|<6.3\cdot10^{-36}$, $\pi^2<10$, and $b,b_*<1$ gives
\[
\begin{aligned}
|Q_m^\#-Q_m|
&<\max\{(4/3)(6.3+10),\,6.3/3\}10^{-36}<3\cdot10^{-35},\\
|D_m^\#-D_m|
&<\max\{(5+44/7)/18,\,29/9\}10^{-36}<4\cdot10^{-36}.
\end{aligned}
\]
The exact rational lower bound used
after~\eqref{ind:pair-stability} therefore shows that both denominators
$Q_m,Q_m^\#$ exceed $7/4$. The moments have modulus below $10^4$, so the
additional division and $D_m$ errors are less than $4\cdot10^{-31}$.
Together with \eqref{cert:rational-small-quadrature}, these bounds
give an error less than $7\cdot10^{-14}<10^{-12}$ in either matrix
entry, as asserted.
\end{proof}
If every entry of this $6\times6$ matrix has error below
$\varepsilon=10^{-12}$, its induced column-norm error is below
$6\varepsilon$. The program accordingly enlarges the computed
defect norm by $\|W_\eta\|_1\,6\varepsilon$, then verifies that the
result is below $.003$ and that the resulting inverse bound is below
$32$. This is the same perturbation-of-identity argument as in
Section~\ref{cert:matrix-bounds}. It checks both small inverse
certificates without asserting that the full $N=256$ rational
prescription has been executed.
==== END A-sharp-Fourier-certificate-for-planar-circle-packing//rational-check.tex ====
==== BEGIN A-sharp-Fourier-certificate-for-planar-circle-packing//signs.tex ====
\section{Global signs}\label{sec:signs}
We now prove the global signs for the exact functions supplied by
Proposition~\ref{interpolation:solution}. The strict bounds will also
identify their zeros within the respective sign domains. Write
\[
\sigma_1=-1,\qquad \sigma_2=1,\qquad
\delta=3\cdot10^{-5}.
\]
The difference between each exact coefficient list and its tabulated list has
norm less than $\delta$. Their fixed coefficients at zero and their constants
$C_i$ agree. Tildes will always denote the functions constructed from the
finite tabulated lists, so that
$\widetilde H_i=\widetilde p_i+K_{\widetilde p_{3-i}}$.
We use the spectral envelopes of Section~\ref{sec:interpolation} throughout.
\subsection{Taylor quotients on half-gaps}
Set $\mathcal N_0=\{0\}\cup\mathcal N$. For each
$m\in\mathcal N_0\cap[0,40]$, let $m^+$ be its successor in
$\mathcal N_0$ and include
\[
y=\frac{m^+-m}{2}.
\]
If $m>0$, also include $y=(m^--m)/2$, where $m^-$ is its predecessor.
Thus $y$ is signed, and the associated half-gap is
\[
I_{m,y}=\{m+yv:0\le v\le1\}.
\]
For function index $2$ we use every such interval, with order
$\nu=0$ at $m=0$ and $\nu=2$ otherwise. For function index $1$ we use only
those intervals contained in $[1,\infty)$, with $\nu=1$ at $m=1$ and
$\nu=2$ otherwise. In particular, the only odd order occurs on the
rightward half-gap from $1$.
There are $22$ centers from $0$ through $40$. They give $43$ half-gaps for
index $2$ and, after removing the two intervals not contained in $[1,\infty)$, $41$ for index
$1$, for a total of $84$ certified intervals.
Fix one of the specified triples $(i,m,y)$. For $s=m+u$ in its interval,
we will prove positivity of
\[
\frac{\sigma_i H_i(m+u)}{u^\nu}.
\]
Here $u^0=1$, and for $\nu>0$ the quotient at $u=0$ is interpreted by
continuity. The interpolation conditions make the derivatives of $H_i$
of orders below $\nu$ vanish at $m$, so this extension exists. Positivity
of the quotient gives the required sign: every order is even except for
the rightward half-gap from $1$, where $u\ge0$.
Figure~\ref{fig:half-gaps} shows the node-to-midpoint intervals on which
these quotients are studied.
\input{figures/half-gaps}
The tabulated function $\widetilde H_i$ approximates $H_i$ but need not
reproduce its exact zero and slope at a node. We therefore compare the
exact quotient with the Taylor series of $\widetilde H_i$ after deleting
the terms of orders below $\nu$. Put $\Delta_i=H_i-\widetilde H_i$.
Since $\widetilde H_i$ is entire and the lower derivatives of $H_i$ vanish,
this Taylor series gives
\begin{align}\label{signs:deleted-jets}
&\frac{\sigma_i H_i(m+u)}{u^\nu}
-\sigma_i\sum_{r=0}^{\infty}
\frac{\widetilde H_i^{(r+\nu)}(m)}{(r+\nu)!}u^r
\nonumber\\
&\hspace{1cm}=
\frac{\sigma_i}{u^\nu}
\left(\Delta_i(m+u)-
\sum_{j=0}^{\nu-1}\frac{\Delta_i^{(j)}(m)}{j!}u^j\right).
\end{align}
For $\nu=0$, the sum on the second line is empty and the right-hand side is
simply $\sigma_i\Delta_i(m+u)$. For $\nu>0$, Taylor's integral remainder
gives the exact identity
\[
\frac{\Delta_i(m+u)-\sum_{j=0}^{\nu-1}
\Delta_i^{(j)}(m)u^j/j!}{u^\nu}
=\frac{1}{(\nu-1)!}\int_0^1(1-t)^{\nu-1}
\Delta_i^{(\nu)}(m+tu)\,\dd t.
\]
It holds for either sign of $u$ and extends continuously to $u=0$.
Since $m+tu$ stays in the same half-gap, the absolute value of the
right-hand side of \eqref{signs:deleted-jets} is at most
\begin{equation}\label{signs:coefficient-error}
\frac{1}{\nu!}\sup_{s\in I_{m,y}}
\abs{\Delta_i^{(\nu)}(s)}.
\end{equation}
The same bound applies at $u=0$ by continuity and for either orientation
of the half-gap. Deleting the lower Taylor terms is essential here:
dividing the raw table error by $u^\nu$ could be singular at the node.
After writing $u=yv$, the first $29$ terms of the tabulated series in
\eqref{signs:deleted-jets} form a polynomial of degree at most $28$ in
$v$. We bound it on $[0,1]$ using the standard range enclosure by Bernstein
coefficients \cite[Section~3]{TitiGarloff2019}, retaining its elementary
proof below.
\begin{lemma}[Finite Bernstein certificate]\label{signs:bernstein}
For every specified triple $(i,m,y)$, put
\begin{equation}\label{signs:power-coefficients}
T_r=\sigma_i y^r
\frac{\widetilde H_i^{(r+\nu)}(m)}{(r+\nu)!},
\qquad 0\le r\le28.
\end{equation}
Then, for $0\le k\le28$,
\begin{equation}\label{signs:bernstein-bound}
B_k:=\sum_{r=0}^k
\frac{\binom{k}{r}}{\binom{28}{r}}T_r
>\begin{cases}
0.74,&m=0,\\
0.18,&m=1,\\
0.02,&m=3,4,\\
0.009,&4<m\le40.
\end{cases}
\end{equation}
\end{lemma}
\begin{proof}
The primary interval calculation in Appendix~\ref{sec:certificates} encloses
the exact derivatives of the tabulated functions in
\eqref{signs:power-coefficients}, using the quadrature assertion of
Lemma~\ref{cert:moment-error} for the analytic moment remainder.
Its exact grouped bounds~\eqref{cert:bernstein-table} imply the stated
thresholds.
\end{proof}
The relevance of this finite statement is the identity
\begin{equation}\label{signs:bernstein-identity}
\sum_{r=0}^{28}T_r v^r
=\sum_{k=0}^{28} B_k\binom{28}{k}v^k(1-v)^{28-k}.
\end{equation}
Indeed, for $0\le r\le k\le d=28$, the binomial identity
\[
\frac{\binom{k}{r}}{\binom{d}{r}}\binom{d}{k}
=\binom{d-r}{k-r}
\]
gives
\[
\sum_{k=r}^{d}\frac{\binom{k}{r}}{\binom{d}{r}}
\binom{d}{k}v^k(1-v)^{d-k}
=v^r\sum_{l=0}^{d-r}\binom{d-r}{l}v^l(1-v)^{d-r-l}=v^r.
\]
Multiply by $T_r$, sum over $r$, and interchange the two finite sums to
obtain \eqref{signs:bernstein-identity}.
For $0\le v\le1$, the basis functions on the right are nonnegative and
sum to one. Hence the lower bounds in
\eqref{signs:bernstein-bound} hold on the whole interval for the polynomial
on the left. To obtain a lower bound for the exact quotient, it remains
to control the coefficient error in \eqref{signs:coefficient-error} and the
terms beyond degree $28$ in the tabulated series.
\subsection{Bounds for the two errors}
The density bounds and \eqref{interpolation:moment-bound}, together with
the separate coefficient error bounds $\delta$ for the two functions,
give, for $s\ge s_*\ge0$,
\begin{equation}\label{signs:perturbation-envelope}
\abs{\Delta_i^{(\nu)}(s)}
\le \delta\sum_{j=1}^6 M_j
\left(\frac{(\pi t_j)^\nu}{6}
+L_j(s_*)(\pi a_{j-1}^*)^\nu\right).
\end{equation}
There are no atom contributions because the constants $C_i$ have not
changed. Substitution of the envelopes proves the following uniform
bounds for \eqref{signs:coefficient-error}:
\begin{equation}\label{signs:perturbation-caps}
\frac{1}{\nu!}\sup_{I_{m,y}}\abs{\Delta_i^{(\nu)}}
<\begin{cases}
0.004,&m=0,\\
0.011,&m=1,\\
0.002,&m\ge3.
\end{cases}
\end{equation}
Here one uses $s_*=0$ for $m=0$, $s_*=1/2$ for $m=1$ and both possible
orders, and $s_*=2$ for $m\ge3$. For example the first bound is
$\delta\sum_{j=1}^6(M_j/6+M_jL_j(0))<0.004$.
All these are bounds on intervals, obtained from the derivative envelopes.
It remains to truncate the infinite series in
\eqref{signs:deleted-jets}. Put $Y=\abs{y}$ and $q=29+\nu$.
The total variation of the direct spectral measure of either tabulated
function is less than $80$, since
\[
2.29\sum_{j=1}^6\frac{M_j}{6}
+0.017\sum_{j=0}^6 m_j^*<80.
\]
For the transformed part, use coefficient norm less than $2.3$ and
$\abs{C_i}\le0.017$. Expanding each spectral exponential at $m$ gives
the following bound for the omitted quotient terms:
\begin{align}\label{signs:taylor-bound}
\mathcal E(m,Y,\nu)
:= {}&\frac{Y^{-\nu}}{q!}
\frac{1}{1-\pi a_0^*Y/(q+1)}
\Bigg(80(\pi Y)^q
\nonumber\\[-2pt]
&\qquad+2.3\sum_{j=1}^6
M_jL_j(m)(\pi a_{j-1}^*Y)^q
+0.017\sum_{j=0}^6
L_j^P(m)(\pi a_j^*Y)^q\Bigg).
\end{align}
Indeed the first omitted exponential term has degree $q$ and subsequent
absolute-value ratios are at most $\pi a_0^*Y/(q+1)<1$.
After division by $u^\nu$, all powers still have degree at least $29$,
so their maximum for $\abs{u}\le Y$ is bounded by evaluating these
absolute-value majorants at $Y$. This proves
\eqref{signs:taylor-bound} also when $y<0$.
The five representative substitutions are
\begin{equation}\label{signs:taylor-representatives}
(m,Y,\nu)=(0,\tfrac12,0),\ (1,1,1),\ (1,1,2),\ (3,1,2),\
(4,\tfrac32,2).
\end{equation}
In each case \eqref{signs:taylor-bound} is less than $0.001$.
Increasing $m$ or decreasing $Y$ only decreases the bound. The first
three cases handle centers $0$ and $1$, the fourth handles $m=3$, and the
last handles every center $m\ge4$. Thus the truncation error is less
than $0.001$ on every required half-gap.
Combining \eqref{signs:bernstein-identity},
\eqref{signs:deleted-jets}, \eqref{signs:perturbation-caps}, and
\eqref{signs:taylor-bound}, we obtain
\begin{equation}\label{signs:finite-quotient}
\frac{\sigma_i H_i(s)}{(s-m)^\nu}
>\begin{cases}
0.735,&m=0,\\
0.168,&m=1,\\
0.017,&m=3,4,\\
0.006,&4<m\le40,
\end{cases}
\qquad s\in I_{m,y},
\end{equation}
with endpoint values interpreted by continuity. Since the only odd
order has $s-m\ge0$, these inequalities imply the desired signs on all
of these intervals. Their union is $[0,41.5]$ for index $2$ and
$[1,41.5]$ for index $1$: adjacent half-gaps meet at their common
midpoint, and the last interval runs from $40$ halfway to $43$.
\subsection{The unbounded region}
For $s\ge s_0:=41.5$, split the exact direct part as
\begin{equation}\label{signs:tail-split}
\begin{aligned}
p_i(s)&=P(s)R_{i,0}(s)+p_{i,1}(s),\\
R_{i,0}(s)&=C_i+
\sum_{\substack{n\in\mathcal N_0\\n\le40}}
\left(\frac{c_{i,n}}{(s-n)^2}+\frac{d_{i,n}}{s-n}\right),\\
\widetilde R_{i,0}(s)&=C_i+
\sum_{\substack{n\in\mathcal N_0\\n\le40}}
\left(\frac{\widetilde c_{i,n}}{(s-n)^2}
+\frac{\widetilde d_{i,n}}{s-n}\right).
\end{aligned}
\end{equation}
Here $p_{i,1}$ is the cardinal contribution from exact coefficients at
nodes greater than $40$. The function $R_{i,0}$ uses the exact coefficients,
whereas
$\widetilde R_{i,0}$ uses the finite table; their constants and node-zero
coefficients agree. All denominators in these two rational functions are
at least $1.5$. Write the corresponding decomposition of the exact
function as
\[
E_i:=p_{i,1}+K_{p_{3-i}},\qquad H_i=P R_{i,0}+E_i.
\]
For each $s\ge s_0$, the final comparison will use a nearest node $m\ge43$.
At such a node, the exact interpolation conditions and the double zero of
$P$ make $E_i$ have zero value and slope; we verify this at the comparison
step below. The useful scale for both terms is therefore $(s-m)^2$. We will bound
$\sigma_iR_{i,0}$ below, bound $P(s)$ below by a multiple of $(s-m)^2$,
and bound $|E_i(s)|$ above at the same scale through a uniform
second-derivative estimate.
We first show
\begin{equation}\label{signs:rational-tail}
\sigma_i R_{i,0}(s)>0.011\qquad(s\ge s_0,\ i=1,2).
\end{equation}
For the tabulated coefficients put
$u_n^*=\sigma_i\widetilde c_{i,n}$,
$v_n^*=\sigma_i\widetilde d_{i,n}$, and $[w]_-:=\min\{w,0\}$.
Using
\[
\frac{1}{s-n}=\frac1s+\frac{n}{s(s-n)}
\]
and discarding positive terms gives the following lower bound for
$\sigma_i\widetilde R_{i,0}(s)$, uniformly for $s\ge s_0$:
\begin{equation}\label{signs:rational-lower-sum}
\sigma_i C_i+
\frac{[\sum_{n\le40}v_n^*]_-}{s_0}
+\sum_{n\le40}
\left(\frac{[u_n^*]_-}{(s_0-n)^2}
+\frac{[v_n^*]_-n}{s_0(s_0-n)}\right),
\end{equation}
where both sums include $n=0$. Each negative reciprocal term is smallest
at $s=s_0$. The exact rational sums in
\eqref{signs:rational-lower-sum} exceed $0.012$ for index $1$ and $0.016$
for index $2$. To pass to the exact rational part, use $s-n\ge1.5>1$ and
the coefficient enclosure to obtain
\[
\begin{split}
\abs{R_{i,0}(s)-\widetilde R_{i,0}(s)}
&\le\sum_{n\le40}\left(
\frac{\abs{c_{i,n}-\widetilde c_{i,n}}}{(s-n)^2}
+\frac{\abs{d_{i,n}-\widetilde d_{i,n}}}{s-n}\right)\\
&<\frac{\delta}{s_0-40}=0.00002.
\end{split}
\]
The finite lower sums therefore imply the asserted exact bound
\eqref{signs:rational-tail}.
The coefficient norm of the tabulated part with $n>40$ is less than
$0.00002$ for either function. Hence the exact coefficient norm for
$p_{i,1}$ is less than $0.00005$. The derivative envelopes imply
\begin{align}\label{signs:tail-second-derivative}
\sup_{s\ge s_0}\abs{E_i''(s)}
\le {}&0.00005\cdot77
+2.3\sum_{j=1}^6 M_jL_j(s_0)(\pi a_{j-1}^*)^2
\nonumber\\
&+0.017\sum_{j=0}^6 L_j^P(s_0)(\pi a_j^*)^2
<0.006.
\end{align}
For the direct part we used
$\sum_{j=1}^6 M_j(\pi t_j)^2/6<77$; for the transformed part we used the
all-node coefficient norm less than $2.3$ from
Proposition~\ref{interpolation:solution}.
It remains to give the quadratic lower bound for $P$.
\begin{lemma}\label{signs:P-lower-bound}
If $m\in12\ZZ+A$ is a nearest zero of $P$ to a real number $s$, then
\[
P(s)\ge0.68(s-m)^2.
\]
\end{lemma}
\begin{proof}
On either half-gap from $m$, consider $P(m+u)/u^2$, with value $Q_m$ at
$u=0$. Put $a=\pi/12$. For a sine factor whose residue is $a'\ne m$
modulo $12$, the second derivative of its logarithm is
$-2a^2\csc^2(a(m+u-a'))<0$. These factors stay nonzero on the half-gap.
The factor vanishing at $m$, after division by $u^2$, has the same
logarithmic concavity because
\[
\frac{\dd^2}{\dd u^2}\log\abs{\frac{\sin(au)}{u}}
=\frac1{u^2}-\frac{a^2}{\sin^2(au)}\le0.
\]
Here $\abs{\sin(au)}\le\abs{au}$, and the logarithm has a smooth limiting
value at $u=0$. Squaring the divided sine only multiplies this second
derivative by two.
The quotient is therefore logarithmically concave on each half-gap,
including its continuous value at the center. Its minimum is bounded
below by the smaller of its endpoint values. At the center
$Q_m>1.75$. At the gap midpoints, direct substitution gives the values
listed below; periodicity reduces the calculation to one period.
To evaluate them, put $\theta=\pi(s-2)/12$ and $c=\cos(2\theta)$.
The identity $4\sin v\sin w=2(\cos(v-w)-\cos(v+w))$ shows that the paired
unsquared factors for residues $(0,4)$, $(1,3)$, and $(7,9)$ are,
respectively,
\[
1-2c,\qquad \sqrt3-2c,\qquad \sqrt3+2c.
\]
Consequently
\[
P(s)=(4c^2-3)^2(1-2c)^2,
\qquad c=\cos\frac{\pi(s-2)}6,
\]
which evaluates every entry directly.
\begin{center}
\begin{tabular}{c|c}
Gap between consecutive zeros & Midpoint quotient $P(s)/(s-m)^2$\\
\hline
$[0,1],\ [3,4]$ & $12-8\sqrt2$\\
$[1,3]$ & $1$\\
$[7,9]$ & $9$\\
$[4,7],\ [9,12]$ & $\dfrac89(3+\sqrt2+\sqrt3+\sqrt6)$
\end{tabular}
\end{center}
The distances to the two ends of a gap agree at its midpoint, so this
table supplies the endpoint value for both half-gaps. The smallest value
is $12-8\sqrt2>0.68$: indeed $(283/200)^2>2$, so
$\sqrt2<283/200$ and $12-8(283/200)=0.68$. The other three values are
at least $1$. This proves the claim.
\end{proof}
For any $s\ge41.5$, choose a nearest node $m\ge43$; at $s=41.5$, choose
$43$ rather than the tied node $40$. Both $H_i$ and $P R_{i,0}$ vanish
to order at least two at $m$, by \eqref{fourier:interpolation} and the absence of a
pole of $R_{i,0}$ there. Thus
$E_i=H_i-P R_{i,0}$ vanishes to order at least two at $m$. The segment between
$s$ and $m$ lies in $[41.5,\infty)$, and
\eqref{signs:tail-second-derivative} gives
\[
\abs{E_i(s)}\le\frac{0.006}{2}(s-m)^2.
\]
Using \eqref{signs:rational-tail} and
Lemma~\ref{signs:P-lower-bound}, we obtain on the entire unbounded region
\begin{equation}\label{signs:tail-margin}
\sigma_iH_i(s)
\ge\left(0.011\cdot0.68-\frac{0.006}{2}\right)(s-m)^2
=0.00448(s-m)^2\ge0.
\end{equation}
\begin{proposition}[Global signs and zeros]\label{signs:global}
The exact functions have the zero sets
\begin{equation}\label{signs:zero-sets}
\{s\ge1:H_1(s)=0\}=\mathcal N,\qquad
\{s\ge0:H_2(s)=0\}=\mathcal N,
\end{equation}
and satisfy the strict inequalities
\[
H_1(s)<0\quad(s\ge1,\ s\notin\mathcal N),\qquad
H_2(s)>0\quad(s\ge0,\ s\notin\mathcal N).
\]
The zero of $H_1$ at $1$ is simple. Its zeros at $n\in\mathcal N$ with
$n>1$, and the zeros of $H_2$ at every $n\in\mathcal N$, have order exactly
two as zeros of the entire functions of $s$.
Consequently $F(x)\le0$ when $b\abs{x}^2\ge1$, and
$\widehat F(\xi)\ge0$ for every $\xi\in\RR^2$.
\end{proposition}
\begin{proof}
At every node $n\in\mathcal N$, the exact interpolation conditions give
$H_1(n)=H_2(n)=0$. Away from the nodes, the quotient bounds
\eqref{signs:finite-quotient} are strict on the finite region, and
\eqref{signs:tail-margin} is strict on the remaining region. This proves
\eqref{signs:zero-sets} and the asserted strict signs.
The interpolation slope $H_1'(1)=-Q_1\ne0$ makes its zero at $1$ simple.
For $H_2$ at every node in $\mathcal N\cap[1,40]$, and for $H_1$ at those
nodes greater than $1$, the chosen order is $\nu=2$.
The endpoint value in \eqref{signs:finite-quotient} is then
$\sigma_iH_i''(m)/2>0$, so the zero has order exactly two.
For a node $m\in\mathcal N$ with $m\ge43$, take $u$ sufficiently small that
$m$ remains the nearest node to $m+u$ and $m+u\ge41.5$. The exact value and first derivative
at $m$ vanish, so \eqref{signs:tail-margin} yields
\[
\frac{\sigma_iH_i''(m)}2
=\lim_{u\to0}\frac{\sigma_iH_i(m+u)}{u^2}
\ge0.00448>0.
\]
These zeros also have order exactly two. Finally, the factors
$e^{-\pi h s}$ in $F$ and $\widehat F$ are positive, which gives their
non-strict sign inequalities.
\end{proof}
==== END A-sharp-Fourier-certificate-for-planar-circle-packing//signs.tex ====
==== BEGIN A-sharp-Fourier-certificate-for-planar-circle-packing//references.bib ====
@article{CohnElkies2003,
author = {Cohn, Henry and Elkies, Noam},
title = {New upper bounds on sphere packings {I}},
journal = {Annals of Mathematics},
volume = {157},
number = {2},
year = {2003},
pages = {689--714},
doi = {10.4007/annals.2003.157.689},
url = {https://annals.math.princeton.edu/2003/157-2/p09}
}
@article{Hales2000,
author = {Hales, Thomas C.},
title = {Cannonballs and honeycombs},
journal = {Notices of the American Mathematical Society},
volume = {47},
number = {4},
year = {2000},
pages = {440--449},
url = {https://www.ams.org/notices/200004/fea-hales.pdf}
}
@article{Viazovska2017,
author = {Viazovska, Maryna S.},
title = {The sphere packing problem in dimension 8},
journal = {Annals of Mathematics},
volume = {185},
number = {3},
year = {2017},
pages = {991--1015},
doi = {10.4007/annals.2017.185.3.7},
url = {https://annals.math.princeton.edu/2017/185-3/p07}
}
@article{CKMRV2017,
author = {Cohn, Henry and Kumar, Abhinav and Miller, Stephen D. and
Radchenko, Danylo and Viazovska, Maryna},
title = {The sphere packing problem in dimension 24},
journal = {Annals of Mathematics},
volume = {185},
number = {3},
year = {2017},
pages = {1017--1033},
doi = {10.4007/annals.2017.185.3.8},
url = {https://annals.math.princeton.edu/2017/185-3/p08}
}
@article{Johansson2017,
author = {Johansson, Fredrik},
title = {{Arb}: Efficient arbitrary-precision midpoint-radius interval arithmetic},
journal = {IEEE Transactions on Computers},
volume = {66},
number = {8},
year = {2017},
pages = {1281--1292},
doi = {10.1109/TC.2017.2690633},
url = {https://arxiv.org/abs/1611.02831v1}
}
@article{RadchenkoViazovska2019,
author = {Radchenko, Danylo and Viazovska, Maryna},
title = {Fourier interpolation on the real line},
journal = {Publications Math{\'e}matiques de l'IH{\'E}S},
volume = {129}, pages = {51--81}, year = {2019},
doi = {10.1007/s10240-018-0101-z},
note = {The cited theorem numbering follows arXiv:1701.00265v2},
url = {https://arxiv.org/abs/1701.00265v2}
}
@article{CKMRV2022,
author = {Cohn, Henry and Kumar, Abhinav and Miller, Stephen D. and
Radchenko, Danylo and Viazovska, Maryna},
title = {Universal optimality of the {$E_8$} and {Leech} lattices and interpolation formulas},
journal = {Annals of Mathematics},
volume = {196}, number = {3}, pages = {983--1082}, year = {2022},
doi = {10.4007/annals.2022.196.3.3},
note = {The cited theorem numbering follows arXiv:1902.05438v3},
url = {https://arxiv.org/abs/1902.05438v3}
}
@misc{Sardari2021,
author = {Talebizadeh Sardari, Naser},
title = {Higher {Fourier} interpolation on the plane},
year = {2021},
eprint = {2102.08753}, archivePrefix = {arXiv}, primaryClass = {math.NT},
note = {Version 2, May 4, 2021},
url = {https://arxiv.org/abs/2102.08753v2}
}
@article{TitiGarloff2019,
author = {Titi, Jihad and Garloff, J{\"u}rgen},
title = {Matrix methods for the tensorial {Bernstein} form},
journal = {Applied Mathematics and Computation},
volume = {346}, pages = {254--271}, year = {2019},
doi = {10.1016/j.amc.2018.08.049},
url = {https://www-home.htwg-konstanz.de/~garloff/AMC-2.pdf}
}
@misc{Zhitniaia2026,
author = {Zhitniaia, Daria},
title = {Linear Programming Bound for Sphere Packing in Dimension 2},
year = {2026},
howpublished = {Bachelor's thesis, Faculty of Mathematics, HSE University},
note = {Institutional abstract; full thesis unavailable in the inspected record},
url = {https://www.hse.ru/en/edu/vkr/1163516275}
}
@article{CarneiroLittmannVaaler2013,
author = {Carneiro, Emanuel and Littmann, Friedrich and Vaaler, Jeffrey D.},
title = {Gaussian subordination for the {Beurling--Selberg} extremal problem},
journal = {Transactions of the American Mathematical Society},
volume = {365}, number = {7}, pages = {3493--3534}, year = {2013},
doi = {10.1090/S0002-9947-2013-05716-9},
url = {https://arxiv.org/abs/1008.4969v2}
}
@article{Waldvogel2006,
author = {Waldvogel, J{\"o}rg},
title = {{Fast Construction of the Fej{\'e}r and Clenshaw--Curtis Quadrature Rules}},
journal = {BIT Numerical Mathematics},
year = {2006},
volume = {46},
number = {1},
pages = {195--202},
doi = {10.1007/s10543-006-0045-4},
url = {https://people.math.ethz.ch/~waldvoge/Papers/fejer.pdf}
}
==== END A-sharp-Fourier-certificate-for-planar-circle-packing//references.bib ====