SIMD Intrinsics: Vectorizing with AVX2 and NEON

SIMD — Single Instruction, Multiple Data — lets one instruction operate on 4, 8, or 16 values simultaneously. Every modern CPU has wide vector registers: AVX2 processes 8 floats or 32 chars in a single instruction; ARM NEON processes 4 floats or 16 bytes. When you write float c = a + b, you're leaving 7/8 of the hardware idle.

This page covers x86 AVX2 intrinsics in depth (the 95% case on modern desktop/server hardware), NEON patterns for ARM, the GCC vector extension as a cleaner alternative, and std::simd (C++26) for portable code.

When to use SIMD

SIMD is the right tool when:

You have a tight loop over a large contiguous array
The bottleneck is compute (not memory bandwidth or latency)
The auto-vectorizer fails (inspect assembly: SIMD instructions start with v on x86)
You need a specific 3–12× speedup that other optimizations cannot achieve

Profile first. SIMD code is hard to read, write, and maintain. Auto-vectorization (enabled by -O2) already handles many simple cases. Only write intrinsics when you've confirmed the compiler missed the opportunity and the hotspot matters.

x86 register width and type naming

cpp

Godbolt

SSE (1999):  16 × 128-bit registers: xmm0–xmm15
AVX (2011):  16 × 256-bit registers: ymm0–ymm15
AVX-512:     32 × 512-bit registers: zmm0–zmm31

C/C++ vector types that map to these registers:

Type	Width	Holds
`__m128`	128-bit	4 × float
`__m128d`	128-bit	2 × double
`__m128i`	128-bit	16 × int8, 8 × int16, 4 × int32, 2 × int64
`__m256`	256-bit	8 × float
`__m256d`	256-bit	4 × double
`__m256i`	256-bit	32 × int8, 16 × int16, 8 × int32, 4 × int64

Intrinsic naming convention: _mm<width>_<operation>_<type>

cpp

Godbolt

_mm256_add_ps   → add 8 floats (ps = packed single)
_mm256_add_pd   → add 4 doubles (pd = packed double)
_mm256_add_epi32 → add 8 int32s (epi32 = extended packed int32)
_mm256_mul_ps   → multiply 8 floats
_mm256_fmadd_ps → fused multiply-add: a*b + c (FMA3, 2013+)
_mm256_loadu_ps → load 8 floats (unaligned)
_mm256_storeu_ps → store 8 floats (unaligned)

Setup

cpp

Godbolt

#include <immintrin.h>  // includes everything: SSE + AVX + AVX2 + FMA

// Tell the compiler we're targeting AVX2 (alternative to -mavx2 flag):
#pragma GCC target("avx2,fma")

// Runtime CPU feature detection (no segfault on CPUs without AVX2):
bool has_avx2 = __builtin_cpu_supports("avx2");
bool has_avx512 = __builtin_cpu_supports("avx512f");
bool has_neon = false;  // compile-time, not runtime, for ARM

cmake

# Tell the compiler to enable AVX2 for the whole target:
target_compile_options(mylib PRIVATE -mavx2 -mfma)

# Or use modern CMake ISA flags:
include(CheckCXXCompilerFlag)
check_cxx_compiler_flag("-mavx2" COMPILER_SUPPORTS_AVX2)
if(COMPILER_SUPPORTS_AVX2)
    target_compile_options(mylib PRIVATE -mavx2 -mfma)
endif()

Load, compute, store — the basic pattern

cpp

Godbolt

#include <immintrin.h>
#include <vector>

// Scalar baseline: a[i] + b[i]
void add_scalar(const float* a, const float* b, float* c, int n) {
    for (int i = 0; i < n; ++i)
        c[i] = a[i] + b[i];
}

// AVX2: adds 8 floats per iteration instead of 1
void add_avx2(const float* a, const float* b, float* c, int n) {
    int i = 0;
    // Main loop: process 8 elements at a time
    for (; i + 8 <= n; i += 8) {
        __m256 va = _mm256_loadu_ps(a + i);  // load 8 floats from a
        __m256 vb = _mm256_loadu_ps(b + i);  // load 8 floats from b
        __m256 vc = _mm256_add_ps(va, vb);   // add all 8 pairs
        _mm256_storeu_ps(c + i, vc);         // store 8 results
    }
    // Scalar tail: handle remainder when n is not a multiple of 8
    for (; i < n; ++i)
        c[i] = a[i] + b[i];
}

Aligned vs unaligned loads

cpp

Godbolt

// _mm256_loadu_ps — unaligned, works anywhere (may be 1 cycle slower on some CPUs)
__m256 v = _mm256_loadu_ps(ptr);

// _mm256_load_ps — requires ptr to be 32-byte aligned, slightly faster
alignas(32) float buf[8];
__m256 v = _mm256_load_ps(buf);  // guaranteed aligned

// Use std::aligned_alloc or alignas for aligned heap/stack buffers
float* data = static_cast<float*>(std::aligned_alloc(32, n * sizeof(float)));

Dot product — data dependency and multiple accumulators

A single accumulator creates a data dependency: each FMA must wait 4–5 cycles for the previous result before it can start. The fix is multiple independent accumulators that drain into a horizontal sum at the end.

cpp

Godbolt

#include <immintrin.h>
#include <numeric>

// Single accumulator (slow — bottlenecked by 5-cycle FMA latency)
float dot_single(const float* a, const float* b, int n) {
    __m256 acc = _mm256_setzero_ps();
    int i = 0;
    for (; i + 8 <= n; i += 8) {
        __m256 va = _mm256_loadu_ps(a + i);
        __m256 vb = _mm256_loadu_ps(b + i);
        acc = _mm256_fmadd_ps(va, vb, acc);  // acc += va * vb
    }
    // Horizontal sum of 8 floats in acc
    __m128 lo = _mm256_castps256_ps128(acc);
    __m128 hi = _mm256_extractf128_ps(acc, 1);
    __m128 sum128 = _mm_add_ps(lo, hi);
    __m128 shuf = _mm_movehdup_ps(sum128);
    __m128 sums = _mm_add_ps(sum128, shuf);
    shuf = _mm_movehl_ps(shuf, sums);
    sums = _mm_add_ss(sums, shuf);
    float result = _mm_cvtss_f32(sums);
    // Scalar tail
    for (; i < n; ++i) result += a[i] * b[i];
    return result;
}

// Four accumulators (2-4x faster — hides FMA latency)
float dot_four_accumulators(const float* a, const float* b, int n) {
    __m256 acc0 = _mm256_setzero_ps();
    __m256 acc1 = _mm256_setzero_ps();
    __m256 acc2 = _mm256_setzero_ps();
    __m256 acc3 = _mm256_setzero_ps();

    int i = 0;
    for (; i + 32 <= n; i += 32) {
        acc0 = _mm256_fmadd_ps(_mm256_loadu_ps(a+i),    _mm256_loadu_ps(b+i),    acc0);
        acc1 = _mm256_fmadd_ps(_mm256_loadu_ps(a+i+8),  _mm256_loadu_ps(b+i+8),  acc1);
        acc2 = _mm256_fmadd_ps(_mm256_loadu_ps(a+i+16), _mm256_loadu_ps(b+i+16), acc2);
        acc3 = _mm256_fmadd_ps(_mm256_loadu_ps(a+i+24), _mm256_loadu_ps(b+i+24), acc3);
    }

    // Combine accumulators
    __m256 acc = _mm256_add_ps(_mm256_add_ps(acc0, acc1),
                               _mm256_add_ps(acc2, acc3));

    // Horizontal sum + scalar tail (same as above)
    // ...
    float result = /* hadd */ 0.0f;
    for (; i < n; ++i) result += a[i] * b[i];
    return result;
}

Rule: any tight loop with a single reduction accumulator should use 4+ independent accumulators to fill the CPU's pipeline depth.

Practical example: RGB to grayscale

Converting 24-bit RGB to 8-bit grayscale using 16-bit fixed-point arithmetic — processing 16 pixels at once with SSE2, vs 8 with AVX2 floats:

cpp

Godbolt

#include <immintrin.h>
#include <cstdint>

// Luminance weights: Y = 0.2126*R + 0.7152*G + 0.0722*B
// Fixed-point Q15: multiply by 32768 (1 << 15)
static constexpr int16_t W_R = static_cast<int16_t>(0.2126f * 32768);
static constexpr int16_t W_G = static_cast<int16_t>(0.7152f * 32768);
static constexpr int16_t W_B = static_cast<int16_t>(0.0722f * 32768);

void rgb_to_gray_sse2(const uint8_t* rgb, uint8_t* gray, int pixel_count) {
    const __m128i wr = _mm_set1_epi16(W_R);
    const __m128i wg = _mm_set1_epi16(W_G);
    const __m128i wb = _mm_set1_epi16(W_B);
    const __m128i zero = _mm_setzero_si128();

    int i = 0;
    for (; i + 16 <= pixel_count; i += 16) {
        // Load 16 pixels (48 bytes of interleaved RGB)
        __m128i p0 = _mm_loadu_si128(reinterpret_cast<const __m128i*>(rgb + i*3));
        __m128i p1 = _mm_loadu_si128(reinterpret_cast<const __m128i*>(rgb + i*3 + 16));
        __m128i p2 = _mm_loadu_si128(reinterpret_cast<const __m128i*>(rgb + i*3 + 32));

        // Deinterleave: separate R, G, B channels (simplified — real code uses shuffles)
        // ... (omitted for brevity; use _mm_shuffle_epi8 with a deinterleave mask)

        // Fixed-point multiply-high: (r * W_R) >> 15
        // _mm_mulhi_epu16 returns the high 16 bits of the 32-bit product
        __m128i yr = _mm_mulhi_epu16(r16, wr);
        __m128i yg = _mm_mulhi_epu16(g16, wg);
        __m128i yb = _mm_mulhi_epu16(b16, wb);

        // Sum channels and pack 16-bit back to 8-bit
        __m128i y16 = _mm_add_epi16(_mm_add_epi16(yr, yg), yb);
        __m128i y8  = _mm_packus_epi16(y16, zero);  // saturate to [0, 255]
        _mm_storel_epi64(reinterpret_cast<__m128i*>(gray + i), y8);
    }
    // Scalar tail
    for (; i < pixel_count; ++i) {
        int r = rgb[i*3], g = rgb[i*3+1], b = rgb[i*3+2];
        gray[i] = static_cast<uint8_t>(0.2126f*r + 0.7152f*g + 0.0722f*b);
    }
}

Fixed-point 16-bit fits 2× more values per register than 32-bit float, and _mm_mulhi_epu16 has the same throughput as _mm_mul_ps. This is why real-world image processing typically runs on 16-bit integer SIMD rather than floats.

GCC vector extensions — cleaner syntax

GCC (and Clang) provide a way to define SIMD types with operator overloading, avoiding the intrinsic naming chaos:

cpp

Godbolt

// Define vector types with __attribute__((vector_size))
typedef float  v8f  __attribute__((vector_size(32)));  // 8 × float (256-bit)
typedef int    v8i  __attribute__((vector_size(32)));  // 8 × int   (256-bit)
typedef double v4d  __attribute__((vector_size(32)));  // 4 × double (256-bit)

v8f a = {1, 2, 3, 4, 5, 6, 7, 8};
v8f b = {8, 7, 6, 5, 4, 3, 2, 1};

v8f c = a + b;   // calls _mm256_add_ps under the hood
v8f d = a * b;   // _mm256_mul_ps
v8f e = a / b;   // _mm256_div_ps
v8i mask = a > b; // compare → mask (all-ones where true, zero where false)

// Mix with intrinsics when needed (zero-cost conversion):
int movemask = _mm256_movemask_ps((__m256)a);

Vector extensions are much more readable for simple arithmetic. Use intrinsics when you need specific operations (shuffles, horizontal sums, blends) that have no C++ operator equivalent.

ARM NEON — same ideas, different syntax

On ARM (Apple M-series, Cortex-A, Raspberry Pi), use <arm_neon.h>:

cpp

Godbolt

#include <arm_neon.h>

// NEON: 128-bit registers, process 4 floats at once
void add_neon(const float* a, const float* b, float* c, int n) {
    int i = 0;
    for (; i + 4 <= n; i += 4) {
        float32x4_t va = vld1q_f32(a + i);  // load 4 floats
        float32x4_t vb = vld1q_f32(b + i);
        float32x4_t vc = vaddq_f32(va, vb); // add 4 pairs
        vst1q_f32(c + i, vc);               // store 4 results
    }
    for (; i < n; ++i) c[i] = a[i] + b[i];
}

// FMA: a * b + c (very common in ML inference)
float32x4_t result = vmlaq_f32(c, a, b);  // c + a*b

x86 AVX2	ARM NEON	Width	Operation
`_mm256_loadu_ps`	`vld1q_f32`	128-bit	load floats
`_mm256_add_ps`	`vaddq_f32`	128-bit	add floats
`_mm256_mul_ps`	`vmulq_f32`	128-bit	multiply floats
`_mm256_fmadd_ps`	`vmlaq_f32`	128-bit	FMA
`_mm256_storeu_ps`	`vst1q_f32`	128-bit	store floats

std::simd — portable SIMD (C++26)

std::simd (standardized in C++26, in <experimental/simd> in GCC 11+ as std::experimental::simd) lets you write portable SIMD code that compiles to the widest register the target supports:

cpp

Godbolt

#include <experimental/simd>
namespace stdx = std::experimental;

using floatv = stdx::native_simd<float>;  // width = hardware max (4/8/16 floats)

void add_portable(const float* a, const float* b, float* c, int n) {
    constexpr int W = floatv::size();  // 8 on AVX2, 4 on SSE/NEON

    int i = 0;
    for (; i + W <= n; i += W) {
        floatv va(a + i, stdx::element_aligned);
        floatv vb(b + i, stdx::element_aligned);
        floatv vc = va + vb;
        vc.copy_to(c + i, stdx::element_aligned);
    }
    for (; i < n; ++i) c[i] = a[i] + b[i];
}

// Reduction
float dot_portable(const float* a, const float* b, int n) {
    constexpr int W = floatv::size();
    floatv acc = 0.0f;
    int i = 0;
    for (; i + W <= n; i += W) {
        floatv va(a + i, stdx::element_aligned);
        floatv vb(b + i, stdx::element_aligned);
        acc += va * vb;
    }
    float result = stdx::reduce(acc);  // horizontal sum
    for (; i < n; ++i) result += a[i] * b[i];
    return result;
}

cmake

target_compile_options(myapp PRIVATE -std=c++23 -mavx2)
# Or for C++26 once widely supported:
target_compile_options(myapp PRIVATE -std=c++26 -mavx2)

Status: std::simd was merged into the C++26 standard (P1928). As of 2026, it is available via -std=c++26 in GCC 14+ and Clang 18+.

Performance pitfalls

Gather/scatter — the latency killer

SIMD loads assume contiguous data. Gathering non-contiguous elements is 5–30× slower than a sequential load:

cpp

Godbolt

// Gather (slow — vgatherdps has 20+ cycle latency)
__m256i indices = _mm256_set_epi32(7, 6, 5, 4, 3, 2, 1, 0);
__m256 gathered = _mm256_i32gather_ps(base_ptr, indices, 4);  // stride in bytes

// Solution: reorganize your data to be SoA (struct of arrays), not AoS
// AoS (bad for SIMD):   [x,y,z,x,y,z,x,y,z,...]
// SoA (good for SIMD):  [x,x,x,...][y,y,y,...][z,z,z,...]

Alignment matters

cpp

Godbolt

// Misaligned 32-byte reads can split a cache line, adding latency:
float* data = static_cast<float*>(std::aligned_alloc(32, n * sizeof(float)));

// Or use alignas for stack/global:
alignas(32) float buf[256];

// _mm256_load_ps requires 32-byte alignment; _mm256_loadu_ps doesn't:
__m256 v = _mm256_load_ps(buf);     // safe: buf is 32-byte aligned
__m256 v = _mm256_loadu_ps(ptr);    // safe: ptr can be any alignment

Horizontal reductions are expensive

Avoiding them is often the key to performance:

cpp

Godbolt

// Bad: horizontal sum in every loop iteration
for (int i = 0; i < n; ++i) {
    __m128 v = _mm_loadu_ps(a + i);
    result += _mm_cvtss_f32(_mm_hadd_ps(_mm_hadd_ps(v, v), v));
}

// Good: accumulate in a vector register, horizontal sum once at the end
__m256 acc = _mm256_setzero_ps();
for (int i = 0; i + 8 <= n; i += 8)
    acc = _mm256_add_ps(acc, _mm256_loadu_ps(a + i));
// One horizontal sum at the end

Library alternatives

For most use cases, a SIMD library gives 90% of the performance with 10% of the code complexity:

Library	Approach	Platform	Notes
Highway (Google)	Portable vectors, runtime dispatch	x86/ARM/WASM	Recommended for new code
xsimd	Portable SIMD API	x86/ARM	Header-only, good Eigen integration
VCL (Agner Fog)	x86-specific, very complete	x86 only	Extremely well-documented
std::simd	Standard (C++26)	Compiler-dependent	The future; limited today
Eigen	Auto-vectorized matrices	All	Handles SIMD transparently

Google Highway in particular provides automatic architecture dispatch: one binary, multiple code paths selected at runtime based on the CPU's SIMD capabilities.

Inspecting what the compiler generates

bash

# Build with AVX2 and look at the assembly
g++ -O3 -mavx2 -mfma -S -o out.s myfile.cpp

# Check if the compiler auto-vectorized a loop (look for ymm registers)
grep -i "ymm\|vzero\|vbroadcast" out.s | head -20

# Use Godbolt (compiler explorer) online for quick checks
# Compiler options: -O3 -mavx2 -mfma

# Check what intrinsics map to with -fverbose-asm
g++ -O3 -mavx2 -fverbose-asm -S myfile.cpp

If the compiler already generates ymm instructions on a plain scalar loop, you don't need intrinsics — the auto-vectorizer did the work.