std::simd: шаблоны векторизации без intrinsics в C++ / Хабр

Привет, Хабр!

Векторизация в C++ давно живёт на двух этажах. Внизу автоворожбы компилятора: достаточно аккуратно написать цикл, и при нужных флагах он соберёт SIMD‑инструкции сам. Наверху низкоуровневые intrinsics, где вы контролируете каждый shuffle и predication, но платите за это портируемостью и временем на поддержку. Между ними появился удобный этаж: std::simd. Он даёт вам явные векторные типы и операции без прыжков по AVX/NEON‑интринсикам и при этом остаётся переносимым. В 2025 году картина такая: полноценный std::simd принят в стандарт C++26, а использовать в проде уже сейчас можно std::experimental::simd из Parallelism TS v2, которая есть в libstdc++ (GCC 11+) и постепенно доезжает до остальных реализаций. Для чтателя это означает простую вещь: сегодня пишем под <experimental/simd>, а миграция на <simd> будет механической.

Коротко о модели std::simd

Библиотека даёт типы вида simd<T, Abi> и маски simd_mask<T, Abi>. ABI‑теги задают ширину и представление вектора: чаще всего вам нужен simd_abi::native<T> для нативной ширины под текущие флаги компиляции, или simd_abi::fixed_size<N> если вы хотите фиксированную ширину, одинаковую на x86 и ARM. Есть теги выравнивания для загрузок/выгрузок: element_aligned, vector_aligned, overaligned<Align>. Все арифметические и логические операции действуют поэлементно.

Важно про статус API. В C++26 это уже <simd> и перегрузки математических функций для SIMD есть в стандарте. В TS v2 (то, что в libstdc++) интерфейс находится в <experimental/simd>, и именование чуть отличается, но семантика та же. Если вы хотите «нативную» ширину SIMD сейчас, берите std::experimental::native_simd<T> (алиас поверх simd<T, simd_abi::native<T>>).

Минимальные флаги для релиза и автопреобразований:

GCC/Clang: -O3 -march=native -fno-math-errno -fno-trapping-math
-ffast-math даёт скорость, но будьте осторожны в численных ядрах вроде softmax.
Проверять автогенерацию векторализаторов удобно через отчёты:
GCC: -fopt-info-vec семейство, LLVM.

Базовые паттерны: скелет цикла, загрузки, маски, хвосты

Рассмотрим безопасный шаблон для прохода по массиву float:

#include <experimental/simd>
#include <cstddef>
#include <cstring>

namespace stdx = std::experimental;

template<class T>
using vnat = stdx::native_simd<T>;                     // нативная ширина
template<class T>
using mnat = stdx::simd_mask<T, typename vnat<T>::abi_type>;

void add_inplace(float* __restrict dst,
                 const float* __restrict src,
                 std::size_t n)
{
    using V = vnat<float>;
    using M = mnat<float>;
    constexpr std::size_t W = V::size();

    std::size_t i = 0;

    // основной векторный проход
    for (; i + W <= n; i += W) {
        // безопасно: contiguous, элементное выравнивание
        V a; a.copy_from(dst + i, stdx::element_aligned);
        V b; b.copy_from(src + i, stdx::element_aligned);
        a += b;
        a.copy_to(dst + i, stdx::element_aligned);
    }

    // хвост через маску
    if (i < n) {
        const std::size_t tail = n - i;
        // маска вида [1..tail, 0..]
        M mask = stdx::where(stdx::simd<std::size_t, typename V::abi_type>(
                                 [](auto j){ return j < tail ? 1u : 0u; }
                             ) != 0, M(true));
        V a, b;
        a.copy_from(dst + i, mask, stdx::element_aligned);
        b.copy_from(src + i, mask, stdx::element_aligned);
        a += b;
        a.copy_to(dst + i, mask, stdx::element_aligned);
    }
}

copy_from/copy_to с тегами выравнивания. Если памяти гарантированно хватает и она выровнена под «вектор», используйте vector_aligned или overaligned<N>.

«Хвост» лучше делать через маску, а не отдельным скалярным циклом..

Про выравнивание контейнеров. Дефолтный std::vector не гарантирует векторное выравнивание. Если хотите NV‑friendly 32/64 байта, используйте кастомный аллокатор или специализированные контейнеры. Вариант с аллокатором:

template<class T, std::size_t Align>
struct aligned_allocator {
  using value_type = T;
  T* allocate(std::size_t n) {
    void* p = nullptr;
    if (posix_memalign(&p, Align, n * sizeof(T)) != 0) throw std::bad_alloc();
    return static_cast<T*>(p);
  }
  void deallocate(T* p, std::size_t) noexcept { free(p); }
};

Или хотя бы сообщите компилятору про выравнивание указателя там, где это корректно: std::assume_aligned<32>(ptr). Но это обещание, а не проверка: если оно ложное, будет UB.

Микро-ядро 1: редукция суммы

Базовый скаляр:

float sum_scalar(const float* a, std::size_t n) noexcept {
    float s = 0.f;
    for (std::size_t i = 0; i < n; ++i) s += a[i];
    return s;
}

Векторный вариант:

float sum_simd(const float* a, std::size_t n) noexcept {
    using V = vnat<float>;
    constexpr std::size_t W = V::size();

    V acc = 0.f;
    std::size_t i = 0;
    for (; i + W <= n; i += W) {
        V v; v.copy_from(a + i, stdx::element_aligned);
        acc += v;
    }

    // горизонтальная редукция + хвост
    float s = stdx::reduce(acc); // sum по всем lane
    for (; i < n; ++i) s += a[i];
    return s;
}

std::experimental::reduce определён как горизонтальная операция, есть также hmin/hmax. Важно помнить, что ассоциативность у вас должна быть допустима для произвольной группировки. Для суммы float это обычно нормально, но для строгой воспроизводимости включайте Kahan или задавайте порядок.

А/Б‑каркас для измерений:

#include <chrono>
#include <random>
#include <iostream>

template<class F>
static double bench_ms(F&& f, int iters = 10) {
    using clk = std::chrono::steady_clock;
    double best = 1e300;
    for (int i = 0; i < iters; ++i) {
        auto t0 = clk::now();
        auto volatile sink = f();
        auto t1 = clk::now();
        double ms = std::chrono::duration<double, std::milli>(t1 - t0).count();
        if (ms < best) best = ms;
    }
    return best;
}

int main() {
    const std::size_t N = 1 << 24;
    std::vector<float, aligned_allocator<float, 64>> a(N);
    std::mt19937 rng(123);
    std::uniform_real_distribution<float> dist(0.f, 1.f);
    for (auto& x : a) x = dist(rng);

    double t_scalar = bench_ms([&]{ return sum_scalar(a.data(), N); });
    double t_simd   = bench_ms([&]{ return sum_simd  (a.data(), N); });

    std::cout << "scalar " << t_scalar << " ms\n";
    std::cout << "simd   " << t_simd   << " ms\n";
}

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

Микро-ядро 2: softmax с численной устойчивостью

Стабильный softmax по вектору x состоит из трёх этапов: найти максимум, вычесть его, посчитать exp, просуммировать, поделить. Чисто скалярно это выглядит очевидно. Сделаем SIMD‑версию с аккуратной обработкой хвоста и двумя реализациями экспоненты: а) через SIMD‑перегрузки std::exp (когда у вас уже <simd> из C++26), б) через приближение полиномом для TS v2.

struct SoftmaxResult {
    float sum;
    float maxv;
};

template<class V>
static inline V exp_poly(V z) noexcept {
    // Простое полиномиальное приближение exp на ограниченном диапазоне.
    // Для продакшена ограничьте вход, чтобы |z| ≤ ~10.
    const V c1 = 1.0f;
    const V c2 = 1.0f;
    const V c3 = 0.5f;
    const V c4 = 1.0f/6.0f;
    const V c5 = 1.0f/24.0f;
    const V c6 = 1.0f/120.0f;
    return c1 + z*(c2 + z*(c3 + z*(c4 + z*(c5 + z*c6))));
}

SoftmaxResult softmax_pass1_findmax(const float* x, std::size_t n) {
    using V = vnat<float>;
    constexpr std::size_t W = V::size();
    V vmax = std::numeric_limits<float>::lowest();
    std::size_t i = 0;
    for (; i + W <= n; i += W) {
        V v; v.copy_from(x + i, stdx::element_aligned);
        vmax = stdx::max(vmax, v);
    }
    float m = stdx::hmax(vmax);
    for (; i < n; ++i) m = std::max(m, x[i]);
    return { /*sum=*/0.f, /*maxv=*/m };
}

void softmax_inplace(float* x, std::size_t n) {
    using V = vnat<float>;
    using M = mnat<float>;
    constexpr std::size_t W = V::size();

    auto [_, maxv] = softmax_pass1_findmax(x, n);
    const V vmax(maxv);
    V vsum = 0.f;

    std::size_t i = 0;
    for (; i + W <= n; i += W) {
        V v; v.copy_from(x + i, stdx::element_aligned);
        v -= vmax;

        // Вариант A: когда есть SIMD-перегрузки std::exp (C++26 <simd>)
        // v = std::exp(v);

        // Вариант B: TS v2 — приближение полиномом
        v = exp_poly(v);

        v.copy_to(x + i, stdx::element_aligned);
        vsum += v;
    }

    float sum = stdx::reduce(vsum);
    for (; i < n; ++i) {
        float z = std::exp(x[i] - maxv); // хвост можно посчитать точно
        x[i] = z;
        sum += z;
    }

    const V vsumv(sum);
    i = 0;
    for (; i + W <= n; i += W) {
        V v; v.copy_from(x + i, stdx::element_aligned);
        v /= vsumv;
        v.copy_to(x + i, stdx::element_aligned);
    }
    for (; i < n; ++i) x[i] /= sum;
}

В C++26 перегрузки математических функций для SIMD типов уже есть, std::exp(v) работает поэлементно. В TS v2 этого может не быть в вашей либе; тогда либо используйте аккуратное приближение, либо делегируйте в специализированные векторные math‑библиотеки.

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

Микро-ядро 3: горизонтальный фильтр изображений 3×3

Разобьём классическое 3×3 на две одно‑мерные свёртки: горизонтальную и вертикальную. В горизонтали данные лежат подряд, SIMD заходит «в лоб». Границы и правый хвост закроем масками.

struct Image {
    int w, h, stride; // stride в пикселях
    float* data;      // одиночный канал для простоты
};

void box3x3_horizontal(const Image& in, Image& tmp) {
    using V = vnat<float>;
    using M = mnat<float>;
    constexpr std::size_t W = V::size();

    const float k[3] = { 1.f/3, 1.f/3, 1.f/3 };
    const V k0(k[0]), k1(k[1]), k2(k[2]);

    for (int y = 0; y < in.h; ++y) {
        const float* src = in.data + y * in.stride;
        float* dst = tmp.data + y * tmp.stride;

        // левая кромка скаляром
        if (in.w >= 1) dst[0] = (src[0] + src[1]) * 0.5f;
        if (in.w >= 2) dst[1] = (src[0] + src[1] + src[2]) / 3.f;

        int x = 2;
        for (; x + (int)W + 1 <= in.w; x += (int)W) {
            V a, b, c;
            a.copy_from(src + x - 1, stdx::element_aligned);
            b.copy_from(src + x,     stdx::element_aligned);
            c.copy_from(src + x + 1, stdx::element_aligned);

            V r = a * k0 + b * k1 + c * k2;
            r.copy_to(dst + x, stdx::element_aligned);
        }

        // хвост: x в [in.w - W - 1, in.w)
        if (x < in.w) {
            const int rem = in.w - x;
            M mask([&](auto i){ return i < (std::size_t)rem; });

            V a, b, c;
            a.copy_from(src + x - 1, mask, stdx::element_aligned);
            b.copy_from(src + x,     mask, stdx::element_aligned);
            // Для c у последнего элемента читаем src[in.w-1]
            // Сформируем охранную копию:
            alignas(64) float guard[64];
            std::memcpy(guard, src + x + 1, rem * sizeof(float));
            for (int t = rem; t < (int)W; ++t) guard[t] = src[in.w - 1];
            c.copy_from(guard, stdx::element_aligned);

            V r = a * k0 + b * k1 + c * k2;
            r.copy_to(dst + x, mask, stdx::element_aligned);
        }

        // правая кромка скаляром
        if (in.w >= 2) dst[in.w - 2] = (src[in.w - 3] + src[in.w - 2] + src[in.w - 1]) / 3.f;
        if (in.w >= 1) dst[in.w - 1] = (src[in.w - 2] + src[in.w - 1]) * 0.5f;
    }
}

Вертикальный проход можно сделать блоками по нескольку строк, чтобы держать рабочее окно в L1/L2. Хорошая эвристика: высота блока такая, чтобы три строки источника и одна строка назначения плюс пару строк «запаса» умещались в L1, с учётом ширины в байтах. Классическая литература по кеш‑паттернам напоминает: стриминговые проходы, предсказуемые адреса, prefetch при больших шагах.

Где хватит авто-векторизации, а где нет

Когда можно положиться на -O3 и std::execution::unseq:

Простые линейные проходы со сведением зависимостей к редукции.
Отсутствие ветвлений внутри горячего цикла.
Нет «непонятных» компилятору вызовов, или он умеет их мэппить на инструкции (пример из LLVM docs: floorf превращается в roundps при наличии SSE4.1).

Когда std::simd даёт выигрыши и стабильность:

Сложные хвосты, которые автопроход ломает на контроле потока.
Маскированные операции и раннее отбрасывание элементов.
Тяжёлая математика внутри итерации, например softmax или активации, где компилятор боится лезть из‑за точности и исключений.
Специфичные паттерны загрузок/выгрузок и строгие требования к выравниванию, которые нужно выписать явно.

x86 против ARM: ширина, предикация, SVE

На x86 вы обычно живёте в фиксированных ширинах: SSE 128, AVX2 256, AVX-512 512. На ARM исторически NEON 128, а в серверном мире появился SVE/SVE2 со скалируемой длиной вектора, которую выбирает конкретное железо. Это меняет стиль: под SVE хочется писать «длино‑агностичный» код. std::simd частично закрывает эту историю: native<T> подбирает «правильную» ширину под флаги компилятора, а libstdc++ уже втащил поддержку SVE в std::experimental::simd. Если вам нужна фиксированная, берите fixed_size<N>, но учитывайте компромисс с SVE.

Cледствия:

Не пришивать ширину константой в коде, если не вынуждены.
Тестировать на ARM не только NEON, но и SVE. В ряде окружений SVE и NEON сосуществуют, и существуют удобные мосты между ними на уровне компиляторов и заголовков.

Маски и where: идиомы

Маска ― это тип, соответствующий ABI вашего вектора. С ним удобно писать условные обновления без ветвлений.

template<class T>
void clamp_above(T* a, std::size_t n, T hi) {
    using V = vnat<T>;
    using M = stdx::simd_mask<T, typename V::abi_type>;
    constexpr std::size_t W = V::size();

    std::size_t i = 0;
    for (; i + W <= n; i += W) {
        V v; v.copy_from(a + i, stdx::element_aligned);
        M m = v > V(hi);
        stdx::where(m, v) = V(hi); // только там, где m=true
        v.copy_to(a + i, stdx::element_aligned);
    }
    for (; i < n; ++i) if (a[i] > hi) a[i] = hi;
}

where работает и как «маска‑загрузка/выгрузка», и как «маска‑обновление». Это компилируется в предицированные инструкции на архитектурах, где они есть, либо в blend/shuffle там, где предикации нет

Ещё два микро-паттерна для реального кода

Пороговый фильтр с подсчётом (mask + reduce)

std::size_t count_greater_than(const float* a, std::size_t n, float thr) {
    using V = vnat<float>;
    using M = mnat<float>;
    constexpr std::size_t W = V::size();
    std::size_t i = 0;
    std::size_t cnt = 0;

    for (; i + W <= n; i += W) {
        V v; v.copy_from(a + i, stdx::element_aligned);
        M m = v > V(thr);
        cnt += stdx::popcount(m); // количество true в маске
    }
    for (; i < n; ++i) if (a[i] > thr) ++cnt;
    return cnt;
}

popcount(mask) — очень удобная штука для всяких фильтров и QPS‑счётчиков.

Свёртка 1D фиксированного ядра (расписываем FMA):

template<int K>
void conv1d_valid(const float* x, const float* k, float* y, std::size_t n) {
    using V = vnat<float>;
    constexpr std::size_t W = V::size();
    std::size_t end = n - K + 1;

    // загрузим коэффициенты как скаляры и будем вещать в V
    V kv[K];
    for (int t = 0; t < K; ++t) kv[t] = V(k[t]);

    std::size_t i = 0;
    for (; i + W <= end; i += W) {
        V acc = 0.f;
        for (int t = 0; t < K; ++t) {
            V v; v.copy_from(x + i + t, stdx::element_aligned);
            acc = stdx::fma(v, kv[t], acc); // если реализация подтянет FMA
        }
        acc.copy_to(y + i, stdx::element_aligned);
    }
    for (; i < end; ++i) {
        float s = 0.f;
        for (int t = 0; t < K; ++t) s += x[i + t] * k[t];
        y[i] = s;
    }
}

Некоторые ошибочки, которые можно допустить

Игнорировать выравнивание. Если у вас большие массивы под горячие ядра, позаботьтесь о 32/64-байтовом выравнивании и сообщайте об этом библиотеке copy_from/copy_to тегами.

«Хвост» через отдельный второй цикл, который ломает предсказание ветвлений. Маска понятнее и часто быстрее.

Полагаться на «магический» -ffast-math в численно хрупких местах. Выигрыш может оказаться ложным. Смотрите рекомендации LLVM/ARM по авто‑векторизации и режимам математики.

Жёстко кодировать ширину под AVX2, а потом пытаться гонять то же на NEON/SVE. Используйте native<T> или слоями выделяйте fixed_size<N> там, где это действительно нужно.

Итог

Если коротко: std::simd ― это тот случай, когда «средний этаж» наконец удобен. Он закрывает 80 процентов задач без боли intrinsics и без молитв компилятору, что он сам догадается.

Если вам интересно разобраться, как современные стандарты C++ меняют подход к работе с памятью, многопоточностью и производительностью, загляните на курс C++ Developer. Professional. Программа построена вокруг практики: стандартные библиотеки, асинхронность, профилирование и оптимизация. Пройдите бесплатное тестирование по курсу, чтобы оценить свои знания и навыки.

Отзыв студента курса C++ Developer. Professional — *Отзыв студента курса* ***C++ Developer. Professional***

Рост в IT быстрее с Подпиской — дает доступ к 3-м курсам в месяц по цене одного. Подробнее

std::simd: шаблоны векторизации без intrinsics в C++