What Zero-Cost Abstractions Are Worth

Writing the fastest linear algebra, counting cycles, chasing a regression

Artem Zyktin

tenvex repo

TL;DR. I wrote a small SIMD library for vec4 and checked with a disassembler whether C++ zero-cost abstractions really cost nothing. Answer: byte-to-byte identical to hand-written intrinsics. Then I chased a +39% regression on arm64, convinced I'd caught an abstraction leak, tested several hypotheses, and smashed each one against the disassembler. What it actually was, I figured out, but not immediately and not where I was looking. Along the way, it became clear that "costs nothing" and "beats naive code" are two different things: sometimes the compiler wins, sometimes we do. About registers, disassemblers, and how a beautiful hypothesis loses to facts.

One night, staying up late, I decided: why not write the fastest linear algebra library. Just kidding, of course that's too ambitious for one developer. Entire teams have been fighting over this for decades (e.g. DirectXMath, GLM, Eigen) hardly I could compete with them (well, we'll see). In reality, this thought came to me for several much more down-to-earth reasons.

First, it's a great excuse to brush up on linear algebra. Especially in the context of computer graphics, where I've always had gaps. Second, I'd never dealt with SIMD operations and wide registers (SSE, AVX, and so on) closely, and this is an excellent topic to get hands-on with such things. Third, I wanted to play with optimizations. And templates: I love C++ for its flexibility, and I'm good at metaprogramming and templates, and I enjoy tinkering with them.

Also, I wanted to test a promise that C++ has been repeating for decades, about zero-cost abstractions: convenient, readable code that costs nothing at runtime, for real, with benchmarks, tests, and a disassembler. Does it really cost nothing? That's what I wanted to feel. Let's go.

We'll start with the very basics, the vector. As it turns out, this is quite enough to cover a wide spectrum of low-level computing programming topics. We won't touch linear algebra theory; we're interested not in the mathematics, but in how it maps to hardware. Can an expression like dot(A,B)*t*C compile into the same instructions I'd write by hand with intrinsics (built-in functions that translate directly into low-level instructions)? Not "comparably fast," but literally the same disassembly. If yes, then I have convenient math without paying for it. If not, then I'm paying for beauty, and I need to know how much. The discussion will mainly revolve around several things: the hardware we compute on, and low-level memory work (even lower than you might imagine). True speed lives in registers. Everything else in this article is a consequence of that fact.

Part I. The Promise of Zero-Cost Abstractions

Where Fast Code Lives: In Registers, Not in Memory

When people talk about optimization in code, they usually mean the right data layout. The main bottleneck is memory, it's the slowest, and good data layout is the first thing to think about. Any computing library will push you toward this: pack data densely, uniformly, and process it pipelined to fit in cache and minimize trips to RAM. We'll deliberately skip this layer now and look one level lower, at computing individual expressions.

And here, essentially, it's the same story. Against the computations themselves (1-4 cycles), an L1 cache access is a comparable cost (around 4-5 cycles). So "just" going to cache is already expensive for us. Therefore, the main task is to keep the expression in registers. And things like spills (unloading a value from a register to memory; this is expensive) already change the performance mode, tens and hundreds of cycles if it reaches RAM.

Table 1: Operation latency for Ice Lake microarchitecture (uops.info, uops.info/cache). Order of magnitude; comparable microarchitectures are similar.

#include <immintrin.h>

static inline __m128 kernel(__m128 x, __m128 y)
{
    return _mm_add_ps(_mm_mul_ps(x, y), y);
}

__m128 clean(__m128 a, __m128 b, __m128 c)
{
    __m128 r = kernel(a, b);
    return kernel(r, c);
}

__m128 opaque_kernel(__m128 x, __m128 y);

__m128 spilled(__m128 a, __m128 b, __m128 c)
{
    __m128 r = opaque_kernel(a, b);
    return opaque_kernel(r, c);
}

_Z5cleanDv4_fS_S_ :
.LFB7302 :
    .cfi_startproc
    endbr64
    mulps xmm0, xmm1
    addps xmm0, xmm1
    mulps xmm0, xmm2
    addps xmm0, xmm2
    ret
    .cfi_endproc

_Z7spilledDv4_fS_S_ :
.LFB7303 :
    .cfi_startproc
    endbr64
    sub rsp, 24
    .cfi_def_cfa_offset 32
    movaps XMMWORD PTR[rsp], xmm2
    call _Z13opaque_kernelDv4_fS_@PLT
    movaps xmm1, XMMWORD PTR[rsp]
    add rsp, 24
    .cfi_def_cfa_offset 8
    jmp _Z13opaque_kernelDv4_fS_@PLT
    .cfi_endproc

Listing 1: C++ code and its assembly: clean - pure function without spills to memory, spilled - additional calls and unloading temporaries to memory.

Here I need to be honest: registers amount is finite, and besides our expression they hold other computations and service code - this is called register pressure. Sometimes a spill is a forced and absolutely correct step. So going forward I'll carefully distinguish two things: a forced spill due to register shortage, and an extra spill because optimization didn't squeeze hard enough. Only the second is worth fighting.

Why the Compiler Is So Good at Scalar Code and Whether It Can Be Caught

On scalar expressions, the compiler is almost unbeatable. Decades of optimizations tuned for specific architectures behind it, it excellently recognizes patterns, reorders code, schedules it for a specific pipeline, allocates registers, and picks instructions. In a typical scalar expression, you can't outplay it by hand.

What's important to accept right away. The goal of our manual code is not to beat the compiler on scalars. There we'll play to a draw at best, and more likely lose. But there's a subtlety: the compiler is strong as long as it has guarantees. Where it loses information about semantics - aliasing (when two pointers may point to the same memory area or overlapping areas), escape analysis (compiler analysis determining whether a variable's address can escape outside a function; if not, the variable can be kept in a register, but if it can, it needs to be spilled to memory), subtleties of floating-point arithmetic (e.g., handling NaN situations) - it is forced to play it safe and be conservative.

float dot3(const float* a, const float* b)
{
    return a[0] * b[0] + a[1] * b[1] + a[2] * b[2];
}

void axpy3(float* dst, const float* src, float k)
{
    dst[0] += k * src[0];
    dst[1] += k * src[1];
    dst[2] += k * src[2];
}

void axpy3_r(float* __restrict dst, const float* __restrict src, float k)
{
    dst[0] += k * src[0];
    dst[1] += k * src[1];
    dst[2] += k * src[2];
}

_Z4dot3PKfS0_:
.LFB0 :
    .cfi_startproc
    endbr64
    vmovss xmm0, DWORD PTR 4[rdi]
    vmulss xmm0, xmm0, DWORD PTR 4[rsi]
    vmovss xmm1, DWORD PTR [rdi]
    vfmadd231ss xmm0, xmm1, DWORD PTR [rsi]
    vmovss xmm2, DWORD PTR 8[rdi]
    vfmadd231ss xmm0, xmm2, DWORD PTR 8[rsi]
    ret
    .cfi_endproc

_Z5axpy3PfPKff :
.LFB1 :
    .cfi_startproc
    endbr64
    vmovss xmm1, DWORD PTR [rsi]
    vfmadd213ss xmm1, xmm0, DWORD PTR [rdi]
    vmovss xmm2, DWORD PTR 8[rdi]
    vmovss DWORD PTR [rdi], xmm1
    vmovss xmm1, DWORD PTR 4[rsi]
    vfmadd213ss xmm1, xmm0, DWORD PTR 4[rdi]
    vmovss DWORD PTR 4[rdi], xmm1
    vfmadd132ss xmm0, xmm2, DWORD PTR 8[rsi]
    vmovss DWORD PTR 8[rdi], xmm0
    ret
    .cfi_endproc

type _Z7axpy3_rPfPKff, @function
_Z7axpy3_rPfPKff :
.LFB2 :
    .cfi_startproc
    endbr64
    vmovsldup xmm2, xmm0
    vmovq xmm3, QWORD PTR [rdi]
    vmovq xmm1, QWORD PTR [rsi]
    vmovss xmm4, DWORD PTR 8[rdi]
    vfmadd132ss xmm0, xmm4, DWORD PTR 8[rsi]
    vmovq xmm2, xmm2
    vfmadd132ps xmm1, xmm3, xmm2
    vmovss DWORD PTR 8[rdi], xmm0
    vmovlps QWORD PTR [rdi], xmm1
    ret
    .cfi_endproc

Listing 2.

Look at Listing 2, dot3. Here the compiler didn't vectorize the computation - all operations remained scalar. But it did something else: it merged several multiplications and additions into one multiplication and two FMAs - this is FMA contraction (Fused Multiply-Add - an instruction that does multiplication and addition in one operation: a*b+c; instead of two cycles - one). Now look at axpy3 and axpy3_r. Seems like identical expressions, yet the __restrict annotation gives so much: instead of loading data sequentially and computing them one by one (ss - scalar single-precision), the compiler loads two components at once and computes them with a packed instruction (ps - packed single-precision, not ss!), i.e., it starts vectorizing on its own. Without the guarantee, it must assume that writing to dst[0] could have touched src[1] (the pointers might overlap), so it keeps everything sequential and scalar. With __restrict, the fear is lifted, and it voluntarily packs a pair of arrays into SIMD. Note that it's specifically the pair (x, y) that gets packed with one ps instruction, while the third component (z) is finished off scalarly (ss) - an odd tail doesn't fit into SIMD. This is where our bread is buttered. Not "how to beat the compiler," but "where it stops on its own." It stops not because it's stupid, but because it lacks proof. And part of our job is to provide that proof. But more on that later.

Can We Catch Hand-Written Intrinsics?

Yes - and that's the bar we're aiming for. The goal is for the abstraction to output byte-to-byte the same code a human would write with intrinsics. Not "close," not "almost," but identical machine code. This, by the way, is a convenient criterion: we compare disassemblies, not feelings.

Why is this even possible? Intrinsics are the thinnest wrapper around processor instructions. If an abstraction compiles to the same instructions, then by definition there's no runtime difference between it and hand-written code. What gets in the way? Every layer is expression nodes, wrappers, references - that's what the optimizer must untangle back to the same intrinsics. Untangled - we get zero-cost: convenient code that's pleasant to read and write, and that doesn't lose on speed. Didn't untangle - we get a leak, an abstraction we pay for.

#include <immintrin.h>

static inline __m128 dot3(__m128 l, __m128 r)
{
    __m128 t = _mm_mul_ps(l, r);
    __m128 ty = _mm_shuffle_ps(t, t, _MM_SHUFFLE(1, 1, 1, 1));
    __m128 xy = _mm_add_ss(t, ty);
    __m128 tz = _mm_movehl_ps(t, t);
    __m128 xyz = _mm_add_ss(xy, tz);
    return _mm_shuffle_ps(xyz, xyz, _MM_SHUFFLE(0, 0, 0, 0));
}

static inline __m128 norm3(__m128 v)
{
    __m128 len = _mm_sqrt_ps(dot3(v, v));
    return _mm_blend_ps(_mm_div_ps(v, len), v, 0b1000);
}

__m128 expr(__m128 a, __m128 b, __m128 c)
{
    __m128 d = dot3(c, b);
    __m128 v = _mm_add_ps(a, _mm_mul_ps(b, _mm_set1_ps(2.0f)));
    return _mm_add_ps(
        _mm_mul_ps(norm3(v), d),
        _mm_mul_ps(c, _mm_set1_ps(3.0f))
    );
}

#include <tenvex.h>

tnvx::vec4 expr(tnvx::vec4 a, tnvx::vec4 b, tnvx::vec4 c)
{
    using namespace tnvx;
    return norm3(a + b * 2.0f) * dot3(b, c) + c * 3.0f;
}

g++ -std=c++20 -O3 -msse4.1 -mfma -S -masm=intel -fno-asynchronous-unwind-tables sample3.cpp -o sample3.asm

g++ -std=c++20 -O3 -msse4.1 -mfma -S -masm=intel -fno-asynchronous-unwind-tables -I./projects/tenvex/src/tenvex sample3_tenvex.cpp -o sample3_tenvex.asm

_Z4exprDv4_fS_S_:
.LFB7301 :
    .cfi_startproc
    endbr64
    vmovaps xmm5, xmm0
    vmulps xmm0, xmm2, xmm1
    vpermilps xmm3, xmm0, 85
    vaddss xmm3, xmm0, xmm3
    vmovhlps xmm0, xmm0, xmm0
    vaddss xmm3, xmm3, xmm0
    vbroadcastss xmm0, DWORD PTR .LC1[rip]
    vfmadd231ps xmm5, xmm1, xmm0
    vbroadcastss xmm0, DWORD PTR .LC3[rip]
    vpermilps xmm4, xmm3, 0
    vmulps xmm1, xmm5, xmm5
    vpermilps xmm3, xmm1, 85
    vaddss xmm3, xmm1, xmm3
    vmovhlps xmm1, xmm1, xmm1
    vaddss xmm1, xmm3, xmm1
    vpermilps xmm1, xmm1, 0
    vsqrtps xmm1, xmm1
    vdivps xmm1, xmm5, xmm1
    vblendps xmm1, xmm1, xmm5, 8
    vmulps xmm4, xmm4, xmm1
    vfmadd132ps xmm0, xmm4, xmm2
    ret
    .cfi_endproc

_Z4exprN4tnvx4vec4ES0_S0_:
.LFB800 :
    .cfi_startproc
    endbr64
    vmovaps xmm5, xmm0
    vmulps xmm0, xmm2, xmm1
    vpermilps xmm3, xmm0, 85
    vaddss xmm3, xmm0, xmm3
    vmovhlps xmm0, xmm0, xmm0
    vaddss xmm3, xmm3, xmm0
    vbroadcastss xmm0, DWORD PTR .LC1[rip]
    vfmadd231ps xmm5, xmm1, xmm0
    vbroadcastss xmm0, DWORD PTR .LC3[rip]
    vpermilps xmm4, xmm3, 0
    vmulps xmm1, xmm5, xmm5
    vpermilps xmm3, xmm1, 85
    vaddss xmm3, xmm1, xmm3
    vmovhlps xmm1, xmm1, xmm1
    vaddss xmm1, xmm3, xmm1
    vpermilps xmm1, xmm1, 0
    vsqrtps xmm1, xmm1
    vdivps xmm1, xmm5, xmm1
    vblendps xmm1, xmm1, xmm5, 8
    vmulps xmm4, xmm4, xmm1
    vfmadd132ps xmm0, xmm4, xmm2
    ret
    .cfi_endproc

Listing 3: With a different expression writing order, the compiler will arrange instructions differently - but it will be an equivalent assembly, with the same instruction set. Here both versions are arranged so that the disassemblies match literally.

So, the bar is taken, I got zero-cost abstraction, byte-to-byte like hand-written intrinsics code (Listing 3). And without spilling temporaries to memory. We'll discuss benchmarks below.

Why Templates Are an Abstraction Without Overhead

The mechanism is this. A regular abstraction is virtual calls, function pointers, type erasure - we pay at runtime for dispatching and inability to inline. A template, by its nature, is closer to substitution: for each specific type, the compiler generates separate, specialized code that it sees in full. Instead of abstractly written template, the output is concrete code for a concrete task. The key point is that the computation structure (which operations and in what order) becomes part of the type. The compiler sees the entire tree statically, and therefore can inline and fold it all the way through. No runtime decisions, nothing to dispatch.

And here we come to the main technique in my math - expression templates (ET). This is an abstraction over abstraction, a template that builds not a value, but a type-expression. dot(A,B)*t*C is not evaluated on the spot, but assembled into a type like Mul< Mul< Dot<vec4,vec4>, float >, vec4 >. The computation itself is the compiler unfolding this type.

// --- traits: mark what is an expression ---
template<typename> inline constexpr bool is_vec_expr = false;
template<typename> inline constexpr bool is_scalar_expr = false;

// node is stored by value if small and trivial; otherwise - by reference
template<typename T>
inline constexpr bool is_stored_by_value_v =
    std::is_trivially_copyable_v<T>
    && std::is_trivially_destructible_v<T>
    && sizeof(T) <= 16;

template<typename T>
using tnvx_ref_or_value_t =
    std::conditional_t<is_stored_by_value_v<T>, T, const T&>;

// --- CRTP base ---
template<typename Derived>
struct Expr
{
    const Derived& self() const noexcept
    {
        return static_cast<const Derived&>(*this);
    }
};

// --- expression concept: can eval() into a register and knows its type ---
template<typename T>
concept expression = requires(T t)
{
    { t.eval() } -> std::same_as<vf4>;
    { t.self() } -> std::same_as<const T&>;
} && std::is_base_of_v<Expr<T>, T>;

// --- multiplication node ---
template<expression L, expression R>
struct Mul : Expr<Mul<L, R>>
{
    Mul(const L& l, const R& r) noexcept : _l(l), _r(r) {}

    vf4 eval() const noexcept
    {
        return _mm_mul_ps(_l.eval(), _r.eval());
    }

private:
    tnvx_ref_or_value_t<L> _l;
    tnvx_ref_or_value_t<R> _r;
};

Listing 4.

In the end, the expression assembly scheme looks like this:

Figure 1: expression tree building scheme in tenvex

The cost, of course, doesn't disappear, but moves to compile-time: compilation time, type bloat, diagnostic wall-of-text. Zero-cost at runtime is not zero-cost for the developer, and the library author pays for the user's convenience.

What Zero-Cost Abstraction Is and How to Consume It

Zero-cost is a promise: you don't pay for what you don't use, and for what you do use, you don't pay more than you would for a manual implementation. And this promise is quite conditional. Because zero-cost is a property not of the code itself, but of the pair code + compiler (more precisely, of the specific backend - the part of the compiler generating code for a specific architecture, be it x86_64 or ARM64). In our task, it consists of several independent things: function inlining, SROA (Scalar Replacement of Aggregates - an optimization that breaks down a complex structure into separate variables and places them in registers instead of memory), register promotion (moving values from memory to registers), and DSE (Dead Store Elimination - an optimization that removes writes to memory if they are never read; cleans up garbage). If even one fails to come together - zero-cost is gone, the abstraction leaks.

Table 2: median values over 30 benchmarks, including throughput benchmarks where the next computation didn't depend on the previous one, and latency benchmarks where each computation depended on the previous one. manual - hand-built from the same kernels that tenvex uses; overhead = tenvex/manual.

And since the promise lives in the optimizer, it can be verified and it can be broken. Hence begins the second part of the story - the fight with compilers. Essentially, it's exactly about either giving the optimizer those missing guarantees, or discovering that a specific backend, one of the passes, doesn't pull through.

Looking at the comparison table, we see that the zero-cost promise is cracking, but we'll deal with that in the second part of the article.

Zero-cost abstraction is not what you wrote, but what the optimizer agreed to implement.

Part II. Fighting Compilers

A brief reminder: zero-cost is a property of the pair code + compiler, and it holds on the combination of optimizations (inline, SROA, register promotion, DSE). If even one fails to come together, the abstraction leaks. The comparison table showed one outlier cell: GCC on aarch64 - +39% against hand-written intrinsics code. That's the one we'll figure out.

At first I decided it was a hole in the library. It turned out to be more complicated, and along the way I made several mistakes. These mistakes are worth showing, because they're more instructive than the fix itself. Each is about how a convenient hypothesis loses to the disassembler.

Investigation Tools

Before chasing the bug, let's talk about tools. There are essentially four, and without the first three, any reasoning about performance hangs in the air. The first is Google Benchmark, and here it's important to immediately distinguish two modes: a latency benchmark runs a serial chain where each computation depends on the previous one - it measures the length of the critical path; a throughput benchmark runs independent iterations that the processor overlaps in the pipeline - it measures throughput. These are different numbers about different things, and confusing them is a mistake; what exactly you're measuring is decided by how the result leaks out through a sink like DoNotOptimize.

The second tool is the disassembler. Everything I claim going forward is verified not by a feeling of "got faster," but by comparing disassemblies: is there a stack access, is there a dead store left, has a node been simplified.

The third is case reduction: the ability to step-by-step reduce the real spilling translation unit to a minimum without losing the bug along the way. More on that in the unfolding section when it's needed.

And the fourth tool is Claude AI. I'm not too strong in asm - I know the very basics, but making sense of wall-of-text disassemblies and searching for specific facts in them was helped by Claude. It also helped me prepare benchmarks and specific code to test various hypotheses. And it took on writing tests for tenvex based on my requirements description.

static void BM_Dot3_Latency(benchmark::State& state)
{
    vec4 a { 1.0f, 2.0f, 3.0f, 0.0f };
    vec4 b { 0.5f, 0.5f, 0.5f, 0.0f };
    float acc = 0.0f;

    for (auto _ : state)
    {
        vec4 bb = b + vec4{ acc, 0.0f, 0.0f, 0.0f };
        float d = dot3(a, bb);
        acc = d * 1e-7f;
        benchmark::DoNotOptimize(acc);
    }
}

BENCHMARK(BM_Dot3_Latency);

static void BM_Dot3_Throughput(benchmark::State& state)
{
    const auto a = make_vecs(1024, 1);
    const auto b = make_vecs(1024, 2);
    std::size_t i = 0;

    for (auto _ : state)
    {
        float d = dot3(a[i], b[i]);
        benchmark::DoNotOptimize(d);
        i = (i + 1) & 1023;
    }
}

BENCHMARK(BM_Dot3_Throughput);

Listing 5: examples of tenvex benchmarks on Google Benchmark.

These two benchmarks measure two physically different characteristics of the same operation, and confusing them is a mistake that later makes numbers not match. In short: latency answers the question of how long one value takes to compute if the next can't start without the previous one, while throughput answers how many independent such operations I can push per unit time when there are no dependencies between them. Where does the difference come from at the hardware level? Take a regular mulps: its latency is ~3-4 cycles, but throughput is ~0.5 cycles, meaning the execution block is pipelined and can accept a new mul every cycle (or even two per cycle), though each individual result is only ready after 4. If you measure only latency, you'll see "4" and decide the operation is expensive. If only throughput, you'll see "0.5" and decide it's almost free. Both numbers are correct, they just answer different questions.

Miss #1: "There Are Spills in the Hot Path"

The first specific hypothesis was straightforward: +39% means spills in the computation chain - str/ldr shuttling values to the stack and back. Simple logic: overhead is extra memory accesses. I opened the asm, already knowing what I'd see. And I didn't. The computation chain turned out to be clean, entirely in vector registers, not a single load back from the stack. Instead of expected spills, there was a dead stack frame: writes to [sp] that no one then reads, several stp q, q, [sp] per iteration. This isn't a spill in the sense of "ran out of registers," values calmly live in registers. These are dead store that the DSE pass didn't squeeze and didn't delete. A different beast, and it heals differently.

#include <tenvex.h>
#include <cstddef>

using namespace tnvx;

__attribute__((noinline))
void loop_et(const vec4* va, const vec4* vb, const vec4* vc, vf4* out, std::size_t n)
{
    for (std::size_t i = 0; i < n; ++i)
        out[i] =
            (norm3(va[i] + vb[i] * 2.0f) *
             dot3(vb[i], vc[i]) +
             vc[i] * 3.0f)
            .eval();
}

g++ -std=c++20 -O3 -I./tenvex/src/tenvex -S sample4.cpp -o sample4.asm

_Z7loop_etPKN4tnvx4vec4ES2_S2_P13__Float32x4_tm:
.LFB4063 :
.cfi_startproc
    cbz x4, .L7
    sub sp, sp, #128
    .cfi_def_cfa_offset 128
    mov x5, 0
    fmov v31.4s, 2.0e+0
    add x8, sp, 96
    fmov v30.4s, 3.0e+0
    add x9, sp, 32
    mov x7, sp
    mov x6, 0
    .p2align 5, , 15
.L2:
    ldr q29, [x1, x5]
    str q30, [sp, 80]
    adrp x10, .LC0
    ldr q28, [x0, x5]
    str q29, [x8]
    add x6, x6, 1
    ldr q27, [x2, x5]
    fmla v28.4s, v29.4s, v31.4s
    ldr q25, [x10, #:lo12:.LC0]
    stp q29, q27, [x9]
    fmul v26.4s, v28.4s, v28.4s
    fmul v29.4s, v29.4s, v27.4s
    stp q27, q25, [x7]
    movi v24.4s, 0x40, lsl 24
    str q24, [x8, 16]
    ins v26.s[3], wzr
    faddp v26.4s, v26.4s, v26.4s
    ins v29.s[3], wzr
    faddp v29.4s, v29.4s, v29.4s
    fmul v27.4s, v27.4s, v30.4s
    faddp v26.4s, v26.4s, v26.4s
    faddp v29.4s, v29.4s, v29.4s
    dup v26.4s, v26.s[0]
    fsqrt v26.4s, v26.4s
    fdiv v26.4s, v28.4s, v26.4s
    ins v26.s[3], v28.s[3]
    fmla v27.4s, v26.4s, v29.s[0]
    str q27, [x3, x5]
    add x5, x5, 16
    cmp x4, x6
    bne .L2
    add sp, sp, 128
    .cfi_def_cfa_offset 0
    ret
.L7:
    ret
    .cfi_endproc

Listing 6

And here my hypothesis diverged from reality. There were no spills in the sense I expected: "computed, ran out of registers, spilled, then loaded back." The computation itself is clean, values live in registers, nothing evicts them. But at the same time, the function for some reason allocates 128 bytes on the stack and writes there six times, two of them - stp, two registers per instruction. I searched for each of these writes on the read side and found none: not a single ldr back from these addresses. Values go to the stack and die there, no one reads them. This isn't a spill from register pressure, these are dead writes that for some reason weren't deleted. And DoNotOptimize, and explicit .eval() in vf4 before writing didn't remove the frame, the result consumer was not at fault. Why the compiler left them, I didn't yet understand at this stage.

"Slow" and "extra memory accesses" are not synonyms, and the predicted mechanism may not be the one. I expected spills from register shortage, but found dead writes during clean computation. Look at the disassembly, not at the model in your head.

Miss #2: "Storage Policy Is the Lever"

This is the most instructive of the mistakes. ET nodes store operands via tnvx_ref_or_value_t: small leaves are stored by value, large subexpressions - by const&. The hypothesis suggested itself: it's the const& storage that forces nodes to hold references to temporary nodes that need to be materialized (spilled to memory); switch everything to by-value and the frame will collapse. I assembled a toy example entirely on by-value, and it confirmed it: clean asm, no frame. Elated, I switched the real library to by-value.

template<typename T>
inline constexpr bool is_stored_by_value_v = std::is_trivially_copyable_v<T>
&& std::is_trivially_destructible_v<T>
&& sizeof(T) <= 16;

template<typename T>
using tnvx_ref_or_value_t =
    std::conditional_t<is_stored_by_value_v<T>, T, const T&>;

I changed to

template<typename T>
inline constexpr bool force_value_storage_v = true;

template<typename T>
using tnvx_ref_or_value_t =
    std::conditional_t<force_value_storage_v<T>, T, const T&>;

Listing 7.

Easy fix, a good point for changes, change in one place, and if something needs to be done with this in the future, we know what and where.

And this was a footgun... Backfire. +39% turned into +147-170%. The frame didn't collapse, it grew: from ~128 bytes to ~400, and instead of four dead store there were fifteen. The clean toy lied: it was so sterile that it couldn't show the main thing - real tenvex on by-value was not better, but really worse. Storage policy wasn't the lever at all, I was kneading the wrong dough.

Table 3: median values over 30 benchmarks, including throughput benchmarks where the next computation didn't depend on the previous one, latency benchmarks where each computation depended on the previous one. manual - hand-built from the same kernels that tenvex uses; overhead = tenvex/manual.

#include <tenvex.h>
#include <cstddef>

using namespace tnvx;

__attribute__((noinline))
void loop_et(const vec4* va, const vec4* vb, const vec4* vc, vf4* out, std::size_t n)
{
    for (std::size_t i = 0; i < n; ++i)
        out[i] =
            (norm3(va[i] + vb[i] * 2.0f) *
             dot3(vb[i], vc[i]) +
             vc[i] * 3.0f)
            .eval();
}

g++ -std=c++20 -O3 -I./tenvex/src/tenvex -S sample4.cpp -o sample5.asm

_Z7loop_etPKN4tnvx4vec4ES2_S2_P13__Float32x4_tm:
.LFB4063 :
.cfi_startproc
    cbz x4, .L7
    sub sp, sp, #400
    .cfi_def_cfa_offset 400
    mov x5, 0
    fmov v31.4s, 2.0e+0
    add x15, sp, 128
    fmov v30.4s, 3.0e+0
    mov x14, sp
    add x13, sp, 352
    add x8, sp, 160
    add x7, sp, 256
    add x6, sp, 208
    add x12, sp, 96
    add x11, sp, 32
    add x10, sp, 64
    mov x9, 0
    .p2align 5, , 15
.L2:
    ldr q29, [x1, x5]
    adrp x16, .LC0
    ldr q27, [x0, x5]
    add x9, x9, 1
    ldr q0, [x2, x5]
    mov v1.16b, v29.16b
    stp q27, q29, [x6]
    fmul v26.4s, v0.4s, v30.4s
    stp q27, q29, [x7]
    movi v23.4s, 0x40, lsl 24
    stp q27, q29, [x8]
    str q29, [x11]
    str q0, [x10]
    stp q29, q0, [x12]
    fmul v29.4s, v0.4s, v29.4s
    str q27, [x13]
    fmla v27.4s, v31.4s, v1.4s
    str q1, [x14]
    ins v29.s[3], wzr
    faddp v29.4s, v29.4s, v29.4s
    str q1, [x15]
    fmul v28.4s, v27.4s, v27.4s
    stp q0, q30, [sp, 304]
    faddp v29.4s, v29.4s, v29.4s
    str q31, [sp, 336]
    ins v28.s[3], wzr
    faddp v28.4s, v28.4s, v28.4s
    stp q1, q23, [x13, 16]
    ldr q22, [x16, #:lo12:.LC0]
    str q0, [x11, 16]
    faddp v28.4s, v28.4s, v28.4s
    str q23, [x14, 16]
    str q22, [x10, 16]
    dup v28.4s, v28.s[0]
    str q23, [x15, 16]
    str q23, [x6, 32]
    fsqrt v28.4s, v28.4s
    str q23, [x7, 32]
    str q23, [x8, 32]
    fdiv v28.4s, v27.4s, v28.4s
    ins v28.s[3], v27.s[3]
    fmla v26.4s, v28.4s, v29.s[0]
    str q26, [x3, x5]
    add x5, x5, 16
    cmp x4, x9
    bne .L2
    add sp, sp, 400
    .cfi_def_cfa_offset 0
    ret
.L7:
ret
.cfi_endproc

Listing 8.

A too-clean toy repro can confirm a false hypothesis, precisely because it's clean. It doesn't contain the difference from reality where the bug lives. The toy confirmed "by-value saves", but reality responded with doubled regression.

Method Reversal: A Repro That Didn't Reproduce

After the by-value failure, I went the other way: building a minimal repro from the bottom up. From a clean toy, adding real tenvex traits one by one. Added tnvx_ref_or_value_t. Added empty CRTP base Expr<Derived>. Added nested aggregates. Added layout with a vector member. The bug didn't appear at any step - not from any single trait, nor from their combination. The minimal repro remained sterile no matter what I added.

(I decided not to insert the listing of these bug-hunting attempts here - it takes over 200 lines of code)

And this was the turning point. Not because it found the bug, but because it told me: the method is wrong. Building up a clean toy, you can't catch a bug that lives in the difference between the toy and reality, because that's exactly the difference I was discarding each time. The right move was the opposite: search downward from the real spilling translation unit, replacing real components with toy ones one by one, until the frame disappears. The frame disappearing at some step would then point to the culprit.

When a minimal repro doesn't reproduce the bug, stop building clean. Search downward from the real spilling TU, not upward from a toy.

The Reveal: It Was a Benchmark Bug

When I finally started searching downward, the cause turned out to be offensively mundane. It wasn't in the library at all. It was in the benchmark: the compound expression used the same array elements - vb[i], vc[i] - twice, right inside the ET expression. Repeated array access by index in the role of a leaf.

On GCC/aarch64 this breaks SROA: the optimizer can't prove that the materialized node's address doesn't escape, and just in case leaves the node in memory, correctness trumps speed. Hence the dead ~400-byte frame and those stp ..., [sp] that no one reads. And why did toys never reproduce this? Because they by coincidence always bound leaves to local variables - I intuitively wrote this in toys, but didn't do it in benchmarks. Constructors like Leaf(va[i]) first fetched the element into a local variable. The very thing that turned out to be the fix, toys did by accident. That's why too-clean confirmed the false hypothesis time after time.

The fix is one line in user code: bind leaves to local variables before the expression (vec4 a = va[i], b = vb[i], c = vc[i];), and write the expression over a, b, c. The dead frame collapses to zero. I confirmed this in five reduction steps: literally the second step with leaf binding removes the frame. It turned out the same artifact had infected other benchmarks too. After fixing, 6 out of 6 cells were zero-cost, within measurement noise.

static void BM_Compound_tenvex(benchmark::State& state)
{
    const auto va = make_vecs(1024, 1);
    const auto vb = make_vecs(1024, 2);
    const auto vc = make_vecs(1024, 3);
    std::size_t i = 0;
    for (auto _ : state)
    {
        vec4 r =
            norm3(va[i] + vb[i] * 2.0f) *
            dot3(vb[i], vc[i]) +
            vc[i] * 3.0f;
        benchmark::DoNotOptimize(r);
        i = (i + 1) & 1023;
    }
}

static void BM_Compound_tenvex(benchmark::State& state)
{
    const auto va = make_vecs(1024, 1);
    const auto vb = make_vecs(1024, 2);
    const auto vc = make_vecs(1024, 3);
    std::size_t i = 0;

    for (auto _ : state)
    {
        vec4 a = va[i], b = vb[i], c = vc[i];

        vec4 r =
            norm3(a + b * 2.0f) *
            dot3(b, c) +
            c * 3.0f;

        benchmark::DoNotOptimize(r);
        i = (i + 1) & 1023;
    }
}

static void BM_Lerp_Throughput(benchmark::State& state)
{
    const auto a = make_vecs(1024, 1);
    const auto b = make_vecs(1024, 2);
    Scalar t { 0.35f };
    std::size_t i = 0;

    for (auto _ : state)
    {
        vec4 r = lerp(a[i], b[i], t);
        benchmark::DoNotOptimize(r);
        i = (i + 1) & 1023;
    }
}

static void BM_Lerp_Throughput(benchmark::State& state)
{
    const auto a = make_vecs(1024, 1);
    const auto b = make_vecs(1024, 2);
    Scalar t { 0.35f };
    std::size_t i = 0;

    for (auto _ : state)
    {
        vec4 av = a[i], bv = b[i];
        vec4 r = lerp(av, bv, t);
        benchmark::DoNotOptimize(r);
        i = (i + 1) & 1023;
    }
}

static void BM_Reflect_Throughput(benchmark::State& state)
{
    const auto v = make_vecs(1024, 1);
    const auto n = make_vecs(1024, 2);
    std::size_t i = 0;
    for (auto _ : state)
    {
        vec4 r = reflect(v[i], n[i]);
        benchmark::DoNotOptimize(r);
        i = (i + 1) & 1023;
    }
}

static void BM_Reflect_Throughput(benchmark::State& state)
{
    const auto v = make_vecs(1024, 1);
    const auto n = make_vecs(1024, 2);
    std::size_t i = 0;

    for (auto _ : state)
    {
        vec4 vv = v[i], nv = n[i];
        vec4 r = reflect(vv, nv);
        benchmark::DoNotOptimize(r);
        i = (i + 1) & 1023;
    }
}

Listing 9.

And here's the same calculation when leaves are bound to local variables. Compare line by line with the previous listing: the same fmla on a + b*2, the same fmuls on products, the same horizontal reduction (ins into w-lane plus a pair of faddp), fsqrt, fdiv, final fmla. The computations match instruction to instruction - it's the same norm, the same dot, the same c*3. The difference is in what's not here. There's no sub sp: the function doesn't open a frame at all. Memory is touched at exactly four places: three ldrs on entry (va[i], vb[i], vc[i]) and one str on exit (out[i]), i.e., only reading arguments and writing the result. Intermediate values don't leave registers: not a single str to [sp], not a single address computed for spilling. The entire dead frame rigging from loop_et - sub/add sp, three service addresses, six dead writes - is simply absent. This is the proof: the extra frame is not a payment for the math itself, that matches byte for byte, but an artifact of the form in which the math is written. Bind array elements to local variables, and each address becomes provably non-escaping, SROA breaks the expression into registers, the frame collapses to zero. That very one-line fix.

#include <tenvex.h>
#include <cstddef>

using namespace tnvx;

__attribute__((noinline))
void loop_manual(const vec4* va, const vec4* vb, const vec4* vc, vf4* out, std::size_t n)
{
    using namespace tnvx::detail;

    for (std::size_t i = 0; i < n; ++i)
    {
        vf4 a = va[i].eval(),
             b = vb[i].eval(),
             c = vc[i].eval();

        vf4 t = norm3(add(a, mul(b, set_all(2.0f))));
        vf4 d = dot3(b, c);
        out[i] = add(mul(t, d), mul(c, set_all(3.0f)));
    }
}

g++ -std=c++20 -O3 -I./tenvex/src/tenvex -S sample6.cpp -o sample6.asm

_Z11loop_manualPKN4tnvx4vec4ES2_S2_P13__Float32x4_tm :
.LFB4064 :
.cfi_startproc
    cbz x4, .L11
    fmov v31.4s, 2.0e+0
    mov x5, 0
    fmov v30.4s, 3.0e+0
    mov x6, 0
    .p2align 5, , 15
.L13:
    ldr q29, [x1, x5]
    add x6, x6, 1
    ldr q27, [x0, x5]
    ldr q28, [x2, x5]
    fmla v27.4s, v29.4s, v31.4s
    fmul v29.4s, v28.4s, v29.4s
    fmul v28.4s, v28.4s, v30.4s
    fmul v26.4s, v27.4s, v27.4s
    ins v29.s[3], wzr
    faddp v29.4s, v29.4s, v29.4s
    ins v26.s[3], wzr
    faddp v26.4s, v26.4s, v26.4s
    faddp v29.4s, v29.4s, v29.4s
    faddp v26.4s, v26.4s, v26.4s
    dup v26.4s, v26.s[0]
    fsqrt v26.4s, v26.4s
    fdiv v26.4s, v27.4s, v26.4s
    ins v26.s[3], v27.s[3]
    fmla v28.4s, v26.4s, v29.s[0]
    str q28, [x3, x5]
    add x5, x5, 16
    cmp x4, x6
    bne .L13
.L11:
    ret
.cfi_endproc

Listing 10

Table 4: median values over 30 benchmarks, including throughput benchmarks where the next computation didn't depend on the previous one, latency benchmarks where each computation depended on the previous one. manual - hand-built from the same kernels that tenvex uses.

A microbenchmark must mirror the real access pattern. A too-clean benchmark is not neutral, it confirms a lie: my bench measured an artifact that doesn't exist in normal code.

Six out of six cells - zero-cost. The only one that stood out was held not by library properties, but by one line in the bench: repeated vb[i] inside the expression. Remove the repetition, and tenvex and hand-written code converge (within measurement noise). On GCC/arm64 this is literally the same machine code. It's worth spelling out what this means and what it doesn't. The zero-cost promise held, the ET abstraction costs nothing at runtime relative to the same expression hand-built from the same kernels. But "costs nothing" here is about abstraction overhead, not that tenvex is the fastest way to compute this expression in principle. All this time we were answering one question: does the abstraction leak? Answer: no, it doesn't (and where it almost did, the bench was at fault). The question "and where does it actually win, and by how much?" is a different one, and we'll get to it soon.

Where the Dead Frame Comes From

Now that the cause is clear, the theory falls into place cleanly. A local variable's address is provably non-escaping; an array element's - is. When the same element appears twice as a leaf, GCC/aarch64's escape analysis can't prove that the materialized node doesn't escape, and leaves it in memory. And the DSE pass on this backend doesn't squeeze the resulting dead writes, hence the frame. Bind the leaf to a local variable, the address becomes provably non-escaping, SROA breaks down the node, and there's no frame.

It remains to show the flip side - where the abstraction not only doesn't leak, but actually wins. To that we now turn.

Part III. Where We Win and Where We Don't

We've established that the abstraction doesn't leak. The reverse question remains: what does it actually win? There are two competitors: you can win against a naive implementation - against a person who wrote the same thing straightforwardly; you can measure yourself against the compiler - how far it will stretch naive code with optimizations on its own. I warn you right away, there will be no parade of victories. In some places tenvex loses, and fairly. Those places will be the most interesting.

Summary and How to Read It

First, raw numbers: naive scalar implementation, tenvex, and their ratio by operation.

Table 5.

One caveat about reading. Numbers vary noticeably from compiler to compiler, this is normal and interesting in its own way, but it's a separate story, and I won't get into it here. Latency and throughput are numbers about different things: critical path length versus throughput capacity.

Three Classes of Operations

Operations fall into three classes, and each has its own story.

The first is element-wise, which the compiler vectorizes on its own. Naive code and tenvex compile to the same thing. Zero win, the work was already done for us, nothing to compete with.

#include <tenvex.h>
#include <naive_vec4.h>
#include <cstddef>

template <class T>
static inline void DoNotOptimize(T const& v)
{
    asm volatile("" : : "x,r,m"(v) : "memory");
}

__attribute__((noinline))
void LerpNaive(const naive::vec4* a, const naive::vec4* b, float t, std::size_t reps)
{
    std::size_t i = 0;

    for (std::size_t k = 0; k < reps; ++k)
    {
        naive::vec4 r = naive::lerp(a[i], b[i], t);
        DoNotOptimize(r);
        i = (i + 1) & 1023;
    }
}

__attribute__((noinline))
void LerpTnvx(const tnvx::vec4* a, const tnvx::vec4* b, float t, std::size_t reps)
{
    std::size_t i = 0;

    for (std::size_t k = 0; k < reps; ++k)
    {
        tnvx::vec4 r = tnvx::lerp(a[i], b[i], t);
        DoNotOptimize(r);
        i = (i + 1) & 1023;
    }
}

g++ -std=c++20 -O3 -msse4.1 -mfma -S -masm=intel -fno-asynchronous-unwind-tables -I./tenvex/src/naive -I./tenvex/src/tenvex sample_lerp.cpp -o sample_lerp.asm

_Z9LerpNaivePKN5naive4vec4ES2_fm :
.LFB2664 :
.cfi_startproc
    endbr64
    mov r8, rsi
    mov rsi, rdx
    test rdx, rdx
    je .L8
    xor edx, edx
    xor eax, eax
    vshufps xmm0, xmm0, xmm0, 0
    .p2align 4
    .p2align 3
.L3:
    mov rcx, rax
    sal rcx, 4
    vmovups xmm2, [rdi + rcx]
    vmovups xmm1, [r8 + rcx]
    vsubps xmm1, xmm1, xmm2
    vfmadd132ps xmm1, xmm2, xmm0
    add rax, 1
    add rdx, 1
    and eax, 1023
    cmp rsi, rdx
    jne .L3
.L8:
    ret
.cfi_endproc

Z8LerpTnvxPKN4tnvx4vec4ES2_fm :
.LFB2665 :
.cfi_startproc
    endbr64
    mov r8, rdi
    mov rdi, rsi
    mov rsi, rdx
    xor eax, eax
    xor edx, edx
    vshufps xmm0, xmm0, xmm0, 0
    test rsi, rsi
    je .L16
.L12:
    mov rcx, rax
    sal rcx, 4
    vmovaps xmm2, [r8 + rcx]
    vmovaps xmm1, [rdi + rcx]
    vsubps xmm1, xmm1, xmm2
    vfmadd132ps xmm1, xmm2, xmm0
    add rax, 1
    add rdx, 1
    and eax, 1023
    cmp rsi, rdx
    jne .L12
.L16:
    ret
.cfi_endproc

Listing 11: disassembly for tenvex and naive implementations is equivalent in the computational kernel; differences reduce to register allocation chosen by the compiler, and using vmovaps instead of vmovups (tnvx::vec4 is 16-byte aligned).

Second: element-wise where the auto-vectorizer gives up, this is a win. For saturate (clamp to [0,1] - min(max(v,0),1)) both versions turn out vector, but the naive one spends an extra instruction: the compiler folds std::min/std::max into cmp+blendvps+minps, while tenvex gets by with cmp+andps+minps. Same for floor/ceil/frac - tenvex takes one packed operation and wins consistently (frac ~1.5×, reflect ~1.3×). The advantage is exactly where the compiler couldn't fold as tightly.

#include <tenvex.h>
#include <naive_vec4.h>

naive::vec4 TestSaturateNaive(naive::vec4 v)
{
    return naive::saturate(v);
}

tnvx::vec4 TestSaturateTnvx(tnvx::vec4 v)
{
    return tnvx::saturate(v);
}

g++ -std=c++20 -O3 -msse4.1 -mfma -S -masm=intel -fno-asynchronous-unwind-tables -I./tenvex/src/naive -I./tenvex/src/tenvex sample_saturate.cpp -o sample_saturate.asm

_Z17TestSaturateNaiveN5naive4vec4E :
.LFB2651 :
.cfi_startproc
    endbr64
    vpunpcklqdq xmm1, xmm0, xmm1
    vxorps xmm0, xmm0, xmm0
    vcmpltps xmm2, xmm1, xmm0
    vblendvps xmm1, xmm1, xmm0, xmm2
    vbroadcastss xmm0, DWORD PTR .LC1[rip]
    vminps xmm3, xmm0, xmm1
    vmovaps XMMWORD PTR -24[rsp], xmm3
    vmovq xmm1, QWORD PTR -16[rsp]
    vmovq xmm0, QWORD PTR -24[rsp]
    ret
.cfi_endproc

_Z16TestSaturateTnvxN4tnvx4vec4E :
.LFB2652 :
.cfi_startproc
    endbr64
    vmovaps xmm1, xmm0
    vxorps xmm0, xmm0, xmm0
    vcmpltps xmm0, xmm0, xmm1
    vandps xmm0, xmm0, xmm1
    vbroadcastss xmm1, DWORD PTR .LC1[rip]
    vminps xmm0, xmm0, xmm1
    ret
.cfi_endproc

Listing 12.

The issue here is std::min/std::max and NaN handling. std::min(a, b) is defined as b < a ? b : a, on equality and on NaN it has its own operand return semantics. And minps/maxps behave differently on NaN: with NaN in any operand they return the second operand, and the argument order is baked into the instruction itself. That is, minps is not a bitwise equivalent of four std::mins: on NaN inputs they'll give different results. The compiler is obliged to preserve observable program behavior. Since minps can diverge from std::min on NaN, it has no right to replace one with the other and, lacking a guarantee that NaN won't come, plays it safe and leaves scalar code. This is not optimizer stupidity, but its honesty - it doesn't know NaN won't be in the data.

Third, slightly subtler. Several latency-linked cells on ARM where tenvex loses for real. AABB by latency: naive min/max the compiler folds into fcmgt+bit (compare plus bitselect), while tenvex takes fmin/fmax, and on Cortex-A76 in a dependent chain they're more expensive due to NaN semantics. Where each iteration's result waits for the previous one, the cheaper-by-latency compare-select of the naive wins. On x86_64 this doesn't happen, there min/max instructions are cheap, and the same operations are ahead for tenvex. This is the price of microarchitecture.

#include <tenvex.h>
#include <naive_vec4.h>
#include <cstddef>

template <class T>
static inline void DoNotOptimize(T const& v)
{
    asm volatile("" : : "x,r,m"(v) : "memory");
}

__attribute__((noinline))
void AABBNaive(const naive::vec4* pts, std::size_t size)
{
    naive::vec4 lo { 1e30f,  1e30f,  1e30f, 0.0f };
    naive::vec4 hi { -1e30f, -1e30f, -1e30f, 0.0f };

    for (std::size_t i = 0; i < size; ++i)
    {
        naive::vec4 vertex = pts[i];
        lo = naive::min(lo, vertex);
        hi = naive::max(hi, vertex);
    }

    DoNotOptimize(lo);
    DoNotOptimize(hi);
}

__attribute__((noinline))
void AABBTnvx(const tnvx::vec4* pts, std::size_t size)
{
    tnvx::vec4 lo { 1e30f,  1e30f,  1e30f, 0.0f };
    tnvx::vec4 hi { -1e30f, -1e30f, -1e30f, 0.0f };

    for (std::size_t i = 0; i < size; ++i)
    {
        tnvx::vec4 vertex = pts[i];
        lo = tnvx::min(lo, vertex);
        hi = tnvx::max(hi, vertex);
    }

    DoNotOptimize(lo);
    DoNotOptimize(hi);
}

g++ -std=c++20 -O3 -msse4.1 -mfma -S -masm=intel -fno-asynchronous-unwind-tables -I./tenvex/src/naive -I./tenvex/src/tenvex sample_aabb.cpp -o sample_aabb.asm

_Z9AABBNaivePKN5naive4vec4Em :
.LFB2664 :
.cfi_startproc
    endbr64
    test rsi, rsi
    je .L4
    sal rsi, 4
    vmovaps xmm1, XMMWORD PTR .LC0[rip]
    vmovaps xmm2, XMMWORD PTR .LC1[rip]
    add rsi, rdi

    .p2align 5
    .p2align 4
    .p2align 3

.L3:
    vmovups xmm0, [rdi]
    add rdi, 16
    vmaxps xmm1, xmm0, xmm1
    vminps xmm2, xmm0, xmm2
    cmp rdi, rsi
    jne .L3
    ret

.L4:
    vmovaps xmm1, XMMWORD PTR .LC0[rip]
    vmovaps xmm2, XMMWORD PTR .LC1[rip]
    ret
.cfi_endproc

_Z8AABBTnvxPKN4tnvx4vec4Em :
.LFB2665 :
.cfi_startproc
    endbr64
    test rsi, rsi
    je .L10
    sal rsi, 4
    vmovaps xmm2, XMMWORD PTR .LC1[rip]
    vmovaps xmm1, XMMWORD PTR .LC0[rip]
    add rsi, rdi

    .p2align 5
    .p2align 4
    .p2align 3

.L9:
    vmovaps xmm0, [rdi]
    add rdi, 16
    vminps xmm2, xmm2, xmm0
    vmaxps xmm1, xmm1, xmm0
    cmp rsi, rdi
    jne .L9
    ret

.L10:
    vmovaps xmm2, XMMWORD PTR .LC1[rip]
    vmovaps xmm1, XMMWORD PTR .LC0[rip]
    ret
.cfi_endproc

g++ -std=c++20 -O3 -S -fno-asynchronous-unwind-tables -I./tenvex/src/naive -I./tenvex/src/tenvex sample_aabb.cpp -o sample_aabb.asm

_Z9AABBNaivePKN5naive4vec4Em :
.LFB5938 :
.cfi_startproc
    cbz x1, .L4
    adrp x2, .LC0
    add x1, x0, x1, lsl 4
    ldr q7, [x2, #:lo12:.LC0]
    adrp x2, .LC1
    ldr q6, [x2, #:lo12:.LC1]

    .p2align 5

.L3:
    ldr q5, [x0], 16
    fcmgt v0.4s, v5.4s, v7.4s
    fcmgt v4.4s, v6.4s, v5.4s
    bit v7.16b, v5.16b, v0.16b
    bit v6.16b, v5.16b, v4.16b
    cmp x0, x1
    bne .L3
    ret

.L4:
    adrp x0, .LC0
    ldr q7, [x0, #:lo12:.LC0]
    adrp x0, .LC1
    ldr q6, [x0, #:lo12:.LC1]
    ret
.cfi_endproc

_Z8AABBTnvxPKN4tnvx4vec4Em :
.LFB5939 :
.cfi_startproc
    cbz x1, .L10
    adrp x2, .LC0
    add x1, x0, x1, lsl 4
    ldr q7, [x2, #:lo12:.LC0]
    adrp x2, .LC1
    ldr q6, [x2, #:lo12:.LC1]

    .p2align 5

.L9:
    ldr q31, [x0], 16
    fmin v6.4s, v6.4s, v31.4s
    fmax v7.4s, v7.4s, v31.4s
    cmp x0, x1
    bne .L9
    ret

.L10:
    adrp x0, .LC0
    ldr q7, [x0, #:lo12:.LC0]
    adrp x0, .LC1
    ldr q6, [x0, #:lo12:.LC1]
    ret
.cfi_endproc

Listing 13.

It would seem tenvex should win on arm64 - shorter by instruction count, two against four. But it loses by latency, because what matters is not quantity, but latency on the critical path.

Horizontal Operations on Wide Registers Are Expensive by Latency

The most interesting case, dot3/dot4. By latency dot loses everywhere: on ARM 0.44× for dot3, 0.53× for dot4. This is the honest price of horizontal reduction (an operation that collapses 4 numbers in a wide register into one) on the critical path, and the reason is visible in asm: long shuffle-chain against short FMA-scalar. Horizontal reduction is expensive by latency. In throughput the cost of reduction partially hides behind overlapping independent iterations: the critical path is no longer visible, port occupancy matters. But hides doesn't mean disappears: where tenvex has port headroom, the gap flips (dot4 on ARM pulls ahead, up to 1.81×), and where shuffles are still in the bottleneck, it stays around parity (dot3 - weak plus, balancing at zero). The confident flip happens where ports suffice, at dot4.

#include <tenvex.h>
#include <naive_vec4.h>

float TestDot3Naive(naive::vec4 a, naive::vec4 b)
{
    return naive::dot3(a, b);
}

float TestDot3Tnvx(tnvx::vec4 a, tnvx::vec4 b)
{
    return tnvx::dot3(a, b);
}

float TestDot4Naive(naive::vec4 a, naive::vec4 b)
{
    return naive::dot4(a, b);
}

float TestDot4Tnvx(tnvx::vec4 a, tnvx::vec4 b)
{
    return tnvx::dot4(a, b);
}

g++ -std=c++20 -O3 -msse4.1 -mfma -S -masm=intel -fno-asynchronous-unwind-tables -I./tenvex/src/naive -I./tenvex/src/tenvex sample_dots.cpp -o sample_dots.asm

_Z13TestDot3NaiveN5naive4vec4ES0_ :
.LFB2651 :
.cfi_startproc
    endbr64
    vmovdqa xmm5, xmm2
    vmovdqa xmm4, xmm0
    vshufps xmm2, xmm2, xmm2, 85
    vshufps xmm0, xmm0, xmm0, 85
    vmulss xmm2, xmm0, xmm2
    vfmadd132ss xmm4, xmm2, xmm5
    vfmadd132ss xmm1, xmm4, xmm3
    vmovaps xmm0, xmm1
    ret
.cfi_endproc

_Z12TestDot3TnvxN4tnvx4vec4ES0_ :
.LFB2652 :
.cfi_startproc
    endbr64
    vmulps xmm1, xmm1, xmm0
    vpermilps xmm0, xmm1, 85
    vaddss xmm0, xmm1, xmm0
    vmovhlps xmm1, xmm1, xmm1
    vaddss xmm0, xmm0, xmm1
    ret
.cfi_endproc

Z13TestDot4NaiveN5naive4vec4ES0_ :
.LFB2653 :
.cfi_startproc
    endbr64
    vmovdqa xmm7, xmm2
    vmovdqa xmm5, xmm0
    vshufps xmm2, xmm2, xmm2, 85
    vmovdqa xmm6, xmm3
    vshufps xmm0, xmm0, xmm0, 85
    vmulss xmm0, xmm0, xmm2
    vmovdqa xmm4, xmm1
    vshufps xmm3, xmm3, xmm3, 85
    vshufps xmm1, xmm1, xmm1, 85
    vfmadd132ss xmm5, xmm0, xmm7
    vmovaps xmm0, xmm3
    vfmadd132ss xmm4, xmm5, xmm6
    vfmadd132ss xmm0, xmm4, xmm1
    ret
.cfi_endproc

_Z12TestDot4TnvxN4tnvx4vec4ES0_ :
.LFB2654 :
.cfi_startproc
    endbr64
    vmulps xmm0, xmm1, xmm0
    vhaddps xmm0, xmm0, xmm0
    vhaddps xmm0, xmm0, xmm0
    ret
.cfi_endproc

Listing 14.

And again fewer instructions, but higher latency. In latency what matters is critical path length. Naive dot3: the compiler pulls up components with a pair of shufps, then mul and two merged fmas - short tail on FMA ports. SIMD-dot3 drags horizontal reduction: mul, then shuffle to merge lanes, then dependent additions. The chain is longer, and two shuffles sit on the latent shuffle-port. In throughput, when there are many dots and they're independent, the critical path hides behind iteration overlap: one vector load plus reduction at tenvex overlaps better than scattered scalar reads of the naive, but shuffles are still occupied, and dot3 in throughput only pulls up to parity, not ahead. Confident flip happens where ports suffice, at dot4.

How do grown-up libraries live with this? The same way. DirectXMath does the same horizontal reduction (mul → shuffle → add → shuffle → add → splat in baseline-SSE), branching by available instructions: scalar without intrinsics, vadd+vpadd on NEON, dpps on SSE4.1, two hadds on SSE3, shuffle-reduction on baseline SSE. GLM - the same for vec4: dpps on SSE4.1, two hadds on SSE3, shuffle+add on baseline. Both have several implementations for different instruction sets, and the very fact of these #if trees proves there's no single fast way for a horizontal operation.

inline XMVECTOR XM_CALLCONV XMVector4Dot
(
    FXMVECTOR V1,
    FXMVECTOR V2
) noexcept
{
#if defined(_XM_NO_INTRINSICS_)
    XMVECTORF32 Result;
    Result.f[0] =
        Result.f[1] =
        Result.f[2] =
        Result.f[3] =
        V1.vector4_f32[0] * V2.vector4_f32[0] +
        V1.vector4_f32[1] * V2.vector4_f32[1] +
        V1.vector4_f32[2] * V2.vector4_f32[2] +
        V1.vector4_f32[3] * V2.vector4_f32[3];
    return Result.v;

#elif defined(_XM_ARM_NEON_INTRINSICS_)
    float32x4_t vTemp = vmulq_f32(V1, V2);
    float32x2_t v1 = vget_low_f32(vTemp);
    float32x2_t v2 = vget_high_f32(vTemp);
    v1 = vadd_f32(v1, v2);
    v1 = vpadd_f32(v1, v1);
    return vcombine_f32(v1, v1);

#elif defined(_XM_SSE4_INTRINSICS_)
    return _mm_dp_ps(V1, V2, 0xff);

#elif defined(_XM_SSE3_INTRINSICS_)
    XMVECTOR vTemp = _mm_mul_ps(V1, V2);
    vTemp = _mm_hadd_ps(vTemp, vTemp);
    return _mm_hadd_ps(vTemp, vTemp);

#elif defined(_XM_SSE_INTRINSICS_)
    XMVECTOR vTemp2 = V2;
    XMVECTOR vTemp = _mm_mul_ps(V1, vTemp2);
    vTemp2 = _mm_shuffle_ps(vTemp2, vTemp, _MM_SHUFFLE(1, 0, 0, 0));
    vTemp2 = _mm_add_ps(vTemp2, vTemp);
    vTemp = _mm_shuffle_ps(vTemp, vTemp2, _MM_SHUFFLE(0, 3, 0, 0));
    vTemp = _mm_add_ps(vTemp, vTemp2);
    return XM_PERMUTE_PS(vTemp, _MM_SHUFFLE(2, 2, 2, 2));
#endif
}

GLM_FUNC_QUALIFIER glm_vec4 glm_vec4_dot(glm_vec4 v1, glm_vec4 v2)
{
#       if GLM_ARCH & GLM_ARCH_SSE41_BIT
    return _mm_dp_ps(v1, v2, 0xff);
#       elif GLM_ARCH & GLM_ARCH_SSE3_BIT
    glm_vec4 const mul0 = _mm_mul_ps(v1, v2);
    glm_vec4 const hadd0 = _mm_hadd_ps(mul0, mul0);
    glm_vec4 const hadd1 = _mm_hadd_ps(hadd0, hadd0);
    return hadd1;
#       else
    glm_vec4 const mul0 = _mm_mul_ps(v1, v2);
    glm_vec4 const swp0 = _mm_shuffle_ps(mul0, mul0, _MM_SHUFFLE(2, 3, 0, 1));
    glm_vec4 const add0 = _mm_add_ps(mul0, swp0);
    glm_vec4 const swp1 = _mm_shuffle_ps(add0, add0, _MM_SHUFFLE(0, 1, 2, 3));
    glm_vec4 const add1 = _mm_add_ps(add0, swp1);
    return add1;
#       endif
}

Listing 15: dot implementations in DirectXMath (left) and GLM (right)

Honest conclusion about latent dot: there's not much room to optimize it, horizontal operations are baked into the task of "one vec4 → one scalar." The way out: if there are many dots and latency specifically matters, then you need to change not the operation, but the data layout. SoA (Structure of Arrays): instead of an array of vectors, hold a structure of arrays, and then horizontal becomes vertical, dot is computed as x*x + y*y + z*z by lanes, without a single shuffle. tenvex doesn't have SoA yet; this is a direction that logic leads to and where the industry has long gone (DXM streaming functions, SoA in gamedev), not a result I measured.

Then Why ET?

Here's a question that suggests itself. DirectXMath gets by with ordinary inline functions, no expression templates whatsoever. After inlining, the compiler sees the whole expression anyway. Why then all this ET machinery?

For the speed of the optimizer-generated code itself - not needed, and this should be said straight. After inlining, the optimizer sees the same instruction stream whether from a chain of functions or from a tree of nodes, and folds intermediates the same way. Remember loop_manual? - flat kernel calls without any ET? It compiled no worse than the ET version. The same expression ET won't turn into faster machine code than well-inlined functions.

But there's a second lever, and this one is about real optimizations. ET can change not the code quality, but the amount of computation itself, compute not the same thing faster, but simply less. A function is obliged to return its full result, it doesn't know how it'll be used. An ET node knows context, it's part of a tree, and the top eval() can rewrite the subtree. magnitude3(a) < magnitude3(b) - both sqrts are unnecessary, comparing squares suffices: and this is unavailable to the compiler from ordinary functions, because it requires algebraic knowledge that sqrt is monotonic, not just dead code elimination. There are also simpler cases, laziness: took only .x() from cross3 - the other components can be skipped (though a good compiler sometimes cleans this itself via DCE). Plus ergonomics: a + b*2 instead of XMVectorAdd(a, XMVectorScale(b, 2)). And ergonomics is not nothing. We all know from experience that code that reads well is easier to debug (which is exactly why we sweep the complexity of debugging expression templates under the rug).

namespace tnvx
{

// Generic case: both sides collapse to float, compare as usual.
template<vec_expr L, vec_expr R>
[[nodiscard]] TNVX_INLINE
bool operator<(const L& l, const R& r) noexcept
{
    return static_cast<float>(l) < static_cast<float>(r);
}

// Rewrite case: |a| < |b|  ->  |a|^2 < |b|^2  (no sqrt on either side)
template<vec_expr A, vec_expr B>
[[nodiscard]] TNVX_INLINE
bool operator<(const Magn3<A>& a, const Magn3<B>& b) noexcept
{
    // peel the argument back out of each Magn3 node and re-wrap as Magn3Sq
    return static_cast<float>(magnitude3_sq(a.arg())) <
           static_cast<float>(magnitude3_sq(b.arg()));
}

}

Listing 16: Sketch: how one expression-level rewrite might look. tenvex doesn't have this yet, magnitude3_sq is chosen by hand now. Shown to make ET potential tangible, not abstract.

Conclusion

I started with a half-joke about the fastest linear algebra library and an honest question whether zero-cost abstraction really costs nothing. The answer exists, and it's conditional. Zero-cost is a property not of the code and not of the compiler separately, but of their pair: the promise holds when inlining, SROA, promotion to registers, and dead store elimination all come together. If they come together, then an ET expression compiles to the same code as hand-written intrinsics. If even one link fails, the abstraction leaks. But this doesn't mean the abstraction isn't worth making. It leaks rarely and predictably, and fixing it is usually simple, one line (though finding the problem may not be entirely trivial). Writing intrinsics by hand across all code for the sake of that one percent where it would leak is clearly a bad trade. On all six backends, when the bench finally became honest, the promise held.

The main lesson is about an honest bench. For most of the journey, my opponent was not the compiler, but my own measurements. Those very +39% turned out to be not an abstraction leak, but a dead frame from repeated index access in my own benchmark, and this bug sat in half the bench, giving imaginary regressions where there were none. Before blaming the backend, check what you're measuring, and look at the disassembly, not at the model in your head of what's happening. Writing a good benchmark is a whole art.

What's also important is what zero-cost doesn't touch: it never promised a lead over naive code. Where the auto-vectorizer gave up, the wide operation wins; where the compiler is strong on its own - draw. And somewhere you honestly lose: horizontal operations are expensive by latency, fmin is expensive in a dependent chain. And separately about expression templates, for which all this was started - for machine code speed they're not needed, after inlining the compiler folds the same expression the same way. Their real niche is to do less work, skip an extra sqrt, not compute an unnecessary component, but that's potential so far, not implemented. So today ET gives ergonomics, and speed lives in registers, in data layout, and in whether the compiler had enough guarantees. Was all this fuss worth one vector? By the original count, yes. I wanted to get hands-on with SIMD, play with optimizations, and test the zero-cost promise for strength, and the understanding of what speed is made of turned out more honest and interesting.