A sharp Fourier certificate for planar circle packing

1 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

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 [5, pp. 442–443].

The Fourier-analytic approach of Cohn and Elkies [2] 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

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.

Theorem 1.1 gives an affirmative answer to the two-dimensional case of Cohn and Elkies’s Conjecture 7.3 [2], 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

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

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 5.2). 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 2 gives an explicit extra node, and Section 5 records the zero sets in their respective sign domains and final radial coordinates.

1.1 Fourier certificates and interpolation

Cohn and Elkies formulated their sharpness conjecture in dimensions 2, 8, and 24 [2, Conjecture 7.3]. 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 [2, Section 7, pp. 704–705]. Exact attainment of the conjectured sharp bound remained unresolved in that approach. Viazovska constructed the required function in dimension 8 [10, Theorems 3 and 4], and Cohn, Kumar, Miller, Radchenko, and Viazovska constructed its dimension-24 counterpart [3, Section 4]. 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 [2, p. 694].

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 [7, Theorem 1]. 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 [4, Theorems 1.4, 1.7, and 1.9].

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 [4, Conjecture 7.5]. Sardari proved this nonuniqueness by constructing infinitely many linearly independent radial Schwartz functions with double zeros, before and after Fourier transformation, on those shells [8, Theorem 1.7]. 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 [12] 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 1.1 remains to be determined from the full thesis.

1.2 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 [2, Section 5]. 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 1.1. 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 [8, Section 6.1]; 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 [1, formulas (1.2)–(1.3)]; 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 2 develops the cardinal functions and their Fourier transforms. Section 3 solves the interpolation equations, and Section 4 proves the global signs and identifies their zeros. Section 5 establishes the exact normalization, performs the final dilation, and derives the packing consequences. Appendix A gives the finite data and the primary interval certificate, including the analytic integration error that connects its finite sums to exact moments. Appendix B gathers the rational alternatives and independent checks, with their precise mathematical and computational scopes.

2 Gaussian cardinal functions

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

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

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

We impose interpolation conditions at the positive nodes of a periodic set:

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 [8, Section 6.1]; the Gaussian transformation and the coupled inversion needed here will be proved directly.

2.1 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 [1, formulas (1.2)–(1.3)]. The nodes and simultaneous Fourier constraints here are different.

For -6\le j\le6, put t_j=j/6, and define

Here P_{-j}=\overline{P_j}, and

For example, with w=e^{i\pi s/6} and \rho=e^{i\pi/6}, each factor is

The six minus signs cancel, and \rho^{-\sum_{a\in A}a}=\rho^{-24}=1. Thus the product in (2.5) equals

Expanding this Laurent polynomial gives (2.6).

Every zero n\in12\ZZ+A of P is double. Define Q_n,D_n by

If a\in A is the residue of n modulo 12, then

Thus Q_n,D_n are periodic in n, and their exact values are

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,

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

Away from the nodes, define

We next show that the apparent singularities are removable, with uniform control of the resulting entire function.

In particular, for every integer k\ge0,

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 (2.7) gives the cardinal identities

The other summands in (2.11) vanish to second order at n, which explains why these identities involve only its own pair of coefficients.

2.2 Fourier transformation

For a cardinal function p, let

Define, for real -1\le t\le1,

and, using the measure from Lemma 2.1, put

No cardinal representation of K_p is assumed. Only its entireness, its values and derivatives, and the estimates obtained from (2.19) 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.

2.3 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

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

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 (2.23) are real, radial, and Schwartz for every pair of admissible coefficient lists.

The required interpolation conditions are

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.

3 Solving the interpolation equations

The equations (2.25) 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

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 (2.11) with these positive-node coefficients and with C=c_0=d_0=0. Define the real linear operator S by

The bounds proved below show, in particular, that S is a bounded operator on X.

We use two successive decompositions of this operator. Set

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.

3.1 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.

At each endpoint t_j=j/6,

On the interval t_{j-1}\le t\le t_j, the corresponding bounds are

The density bounds hold uniformly for n\in\{0\}\cup\mathcal N.

Here are direct ways to verify these estimates. The identities

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 (3.3) 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 jth interval is

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

For s\ge0, set

When \beta_js=0, the minimum in the first line is understood to be 1/6. Integrating the exponential envelope in (3.3) 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 (2.11),

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,

because the minimum factor in (3.5) is nonincreasing. For every real cutoff x\ge0, the value and derivative terms in (3.1), together with (3.6), therefore imply

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 2.1; 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 (3.1), not merely the coordinates of an auxiliary matrix limit.

For x=0,4,16,100, respectively, the right side of (3.7) is less than

We shall use the slightly weaker bounds

For the input C=1, the analogue of (3.7) 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

The scalar substitutions in these envelope estimates use only the explicit constants and the node data. Section A.6 records the finite scalar bounds and the precise coverage of their supplied verifier.

3.2 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 2, with the indicated factor 10^{-10}, and are zero at positive nodes beyond 100. Append the fixed data (2.22) to obtain \widetilde p_i, and write \widetilde K_i=K_{\widetilde p_i}.

3.3 Inversion on the whole coefficient space

We first prove an inverse bound that keeps finite and tail residuals separate.

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 [8, Theorem 1.7].

4 Global signs

We now prove the global signs for the exact functions supplied by Proposition 3.4. The strict bounds will also identify their zeros within the respective sign domains. Write

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 3 throughout.

4.1 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

If m>0, also include y=(m^--m)/2, where m^- is its predecessor. Thus y is signed, and the associated half-gap is

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

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 1 shows the node-to-midpoint intervals on which these quotients are studied.

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

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

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 () is at most

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 () form a polynomial of degree at most 28 in v. We bound it on [0,1] using the standard range enclosure by Bernstein coefficients [9, Section 3], retaining its elementary proof below.

The relevance of this finite statement is the identity

Indeed, for 0\le r\le k\le d=28, the binomial identity

gives

Multiply by T_r, sum over r, and interchange the two finite sums to obtain (4.5). For 0\le v\le1, the basis functions on the right are nonnegative and sum to one. Hence the lower bounds in (4.4) 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 (4.2) and the terms beyond degree 28 in the tabulated series.

4.2 Bounds for the two errors

The density bounds and (3.6), together with the separate coefficient error bounds \delta for the two functions, give, for s\ge s_*\ge0,

There are no atom contributions because the constants C_i have not changed. Substitution of the envelopes proves the following uniform bounds for (4.2):

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 (). 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

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:

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 () also when y<0.

The five representative substitutions are

In each case () 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 (4.5), (), (4.7), and (), we obtain

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.

4.3 The unbounded region

For s\ge s_0:=41.5, split the exact direct part as

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

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

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

and discarding positive terms gives the following lower bound for \sigma_i\widetilde R_{i,0}(s), uniformly for s\ge s_0:

where both sums include n=0. Each negative reciprocal term is smallest at s=s_0. The exact rational sums in (4.13) 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

The finite lower sums therefore imply the asserted exact bound (4.12).

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

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 3.4.

It remains to give the quadratic lower bound for P.

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 (2.24) 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 () gives

Using (4.12) and Lemma 4.2, we obtain on the entire unbounded region

5 Poisson 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 3.4, and put

Section 2 proves that F is real, radial, and Schwartz, with \widehat F(\xi)=G_2(b|\xi|^2). The interpolation conditions give

Indeed, G_i'(n)=e^{-\pi hn}H_i'(n) because H_i(n)=0.

5.1 The dilation identity

Use the lattice \Lambda of covolume one and the shell coordinate s_a from (2.2)–(2.3). 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

As shown in Section 2, \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,

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 (5.2).

5.2 The final dilation

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

Proposition 4.3 therefore gives the exact zero sets in the respective sign domains:

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

5.3 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

Poisson summation on each translated lattice gives

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 [2]. 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 [2, Section 8].

5.4 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.

A Finite certificates

This appendix completes the finite evidence used in Lemma 3.2 and Lemma 4.1. 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 B gives the rational alternatives and their separately stated scope.

A.1 Exact input lists and matrices

Table 2 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,

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

Here and below an index n\ge0 in a coefficient sum ranges over \{0\}\cup\mathcal N.

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

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 A.3.

A.2 Moments and the integration rule

For a finite input p=PR supported on nodes at most 100, put

The domain for the quadrature error estimate is stated explicitly in Lemma A.1.

We use the first Fejér quadrature rule, whose cosine-sum weights are described in [11, Section 2, formula (2.3)]. 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

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 (A.6), so the rule is exact through degree 255. Moreover,

Thus the six pieces together have total weight one.

For clarity, the complete density on piece j used in this appendix is

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

The real part is taken after summation.

This domain contains all basis inputs and the two tabulated inputs used below.

A.3 Interval evaluation

For a finite input p in the real finite-input domain of Lemma A.1 and a node row m\in I_f, form the two real numbers

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 4.1; 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

The signed scaling is exact. Definition (A.5) gives U_{p,l}(m,y)=y^lU_{p,l}(m,1) and the same identity for X. Therefore, for positive or negative y,

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 (A.6), (A.7), (A.8), (A.11), and (A.12) 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.

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

Here stack means all first entries followed by all second entries. The recurrence for a equals \binom{k}{r}/\binom{28}{r}. 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 (A.12). Each unscaled real moment enclosure is enlarged by [-10^{-22},10^{-22}], as justified by Lemma A.1. Consequently the analytic-remainder contribution to the enclosure radius of a Bernstein coefficient is bounded by

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 [6]. 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 verification/verify_all.py runs two verifiers in the same directory, interval_certificate.py and interval_bounds.py. Their output records, finite.json and scalar.json, and the aggregate record summary.json are placed in its results/ subdirectory. The first verifier checks the exact bounds in (A.13), (A.14), (A.15), and (A.16), as well as all 2436 individual Bernstein inequalities in Lemma 4.1. 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.

A.4 Matrix and residual certificates

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 (A.11).

The unmultiplied exterior block also satisfies

The certified enclosure calculation described in Section A.3 proves these inequalities. For example, the two defect norms are enclosed respectively in

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 (A.13) and (A.14) 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

It follows from (A.13) that

Equation (A.14) gives the remaining bound <.09 in Lemma 3.2.

For the two full tabulated inputs, including (A.1), the residual certificate is

The primary evaluation in Section A.3 proves (A.15) for these exact residuals. Since 3\cdot10^{-10}<10^{-9}, these inequalities, the preceding primary matrix conclusion, and (A.2) prove every finite assertion of Lemma 3.2.

A.5 Bernstein coefficient certificates

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

The primary evaluation in Section A.3 proves the grouped bounds (A.16) for the exact derivatives of the tabulated functions. Each printed lower bound is stronger than its corresponding threshold in (4.4). Therefore these exact grouped bounds imply Lemma 4.1.

A.6 Scalar envelope and sign bounds

The remaining finite substitutions involve elementary real expressions. For the endpoint and density envelopes in (3.3), 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 3.

For the row envelope in (3.7), certified outward bounds are

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 3.

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.

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 (4.9), the corresponding bounds for \mathcal E(m,Y,\nu) in () are

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 3. 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 4.2 supply the global lower bound for P.

B Alternative certification methods

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.

B.1 Rational exponential primitives

The fixed rational procedures below use rational arithmetic alone, with a complex number represented by a pair of rationals. Put

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

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:

Real and imaginary parts of e_*(u) give the cosine and sine replacements when u is real.

Positive exponential upper bounds. For x\ge0 and an integer n\ge0 with x<n+2, define

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.

B.2 A fixed rational evaluation

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 A.1. After defining the rounded sums, we prove their uniform approximation error. No matrix, residual, or Bernstein values are asserted by that error theorem.

In (A.6), 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 (A.7) 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

and set

The rational moment sums are therefore

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 A.3, 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.

B.3 Propagation to finite certificates

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 (A.11), 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 (A.11) gives

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 (2.9) exceeds 7/4 by 1268/450000, so \eta_Q\le10^{-25} also ensures that the approximate Q_m exceeds q.

Analytic pair-error allowance. For an input p in the rational finite-input domain of Lemma B.2 and a row m\in I_f, define the rational version \mathcal S_m^\# by using the quantities with superscript \# throughout (A.11). We record the division estimate explicitly for both basis and full-table inputs. By (3.3) and the atom masses there, the total variation of the folded measure for any input in the real finite-input domain of Lemma A.1 is at most

On the real segment, |\lambda|,|z|<3 and |e^{i\pi z m}|\le1 for m\ge0. Thus, using \pi<22/7,

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 B.2, the errors in X_0,X_1 are less than \varepsilon and those in Q,D are less than \delta. Equation (2.9) 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^{\#},

The two quotient differences are consequently bounded by

In particular, the uniform error in each entry of \mathcal S_m^{\#} is strictly less than 10^{-12}; we retain this looser allowance below.

B.3.1 Bernstein multipliers

For every signed half-gap in Lemma 4.1, y\ne0 and |y|^{-\nu}\le4. Directly from (A.5),

This identity also holds for negative y and gives exactly the power coefficients of (A.12); 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

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}.

B.4 An independent rational check of the first block

The optional program verification/check_sjj_rational.py checks the 6\times6 block on J=\{1,3,4\} by exact rational arithmetic. It uses the cosine rule (A.6) 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.

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 A.4. It checks both small inverse certificates without asserting that the full N=256 rational prescription has been executed.

B.5 Comparison of first-block majorants

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:

The separated route uses the explicit density bound D_0 from Section B.4. The two common routes use density and exponential bounds (1800,2.8) and (2500,3), respectively.

B.5.1 A common kernel majorant

The rational comparisons

follow by clearing positive denominators in (B.2); 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

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,

The factor two here compares the integral and quadrature averages; a direct Taylor average has only one tail, as in Lemma B.4.

The rounding proof in Section B.4 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 (B.6), the additional Q,D contribution is less than

Consequently

B.5.2 A separate rounding majorant

The weaker bounds e^{4\pi/6}<9 and e^{4\pi(.081)}<3 also follow from (B.10). 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

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 (B.6), now with d=4.3, bounds their additional contribution by 6\cdot10^{-31}. Therefore

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.

B.6 Rational-center Taylor integration

For a finite input in the real moment domain of Lemma A.1, let F^U_{j,\ell} and F^X_{j,\ell} be the density (A.7) multiplied by the respective factors in (A.8). 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

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.

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:

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 (B.2). 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

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 (B.2) 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,

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

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

The second line follows by comparing coefficients in E'=g'E. Use (B.19) for E_0 and () 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

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

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

The constant exponentials in L_{n,j},H_{n,j},\phi_{n,0} are enclosed by (B.19). For a series v(x) equal to c_j+x or z(c_j+x), the divided powers obey

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 (A.8) 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

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 (2.8) 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 (B.20)–(B.24). Their radius formulas tend to zero when their input radii tend to zero, while (B.19) 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

on each piece and moment imply that the weighted coefficient-radius sum in (B.16) is at most \delta.

We now combine these integration budgets with the pair estimate (B.6) and the Bernstein estimate from Section B.3.1.

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 (A.15) are exact rationals, so its error is just the corresponding pair error. Finally, for the signed half-gaps of Lemma 4.1, 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 B.3.1. 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 A.

B.7 Direct node-Taylor representation

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 (2.11). Write \mathcal N_0=\{0\}\cup\mathcal N and, for m\in\mathcal N_0, put

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.

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 (2.19), and this local identity neither places K_p in the cardinal class nor replaces the exact-derivative certificates for the tabulated functions.