Take One, Measure Six: What Go's SIMD Experiment Gave Mere

In September the Go team wrote up their SIMD experiment. Go 1.26 had shipped SIMD for amd64 behind GOEXPERIMENT=simd, with the API shaped by the architecture. Go 1.27 turned it into a platform-independent simd package that also covers arm64 NEON and WebAssembly. It is still an experiment, and the post is a good read on how they reasoned about it.

Mere is the small language I write for myself. It has an interpreter and four compiled backends – C, LLVM IR, WebAssembly, and RISC-V machine code – and it already had SIMD: a compiler pass that lets clang vectorize ordinary loops, and three 128-bit types, f64x2, f32x4 and u8x16 (the earlier note has the numbers). So the question the Go post raised was not “should Mere have SIMD” but “which of Go’s decisions would make Mere’s SIMD better”.

I went through Go’s design one feature at a time and measured each before deciding. One came across. This note is about that one, the six that did not, and what the checking turned up. Every number was measured on arm64 (Apple clang 21, node 24).

The table first

Go’s feature	Mere	What decided it
MulAdd, fused multiply-add	taken: fma, f64x2_fma	a vertex transform in lanes went from 140 ms to 100 ms
LoadVPart / StorePart for loop tails	not taken	measured upper bound on the one program that needs it: 0.16 s to 0.15 s
a mask type per lane width, IfElse	not taken yet	nothing needs a select, and the cost can be removed without a new type
emulation, so the code always runs	not a rule	Mere computes in software only what it can reproduce exactly – so far, fma on Wasm
vector width removed from the type	not taken	every Mere target is 128 bits wide
feature variants (planned for 1.28)	not taken	a runtime-detection problem an ahead-of-time build does not have
reductions, missing in Go 1.27	already there	reduce_add, any_true, first_true

My first recommendation was wrong

Reading the post, I would have put partial loads first. Go’s LoadFloat32sPart(x[i:]) loads the last, short vector of a loop, fills the rest with zeros, and tells you how many elements it read, so the tail of the loop stays in lanes. I assumed Mere’s byte-lane programs duplicate their loops in scalar code for the last few bytes, and that a partial load would remove that duplication.

I checked before building anything, and the premise fell over in one file. The UTF-8 validator already pads its tail into one zero-filled block in user code. Zero is ASCII, so the padding changes nothing about the answer. Only one program, the literal prefilter in mgrep (a grep written in Mere), does its tail, and lines shorter than 16 bytes, in scalar code.

So I measured an upper bound on it: 81 MB of source-code lines, 1.7 million lines averaging 49 bytes, 14% of the bytes going through the scalar path.

mgrep variant	first byte rare	first byte common
as it is (scalar tail)	0.16 s	0.31 s
tail padded in user code	0.26 s	0.38 s
tail skipped entirely – wrong answers, the ceiling	0.15 s	0.27 s

A built-in partial load could at best reach the last row, and that row flatters it, because it also skips the matches in the tail. So the feature I would have ranked first is worth about 6% to the only program that could use it.

The measurement did show where mgrep actually loses. Mere’s SIMD loads read bytes, and a line is a str, so every line is copied first. That copy costs about 40% of the available gain at ordinary line lengths. It is Mere’s own problem, not one Go’s design speaks to, and it stays open under its own name.

Masks and emulation, in brief

In Mere, comparing two u8x16 gives another u8x16 with 0xFF in the matching lanes, which is the NEON and WebAssembly convention. On RISC-V the comparison produces a bit mask in a mask register, and turning it back into bytes costs two instructions. On WebAssembly, asking for the first matching lane costs two more than it would with a real mask type. Against that, the UTF-8 validator uses a comparison as data: it ANDs the result with 1 to count continuation bytes. A separate mask type would need conversions there, and no program needs a select yet. Pattern-matching the common shapes in the backend would remove the extra instructions without a new type. That waits for a JSON scanner or a base64 codec that needs a select.

Go emulates every SIMD operation where the hardware lacks it, so a Go program using simd runs everywhere. Mere refuses an operation by name on a backend that cannot lower it; for example, the RISC-V backend refuses f64x2. It also has an interpreter that every compiled backend is compared against, which is what Go’s GODEBUG=simd=0 provides as a testing tool. Emulation as a rule is worth less to Mere. Emulation as a case-by-case decision turned out to matter, and the case was the one feature Mere did take.

Fused multiply-add

fma a b c is a * b + c rounded once. Written out, a * b + c rounds twice, once after the product and once after the sum. The two answers can differ in the last bits, and they differ most when the product and c nearly cancel: fma 0.1 10.0 (-1.0) is 5.55e-17, while 0.1 * 10.0 - 1.0 is exactly 0.0.

Mere had deliberately gone without it. The C backend emits #pragma STDC FP_CONTRACT OFF, because clang on arm64 fuses a * b + c into an fma instruction at -O2 on its own. When that happened, a dot product printed different bits under the interpreter and under the optimized binary, and Mere’s rule is that the optimizer may not change the answer. Switching fusion off kept that rule, and cost the one instruction that makes a chain of multiply-adds fast.

The Go post suggested asking a different question: could a program ask for the fused answer by name? That keeps the rule. a * b + c still rounds twice everywhere, and fma rounds once everywhere.

Before writing a line, I measured what it could be worth by turning the pragma on in the C that Mere emits for a vertex transform (20,000 vertices through a 4x4 matrix, 5,000 passes):

kernel	pragma off (as shipped)	pragma on	fused instructions
scalar	0.20 s	0.15 s	0 -> 348
f64x2 lanes	0.14 s	0.10 s	0 -> 72

It is 1.4x on both. The answer changed, as it had to, but the fused scalar program and the fused lane program printed the same bits as each other. fma is deterministic, so a fused answer is a well-defined answer and not a rounding accident. With the builtin written, the lane kernel spelled with f64x2_fma runs in 100 ms against 140 ms through the C backend. clang turns each pair of lane fmas into one fmla.2d.

The backend without an fma instruction

Four of the five backends were straightforward. The interpreter uses OCaml’s Float.fma, C calls fma(3), and LLVM calls llvm.fma. The RISC-V backends have no double-precision unit and refuse the name. WebAssembly was the real decision. It has no fma instruction, and its JavaScript host has no Math.fma.

The two options were to refuse it there, which is what Mere usually does, or to compute it in software, which is what Go does. I chose software for one reason: fma is one of the operations IEEE-754 requires to be correctly rounded, like +, -, *, / and sqrt. There is exactly one right answer and it can be reproduced. A refusal would not have spared anyone the problem. A program that wanted fma would write a * b + c for its WebAssembly build and get a different answer there, which is the very thing the pragma exists to prevent.

My first software fma was about fifteen lines of Mere. It uses a trick due to Boldo and Melquiond: compute the exact product as two doubles, add, and round to odd once in the middle. It matched the hardware on 30,000 cases on all four backends. But it holds only while nothing in the middle overflows or underflows – I had checked it on operands between about 2^-100 and 2^100 and nowhere else – and on WebAssembly it allocated per call and ran out of memory at about two million calls. What shipped instead is written in integers inside the WebAssembly runtime, with the structure of musl’s fma.c:

1. The three significands are normalized to 54 bits.
2. The 108-bit product is built from 32-bit halves, because WebAssembly has no 64-by-64 multiply that returns 128 bits.
3. c is aligned against the product, and anything shifted out becomes a sticky bit.
4. The two are added or subtracted.
5. The result is rounded exactly once: the conversion of the top 63 bits to a double is the rounding. A result that will be subnormal is pre-rounded at the subnormal’s last place, so the final scaling does not round a second time.

Checking it took longer than writing it:

• In C first, against the hardware, on 250 million cases drawn from special values, subnormals, the edges of the exponent range, cancellation, and ties. No mismatches.
• Each branch broken on purpose. Eight of ten mutations produced mismatches. The other two could not be made to change any answer, and I worked out why: it would take a 53-bit by 53-bit product with 64 zero bits between two set ones, and there is no such product. So those two are argued unobservable in writing, not claimed as tested.
• Then the WAT itself. I cut it out of an emitted module and ran it under node against hardware results for 50 million more cases. No mismatches. The node harness was checked too, by breaking the WAT and watching it fail.

One of those checks nearly fooled me. Dropping the sticky bit for a shifted-out c passed 20 million random cases. Random inputs essentially never land on a tie, the one case where that bit decides the rounding. With a family of constructed ties added – an exact odd 54-bit product, which rounds at exactly half, and a tiny c far below it that decides which way – the same mutation failed 10,272 of 20 million. Ties are in the permanent gate now. If that bit is removed from the compiler’s WebAssembly runtime, the gate fails and names the family: Wasm differs from fma(3) in: ties.

The cost on WebAssembly is about 12 ns a call. Eight nested fmas per iteration took 130 ms where the unfused expression took 30 ms. In a loop that already allocates a float per iteration, the difference disappears into the noise.

What the gate found that was not about fma

The permanent gate runs six families of inputs through the interpreter, C at two optimization levels, LLVM and WebAssembly, and compares each against the C library’s fma(3), bit for bit. Before it ever went green it had found two bugs in the LLVM backend that had nothing to do with fma.

x != x said a NaN was not a NaN. The gate folds every NaN result into one value, because IEEE-754 does not specify a NaN’s payload, and it detects them with r != r. The LLVM leg disagreed with the hardware on exactly one family, the only one that produces NaNs. The backend compiled float != as the ordered comparison, which is false whenever either side is a NaN. The interpreter, C and WebAssembly all answer true. The parity test for float edge cases had asked ==, < and > of a NaN, but never !=.

A local binding of a builtin’s name leaked into the whole program. The parity test for fma checks that a local function named fma belongs to its own scope. On LLVM, it took over every fma in the file. The top-level code was emitted with the table of lifted inner functions left over from whichever function had been emitted last. The same bug turned atan2 1.0 1.0 into 2.0 when some function defined its own atan2. My first regression test for it passed with the bug still in place, because only the last function’s table leaked. The test now puts the shadowing function last.

Looking next to the != bug turned up a third: f_min and f_max with a NaN argument returned NaN on the interpreter and WebAssembly, and on C and LLVM returned the other operand, depending on argument order. That one went into the list of known failures and was fixed shortly after. The fix copies OCaml’s own Float.min into every backend, step for step. My first version used a + b to produce the NaN, and a gate that compares bits caught it: a signalling NaN comes back from OCaml unchanged, and an addition quiets it.

What I would tell myself before starting

• Measure the upper bound before building the feature. Partial loads would have been a new builtin on five backends for 6%, and the variant that skipped the tail entirely – wrong answers and all – settled it in an afternoon.
• The rule was “don’t let the optimizer change the answer”, not “no fma”. Reading why the pragma existed is what turned a 1.4x win from forbidden into available.
• Emulate what you can reproduce exactly, and refuse the rest. That is a per-operation test, not a global policy, and fma was the first operation to pass it.
• Random testing does not reach ties. Build them on purpose.
• A gate’s harvest is often next to its subject. fma was right from the first day. What went red was the NaN test the gate used to fold its results, and the local binding its parity case used to check scoping.

fma has no real user in the tree yet beyond its tests and the measurement above; the natural one is a small 3D engine’s linear algebra, which already checks that every backend prints the same bits for correctly rounded operations. The detailed record, with every decision and its measurement, is docs/simd-go-comparison.md in the Mere repository.