Rigorous numerical DFTs in fixed point

Ball arithmetic [Hoe2009; Joh2017] gives rigorous results for any computation, but for a DFT it is wasteful: every operation propagates a radius, the radii grow with the depth of the transform, and each multiplication costs several floating-point operations. Since all roots of unity have modulus 1, the magnitudes in a DFT are easy to bound in advance, which makes fixed-point arithmetic a better fit. gr_dft runs transforms over two gr contexts of fixed-point numbers, computes a bound for every intermediate value and for the output error from the structure of the plan alone, and uses these bounds to produce enclosures for transforms of acb vectors (gr_dft_acb, which also implements acb_dft). For floating-point FFTs, error bounds of a similar kind are classical [Hig2002, ch. 24; Per2003].

Arithmetic#

With $k$ limbs of $\beta = \mathtt{FLINT\_BITS}$ bits, a real fixed-point number is $t = \pm f\, 2^{-k\beta}$ with an integer $0 \le f < 2^{k\beta}$, so $|t| < 1$, and we write $\mathrm{ulp} = 2^{-k\beta}$. A complex number is a pair of real ones. All error bounds below are in ulp and, for complex values, bound the modulus of the error unless stated otherwise.

The error of a real product is less than 1 ulp, but the analysis allows 2, which covers flint_mpn_mulhigh_n; for $k \le 2$ the products are exactly truncated, and a complex product combines two of them exactly. The three-product formula from 22 limbs computes $(x_r + x_i)(y_r + y_i) - x_r y_r - x_i y_i$ with the sums held in an extra limb, so its imaginary part combines three rounded products. Overflow is detected from the carries, so a computation that returns GR_SUCCESS has kept every value in $(-1, 1)$. The constants $\pm 1$ are not representable and are replaced by $\pm(1 - \mathrm{ulp})$; the transforms never multiply by them, since in complex mode the powers $\pm 1$, $\pm i$ of the root are applied by exchanging and negating parts.

The table of roots#

The plan's table holds approximations $\hat w_j$ of $w^j = e^{-2\pi i j/n}$. The primitive root is computed with arb at $k\beta + 64$ bits and truncated, with error at most 2 ulp per part. The entries up to $L \approx n/8$ (or $n/4$, or $n/2$, if $4 \nmid n$ or $2 \nmid n$) are computed as $\hat w_j = \hat w_{\lfloor j/2 \rfloor} \hat w_{\lceil j/2 \rceil}$, and the others are exact copies with exchanged and negated parts, by the symmetries $w^{n/4 - j} = -i\, \overline{w^j}$, $w^{n/2 - j} = -\overline{w^j}$ and $w^{n - j} = \overline{w^j}$.

The recurrence holds up to second-order terms, because $|\hat w_a \hat w_b - w^a w^b| \le |\hat w_a - w^a| |\hat w_b| + |\hat w_b - w^b|$ and $|\hat w_b| \le 1 + E(b)\,\mathrm{ulp}$; the implementation inflates the bound by a factor $1 + 10^{-7}$ to cover them. The resulting bound $D_w$ for the whole table, about $(2 + K) L$, is stored in the plan and is $O(n)$ ulp. Sub-plans whose roots are powers of the parent root take strided copies of the parent table and its bound.

Propagating bounds through a plan#

The bound computation gr_dft_precomp_nfixed_bound follows the plan and tracks two numbers through the transform: a bound $B$ for the modulus of the current values and a bound $D$ for the modulus of their errors, in ulp. A running maximum, the peak, records the largest modulus reached by any input, intermediate or output value.

Every algorithm is a sequence of these steps, so its bounds are affine functions of the input bounds $(B_0, D_0)$, and for recursive algorithms they are computed with memoised coefficients in $O(\log n)$ operations on doubles, each inflated by a factor $1 + 10^{-10}$ to cover the rounding of the doubles themselves.

The other algorithms are handled in the same way: the split-radix recursion by memoised affine coefficients per depth, the direct prime kernel of Lemma 5 of Discrete Fourier transforms over a generic ring by bounding its $p - 1$ differences of modulus $2B$ and their rotations, the prime factor and four-step algorithms by composing the bounds of their sub-plans, and the naive algorithm as $n$ rotated terms. Since $D_w = O(n)$, the output error of a radix-2 transform is $O(n^2)$ ulp for inputs of modulus about $1/n$, compared with the $O(\log n)$ relative error of a floating-point FFT with accurate roots [Hig2002, ch. 24]; the cost is a few guard bits per doubling of $n$, about 3 in practice.

Bluestein's algorithm#

A generic bound on Bluestein's algorithm would be too weak: the inverse convolution transform doubles the bound at each of its $\log_2 N$ levels, starting from values that are already $N$ times the input. Here the kernel $\hat h$ is scaled by $1/(2N)$, so $|\hat h_k| \le (2n - 1)/(2N) < 1/2$, and the values inside the inverse transform are bounded directly: each intermediate value of the unscaled inverse transform at level $s$ is $2^s$ times a sum of $N/2^s$ aliased values of the exact convolution $u = z \ast h$, and $|u_k| \le \|z\|_1 \max|h| \le n B_z/(2N)$. All intermediate values of the inverse transform are therefore at most $n B_z / 2$, and the final doubling and chirp multiplication give outputs bounded by $n B_z$, as for any DFT. The errors follow the generic rules.

From bounds to enclosures#

The bounds are linear in $B_0$, so for unit inputs the plan has a peak $P_1 \ge n$, and inputs of modulus less than $1/P_1$ cannot overflow. gr_dft_acb uses this as follows, for working precision $p$.

1. With $P_1 < 2^{e_1}$, the input scale is $m = 2^{-e_1 - 1}$, and $E$ is the output error bound for inputs bounded by $m$ with errors of 1 ulp (the truncation of the input).
2. The number of limbs $k$ is the least with $k\beta \ge p$ and $E\, 2^{-k\beta} < 2^{-p-2} m$, so that the roundoff is below $2^{-p}$ relative to the largest input.
3. For an input whose midpoint parts are less than $2^{E_M}$ in absolute value, the midpoints are multiplied by $2^{-e}$, $e = E_M + e_1 + 1$, truncated to $k$ limbs and transformed. The outputs are multiplied by $2^e$, rounded to $p$ bits, and get the radius $E\, 2^{e - k\beta} + \rho$ plus their rounding error, where $\rho$ is the sum of all input radii.

The output therefore contains the exact transform of every point of the input balls: the transform of the midpoints is within the roundoff bound, and moving from the midpoints to any point of the balls changes the result by at most $\rho$. The bound $\rho$ is computed once with upward rounding, which is both cheaper and, for inputs with radii of similar size, sharper than propagating radii through $\log n$ levels of ball arithmetic. For the inverse transform, the radius is divided by $n$, exactly by an exponent shift for powers of two and otherwise together with a single rounded division of each output midpoint.

@misc{flint-doc-algorithms-fft-fixed-point-dft,
  author       = {{The FLINT team}},
  title        = {{Rigorous numerical DFTs in fixed point}},
  howpublished = {FLINT: Fast Library for Number Theory, documentation, version 3.7.0-dev},
  year         = {2026},
  url          = {https://flintlib.org/doc/algorithms/fft/fixed-point-dft.html},
  note         = {Generated 2026-11-02}
}