Чего стоят zero-cost-абстракции

Пишу самую быструю линейную алгебру, считаю такты, гоняюсь за регрессией

Артём Зыктин

репозиторий tenvex

TL;DR. Написал маленькую SIMD-библиотеку для vec4 и с дизассемблером проверил, правда ли C++ zero-cost абстракции ничего не стоят. Ответ: байт-в-байт как ручные интринсики. Потом гонялся за регрессией в +39% на arm64, был уверен, что поймал протечку абстракции, перебрал несколько гипотез, и каждую разбил об дизассемблер. В чем было дело, разобрался, но не сразу и не там, где искал. А заодно выяснилось, что "ничего не стоит" и "обгоняет наивный код" - это две разные вещи: где-то компилятор нас побеждает, а где-то мы его. Про регистры, дизассемблер и то, как красивая гипотеза проигрывает фактам.

Однажды, засидевшись ночью, я решил: а почему бы мне не написать самую быструю библиотеку для линейной алгебры. Шучу, конечно, это слишком амбициозная задача для одного разработчика. Есть целые команды, которые бьются над ней десятилетиями - DirectXMath, GLM, Eigen, - вряд ли я смогу с ними тягаться (ну, это мы еще посмотрим). На самом деле эта мысль пришла мне в голову по нескольким куда более приземленным причинам.

Во-первых, это отличный повод напомнить себе линейную алгебру. Особенно в контексте компьютерной графики, где у меня всегда были пробелы. Во-вторых, я никогда не имел дела вплотную с SIMD-операциями и широкими регистрами (SSE, AVX и так далее), а это отличная тема, чтобы пощупать такие вещи руками. В-третьих, мне хотелось поиграть с оптимизациями. Ну и шаблоны: я вообще люблю C++ за его гибкость, а в метапрограммировании и шаблонах я хорош, и мне нравится возиться с этим.

А еще хотелось проверить обещание, которое C++ повторяет уже не первое десятилетие, про zero-cost абстракции: удобный, читаемый код, который при этом ничего не стоит в рантайме, по-настоящему, с бенчмарками, тестами и дизассемблером. Правда ли он ничего не стоит? Вот это и хотелось пощупать. Поехали.

Начнем мы с самой базы, с вектора. Как оказалось, этого вполне достаточно, чтобы покрыть широкий спектр тем программирования низкоуровневых вычислений. При этом теории линейной алгебры мы касаться не будем, нас интересует не математика, а то, как она ложится на железо. Может ли выражение вроде dot(A,B)*t*C скомпилироваться в те же самые инструкции, что я написал бы руками интринсиками (встроенные функции, превращаются напрямую в низкоуровневые инструкции)? Не "сравнимо быстро", а буквально тот же дизассемблер. Если да, то у меня есть удобная математика без расплаты. Если нет, то я плачу за красоту, и надо знать, сколько. Разговор в основном пойдет про несколько вещей: железо, на котором мы что-то считаем и низкоуровневая работа с памятью (даже ниже, чем вы можете себе представить). Настоящая скорость живет в регистрах. Все остальное в этой статье - следствие этого факта.

Часть I. Обещание про zero-cost абстракции

Где живет быстрый код: регистры, не память

Когда говорят про оптимизацию в коде, чаще всего имеют в виду правильное представление данных. Главная точка преткновения - это память, она самая медленная, а хорошая раскладка данных (layout) - это первое, о чем стоит думать. Любая вычислительная библиотека будет подталкивать вас к этому: упакуйте данные плотно, однотипно, а обрабатывайте конвейерно, чтобы пролезть в кэш, и процессор как можно реже ходил в оперативную память. Этот пласт мы сейчас сознательно опустим и посмотрим уровнем ниже, на вычисление отдельных выражений.

И здесь, по сути, та же история. На фоне самих вычислений (1-4 такта) обращение к кэшу L1 - сопоставимая по стоимости операция (порядка 4-5 тактов). То есть "всего лишь" сходить в кэш для нас уже накладно. Поэтому главная задача - это удержать выражение в регистрах. А такие вещи, как спилл (spill, выгрузка значения из регистра в память; это дорого) - это уже смена режима производительности, десятки и сотни тактов, если речь доходит до RAM.

Таблица 1: Латентность операций для микроархитектуры Ice Lake - uops.info,uops.info/cache. Порядок величин; на сопоставимых микроархитектурах похож.

#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

Листинг 1: С++ код и его ассемблер: clean - чистая функция без спиллов в память, spilled - дополнительные вызовы и выгрузка временных переменных в память.

Здесь нужно честно оговориться: регистров конечное число, и кроме нашего выражения в них живут другие вычисления и служебный код - это называется давлением на регистры (register pressure). Иногда спилл - вынужденный и абсолютно корректный шаг. Поэтому дальше я буду аккуратно различать две вещи: вынужденный спилл из-за нехватки регистров и лишний спилл из-за того, что оптимизация не дожала. Бороться имеет смысл только со вторым.

Почему компилятор так хорош на скалярном коде и можно ли его догнать

На скалярных выражениях компилятор почти непобедим. За ним десятилетия оптимизаций, заточенных под конкретные архитектуры, он отлично распознает паттерны, переставляет код, планирует его под конкретный конвейер, распределяет регистры и подбирает инструкции. Руками в типичном скалярном выражении его не переиграть.

Что важно принять сразу. Цель нашего ручного кода не в том, чтобы обогнать компилятор на скаляре. Там мы в лучшем случае сыграем вничью, а скорее проиграем. Но есть тонкость, компилятор силен, пока у него есть гарантии. Там, где он теряет информацию о семантике - алиасинг (ситуация, когда два указателя могут указывать на одну область памяти или пересекающиеся области), escape-анализ (анализ компилятора, определяющий, может ли адрес переменной убежать за пределы функции; если не может, то переменную можно держать в регистре, а если может, то нужно выгружать в память), тонкости арифметики с плавающей точкой (например, обработка ситуаций с NaN), - он вынужден перестраховываться и быть консервативным.

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

Листинг 2.

Посмотрите на Листинг 2, dot3. Здесь компилятор не векторизовал вычисление - все операции остались скалярными. Зато он сделал другое: слил несколько умножений и сложений в одно умножение и два fma - это контракция FMA (Fused Multiply-Add - инструкция, которая делает умножение и сложение за одну операцию: a*b+c; вместо двух тактов - один). Теперь обратим внимание на axpy3 и axpy3_r. Казалось бы, одинаковые выражения, однако указание __restrict дает так много: вместо того чтобы грузить данные по очереди и считать их поодиночке (ss - scalar single-precision), компилятор грузит два компонента сразу и считает их упакованной инструкцией (ps - packed single-precision, не ss!), т.е. начинает сам векторизовать. Без гарантии он обязан считать, что запись в dst[0] могла задеть src[1] (вдруг указатели пересекаются), поэтому держит все последовательным скаляром. С __restrict, страх снят, и он по своей воле пакует пару массивов в SIMD. Обратите внимание, пакуется именно пара (x, y) одной ps-инструкцией, а третий компонент (z) досчитывается скаляром (ss) - нечетный хвост в SIMD не уложить. Вот здесь и зарыт наш хлеб. Не "как обогнать компилятор", а "где он сам останавливается". Останавливается он не потому, что глуп, а потому что ему не хватает доказательств. И часть нашей работы, эти доказательства ему предоставить. Но об этом чуть позже.

Можно ли догнать рукописные интринсики?

Да - и это та самая планка, к которой мы стремимся. Цель в том, чтобы абстракция выдавала байт-в-байт тот же код, что написал бы человек интринсиками. Не "близко", не "почти", а идентичный машинный код. Это, кстати, удобный критерий, сравниваем дизассемблер, а не ощущения.

Почему это вообще возможно? Интринсики - это тончайшая обертка над инструкциями процессора. Если абстракция компилируется в те же инструкции, то рантайм-разницы между ней и ручным кодом нет по определению. Что мешает? Каждый слой - это узлы выражения, обертки, ссылки, - это то, что оптимизатор обязан разобрать обратно, до уровня тех же интринсиков. Разобрал - получаем zero-cost: удобный код, который приятно читать и писать, и который не проигрывает по скорости. Не разобрал - получаем протечку, - абстракцию, за которую мы платим.

#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

Листинг 3: При другом порядке записи выражения компилятор расставит инструкции иначе - но это будет эквивалентный по работе ассемблер, с тем же набором инструкций. Здесь обе версии выстроены так, чтобы дизассемблер совпал буквально.

Итак, планка взята, я получил zero-cost абстракцию, байт-в-байт, как ручной код на интринсиках (Листинг 3). Притом без выгрузки временных объектов в память. Обсудим бенчмарки ниже.

Почему шаблоны - это абстракция без накладных расходов

Механизм такой. Обычная абстракция - это виртуальные вызовы, указатели на функции, стирание типов, - мы платим в рантайме за диспетчеризацию и невозможность инлайна. Шаблон же по своей сути ближе к подстановке: на каждый конкретный тип компилятор генерирует отдельный, специализированный код, который видит целиком. Вместо абстрактно написанного шаблона на выходе оказывается конкретный код под конкретную задачу. Ключевой момент в том, что структура вычисления (какие операции и в каком порядке) становится частью типа. Компилятор видит все дерево статически, а значит, может инлайнить и сворачивать его насквозь. Рантайм-решений нет, диспетчеризировать нечего.

И тут мы приходим к главному приему в моей математике - expression templates (ET). Это абстракция над абстракцией, шаблон, который строит не значение, а тип-выражение. dot(A,B)*t*C не вычисляется на месте, а собирается в тип вроде Mul< Mul< Dot<vec4,vec4>, float >, vec4 >. Само вычисление - это разворачивание этого типа компилятором.

// --- 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;
};

Листинг 4.

В итоге получается такая схема сборки выражения:

Схема 1: построение дерева выражения в tenvex

Цена, разумеется, не исчезает, но она переезжает в compile-time: время компиляции, раздувание типов, простыни в диагностике. Zero-cost в рантайме - это не zero-cost для разработчика, и за удобство пользователя библиотеки расплачивается ее автор.

Что такое zero-cost абстракция и с чем ее едят

Zero-cost - это обещание: не платишь за то, чего не используешь, а за то, что используешь, не платишь больше, чем заплатил бы за ручную реализацию. И обещание это весьма условное. Потому что zero-cost - это свойство не самого кода, а пары код + компилятор (точнее, конкретного бэкенда - части компилятора, генерирующей код для конкретной архитектуры, будь то x86_64 или ARM64). В нашей задаче оно складывается из нескольких независимых вещей: инлайна функций, SROA (Scalar Replacement of Aggregates - оптимизация, которая разбирает сложную структуру на отдельные переменные и помещает их в регистры вместо памяти), register promotion (продвижения значений из памяти в регистры) и DSE (Dead Store Elimination - оптимизация, которая удаляет записи в память, если их никто никогда не читает; убирает мусор). Достаточно чему-то одному не сложиться - и zero-cost пропадает, абстракция течет.

Таблица 2: медианные значения по 30 бенчмаркам, включая throughput - пропускные бенчмарки, где следующее вычисление не зависело от предыдущего, latency - бенчимарки, где каждое вычисление зависит от предыдущего. manual - ручная сборка из тех же ядер, что использует tenvex; overhead = tenvex/manual.

А раз обещание живет в оптимизаторе, значит, его можно проверить и можно нарушить. Отсюда и начинается вторая часть истории - борьба с компиляторами. По сути, она ровно про то, чтобы либо дать оптимизатору те самые недостающие гарантии, либо обнаружить, что конкретный бэкенд, один из проходов не вытягивает.

Глядя в сравнительную таблицу, мы видим, что обещание zero-cost трещит, но разберемся мы с этим во второй части статьи.

Zero-cost абстракция - это не то, что ты написал, а то, что оптимизатор согласился реализовать.

Часть II. Борьба с компиляторами

Краткое напоминание: zero-cost - это свойство пары код + компилятор, и держится оно на совокупности оптимизаций (инлайн, SROA, register promotion, DSE). Стоит одной не сложиться, и абстракция течет. Сравнительная таблица показала одну выпадающую ячейку: GCC на aarch64 - +39% против ручного кода на интрисинках. С ней и разберемся.

Поначалу я решил, что это дыра в библиотеке. Оказалось сложнее, и по дороге я несколько раз ошибся. Эти ошибки стоит показать, потому что они поучительнее самого фикса. Каждая про то, как удобная гипотеза проигрывает дизассемблеру.

Инструменты следствия

Прежде чем гоняться за багом, поговорим об инструментах. Их по сути четыре, и без первых трех любые рассуждения о производительности повисают в воздухе. Первый - это Google Benchmark, и здесь важно сразу развести два режима: latency-бенчмарк гоняет серийную цепочку, где каждое вычисление зависит от предыдущего, - он измеряет длину критического пути; throughput-бенчмарк гоняет независимые итерации, которые процессор перекрывает в конвейере, - он измеряет пропускную способность. Это разные числа про разное, и путать их не стоит; а что именно ты меряешь, решает то, как результат утекает наружу через sink вроде DoNotOptimize.

Второй инструмент - дизассемблер. Все, что я утверждаю дальше, проверяется не ощущением "стало быстрее", а сравнением дизассемблера: видно ли обращение к стеку, остался ли мертвый store (выделенная область памяти, куда что-то записывается, но эти значения никогда не читаются), упростился ли узел.

Третий - редукция кейса: способ пошагово свести реальную спиллящую единицу трансляции к минимуму, не теряя по дороге сам баг. Про него подробно в разворотном разделе, когда он понадобится.

Ну и четвертый инструмент - Claude AI. Я не слишком силен в asm - знаю самую базу, но разобраться в простынках сквозного дизассемблера и искать конкретные факты в них мне помогал именно Claude. Кроме того, он помогал мне готовить бенчмарки и конкретный код для проверок тех или иных гипотез. А также, он занялся написанием тестов для tenvex по моему описанию требований.

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);

Листинг 5: примеры бенчмарков tenvex на Google Benchmark.

Эти два бенчмарка меряют две физически разные характеристики одной и той же операции, и путать их является ошибкой, из-за которой потом не сходятся числа. Если коротко: latency отвечает на вопрос как долго считается одно значение, если следующее не может начаться без предыдущего, а throughput - сколько независимых таких операций я могу продавить в единицу времени, когда между ними нет зависимостей. Откуда берется разница на уровне железа? Возьми обычный mulps: его латентность ~3-4 такта, но throughput ~0.5 такта, то есть исполнительный блок конвейеризован и может принимать новый mul каждый такт (а то и два за такт), хотя каждый отдельный результат готов только спустя 4. Если ты измеряешь только latency, увидишь "4" и решишь, что операция дорогая. Если только throughput - увидишь "0.5" и решишь, что почти бесплатно. Оба числа верны, просто отвечают на разные вопросы.

Промах №1: "там спиллы в горячем пути"

Первая конкретная гипотеза была прямолинейной: раз +39%, значит, в цепочке вычислений спиллы - str/ldr гоняющие значения на стек и обратно. Логика простая, overhead - это лишние обращения к памяти. Я открыл asm, заранее зная, что увижу. И не увидел. Цепочка вычислений оказалась чистой, целиком в векторных регистрах, ни одной загрузки обратно со стека. Вместо ожидаемых спиллов там был мертвый кадр на стеке: записи в [sp] которые потом никто не читает, по несколько stp q, q, [sp] на итерацию. Это не спилл в смысле "не хватило регистров", значения преспокойно живут в регистрах. Это мертвые store, которые проход DSE не дожал и не удалил. Зверь другой, и лечится он иначе.

#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

Листинг 6

И вот тут моя гипотеза разошлась с реальностью. Спиллов в том смысле, в каком я их ждал: "посчитали, не хватило регистров, выгрузили, потом загрузили обратно", - здесь нет. Само вычисление чистое, значения живут в регистрах, их никто не вытесняет. Но при этом функция зачем-то аллоцирует на стеке 128 байт и пишет туда шесть раз, два из них - stp, по два регистра за инструкцию. Я искал каждую из этих записей на стороне чтения и не нашел: ни одного ldr обратно из этих адресов. Значения уходят на стек и там умирают, никто их не читает. Это не спилл из-за давления на регистры, это мертвые записи, которые почему-то не удалились. И DoNotOptimize, и явный .eval() в vf4 до записи кадр не убирали, потребитель результата был ни при чем. Почему компилятор их оставил, на этом этапе я еще не понимал.

"Медленно" и "лишние обращения к памяти" - это не синонимы, и предсказанный механизм может оказаться не тем. Я ждал спиллы из-за нехватки регистров, а нашел мертвые записи при чистом компьюте. Смотри в дизассемблинг, а не в модель у себя в голове.

Промах №2: "политика хранения - вот рычаг"

Это самая поучительная из ошибок. Узлы ET хранят операнды через tnvx_ref_or_value_t: маленькие листья хранятся по значению, крупные подвыражения - по const&. Гипотеза напрашивалась сама, именно const&-хранение заставляет узлы держать ссылки на временные узлы, которые приходится материализовать (выгружать в память); переключим все на by-value и кадр схлопнется. Я собрал игрушечный пример целиком на by-value, и он это подтвердил: чистый asm, никакого кадра. Окрыленный, я перевел на 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&>;

я поменял на

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&>;

Листинг 7.

Легкий фикс, хорошая точка для изменений, меняем в одном месте, а если в будущем придется что-то сделать с этим, мы знаем что и где делать.

И это был выстрел в ногу из… Бэкфайр. +39% превратились в +147-170%. Кадр не схлопнулся, он вырос: со ~128 байт до ~400, и вместо четырех мертвых store их стало пятнадцать. Чистая игрушка соврала: она была настолько стерильной, что не могла показать главного - реальный tenvex на by-value не лучше, а заметно хуже. Политика хранения вообще не была рычагом, я месил не то тесто.

Таблица 3: медианные значения по 30 бенчмаркам, включая throughput - пропускные бенчмарки, где следующее вычисление не зависело от предыдущего, latency - бенчимарки, где каждое вычисление зависит от предыдущего. manual — ручная сборка из тех же ядер, что использует tenvex; 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

Листинг 8.

Слишком чистый игрушечный репро может подтвердить ложную гипотезу, именно потому, что он чист. Он не содержит того отличия от реальности, в котором живет баг. Игрушка подтвердила "by-value спасает", а реальность ответила удвоением регрессии.

Разворот метода: репро, которое не воспроизводило

После провала by-value я пошел в другую сторону: строить минимальный репро снизу вверх. Из чистой игрушки, добавляя черты реального tenvex по одной. Добавил tnvx_ref_or_value_t. Добавил пустую CRTP-базу Expr<Derived>. Добавил вложенные агрегаты. Добавил раскладку с вектор-членом. Баг не появлялся ни на одном шаге - ни от одной черты по отдельности, ни от их комбинации. Минимальный репро оставался стерильно чистым, что бы я в него не досыпал.

(Я решил не вставлять тут листинг этих попыток найти баг - он занимает более 200 строк кода)

И вот это был поворотный момент. Не потому, что он нашел баг, а потому, что он сказал мне: метод неверен. Достраивая чистую игрушку, не получается поймать баг, который живет в разнице между игрушкой и реальностью, ведь именно эту разницу я каждый раз и выбрасывал. Правильный ход был противоположный: искать вниз из реальной спиллящей единицы трансляции, по одному заменяя настоящие компоненты на игрушечные, пока кадр не исчезнет. Тогда исчезнувший на каком-то шаге кадр и укажет на виновника.

Когда минимальный репро не воспроизводит баг, перестань достраивать чистое. Ищи вниз из реальной спиллящей TU, а не вверх из игрушки.

Развязка: это был баг бенчмарка

Когда я наконец стал искать вниз, причина оказалась до обидного приземленной. Она вообще была не в библиотеке. Она была в бенчмарке: составное выражение использовало одни и те же элементы массива - vb[i], vc[i] - дважды, прямо внутри ET-выражения. Повторяющееся обращение к массиву по индексу в роли листа.

На GCC/aarch64 это рушит SROA: оптимизатор не может доказать, что адрес материализованного узла не убегает, и на всякий случай оставляет узел в памяти, корректность важнее скорости. Отсюда мертвый ~400-байтный кадр и те самые stp ..., [sp] которые никто не читает. А почему игрушки никогда это не воспроизводили? Потому что они по совпадению всегда биндили листья в локальные переменные - я по-наитию написал это в игрушках, но не сделал этого в бенчмарках. Конструкторы вроде Leaf(va[i]) сначала забирали элемент в локальную переменную. То самое, что оказалось фиксом, игрушки делали по-случайности. Вот почему слишком чисто раз за разом подтверждало ложную гипотезу.

Фикс - одна строка в пользовательском коде, связать листья в локальные переменные до выражения (vec4 a = va[i], b = vb[i], c = vc[i];), а выражение писать уже над a, b, c. Мертвый кадр схлопывается в ноль. Я подтвердил это за пять шагов редукции: буквально второй шаг с биндингом листьев убирает кадр. Тем же артефактом, как выяснилось, были заражены еще и другие бенчмарки. После починки 6 из 6 ячеек zero-cost, в пределах шума измерений.

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;
    }
}

Листинг 9.

А вот тот же расчет, когда листья связаны в локальные переменные. Сверь с предыдущим листингом по арифметике, построчно: тот же fmla на a + b*2, те же fmul на произведения, та же горизонтальная редукция (ins в w-дорожку плюс пара faddp), fsqrt, fdiv, финальный fmla. Вычисления совпадают инструкция в инструкцию - это та же norm, тот же dot, тот же c*3. Математика идентична. Разница в том, чего здесь нет. Нет sub sp: функция вообще не открывает кадр. Из памяти трогаются ровно четыре места: три ldr на входе (va[i], vb[i], vc[i]) и один str на выходе (out[i]), то есть только чтение аргументов и запись результата. Промежуточные значения не покидают регистров: ни одного str в [sp], ни одного адреса, вычисленного под выгрузку. Вся обвязка мертвого кадра из loop_et - sub/add sp, три служебных адреса, шесть мертвых записей здесь просто отсутствует. Это и есть доказательство: лишний кадр - не плата за саму математику, она-то совпадает байт в байт, а артефакт формы, в которой математика записана. Стоит связать элементы массива в локальные переменные, и адрес каждого становится доказуемо неубегающим, SROA разбирает выражение на регистры, кадр схлопывается в ноль. То самое исправление в одну строку.

#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)));
    }
}

_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

Листинг 10.

Таблица 4: медианные значения по 30 бенчмаркам, включая throughput - пропускные бенчмарки, где следующее вычисление не зависело от предыдущего, latency - бенчимарки, где каждое вычисление зависит от предыдущего. manual — ручная сборка из тех же ядер, что использует tenvex;

Микробенчмарк обязан зеркалить реальный паттерн доступа. Слишком чистый бенчмарк не нейтрален, он подтверждает ложь: мой стенд мерил артефакт, которого в нормальном коде нет.

Шесть ячеек из шести - zero-cost. Та единственная, что выбивалась, держалась не на свойствах библиотеки, а на одной строке в стенде: повторный vb[i] внутри выражения. Убрали повтор и tenvex и ручной код сошлись (в пределах шума измерений). На GCC/arm64 это буквально один и тот же машинный код. Стоит проговорить, что именно это означает, и чего не означает. Обещание zero-cost устояло, абстракция ET не стоит ничего в рантайме относительно того же выражения, собранного руками из тех же ядер. Но "не стоит ничего" здесь - это про накладные расходы абстракции, а не про то, что tenvex - это самый быстрый способ посчитать это выражение в принципе. Мы все это время отвечали на один вопрос: течет ли абстракция? Ответ: нет, не течет (а где почти потекла, виноват был стенд). Вопрос "а где она вообще дает выигрыш и насколько?" другой, и к нему мы скоро перейдем.

Откуда берется мертвый кадр

Теперь, когда причина ясна, теория складывается чисто. У локальной переменной адрес доказуемо не убегает; у элемента массива - да. Когда один и тот же элемент появляется дважды как лист, escape-анализ GCC/aarch64 не может доказать, что материализованный узел не убежит, и оставляет его в памяти. А проход DSE на этом бэкенде не дожимает образовавшиеся мертвые записи, отсюда кадр. Свяжешь лист в локальную переменную и адрес становится доказуемо неубегающим, SROA разбирает узел, и кадра нет.

Осталось показать обратную сторону - где абстракция не просто не протекает, а реально выигрывает. К этому и перейдем.

Часть III. Где мы выигрываем и где нет

Мы выяснили, что абстракция не течет. Остался обратный вопрос, что она вообще выигрывает? Соперников два: можно выигрывать у наивной реализации - у человека, написавшего то же самое в лоб, можно мериться силами с компилятором - с тем, насколько далеко он сам дотянет наивный код оптимизациями. Сразу предупрежу, парада побед не будет. Кое-где tenvex проигрывает и по делу. Эти места и окажутся самыми интересными.

Сводка и как ее читать

Сначала сырые числа: наивная скалярная реализация, tenvex и их отношение по операциям.

Таблица 5.

Одна оговорка про чтение. Числа заметно разнятся от компилятора к компилятору, это нормально и по-своему интересно, но это отдельная история, и здесь я в нее лезть не буду. Latency и throughput - это числа про разное, длина критического пути против пропускной способности.

Три класса операций

Операции распадаются на три класса, и в каждом своя история.

Первый - element-wise, которые компилятор и сам векторизует. Наивный код и tenvex компилируются в одно и то же. Выигрыша ноль, работу уже сделали за нас, конкурировать не с чем.

#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

Листинг 11: дизассемблинг для tenvex и наивной реализаций эквивалентен по вычислительному ядру; различия сводятся к регистровой раскладке, выбранной компилятором, и использованию vmovaps вместо vmovups (tnvx::vec4 - выровнен по 16 байт).

Вторая: element-wise, где автовекторизатор сдается, это победа. У saturate (зажать в [0,1] - min(max(v,0),1)) обе версии получаются векторными, но наивная тратит лишнюю инструкцию: компилятор сворачивает std::min/std::max в cmp+blendvps+minps, тогда как tenvex обходится cmp+andps+minps. То же и у floor/ceil/frac - tenvex берет одну упакованную операцию и выигрывает стабильно (frac ~1.5×, reflect ~1.3×). Преимущество ровно там, где компилятор не смог свернуть так же плотно.

#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

Листинг 12.

Дело тут в std::min/std::max и обработке NaN. std::min(a, b) определен как b < a ? b : a, при равенстве и при NaN у него своя семантика возврата операнда. А minps/maxps на NaN ведут себя иначе: при NaN в любом операнде возвращают второй операнд, и порядок аргументов у них вшит в саму инструкцию. То есть minps это не побитовый эквивалент четырех std::min: на NaN-входах они дадут разный результат. Компилятор обязан сохранять наблюдаемое поведение программы. Раз minps может разойтись с std::min на NaN, он не имеет права заменить одно другим и, не имея гарантии, что NaN не придет, перестраховывается и оставляет скалярный код. Это не глупость оптимизатора, а его честность, он не знает, что в данных NaN не будет.

Третья, чуть более тонкая. Несколько latency-связанных ячеек на arm, где tenvex проигрывает по-настоящему. AABB по latency: наивные min/max компилятор сворачивает в fcmgt+bit (сравнение плюс битселект), а tenvex берет fmin/fmax, а они на Cortex-A76 в зависимой цепочке дороже из-за NaN-семантики. Где результат каждой итерации ждет предыдущей, более дешевый по латентности compare-select наива выигрывает. На x86_64 этого нет, там min/max-инструкции дешевы, и те же операции у tenvex впереди. Это цена микроархитектуры.

#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

Листинг 13.

Казалось бы, tenvex должен выигрывать на arm64 - короче по числу инструкций, две против четырех. Но проигрывает по latency, потому что дело не в количестве, а в латентности на критическом пути.

Горизонтальные операции над широкими регистрами дороги по latency

Самый интересный случай, dot3/dot4. По latency dot lпроигрывает везде: на arm 0.44× у dot3, 0.53× у dot4. Это честная цена горизонтальной свертки (операции, которая схлопывает 4 числа в широком регистре в одно) на критическом пути, и причина видна в asm: длинная shuffle-цепочка против короткого FMA-скаляра. Горизонтальная свертка дорога по латентности. В throughput стоимость свертки частично прячется за перекрытием независимых итераций: критический путь больше не виден, важна занятость портов. Но прячется не значит исчезает: где у tenvex есть запас по портам, разрыв переворачивается (dot4 на arm выходит вперед, до 1.81×), a где шаффлы все еще в горлышке, остается около ничьей (dot3 - слабый плюс, балансирующий на нуле). Опираться стоит на устойчивое: проигрыш по latency.

#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

Листинг 14.

И опять инструкций меньше, а латентность выше. В latency важна длина критического пути. Наивный dot3: компилятор подтягивает компоненты парой shufps, затем mul и два слитых fma - короткий хвост на FMA-портах. SIMD-dot3 тащит горизонтальную свертку: mul, затем shuffle чтобы свести дорожки, затем зависимые сложения. Цепочка длиннее, и два шаффла садятся на латентный shuffle-порт. В throughput же, когда дотов много и они независимы, критический путь прячется за перекрытием итераций: одна векторная загрузка плюс редукция у tenvex перекрываются лучше разрозненных скалярных чтений наива, но шаффлы все еще заняты, и dot3 в throughput лишь подтягивается к ничьей, а не выигрывает. Уверенный переворот случается там, где портов хватает, у dot4.

Как с этим живут взрослые библиотеки? Так же. DirectXMath делает ту же горизонтальную свертку (mul → shuffle → add → shuffle → add → splat в baseline-SSE), ветвясь по доступным инструкциям: скаляр без интринсик, vadd+vpadd на NEON, dpps на SSE3, два hadd на SSE3, shuffle-редукция на базовом SSE. GLM - то же самое для vec4: dpps на SSE4.1, два hadd на SSE3, shuffle+add на базовом. У обеих по несколько реализаций под разные наборы инструкций, и сам факт этих #if-деревьев доказывает, что единого быстрого способа для горизонтальной операции нет.

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
}

Листинг 15: реализации dot в DirectXMath (слева) и GLM (справа)

Честный вывод про латентный dot: оптимизировать его особо некуда, горизонтальные операции вшиты в задачу "один vec4 → один скаляр". Выход: если dot много и важна именно латентность, то нужно менять не операцию, а раскладку данных. SoA (Structure of Arrays): вместо массива векторов держать структуру из массивов, и тогда горизонталь становится вертикалью, dot считается как x*x + y*y + z*z по дорожкам, без единого шаффла. В tenvex SoA пока нет; - это направление, к которому ведет логика и куда давно ушла индустрия (потоковые функции DXM, SoA в геймдеве), а не измеренный мной результат.

Зачем тогда ET?

Тут уместен вопрос, который напрашивается сам. DirectXMath обходится обычными инлайн-функциями, никаких expression templates. После инлайна компилятор и так видит выражение целиком. Зачем тогда вся ET-машинерия?

Для скорости самого сгенерированного оптимизатором кода незачем, и это стоит сказать прямо. После инлайна оптимизатор видит один и тот же поток инструкций, что у цепочки функций, что у дерева узлов, и сворачивает промежуточные значения одинаково. Помните loop_manual? - плоские вызовы ядер без всяких ET? Он компилировался не хуже ET-версии. Одно и то же выражение ET не превратят в более быстрый машинный код, чем хорошо заинлайненные функции.

Но есть второй рычаг, а вот он уже про оптимизации по-настоящему. ET могут менять не качество кода, а сам объем вычислений, считать не то же самое быстрее, а попросту меньше. Функция обязана вернуть свой полный результат, она не знает, как он будет использован. Узел ET знает контекст, он часть дерева, и верхний eval() может переписать поддерево. magnitude3(a) < magnitude3(b) - оба sqrt лишние, достаточно сравнить квадраты: и это компилятору из обычных функций недоступно, потому что требует алгебраического знания, что sqrt монотонна, а не просто выкидывания мертвого кода. Бывают и более простые случаи, ленивость: взяли только .x() от cross3 - остальные компоненты можно не считать (хотя такое хороший компилятор иногда вычищает и сам через DCE). Плюс эргономика: a + b*2 вместо XMVectorAdd(a, XMVectorScale(b, 2)). И эргономика - это не мало. Все мы на своем опыте знаем, что код, который хорошо читается, проще отлаживать (именно для этого мы заметаем под ковёр сложность отладки expression templates).

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()));
}

}

Листинг 16: Эскиз: как мог бы выглядеть один expression-level rewrite. В tenvex его нет, сейчас magnitude3_sq выбирается руками. Показано, чтобы потенциал ET стал осязаемым, а не абстрактным.

Итог

Я начинал с полушутки про самую быструю библиотеку линейной алгебры и честного вопроса, правда ли zero-cost абстракция ничего не стоит. Ответ есть, и он условный. Zero-cost - свойство не кода и не компилятора по отдельности, а их пары: обещание держится, когда складываются инлайн, SROA, продвижение в регистры и удаление мертвых записей. Если сложились, то ET-выражение компилируется в тот же код, что и ручные интринсики. Не сложилось хоть одно звено, и абстракция течет. Но это не значит, что абстракцию не стоит делать. Течет она редко и предсказуемо, а чинится чаще всего просто, одной строкой (хотя найти проблему может быть не совсем тривиально). Писать же интринсики руками на всем коде ради того процента, где она дала течь, размен явно невыгодный. На всех шести бекендах, когда стенд наконец стал честным, обещание устояло.

Главный урок - про честный стенд. Большую часть пути моим противником был не компилятор, а мои собственные измерения. Те самые +39% оказались не протечкой абстракции, а мертвым кадром от повторного обращения по индексу в моем же бенчмарке, и сидел этот баг в половине стенда, выдавая мнимые регрессии там, где их не было. Прежде чем винить бэкенд, проверь, что измеряешь, и смотри в дизассемблинг, а не в свою модель происходящего, которая у тебя в голове. Написать же хороший бенчмарк - это целое искусство.

Важно и то, чего zero-cost не касается: форы над наивным кодом он не обещал. Там, где автовекторизатор сдался, выигрывает широкая операция, где компилятор силен сам - ничья. А кое-где честно проигрываешь: горизонтальные операции дороги по латентности, fmin дорог в зависимой цепочке. И отдельно про expression templates, ради которых все затевалось - для скорости машинного кода они не нужны, после инлайна компилятор сворачивает то же выражение так же. Их реальная ниша - сделать меньше работы, выкинуть лишний sqrt, не считать ненужную компоненту, но это пока задел, а не реализованное. Так что сегодня ET дают эргономику, а скорость живет в регистрах, в раскладке данных и в том, насколько компилятору хватило гарантий. Стоила ли вся эта возня одного вектора? По исходному счету да. Я хотел пощупать руками SIMD, поиграть с оптимизациями и проверить zero-cost обещание на прочность, и понимание, из чего складывается скорость, вышло честнее и интереснее.