Most numerical code is trusted the way most bridges are trusted: it hasn’t fallen down yet. You run it, the answer looks plausible, the plot has the right shape, and you move on. For a great deal of engineering that is entirely rational — the cost of being wrong is small and the feedback is fast. But there is a narrower class of code where „it looks right“ is not good enough: a reference implementation, the thing other people will check their own work against. A reference is only worth the name if its correctness is earned and auditable, not merely un-refuted. This is a short account of what that costs, told through one of the most-implemented algorithms in computing — the discrete Fourier transform — and a small project that set out to implement it to that standard in four languages at once.
Why the DFT, and why four times
The discrete Fourier transform turns a sequence of samples into a sum of pure frequencies. It is the workhorse behind spectrum analysers, image compression, the numerical solution of partial differential equations, radio astronomy, and the detection of gravitational waves. Its fast form, the FFT, is arguably the most consequential numerical algorithm of the twentieth century. It is also, precisely because it is everywhere, a place where errors hide comfortably: everyone assumes the DFT is correct, so nobody checks the DFT.
Implementing it once proves you can get an answer. Implementing it independently in several languages and demanding that they agree proves something stronger — that the answer is a property of the mathematics, not an accident of one compiler, one library, or one programmer’s convention. The fourier project implements the DFT/FFT as a hand-written kernel in Fortran, C++, Pascal, and Rust, and then holds those four renderings to a single question: do they produce the same bits?
The trap of the plausible average
The first and most important lesson is about how you compare. The instinct is to reduce: compute the total energy of the spectrum, or its norm, and check that the two implementations agree to some tolerance. This is exactly the comparison that fails silently.
A Fourier spectrum’s total energy is invariant under a sign flip or a complex-conjugation of individual bins (Parseval’s theorem sees only magnitudes). So a backend with a sign-convention error in its twiddle factors — a genuinely wrong transform — will still report the correct total energy. Across the sibling project kh-sim, four independent hand-rolled FFTs were each carrying the same „conjugate-of-standard“ convention drift, and every scalar check passed. The bug was caught, four times over, only by an element-wise comparison: the full complex field, bin by bin, real and imaginary parts separately. The rule that falls out of this is blunt and load-bearing: reductions are lossy by design; compare the whole field, and make the scalar shortcut something you have to ask for explicitly.
Tolerance is not machine epsilon
The second lesson is about the number you compare to. Floating-point arithmetic is not exact, so „identical“ needs a definition. Two definitions are wrong in opposite directions. Demanding equality to machine epsilon (~2×10⁻¹⁶) is too strict — an FFT accumulates rounding over its butterflies, and the honest error grows roughly as ε·√N. Picking a tolerance loose enough to pass is worse: it launders the very errors the test exists to catch.
The project uses a single, defensible rule — Higham’s backward-error bound, in the form tol = factor · scale · ulp(precision) · √(n_ops), with the scale floored per observable class. The per-backend factor lives with the backend, never in the gate, so the gate stays platform-neutral. Crucially, when a test cannot meet its gate, the answer is never to widen the number. The failing test stays failing, with the real cause recorded, and the decision is escalated. A green light bought by relaxing a gate past the source’s own reported precision is not a pass; it is a regression wearing a pass’s clothing.
Bit-identity, and the one honest exception
Held to element-wise comparison at a principled tolerance, the four languages converge on something quietly remarkable: for the same inputs at the same precision, the hand-written Fortran, C++ and Pascal kernels produce bit-for-bit identical IEEE-754 results. Not „close“ — identical. Three people, three languages, three compilers, one bit pattern. That is what it means for a result to be a property of the mathematics rather than of the tooling.
The exception is instructive rather than embarrassing. The Rust backend diverges from the other three by about 4 units-in-the-last-place — not from any error, but because it reaches for the language’s intrinsic transcendental functions where the others call the same underlying libm. The three libm-backed backends are bit-identical to each other; Rust is bit-identical-within-locked-tolerance. Documenting that boundary — libm-versus-intrinsic — is itself part of the canonical claim. A reference implementation that couldn’t explain its own last four bits would not deserve the adjective.
The case for keeping the hard way
There is an obvious objection: why hand-write an FFT at all, when battle-tested libraries exist? The project’s answer is a deliberate inversion. The library-backed implementations were correct from the first run and stayed correct. Every hand-rolled one had the convention bug — and that is the argument for them, not against. The bug was findable precisely because the transform was four lines of visible arithmetic rather than an opaque call. A reference exists to be read and checked; the languages that keep the mathematics on the surface — Fortran and Pascal especially, with their long lineage in Numerical Recipes, whose 1986 Pascal edition is a direct historical anchor for this work — earn their place by staying reviewable. The library versions are the trusted units; the hand-written ones are the ones that teach.
The philosophical residue
Strip away the specifics and what remains is a small epistemology of computation. We tend to treat a number produced by a program as a fact. It is not; it is a claim, and like any claim it is only as good as the procedure that could have shown it false. The discipline sketched here — compare everything element-wise, gate against the precision the source actually reports, never relax a gate to pass, seek agreement across independent renderings, and keep the mathematics visible enough to audit — is really just the experimental method applied to arithmetic. Reproducibility across independent implementations is to numerical software what reproducibility across independent laboratories is to physics: not a nicety, but the thing that turns a result into knowledge.
Four languages, one transform, identical bits. The bits are the easy part to state and the hard part to earn. What the earning buys is the right to say canonical and mean it.
—
Source: github.com/lege-artis/fourier (v0.3.1). Code Apache-2.0; this article CC-BY-SA-4.0. Part of the lege-artis canonical-reference series.
