Roadmap

go-fft/fft is a pure-Go (cgo-free) FFT library — the numpy.fft / scipy.fft equivalent for Go, with no dependency on the native FFTW3 C library.

Motivation

  • Ruby has no cgo-free FFT. Every Ruby option binds the FFTW3 C library, so an embedded-Ruby runtime that wants to stay CGO_ENABLED=0 cannot use them. go-fft is the FFT backend for go-embedded-ruby (shipped in Phase 5, see below).
  • gonum/dsp/fourier is pure Go but its optimized assembly is amd64-only. go-fft targets a portable core first, on all six 64-bit Go targets (amd64, arm64, riscv64, loong64, ppc64le, s390x), then SIMD kernels generated by go-asmgen wherever they are bit-identical and measured faster. Today that is amd64 (see Phase 4).

Conventions (org HARD RULES)

  • Pure Go, CGO_ENABLED=0, no vendoring, build from source.
  • 100% statement coverage, enforced in CI (gate >= 100.0); the test suite must also exit 0, not merely report coverage.
  • Validated on all six 64-bit targets in CI (amd64/arm64 native, the rest under qemu-user).
  • English only.

Phase 0 — correct pure-Go FFT — DONE

  • Complex FFT/IFFT over []complex128, any length.
    • Radix-2 iterative Cooley–Tukey for power-of-two lengths (bit-reversal permutation + in-place butterfly stages).
    • Bluestein’s chirp-z algorithm for arbitrary N, reducing the DFT to a power-of-two convolution.
    • Unnormalized forward convention; IFFT divides by N so IFFT(FFT(x)) ≈ x.
  • FFTReal(x []float64) []complex128 — full complex spectrum.
  • Edge cases: empty input, length 1; inputs are never mutated.
  • Butterfly/bit-reversal loops isolated in internal/kernels so SIMD variants can drop in behind the same signatures.
  • Tests: round-trip (power-of-two and non-power-of-two), cross-check against a naive O(N²) reference DFT, impulse → flat spectrum, constant → DC spike, sinusoid → single bin spike, linearity, no-mutation.

Phase 1 — real-optimized RFFT — DONE

  • RFFT(x []float64) []complex128 returns only the non-redundant N/2+1 bins (numpy.fft.rfft convention); IRFFT(spectrum, n) []float64 reconstructs the length-n real signal (numpy.fft.irfft), normalized by n so IRFFT(RFFT(x), len(x)) ≈ x.
  • Even lengths use the half-length complex packing trick: pack the real signal into a length-N/2 complex array, take one size-N/2 FFT, then untangle — ~2× throughput over FFTReal. Odd lengths fall back to the full complex transform and slice the kept bins.
  • First benchmark suite added (FFT/RFFT/FFTReal/IRFFT vs the naive O(N²) DFT baseline) to anchor the Phase 4 SIMD work.
  • Tests: cross-check against the full complex spectrum and the naive DFT, conjugate symmetry of the dropped half, real DC/Nyquist bins, an analytic cosine spike, round-trip for even and odd N, explicit Hermitian-mirror reconstruction, edge cases, and no-mutation.

Phase 2 — multi-dimensional transforms — DONE

  • FFTN(data, shape)/IFFTN(data, shape) over a flat row-major (C-order) []complex128 paired with an explicit []int shape: the separable N-dimensional DFT, applying the 1-D FFT along each axis in turn (numpy.fft.fftn/ifftn). Forward is unnormalized; the inverse divides by the product of the transformed axis lengths so IFFTN(FFTN(x)) ≈ x. A mismatched or non-positive shape panics, mirroring numpy.
  • FFT2(data, [2]int)/IFFT2(data, [2]int) — the 2-D specialization (numpy.fft.fft2/ifft2), a thin wrapper over FFTN.
  • Real 2-D variants RFFT2(data, [2]int)/IRFFT2(spectrum, [2]int) for image-style workloads (numpy.fft.rfft2/irfft2): the real transform runs along the last axis (keeping cols/2+1 bins per row), the full complex transform down the columns; IRFFT2 inverts and is normalized by rows*cols. The caller passes the target shape explicitly.
  • None of these mutate the caller’s input.
  • Tests: a separable naive-DFT N-D oracle, round-trip, FFT2≡FFTN, explicit row-then-column separability, impulse→flat, constant→DC spike, linearity, RFFT2 vs the kept bins of FFT2, RFFT2/IRFFT2 round-trip, zero-padded short spectra, empty/length-1 shapes, shape-validation panics, and no-mutation.

Phase 3 — windowing and spectral helpers — DONE

  • Symmetric window functions returning []float64 of a given length: Hann, Hamming, Blackman (matching numpy.hanning/hamming/blackman), Bartlett (numpy.bartlett), and the 4-term BlackmanHarris (scipy.signal.windows.blackmanharris). n <= 0 returns an empty slice, n == 1 returns [1].
  • FFTFreq(n, d) and RFFTFreq(n, d) bin-frequency helpers matching numpy.fft.fftfreq/rfftfreq exactly — including the [0, …, (n-1)//2, -(n//2), …, -1] / (d·n) ordering for fftfreq and the n//2+1 non-negative bins for rfftfreq.
  • PSD(x, d) one-sided power spectral density (periodogram, density scaling with the negative-half folding, matching scipy.signal.periodogram), and Spectrogram(x, segment, overlap, window, d) which windows successive overlapping segments and returns a slice of one-sided PSD frames (STFT convention: trailing partial segments dropped).
  • Tests assert the EXACT numpy/scipy window and frequency values computed by hand, window symmetry, a Parseval check tying PSD back to the signal mean square (even and odd N, exercising the Nyquist-doubling branch), spectrogram framing against an independent per-segment PSD, edge cases, validation panics, and no-mutation.

Phase 4 — SIMD kernels via go-asmgen — DONE (routed on amd64)

The hot kernel chosen first is the pointwise (Hadamard) complex multiply, a[i] *= b[i] over []complex128, which drives Bluestein’s spectral product. The scalar CMulScalar is the portable correctness oracle; CMul is the stable dispatch seam the FFT core calls; per-arch cmulSIMD wrappers (over go-asmgen assembly) are validated bit-for-bit against the oracle.

Generated kernels — amd64 (SSE2), arm64 (NEON), s390x (vector), riscv64 (RVV), all bit-identical

  • internal/kernels/asmgen/amd64/gen.go drives go-asmgen to emit cmul_amd64.s: an SSE2 loop, one complex128 per iteration, computing (ar·br − ai·bi, ar·bi + ai·br) with MULPD/SHUFPD/XORPD/ADDPD and a low-lane sign mask built in-register (no read-only data table). It is bit-for-bit identical to the scalar oracle because every product and the final add are separately rounded exactly as the scalar code rounds them — no fused multiply-add (GOAMD64=v1 does not fuse).
  • internal/kernels/asmgen/arm64/gen.go drives go-asmgen to emit cmul_arm64.s: a 2-wide NEON loop processing two complex128 per iteration. VLD2 deinterleaves a pair into a vector of reals and a vector of imaginaries, the arithmetic runs two-lane, and VST2 re-interleaves on store; an odd final element is handled by a scalar tail. It is bit-for-bit identical to the scalar oracle — see the fusion note below.
  • internal/kernels/asmgen/s390x/gen.go drives go-asmgen to emit cmul_s390x.s: a 2-wide z/Architecture vector-facility loop processing two complex128 per iteration. VL loads the pairs, VMRHG/VMRLG deinterleave them into a vector of reals and a vector of imaginaries (big-endian lane order), VFMDB/VFMSDB/VFMADB run the fused arithmetic two-lane, and VMRHG/VMRLG + VST re-interleave on store; an odd final element is handled by a scalar tail. It is bit-for-bit identical to the scalar oracle on the project’s only big-endian target — see the fusion note below.
  • internal/kernels/asmgen/riscv64/gen.go drives go-asmgen to emit cmul_riscv64.s: a strip-mined RVV 1.0 loop. Each iteration asks the hardware for vl = min(remaining, VLMAX) via VSETVLI (E64, M1), then VLSEG2E64V deinterleaves a run of {re, im} pairs into a vector of reals and a vector of imaginaries, the fused arithmetic runs vl-wide (VFMUL.VV/VFNMSAC.VV/VFMACC.VV), and VSSEG2E64V re-interleaves on store. Strip-mining folds the tail (the last iteration just gets a smaller vl), so there is no separate scalar tail. It is bit-for-bit identical to the scalar oracle, hardware-validated on cfarm95 — see the fusion and soundness notes below.
  • The generators live in a separate Go module (internal/kernels/asmgen/go.mod, requiring go-asmgen), so the FFT library itself keeps zero third-party dependencies. go generate (or the per-arch CI job) regenerates the committed .s and CI fails if it is stale.

Bit-identity decides which arches ship a kernel (correctness over coverage)

The contract is bit-for-bit equality with the scalar oracle, and which floating-point form the gc compiler gives the oracle differs by arch:

  • amd64: GOAMD64=v1 does not fuse, so the oracle is separately-rounded multiplies + a rounded add/sub. SSE2 MULPD/ADDPD reproduce that exactly → bit-identical.

  • arm64: FMA is baseline, so the gc compiler fuses the oracle itself (FMSUBD for ar·br − ai·bi, FMADDD for ar·bi + ai·br; verified by disassembling CMulScalar with -gcflags=-S). A NEON kernel is bit-identical only if it reproduces that same fusion — which it does: each product is accumulated into a zeroed register with VFMLA (adding +0.0 leaves the product’s rounding unchanged, matching the oracle’s FMULD), and the second term is folded in with VFMLS/VFMLA (the fused subtract/add the oracle performs). The lesson is fusion is matched, not avoided — a non-fused NEON kernel, or one whose fusion form differed, diverges by up to 1 ULP, caught by the random SIMD-vs-scalar test (Go’s arm64 assembler exposes no non-fused vector multiply — there is no VFMUL — so matching the fused form is the only sound option anyway).

  • s390x: FMA is baseline and — uniquely among the four qemu-only targets — the Go assembler exposes 2-wide double-precision vector FMA. The gc compiler fuses the oracle (FMUL then FMSUB/FMADD; verified by disassembly), and the vector-facility kernel reproduces that fusion two-lane (VFMDB for the exactly-rounded product, VFMSDB/VFMADB for the fused subtract/add), with VL/VMRHG/VMRLG/VST deinterleaving and re-interleaving the complex pairs. s390x is the project’s only big-endian target, so its per-arch qemu job proving the kernel bit-identical (correct lane order included) is the highest-value validation in the suite. Ships a kernel. ✅

  • riscv64: FMA is baseline, so the gc compiler fuses the oracle (FMULD then FNMSUBD for ar·br − ai·bi; FMULD then FMADDD for ar·bi + ai·br; verified by disassembling CMulScalar with GOARCH=riscv64 -gcflags=-S). The RVV 1.0 kernel reproduces that fusion vl-wide — VFMUL.VV for each exactly-rounded product, then VFNMSAC.VV (vd = −(vs1·vs2)+vd) and VFMACC.VV (vd = +(vs1·vs2)+vd) for the fused subtract/add — the precise RVV analogue of the arm64 NEON kernel. Crucially, the V extension is OPTIONAL on riscv64, so unlike the other three SIMD arches cmulSIMD detects V at run time (parsing /proc/cpuinfo’s isa string) and calls the RVV kernel only when V is present; on a non-V CPU it falls back to scalar (calling VSETVLI there would SIGILL). The kernel is hardware-validated bit-identical on cfarm95 (GCC Compile Farm: a SpacemiT X60, RVA22 + RVV 1.0), where the SIMD-vs-scalar test runs and passes on structured and random inputs across a length sweep. Under CI’s non-V qemu-riscv64 the scalar path is taken and the SIMD test skips its RVV comparison (t.Skip("no RVV")), so CI stays green and the RVV instructions never execute there. Ships a kernel (hardware-validated). ✅

  • loong64 / ppc64le: keep the validated scalar path, each for a concrete, checked reason — NOT for lack of trying:

    • loong64: the gc compiler fuses the oracle (MULD, then FMSUBD / FMADDD; verified by disassembly), so a bit-identical kernel needs a vector fused multiply-add. As of Go 1.26 the loong64 assembler has vector double add, subtract and multiply (VADDD/VSUBD/VMULD and their 256-bit XV forms) but no vector FMA, which would have to be encoded by hand. (Until Go 1.26 it exposed vector floating point only as unary ops.)
    • ppc64le: the gc compiler fuses the oracle (FMUL, then FMSUB / FMADD), and the Go ppc64 assembler exposes no vector double arithmetic (as of Go 1.26 the VSX surface is loads/stores, logicals, permutes and conversions; XVMADDDP/XVMULDP/XVADDDP are not assemblable), so a bit-identical vector kernel would have to be encoded by hand entirely.

    Per “on n’a pas le droit de se tromper”, a 1-ULP-different kernel is a correctness divergence the project rejects; partial-arch SIMD with identical results beats full-arch SIMD with a divergence.

Measured perf (honest delta)

On Apple-silicon arm64 the Go compiler’s scalar loop reaches ~48 GB/s (L1) and out-runs the hand-written 2-wide NEON kernel at this width (measured ~3× slower; the VLD2/VST2 deinterleave dominates). amd64 SSE2 is similar. So even where a kernel is bit-identical, routing the hot path through it would be a regression. The dispatch (CMul) therefore keeps the scalar implementation on the hot path; each SIMD kernel is retained as a per-arch-validated artifact and the reference for future wider kernels (AVX2/AVX-512, or arm64 SVE / unrolled NEON), where a win is plausible. Benchmarks: BenchmarkCMulScalar* vs BenchmarkCMulSIMD*.

Split CI + coverage policy

CI is split into a pure-Go coverage gate and per-arch execution jobs, mirroring the go-asmgen org:

  • The test gate builds, vets, race-tests, and enforces 100% statement coverage of the portable core. It runs on amd64 (ubuntu/windows) and arm64 (macos); on each host exactly one cmulSIMD wrapper compiles and is exercised by the SIMD-vs-scalar test, so the gate reaches 100% without excluding any reachable line.
  • The generated .s carries no Go statements (it cannot be statement- covered), so it is validated by the per-arch jobs: arch-native (amd64/arm64) regenerates the committed assembly and fails if stale, then vet (asmdecl cross-checks name+offset(FP)), build (cmd/asm encodes the mnemonics), and runs the suite whose TestCMul* assert SIMD == scalar bit-for-bit; arch-qemu runs riscv64/loong64/ppc64le/s390x under qemu-user (s390x is the big-endian guard; riscv64’s RVV .s is regen-checked there but runs scalar, since that qemu CPU lacks V — the RVV bit-identity proof is run on real hardware, cfarm95).

Every line is either statement-covered by the gate or executed by a per-arch job; the gate is never lowered and there is no coverage-gaming knob.

Progress

  • amd64 (SSE2) — bit-identical SIMD kernel, validated by the native job. ✅
  • arm64 (NEON) — bit-identical 2-wide SIMD kernel, validated by the native job (reproduces the oracle’s FMSUBD/FMADDD fusion). ✅
  • s390x (vector facility) — bit-identical 2-wide SIMD kernel (VFMDB/VFMSDB/VFMADB reproducing the oracle’s fused FMSUB/FMADD), validated by the qemu-s390x job on the project’s only big-endian target. ✅
  • riscv64 (RVV 1.0) — bit-identical strip-mined SIMD kernel (VLSEG2E64V + VFMUL.VV/VFNMSAC.VV/VFMACC.VV reproducing the oracle’s fused FNMSUBD/FMADDD), with run-time V detection, hardware-validated on cfarm95 (SpacemiT X60, RVV 1.0); under CI’s non-V qemu it runs scalar and the SIMD test skips. ✅
  • loong64 / ppc64le — validated scalar fallback (qemu jobs); no SIMD kernel, each for a concrete checked reason (no assemblable vector FMA, and on ppc64le no vector double arithmetic at all — see above). 🔒

What the transforms run today (v0.19.0)

The pointwise multiply above was the first kernel; it stayed off the hot path everywhere, because the compiler’s scalar loop measured as fast. The kernels that do run came later, and they are whole passes. Each is bit-identical to the Go code it replaces and tested bit for bit on generic, ±0 and ∞ inputs.

  • amd64, AVX2: the Stockham passes of radix 2/3/4/5/8, plus a final pass vectorized across blocks. Every power of two runs entirely on them, and so does every other 5-smooth length; radix-7, 11 and 13 passes stay in Go.
  • amd64, AVX-512: the same passes four points per register, used for powers of two of 256 points or more. Elsewhere, mixing widths in a transform measured slower.
  • amd64, AVX2: the real-FFT untangle and its inverse.
  • amd64, AVX2, N-D: an axis other than the last runs as batched passes over strips of 16 or 32 neighbouring lines, in place (v0.8.0). The radix order puts 3 and 5 first, then 4 before 8, then 7/11/13 (measured over every ordering of 12 lengths).
  • amd64, AVX2 without AVX-512: a radix-16 pass for 128, 256 and 1024 points, and 2048 on Intel (v0.10.0). From 4096 points its 16 streams fall into one L1 set and it loses, so larger sizes keep radix 8/4.
  • amd64, AVX-512 on Intel, from 65536 points: a blocked schedule (v0.12.0). The first passes run over the whole array (the first two fused into one sweep), the rest in 8192-point groups between two L2-sized buffers, and the final pass writes each group straight into the output. Beyond L2 every pass runs at copy speed; this crosses the L2 boundary twice instead of six or seven times. AMD keeps the plain schedule (its L3 holds the arrays).
  • amd64, AVX2 on AMD, powers of two: the passes keep their data split (four real parts, then four imaginary parts, per 64 bytes) between radix-4/8 AVX2 passes, removing every shuffle and sign flip (v0.14.0).
  • amd64, Intel AVX2, composites up to 16384: a radix-12 prime-factor pass (Good–Thomas 3·4) and radix-16 tails after the odd passes (v0.15.0).
  • amd64, AMD, composites: radix-10/12/15/20 prime-factor passes with stride limits (v0.18.0).
  • amd64, Intel AVX-512, powers of two from 128: 512-bit split kernels (eight points per register) from 256 points, radix 16 at 128 (v0.19.0).
  • amd64 without AVX2: SSE2 butterfly stage kernels for the iterative pow2 kernel. It also had AVX2 stages until v0.1.8; since v0.1.3 an AVX2 machine routes every power of two to the Stockham engine, so they ran only in tests and were removed.
  • arm64, NEON: the Stockham passes of radix 2/3/4/5/8 at any stride, and the radix-2 to 5 final passes over an even number of blocks (v0.7.0). gc fuses multiply-adds on arm64, so each kernel reproduces, pass by pass, which product the Go code rounds and which it fuses; results equal the Go code’s bit for bit except the sign of a NaN produced from infinite inputs. Radix 7, the general radix and odd final passes stay in Go. Powers of two from 16 points keep their data block-split between passes (real and imaginary parts apart within each 32 bytes), which removes the per-pass deinterleave (v0.11.0), and the columns of N-D plans run as NEON batched passes. Since v0.15.0 composites run odd radices first and split too (radix-3/5 split kernels), and the final passes take any block count. The real-FFT untangle runs as a NEON kernel, and small N-D plans run each axis as one pass sequence over all lines (v0.17.0).
  • float32, AVX2 and NEON: the single-precision Stockham passes (v0.13.0), the real-FFT untangle and retangle, and the N-D column passes (v0.16.0) have their own kernels, bit-identical to the float32 Go code (except NaN bits from infinite inputs).
  • ppc64le, riscv64, loong64, s390x: the Go passes, compiled to scalar instructions. gc does not vectorize, and fuses multiply-adds on all four.

Every routing and radix rule behind this was measured on real hardware; the dated rounds in BENCHMARKS.md give the numbers, including what was tried and dropped:

  • a first NEON Stockham pass (WORD-encoded, interleaved twiddles): 1.13–1.24× on Neoverse-N1 but slower on large strides; the shipped one stores twiddles pre-split and wins at every stride there (Round 18);
  • the six-step algorithm;
  • a fan-out threshold that depends on the machine.

Remaining

  • Small powers of two on amd64 (v0.27.0, see Performance): complex 256 1.10× FFTW on Zen 3 and 1.052× on Cascade Lake, complex 1,024 1.13× and 1.19×, and on Zen 3 RFFT 256 and 4,096 (1.051×, 1.09×). FFTW does complex 256 and 1,024 and RFFT 4,096 in two passes with 32- and 64-point codelets. Round 38 built a bit-identical radix-32 pass for amd64; it lost to split radix 8 wherever it was timed and was not shipped. Neoverse-N1 is faster than FFTW on all 24 rows since v0.24.0 (re-measured at v0.26.1: RFFT 256 0.998×, RFFT 1,080 0.994×).
  • Zen 3: complex and RFFT 1,080 (1.08×, 1.06×). Round 38 re-measured its RFFT rows against FFTW at v0.26.1, after v0.25.0’s untangle tail jump: they moved within that measurement’s spread.
  • Cascade Lake: the composites 1,000–1,920 (1.10–1.13×), complex and RFFT 65,536 (1.08×, 1.09×). Every 2-D row is faster than FFTW since v0.27.0 (Round 39). The 2-D 1024² regression of v0.19.0 is fixed in v0.22.0 (Round 31).
  • Haswell (Intel AVX2-only) was measured in Rounds 29 and 37 (v0.21.0, v0.24.0), but it has no FFTW, so its ratio to FFTW is unknown; v0.25.0’s untangle tail jump, which addresses the ~1.5% Round 37 measured on RFFT 256 there, has not been re-measured on it. Apple M4 has not been measured with the current engines.
  • float32: AVX-512 kernels, and a comparison with FFTW’s single-precision library (not built on the hosts).
  • loong64. Never measured: its only host stayed at load 150 or unreachable.
  • An RVV-capable CI runner, so the riscv64 bit-identity proof also runs in CI. Today it runs scalar there, and the proof is on real hardware (cfarm95).

Phase 5 — Ruby binding — DONE

  • Exposed through go-embedded-ruby as an FFT module, so Ruby gets a cgo-free FFT matching the numpy.fft / scipy.fft surface: fft/ifft/rfft/irfft, the N-D and 2-D transforms (fftn/ifftn/fft2/ifft2), the bin-frequency helpers (fftfreq/rfftfreq), the window functions (hann/hamming/blackman/blackman_harris/bartlett), and the spectral helpers (psd/spectrogram). Spectra are returned as Ruby Complex values. The binding required a Ruby Complex type and the :: scope operator, both added to go-embedded-ruby; it is differentially validated against the numpy.fft conventions and held to the same 100%-coverage gate.

Phase 6 — numpy.fft / scipy.fft surface — DONE (v0.3.0–v0.6.0)

  • numpy’s optional arguments (v0.4.0): n, norm and axes through ...With(Options) variants whose zero value is numpy’s defaults; HFFT/ IHFFT, RFFTN/IRFFTN, FFTShift/IFFTShift (N-D too), NextFastLen, batched 1-D transforms, and a Norm argument on every plan.
  • DCT and DST of types I–IV, forward, inverse and N-D, as scipy.fft defines them (v0.5.0), with reusable plans.
  • Single precision (v0.6.0): Plan32, RealPlan32, FFT32, RFFT32; then N-D and 2-D (FFTN32, RFFTN32, FFT2_32, RFFT2_32, PlanN32, RealPlan2_32), the Options variants (FFT32With, …, HFFT32) and the DCT/DST (DCT32, DSTN32, DCTPlan32, …), checked against scipy’s float32 results.
  • Not planned: numpy’s s (per-axis crop/pad) for N-D, orthogonalize=False for the DCT/DST.