SIMD and auto-vectorisation

A modern core can add four floats in the time it takes to add one — if the data and the code are shaped for it. SIMD (single instruction, multiple data) instructions operate on short fixed-width vectors of lanes. This skill covers when the compiler does that for you, why it often refuses, and how to write it yourself, portably, with exact handling of the awkward parts: the tail, alignment, masks and reductions.

1. Lanes and registers

A plain g++ -O2 build for x86-64 may only assume SSE2. AVX needs -mavx2 (or -march=native, which also makes the binary unportable to older CPUs). That's why this app's grader and your Compiler tab default to 128-bit vectors. The AArch64 track will come back to NEON and SVE.

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

int main() {
    // native_simd<T>: as many lanes as the target's native register holds.
    std::println("float lanes: {}", stdx::native_simd<float>::size());
    std::println("int32 lanes: {}", stdx::native_simd<std::int32_t>::size());
    std::println("byte lanes:  {}", stdx::native_simd<std::uint8_t>::size());
}

On baseline x86-64 that prints 4, 4 and 16.

2. Auto-vectorisation: what the compiler needs

GCC and Clang vectorise loops on their own (GCC from -O2 with a cautious cost model, more aggressively at -O3). Ask GCC what it did with -fopt-info-vec-optimized and what it refused with -fopt-info-vec-missed:

$ g++ -O3 -c kernels.cpp -fopt-info-vec-optimized
kernels.cpp:2: optimized: loop vectorized using 16 byte vectors and unroll factor 4
kernels.cpp:2: optimized:  loop versioned for vectorization because of possible aliasing
kernels.cpp:7: optimized: loop vectorized using 16 byte vectors and unroll factor 16
$ g++ -O3 -c kernels.cpp -fopt-info-vec-missed
kernels.cpp:6: missed: not vectorized, possible dependence between data-refs

A loop vectorises when:

• The trip count is computable before the loop starts. for (i = 0; i < n; ++i) is fine even if n is only known at run time. A loop with an early return or break (a search) usually isn't vectorised by GCC on x86-64.
• Iterations are independent. a[i] += a[i-1] (a prefix sum) needs the previous iteration's result: a loop-carried dependency.
• Pointers don't alias, or the compiler can check. For a[i] += b[i] with int* a, const int* b, a and b might overlap. At -O3 GCC versions the loop (a runtime overlap check, then a vector or a scalar path); at -O2 it gives up. __restrict (a common extension) or separate arrays remove the question.
• Access is contiguous. Unit stride is cheap; strided or indirect access (a[idx[i]]) needs gathers, which are slow or absent on baseline x86-64.
• Reductions can be reordered. See section 3.

Use the Compiler tab: paste a loop, pick -O3, and look for packed instructions — paddd (4 × int32 add), addps (4 × float add), pcmpeqb (16 byte comparisons), vpaddd/vaddps with AVX. Scalar code shows addss, add, cmp per element instead.

3. Reductions and floating point

Vectorising a sum means keeping V::size() partial sums and combining them at the end. That reorders the additions. For integers that's fine — wrapping addition is associative. For floating point it isn't:

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

int main() {
    float xs[4] = {1e8f, 1.0f, -1e8f, 1.0f};

    float s = 0;
    for (float x : xs) s += x;        // ((1e8 + 1) - 1e8) + 1

    using V = stdx::fixed_size_simd<float, 2>;
    V acc = 0.0f;                      // lanes: (1e8 + -1e8), (1 + 1)
    for (int i = 0; i < 4; i += 2) acc += V(xs + i, stdx::element_aligned);

    std::println("sequential {}, lane-wise {}", s, stdx::reduce(acc));
}

That prints sequential 1, lane-wise 2: 1e8 + 1 isn't representable as a float and rounds back to 1e8. So GCC at -O3 loads a float sum in vectors but still adds the elements one at a time (addss) to keep the source's order. You grant permission to reorder with -ffast-math (which also assumes no NaNs or infinities — a bigger promise than it sounds) or the narrower -fassociative-math, or by writing the lane-wise sum yourself. Either way, results can change in the last bits: compare float SIMD results with a tolerance, never ==.

4. Explicit SIMD with portable types

When the compiler won't (or you want guarantees), write the vector code yourself. Two options:

• Intrinsics — one function per instruction, tied to an instruction set. Precise, unportable, verbose:

  #include <emmintrin.h>   // SSE2, always available on x86-64
  #include <print>
  
  int main() {
      alignas(16) int a[4] = {1, 2, 3, 4}, b[4] = {10, 20, 30, 40}, c[4];
      __m128i va = _mm_load_si128(reinterpret_cast<const __m128i*>(a));
      __m128i vb = _mm_load_si128(reinterpret_cast<const __m128i*>(b));
      _mm_store_si128(reinterpret_cast<__m128i*>(c), _mm_add_epi32(va, vb));  // paddd
      std::println("{} {} {} {}", c[0], c[1], c[2], c[3]);
  }
  

• Portable SIMD types — a value type with one element per lane and ordinary operators, compiled to whatever the target has. <experimental/simd> (the Parallelism TS v2) is the TS ancestor of C++26 std::simd: same model, slightly different names, and it works in both GCC and Clang here. It's what the build task uses.

Worked example: a float sum with explicit lanes

#include <experimental/simd>
#include <cstddef>
#include <print>
#include <span>
#include <vector>
namespace stdx = std::experimental;

float sum(std::span<const float> x) {
    using V = stdx::native_simd<float>;
    V acc = 0.0f;                                    // broadcast: every lane 0
    std::size_t i = 0;
    for (; i + V::size() <= x.size(); i += V::size())
        acc += V(x.data() + i, stdx::element_aligned);   // load one chunk
    float s = stdx::reduce(acc);                     // combine the lanes
    for (; i < x.size(); ++i) s += x[i];             // the tail
    return s;
}

int main() {
    std::vector<float> v(10);
    for (std::size_t i = 0; i < v.size(); ++i) v[i] = float(i) * 0.5f;
    std::println("{}", sum(v));   // 0 + 0.5 + ... + 4.5 = 22.5
}

Every explicit SIMD loop has this shape:

1. Broadcast scalars you need in every lane (V acc = 0.0f;).
2. Chunk loop while a whole chunk fits: i + V::size() <= n. Using i < n would load past the end on the last iteration — undefined behaviour that ASan reports and that crashes for real when the buffer ends at a page boundary.
3. Load with element_aligned (any address) unless you know the address is aligned to stdx::memory_alignment_v<V>, in which case vector_aligned may be faster. A subspan at an arbitrary offset is never known to be aligned.
4. Store results back explicitly (v.copy_to(p + i, element_aligned)): a simd value is a copy, not a view.
5. Reduce lane-wise accumulators once, after the loop.
6. Tail: finish the last 0 … V::size()-1 elements with a scalar loop.

What it buys, measured on this machine with g++ -O2 (1 Mi floats, best of 20 runs — your numbers will differ, run it yourself):

float sum, 1M floats: scalar 846 us, simd 212 us (4.0x)

The scalar loop can't be vectorised without permission (section 3); the explicit version adds four lanes per instruction.

5. Masks and branch-free code

Comparisons between simd values produce a mask — one bool per lane. Masks replace branches:

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

int main() {
    stdx::fixed_size_simd<int, 4> v([](int i) { return i * 3; });   // 0 3 6 9
    auto big = v > 4;                                                // F F T T

    std::println("count {}", stdx::popcount(big));                  // 2
    std::println("first {}", stdx::find_first_set(big));            // 2
    std::println("any {} all {}", stdx::any_of(big), stdx::all_of(big));

    stdx::where(big, v) = 0;               // masked assignment: only the T lanes
    std::println("sum {}", stdx::reduce(v)); // 0 + 3 + 0 + 0 = 3

    auto c = stdx::min(stdx::max(v, decltype(v)(1)), decltype(v)(2)); // clamp to [1, 2]
    std::println("clamped {} {} {} {}", int(c[0]), int(c[1]), int(c[2]), int(c[3]));
}

• popcount(mask) counts matches without a branch per element.
• where(mask, v) = x (and += etc.) is a blend: lanes outside the mask keep their value.
• min/max compile to minps/maxps — a branch-free clamp. Unpredictable data makes branches expensive (see the profiling skill); masks cost the same whatever the data.

Masked tails. A tempting alternative to the scalar tail is "load the last partial chunk padded with zeros". The padding then takes part in every comparison: counting x > -1 counts each zero. Either pad with a value that can't affect the result, or mask out the padding lanes.

Faded example: counting values above a threshold

Complete it (the answer is below):

std::size_t count_above(std::span<const int> xs, int t) {
    using V = stdx::native_simd<int>;
    const V vt = /* ??? */;
    std::size_t n = 0, i = 0;
    for (; /* ??? */; i += V::size()) {
        V v(xs.data() + i, stdx::element_aligned);
        n += /* ??? */;
    }
    for (; i < xs.size(); ++i) n += /* ??? */;
    return n;
}

6. Counting in lanes: watch the width

Byte kernels (counting newlines, character classes) get 16 lanes per SSE register. If you count matches in byte lanes — where(v == '\n', counts) += 1 — each counter is an unsigned char and wraps after 255. Either add popcount(mask) into a std::size_t every chunk, or flush the byte counters into a wider total at least every 255 chunks. The same applies to any narrow accumulator lane (int16, float precision).

7. Data layout for SIMD

SIMD wants each lane's data next to each other. An array of structs {x, y, z, mass} interleaves fields, so loading "4 x values" needs a gather or shuffles. A struct of arrays (xs[], ys[], …) makes each field one contiguous stream — load, compute, store. That's the AoS → SoA transformation from the memory-layout skill; data-oriented design takes it further.

8. C++26 std::simd

C++26 standardises the same model in <simd>, under namespace std::simd. GCC 16 already ships an early implementation (it doesn't define __cpp_lib_simd yet, so treat it as experimental; Clang 23's library doesn't have it):

#include <simd>
#include <cstddef>
#include <print>
#include <span>
#include <vector>
namespace simd = std::simd;

std::size_t count_above(std::span<const int> xs, int t) {
    using V = simd::vec<int>;                        // TS: native_simd<int>
    std::size_t n = 0, i = 0;
    for (; i + V::size() <= xs.size(); i += V::size())
        n += simd::reduce_count(simd::unchecked_load<V>(xs.subspan(i, V::size())) > t);
    for (; i < xs.size(); ++i) n += xs[i] > t;
    return n;
}

int main() {
    std::vector<int> v{1, 9, 3, 12, 5, 7, 20, 2, 11};
    std::println("{}", count_above(v, 6));   // 5
}

Summary

• SIMD lanes: 4 floats/ints or 16 bytes per SSE/NEON register; more with AVX.
• Compilers vectorise independent, countable, contiguous loops; aliasing, loop-carried dependencies, early exits and float reductions stop them. Check with -fopt-info-vec-missed or the Compiler tab.
• Float sums reorder additions when vectorised: only with permission, and the result changes in the last bits.
• Explicit SIMD shape: broadcast → chunk loop while a whole chunk fits → element_aligned loads → store back → reduce → scalar tail.
• Masks (popcount, where, min/max) replace branches; padding and narrow lane counters are the classic correctness traps.
• <experimental/simd> today; C++26 std::simd has the same model.