Skip to content

go-fft 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 (planned binding, Phase 4).
  • gonum/dsp/fourier is pure Go but its optimized assembly is amd64-only. go-fft targets a portable scalar core first, then SIMD kernels generated by go-asmgen across all six 64-bit Go targets (amd64, arm64, riscv64, loong64, ppc64le, s390x).

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 — IN PROGRESS

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 Go loong64 assembler exposes LSX vector floating-point only as unary ops (VFSQRTD, VFRINTD, VFRECIPD, …) — there is no vector float add/multiply/FMA to build a complex product from, so a bit-identical vector kernel is not expressible.
  • ppc64le: the Go ppc64 assembler exposes no vector double arithmetic (the VSX surface is loads/stores, logicals, permutes, conversions — XVMADDDP/XVMULDP/XVADDDP are not assemblable), so likewise no bit-identical vector kernel is expressible.

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 double arithmetic in the Go assembler — see above). 🔒

Remaining

  • loong64/ppc64le SIMD await Go-assembler support for vector double arithmetic (or a different, still bit-identical formulation). riscv64 RVV is done — hardware-validated on cfarm95; a future improvement would be an RVV-capable CI runner so the bit-identity proof also runs in CI (today it runs scalar there and the proof is on real hardware).
  • Wider kernels (amd64 AVX2/AVX-512, arm64 SVE / unrolled NEON, multiple complex per iteration) and a measurement to decide if/when to route CMul through SIMD on the hot path.
  • A bit-identical vectorized butterfly (Radix2 inner loop) once the pointwise multiply is solid across arches.

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.