---
title: Kulisch Fixed-Point Complete Arithmetic
url: https://www.emergentmind.com/topics/kulisch-fixed-point-complete-arithmetic
type: topic
---

# Kulisch Fixed-Point Complete Arithmetic

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 [2606.23698].

## 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 [2606.23698].

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 \(\to\) length-\(p\) column FFTs \(\to\) twiddle multiplication \(\to\) transpose \(\to\) length-\(q\) column FFTs \(\to\) 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 \(N=1024\), the natural factorization is \(p=q=32\) [2606.23698].

This factorization creates the key difficulty. In the Tensor–Memory Equilibrium model, the emulated execution time is
\[
T_{\text{emu}}=\max\!\left(\frac{\alpha W}{P_{\text{low}}},\frac{\beta Q}{B_{\mathrm{HBM}}}\right)+\gamma n_{\text{out}},
\]
where \(W\) is the FP64-equivalent work, \(Q\) is memory traffic, \(P_{\text{low}}\) is low-precision tensor throughput, \(\alpha\) is the compute expansion under emulation, \(\beta\) is the bandwidth multiplier, \(\gamma\) is reconstruction latency per output, and \(n_{\text{out}}\) is the number of outputs reconstructed.

For the Bailey FFT, the small inner factor is \(k\approx \sqrt{N}\), hence \(k=32\) for \(N=1024\). With recommended \(r=12\) CRT moduli, one has \(r^2=144\), so \(k=32\ll r^2\). In the companion model, reconstruction latency \(\gamma\) had amortized when \(k\gg r^2\); 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:
\[
v=v_1+m_1(v_2+m_2(v_3+\cdots +m_{r-1}v_r)),
\qquad v_k\in [0,m_k),
\]
with digits
\[
v_k=\left(v'_k-\sum_{j<k} v_j\prod_{i<j} m_i\right)\cdot \left(\prod_{i<k} m_i\right)^{-1}\pmod{m_k}.
\]
This reconstruction costs \(O(r^2)\) small-modulus integer work per output. The stated estimate is
\[
\text{ops/output}\approx 2.5r^2,
\]
which for \(r=12\) is about \(360\) INT32 ops per output. Over all Bailey outputs for a \(1024^3\) FFT, the paper gives
\[
W_{\text{Garner}}\approx 6.8\times 10^{12}\ \text{INT32 ops},
\]
and with effective Barrett-reduced INT32 throughput of \(\sim 25\) TOPS on B300, the projected wall time is near \(260\) ms. This is identified as the first \(\gamma\)-roof bottleneck [2606.23698].

To reduce that cost, the reconstruction is reformulated through the forward CRT identity
\[
C \equiv \sum_{k=1}^r v'_k\cdot u_k \pmod M,\qquad
u_k=\frac{M}{m_k}\cdot \left(\frac{M}{m_k}\right)^{-1}_{m_k},\qquad
M=\prod_i m_i.
\]
The coefficients \(u_k\) are fixed for a given modulus set, so they can be byte-sliced. Let
\[
S=\left\lceil \frac{7r}{8}\right\rceil.
\]
For \(r=12\), \(S=11\), and one writes
\[
u_k=\sum_s 256^s u_k^{(s)},\qquad u_k^{(s)}\in [0,256),
\]
so that
\[
\sum_k v'_k u_k=\sum_{s=0}^{S-1}256^s \underbrace{\sum_{k=1}^r v'_k u_k^{(s)}}_{=:P_s}.
\]

This yields a two-phase decomposition. Phase A computes the \(P_s\) as small inner products on FP8/INT8 tensor cores. Phase B reduces
\[
\sum_s P_s 256^s
\]
to the final reconstructed value. Phase A is cheap; for a \(1024^3\) FFT on B300 it is estimated at \(\sim 0.5\)–\(1\) ms. However, without Kulisch arithmetic, Phase B remains a per-output FP64-weighted sum
\[
\sum_s P_s\cdot c_s,\qquad c_s=\frac{256^s}{s_A s_B},
\]
which moves the bottleneck onto the weak FP64 vector pipe. The paper estimates
\[
T_{\text{Phase B}}=\frac{3\cdot 1024^2\cdot 6\cdot 1024\cdot S}{\eta_{\text{FP64-vec}}}\approx 163\ \text{ms}
\]
on B300. In other words, tensor-core Garner removes the \(O(r^2)\) 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:
\[
y=(s_A s_B)^{-1}\sum_{s=0}^{S-1} P_s\cdot 256^s,\qquad P_s\in \mathbb{Z},\ \log_2 P_s\lesssim 20,\ S=11.
\]
Because \(256^s=2^{8s}\), 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 [2606.23698].

The dynamic-range estimate is
\[
\log_2\!\left(S\cdot \max P_s\cdot 256^{S-1}\right)\approx 104\ \text{bits}.
\]
A \(160\)-bit accumulator, implemented as \(5\times\)INT32, is therefore stated to be sufficient with margin. The compile-time constants are
\[
\mathit{word}[s]=\left\lfloor \frac{8s}{32}\right\rfloor,\qquad
\mathit{shift}[s]=(8s)\bmod 32.
\]

```text
Algorithm: Kulisch Phase B for one Bailey-FFT output
W[0..4] <- 0      // five INT32 registers = 160-bit accumulator
for s = 0,1,...,10:
    w <- floor(8s / 32)
    b <- (8s) mod 32
    lo <- P_s << b
    hi <- P_s >> (32 - b) if b > 0 else 0
    (W[w], c) <- addcarry(W[w], lo)
    W[w+1] <- W[w+1] + hi + c
return fp64(W) / (s_A s_B)
```

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 \(S=11\), this gives
\[
cS\approx 4\cdot 11=44\ \text{INT32 ops/output}.
\]

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 \(P_s256^s\) 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 \(\le 1\) ulp and the subsequent FP64 divide by \(s_A s_B\) contributes another \(\le 1\) ulp, for a total claimed output error of
\[
\le 2\ \text{ulps}.
\]

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
\[
Q=96N^3\ \text{bytes},
\]
so for \(N=1024\),
\[
Q\approx 103\ \text{GB}.
\]
At \(B_{\mathrm{HBM}}=8\) TB/s, the memory roof is
\[
T_{\text{mem}}=\frac{Q}{B_{\mathrm{HBM}}}\approx 12.9\ \text{ms}.
\]

The native operational intensity is
\[
\mathrm{OI}_{\text{FFT}}=\frac{W_{\text{FFT}}}{Q_{\text{FFT}}}
=\frac{15\log_2 N}{96}.
\]
For \(N=1024\), \(\log_2 N=10\), so
\[
\mathrm{OI}_{\text{FFT}}=\frac{150}{96}=1.5625.
\]

The key parity floors are summarized below.

| Quantity | Expression | Value stated for \(N=1024\), \(S=11\), \(c=4\) |
|---|---|---|
| Native FP64 floor | \(\eta_{\text{FP64}^{\text{native}}}\approx 1.56\,B_{\mathrm{HBM}}\) | \(12.5\) TFLOPS at \(8\) TB/s |
| Naive Ozaki Phase-B FP64 floor | \(\eta_{\text{FP64}^{\text{naive-OBFFT}}}\approx 2.06\,B_{\mathrm{HBM}}\) | Derived for FP64 weighted Phase B |
| Kulisch INT32 sub-floor | \(\eta_{\text{INT32}^{\text{Kulisch}}}\approx 8.25\,B_{\mathrm{HBM}}\) | \(66\) TOPS at \(8\) TB/s |
| Kulisch INT32 sub-floor with overlap | \(\eta_{\text{INT32}^{\text{Kulisch, overlap}}}\approx 4.125\,B_{\mathrm{HBM}}\) | Factor-of-two relaxation |
| FP8 Phase-A floor | \(\eta_{\mathrm{FP8}^{\text{Phase A}}}\approx 170\,B_{\mathrm{HBM}}\) | \(1.36\) PFLOPS at \(8\) TB/s |

The Kulisch sub-floor is obtained by equating total INT32 work \(18N^3\cdot cS\) to the memory-roof time \(96N^3/B_{\mathrm{HBM}}\), yielding
\[
\eta_{\text{INT32}^{\text{Kulisch}}}=\frac{3cS}{16}B_{\mathrm{HBM}}\approx 8.25\,B_{\mathrm{HBM}}.
\]
At \(8\) TB/s, this becomes \(66\) TOPS. The paper states that B300 provides \(\sim 75\) TOPS INT32 vector throughput, exceeding the requirement by about \(14\%\), while its FP8 capability of \(5\) PFLOPS also exceeds the \(1.36\) PFLOPS Phase-A floor [2606.23698].

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 \(1024^3\) FP64 3-D FFT on B300, the paper gives: memory roof \(12.9\) ms, native FP64 path \(\sim 124\) ms, recursive Garner Ozaki-Bailey \(\sim 260\) ms, tensor-core Garner plus naive FP64 sum \(\sim 163\) ms, and tensor-core Garner plus Kulisch Phase B projected at \(\sim 18\) ms. The Kulisch Phase B component itself is estimated as
\[
T_{\text{Kulisch}}^{\text{B300}}\approx \frac{18N^3\cdot 44}{\eta_{\text{INT32}}}\approx 13\ \text{ms at peak},
\]
or about \(18\)–\(25\) 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 \(s_A s_B\) 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 [2606.23698].

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 \(P_s\) 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 \(5\times\)INT32 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; \(P_s\) 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
\[
y_i=\sum_{j=1}^{W} P_{ij}w_j,
\]
where \(P_{ij}\) are small bounded integers and \(w_j\) are fixed positional weights, and derives the generalized sub-floor
\[
\eta_{\text{INT32}^{\text{Kulisch}}}=c\cdot \mathrm{OI}_{\text{red}}\cdot B_{\mathrm{HBM}}.
\]
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.

Source: https://www.emergentmind.com/topics/kulisch-fixed-point-complete-arithmetic