In the original MFFT post the polynomials were written down to single binary digits. The point of that choice was not aesthetics: once every coefficient is \(0\) or \(1\), a pointwise matrix product

\[ C_{ij} = \sum_k A_{ik}\,B_{kj} \]

stops being a general multiply–add. It is the count of positions where both factors are one:

\[ C_{ij} = \#{k : A_{ik}=B_{kj}=1} = \operatorname{popcount}!\bigl(\mathrm{row}_i(A)\ \mathrm{AND}\ \mathrm{col}_j(B)\bigr). \]

No scalar multiply in the leaf. Only bitwise AND and a population count along the contraction index \(k\).

mfft-bench now has that leaf on CPU and GPU. This note is the report: how we got there, what the complexity claim actually means, and the numbers on an RTX 5070 Ti.

---

Why this sits next to MFFT

Classical integer multiplication (Karatsuba, Toom–Cook, Schönhage–Strassen) reduces a wide product to many digit products. MFFT does the same for matrices of wide integers: each digit plane is an \(n\times n\) integer matrix, and each plane product is a GEMM.

If the digit base is \(2\) — one bit per limb — those GEMMs are Boolean. The cost model people care about then changes:

leaf	rough cost for width-\(m\) entries
schoolbook digits + general GEMM	\(\sim m^2\,n^\omega\)
binary planes + AND/popcount GEMM	\(\sim m\,n^3/w\)

where \(w\) is the machine word width (64 on CPU, 32 in the GPU pack). Scaling in bit width \(m\) is linear instead of quadratic. That is the advantage the base-2 story was always aiming at. Matrix-Strassen still helps the \(n\) exponent, but it does not change the \(m^2\) digit factor.

The trap: Karatsuba breaks \({0,1}\)

A first attempt to measure boolpack under --only karatsuba with LIMB_BITS=1 looked disappointing and inconsistent. The reason is algebraic, not implementation:

Karatsuba forms sums \(A_0+A_1\). On binary limbs that sum lives in \({0,1,2}\). Treating a \(2\) as a set bit is wrong; checking is_01 and falling back to schoolbook is correct — and then most of the timed work is no longer Boolean.

Fair measurement uses schoolbook on the limb axis (--only limbplane) so every plane stays in \({0,1}\).

CPU: limbplane + boolpack

With make LIMB_BITS=1 and 64-bit entries (\(L=64\) planes), the CPU leaf that packs \(k\) into uint64_t and runs popcountll wins cleanly:

\(n\)	ikj (s)	strassen (s)	boolpack (s)	vs ikj
64	0.26	0.42	0.097	2.7×
128	2.04	3.29	0.61	3.3×
256	15.5	23.1	3.48	4.5×

The gap widens with \(n\). Matrix-Strassen is a poor fit for pure \(0/1\) planes here: recursion builds denser intermediate panels while the Boolean leaf is already doing word-parallel work.

SIMD helps further. On hosts with AVX-512 VPOPCNTDQ the inner loop processes eight columns per step (_mm512_popcnt_epi64 on A_row AND B_cols). Scalar POPCNT remains the AVX2 fallback.

GPU: ballot pack + popc + tiles

The CUDA side lives in cuda/boolpack_gpu.cuh and shows up as four rows in gemm_bench:

row	timed work
boolpack-gpu	scalar pack + GEMM
boolpack-ballot	__ballot_sync pack + GEMM
boolpack-gemm	GEMM only (pack outside the timer)
boolpack-tiled	shared-memory tiled GEMM

Packing one row of bits: each warp lane reads one \(k\), then __ballot_sync(0xffffffff, pred) writes a 32-bit word. The GEMM is

acc += __popc(A_row_words[w] & B_col_words[w]);

over those words. The tiled variant stages blocks of words in shared memory.

On an NVIDIA GeForce RTX 5070 Ti (sm_120), \(n=1024\), best of three:

method	ms	TFLOP/s	vs sgemm
cublas-sgemm	0.230	9.4	1.0×
int8-dp4a	0.054	39.6	4.3×
boolpack-tiled	0.034	63.6	6.8×
boolpack-gemm	0.096	22.4	2.4×
boolpack-ballot	0.120	17.8	1.9×
boolpack-gpu	0.295	7.3	0.8×

All four boolpack rows are exact against a host reference on \(n\le 512\).

Two practical notes:

1. Packing is real work. Scalar pack + GEMM loses to sgemm; ballot pack recovers most of it; measuring GEMM alone (or tiling) is what beats the baselines.
2. Compare apples to apples. These rows multiply random 0-1 matrices, not the float embeddings used by cublas-sgemm in the same table. The comparison answers “how fast is exact binary plane GEMM on this GPU?”, which is the leaf cost MFFT pays once limbs are bits.

At \(n=256\) launch overhead dominates every method; the interesting regime starts around \(n=1024\).

What this does not claim

• It does not make full fp32/fp64 MFFT faster than cuBLAS on dense floats. The float embedding still expands into many non-binary limbs; int8 schoolbook and Ozaki slices remain the practical exact paths there.
• It does not replace matrix-Strassen for general integer panels. Once values leave \({0,1}\), boolpack must fall back.
• TFLOP/s on the boolpack rows is a normalized \(2n^3/\mathrm{time}\) figure for fair table layout, not IEEE floating-point ops.

How to reproduce

# CPU binary planes
make clean && make WITH_OPENMP=1 LIMB_BITS=1
./mfft-bench --n 256 --bits 64 --reps 2 --no-naive --only limbplane

# GPU Boolean GEMM (rows appear in the main table)
make -C cuda
./cuda/gemm_bench --n 1024 --reps 3 --check

Code: github.com/hadilq/mfft-bench. Earlier measurement context: MFFT, measured.

Takeaway

The original post’s base-2 model is not only a proof device. On hardware that can pack bits and popcount quickly, the pointwise step of a binary-limb matmul is a Boolean GEMM, and that leaf can be faster than sgemm while staying bit-exact. The open engineering problem is no longer “is AND+popcount real?” — it is fusing that leaf into the limb methods (and MFFT) so that every plane that is still \(0/1\) takes the fast path, and every plane that is not falls back cleanly.

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},
}

@misc{hadilq2026mfftboolpack,
    author = {{Hadi Lashkari Ghouchani}},
    note = {Published electronically; source at \url{https://github.com/hadilq/mfft-boolpack}},
    title = {MFFT, boolpack: Boolean GEMM on the limb axis: when AND+popcount beats sgemm},
    year = {2026},
}