How I write SIMD kernels_

I widened two neighboring dot-product kernels. The E5M2 kernel became 1.5% faster at p99, so the change stayed. The f32 kernel did not move, so I removed it.

Both changes used the same idea: carry more independent accumulators through the loop. Their different results came from the work around those accumulators. The compact kernel converts bytes while it scores them; the f32 kernel mostly loads and multiplies values that are already in the right format.

This article follows that difference from source code to emitted instructions, then checks it against recall and latency rather than treating a plausible SIMD argument as proof.

The targets come first

Before the compact lanes existed I fixed two targets: at least 99% recall and add p99 below half a millisecond. The acceptance run enforces both, so a kernel change cannot trade quality for speed without saying so. I read each result in the same order: recall, p95 and p99, then p50. The average is not part of the decision because it can hide the one-in-a-hundred insert the p99 budget exists to bound.

One piece of standing context before the kernels: every measurement in this post comes from one machine, an Intel Core Ultra 7 265 (20 cores), GCC 13.3.0, Ubuntu 24.04. The numbers are specific to that machine and toolchain; on different hardware the same changes would need to be measured again.

What the dot kernel does

The read kernel operates on E5M2, one byte with 1 sign bit, 5 exponent bits, and 2 mantissa bits. The format is compact and coarse, so the byte scores route the graph search and the survivors are reranked against their original f32 vectors. The conversion post covers the bit layout and recall trade.

The kernel problem starts after storage: x86 has no instruction that loads an E5M2 row and scores it directly. The loop has to turn bytes into f32 lanes while it computes the dot product, and any kernel change still has to preserve enough approximate ordering to pass the recall gate.

So the dot kernel converts while it scores. It loads 8 bytes and zero-extends each into its own 32-bit lane, giving every byte room to move. It shifts each lane left by 8, which parks the E5M2 bits exactly where a half-float keeps its top bits; the byte for 1.0, shifted left by 8, is bit-for-bit the half-float for 1.0. It packs the lanes down to 8 half words, and one hardware half-to-float instruction turns those into 8 f32 lanes, ready for ordinary FMA work against the query. It works eight values at a time because eight is the width the converter takes. Finding that exact sequence took five attempts in one day; the conversion tournament is a separate post.

The part of the sequence this post needs is the split it creates. Every iteration of the kernel now carries two kinds of work: conversion work, which manufactures f32 lanes out of bytes, and FMA work, which folds those lanes into a running score. Converting the next group of bytes does not depend on any FMA result, so the two streams are independent of each other. How much that independence is worth depends on how the loop is organized, and that is the accumulator question.

Accumulators

Set the conversion aside for a moment, because the accumulator question already exists in a plain f32 dot product. The naive loop keeps one running sum:

// one running sum: every fma depends on the one before it
float dot(const float* a, const float* b, int n) {
    float s = 0.0f;
    for (int i = 0; i < n; ++i)
        s = std::fma(a[i], b[i], s);
    return s;
}

Every FMA reads the sum the previous FMA produced, so the whole loop is one dependency chain, and it cannot run faster than one FMA latency per element no matter how wide the core is. The gap between two numbers is the whole story here: an FMA takes a few cycles to finish, but the core can start new, independent FMAs every cycle. Call the latency four cycles and the issue rate two per cycle: a single chain then delivers one FMA every four cycles on a machine built to start eight in that time. (Illustrative arithmetic to show the shape, not this core's exact numbers.)

The fix is independent partial sums:

// four independent chains, folded once at the end (n a multiple of 4)
float dot4(const float* a, const float* b, int n) {
    float s0 = 0, s1 = 0, s2 = 0, s3 = 0;
    for (int i = 0; i < n; i += 4) {
        s0 = std::fma(a[i+0], b[i+0], s0);
        s1 = std::fma(a[i+1], b[i+1], s1);
        s2 = std::fma(a[i+2], b[i+2], s2);
        s3 = std::fma(a[i+3], b[i+3], s3);
    }
    return (s0 + s1) + (s2 + s3);
}

Now there are four chains and each waits only on itself. Nothing about the FMA latency changed; the four latencies now overlap instead of queueing one behind the next. The chains run side by side through the whole loop and collapse into a single number only at the end, in the pairwise reduction; even that line is shaped by the same logic, because (s0 + s1) + (s2 + s3) is two levels of dependency where a left-to-right sum would be three.

The same arithmetic says the useful count is finite. With latency four and issue two, eight independent chains would saturate the FMA units on their own, but the loop also has loads to issue, and this loop will have conversion work too, so the practical number sits below the paper ceiling, and finding it takes a measurement.

I also ruled out waiting for the compiler to do it: it will not, at least not here. Splitting one float sum into four changes the order of the additions, float addition is not associative, and without fast-math flags the compiler is not allowed to reassociate on its own. The four sums have to be written out in the source. Fast-math flags would license the compiler to reassociate, but then the accumulator count changes with compiler versions and flag sets instead of with measurements. So the accumulator count is a decision the kernel author makes, and it can be wrong in both directions.

None of this needs SIMD, and in astral it predates it: the scalar E5M2 dot kernel had exactly this shape from its very first version, four partial sums at stride four, summed at the end, before a single vector instruction existed in the lane. When the kernel got its AVX2 version, each partial sum became a vector of 8 f32 lanes and the loop carried two vector accumulators, 16 values per iteration. Two was the shape the vector code arrived with, not a measured answer, and the E5M2 kernel had a specific reason to suspect two was too few.

Here is the two-accumulator kernel, and the loop clang-18 actually emits for it on x86-64. The widen is the five-instruction sequence the search settled on, zero-extend the bytes into 32-bit lanes, shift, pack the lanes back down to halves, then convert, and a single vfmadd231ps folds the widened lane into its accumulator:

// widen 8 E5M2 bytes to 8 f32 through the F16C unit
static inline __m256 widen_e5m2_x8(const uint8_t* p) {
    __m256i w = _mm256_slli_epi32(_mm256_cvtepu8_epi32(       // byte -> own 32-bit lane
                    _mm_loadl_epi64((const __m128i*)p)), 8);  // << 8 into half position
    __m128i h = _mm_packus_epi32(_mm256_castsi256_si128(w),
                                 _mm256_extracti128_si256(w, 1)); // pack lanes -> 8 halves
    return _mm256_cvtph_ps(h);                                   // 8 half -> 8 f32 (F16C)
}
// the E5M2 dot kernel, two vector accumulators (n a multiple of 16)
float dot_e5m2_2acc(const uint8_t* a, const float* b, int n) {
    __m256 s0 = _mm256_setzero_ps(), s1 = _mm256_setzero_ps();
    for (int i = 0; i < n; i += 16) {
        s0 = _mm256_fmadd_ps(widen_e5m2_x8(a + i),     _mm256_loadu_ps(b + i),     s0);
        s1 = _mm256_fmadd_ps(widen_e5m2_x8(a + i + 8), _mm256_loadu_ps(b + i + 8), s1);
    }
    return hsum256(_mm256_add_ps(s0, s1));   // hsum256 = extract-add + two hadds
}

; clang-18 -O2 -mavx2 -mfma -mf16c (x86-64), dot_e5m2_2acc inner loop
.LBB0_5:
    vpmovzxbd    ymm2, qword ptr [rdi + rcx]             ; 8 E5M2 bytes -> 8x u32
    vpslld       ymm2, ymm2, 8                           ; byte << 8, into half position
    vextracti128 xmm3, ymm2, 1                           ; high 4 lanes
    vpackusdw    xmm2, xmm2, xmm3                        ; pack 8x u32 -> 8x u16 (halves)
    vcvtph2ps    ymm2, xmm2                              ; 8 half -> 8 f32   (F16C)
    vfmadd231ps  ymm1, ymm2, ymmword ptr [rsi + 4*rcx]   ; acc0 += widened * b
    vpmovzxbd    ymm2, qword ptr [rdi + rcx + 8]
    vpslld       ymm2, ymm2, 8
    vextracti128 xmm3, ymm2, 1
    vpackusdw    xmm2, xmm2, xmm3
    vcvtph2ps    ymm2, xmm2
    vfmadd231ps  ymm0, ymm2, ymmword ptr [rsi + 4*rcx + 32] ; acc1 += widened * b
    add          rcx, 16
    cmp          rcx, rax
    jb           .LBB0_5

Two accumulators means two vfmadd231ps instructions into ymm0 and ymm1, each with its own widen stream feeding it. Going to four duplicates the complete widen-plus-FMA block twice more and advances by 32 components. The excerpt below shows only the four accumulator destinations; the full loop has four copies of the conversion sequence shown above:

; accumulator destinations in the 4-accumulator loop:
    vfmadd231ps ymm0, ymm4, ymmword ptr [rsi + 4*rcx]      ; acc0
    vfmadd231ps ymm1, ymm4, ymmword ptr [rsi + 4*rcx + 32] ; acc1
    vfmadd231ps ymm3, ymm4, ymmword ptr [rsi + 4*rcx + 64] ; acc2
    vfmadd231ps ymm2, ymm4, ymmword ptr [rsi + 4*rcx + 96] ; acc3

A standalone Clang 18 reproduction kept all four sums and the two widening temporaries in registers, with no stack reference or vector spill in the loop. Register pressure is the next boundary to inspect, not a measured explanation for where this kernel stops.

Two accumulators to four

The reason is the conversion stream. A plain f32 dot kernel fills the space between FMA completions with loads and not much else. The E5M2 kernel issues real instructions in that space, the zero-extend, the shift, the pack, the convert, and none of them wait on any FMA chain. With two accumulators, once both chains are stalled on their previous results, the conversion stream is the only independent work in flight. With four, the scheduler holds four FMA chains plus the conversion stream at once, so more of every chain's latency gets covered by useful work before the loop finally collapses into its reduction. The total arithmetic is unchanged; more of it is independent at any given moment.

So the evening after the conversion path settled, I widened both E5M2 dot kernels, the one that scores stored bytes against an f32 query and the one that scores stored bytes against other stored bytes, both of which the search needs, from two vector accumulators to four: 32 values per loop iteration, four sums, pairwise reduce at the end. The q8 lane's dot kernel had taken the same widening the day before, so the widening was already the house pattern; the dot kernels across lanes share one shape, same stride, same reduction, a different conversion in the middle, and the E5M2 lane was the lane with the most latency to hide.

I judged the widening through the insert benchmark, because inserting a vector into the graph means searching the graph for the new vector's neighbors, and the search's time goes into exactly these kernels. The fixture inserts 100k vectors and reads the latency percentiles across the individual inserts. A dot kernel runs many times inside each insert, so a kernel change arrives at the percentiles shifted and smoothed rather than raw; small kernel wins show up as small tail wins. Measured there: p95 went from 1,926,399 to 1,899,105 ns, p99 went from 2,079,024 to 2,048,914 ns, p50 stayed basically flat, and recall did not move.

The widening removed about 1.5% from p99 and nothing from the median or recall, so it stayed. The run did not say where additional accumulators would stop helping.

Where the widening stops

The latency-hiding argument has no built-in stopping point. If four independent chains beat two, the same reasoning argues for eight, so the stopping point has to come from somewhere else. Every extra accumulator occupies a vector register for the entire loop, the conversion sequence needs scratch registers of its own, and AVX2 gives the two of them sixteen register names to share; the moment sums start spilling to the stack, the loop is doing memory traffic just to keep its own accumulators alive. Every doubling of the stride also grows the tail of leftover elements the loop has to finish scalar. (Those are the mechanics in general, not a measured diagnosis of where this kernel's ceiling sits.)

Two measurements placed the line. Widening to four helped the E5M2 dots, as above. Around the same period I tried a wider unroll on the plain f32 dot kernel, and the measurement rejected it, so it came back out and the f32 kernel kept its width. That result fits the story, because a plain f32 loop has no conversion stream to overlap and so has the least to gain, but I reverted it because of the number, not because it fit a story. And the widening never hardened into a policy: the L2 distance kernels that sit beside the dot kernels still run two accumulators wide. Each kernel has its own mix of conversion work, arithmetic, and loads, and that mix is why four helped in the E5M2 dots.

So the widths across neighboring kernels ended up uneven: four for the E5M2 dots, the old width for f32 after a measured attempt, two for the L2 distances. I cannot say why the L2 kernels are still at two; the record does not say whether a wider version lost a benchmark or was never tried. That missing result is a reason to measure the next change, not permission to invent the old one.

The read path now overlaps conversion with four independent sums. The write path has a different problem: data-dependent branches.

Branchless on the write path

Vectors enter a lane through its quantizer: insert takes f32 components in and produces E5M2 bytes. Insert is the gated path, the half-millisecond aim and the p99 ceiling live here, and the quantizer touches every component of every inserted vector. A vector is 384 components in my fixtures, so every insert is hundreds of quantizations before the search even starts. So it is a hot kernel, not one-time setup work.

Quantizing a float down to eight bits is a pile of cases. A component can exceed the format's largest finite value and need saturating. The sign has to come off and go back on. An in-range value lands between two of the format's four steps per power of two and needs rounding. Scalar code writes each case as an if and thinks nothing of it. A SIMD kernel cannot, because the eight components in a vector register share one instruction stream: a branch on one component's value means either breaking the vector apart or branching on a summary of all eight lanes. Worse, these branches depend on the data. Whether a component needs saturating is decided by whatever embedding happens to arrive, so the predictor has nothing stable to learn, and a mispredicted branch discards a pipeline's worth of work. On the gated path that shows up as insert latency that varies with the data, which is the pattern the p99 ceiling is meant to catch.

So the quantizer computes every case and selects, compare-and-select instead of branches. The shape, simplified to the idea:

// branchless saturate: round every lane, then overwrite the overflow lanes
__m256 saturate_select_x8(__m256 x, __m256 rounded_code) {
    __m256 over = _mm256_cmp_ps(x, _mm256_set1_ps(57344.0f), _CMP_GT_OQ);
    return _mm256_blendv_ps(rounded_code, _mm256_set1_ps(123.0f), over); // 0x7B
}

; clang-18 -O2 -mavx2 (x86-64), saturate_select_x8
    vbroadcastss ymm2, dword ptr [rip + .LCPI2_0]   ; 57344.0 = E5M2 max finite
    vbroadcastss ymm3, dword ptr [rip + .LCPI2_1]   ; 123.0   = 0x7B, top finite code
    vcmpltps     ymm0, ymm2, ymm0                   ; mask = (57344 < x), i.e. x > max
    vblendvps    ymm0, ymm1, ymm3, ymm0             ; pick 0x7B where over, else rounded
    ret

Every lane computes both outcomes of every case. That sounds wasteful, but the compare and the select are one instruction each, and the kernel's running time stops depending on the data entirely. The sign case gets the same treatment: strip the sign bits with a mask up front, quantize the magnitudes, put the sign bits back with an or at the end. Astral's f32-to-E5M2 quantizer is an AVX2 kernel of exactly this shape, compare and select throughout.

The saturation constant is 0x7B, the format's top finite code, the 57344.0 value, and everything above it saturates there rather than rolling into the four codes past it, whose half-float shapes decode to infinity and NaN. Those four codes are the one place astral's two decoders disagree: the hardware converter would turn them into inf or NaN, while the fallback decoder, a 256-entry lookup table kept for machines without the converter, clamps them to plus or minus 57344. Because the quantizer never emits the four codes, the disagreement is unreachable through astral's own write path. That argument sits in a comment next to the code, because a divergence that no input can currently reach is exactly the kind of fact that gets rediscovered as a bug later if nobody wrote down why it is safe.

So the write path, end to end, is a handful of masked steps: strip the signs, round every lane, saturate anything past 0x7B, put the signs back, and never emit a code the two decoders disagree on. Every case resolves with a mask instead of a branch, and the edge codes are documented next to the code they constrain.

Rounding by band

E5M2 is not the index's only compact format. Its smaller sibling E3M2 uses six bits per component, 1 sign bit, 3 exponent bits, 2 mantissa bits: the same rough four-step mantissa with a quarter of the exponent values. The same write-path job comes with it: round every f32 component to the nearest value the format can represent.

The first version of that rounding searched for its answer. Every magnitude an E3M2 code can represent fits in a table of 31 or 32 entries, and the quantizer scanned that table for each component and kept the nearest entry. Thirty-odd compares per component would be fine on a cold path, and this path is not cold.

Inserts quantize whole vectors, and at query time the query itself gets quantized to E3M2, so that the search can score stored components against it in pure integer arithmetic, int16 against int16, with no float conversion inside the search. That is the point of quantizing the query at all: one conversion up front instead of one at every candidate the search visits. My fixtures run 384-dimensional vectors, so quantizing one query by scan is up to 384 by 32 comparisons, roughly twelve thousand, before the search begins. (Sizing arithmetic, not a measurement.)

The replacement rounds by band. Three exponent bits split the representable magnitudes into bands, one band per exponent value, four mantissa steps inside each: a band for values between 1.0 and 2.0, a band for values between 2.0 and 4.0, and so on across the format's range. An incoming f32 already carries its band, because its own exponent bits are the band index. So the kernel reads the band straight off the input's exponent, rounds two bits of mantissa within that band, and is done: constant work per component. The values it rounds to are the same nearest values the scan chose. What changed is how they get found: the kernel now reads the band index that was already in the input, instead of searching for it.

The E5M2 quantizer resolves its cases with compares and selects; the E3M2 one reads its rounding band out of the input's own exponent. Neither branches on data or reads a table on the hot path, and both landed through the same acceptance gates as everything else that month.

Reading what the compiler emitted

A benchmark measures the binary that exists, not the path I meant to compile. The first E5M2 run used the fallback because the F16C build flag was missing; every answer was still correct. Reading the loop found no vcvtph2ps, so the flag was fixed and generated code became part of the check.

For each kernel I now verify the few instructions the argument depends on: the conversion is present, the multiply and add fused, and compare-and-select did not turn back into a branch. The compiler-specific F16C guard and the separate runtime CPUID/XCR0 check belong to the conversion post; they do not change the two accumulator verdicts here.

The two changes from the opening looked equally plausible before the runs and ended in opposite places. Four accumulators remain in the E5M2 dot kernels because p99 fell 1.5% with recall and the median unchanged. The wider f32 loop remains reverted because its measurements did not move. The unresolved item is narrower: the L2 kernels still use two accumulators, and the record does not say whether four lost a benchmark or was never tried.