Ozaki-Bailey FFT Methods

• Ozaki-Bailey FFT is a class of reconstruction methods that reformulate problems to recover coefficients or transform values via fast Fourier techniques.
• Techniques include coefficient extraction from roots using FFT inversion, binary-extension-field transforms for Reed-Solomon codes, and Bluestein-based exact cyclic convolutions.
• Recent GPU implementations combine Bailey’s six-step decomposition with tensor-core CRT and Kulisch fixed-point accumulation to optimize performance.

Ozaki-Bailey FFT denotes a family of FFT-mediated reconstruction strategies in which a target object—polynomial coefficients, transform values over binary extension fields, or full target-precision Fourier outputs—is recovered from structured evaluations or split subcomputations at transform-friendly nodes. In the recent GPU literature, the term is used explicitly for an emulated 3-D FFT that combines Ozaki Scheme II with the Bailey six-step decomposition (Matsuoka, 28 May 2026). Closely related constructions use the same underlying pattern for coefficient extraction from roots (Thomann, 2016), for FFTs over $\mathbb{F}_{2^m}$ in Reed-Solomon coding (Lin et al., 2015), and for target-precision FFTs obtained from lower-precision exact cyclic convolutions via Bluestein, NTTs, and CRT (Kawakami et al., 31 Mar 2026).

1. Structured evaluation and reconstruction pattern

A recurring feature of Ozaki-Bailey-style methods is that the Fourier or coefficient problem is not attacked in its original form. Instead, the computation is first recast so that a structured transform becomes available. In the polynomial-coefficient setting, the coefficients are recovered by inverse discrete Fourier transform from polynomial values sampled at roots of unity. In the binary-extension-field setting, FFT and IFFT are defined over affine subspaces $V_k+\beta$ using a basis built from subspace polynomials. In the target-precision setting, Bluestein converts the DFT into a cyclic convolution, after which Ozaki-style splitting, NTT evaluation, and CRT reconstruction are applied. In the GPU setting, Bailey’s six-step FFT decomposes the transform into small dense DFT stages that can be expressed as GEMMs and then emulated through FP8 tensor cores with CRT/Garner reconstruction (Thomann, 2016, Lin et al., 2015, Kawakami et al., 31 Mar 2026, Matsuoka, 28 May 2026).

This shared structure does not imply a single implementation template. The evaluation nodes, arithmetic domains, and reconstruction mechanisms differ substantially. Roots of unity are used in coefficient recovery; affine subspaces and linearized polynomials are used over $\mathbb{F}_{2^m}$; cyclic convolution is induced by Bluestein in the lower-precision scheme; and modular tensor-core residue computations with Garner or Kulisch reduction are used in the GPU formulation. This suggests that the unifying content of the term is architectural rather than notational: a problem is reformulated so that a fast transform or transform-like subroutine becomes the central recovery mechanism.

2. Coefficient extraction from roots and symmetric functions

For the polynomial

$$P(\omega)=\prod_{k=1}^n(\omega-\omega_k)=\sum_{m=0}^n a_m\omega^m,$$

the FFT-based coefficient algorithm of Katsifarakis and Papageorgiou converts “coefficients from roots” into a discrete Fourier inversion problem. One chooses an integer $N$ with $n<N\le 2n$, evaluates

$$p_j=P(\omega_N^{-j})=\prod_{k=1}^n(\omega_N^{-j}-\omega_k), \qquad j=0,\dots,N-1,$$

and then recovers the coefficients by

$$a_m=\frac{1}{N}\sum_{j=0}^{N-1}P(\omega_N^{-j})\,\omega_N^{jm}, \qquad 0\le m\le n.$$

This is the paper’s Algorithm $P$: choose $N$, compute the $V_k+\beta$0, apply the FFT, and output the first $V_k+\beta$1 Fourier coefficients. The implementation uses radix-2 Cooley-Tukey FFT, while Mathematica’s FFT can allow $V_k+\beta$2. The product stage costs about $V_k+\beta$3 complex multiplications/additions, the FFT costs $V_k+\beta$4, and with radix-2 and $V_k+\beta$5 the leading complexity is $V_k+\beta$6 floating-point operations with $V_k+\beta$7 storage for tasks other than Vandermondian inversion. For complex roots the paper states roughly between $V_k+\beta$8 and $V_k+\beta$9 floating-point operations, and for real roots roughly between $\mathbb{F}_{2^m}$0 and $\mathbb{F}_{2^m}$1. A notable property is that the method is insensitive to root ordering (Thomann, 2016).

The same paper contrasts Algorithm $\mathbb{F}_{2^m}$2 with the recursive Algorithm $\mathbb{F}_{2^m}$3,

$\mathbb{F}_{2^m}$4

and with the Leja-enhanced variant $\mathbb{F}_{2^m}$5. Algorithm $\mathbb{F}_{2^m}$6 also has $\mathbb{F}_{2^m}$7 time and $\mathbb{F}_{2^m}$8 storage, with leading floating-point cost about $\mathbb{F}_{2^m}$9 for complex roots and $$P(\omega)=\prod_{k=1}^n(\omega-\omega_k)=\sum_{m=0}^n a_m\omega^m,$$0 for real roots, but it is very sensitive to root ordering. Leja ordering is introduced through

$$P(\omega)=\prod_{k=1}^n(\omega-\omega_k)=\sum_{m=0}^n a_m\omega^m,$$1

The main empirical conclusion is that for roots on or near circles in the complex plane, Algorithm $$P(\omega)=\prod_{k=1}^n(\omega-\omega_k)=\sum_{m=0}^n a_m\omega^m,$$2 is generally superior to $$P(\omega)=\prod_{k=1}^n(\omega-\omega_k)=\sum_{m=0}^n a_m\omega^m,$$3, especially for larger $$P(\omega)=\prod_{k=1}^n(\omega-\omega_k)=\sum_{m=0}^n a_m\omega^m,$$4 and in the reported tests up to and beyond $$P(\omega)=\prod_{k=1}^n(\omega-\omega_k)=\sum_{m=0}^n a_m\omega^m,$$5. The paper repeatedly reports that $$P(\omega)=\prod_{k=1}^n(\omega-\omega_k)=\sum_{m=0}^n a_m\omega^m,$$6 is robust for roots on circles, including unit roots; that $$P(\omega)=\prod_{k=1}^n(\omega-\omega_k)=\sum_{m=0}^n a_m\omega^m,$$7 can be very accurate for small problems but tends to deteriorate quickly as $$P(\omega)=\prod_{k=1}^n(\omega-\omega_k)=\sum_{m=0}^n a_m\omega^m,$$8 increases; and that on the real line both methods degrade, though $$P(\omega)=\prod_{k=1}^n(\omega-\omega_k)=\sum_{m=0}^n a_m\omega^m,$$9 may sometimes be slightly better. Scaling by $N$0 changes coefficients by $N$1, and for roots inside a small circle of radius $N$2, choosing $N$3 can improve the accuracy of $N$4 by $N$5–$N$6 orders of magnitude, whereas scaling is generally not helpful for $N$7.

Because the coefficients are, up to sign, the elementary symmetric functions of the roots, the same method computes symmetric functions efficiently. The paper also applies the coefficient stage to reduced polynomials, Vandermonde inversion, and interpolation. In Vandermonde inversion, replacing the classical recursive coefficient stage by Algorithm $N$8 yields the method $N$9, and the paper reports that $n<N\le 2n$0 can accurately invert Vandermonde matrices of dimension $n<N\le 2n$1 and more when roots lie on the unit circle. In interpolation, the alternative algorithm $n<N\le 2n$2 first forms the barycentric weights

$n<N\le 2n$3

then evaluates the interpolation polynomial at the $n<N\le 2n$4-th roots of unity and finally recovers the coefficients using the FFT exactly as in Algorithm $n<N\le 2n$5. The reported complexity is about $n<N\le 2n$6 complex additions/multiplications with $n<N\le 2n$7 storage, and $n<N\le 2n$8 is described as superior to the conventional Vandermonde-solve approach $n<N\le 2n$9 in numerical performance.

3. Binary-extension-field transforms and Reed-Solomon decoding

A distinct Ozaki-Bailey-style formulation appears in the FFT over binary extension fields developed for systematic $$p_j=P(\omega_N^{-j})=\prod_{k=1}^n(\omega_N^{-j}-\omega_k), \qquad j=0,\dots,N-1,$$0 Reed-Solomon codes over $$p_j=P(\omega_N^{-j})=\prod_{k=1}^n(\omega_N^{-j}-\omega_k), \qquad j=0,\dots,N-1,$$1, with $$p_j=P(\omega_N^{-j})=\prod_{k=1}^n(\omega_N^{-j}-\omega_k), \qquad j=0,\dots,N-1,$$2 a power of two. The construction uses a nonstandard basis built from subspace polynomials over the nested chain

$$p_j=P(\omega_N^{-j})=\prod_{k=1}^n(\omega_N^{-j}-\omega_k), \qquad j=0,\dots,N-1,$$3

with

$$p_j=P(\omega_N^{-j})=\prod_{k=1}^n(\omega_N^{-j}-\omega_k), \qquad j=0,\dots,N-1,$$4

and

$$p_j=P(\omega_N^{-j})=\prod_{k=1}^n(\omega_N^{-j}-\omega_k), \qquad j=0,\dots,N-1,$$5

The basis is defined by

$$p_j=P(\omega_N^{-j})=\prod_{k=1}^n(\omega_N^{-j}-\omega_k), \qquad j=0,\dots,N-1,$$6

and each basis element satisfies $$p_j=P(\omega_N^{-j})=\prod_{k=1}^n(\omega_N^{-j}-\omega_k), \qquad j=0,\dots,N-1,$$7. This degree-ordered property is what enables a divide-and-conquer FFT over affine subspaces $$p_j=P(\omega_N^{-j})=\prod_{k=1}^n(\omega_N^{-j}-\omega_k), \qquad j=0,\dots,N-1,$$8 (Lin et al., 2015).

The transform recursion exploits

$$p_j=P(\omega_N^{-j})=\prod_{k=1}^n(\omega_N^{-j}-\omega_k), \qquad j=0,\dots,N-1,$$9

so that a length-$$a_m=\frac{1}{N}\sum_{j=0}^{N-1}P(\omega_N^{-j})\,\omega_N^{jm}, \qquad 0\le m\le n.$$0 transform reduces to two length-$$a_m=\frac{1}{N}\sum_{j=0}^{N-1}P(\omega_N^{-j})\,\omega_N^{jm}, \qquad 0\le m\le n.$$1 transforms on shifted evaluation sets. For $$a_m=\frac{1}{N}\sum_{j=0}^{N-1}P(\omega_N^{-j})\,\omega_N^{jm}, \qquad 0\le m\le n.$$2, the addition and multiplication recurrences are

$$a_m=\frac{1}{N}\sum_{j=0}^{N-1}P(\omega_N^{-j})\,\omega_N^{jm}, \qquad 0\le m\le n.$$3

hence

$$a_m=\frac{1}{N}\sum_{j=0}^{N-1}P(\omega_N^{-j})\,\omega_N^{jm}, \qquad 0\le m\le n.$$4

and both FFT and IFFT run in $$a_m=\frac{1}{N}\sum_{j=0}^{N-1}P(\omega_N^{-j})\,\omega_N^{jm}, \qquad 0\le m\le n.$$5.

The coding application is central. With $$a_m=\frac{1}{N}\sum_{j=0}^{N-1}P(\omega_N^{-j})\,\omega_N^{jm}, \qquad 0\le m\le n.$$6, encoding and syndrome computation are reorganized so that only transforms of length $$a_m=\frac{1}{N}\sum_{j=0}^{N-1}P(\omega_N^{-j})\,\omega_N^{jm}, \qquad 0\le m\le n.$$7 are needed for the parity or syndrome part, avoiding a full length-$$a_m=\frac{1}{N}\sum_{j=0}^{N-1}P(\omega_N^{-j})\,\omega_N^{jm}, \qquad 0\le m\le n.$$8 transform in those stages. The syndrome polynomial is extracted from blockwise inverse transforms, while the error-locator polynomial is obtained through a fast half-GCD algorithm on degree-$$a_m=\frac{1}{N}\sum_{j=0}^{N-1}P(\omega_N^{-j})\,\omega_N^{jm}, \qquad 0\le m\le n.$$9 polynomials with complexity $P$0. The overall decoding pipeline has complexity

$P$1

improving upon the cited $P$2 benchmark and reaching the best known complexity bound that was established by Justesen in 1976, but now over the practically important family $P$3. The authors contrast this with Cooley-Tukey FFTs, which apply only when the field contains enough roots of unity, and with arbitrary-field FFTs of Schönhage/Cantor type, which incur extra $P$4 factors. Computer simulations report that the proposed decoder is $P$5 times faster than the conventional one for the $P$6 Reed-Solomon code over $P$7.

4. Target-precision FFTs from lower-precision exact convolutions

Another major development formulates an Ozaki-Bailey-style FFT for target precision by applying the Ozaki scheme to the cyclic convolution inside Bluestein’s FFT. For a length-$P$8 DFT,

$P$9

Bluestein rewrites the transform as

$N$0

The method then treats the induced cyclic convolution as the computational core, splits the inputs into lower-precision pieces in Ozaki fashion,

$N$1

and computes every split convolution exactly using 32-bit NTTs together with CRT reconstruction rather than floating-point FFTs (Kawakami et al., 31 Mar 2026).

The exactness condition is what distinguishes this scheme from ordinary floating-point convolution. For a single modulus $N$2, exact recovery requires

$N$3

which yields

$N$4

With two moduli $N$5, the split width becomes

$N$6

The implementation uses the 32-bit primes

$N$7

This larger admissible $N$8 is identified as the key reason that the NTT/CRT route yields fewer splits than a floating-point exact-rounding route based on FFT error bounds.

The paper also introduces NTT-domain accumulation. If two split-pair products share the same scale, their transform-domain products can be accumulated before the inverse NTT, reducing the inverse-transform count. With capped splits $N$9, the total NTT plus inverse NTT count is at most $V_k+\beta$00; with $V_k+\beta$01, it becomes $V_k+\beta$02 for $V_k+\beta$03 and $V_k+\beta$04 across all tested lengths and matches the stated theoretical minimum $V_k+\beta$05 when $V_k+\beta$06. The method is implemented as a double-precision Bluestein FFT using TS arithmetic for outer complex arithmetic and 32-bit NTTs for exact split convolutions.

The numerical profile is mixed in a characteristic way. Against FFTW double precision, Stockham double precision, Bluestein double precision, TS Stockham FFT, and TS Bluestein FFT, the proposed method usually has the smallest relative error; it does not show a strong growth trend with $V_k+\beta$07; and it substantially improves over TS Bluestein, indicating that the exact split convolutions suppress error accumulation. At the same time, it does not always give the smallest maximum relative error. On an Intel Xeon Platinum 8468 for lengths $V_k+\beta$08–$V_k+\beta$09, the execution time is approximately $V_k+\beta$10–$V_k+\beta$11 that of FFTW’s double-precision FFT, and the paper states that 32-bit NTTs account for approximately $V_k+\beta$12 of the total time.

5. Bailey six-step decomposition, tensor-core CRT, and the Kulisch route

In the 2026 GPU formulation, Ozaki-Bailey FFT is defined explicitly as an emulated 3-D FFT obtained by combining Ozaki Scheme II with Bailey’s six-step decomposition. The hardware motivation is that NVIDIA’s Blackwell Ultra B300 cuts FP64 vector throughput to about $V_k+\beta$13 TFLOPS per GPU, roughly $V_k+\beta$14 below B200, while FFT is treated as a bandwidth-limited primitive that should run near the memory roof on a well-balanced GPU. For a $V_k+\beta$15 3-D FFT, the paper computes memory traffic

$V_k+\beta$16

for $V_k+\beta$17, giving a memory-roof time of about $V_k+\beta$18 ms at $V_k+\beta$19 TB/s HBM bandwidth. Ozaki Scheme II is then applied to the GEMM-represented inner DFT stages of Bailey’s decomposition, routing them through FP8 tensor cores with mantissa-sliced Chinese-remainder reconstruction (Matsuoka, 28 May 2026).

The Bailey six-step factorization uses $V_k+\beta$20; for $V_k+\beta$21, the natural choice is $V_k+\beta$22. The two 1-D FFT stages are implemented as dense DFT matrix multiplications, so in the 3-D case the transform reduces to repeated small GEMM-like kernels. For B300-class hardware, the paper recommends $V_k+\beta$23 moduli for FP64-equivalent accuracy, with tensor-core cost multiplier

$V_k+\beta$24

The crucial observation is that Bailey’s inner dimension $V_k+\beta$25 gives $V_k+\beta$26 for $V_k+\beta$27, while $V_k+\beta$28. Hence the regime is $V_k+\beta$29, so reconstruction overhead no longer amortizes and the third Tensor-Memory Equilibrium parameter, $V_k+\beta$30, becomes the binding cost.

Standard Garner reconstruction is estimated at roughly $V_k+\beta$31 INT32 operations per output, about $V_k+\beta$32 INT32 operations per output for $V_k+\beta$33, and about $V_k+\beta$34 INT32 operations over all FFT outputs. On B300, the effective INT32 throughput for this modular work is estimated at about $V_k+\beta$35 TOPS after Barrett-reduction overhead, yielding a reconstruction time of about $V_k+\beta$36 ms, far above the $V_k+\beta$37 ms memory roof. The paper therefore reformulates CRT reconstruction as a forward CRT sum,

$V_k+\beta$38

and slices each $V_k+\beta$39 into $V_k+\beta$40 8-bit chunks. This creates Phase A, an $V_k+\beta$41 inner-product GEMM on tensor cores, and Phase B, the per-output reduction. For $V_k+\beta$42, the main-text choice is $V_k+\beta$43. Phase A is reported as around $V_k+\beta$44–$V_k+\beta$45 ms; Phase B is the bottleneck if executed naively in FP64.

The paper’s “Kulisch escape route” replaces Phase B by exact fixed-point accumulation on the INT32 SIMT pipe. The reconstruction is written as

$V_k+\beta$46

and because $V_k+\beta$47, each term lands at a fixed bit position. For $V_k+\beta$48, the running sum needs about $V_k+\beta$49 bits, so a $V_k+\beta$50-bit accumulator made from five INT32 registers is sufficient. The implementation model uses five INT32 registers, about four INT32 operations per slice, and roughly $V_k+\beta$51 INT32 operations per output for $V_k+\beta$52. The reported B300 estimate is about $V_k+\beta$53 ms at peak INT32 throughput and $V_k+\beta$54–$V_k+\beta$55 ms at realistic $V_k+\beta$56–$V_k+\beta$57 efficiency.

The resulting parity analysis is expressed as closed-form floors. For FFT, the operational intensity is

$V_k+\beta$58

which for $V_k+\beta$59 gives $V_k+\beta$60, hence the native FP64 floor $V_k+\beta$61. The naive Ozaki-Bailey FP64 floor is

$V_k+\beta$62

which is $V_k+\beta$63 for $V_k+\beta$64. The Kulisch INT32 sub-floor is

$V_k+\beta$65

with an overlap-optimized version reducing this to roughly $V_k+\beta$66. The FP8 floor is $V_k+\beta$67. The projected performance for a $V_k+\beta$68 3-D FFT is then: memory roof $V_k+\beta$69 ms, native B300 FP64 about $V_k+\beta$70 ms, recursive Garner Ozaki-Bailey about $V_k+\beta$71 ms, tensor-core Garner plus naive FP64 Phase B about $V_k+\beta$72 ms, and tensor-core Garner plus Kulisch Phase B about $V_k+\beta$73 ms. The paper is explicit that these are projections, not measured results.

6. Performance regimes, applications, and misconceptions

The various Ozaki-Bailey formulations share a common design logic, but they behave differently because different stages dominate in different regimes. In polynomial coefficient extraction, the FFT is only the final recovery step; the expensive part is the product evaluation

$V_k+\beta$74

and the method is most effective for roots on or near circles (Thomann, 2016). In the binary-extension-field transform, the main gain comes from a basis and recursion that make FFT, IFFT, multiplication, division, and half-GCD compatible over $V_k+\beta$75, allowing Reed-Solomon decoding complexity $V_k+\beta$76 (Lin et al., 2015). In the lower-precision Bluestein scheme, exactness is transferred to the split cyclic convolutions, so performance is governed by NTT counts and NTT efficiency rather than by a conventional FFT kernel (Kawakami et al., 31 Mar 2026). In the GPU formulation, the central issue is not the Bailey decomposition itself but the per-output reconstruction latency: FFT lies in the $V_k+\beta$77 regime, so $V_k+\beta$78 binds unless Phase B is reformulated through Kulisch accumulation (Matsuoka, 28 May 2026).

Several common misconceptions are contradicted by these results. One is that FFT-based recovery is determined only by asymptotic complexity. The coefficient-from-roots paper states that algorithm quality is not determined only by condition number and emphasizes the role of the evaluation-point structure; circle-like distributions are favorable, whereas line-like distributions are much harder (Thomann, 2016). Another is that lower-precision subcomputations necessarily degrade accuracy. The Bluestein-NTT method obtains a double-precision FFT from 32-bit NTTs and reports reduced relative error compared with FFTW double precision and Triple-Single precision arithmetic, while the GPU Kulisch route is presented as preserving full FP64 accuracy through exact fixed-point accumulation with one final conversion (Kawakami et al., 31 Mar 2026, Matsuoka, 28 May 2026). A further misconception is that Cooley-Tukey is the only viable FFT template in coding applications. The binary-extension-field work explicitly constructs a true $V_k+\beta$79 FFT over $V_k+\beta$80 without requiring multiplicative roots of unity, precisely to avoid that restriction (Lin et al., 2015).

Taken together, these results indicate that “Ozaki-Bailey FFT” is best understood as a class of transform-centric reconstruction methods whose practicality depends on where the reformulation places the numerical and architectural burden. In some settings the decisive factor is node geometry and scaling; in others it is basis design over finite fields; in others it is exact modular convolution; and in current GPU work it is whether CRT reduction can be moved off the collapsed FP64 vector path.