Kulisch Fixed-Point Complete Arithmetic

• Kulisch fixed-point complete arithmetic is a method that uses wide fixed-point registers to achieve exact accumulation with full FP64 accuracy in the FFT reconstruction phase.
• It reformulates Bailey decomposition reconstruction by leveraging byte-sliced integer arithmetic and carry propagation, thereby mitigating the FP64 pipeline bottleneck.
• The technique transfers critical reductions to high-throughput INT32 SIMT units, ensuring efficient 3-D FFT processing with error bounds limited to about 2 ulps.

Kulisch fixed-point complete arithmetic, in the formulation developed for the Ozaki-Bailey FFT on Blackwell-class GPUs, is a wide fixed-point exact-accumulation method used exclusively in the reconstruction stage of an emulated FP64 3-D FFT. Its role is not to accelerate FFTs in general, but to eliminate a specific reconstruction bottleneck that appears when an Ozaki Scheme II implementation of a Bailey six-step FFT is mapped onto hardware whose native FP64 vector throughput is too weak to sustain memory-roof execution. In that setting, Kulisch arithmetic becomes a “Kulisch escape route”: it reroutes the critical Phase B reduction away from the collapsed FP64 vector pipe and onto the INT32 SIMT pipeline while preserving exact accumulation and a single final rounding at readout (Matsuoka, 28 May 2026).

1. Conceptual definition and problem setting

In the relevant use case, the computational target is a full-FP64-equivalent $1024^3$ 3-D FFT on a GPU such as NVIDIA’s B300. The motivation is architectural: B300 is described as having native FP64 vector throughput of only $\sim 1.3$ TFLOPS/GPU with HBM bandwidth of $\sim 8$ TB/s, whereas the native FP64 bandwidth-parity floor for FFT is

$$\eta_{\text{FP64}^{\text{native}}}=\mathrm{OI}_{\text{FFT}}\cdot B_{\mathrm{HBM}}\approx 1.56\,B_{\mathrm{HBM}},$$

which evaluates to $12.5$ TFLOPS at $8$ TB/s. In the paper’s framing, B300 is therefore about $10\times$ below the native floor, so a native FP64 FFT cannot remain memory-bound (Matsuoka, 28 May 2026).

The recovery path is Ozaki Scheme II. In this framework, FP64 GEMM-like work is replaced by modular integer arithmetic over several small primes, executed on lower-precision fast hardware and reconstructed afterward. The summary equation is

$$C=\mathrm{Garner}(\tilde C^{(1)},\dots,\tilde C^{(r)};\{m_i\})/(s_A s_B),$$

with scaled integer inputs

$$\tilde A=\mathrm{round}(s_A\cdot A),\qquad \tilde B=\mathrm{round}(s_B\cdot B).$$

The significance of Kulisch arithmetic emerges only after this modular decomposition has already been adopted: it is a reconstruction primitive within Ozaki-II, not a substitute for the overall emulation strategy.

The same source is explicit that the intended outcome is not strict bitwise exactness of the entire FFT pipeline, but full FP64 accuracy in the sense of exact accumulation followed by final-rounding effects bounded at the output. A recurring misconception is therefore that Kulisch arithmetic here replaces all floating-point computation. It does not. Rather, it isolates one reduction phase whose structure permits exact fixed-point accumulation and uses that structure to avoid the weak FP64 vector units.

2. Bailey decomposition and the origin of the reconstruction bottleneck

The FFT embedding is based on a Bailey six-step decomposition. For a length-$N=pq$ 1-D FFT, the decomposition is transpose $\sim 1.3$0 length-$\sim 1.3$1 column FFTs $\sim 1.3$2 twiddle multiplication $\sim 1.3$3 transpose $\sim 1.3$4 length-$\sim 1.3$5 column FFTs $\sim 1.3$6 permutation. The implementation consequence is that the two small 1-D FFT stages can be expressed as dense GEMMs with DFT matrices, which makes them amenable to Ozaki/tensor-core acceleration. For $\sim 1.3$7, the natural factorization is $\sim 1.3$8 (Matsuoka, 28 May 2026).

This factorization creates the key difficulty. In the Tensor–Memory Equilibrium model, the emulated execution time is

$\sim 1.3$9

where $\sim 8$0 is the FP64-equivalent work, $\sim 8$1 is memory traffic, $\sim 8$2 is low-precision tensor throughput, $\sim 8$3 is the compute expansion under emulation, $\sim 8$4 is the bandwidth multiplier, $\sim 8$5 is reconstruction latency per output, and $\sim 8$6 is the number of outputs reconstructed.

For the Bailey FFT, the small inner factor is $\sim 8$7, hence $\sim 8$8 for $\sim 8$9. With recommended $$\eta_{\text{FP64}^{\text{native}}}=\mathrm{OI}_{\text{FFT}}\cdot B_{\mathrm{HBM}}\approx 1.56\,B_{\mathrm{HBM}},$$0 CRT moduli, one has $$\eta_{\text{FP64}^{\text{native}}}=\mathrm{OI}_{\text{FFT}}\cdot B_{\mathrm{HBM}}\approx 1.56\,B_{\mathrm{HBM}},$$1, so $$\eta_{\text{FP64}^{\text{native}}}=\mathrm{OI}_{\text{FFT}}\cdot B_{\mathrm{HBM}}\approx 1.56\,B_{\mathrm{HBM}},$$2. In the companion model, reconstruction latency $$\eta_{\text{FP64}^{\text{native}}}=\mathrm{OI}_{\text{FFT}}\cdot B_{\mathrm{HBM}}\approx 1.56\,B_{\mathrm{HBM}},$$3 had amortized when $$\eta_{\text{FP64}^{\text{native}}}=\mathrm{OI}_{\text{FFT}}\cdot B_{\mathrm{HBM}}\approx 1.56\,B_{\mathrm{HBM}},$$4; here it does not. The reconstruction term becomes dominant. Kulisch arithmetic becomes relevant precisely because the Bailey factorization places the FFT in a regime where the cost of reconstruction no longer disappears into the tensor-core work.

This suggests that the importance of Kulisch arithmetic is contingent rather than universal. The method is most compelling when the decomposition generates many outputs with small inner dimensions and the hardware asymmetry between weak FP64 and strong INT32 is sufficiently pronounced.

3. From recursive Garner to tensor-core Garner

The original reconstruction is recursive mixed-radix Garner: $$\eta_{\text{FP64}^{\text{native}}}=\mathrm{OI}_{\text{FFT}}\cdot B_{\mathrm{HBM}}\approx 1.56\,B_{\mathrm{HBM}},$$5 with digits

$$\eta_{\text{FP64}^{\text{native}}}=\mathrm{OI}_{\text{FFT}}\cdot B_{\mathrm{HBM}}\approx 1.56\,B_{\mathrm{HBM}},$$6

This reconstruction costs $$\eta_{\text{FP64}^{\text{native}}}=\mathrm{OI}_{\text{FFT}}\cdot B_{\mathrm{HBM}}\approx 1.56\,B_{\mathrm{HBM}},$$7 small-modulus integer work per output. The stated estimate is

$$\eta_{\text{FP64}^{\text{native}}}=\mathrm{OI}_{\text{FFT}}\cdot B_{\mathrm{HBM}}\approx 1.56\,B_{\mathrm{HBM}},$$8

which for $$\eta_{\text{FP64}^{\text{native}}}=\mathrm{OI}_{\text{FFT}}\cdot B_{\mathrm{HBM}}\approx 1.56\,B_{\mathrm{HBM}},$$9 is about $12.5$0 INT32 ops per output. Over all Bailey outputs for a $12.5$1 FFT, the paper gives

$12.5$2

and with effective Barrett-reduced INT32 throughput of $12.5$3 TOPS on B300, the projected wall time is near $12.5$4 ms. This is identified as the first $12.5$5-roof bottleneck (Matsuoka, 28 May 2026).

To reduce that cost, the reconstruction is reformulated through the forward CRT identity

$12.5$6

The coefficients $12.5$7 are fixed for a given modulus set, so they can be byte-sliced. Let

$12.5$8

For $12.5$9, $8$0, and one writes

$8$1

so that

$8$2

This yields a two-phase decomposition. Phase A computes the $8$3 as small inner products on FP8/INT8 tensor cores. Phase B reduces

$8$4

to the final reconstructed value. Phase A is cheap; for a $8$5 FFT on B300 it is estimated at $8$6–$8$7 ms. However, without Kulisch arithmetic, Phase B remains a per-output FP64-weighted sum

$8$8

which moves the bottleneck onto the weak FP64 vector pipe. The paper estimates

$8$9

on B300. In other words, tensor-core Garner removes the $10\times$0 modular bottleneck but, absent Kulisch arithmetic, simply replaces it with an FP64 summation bottleneck.

4. Kulisch Phase B as complete arithmetic

The Kulisch reformulation observes that Phase B has an exact fixed-point structure: $10\times$1 Because $10\times$2, each term is an integer shifted to a known bit position. The accumulation can therefore be performed as an integer bit deposit plus carry propagation, with no floating-point operations in the inner loop (Matsuoka, 28 May 2026).

The dynamic-range estimate is

$10\times$3

A $10\times$4-bit accumulator, implemented as $10\times$5INT32, is therefore stated to be sufficient with margin. The compile-time constants are

$10\times$6

$\sim 1.3$23

The stated average cost is about 4 INT32 operations per slice: one shift-left, one optional shift-right, one addcarry, and one add-with-carry into the next word. With $10\times$7, this gives

$10\times$8

In the paper’s terminology, this is complete arithmetic in the Kulisch sense: the accumulation is exact and the only rounding occurs at the end. The argument is three-part. First, each term $10\times$9 is exactly representable as a bit pattern in the fixed-point accumulator. Second, the additions are integer additions and hence exact. Third, only the final integer-to-FP64 conversion introduces rounding. The paper then adds that the final conversion contributes error $$C=\mathrm{Garner}(\tilde C^{(1)},\dots,\tilde C^{(r)};\{m_i\})/(s_A s_B),$$0 ulp and the subsequent FP64 divide by $$C=\mathrm{Garner}(\tilde C^{(1)},\dots,\tilde C^{(r)};\{m_i\})/(s_A s_B),$$1 contributes another $$C=\mathrm{Garner}(\tilde C^{(1)},\dots,\tilde C^{(r)};\{m_i\})/(s_A s_B),$$2 ulp, for a total claimed output error of

$$C=\mathrm{Garner}(\tilde C^{(1)},\dots,\tilde C^{(r)};\{m_i\})/(s_A s_B),$$3

A central conceptual distinction follows. Ordinary FP64 accumulation rounds at every add; Kulisch complete arithmetic aligns all terms in a wide fixed-point register and rounds only once at readout. In this setting, “complete arithmetic” therefore refers to exact reduction, not to elimination of all floating-point stages from the full FFT pipeline.

5. Performance model, parity floors, and hardware mapping

The Kulisch reformulation is designed specifically for the INT32 SIMT pipeline. That architectural placement is decisive: the weak resource on B300 is the FP64 vector pipe, whereas INT32 vector throughput remains high. The total 3-D FFT memory traffic is given as

$$C=\mathrm{Garner}(\tilde C^{(1)},\dots,\tilde C^{(r)};\{m_i\})/(s_A s_B),$$4

so for $$C=\mathrm{Garner}(\tilde C^{(1)},\dots,\tilde C^{(r)};\{m_i\})/(s_A s_B),$$5,

$$C=\mathrm{Garner}(\tilde C^{(1)},\dots,\tilde C^{(r)};\{m_i\})/(s_A s_B),$$6

At $$C=\mathrm{Garner}(\tilde C^{(1)},\dots,\tilde C^{(r)};\{m_i\})/(s_A s_B),$$7 TB/s, the memory roof is

$$C=\mathrm{Garner}(\tilde C^{(1)},\dots,\tilde C^{(r)};\{m_i\})/(s_A s_B),$$8

The native operational intensity is

$$C=\mathrm{Garner}(\tilde C^{(1)},\dots,\tilde C^{(r)};\{m_i\})/(s_A s_B),$$9

For $$\tilde A=\mathrm{round}(s_A\cdot A),\qquad \tilde B=\mathrm{round}(s_B\cdot B).$$0, $$\tilde A=\mathrm{round}(s_A\cdot A),\qquad \tilde B=\mathrm{round}(s_B\cdot B).$$1, so

$$\tilde A=\mathrm{round}(s_A\cdot A),\qquad \tilde B=\mathrm{round}(s_B\cdot B).$$2

The key parity floors are summarized below.

Quantity	Expression	Value stated for $$\tilde A=\mathrm{round}(s_A\cdot A),\qquad \tilde B=\mathrm{round}(s_B\cdot B).$$3, $$\tilde A=\mathrm{round}(s_A\cdot A),\qquad \tilde B=\mathrm{round}(s_B\cdot B).$$4, $$\tilde A=\mathrm{round}(s_A\cdot A),\qquad \tilde B=\mathrm{round}(s_B\cdot B).$$5
Native FP64 floor	$$\tilde A=\mathrm{round}(s_A\cdot A),\qquad \tilde B=\mathrm{round}(s_B\cdot B).$$6	$$\tilde A=\mathrm{round}(s_A\cdot A),\qquad \tilde B=\mathrm{round}(s_B\cdot B).$$7 TFLOPS at $$\tilde A=\mathrm{round}(s_A\cdot A),\qquad \tilde B=\mathrm{round}(s_B\cdot B).$$8 TB/s
Naive Ozaki Phase-B FP64 floor	$$\tilde A=\mathrm{round}(s_A\cdot A),\qquad \tilde B=\mathrm{round}(s_B\cdot B).$$9	Derived for FP64 weighted Phase B
Kulisch INT32 sub-floor	$N=pq$0	$N=pq$1 TOPS at $N=pq$2 TB/s
Kulisch INT32 sub-floor with overlap	$N=pq$3	Factor-of-two relaxation
FP8 Phase-A floor	$N=pq$4	$N=pq$5 PFLOPS at $N=pq$6 TB/s

The Kulisch sub-floor is obtained by equating total INT32 work $N=pq$7 to the memory-roof time $N=pq$8, yielding

$N=pq$9

At $\sim 1.3$00 TB/s, this becomes $\sim 1.3$01 TOPS. The paper states that B300 provides $\sim 1.3$02 TOPS INT32 vector throughput, exceeding the requirement by about $\sim 1.3$03, while its FP8 capability of $\sim 1.3$04 PFLOPS also exceeds the $\sim 1.3$05 PFLOPS Phase-A floor (Matsuoka, 28 May 2026).

The hardware interpretation is binary. A GPU reaches memory-roof FFT parity either by satisfying the native FP64 floor or, if it fails that condition, by satisfying both the FP8 Phase-A floor and the Kulisch INT32 sub-floor. The source names this second alternative the Ozaki-Bailey-Kulisch path.

The projected timing changes accordingly. For $\sim 1.3$06 FP64 3-D FFT on B300, the paper gives: memory roof $\sim 1.3$07 ms, native FP64 path $\sim 1.3$08 ms, recursive Garner Ozaki-Bailey $\sim 1.3$09 ms, tensor-core Garner plus naive FP64 sum $\sim 1.3$10 ms, and tensor-core Garner plus Kulisch Phase B projected at $\sim 1.3$11 ms. The Kulisch Phase B component itself is estimated as

$\sim 1.3$12

or about $\sim 1.3$13–$\sim 1.3$14 ms at 50–70% realistic INT32 efficiency.

6. Numerical guarantees, limitations, and broader significance

The strongest numerical claim is local to Phase B. The accumulation in the wide fixed-point register is exact; only the final integer-to-FP64 conversion and the final FP64 division by $\sim 1.3$15 round, producing a claimed total error of at most a couple of ulps. The paper characterizes this as “full fp64 accuracy,” “full fp64 precision,” and “FP64-equivalent” output, and describes it as “strictly better accuracy than naive fp64 sum.” At the same time, the source does not provide a full formal proof of bitwise equivalence to the exact-real result rounded once through all FFT stages (Matsuoka, 28 May 2026).

Several caveats delimit the claim. The wall-time result is projected rather than measured. No production library yet implements the Kulisch Phase B kernel. Cumulative twiddle roundoff at production sizes still requires empirical measurement. Ill-conditioned inputs requiring adaptive scaling are not analyzed in detail. The explicit pseudocode and range discussion are given for non-negative $\sim 1.3$16 slices, and the paper does not present a complete signed fixed-point derivation beyond that setting. A further misconception is therefore that the source provides a full signed fixed-point CRT reconstruction theory for every FFT case; it does not.

The implementation stance is also specific. No special Kulisch hardware is required. The method is presented as implementable entirely in software on standard GPU integer units using a thread-local $\sim 1.3$17INT32 accumulator, fixed compile-time shifts and word indices, add-with-carry idioms, and, where necessary, warp-level carry propagation. Register pressure is stated as feasible on Blackwell without spilling; $\sim 1.3$18 slices may stream from shared memory; and if multiple threads contribute to one output, their 160-bit accumulators can be combined by a 5-word warp.reduce.add, with roughly 25 ops/warp/output included in the 50–70% efficiency estimate.

The broader significance is programmatic but constrained. The paper generalizes the same pattern to reductions of the form

$\sim 1.3$19

where $\sim 1.3$20 are small bounded integers and $\sim 1.3$21 are fixed positional weights, and derives the generalized sub-floor

$\sim 1.3$22

It suggests potential relevance for Ozaki-II SpMV, ReproBLAS-like reductions, and some stencils. A plausible implication is that Kulisch complete arithmetic is being repositioned as a software primitive for accelerator-era mixed-precision emulation rather than as a requirement for specialized arithmetic hardware. Yet the fully developed use case remains narrow: the Phase B reconstruction primitive in Ozaki-Bailey FFT.

In that sense, the designation “Kulisch escape route” is exact. The method neither changes the FFT decomposition nor supersedes native FP64 where native FP64 already satisfies the bandwidth-parity condition. It is a software reformulation of a single bottlenecked reduction, chosen because the byte-sliced forward-CRT terms align naturally to fixed bit positions and because the surviving INT32 SIMT pipeline on B300 can absorb the work that the FP64 vector pipe cannot.