Fast Fourier transforms at awkward lengths

Every FFT library is fast at 1024. The interesting question is what happens at 1009.

algo-fft — the FFT library in my algo-* suite — is at v0.7.3. This post is a full measurement of where it stands: against FFTW3, against the other Go FFT libraries, at powers of two and at the lengths that are actually awkward, with and without its assembly.

How this was measured

The numbers come from go-fft-bench, a harness that runs five FFT libraries over the same lengths in one process:

• algo-fft — this library, both precisions.
• FFTW3, through cgo bindings. The reference implementation of “fast”, and the baseline everything else is divided by. It is C with generated codelets and decades of tuning; the distance to it is the honest measure.
• gonum (dsp/fourier) — the pure-Go library most Go projects actually reach for.
• go-dsp — pure Go, arbitrary lengths via Bluestein, no reusable plan.
• takatoh/fft — a compact radix-2 implementation, power-of-two only.

Five methodological points, because a benchmark you cannot criticise is a benchmark you cannot trust:

1. Steady state, plan reused. Every library that has a plan object gets one, built outside the timing loop. go-dsp and takatoh have no such API and allocate on every call — that is their design, and it shows up as a fixed tax rather than as a measurement artefact.
2. Libraries are interleaved per length, not run in five long blocks. This is a 15 W laptop part that throttles under sustained load, so absolute numbers drift downward over a sweep. Interleaving puts that drift on every library at a given length equally, which is what makes the ratios trustworthy even where the absolute figures are a laptop’s.
3. The runner refuses to start on a busy machine. It reads the load average and declines above a threshold derived from the core count. Nothing downstream can distinguish a throttled result from a slow library, so the check has to happen before the run rather than in a caveat afterwards.
4. Correctness is a precondition, not a footnote. Every non-power-of-two length is cross-checked against gonum’s transform of the same input before it is allowed into the benchmark. A fast wrong answer is not a result, and this is exactly where FFT libraries quietly fail: hand a radix-2 implementation n = 1000 and it will cheerfully return something.
5. The y-axis on the throughput charts is pseudo-MFLOPS, the 5·n·log₂n / t convention from the FFTW benchmarks. Nobody executes exactly that many operations. Its only job is to divide out the O(n log n) growth so that n = 8 and n = 32768 can share an axis: a flat curve means “scaling as the algorithm says it should”, and a dip means a length fell off a fast path.

Everything below is one machine — an i7-1255U, Go 1.26.1, Linux, GOAMD64=v3 — and one machine is not a benchmark suite. Treat the ratios as real and the absolute figures as indicative.

Powers of two

Two groups, and the gap between them is the whole story of writing numerical code in Go. algo-fft and FFTW3 are in one; gonum, go-dsp and takatoh are in the other, a factor of ten to fifty below. As a ratio against FFTW3 across all thirteen lengths, geometric mean: algo-fft 1.37×, gonum 0.09×, go-dsp 0.02×, takatoh 0.03×.

algo-fft is at or above FFTW3 across the whole power-of-two ladder. It wins by a lot below n = 64 — 5.3× at n = 8, where the transform is 6 ns and FFTW3’s 34 ns is mostly the cost of dispatching through a plan at all — and then tracks it within about ±10% from 128 upward, dipping to 0.91× at n = 1024 and rising to 1.27× at 8192.

Against the rest of the Go field it was never close. At n = 1024, one forward complex128 transform costs 1.35 µs in algo-fft, 22.0 µs in gonum, 45.0 µs in go-dsp and 76.9 µs in takatoh, against FFTW3’s 1.24 µs.

What the assembly buys

algo-fft ships 229 registered codelets — kernels bound to one specific transform length, so that a length like 1024 runs straight into a routine that knows it is 1024 and skips every dispatch decision. They are distributed unevenly on purpose:

Tier	Build	complex64	complex128
generic Go	every build	32	32
AVX2	amd64, default	33	27
SSE3 / SSE2	amd64, default	22	22
NEON	arm64, default	22	22
AVX-512	amd64, default	11	6

The SSE tier is split by precision: SSE3 carries complex64 and SSE2 carries complex128, because the complex64 multiply idiom wants ADDSUBPS, which is SSE3. Machines with SSE2 but not SSE3 stopped being made around 2005, so they keep the generic Go path, deliberately. The NEON ladder is complete but untuned — every power of two from 4 to 32768 has an arm64 codelet, verified for correctness under emulation, with priorities mirrored from the amd64 ordering rather than measured, because QEMU timings are not performance data. And the AVX-512 row reaches only n = 128 in double precision, because the machine everything here was measured on has no AVX-512 at all; that tier lives on a Xeon I can borrow at weekends, which is not a tuning loop.

Two things that are not switches, since both get assumed. There is no asm build tag — SIMD is part of the default build, selected at runtime by CPU feature detection. And GOAMD64 gates none of it: that variable controls what Go’s own compiler emits for the portable code, not which hand-written kernel the registry hands you. The single opt-out is -tags purego, which removes every assembly package and every SIMD codelet registrar and leaves the portable engine.

Which makes “how much of this is algorithms and how much is assembly?” a question you can answer by building twice.

The assembly is worth 2.95× on geometric mean for complex128 and 4.78× for complex64, peaking at 8.15× at n = 256. Single precision gains more because a 256-bit register holds eight complex64 lanes against four complex128, and because the complex64 codelet ladder is the more complete of the two. Arbitrary lengths give 2.45×, with a few where it buys nothing at all — more on those below.

Strip the assembly out entirely and algo-fft is still 4.97× ahead of gonum, 19.0× ahead of go-dsp and 14.3× ahead of takatoh, and it does not lose at any length — the narrowest margin is 3.35× at n = 16384. That surprised me, because the plan machinery and the strategy ladder are the parts of this library I have spent the least time defending.

A note on what “AVX2” is worth

One caveat on those figures, because the label is easy to trust too far. algo-fft’s power-of-two workhorses used to be per-size AVX2 radix-4 codelets: VEX-encoded, dispatched only on machines reporting AVX2 — and loading their operands with VMOVSD, 64 bits at a time, with every butterfly done in X registers. They were scalar kernels in VEX clothing, with the three-operand encoding and none of the width. A grep settles it: YMM arithmetic across the whole family came to 0–2 instructions per file, against 100–560 for the kernels that are genuinely vectorized.

Replacing them with one size-generic kernel per precision — covering 16 to 65536, further than any per-size codelet reached — was worth 2.2× to 4.0×. The surprise inside that work was where the time had been going: at n = 16384 the bit-reversal permutation, which does no arithmetic whatsoever, was a third of the kernel’s runtime. So “has an AVX2 codelet” and “is fast” turn out to be almost independent claims.

Awkward lengths

A transform length is easy or hard depending on one thing: whether it factors into radices you have butterflies for. algo-fft has three answers, and choosing between them is most of the planner’s job.

• The mixed-radix engine factors n into 2, 3, 5, 7 and 11 and runs it exactly. Nothing is padded, nothing is wasted.
• Rader’s algorithm turns a prime-length transform into a cyclic convolution of length p−1. A good trade only when p−1 itself factors nicely, which for any given prime is a coin flip you can check at plan time.
• Bluestein’s algorithm handles anything at all by re-expressing the DFT as a convolution and padding to a convenient length of at least 2n−1. It is the universal fallback, and it costs roughly three transforms of a larger length to do the work of one.

Neither Rader nor mixed-radix is a strategy you can ask for by name — both are decisions the planner makes underneath the Bluestein plumbing, each gated by a measured win condition rather than by the mathematical test alone. You can ask a plan what it decided, which is how the benchmark labels its own results:

    97  Bluestein  rader          1153  Bluestein  rader
   257  Bluestein  rader          1009  Bluestein  bluestein
   448  DIT        dit_fallback   2205  Stockham   stockham
  1152  Stockham   stockham      44100  Stockham   stockham

This chart is the most useful thing in the post, and it is backwards from the naive expectation. The class where algo-fft does best is the primes — 1.05× of FFTW3 for primes with smooth p−1, which is to say it wins that class outright — against 0.38× for 5-smooth lengths and 0.48× for the 7/11-smooth ones. The composite lengths ought to be the comfortable case and the primes the hard one. They are not.

Per length, Rader reaches 1.58× of FFTW3 at n = 12289, 1.33× at 641 and 1.06× at 4001 — for prime-length transforms, in Go. The composite side runs 0.2× to 0.6×, and the worst single result in the sweep is n = 2205 at 0.16×.

Three lengths — 257, 1009 and 2003 — are where the SIMD build has no advantage at all over the pure-Go one (0.94×, 1.02×, 0.95×). All three route through Bluestein or Rader, whose padded sub-transforms are powers of two and therefore run in the fastest assembly in the library. Since the totals do not benefit, the time cannot be in the sub-transforms; it is in the chirp construction and the pointwise products around them, which are still ordinary scalar Go. That is a concrete lead rather than a mystery, and it is the next thing on that path.

The one length I want to call out by name is 44100, the canonical audio sample rate. It is 2²·3²·5²·7², it factors perfectly, and it takes 745 µs — 0.25× of FFTW3, and 3.1× faster than gonum. Not a victory, but the routing is finally right: a gate in the planner used to exclude lengths whose power-of-two part is only 2 or 4 from the mixed-radix engine, sending 44100 to Bluestein instead. That gate had a measurement behind it, and the measurement was real; what it had actually measured was a defect in the mixed-radix driver, which dispatched every sub-transform larger than a single point to a codelet — including the size-2 and size-4 recursion leaves that have a perfectly good Go butterfly. A profile at n = 1000 found 1.9% of runtime in the codelet assembly and roughly 40% in dispatch overhead. A threshold fitted to a bug is a good way to lose a length for a year.

Single versus double precision

Halving the element width does not halve the time. Single precision buys 1.53× on geometric mean, unevenly: 1.07× at n = 8, where plan dispatch dominates and the element type barely matters, rising to 2.17× at 32768, with the middle of the range around 1.4–1.5×. Worth knowing before you convert a pipeline to float32 on the assumption of a free doubling — on this hardware, for this workload, it is about a half.

The complex64 curve sits above FFTW3’s double-precision line at every length measured, which is the one place in this post where the comparison is not apples to apples: FFTW can be built for single precision too, and this harness does not link that build. Read the green curve against the blue one, and the orange one only as a sense of scale.

One thing worth knowing if you write numerical Go, single-precision or not: Go does not implement scalar complex64 * complex64 in single precision. It widens all four float32 components to float64, multiplies in double precision, and rounds the two results back — twelve instructions where the complex128 multiply needs six. Any FFT stage written as plain Go therefore costs more in the narrower type, which is a genuinely surprising thing to discover in your own library and an easy one to carry for years without noticing. algo-fft now routes those products through a helper that stays in single precision; the price is 3–9% more relative L2 error against a float64 reference, all of it still at float32 epsilon. You are no longer getting accidental double-precision intermediates, and no longer paying for them.

What is not done

An honest list, because a benchmark post that only reports wins is an advertisement:

• Composite lengths are the weak class. 0.38× and 0.48× of FFTW3 against Rader’s 1.05×. The mixed-radix engine is the part of the library furthest from the reference.
• Bluestein’s scalar overhead is unaddressed, which is why the assembly buys nothing at 257, 1009 and 2003.
• arm64 is unmeasured. The NEON codelet ladder is complete and correct, and its selection priorities are guesses copied from the amd64 ordering. It needs an afternoon on real hardware.
• AVX-512 is barely used, and reaches only n = 128 in double precision — the development machine has none, so that tier is tuned on borrowed weekend access.
• 134 functions that use YMM or ZMM registers never execute VZEROUPPER, which leaves the processor in a state where subsequent legacy-SSE instructions pay a transition penalty. The sweep that found them exists; acting on it does not.
• Plan-level overhead now dominates small transforms — a consequence of the kernels getting faster, and a category of work I had not needed to think about before.
• WebAssembly gets no SIMD, a toolchain limit rather than a decision: Go’s SIMD support is still experimental and amd64-first.

Reproducing it

Everything here is a checkout and two commands:

git clone https://github.com/CWBudde/go-fft-bench
cd go-fft-bench
just sweep    # both builds, waits for an idle machine
just plot     # renders every chart in this post

The charts are drawn by a small tool in that repository using matplotlib-go, my Go port of the Matplotlib model, which makes this post a fairly demanding test of that library as well as of this one. Log-base-2 axes, grouped bars with categorical labels and multi-series legends all came out of it without a fight. It found one bug in return, and it was mine: bars on a logarithmic axis run down to zero, and log(0) is not a coordinate.

algo-fft itself is at github.com/CWBudde/algo-fft, and there is a WebAssembly demo if you would rather watch a transform than read about one.