Certified Gaussian gridding for type 1 and type 2

Fix a dimension $d\ge 1$, finite mode bounds $K_r\ge 0$, a finite type-1 source list $(x_j,c_j{)}_{j<N}$ or type-2 data $(x_j{)}_{j<N},(f_k{)}_{k\in \mathcal{K}}$ as in Type-1, type-2, and type-3 nonuniform Fourier transforms, and a nonempty finite ordered list of candidate tuples $\mathcal{C}=(({\tau }_r,M_r,w_r,L_r{)}_{r=1}^d{)}_{\nu =1}^C,$ an interval-arithmetic format $(p,e_{\min },e_{\max })$ with integers $p\ge 2$ and $e_{\min }\le e_{\max }$, and one operationwise arithmetic model from Operationwise real and complex floating point arithmetic with rounded square roots and fused multiply-add with the same binary precision $p$ and normal exponent limits $j_{\min }=e_{\min }$ and $j_{\max }=e_{\max }$. Let $u$ be a valid unit-roundoff bound for this selected model instance; it is not a separate free control, and require $u\le 1\mathrm{/}16$. Use this same model instance for every computed FFT operation. Let $J\ge 1$ be the finite elementary-series truncation index, $0\le {ϵ}_{\mathrm{t}\mathrm{w}}\le 1\mathrm{/}4$ an exact rational maximum allowed twiddle certificate, $0<{\delta }_{\varphi }\le 1\mathrm{/}4$ an exact rational phase-evaluation tolerance, integers $2\le p_{\varphi ,0}\le p_{\varphi ,\max }$ giving the initial and maximum phase precisions, and a finite recursion-depth cap $d_{\mathrm{F}\mathrm{F}\mathrm{T}}\ge 0$. The phase evaluator uses exponent limits $e_{\min },e_{\max }$ and truncation index $J$. Also supply an exact positive dyadic requested absolute tolerance ${\epsilon }_{\mathrm{r}\mathrm{e}\mathrm{q}}>0$. Source coordinates and coefficients are finite dyadics interpreted as exact real and complex inputs. Each candidate must have rational ${\tau }_r>0$, a power-of-two grid length $M_r=2^{q_r}$ with integer $q_r\ge 0$, $M_r>2K_r$, integer width $w_r\ge 1$ with $M_r\ge 2w_r+1$, and integer image cutoff $L_r\ge 1$. A malformed global control or structurally invalid candidate returns invalid-input before arithmetic. Type 3 is not an input mode for this algorithm.

After validating the global controls and candidate shapes, if $N=0$ return status certified immediately: type 1 returns the exact zero output array and zero radii, and type 2 returns the empty output array and zero radius. For $N>0$, process the candidates in supplied order; the first candidate that produces a complete interval certificate meeting ${\epsilon }_{\mathrm{r}\mathrm{e}\mathrm{q}}$ is returned.

1. For each coordinate compute outward enclosures for $a_r=\frac{\pi }{{\tau }_r{M}_r^2},s_r=\exp \text{\hspace{0.17em}⁣}(-\pi {\tau }_r{M}_r(M_r-2K_r)),d_r=\sqrt{{\tau }_r}\exp (-\pi {\tau }_r{K}_r^2),$ $S_r=\frac{2}{1-\exp (-a_r)},T_r=2\exp (-a_r{w}_r^2)(1+\frac{1}{2a_r{w}_r}),R_r=\frac{2\exp (-\pi L_r^2\mathrm{/}{\tau }_r)}{1-\exp (-\pi L_r\mathrm{/}{\tau }_r)}\mathrm{.}$ Use Certified range evaluation for elementary real and complex functions for every real exponential and square-root node in the procedure, including the finite norms, normalization and complex moduli; use format $(p,e_{\min },e_{\max })$ and truncation index $N=J$ on every such call. Form every Machin enclosure of $\pi$ with this same format and index. Use the outward interval operations for the rational DAGs. Set $A={\prod }_{r=1}^d(1+\frac{2s_r}{1-s_r})-1,T_d={\sum }_{r=1}^d{T}_r{\prod }_{s\mathrm{\ne }r}{S}_s,d_{\min }={\prod }_{r=1}^d{d}_r\mathrm{.}$ Store outward upper bounds $A^{+},T_d^{+}$ and $R_r^{+}$, where $R_r^{+}$ bounds the exact image-tail constant $R_r$, and a lower bound $d_{\min }^{}$. If any required enclosure is non-finite, if a required denominator interval in $A,S_r,T_r,R_r$ contains zero, if some $s_r$ enclosure is not strictly below $1$, or if $d_{\min }^{}\le 0$, reject this candidate as inconclusive; no division by an unverified multiplier is performed.

2. For every source define and construct its periodic nearest-point set by the following finite rule. For each coordinate put $a_{j,r}=M_r{x}_{j,r}$, enumerate the integers from $\lceil a_{j,r}-w_r\rceil$ through $\lfloor a_{j,r}+w_r\rfloor$, and take their least nonnegative residues modulo $M_r$ as $W_{j,r}$; then set $W_j={\prod }_r{W}_{j,r}$. This equals the displayed set because its membership condition is equivalent to some integer $m\equiv {\mathrm{\ell }}_r(\mathrm{m}\mathrm{o}\mathrm{d}{M}_r)$ satisfying $\mathrm{\mid }m-a_{j,r}\mathrm{\mid }\le w_r$. Compute the endpoints by exact rational comparisons. At most $2w_r+1$ consecutive integers are enumerated, and $M_r\ge 2w_r+1$ makes their residues distinct. $W_j={\prod }_{r=1}^d\{{\mathrm{\ell }}_r\in \{0,\dots ,M_r-1\}:{\min }_{n\in \mathbb{Z}}\mathrm{\mid }{\mathrm{\ell }}_r\mathrm{/}{M}_r-x_{j,r}-n\mathrm{\mid }\le w_r\mathrm{/}{M}_r\},$ using the selected candidate. For type 1 only, for each $\mathrm{\ell }\in W_j$ reduce each exact dyadic displacement ${\mathrm{\ell }}_r\mathrm{/}{M}_r-x_{j,r}$ modulo one to $z_r\in [0,1)$. Enclose the one-dimensional kernel by evaluating the finite image sum ${\sum }_{n=-L_r}^{L_r}\exp \text{\hspace{0.17em}⁣}(-\pi (z_r+n{)}^2\mathrm{/}{\tau }_r)$ with Certified range evaluation for elementary real and complex functions using $N=J$, then add the symmetric interval $[-R_r^{+},R_r^{+}]$ from step 1. All kernel and coefficient products and sums use outward interval arithmetic. For each grid cell $\mathrm{\ell }$, this gives a complex rectangle enclosing its exact truncated-spread value $g_{\mathrm{\ell }}$. Choose its lower-left corner $u_{\mathrm{\ell }}$ as a dyadic center and an upward dyadic radius $e_{\mathrm{\ell }}$ large enough to enclose $g_{\mathrm{\ell }}-u_{\mathrm{\ell }}$. Let $E_{\mathrm{e}\mathrm{v}\mathrm{a}\mathrm{l}}$ be an outward upper bound for $({\sum }_{\mathrm{\ell }}{e}_{\mathrm{\ell }}^2{)}^{1\mathrm{/}2}$ and let $U$ be an outward upper bound for $\mathrm{\parallel }u{\mathrm{\parallel }}_2=({\sum }_{\mathrm{\ell }}\mathrm{\mid }{u}_{\mathrm{\ell }}{\mathrm{\mid }}^2{)}^{1\mathrm{/}2}$. Before a type-1 FFT call, require every real and imaginary coordinate of every $u_{\mathrm{\ell }}$ to belong to the representable set $\mathbb{F}$ of the selected operationwise model; otherwise reject this candidate with a range-failure status before calling the FFT. Type 2 skips this spreading-grid construction; its source-window intervals are evaluated in step 4 using the sets $W_j$ just defined.

3. For each mode $k\in \mathcal{K}$, define its grid slot by $\kappa (k)=(k_r<mrow><mi\; mathvariant="normal">m</mi><mi\; mathvariant="normal">o</mi><mi\; mathvariant="normal">d</mi></mrow>M_r{)}_{r=1}^d$, using least nonnegative residues; $M_r>2K_r$ makes this map injective. Put $\eta =10u+{ϵ}_{\mathrm{t}\mathrm{w}}$, ${\rho }_r=(1+\eta {)}^{q_r}-1$ for each axis, and $\rho ={\prod }_r(1+{\rho }_r)-1$. For type 1, run Computed unitary radix-two FFT with certified twiddle phases on each length-$M_r$ grid line along axes $1,\dots ,d$ in a fixed lexicographic line order, with direction forward, the selected operationwise model, phase tolerance ${\delta }_{\varphi }$, exponent limits $e_{\min },e_{\max }$, phase precisions $p_{\varphi ,0},p_{\varphi ,\max }$, series cap $J$, and recursion cap $d_{\mathrm{F}\mathrm{F}\mathrm{T}}$. The supplier passes the exact rational phase argument $q=-2m\mathrm{/}{M}_r$ for its twiddle indexed by $m$ and returns a certificate ${ϵ}_{\mathrm{\ell }}$ for each twiddle it uses. Reject this candidate with the returned status if an FFT call is non-finite or inconclusive; otherwise verify for every returned twiddle certificate that ${ϵ}_{\mathrm{\ell }}\le {ϵ}_{\mathrm{t}\mathrm{w}}$. If any certificate exceeds ${ϵ}_{\mathrm{t}\mathrm{w}}$, reject this candidate with status twiddle-bound-not-met and do not use that FFT output. This check applies literally when ${ϵ}_{\mathrm{t}\mathrm{w}}=0$, so every used twiddle must then have certificate zero. After all checks pass, extract the unitary Fourier coordinates at slots $\kappa (k)$ for $k\in \mathcal{K}$ and scale by $P^{-1\mathrm{/}2}$, where $P={\prod }_r{M}_r$. Use the range evaluator with format $(p,e_{\min },e_{\max })$ and $N=J$ for every real elementary node, including this scaling and the mode multipliers; use the same $J$ Machin enclosure of $\pi$. For each required mode evaluate the multiplier interval $d_k={\prod }_r\sqrt{{\tau }_r}\exp (-\pi {\tau }_r{k}_r^2)$; if its lower endpoint is $\le 0$, reject this candidate as inconclusive before division. Otherwise divide by that certified positive interval. Use outward complex-rectangle arithmetic. For each $k\in \mathcal{K}$, choose the lower-left corner $z_k$ of this computed scaling/deconvolution rectangle and an upward dyadic radius $b_k$ containing that rectangle about $z_k$. This rectangle alone concerns the computed FFT output, not the exact unitary transform. Separately form a certified outward upper bound $R_{\mathrm{F}\mathrm{F}\mathrm{T}}$ for $\frac{\rho U}{\sqrt{P}{d}_{\min }^{}}\mathrm{.}$ Here $U$ is the certified norm bound of the actual dyadic input grid from step 2, and $\rho$ uses the same selected model's valid bound $u$ and checked twiddle limit ${ϵ}_{\mathrm{t}\mathrm{w}}$ as the successful FFT calls. Evaluate this nonnegative radius DAG with outward interval operations: integer powers use finite products, and the normalization square root uses the same range format and index $J$. As with the existing model and rational constants, use certified enclosing constants throughout; in particular enclose the selected model's valid unit-roundoff bound. If that model supplies no tighter certified enclosure for $u$, use the already verified enclosure $[0,1\mathrm{/}16]$, converted outward to the selected format; thus no additional oracle or input for $u$ is required. Reject the candidate as inconclusive if any radius enclosure is non-finite or the normalization denominator enclosure is not strictly positive before division. Define the returned coordinate radius by the upward dyadic sum $r_k=\mathrm{ru}(b_k+R_{\mathrm{F}\mathrm{F}\mathrm{T}})$; if this sum overflows, reject the candidate with its range status. Thus every returned radius includes both the scaling/deconvolution rectangle and the separately bounded error of its computed FFT relative to the exact transform of the dyadic center grid. Set ${ϵ}_{\mathrm{a}\mathrm{n}\mathrm{a}}=A^{+}+\frac{T_d^{+}}{P{d}_{\min }^{}},B_{\mathrm{r}\mathrm{o}\mathrm{u}\mathrm{n}\mathrm{d}}={\max }_{k\in \mathcal{K}}{r}_k+\frac{E_{\mathrm{e}\mathrm{v}\mathrm{a}\mathrm{l}}}{\sqrt{P}{d}_{\min }^{}}\mathrm{.}$ The returned error certificate is ${ϵ}_{\mathrm{a}\mathrm{n}\mathrm{a}}$ and $B_{\mathrm{r}\mathrm{o}\mathrm{u}\mathrm{n}\mathrm{d}}$.

4. For type 2, form the conjugate-transpose finite gridding factorization. For each input mode $k$, evaluate its multiplier interval using the same format and truncation setting $N=J$; if its lower endpoint is $\le 0$, reject the candidate as inconclusive before division. Otherwise divide by that certified positive interval. For the resulting mode rectangles choose lower-left dyadic centers $v_k$ and radii $e_k$ enclosing their deviations. Set $E_f$ to an outward upper bound for $({\sum }_k{e}_k^2{)}^{1\mathrm{/}2}$ and $V$ to an outward upper bound for $\mathrm{\parallel }v{\mathrm{\parallel }}_2$. Require every real and imaginary coordinate of the $v_k$ to belong to $\mathbb{F}$; otherwise reject this candidate with a range-failure status before calling the FFT. Embed each $v_k$ at grid slot $\kappa (k)$, zero-fill all unused grid modes with zero centers and zero radii, and apply Computed unitary radix-two FFT with certified twiddle phases in direction inverse on the same axes and in the same lexicographic line order, passing the same operationwise model and phase controls as in type 1. Reject a non-finite or inconclusive FFT call; for a finite return verify every returned twiddle certificate ${ϵ}_{\mathrm{\ell }}\le {ϵ}_{\mathrm{t}\mathrm{w}}$, with the same literal zero-bound rule. Enlarge the unitary-grid output by radius $E_f+\rho V$ and scale by $P^{-1\mathrm{/}2}$. For each source $x_j$, use exactly the grid cells $\mathrm{\ell }\in W_j$ from step 2; reduce each coordinate displacement ${\mathrm{\ell }}_r\mathrm{/}{M}_r-x_{j,r}$ modulo one to $z_r\in [0,1)$, enclose ${\sum }_{n=-L_r}^{L_r}\exp (-\pi (z_r+n{)}^2\mathrm{/}{\tau }_r)$ by range-evaluator calls with format $(p,e_{\min },e_{\max })$ and $N=J$, and add $[-R_r^{+},R_r^{+}]$. Multiply the coordinate window intervals, conjugate the resulting real window (which leaves it unchanged), and multiply and sum against the scaled grid rectangles using outward complex-rectangle arithmetic at every operation. Return the lower-left corner $z_j$ and an outward radius $r_j$ containing each resulting rectangle; the type-2 rounding certificate is $B_{\mathrm{r}\mathrm{o}\mathrm{u}\mathrm{n}\mathrm{d}}^{\ast }={\max }_j{r}_j$. Its analytic certificate is the same ${ϵ}_{\mathrm{a}\mathrm{n}\mathrm{a}}$ from step 3.

5. Enclose $\mathrm{\parallel }c{\mathrm{\parallel }}_1$ or $\mathrm{\parallel }f{\mathrm{\parallel }}_1$ from above by outward modulus and addition intervals. For each candidate whose operations succeed, compare the outward upper bound ${ϵ}_{\mathrm{a}\mathrm{n}\mathrm{a}}\mathrm{\parallel }c{\mathrm{\parallel }}_1+B_{\mathrm{r}\mathrm{o}\mathrm{u}\mathrm{n}\mathrm{d}}$ (type 1), or ${ϵ}_{\mathrm{a}\mathrm{n}\mathrm{a}}\mathrm{\parallel }f{\mathrm{\parallel }}_1+B_{\mathrm{r}\mathrm{o}\mathrm{u}\mathrm{n}\mathrm{d}}^{\ast }$ (type 2), exactly with ${\epsilon }_{\mathrm{r}\mathrm{e}\mathrm{q}}$. If it is at most ${\epsilon }_{\mathrm{r}\mathrm{e}\mathrm{q}}$, return the centers and their certificates with status certified. Otherwise record this candidate's first non-success reason as requested-tolerance-not-met together with its finite bound, then continue to the next candidate. For any candidate rejected by a range, interval, denominator, mode-positivity, FFT, or twiddle check, record its first failure status and location, then continue to the next candidate. If the finite candidate list is exhausted without a certified return, return status inconclusive with the first recorded non-success record in candidate order. In particular, if every candidate had finite valid bounds above ${\epsilon }_{\mathrm{r}\mathrm{e}\mathrm{q}}$, return inconclusive(requested-tolerance-not-met) with the first such bound; no approximation is returned as certified on this branch.

The analytic term consists only of the exact infinite-grid alias contribution and the omitted finite-spread tail. The reported rounding radius contains the finite image-series enclosure, outward arithmetic, and computed FFT error. The signed-mode slot map $\kappa$ is used for both type-1 extraction and type-2 embedding. The transform and FFT axes are fixed; a different normalization or sign convention is a different procedure.

The unitary radix-two FFT error theorem is used only on successful axes with the stated $u$ and twiddle bounds; range, interval, or denominator failures return inconclusive. The global lower bound $d_{\min }^{}>0$ does not replace checking each separately evaluated mode interval before division. The finite candidate list is part of the input; if none meets the requested tolerance the returned status is inconclusive. This procedure handles types 1 and 2 only. Type 3 requires its own frequency-space scale/product and complexity analysis.