A year ago I wrote about Matrix Fast Fourier Transform (MFFT): treat matrix entries as polynomials over a digit base, evaluate them at a family of matrix-valued roots of unity that are really just signed permutations, multiply pointwise, and map back. The hope was to shrink the digit-side cost of matmul from something like \(m^2 O(n^\omega)\) toward \(\tilde O(m),O(n^\omega)\).

Theory is cheap. Code is less so. So I built mfft-bench: a C and CUDA benchmark that implements the post’s transform, puts it next to the methods people actually use, and measures both throughput and relative error against a bit-exact product.

This post is the report. It worth mentioning that this post is written by Grok 4.5. Also the code was initiated by Opus 5, but concluded by Grok 4.5.

---

What we implemented

Three tracks share one codebase:

1. Exact integer track — \(n\times n\) matrices of \(B\)-bit integers. Every method must produce the bit-identical product. This is the setting MFFT is defined in.
2. ML track — \(n\times n\) fp32 (and fp64) matrices, with accuracy reported against the bit-exact product of a fixed-point embedding.
3. GPU track — the same ideas on int8 tensor-core / __dp4a paths, with cuBLAS, Strassen, Ozaki residual slices, and limb-exact / limb-MFFT rows.

The roots of unity from the original post are not complex numbers. Each power \(I_s^e\) is a signed permutation: exactly one \(\pm 1\) per row and column. In the code they are compiled once into a flat op list (u, v, sign, mode) and then applied as fused butterflies — no general multiplies for twiddles, only adds, subtracts, and sign flips. That part of the post checks out; ./mfft-bench --test-roots verifies the recursive construction of \(H_{s,k}\) against the matrix powers.

One correction relative to a naive reading of the post: if you pack one bit per coefficient, intermediate growth forces an enormous transform length. The practical fix is the usual Schönhage–Strassen move — pack \(S\) limbs into each coefficient so the ring length \(K\) and the number of blocks \(NB\) stay manageable. The GPU path does the same with 7-bit signed limbs that fit in int8 accumulators.

The uncomfortable measurement

At machine-learning widths the embedding lands at a handful of limbs:

embedding	typical \(L\) (limbs)
bf16 / int8 / int4	1
fp32	4–8
fp64	8–12

MFFT only overtakes Karatsuba on product count near \(L \approx 512\). Below that it does more pointwise matmuls than schoolbook. A CPU profile at \(n=128\) makes the bottleneck obvious:

fp32→mfft-rec:
  pack        ~0.1%
  build_ops   ~0%
  transform   ~1–3%
  pointwise   ~95%
  fold        ~2%

Almost all the time is the \(NB\cdot K^2\) small integer GEMMs after the transform. Fusing the butterflies harder will not change the ranking when the transform is three percent of the wall clock.

Representative CPU ML numbers (\(n=128\), fp32 inputs, error vs bit-exact product):

method	relative cost	rel error
sgemm (packed)	1×	\(\sim 3\times 10^{-7}\)
fp32→karatsuba / toom3	\(\sim 30{-}40\times\)	exact (\(\sim 3\times 10^{-8}\))
fp32→limb-plane	\(\sim 50\times\)	exact
fp32→mfft	\(\sim 100\times+\)	exact

So on CPU: more accurate than sgemm, and one to two orders of magnitude slower. That is not a tuning failure. Exactness needs \(\gtrsim 24 + \log_2 n\) bits per entry; even a two-limb Karatsuba product is several integer GEMMs against one fp32 GEMM, and integer GEMM is not free.

The GPU flips the arithmetic

A consumer GPU inverts two CPU facts at once: int8 is several times faster than fp32, and fp64 is roughly an order of magnitude slower. Limb decomposition wants exactly that mix.

On an RTX 5070 Ti at \(n=4096\) (promoted or genuine fp64; best-of-reps):

method	GEMMs	vs baseline	note
cublas-sgemm	1	1×	\(\sim 13,\mathrm{ms}\)
cublas-dgemm	1	1× for fp64	\(\sim 170,\mathrm{ms}\)
limb-fp32-exact	\(\sim 56)	\(\sim 9\times\) sgemm	bit-exact fp32 product
limb-fp64-exact	\(\sim 56{-}64\) (skip-zero)	\(\sim 0.7{-}1.7\times\) dgemm	bit-exact
limb-fp32-faithful	9	\(\sim 1.8\times\) sgemm	correctly-rounded quality
limb-fp64-faithful	25	beats dgemm (\(\sim 0.35\times\))	\(\sim 10^{-10}\)–\(10^{-12}\) error
limb-mfft-fp32	128	much slower	exact, too many products at \(L\sim 8\)
ozaki-i8-s2…s7	4–49	competitive	residual slices, not bit-exact

Two takeaways:

1. The limb embedding is useful on GPUs even when MFFT is not. Exact fp64 via int8 planes is in the same ballpark as cublasDgemm; faithful rounding (keep only the product bits that can affect a correctly rounded binary result) can be faster than native fp64 while staying far more accurate than fp32.
2. MFFT on the GPU track still loses at ML limb counts, for the same structural reason as on the CPU: product count, not butterfly speed.

Where MFFT is the right tool

The post’s complexity story assumed \(m\) large enough that the transform’s \(\tilde O(m)\) behaviour matters. That regime is real; it is just not training a transformer in fp16.

Domain	Why MFFT fits
FHE / lattice cryptography (BFV, BGV, CKKS)	Ciphertexts already live in \(\mathbb{Z}_q[y]/(y^K+1)\) — the same negacyclic ring the implementation uses. Coefficients are hundreds to thousands of bits. Matrix-valued coefficients are a natural way to batch linear algebra over ciphertexts.
Multiprecision integer matmul	Bit-widths past a few kilobits: schoolbook is \(O(L^2)\) limb products; MFFT is \(\tilde O(L)\).
Zero-knowledge proofs of linear algebra	Wide integers, exact arithmetic, often already in a cyclotomic ring.

There the signed-permutation roots are a genuine advantage: no modular multiplies for twiddles, only the fused add/sub/sign butterflies. Recursive SSA on the ring dimension (the implementation’s mfft-rec) is the profitable recursion; a second MFFT level on pointwise entries only helps when those entries are themselves multi-limb, which float embeddings are not.

Methods that sat next to MFFT

A fair bench needs company:

• Schoolbook limb planes — the baseline exact product; still the winner at small \(L\).
• Karatsuba and Toom-3 — fewer products than schoolbook once \(L\) grows; Toom-3 is hybridised with Karatsuba so it never loses on product count at power-of-two widths.
• Ozaki residual slices — split each float into a few scaled integer planes; fewer GEMMs when approximate accuracy is enough.
• Strassen — recursive \(O(n^{\log_2 7})\) on the matrix dimension, not the limb dimension; orthogonal to MFFT.
• Faithful high-limb truncation — drop low-significance limbs so the product is not bit-exact but is correctly rounded (or close). This is the row that makes GPU “exact-ish” fp64 attractive.

Skip-zero plane detection and exponent-band buckets help when the data is narrow (magnitudes in \([0.5,1)\)); under uniform \(U(-1,1)\) every band is live and the bookkeeping costs more than it saves.

What I would tell my past self

1. Implement the transform. The signed-permutation roots and the fused op list are the elegant part of the original post, and they work.
2. Measure product count at the \(L\) you actually have. For fp32/fp64 embeddings, \(L\) is single digits. MFFT is the wrong asymptotic there; Karatsuba or schoolbook is enough.
3. Separate the embedding from the transform. Mapping floats onto a shared fixed-point grid and multiplying exactly is independently useful — especially on GPUs where int8 is cheap and fp64 is expensive.
4. Point at the domains that match the math. FHE, multiprecision, and verifiable linear algebra sit at the limb widths where \(\tilde O(L)\) convolution earns its keep.

MFFT is not a better sgemm. It is a fast multiprecision convolution with a particularly friendly root-of-unity representation, and it becomes practical exactly when the coefficients get wide. That is a narrower claim than the one I started with, and a more useful one.

Code

git clone https://github.com/hadilq/mfft-bench
cd mfft-bench
make && make check
./mfft-bench --ml --n 128 --reps 2 --no-naive --profile
make -C cuda
./cuda/gemm_bench --n 4096 --reps 3 --fp64

Details, tables, and caveats live in the repository README.

Cite

If you found this work useful, please consider citing the original note and this measurement report:

@misc{hadilq2024MatMul,
    author = {{Hadi Lashkari Ghouchani}},
    note = {Published electronically at \url{https://hadilq.com/posts/matrix-fast-fourier-transform/}},
    title = {Matrix Fast Fourier transform(MFFT)},
    year = {2024},
}

@misc{hadilq2026mfftbench,
    author = {{Hadi Lashkari Ghouchani}},
    note = {Published electronically; source at \url{https://github.com/hadilq/mfft-bench}},
    title = {MFFT, measured: when the transform helps and when it does not},
    year = {2026},
}