Guidebook · SIMD in Rust

SIMD in Rust: From the Machine to the Stage

A guidebook for a working engineer who wants to understand SIMD well enough to teach it, build it, reason about it, and defend it in front of 500 people.

By the end you will be able to do four things.

1. Explain what SIMD is, what a lane, a mask and a shuffle are, and why eight lanes almost never means eight times faster.
2. Build real kernels in stable Rust with std::arch: byte search, sums, filters, a matrix multiply, a CSV scanner. Each one is checked against a plain scalar version.
3. Decide when SIMD is the right tool, when the compiler already did the job, and when memory is the real limit.
4. Defend your choices with measurements, because every kernel in this book comes with numbers from a real machine.

How to use this guidebook

Every chapter answers three questions, always in the same order:

Do the exercises. Reading is half the learning. Each chapter ends with a Check yourself block. Questions are marked Warm-up, Core or Stretch, so you can pick your level. Hints are hidden under the word Hint, and answers to the hard ones are at the end of each Part.

One more habit: predict before you run. When a section says "predict first", write the number on paper. The gap between your guess and the real number is where the learning happens.

Everything here was run

All code in this book lives in one Rust crate, and the book you are reading was assembled by a script that copies the code and the program output straight from that crate. Nothing was typed in by hand. The tests pass. The numbers are real. Where something could not be run on the machine used for this book (ARM NEON, WebAssembly), the text says Not run: check this on your machine.

The machine used for this book

This is a shared virtual machine, so timings are noisy. Treat every number as "about this size", and expect your numbers to differ. What should not differ is the shape: which kernel is faster than which, and where the cliffs are.

What you need before you start

Skills. You can write Rust, you know what a slice and a Vec are, and you have seen unsafe. You do not need to know any assembly or any SIMD. We build it from nothing.

Tools.

# Rust 1.89 or newer (AVX-512 intrinsics became stable in 1.89; safe target_feature in 1.86)
rustup update stable
rustc --version        # this book: 1.91.1

# Linux x86-64 is assumed. Useful tools:
objdump --version      # to read the machine code the compiler produced

Which parts need which hardware?

Using an ARM laptop? Put #[cfg(target_arch = "x86_64")] around the x86 modules (the crate already does this) and read the x86 code as a worked example of the ideas. Every idea has an ARM twin, and Chapter 14 shows where to look.

Check your CPU:

grep -m1 flags /proc/cpuinfo | tr ' ' '\n' | grep -E '^(avx2|fma|avx512f|avx512bw|avx512vnni|pclmulqdq)$'

Study paths

Table of contents

Part 1. The picture

1. What SIMD is
2. The machine under your loop
3. Lanes, masks and shuffles
4. Numbers inside lanes
5. Memory is the real boss

Part 2. Seeing SIMD shapes in your code

6. Which loops are SIMD-shaped
7. Data layout decides everything
8. Tails, alignment and aliasing

Part 3. Writing it in Rust

9. Let the compiler try first
10. std::arch and the three contracts
11. Byte kernels: find, count, sum, clamp
12. Dispatch: picking the kernel at runtime
13. AVX-512 in Rust
14. Other worlds: NEON, WebAssembly, std::simd

Part 4. Algorithms that real systems use

15. Dot products and matrix multiply
16. Int8 and the saturation trap
17. Filtering and compaction
18. The nibble trick: classifying bytes with a shuffle
19. Bitmasks and carry-less multiply: a CSV scanner

Part 5. Measuring and verifying

20. Benchmark science
21. Proving a kernel correct

Part 6. Why it is this way

22. Why SIMD exists, and the roads not taken
23. Designing and improving SIMD systems

Part 7. SIMD in the wild

24. Reading real code

Back matter

• Twelve capstone projects
• Teaching kit
• What changed from the earlier SIMD guide
• Glossary
• Reading map
• What I could not verify

Part 1. The picture

Before any instruction names or Rust syntax, we need a clear picture in your head. By the end of this Part you will be able to say what a vector register holds, why a chain of dependent operations is slow, how numbers behave inside lanes, and why memory so often decides the final speed. Every later Part leans on these five chapters.

Chapter 1. What SIMD is

In this chapter: you meet the idea in one picture, run your first real SIMD instruction, and learn the one sentence that explains almost everything: one instruction, many values.

You need first: nothing.

The problem first

You write this loop every week:

for i in 0..n {
    c[i] = a[i] + b[i];
}

You think of it as "n additions". And on a simple CPU it is: load one number, load another, add, store. Then again. And again.

But the chip in your machine is not simple. Inside it are wide registers, 256 or even 512 bits across, and the hardware to do many additions at the same moment. If your loop adds one pair at a time, most of that width sits idle. It is like driving a truck to deliver one parcel.

SIMD is the name for using that width.

The model (the "what")

SIMD stands for Single Instruction, Multiple Data. One instruction. Many pieces of data. The data sit side by side in one wide register, and the instruction acts on all of them at once.

Look at the figure first. Four boxes on top are one register. Each box is a lane. A lane is one slot that holds one value.

 a (one 128-bit register, 4 lanes of f32)       b
 +------+------+------+------+                 +------+------+------+------+
 |  1.0 |  2.0 |  3.0 |  4.0 |                 | 10.0 | 20.0 | 30.0 | 40.0 |
 +------+------+------+------+                 +------+------+------+------+
    |       |       |       |                      |       |       |       |
    +---+---+       |       |                      |       |       |       |
        v           v       v                      v       v       v       v
 +------+------+------+------+       one instruction: vaddps
 | 11.0 | 22.0 | 33.0 | 44.0 |       four additions at the same time
 +------+------+------+------+
 c = a + b, lane by lane

Lane 0 only meets lane 0. Lane 1 only meets lane 1. No lane looks at its neighbour. This is the most important rule, and we will lean on it for the whole book: a normal vector operation is lane-wise.

The register's total size is its width. The lane count is the width divided by the size of one value.

An analogy. Imagine a stamp with eight heads, each one holding the same pattern. You press once and eight pages get stamped. That is SIMD: eight pieces of paper, one press.

Where the analogy breaks. Real stamping lets you choose a different stamp for each page. SIMD does not. All eight heads must do the same operation. And pages cannot stamp themselves based on how the page next to them turned out. The lanes are independent. Hold on to those two limits: they explain most of what is hard about SIMD.

A name worth knowing. In 1966, Michael Flynn sorted computers into four boxes by how many instruction streams and data streams they have. SIMD is one of the four. The other three are SISD (one instruction, one data: the simple CPU), MISD (rare) and MIMD (many cores, each doing its own thing). You will hear "SIMD vs MIMD" in talks, and now you know where it comes from.

Build it (the "how")

Here is the smallest real SIMD program. It adds eight floats with one instruction.

// src/bin/hello.rs
use std::arch::x86_64::*;

/// Add two arrays of 8 floats with ONE vector add instruction.
#[target_feature(enable = "avx")]
fn add8(a: [f32; 8], b: [f32; 8]) -> [f32; 8] {
    let mut out = [0.0f32; 8];
    // SAFETY: each array is exactly 8 floats = 32 bytes, and loadu/storeu need no alignment.
    unsafe {
        let va = _mm256_loadu_ps(a.as_ptr()); // 8 floats -> one 256-bit register
        let vb = _mm256_loadu_ps(b.as_ptr());
        let vc = _mm256_add_ps(va, vb); // eight additions, one instruction
        _mm256_storeu_ps(out.as_mut_ptr(), vc);
    }
    out
}

fn main() {
    let a = [1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0];
    let b = [10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0];
    assert!(is_x86_feature_detected!("avx"));
    // SAFETY: we just checked that this CPU supports AVX.
    let c = unsafe { add8(a, b) };
    println!("{:?}", c);
}

Run it:

[11.0, 22.0, 33.0, 44.0, 55.0, 66.0, 77.0, 88.0]

Three things to notice.

First, #[target_feature(enable = "avx")]. This tells the compiler "inside this function, you may use AVX instructions". Without it, the compiler would refuse to emit them, because not every CPU has AVX. We will spend all of Chapter 10 on this.

Second, we asked the CPU at runtime, with is_x86_feature_detected!, whether AVX exists, before calling the function. That check is what makes the unsafe call safe. If you ran this on a CPU without AVX, the program would crash with an "illegal instruction" error.

Third, look at what the compiler produced. I asked objdump for the machine code of add8:

    b1d0:	vmovups ymm0,YMMWORD PTR [rsi]
    b1d4:	vaddps ymm0,ymm0,YMMWORD PTR [rdx]
    b1d8:	vmovups YMMWORD PTR [rdi],ymm0

One load, one vaddps, one store. That vaddps is the "eight additions at once". ps means packed single-precision: many 32-bit floats packed in one register. (The load of b was folded into the add as a memory operand, which is why you see only one vmovups load.)

Why it is this way (the "why")

Why does the hardware have wide registers at all? Why not just make the clock faster?

For about three decades chips got faster mostly by running at higher clock speeds. The reason this worked is a rule of thumb from Robert Dennard and his colleagues (1974): as transistors shrink, you can raise the frequency without using more power per area. Around the mid-2000s that stopped being true. Power became the limit, not the number of transistors. Chip designers had many more transistors than they could run faster, so they spent them on doing more per cycle instead of more cycles per second.

There were three ways to do more per cycle: more cores (MIMD), a wider pipeline that finds independent instructions by itself (we meet this in Chapter 2), and wider data (SIMD). SIMD is the cheapest of the three. One instruction is fetched and decoded once, but it does the work of eight. Fetching and decoding instructions costs a lot of energy, so spreading that cost over many lanes is a good deal.

The price you pay is the rule from the picture: the lanes must do the same thing. Anything that does not fit that rule, such as different branches per element, dependencies between neighbours, or scattered memory, costs you. The rest of this book is mostly about the cost of that rule, and the tricks for getting around it.

What would have to be true for this design to be wrong? If most programs were full of irregular, branchy work and almost none were regular loops over arrays, wide registers would sit empty and designers would have chosen more cores instead. They did not, because images, audio, text, matrices and databases are mostly regular arrays of the same type.

Common mistakes

• Thinking SIMD is "automatic parallelism". It is parallelism inside one core and one thread. It does not use other cores.
• Thinking eight lanes means eight times faster. It means eight operations per instruction. Whether the program gets 8× faster depends on memory, dependencies and overhead. Chapters 2 and 5 are about exactly this.
• Assuming all CPUs have the same widths. They do not. Chapter 12 deals with it.

Say it simply

SIMD means one instruction works on a whole row of values stored side by side in one wide register. Each slot in that row is a lane. Lanes never talk to each other during a normal operation. It is a cheap way for a CPU to do more work per cycle, as long as the work is the same for every lane.

Check yourself

Quick questions

1. (Warm-up, predict) A register is 256 bits wide. How many lanes does it have when the values are u8? u16? f64? Write the three numbers down, then check against the table above.
2. (Warm-up, explain) In four sentences, explain "lane-wise" to a friend. What analogy would you use, and where does it break?
3. (Core, trace) Register a has u8 lanes [250, 3, 7, 9] and register b has [10, 4, 1, 200]. What are the four lanes of a + b with ordinary wrapping addition? (Hint: a u8 holds 0 to 255.)
4. (Core, explain) Why can a SIMD add use eight lanes of f32 but the same register hold only four lanes of f64?

Problems

5. (Warm-up, implement) Copy hello.rs. Change it to add two arrays of 16 floats by calling the 8-lane function twice. Check the answer against a plain loop.
6. (Core, break it) Remove #[target_feature(enable = "avx")] from add8 and compile. Read the error message. Does the compiler catch it, or does the program crash at runtime?

Think harder

7. (Stretch, design and defend) A friend says: "Instead of SIMD, just buy a CPU with more cores." When is your friend right? When are they wrong? Give one workload for each.

Chapter 2. The machine under your loop

In this chapter: you learn the two numbers that decide how fast any instruction can run (latency and throughput), why one chain of dependent operations is slow, and how to measure them on your own machine.

You need first: Chapter 1.

The problem first

Here is a loop that sums floats:

let mut s = 0.0;
for &v in x { s += v; }

You might guess: "an add is one cycle, so this runs at one element per cycle." On the machine used for this book it runs at about one element every two cycles, which is 5.9 GB/s. And a version of the same sum written with four vector accumulators runs at 129 GB/s: more than 20 times faster, on the same data, on the same core. (We measure it in Chapter 4. Chapter 5 shows what limits it next.)

Where does the factor of twenty come from? Not only from using wide registers. A big part is something quieter: waiting.

The model (the "what")

A modern CPU core is not a worker that does one instruction at a time. It is more like a kitchen with several stations. The cook (the scheduler) looks ahead at the next hundred or so instructions, finds ones that are ready, and hands them to free stations. The stations are called execution ports. An FMA (fused multiply-add) can run on two of them on this machine. A load can run on others.

Two numbers describe every instruction.

Latency is how long, in CPU cycles, from when an instruction starts to when its result is ready for the next instruction to use.

Throughput is how many of that instruction the core can start per cycle when there is nothing in the way. (Some books quote the reciprocal: "cycles per instruction". Both mean the same thing.)

These two are very different. Picture a pizza oven that takes 4 minutes to bake one pizza (latency), but has two shelves, so it can accept a new pizza every 30 seconds if you have pizzas ready (throughput). If each pizza needs the previous pizza as an ingredient, you wait the full 4 minutes every time. If you have 8 independent pizzas, the oven is always busy.

A dependent chain (each step needs the last result):

  cycle:  0    4    8    12   16
          |    |    |    |    |
  chain:  FMA--FMA--FMA--FMA--FMA        one FMA every 4 cycles
                                          the other oven shelves sit empty

Eight independent chains (no step needs another chain's result):

  chain 0: FMA----FMA----FMA----
  chain 1:  FMA----FMA----FMA---
  chain 2:   FMA----FMA----FMA--         two FMAs start every cycle once the
  ...                                    pipeline is full
  chain 7:        FMA----FMA---

Reading the picture: the work is identical. Only the arrangement differs. This idea has a name: instruction-level parallelism, or ILP. It is the CPU finding independent instructions and running them together. It is separate from SIMD (which is about wide data), and it multiplies with it.

The cost model. To estimate the time of a chain-heavy loop, use this:

time per iteration  ≈  max( latency of the longest dependency chain,
                            (number of instructions) / (instructions per cycle the ports allow) )

If one chain dominates, you are latency-bound. If there are many independent chains, you hit the port limit and you are throughput-bound. The cure for the first is more independent chains. The cure for the second is fewer instructions, which is where wider SIMD helps.

Real numbers. Here is what I measured on the machine for this book, for a fused multiply-add on 8 floats (256 bits, ymm) and on 16 floats (512 bits, zmm). "Chains" is how many independent accumulators the loop has.

instruction         chains  cycles/instr  instr/cycle  range of 5 runs
vfmadd231ps ymm          1          3.97         0.25    3.92 - 3.98  
vfmadd231ps ymm          2          1.98         0.50    1.95 - 2.00  
vfmadd231ps ymm          4          0.98         1.02    0.96 - 0.99  
vfmadd231ps ymm          8          0.51         1.95    0.51 - 0.53  
vfmadd231ps ymm         12          0.51         1.95    0.51 - 0.53  
vfmadd231ps zmm          1          3.92         0.26    3.81 - 3.97  
vfmadd231ps zmm          4          0.98         1.02    0.96 - 0.99  
vfmadd231ps zmm          8          0.58         1.72    0.58 - 0.59  
vfmadd231ps zmm         12          0.55         1.82    0.55 - 0.55

Read it slowly.

• With one chain, an FMA takes about 4 cycles. That is the latency. The chip can do 0.25 per cycle.
• With two chains, half that: 2 cycles per FMA. With four, 1 cycle. With eight, about 0.5 cycles.
• At 8 and 12 chains it stops improving at 0.5 cycles, which is 2 FMAs per cycle. That is the throughput limit: two ports can do FMAs.
• The table also says the answer to "how many accumulators do I need?" is latency × throughput = 4 × 2 = 8. Memorize this: you need roughly latency times throughput independent chains to keep a unit busy.
• With 512-bit registers the latency is still 4 cycles, and throughput comes out at about 1.8 per cycle. But each instruction carries 16 floats, not 8. So the work per cycle roughly doubles. (In one of three full runs the 512-bit rows were slower, about 0.85 cycles. I think another thread on the same host interfered, but I could not prove it.)

These match the published numbers for recent Intel cores, which you can look up in Agner Fog's instruction tables or on uops.info. You should still measure on your chip, because AMD and older Intel cores differ.

Build it (the "how")

How do you measure "cycles" when the clock speed changes under you?

The first instinct is to read the clock speed (say 3 GHz), time the loop, and multiply. I tried it. On this virtual machine the clock moved between about 1.6 and 3.0 GHz within one program run, so every conversion was wrong by a different amount.

The second instinct is to time a short loop that takes exactly one cycle per step, use it to learn the clock, then time the real loop. This is better, but it still breaks. My first "ruler" was a two-instruction loop. When another thread on the same physical core competed for the front end of the chip, the ruler ran at half speed. All my cycle counts then came out half as large. The table showed things like "3.2 FMAs per cycle", which is impossible, and that is how I knew.

The fix is a ruler that is hard to disturb, used in an alternating pattern:

1. Make the ruler a long chain of dependent integer adds: 32 of them per loop pass. Each add needs the last one's result, and an integer add takes exactly 1 cycle, so one pass is 32 cycles no matter what the clock does. Because it is long, it does not stress the front end.
2. Run the test loop for about half a millisecond. Then run the ruler for about half a millisecond. Repeat 101 times.
3. In each pair, the ratio (time per test instruction) / (time per ruler cycle) is the cycles per instruction. Clock drift cancels, because both halves of a pair see almost the same clock.
4. Take the median of the 101 ratios, and repeat the whole thing 5 times.

As a sanity check, the same ruler idea with a chain of imul instructions (known latency 3) gave imul/add = 3.05, 2.99, 2.83, ... across eight rounds. Close to 3, so the ruler is honest.

Here is the heart of it. First the ruler:

// src/bin/e1_latency.rs
/// Our "ruler": a chain of 32 dependent adds per loop pass. Each add needs the last
/// one's result, and an integer add takes exactly 1 cycle, so one pass = 32 cycles
/// whatever the clock speed is. It is long on purpose: a tiny loop can be slowed
/// down when another thread shares the core and competes for the front end.
fn ruler(passes: u64) -> f64 {
    let t = Instant::now();
    // SAFETY: plain integer asm, no memory access.
    unsafe {
        core::arch::asm!(
            "2:", add8!(), add8!(), add8!(), add8!(),
            "sub {n}, 1", "jnz 2b",
            x = inout(reg) 0u64 => _, y = in(reg) 1u64, n = inout(reg) passes => _,
            options(nomem, nostack)
        );
    }
    t.elapsed().as_nanos() as f64
}

Then the thing we test: a loop with K independent FMA chains. K is a const generic, so the compiler fully unrolls the inner loop and keeps all K accumulators in registers.

// src/bin/e1_latency.rs
#[target_feature(enable = "avx2,fma")]
fn fma256<const K: usize>(iters: u64) -> f32 {
    let (a, b) = (_mm256_set1_ps(black_box(0.999)), _mm256_set1_ps(black_box(0.001)));
    let mut acc = [_mm256_set1_ps(black_box(1.0)); K];
    for _ in 0..iters {
        for q in 0..K {
            acc[q] = _mm256_fmadd_ps(acc[q], a, b); // chain q depends on its own last value
        }
    }
    let mut s = _mm256_setzero_ps();
    for q in 0..K { s = _mm256_add_ps(s, acc[q]); }
    _mm256_cvtss_f32(s)
}

Now the trap that nearly ruined this experiment. My first version of that function had no black_box on the starting values. The program printed "0.00 cycles per FMA" and "171,627 instructions per cycle". That is nonsense.

What happened: the accumulators started at 1.0, and 1.0 × 0.999 + 0.001 rounds to exactly 1.0 in f32. The compiler worked out that the value never changes, so it deleted the whole loop. The benchmark measured nothing.

The fix is std::hint::black_box, which tells the optimizer "treat this value as unknown". That is why the code above wraps the constants in black_box(...). Remember this story. We will return to it in Chapter 20, because a benchmark that the compiler can see through is the most common way to fool yourself.

Why it is this way (the "why")

Why does the chip have both latency and throughput? Why not just make every instruction instant?

Because of physics and cost. An FMA involves multiplying and adding 24-bit numbers, and that takes several steps of logic. Hardware designers chop it into stages (a pipeline), so each stage is simple and fast, and a new operation can enter the first stage while the previous one is in the second. That gives high throughput, but any single operation still has to walk through all the stages. That walk is latency.

So the design bet is: programs usually have lots of independent work around. If that is true, the pipeline stays full and everything is fast. If your code is one long chain, like a running float sum, the bet fails and you pay for the full latency every step.

The lesson for SIMD: wide registers fix the width problem. They do nothing for the chain problem. A vector sum with one accumulator is still one chain, only eight lanes wide. You need both: wide registers and several independent accumulators.

Where the idea shows up elsewhere. The same thinking governs web servers (many requests in flight hide network latency), disks (queue depth), and GPUs (thousands of threads hide memory latency). "Keep enough independent work in flight to cover the delay" is a general law of fast systems.

Common mistakes

• Reading throughput from a chain. If you time a loop with one accumulator and divide, you measure latency and call it throughput. Always ask: how many chains?
• Benchmarking a loop the compiler can delete. If your result is "too good to be true", look at the machine code.
• Trusting the stated clock speed in a VM. Measure cycles with a ruler, or use a hardware counter if you can (perf stat).

Say it simply

Every instruction has a latency (how long until its answer is ready) and a throughput (how many can start per cycle). One long chain of dependent operations runs at latency speed. Many independent chains run at throughput speed. To keep a unit busy you need about latency × throughput independent chains. Wide SIMD registers do not remove the need for this.

Check yourself

Quick questions

1. (Warm-up, predict) An instruction has latency 3 and throughput 2 per cycle. How many independent chains do you need to keep it busy? Write your answer before checking the rule in the text.
2. (Core, trace) A loop does a = a * x + y on one accumulator, 1000 times, with FMA latency 4. About how many cycles? Now rewrite it with 4 accumulators of 250 steps each. About how many cycles? (Ignore the final combine.)
3. (Core, explain) In five sentences: why is the table above's one-chain row so much slower than the eight-chain row, even though the work is the same?
4. (Warm-up, explain) What is the difference between ILP and SIMD? Give a pizza-oven style analogy of your own.

Problems

5. (Core, predict then check) Run cargo run --release --bin e1_latency on your machine. Before you do, predict the cycle count for 1 chain and for 12 chains. Then compare. If your CPU is older or from AMD, what changes?
6. (Core, break it and fix it) Remove the black_box calls from fma256, compile, and run. What do you see? Look at the loop with objdump -d and explain, in two sentences, what the compiler did.
7. (Stretch, change one thing) Make a copy of fma256 for integer addition (_mm256_add_epi32, which avx2 provides). Predict its latency and throughput, then measure. Is the throughput 2 per cycle, or 3?

Think harder

8. (Stretch, design and defend) Your colleague wants to speed up a running sum of f32 by switching from SSE to AVX-512 (4× wider). You suspect it will barely help. Explain why using this chapter's cost model, and describe the cheapest experiment that would prove you right or wrong.

Chapter 3. Lanes, masks and shuffles

In this chapter: you learn the four moves every SIMD program is made of (lane-wise operations, comparisons that make masks, selection with a mask, and shuffles) by building a tiny simulator you can run on any computer.

You need first: Chapter 1.

The problem first

A scalar loop has if. This one clamps negatives to zero (a "ReLU" in neural networks):

for v in x { if *v < 0 { *v = 0; } }

Now eight values travel together in one register. Some are negative, some are not. There is only one instruction stream, so the CPU cannot go left for lane 1 and right for lane 2. What does if even mean here?

The answer SIMD gives is a little strange at first: don't branch. Compute both answers, then pick per lane. To pick per lane you need something that says, lane by lane, "yes" or "no". That something is a mask.

The model (the "what")

There are only four kinds of move. Learn them and you can read any SIMD code.

1. Lane-wise operations. Add, multiply, min, max, and so on. Lane i of the result depends only on lane i of the inputs. (Chapter 1.)

2. Comparisons make masks. x > 0 on eight lanes gives eight answers, not one. A mask holds one yes/no per lane. On x86 a lane that is "yes" has all its bits set to 1, and a lane that is "no" has all bits 0. On AVX-512 the mask lives in a special small mask register (k0 to k7) with one bit per lane instead.

3. Selection (blend). Given a mask and two vectors a and b, build a new vector: take a in lanes where the mask is yes, take b elsewhere. No branch. This is how if becomes data.

4. Shuffles (cross-lane moves). The only operation family where a lane looks at a different lane's value. "Reverse the order", "broadcast lane 3 everywhere", "pull the even lanes together" are all shuffles. Everything that is not lane-wise costs more, and shuffles are the main reason.

There is also a family that goes from lanes to a plain number: horizontal operations. "Sum all lanes" is one. We call ordinary lane-wise work vertical (the vertical direction is "between registers") and cross-lane work horizontal.

Here is a simulator. Nothing in it is fast. It exists so you can see the meaning of each move and run it on any machine.

// src/lanes.rs
//! A lane simulator: a vector is just N values that move together.
//! Nothing here is fast. It exists so you can see the *meaning* of vector operations.

#[derive(Clone, Copy, Debug, PartialEq)]
pub struct Vec8(pub [i32; 8]);

/// One bool per lane. On real hardware this is all-ones or all-zeros in each lane.
#[derive(Clone, Copy, Debug, PartialEq)]
pub struct Mask8(pub [bool; 8]);

impl Vec8 {
    pub fn splat(x: i32) -> Self {
        Vec8([x; 8])
    }
    /// Lane-wise add. Lane i of the result only looks at lane i of the inputs.
    pub fn add(self, o: Vec8) -> Vec8 {
        let mut r = [0; 8];
        for i in 0..8 {
            r[i] = self.0[i].wrapping_add(o.0[i]);
        }
        Vec8(r)
    }
    pub fn mul(self, o: Vec8) -> Vec8 {
        let mut r = [0; 8];
        for i in 0..8 {
            r[i] = self.0[i].wrapping_mul(o.0[i]);
        }
        Vec8(r)
    }
    /// Comparison does not return one bool. It returns a mask: one answer per lane.
    pub fn gt(self, o: Vec8) -> Mask8 {
        let mut r = [false; 8];
        for i in 0..8 {
            r[i] = self.0[i] > o.0[i];
        }
        Mask8(r)
    }
    /// Cross-lane work: lane i of the result is lane idx[i] of the input.
    pub fn shuffle(self, idx: [usize; 8]) -> Vec8 {
        let mut r = [0; 8];
        for i in 0..8 {
            r[i] = self.0[idx[i]];
        }
        Vec8(r)
    }
    /// Horizontal add, done as a tree: (0+1) (2+3) (4+5) (6+7), then pairs of pairs.
    pub fn reduce_add_tree(self) -> i32 {
        let a = self.0;
        let s1 = [a[0] + a[1], a[2] + a[3], a[4] + a[5], a[6] + a[7]];
        let s2 = [s1[0] + s1[1], s1[2] + s1[3]];
        s2[0] + s2[1]
    }
}

impl Mask8 {
    /// Pick lane by lane: where the mask is true take `a`, else take `b`. No branch.
    pub fn select(self, a: Vec8, b: Vec8) -> Vec8 {
        let mut r = [0; 8];
        for i in 0..8 {
            r[i] = if self.0[i] { a.0[i] } else { b.0[i] };
        }
        Vec8(r)
    }
    /// Squash the mask into the low bits of an integer. Lane 0 becomes bit 0.
    pub fn bitmask(self) -> u32 {
        let mut m = 0;
        for i in 0..8 {
            m |= (self.0[i] as u32) << i;
        }
        m
    }
}

Run it:

x + x       = [6, -2, 8, -2, 10, -18, 4, 12]
x > 0       = [true, false, true, false, true, false, true, true]
as bitmask  = 0b11010101  (lane 0 is the lowest bit)
relu(x)     = [3, 0, 4, 0, 5, 0, 2, 6]
reversed    = [6, 2, -9, 5, -1, 4, -1, 3]
sum (tree)  = 9

Look at the third line. x > 0 gave eight answers. The bitmask squashed them into the integer 0b11010101, with lane 0 as the lowest bit. That "mask to integer" move is called movemask on x86, and it is how SIMD code finally gets to ask a normal if: "is any lane a match?" (bitmask != 0) or "which lane was first?" (bitmask.trailing_zeros()). We use that trick in Chapter 11.

The fourth line is ReLU, done with no branch: select(x > 0, x, 0).

An analogy. Think of eight chefs working the same recipe on eight dishes. The head chef calls out one step at a time and everyone does it. To handle "if the sauce is too thin, add flour", the head chef calls "everyone add flour" and chefs whose sauce was fine hold their hands behind their back. Their hands being behind their back is the mask.

Where it breaks. The chefs who hold back still spend the time. This is the cost of mask-based "if": both sides of the if always run, for all lanes. If one side is very expensive and rarely needed, a branch (skip the whole vector when no lane needs it) can still win. Chapter 17 measures this.

Build it (the "how")

Let us map the four moves to real x86 instructions, so the simulator and the machine line up.

Two surprises are hiding in that table.

Surprise one: shuffles are restricted. _mm256_permutevar8x32_epi32 can move any lane to any lane across all 256 bits. But the very common byte shuffle, _mm256_shuffle_epi8 (called pshufb), only moves bytes within each 128-bit half. An AVX2 register is secretly two 128-bit halves glued together, and many instructions never cross the seam. This is the number one reason AVX2 code is harder to write than SSE code. We use pshufb for a very useful trick in Chapter 18.

Surprise two: horizontal sums are slow. There is no "add all lanes" instruction on x86. You build it from smaller steps: fold the top half onto the bottom, then fold pairs, and so on. Here is the real code that does it for 8 floats:

// src/sums.rs
/// Horizontal sum of 8 lanes: fold high half onto low half, then pairs.
#[target_feature(enable = "avx2")]
pub fn hsum256(v: __m256) -> f32 {
    let lo = _mm256_castps256_ps128(v);
    let hi = _mm256_extractf128_ps::<1>(v);
    let s = _mm_add_ps(lo, hi); // 8 lanes -> 4
    let s = _mm_add_ps(s, _mm_movehl_ps(s, s)); // 4 -> 2
    let s = _mm_add_ss(s, _mm_shuffle_ps::<0b01>(s, s)); // 2 -> 1
    _mm_cvtss_f32(s)
}

Eight lanes → four → two → one: three adds and some shuffles. That is why good SIMD code keeps values vertical in a vector for the whole loop and does the horizontal sum once, at the very end. If you do a horizontal sum inside the loop, you probably lost.

Why it is this way (the "why")

Why is a mask "all ones or all zeros" per lane, rather than a plain bool?

Because it makes selection cheap and general. If a mask lane is all ones, you can combine it with data using plain bitwise AND, OR, ANDNOT, with no special instruction: (mask AND a) OR (NOT mask AND b) is a select. The mask is just data, so you can store it, combine masks with AND/OR, and count them. The idea "control flow turned into data" is the central trick of SIMD, and later of GPUs too.

The design choices differ by instruction set, and they are worth knowing, because they show the trade-offs:

What would have to be true for the all-ones design to be wrong? If masks mostly went to memory instructions (to skip loads and stores), a separate one-bit-per-lane register would be clearly better. That is exactly what AVX-512 found, and it is why AVX-512's masked loads (Chapter 13) let it handle a loop's tail with no scalar code at all.

Common mistakes

• Treating a mask like a bool. In C, if (mask) makes no sense for eight lanes. You need movemask first to turn it into one integer.
• Using the wrong width for movemask. _mm256_movemask_ps gives one bit per 32-bit lane (8 bits). _mm256_movemask_epi8 gives one bit per byte (32 bits). Mixing them gives a mask that is "correct" but means something else.
• Forgetting the 128-bit seam. If your 256-bit shuffle gives garbage in the top half, you probably used a within-lane shuffle across the seam.
• Doing the horizontal sum inside the loop.

Say it simply

A vector holds many lanes. Most instructions work on each lane alone. A comparison gives a mask: one yes or no per lane. Selecting with a mask replaces if. Shuffles are the only moves that mix lanes, and they cost more. Keep work vertical and mix lanes as rarely as you can.

Check yourself

Quick questions

1. (Warm-up, trace) With the simulator's x = [3, -1, 4, -1, 5, -9, 2, 6], what is the bitmask of x.gt(Vec8::splat(2))? Do it by hand, then check by running.
2. (Core, trace) Write select for "abs": |x| using only gt, select and a way to negate. Which operations do you need?
3. (Warm-up, explain) Why does a mask hold all-ones or all-zeros per lane on AVX2, but one bit per lane on AVX-512? What can the second design do that the first cannot?
4. (Core, explain) Explain "vertical" and "horizontal" to a friend. Which is cheap, and why?

Problems

5. (Core, implement) Add Vec8::min(self, o) to the simulator, using only gt and select. Test it against i32::min on random data.
6. (Core, implement) Add Vec8::compress_to_front(self, mask) -> (Vec8, usize) to the simulator: put the lanes where the mask is true at the front, return how many. (Chapter 17 does this on real hardware. Doing it here first will make Chapter 17 easy.)
7. (Stretch, change one thing) In reduce_add_tree, the order of additions is a tree. Change it to a left-to-right running sum ((((a0+a1)+a2)+a3)...). For integers, does the answer change? What about for floats? (We return to this in Chapter 4.)

Think harder

8. (Stretch, design and defend) You have a loop over 1,000,000 elements. For 99.9% of them the "expensive branch" is not needed. Would you use select (compute both sides always) or an ordinary branch per vector (skip when no lane needs it)? What would change your answer?

Chapter 4. Numbers inside lanes

In this chapter: you learn what happens to integers and floats when a lane overflows, how to compare unsigned bytes on a CPU that only has signed compares, and why adding floats in a different order gives a different answer.

You need first: Chapters 1 and 3.

The problem first

You add 250 and 10 in a u8 lane. What is the answer?

In Rust, a plain 250u8 + 10 panics in a debug build. In a SIMD register there is no panic, no overflow flag, no check at all. The lane just does what the instruction says. Which instruction you pick is your answer to the question. This chapter is about reading that choice.

The model (the "what")

Integers. For a lane that overflows, hardware gives you separate instructions for separate meanings. Here are real results from running them:

250 + 10 wrapping   (paddb)  = 4
250 + 10 saturating (paddusb) = 255
250 + 10 widened to u16       = 260

is 200 > 100?  pcmpgtb (signed) says: NO, WRONG  (0xFF = yes, 0x00 = no)
after XOR with 0x80 on both sides:     yes, correct

f32: a = 1 + 2^-12
a*a - a*a            = 0e0   (both products rounded the same way, so exactly 0)
fma(a, a, -(a*a))    = 5.9604645e-8   (= the rounding error that the multiply threw away: 2^-24 = 5.9604645e-8)

Look at the first three lines. The same inputs, 250 and 10, give three different answers:

Image code lives on saturation (a pixel of 250 plus 10 is white, not black). Sums live on widening (a sum of a million bytes needs a bigger lane). Hash code likes wrapping. The instruction name already tells you the contract, so choose on purpose.

Signedness has a trap. Look at the second block of that output. The question "is 200 > 100 as unsigned bytes?" returned NO, WRONG. x86 SSE and AVX2 have a signed byte compare (pcmpgtb) and no unsigned one. Byte 200 is 0xC8. Read as a signed byte, that is −56. And −56 is not greater than 100.

The fix is a clean trick: XOR both sides with 0x80 (flip the top bit) first. Flipping the top bit shifts the number line so that unsigned order becomes signed order. Our compare then works, as the last line shows. You will meet this again in Chapter 17, where we filter u32 values with a signed-only 32-bit compare.

(Min and max are different: SSE2 and AVX2 do have unsigned byte min and max, _mm256_min_epu8, so clamp needs no trick. Always check the instruction list for the operation you need before assuming.)

Floats have a deeper trap: order. Floating-point addition is not associative. That means (a + b) + c can differ from a + (b + c). The reason: each add rounds its result to fit in 24 bits, and which numbers you add first changes which small bits get thrown away.

This matters a lot for SIMD, because a vector sum adds numbers in a different order than a plain loop. The compiler knows that, so it will not vectorize a float sum on its own: it is not allowed to change the order and change the answer.

Here is the same data summed five different ways. The first row of data is "nice" (multiples of 0.25, which float holds exactly) and the second is "nasty" (two huge values that cancel, plus many small ones):

kernel (16 KiB in L1)         ns / call         GB/s
naive loop                         2756          5.9
8 accumulators (scalar src)          342         47.9
AVX2, 1 accumulator                 323         50.8
AVX2, 4 accumulators                127        128.9
AVX-512, 4 accumulators              86        191.3
u32 wrapping (autovec)              325         50.4

results on 'nice' data (multiples of 0.25): 6142.5 6142.5 6142.5 6142.5 6142.5
exact (f64) sum of nasty data      : 532.2100
naive                               : 266.1308
8 accumulators                      : 499.09003
AVX2 4 accumulators                 : 524.0199
AVX-512 4 accumulators              : 528.20996

On the nice data all five agree: 6142.5. On the nasty data the exact answer (computed in f64) is 532.21. The five float answers: 266.13, 499.09, 524.02, 528.21. Every one of them is "a correct float sum", and every one is different. The naive loop is worst here, not best, because one running total of size 1e8 swallows the small numbers one at a time. Several accumulators keep each total smaller, so less is lost. The AVX-512 version, with 64 partial sums in flight (4 registers of 16 lanes), is the closest.

Fused multiply-add (FMA) is the other float surprise. fma(a, b, c) computes a×b + c with one rounding at the end, not two. The last block of the demo shows it: for a = 1 + 2⁻¹², a*a − (a*a) is 0 (two roundings that cancel), but fma(a, a, −(a*a)) is 5.96e-8, which is exactly 2⁻²⁴: the tiny error that the multiply had thrown away, now recovered. FMA is both faster (one instruction) and usually more accurate. But it gives different bits than multiply-then-add. If you need bit-exact agreement with a scalar loop that doesn't use FMA, you cannot use FMA.

Build it (the "how")

Time to see the cost of the chain from Chapter 2, now with numbers from the same program (the first table in the output above):

• the naive float loop: 5.9 GB/s
• 8 accumulators written in plain Rust, compiler vectorizes: 47.9 GB/s
• one AVX2 accumulator: 50.8 GB/s
• four AVX2 accumulators: 128.9 GB/s
• four AVX-512 accumulators: 191.3 GB/s
• integer sum (the compiler does vectorize it): 50.4 GB/s

Here are the first and the best version side by side.

// src/sums.rs
/// The obvious loop. Float addition is not associative, so the compiler
/// must add in exactly this order. That makes one long dependency chain.
pub fn sum_f32_naive(x: &[f32]) -> f32 {
    let mut s = 0.0f32;
    for &v in x {
        s += v;
    }
    s
}

// src/sums.rs
/// Four vector accumulators: four independent chains, so adds overlap.
#[target_feature(enable = "avx2")]
pub fn sum_f32_avx2_4acc(x: &[f32]) -> f32 {
    let end = x.len() / 32 * 32;
    let (mut a0, mut a1, mut a2, mut a3) = (
        _mm256_setzero_ps(),
        _mm256_setzero_ps(),
        _mm256_setzero_ps(),
        _mm256_setzero_ps(),
    );
    let mut i = 0;
    while i < end {
        // SAFETY: i + 32 <= end <= len
        unsafe {
            let p = x.as_ptr().add(i);
            a0 = _mm256_add_ps(a0, _mm256_loadu_ps(p));
            a1 = _mm256_add_ps(a1, _mm256_loadu_ps(p.add(8)));
            a2 = _mm256_add_ps(a2, _mm256_loadu_ps(p.add(16)));
            a3 = _mm256_add_ps(a3, _mm256_loadu_ps(p.add(24)));
        }
        i += 32;
    }
    let acc = _mm256_add_ps(_mm256_add_ps(a0, a1), _mm256_add_ps(a2, a3));
    hsum256(acc) + x[end..].iter().sum::<f32>()
}

What to notice. In the 4-accumulator version, the four vectors a0..a3 are four independent chains of the kind we built in Chapter 2. Every _mm256_add_ps on a0 depends on the last a0, but a0 and a1 do not depend on each other. That is why it is 2.5× faster than one accumulator even though both use 256-bit registers.

The first instinct, and why it breaks. The first instinct is "vectorize my float sum by having the compiler do it": write the loop, turn on optimizations, and hope. As the table shows, nothing happens: 5.9 GB/s. The compiler is not being lazy. It is being correct, because changing the order of additions would change your answer, and you did not give it permission. In Rust there is no -ffast-math to give blanket permission (this is on purpose). The honest fix is you choose a new order, and write it down: eight accumulators in plain Rust (sum_f32_8acc in the file) does exactly that, and the compiler can then vectorize it. We use this trick again in Chapter 9.

What this gives us. Choosing the order yourself does two things for us. It makes the result repeatable (the same input always gives the same sum, for this chunk size). And it removes the dependency chain, so the hardware can overlap the adds.

Why it is this way (the "why")

Why does the hardware leave overflow behaviour to the programmer instead of always checking?

Because checking costs. A vector instruction that also tracks "did any of 32 lanes overflow" would need extra hardware and a way to report which lane. SIMD designers chose: give the programmer a menu of simple instructions (wrap, saturate, widen) and let the software choose. The cost of that choice is that safety is now your job.

Why is float addition not associative? Because a float is a number with a fixed number of significant digits. 1e8 + 0.1 has no room for the 0.1: the result is just 1e8. Do 0.1 + 0.1 + ... (many times) first and you get something big enough to survive the addition to 1e8. Same numbers, different order, different answer. This is not a flaw of SIMD. It is a property of fixed-precision numbers, and SIMD simply makes it visible, because SIMD must reorder.

The deeper principle: a parallel program must say what it computes when the order is not fixed. Reductions are the first place where that rule bites. Whenever someone says "my parallel sum is not bit-identical to my serial sum", this is the reason. You can accept it, or fix the order, or use a more accurate method such as Kahan summation. But it must be a decision.

Common mistakes

• Assuming + means the same thing on a vector. It depends on the instruction: wrap, saturate, widen.
• Comparing unsigned bytes with a signed compare. The symptom is "wrong only for bytes above 127".
• Expecting bit-identical float results between scalar and SIMD. Compare with a tolerance, or fix the order and compare exactly, or accumulate in f64 for the check.
• Summing in a narrow lane. A u8 accumulator that sums bytes wraps after about 256 bytes of 255. The fix is to widen, or to flush the narrow accumulator periodically (we show this in Chapter 11).

Say it simply

A lane has a fixed size, so it can overflow. The hardware gives separate instructions for wrapping, saturating and widening, and you pick. x86 only compares signed bytes, so unsigned compares need a top-bit flip. Float addition depends on order, so a vector sum gives a slightly different answer than a scalar one, and the compiler won't reorder floats unless you do it yourself.

Check yourself

Quick questions

1. (Warm-up, predict) What is 200u8 as a signed byte? What is 0x80 ^ 200, and as a signed byte? Write both, then say why XOR with 0x80 fixes the compare.
2. (Core, trace) A u8 lane holds 250. You add 10 with a wrapping add, then 10 again. What is it after each step? Now do the same with a saturating add.
3. (Core, explain) In five sentences: why does the compiler turn the integer sum into SIMD by itself but not the float sum?
4. (Core, predict) In the nasty data, the exact sum is 532.21. Without running anything, which would you expect to be closer: summing in plain order, or summing as 8 independent accumulators? Why? (Use the order of operations to argue.)

Problems

5. (Core, implement) Write a scalar and an SSE2 sum_i16_saturating that adds a slice of i16 with saturation at each step. What goes wrong with the vector version? (Hint: saturating add isn't associative either.)
6. (Core, break it) Run e2_sums with the nasty data changed so the huge values are 1.0e4 instead of 1.0e8. Do the five answers still differ? By how much? What does that say about when you have to worry?
7. (Stretch, implement) Implement Kahan (compensated) summation in scalar Rust. On the nasty data, how close to 532.21 does it get? How much slower is it than sum_f32_naive? Can you vectorize it?

Think harder

8. (Stretch, design and defend) A numerical library promises "bit-identical results on every machine". Your SIMD sum would be 4× faster but differs in the last bits between AVX2 and AVX-512 machines. What do you do? Give at least two options and their costs.

Chapter 5. Memory is the real boss

In this chapter: you learn why a faster kernel often changes nothing, you meet the cache cliffs on a real machine, and you learn two small formulas (the roofline and Amdahl's law) that predict when SIMD will help.

You need first: Chapters 2 and 4.

The problem first

You wrote the fastest sum you know. Four accumulators, AVX-512. It ran at 191 GB/s on a small array. You are happy. Then you run it on your real data, a few hundred megabytes, and it runs at 13 GB/s. Same code. Fourteen times slower.

Nothing about the instructions changed. The data changed: it no longer fits in the fast memory close to the core.

The model (the "what")

Memory is not one flat thing. It is a stack of levels, each bigger and slower than the one above.

   level        size (this machine)       what it feels like
   ---------    -----------------------   -------------------------------
   registers    32 vector registers       instant
   L1 cache     48 KiB                    a few cycles
   L2 cache     2 MiB                     ~15 cycles
   L3 cache     tens of MiB, shared       ~50+ cycles
   main memory  GBs                       ~200+ cycles

(The cycle counts are rough, typical values for recent Intel server chips. I did not measure latency in this book. The bandwidth below I did measure.)

The unit that moves between levels is the cache line: 64 bytes on x86-64 and on most ARM chips (Apple's M-series chips use 128 bytes, as far as I know). Reading one byte pulls in the whole 64-byte line. This is great for sequential reads (the next 63 bytes come free) and terrible for scattered reads (you pull 64 bytes to use 4).

Bandwidth is bytes per second the memory system can deliver. Here is a real measurement of one core summing floats with our best AVX2 kernel, while the array grows:

  array size    ns / pass       GB/s   (AVX2, 4 accumulators, f32 sum)
       4 KiB           29      142.5
      16 KiB          125      130.6
      64 KiB          540      121.5
     256 KiB         3719       70.5
       1 MiB        17039       61.5
       4 MiB       156892       26.7
      16 MiB       632305       26.5
      64 MiB      2539690       26.4
     256 MiB     19795897       13.6

Read it like a landscape. Up to 64 KiB the array fits in L1 and the core reads about 120 to 140 GB/s. At 256 KiB it falls to ~70 GB/s: it now lives in L2. At 1 MiB, ~61. At 4 MiB and beyond it sits at about 26 GB/s: this is the next level out. At 256 MiB it drops to 13.6 GB/s, which is what one core can pull from main memory on this machine.

So the same kernel runs 10× slower when the data is big. SIMD made the arithmetic cheap, and now arithmetic is not the bottleneck at all.

(Caveat: in a virtual machine the sizes the OS reports for L3 may not be real. I'm trusting the shape of the curve, the plateaus and the cliffs, not the labels.)

Arithmetic intensity is the number that tells you which side of this you live on. It is the amount of useful work divided by the bytes you move:

arithmetic intensity  =  floating-point operations  /  bytes moved

The float sum does 1 add for every 4 bytes read: an intensity of 0.25 flop/byte.

The roofline. A simple formula, from Williams, Waterman and Patterson (2009), gives the best you can possibly do:

attainable speed  ≤  min( peak compute ,  memory bandwidth × arithmetic intensity )

Put in this machine's numbers. Peak compute for AVX2 FMA we measure in Chapter 15's experiment: about 96 GFLOP/s. Main-memory bandwidth for one core: about 13.6 GB/s.

float sum from main memory:  13.6 GB/s × 0.25 flop/byte = 3.4 GFLOP/s   (far below 96)

The sum is so far below the compute peak that the arithmetic hardware is idle 97% of the time. Making the arithmetic faster is pointless. Only moving fewer bytes can help. Compare a kernel that reuses data, like matrix multiply: for a block that stays in L1, each loaded number is used many times, so the intensity is high and the compute peak is the limit. We measure exactly that in Chapter 15.

Amdahl's law is the second formula. If a fraction p of your program's time gets s times faster, the whole program speeds up by:

speedup = 1 / ( (1 − p) + p / s )

Example: half the time in a loop you speed up 8×. Then 1 / (0.5 + 0.5/8) = 1.78. Not 8. The unchanged half sets the limit. A profile tells you p; guessing it is how people waste weeks.

Build it (the "how")

How do you know which regime you are in? Run the kernel on arrays of growing size and watch the curve, exactly as the experiment above did. Here is the whole program, short on purpose:

// src/bin/e3_bandwidth.rs
//! Experiment 3: how fast can ONE core read memory, as the array grows?
use simdbook::harness::*;
use simdbook::sums::x86;

fn main() {
    println!("{:>12} {:>12} {:>10}   (AVX2, 4 accumulators, f32 sum)", "array size", "ns / pass", "GB/s");
    let mut kib = 4usize;
    while kib <= 512 * 1024 {
        let n = kib * 1024 / 4;
        let data: Vec<f32> = (0..n).map(|i| (i & 15) as f32).collect();
        let s = measure(|| unsafe { x86::sum_f32_avx2_4acc(&data) });
        let label = if kib >= 1024 { format!("{} MiB", kib / 1024) } else { format!("{} KiB", kib) };
        println!("{:>12} {:>12.0} {:>10.1}", label, s.med_ns, gbps(n * 4, s.med_ns));
        kib *= 4;
    }
}

The recipe:

1. Pick a size that fits in L1 (16 KiB). Measure GB/s. This is the speed when memory is not the problem.
2. Pick a size far above your last-level cache (hundreds of MiB). Measure GB/s again. This is memory's ceiling.
3. If the second number is about the same as what a plain scalar loop gets there, SIMD will not help for that data size. Chapter 9 shows a case of this.

The first instinct is to "improve" the kernel when it runs slowly on big data: more unrolling, more accumulators. It breaks because the instructions were never the limit. The profile will say your core is waiting. You change what you read: fewer bytes (use u8 instead of f32 if the data allows), better reuse (process in blocks that fit in cache), or fuse two passes into one so you read once and do two things.

What that gives us. A habit: before writing any SIMD, compute the arithmetic intensity and compare it to the machine's ratio of compute to bandwidth. That ratio, 96 GFLOP/s ÷ 13.6 GB/s ≈ 7 flop/byte, is how many operations you must do per byte before compute becomes the limit. Almost every "simple loop over an array" is far below 7.

Why it is this way (the "why")

Why are there levels? Why not one big fast memory?

Because fast memory is physically big and power-hungry per bit. A small memory close to the core is fast because wires are short and the structure is simple. A large memory must be far away and uses slower, denser circuits. Every chip designer faces the same trade: speed, size, cost, pick two. The levels are a compromise built on one bet: programs reuse recent data and read nearby data. Computer scientists call these temporal locality (use it again soon) and spatial locality (use its neighbours soon).

Notice what this means for SIMD. SIMD reads wide chunks of neighbours, which is perfect for spatial locality. But SIMD is also so fast that it eats data faster than the next level can supply it. Wide registers made the arithmetic cheap and thereby made the memory wall easier to hit. This is a general pattern: speeding up one part of a system moves the bottleneck to the next part.

Where this shows up elsewhere. The roofline idea works for GPUs (peak FLOPs versus HBM bandwidth), for databases (CPU versus disk), and for networks (compute versus bytes on the wire). If you only learn one performance model, learn this one.

Common mistakes

• Benchmarking with data that fits in L1 and then being surprised in production.
• Expecting wider vectors to help a bandwidth-bound loop. They cannot.
• Forgetting that writes cost too. Storing to memory that is not in cache usually reads the line first (a "write-allocate"), so a copy moves more bytes than you think.
• Treating the L3 size in a VM as truth.

Say it simply

Memory is a stack of levels, from tiny and instant to huge and slow, and data moves in 64-byte lines. SIMD does so much arithmetic per cycle that on big data the core often waits for memory. A loop's arithmetic intensity (work per byte) and the machine's compute-to-bandwidth ratio tell you whether the arithmetic or the memory is the limit. If memory is, only moving fewer bytes helps.

Check yourself

Quick questions

1. (Warm-up, predict) A loop adds two f32 arrays into a third. How many bytes does it move per element (read, read, write, ignoring write-allocate)? How many flops? What is its arithmetic intensity?
2. (Core, trace) Using the roofline formula with 13.6 GB/s and 96 GFLOP/s, what is the best possible speed of the loop in question 1 from main memory? From L1 (use 130 GB/s)?
3. (Core, trace) Amdahl: 30% of your program's time is a loop you make 10× faster. What is the whole-program speedup? What if it's 70%?
4. (Warm-up, explain) In five sentences: what is a cache line, and why does reading one byte cost the same as reading 64?

Problems

5. (Core, predict then check) Run e3_bandwidth on your machine. Before it finishes, predict where the cliffs are (use lscpu for your cache sizes). How close were you?
6. (Core, implement) Write a version of e3_bandwidth that sums u8 instead of f32 using _mm256_sad_epu8 (we build this kernel in Chapter 11). Is the GB/s curve the same shape? Is the elements per second curve?
7. (Stretch, break it and fix it) Make a loop that is memory-bound: a[i] = b[i] * 2.0 over 256 MiB. Time it. Now change it to also compute a[i] = b[i] * 2.0 + c[i] using a third array. Predict the slowdown from the byte count alone. Is the real slowdown close to your prediction?

Think harder

8. (Stretch, design and defend) You must process a 4 GB file of f32 values twice: once to find the maximum, and once to subtract the maximum from every value. Design the passes. Could you do it in fewer trips through memory? What is the minimum number of times each byte must be read, and what limits you?

Part 1 wrap-up

What you can now do

• Draw a SIMD register, label the lanes, and say what "lane-wise" means.
• Explain latency, throughput, and why one chain of operations is slow.
• Say how many independent chains you need to keep a unit busy.
• Use a mask and a select instead of an if.
• Choose between wrapping, saturating and widening on purpose.
• Explain why float sums depend on order, and why the compiler won't vectorize them.
• Read a roofline and tell whether a loop is limited by compute or memory.

Mini project: your own cost card (2 to 3 hours)

Build a one-page "machine card" for your computer:

1. Run e1_latency and record FMA latency and throughput for 256-bit (and 512-bit if you have it).
2. Run e3_bandwidth and record the GB/s at L1, L2, and main memory sizes.
3. Compute your machine's balance: peak GFLOP/s divided by main-memory GB/s.
4. For three loops of your choice (a sum, a dot product, a copy), compute their arithmetic intensity and say which side of the roofline they sit on.

Done when: the card fits on one page, every number has the command that produced it, and you can say which of your three loops SIMD would help most.

Write-up: half a page. What did you expect? What did you see? What surprised you?

Reading list

• Williams, Waterman, Patterson. Roofline: An Insightful Visual Performance Model for Multicore Architectures. Communications of the ACM, 2009. Read this first: it is short and the picture explains everything.
• Agner Fog, Optimizing software in C++ and the Instruction tables (agner.org/optimize). Where the latency and throughput numbers live. Read the microarchitecture chapter for your CPU.
• uops.info (Abel and Reineke, ICS 2019 paper of the same idea). Look up VFMADD231PS for your CPU and compare to our table.
• Flynn, Very high-speed computing systems, Proceedings of the IEEE, 1966. Two pages on where "SIMD" and "MIMD" come from.
• Hennessy and Patterson, Computer Architecture: A Quantitative Approach. The chapters on memory hierarchy and data-level parallelism go deeper on everything in this Part.

Teach it back. In five sentences, explain to a friend why a loop that adds two big arrays gets almost no speedup from AVX-512, and what you would do instead.

Answer notes for Part 1 (try first!)

• Ch 1, Q1: 32, 16, 4. Q3: [4, 7, 8, 209]: 250+10 wraps to 4, 3+4=7, 7+1=8, 9+200=209.
• Ch 2, Q1: 3 × 2 = 6 chains. Q2: one chain: about 1000 × 4 = 4000 cycles. Four chains of 250: about 250 × 4 = 1000 cycles (the four overlap).
• Ch 3, Q1: x > 2 is [T, F, T, F, T, F, F, T], bitmask 0b10010101.
• Ch 4, Q1: 200 as signed is −56. 0x80 ^ 200 = 0x48 = 72, which as signed is +72. Subtracting 128 from every unsigned value maps it to the signed range with the same order, which is what the XOR does.
• Ch 5, Q1: 12 bytes, 1 flop, intensity 1/12 ≈ 0.083 flop/byte. Q2: 13.6 × 0.083 ≈ 1.1 GFLOP/s from main memory; 130 × 0.083 ≈ 10.8 from L1. Q3: 1/(0.7 + 0.03) = 1.37; 1/(0.3 + 0.07) = 2.70.

Part 2. Seeing SIMD shapes in your code

Part 1 gave you the machine. Part 2 teaches your eyes. Before you write a single intrinsic you should be able to look at a loop and say: this will vectorize easily, this needs a trick, or this should not be vectorized at all. You will also learn the single most important design decision: how your data is laid out in memory.

Chapter 6. Which loops are SIMD-shaped

In this chapter: you learn the handful of loop shapes that SIMD handles well, the shapes that fight it, and a decision tree for telling them apart in about ten seconds.

You need first: Chapters 1, 3 and 4.

The problem first

Most people's first idea is "SIMD speeds up any loop over an array." Try it on this one:

// Fragment, not run: an exponential moving average.
for i in 1..n {
    y[i] = 0.9 * y[i - 1] + 0.1 * x[i];
}

It is a loop over an array. It does simple float maths. It looks like a perfect target. And you cannot put four iterations in four lanes, because iteration i needs the result of iteration i − 1. Lane 1 would need lane 0's answer before lane 0 has finished. Lanes do not talk to each other.

This chapter is the skill of seeing that in advance.

The model (the "what")

Almost every loop you will meet falls into one of eight shapes. Learn the names. They are how experienced people talk.

Here is the picture to keep in your head. Each arrow is "this value is needed to compute that one". Look for arrows between neighbours:

 Map (no arrows between iterations):    Recurrence (a chain):

   x0   x1   x2   x3                      x0 -> y0 -> y1 -> y2 -> y3
   |    |    |    |                              ^     ^     ^
   v    v    v    v                              x1    x2    x3
   y0   y1   y2   y3

 lanes are free to run side by side      lanes would have to wait for each other

The test is this: can iteration i + 1 start before iteration i finishes? If yes, the loop is a candidate. This is exactly the question about independent chains from Chapter 2, asked about the loop instead of the instruction.

A loop-carried dependency is the name for "iteration i needs a value that iteration i − 1 produced". Maps have none. Recurrences are made of them. Reductions have one, but it is an associative one (adding is allowed to be regrouped, up to the float caveat), so we can split it into several chains. Scans are reductions that keep every partial result, which is why they are harder.

An analogy. Think of a car assembly line with eight identical stations. If each car is independent, eight cars go through in parallel. If car 2 cannot start until car 1 is finished (because car 2 uses parts that car 1 leaves behind), eight stations are useless. A recurrence is that second line.

Where it breaks. An EMA over one stream is a recurrence. But if you have eight different streams (eight audio channels, eight sensors) and each is its own EMA, then lane k can hold stream k, and they are independent. The recurrence is still there inside each lane, and it is fine, because no lane needs another lane. So the rule is not "recurrences are bad". The rule is "find independent work to put in the lanes". If there is none, SIMD has nothing to hold.

Build it (the "how")

Here is the decision tree. Use it every time you look at a loop.

                          Look at one loop
                                 |
              Does iteration i+1 need a result from iteration i?
                    /                                  \
                  NO                                   YES
                  |                                     |
     Does it write a different number          Is it just an accumulation
     of outputs than inputs?                   (sum, min, max, count)?
        /            \                              /               \
      NO             YES                          YES                NO
      |               |                            |                  |
   MAP / ZIP     FILTER or SEARCH            REDUCTION        Is there independent work
   (easy)        (compare+mask,              (use several     elsewhere to fill the lanes?
                  then compress or            accumulators)       /            \
                  movemask)                                     YES            NO
                                                                 |              |
                                                       Put independent       Probably leave
                                                       streams in lanes      it scalar

Two more questions to ask at every leaf:

• Is there a data-dependent if or early exit? Early exit is fine for search (check a whole vector, then exit). A different if per lane becomes a mask.
• Is memory the limit? (Chapter 5.) A map over a 4 GB array is limited by memory, not by lanes.

Real numbers for four of these shapes, all measured on the machine for this book:

The wins are not "8×". They are 5× to 30×, because several things stack: width, no chain, no branch, better layout.

The first instinct. "I'll vectorize the hot loop by hand." It breaks in two common ways. Either the compiler already vectorized it (check first, Chapter 9), so you spent a week to match what a flag gives you. Or the loop was a recurrence and you built something that is correct but no faster. The check that saves the week is: read the machine code the compiler made before writing any intrinsic.

Why it is this way (the "why")

Why are lanes independent? Why not let lane 1 see lane 0's result in the same instruction?

Because it would make the hardware serial again. A vector add is fast because all lanes finish in the same time, side by side, with no wires between them. If lane 1 had to wait for lane 0, you would have a chain of eight adders and an instruction that takes eight times as long. The designers chose "fast and independent" over "powerful and slow", and offered cross-lane work as separate, costlier instructions (shuffles). Everything else follows.

The first principle: parallel hardware only speeds up work that does not depend on itself. This is not a SIMD rule. It is a rule for GPUs, for multi-core, for distributed systems. Amdahl's law (Chapter 5) is the same fact seen from the other side: the part with a dependency is the part that stays serial.

What would have to be true for this to be wrong? If a clever algorithm could turn recurrences into independent work cheaply, SIMD would help them. Sometimes it can: a prefix sum can be computed in log n parallel steps instead of n serial ones, at the cost of doing more total additions. The trade (more work for shorter chains) is a real design lever, and Part 6 returns to it.

Common mistakes

• Assuming "has a loop" means "vectorizable".
• Missing a hidden dependency. a[i] = a[i-1] + b[i] is a scan. a[i + 1] += a[i] is too. Two loops that look different can be the same chain.
• Missing aliasing. If dst and src overlap, dst[i] = src[i + 1] may be fine for a vector or may not be. (Chapter 8.)
• Vectorizing something tiny. A loop over 12 items has no steady state. Setup eats everything.

Say it simply

Look at a loop and ask whether iteration i + 1 can start before iteration i is done. If it can, SIMD is a candidate. If each iteration needs the last one's result, there is nothing to put in the lanes, unless you can find independent streams. Maps, searches, filters and sums are the good shapes. Recurrences and scans are the hard ones.

Check yourself

Quick questions (these are code fragments for classification, not run)

1. (Warm-up) Name the shape of each: (a) out[i] = (a[i] > 0.0) as u32 (b) total += x[i] * x[i] (c) if x[i] == b'\n' { return i } (d) out[i] = lookup[idx[i] as usize] (e) y[i] = y[i-1] * 0.5 + x[i]
2. (Core, explain) In five sentences, explain what a loop-carried dependency is and why it blocks lanes.
3. (Core, trace) The loop a[i] = a[i] + a[i+1] (processing i from 0 upward) looks like it reads a neighbour. Is it a loop-carried dependency in the sense above? What changes if the loop runs from the end down?
4. (Core, predict) You have 8 independent audio channels, each needing a one-pole filter (y = a*y + b*x). Describe how you would lay the data out so the 8 channels go in the 8 lanes. What is the layout called?

Problems

5. (Core, implement) Write the scalar EMA from the problem above, and write a version that processes 8 independent streams held as [f32; 8] per time step (a plain Rust array of 8, no intrinsics). Run both on the same total number of samples. Compare time. (The compiler will vectorize the second one. Use objdump to prove it.)
6. (Stretch, break it) Take any loop in a project you own. Predict whether the compiler vectorizes it. Check with objdump -d (look for ymm/xmm packed instructions like vaddps or paddd, not addss). What was blocking it if not?

Think harder

7. (Stretch, design and defend) A prefix sum over 1 million u32 takes ~1 ms scalar. A log-step SIMD version does about 3× more additions but has no long chain. Would you expect it to be faster or slower, and what would you measure first to find out?

Chapter 7. Data layout decides everything

In this chapter: you learn the three ways to arrange a collection of records in memory (AoS, SoA, AoSoA), why layout often matters more than the instructions, and you see a 4× to 10× difference come from layout alone.

You need first: Chapters 5 and 6.

The problem first

You have a particle system. A particle has a position, a velocity, a mass and an id. The natural Rust is this:

struct Particle { x: f32, y: f32, z: f32, vx: f32, vy: f32, vz: f32, mass: f32, id: u32 }
let particles: Vec<Particle> = ...;

Now the update step: x += vx * dt for every particle. It is a perfect map. So SIMD should fly. It does not, and the reason is memory.

The model (the "what")

There are three layouts. Look at where the x values sit in each.

AoS  "Array of Structs"  (what Vec<Particle> gives you)

 address: 0    4    8    12   16   20   24   28 | 32   36   40  ...
          x0   y0   z0   vx0  vy0  vz0  m0   id0| x1   y1   z1  ...
          ^                                        ^
          x values are 32 bytes apart. A 32-byte vector load grabs only ONE x.

SoA  "Struct of Arrays"

 x:   x0  x1  x2  x3  x4  x5  x6  x7  x8 ...      one vector load = 8 x's
 y:   y0  y1  y2  ...
 vx:  vx0 vx1 vx2 ...                             every array is dense

AoSoA  "Array of Struct of Arrays"  (blocks of 8)

 block 0:  x0..x7 | y0..y7 | z0..z7 | vx0..vx7 | ...     <- one block of 8 particles
 block 1:  x8..x15| y8..y15| ...

AoS keeps everything about one particle together. That is natural for humans, and good when you work on one particle at a time and need all its fields. But the xs are scattered. To fill one vector register with 8 xs, the CPU has to gather them from 8 different 32-byte records, which is a shuffle-heavy mess. And every cache line also drags along fields you did not need.

SoA keeps each field in its own dense array. Eight consecutive xs are one clean vector load. Nothing unneeded rides along.

AoSoA is the compromise: groups of 8 (or 16) particles stored SoA inside the group, with the groups laid out one after another. You still get clean vector loads, and the fields of one group sit close together (good for the cache).

Here is the measurement. I wrote the update in plain Rust twice, once for each layout, and let the compiler vectorize. "all3" updates x, y, z. "x-only" updates just x, to show the worst case for AoS. Sizes are chosen to land in L1/L2 (4,096 and 65,536 particles) and in main memory (4 million).

   particles |       all3 AoS       all3 SoA |     x-only AoS     x-only SoA   (ns per particle)
  4096 (  0MB) |          1.179          0.416 |          0.661          0.114
 65536 (  2MB) |          1.336          0.633 |          1.020          0.195
4000000 (122MB) |          4.981          1.023 |          2.749          0.284

Read the rows.

• In cache (4,096 particles): AoS 1.18 ns, SoA 0.42 ns. SoA is 2.8× faster, purely from the layout.
• In main memory (4 million, 122 MB), updating all three: AoS 4.98 ns, SoA 1.02 ns: 4.9× faster.
• In main memory, updating only x: AoS 2.75 ns, SoA 0.28 ns: 9.7× faster. The AoS version loads 32 bytes per particle to use 4 of them, so it moves eight times more data than the SoA version. SoA's extra advantage over the ratio of bytes comes from vectorization on top.

So the same arithmetic ran between 2.8× and 9.7× faster. No intrinsic was written. The layout did it.

Here are the two functions. Notice they are both ordinary safe Rust.

// src/layout.rs
pub fn step_aos(p: &mut [ParticleAos], dt: f32) {
    for q in p {
        q.x += q.vx * dt;
        q.y += q.vy * dt;
        q.z += q.vz * dt;
    }
}

// src/layout.rs
pub fn step_soa(p: &mut ParticlesSoa, dt: f32) {
    for (a, v) in p.x.iter_mut().zip(&p.vx) {
        *a += *v * dt;
    }
    for (a, v) in p.y.iter_mut().zip(&p.vy) {
        *a += *v * dt;
    }
    for (a, v) in p.z.iter_mut().zip(&p.vz) {
        *a += *v * dt;
    }
}

Build it (the "how")

Get into the habit of asking: which fields does this loop touch? Draw the cache lines. In AoS with a 32-byte particle, a 64-byte cache line holds 2 particles. Reading only x for 100 particles pulls in 50 cache lines (3,200 bytes) to use 400 bytes. In SoA it pulls 7 lines. That 7× in bytes moved is the 9.7× we saw.

The first instinct is to keep the AoS layout and fix it with intrinsics: use gathers or shuffles to deinterleave. It is possible, and compilers sometimes do it. It works, but you spend instructions on shuffles that SoA never needs, and you still read the unused bytes. The better fix is to change the data. If you cannot (because other code expects Vec<Particle>), a cheap trick is to convert at the boundary: keep AoS for storage and convert to SoA blocks for the hot loop, if the loop does enough work per element to pay for the conversion.

What the change gives us. It makes SIMD loads simple, and it cuts bytes moved. Two wins from one change.

A real cost: SoA is awkward when you need all fields of one record (look up particle 17 and print it: 6 loads from 6 arrays, 6 cache misses instead of 1). It is also more code to push and remove. So it is a trade, not a law. AoS wins for point lookups, SoA for bulk loops.

Why it is this way (the "why")

The deep point is that a data layout is a decision about what "one thing" means. In AoS the unit is a particle. In SoA the unit is a field across all particles. Neither is wrong. They answer different questions.

The same idea appears far outside SIMD:

• Databases. Row stores keep each row together (AoS). Column stores keep each column together (SoA). Column stores are far better for "sum one column across a billion rows" and worse for "fetch one row". The C-Store paper (Stonebraker et al., VLDB 2005) argues exactly this, and nearly every analytics database since follows it.
• Game engines. Entity-component-system designs store each component type in its own dense array, for the same reason.
• GPUs. Memory is fastest when neighbouring threads read neighbouring addresses. That is SoA again.

The principle: the cheapest data movement is the one you do not do, and you do not do it by storing together what is used together.

What would have to be true for AoS to be better? If nearly every operation needs all fields of one record, and records are looked up individually (a lot of random access, little bulk work), AoS keeps one record in one cache line and wins.

Common mistakes

• Writing intrinsics before looking at layout. Fix the layout first, then see if the compiler is already done.
• Measuring only in-cache. The layout win is biggest when data is large (above: 2.8× in cache, 9.7× in memory for x-only).
• Making everything SoA. Point lookups get slower. Choose per hot loop.
• Padding fields carelessly. A 7-field struct padded to 32 bytes wastes memory bandwidth in AoS.

Say it simply

Memory layout decides how many bytes a loop pulls through the cache and whether a vector load can be one clean instruction. Array of Structs is natural but scatters each field. Struct of Arrays packs each field densely, so vector loads are simple and no unused bytes ride along. Pick the layout that matches the loop that runs the most.

Check yourself

Quick questions

1. (Warm-up, predict) A record is 24 bytes. A loop reads only one 4-byte field of 1,000,000 records stored AoS. About how many bytes of cache lines does it pull in? How many does the same loop pull in from an SoA array of that field? Use 64-byte lines.
2. (Core, explain) In five sentences, explain AoS, SoA and AoSoA to a friend. What analogy would you use? (Hint: library shelves sorted by book, or by topic.)
3. (Core, trace) In AoSoA with blocks of 8, where is the vy of particle 19? Give the block number and the offset within the block, assuming 6 fields of 8 f32 each per block.
4. (Core, explain) Why is the x-only AoS case worse relative to SoA than the all-fields case?

Problems

5. (Core, predict then check) Run e6_layout. Predict the AoS/SoA ratio for x-only at 65,536 particles before looking at the table. (L2 resident.) How close were you?
6. (Core, implement) Write step_aosoa: store particles in blocks of 8 (struct Block { x: [f32; 8], vx: [f32; 8], ... }) and write the update as an inner loop over the 8. Does the compiler vectorize it? How fast is it compared with AoS and SoA?
7. (Stretch, break it and fix it) Add a mass-weighted kinetic energy sum: sum += 0.5 * mass * (vx*vx + vy*vy + vz*vz) to both the AoS and SoA versions. Which is faster now? Which fields does each version have to read? Does the gap shrink or grow compared with the plain update?

Think harder

8. (Stretch, design and defend) Your team's code base uses Vec<Particle> everywhere. A profile says one function that updates positions takes 40% of the runtime. Argue for either (a) rewriting the whole code base to SoA, or (b) converting at the boundary of that one function. What would you measure before deciding?

Chapter 8. Tails, alignment and aliasing

In this chapter: you learn the three boring-sounding details that cause most real SIMD bugs: what to do with the leftover elements at the end of a loop, what alignment is and when it matters, and what happens when two slices overlap.

You need first: Chapters 3 and 6.

The problem first

A vector loop processes 32 bytes at a time. Your slice has 100 bytes. 32 + 32 + 32 = 96. What happens to the last 4 bytes?

The beginner answer is "read 32 bytes anyway and ignore the extra". That feels harmless: the extra bytes are just memory. It is a crash waiting to happen, and in Rust it is also undefined behaviour, meaning the compiler is allowed to assume it never happens and may produce wrong code. This chapter is about getting the end of the loop right, which sounds boring and is not.

The model (the "what")

The tail is the leftover elements after the last full vector. With width W and length n, there are n mod W of them (0 to W − 1).

There are four honest ways to handle a tail. Each has a correctness condition:

You can see where the rule "never read past the end of a slice" comes from. A slice &[u8] of length 100 promises 100 readable bytes. Reading byte 100 is outside that promise. It may land in another allocation, in unmapped memory (a crash), or work "by luck" today and fail after an unrelated change. Rust's rules treat it as undefined behaviour even if the hardware would not fault.

Our x86 code uses the first strategy by default (h[end..].iter().position(...)) and the AVX-512 code uses the third. The overlapping strategy is easy to get wrong, so here is a working one. Look at how the last vector is placed:

// src/bytes.rs
/// No scalar tail: the last vector is placed so it ENDS at the end of the slice,
/// which means it overlaps the previous one. That is only legal because clamp is
/// idempotent: clamping an already-clamped byte changes nothing.
/// (For a sum, the overlap would count bytes twice. Never copy this trick to `sum`.)
#[target_feature(enable = "avx2")]
pub fn clamp_u8_avx2_overlap(v: &mut [u8], lo: u8, hi: u8) {
    if v.len() < 32 {
        super::clamp_u8_scalar(v, lo, hi);
        return;
    }
    let (lo_v, hi_v) = (_mm256_set1_epi8(lo as i8), _mm256_set1_epi8(hi as i8));
    let one = |p: *mut u8| {
        // SAFETY: the caller passes a pointer with 32 valid bytes after it.
        unsafe {
            let b = _mm256_loadu_si256(p.cast());
            _mm256_storeu_si256(p.cast(), _mm256_min_epu8(_mm256_max_epu8(b, lo_v), hi_v));
        }
    };
    let last_start = v.len() - 32;
    let mut i = 0;
    while i < last_start {
        one(v[i..].as_mut_ptr());
        i += 32;
    }
    one(v[last_start..].as_mut_ptr()); // the overlapping final vector
}

 100 bytes, vectors of 32:
 [ 0 .. 32 )  [ 32 .. 64 )  [ 64 .. 96 )
                                    [ 68 .. 100 )     <- the final vector ENDS at 100
                                     ^^^^^^^^
                         bytes 68..96 are processed twice

For clamp this is fine: clamping an already-clamped byte changes nothing. For sum it is a bug: bytes 68 to 95 would be added twice. The overlap trick works only for idempotent operations. The same applies to anything that is not "set to a function of the input", for example x[i] += 1 (not idempotent) and x[i] = min(x[i], 9) (idempotent).

Alignment is about the address of the first byte. A memory address is aligned to 32 if it is a multiple of 32. Older CPUs had aligned loads that were fast and unaligned loads that were slow, and some instructions (_mm256_load_si256, without the u) crash on an unaligned address. On modern x86 chips (about the last ten years), an unaligned load that stays inside one cache line is as fast as an aligned one. The real cost is a load that straddles two cache lines: a 32-byte load starting at byte 48 of a line crosses into the next line and needs two cache accesses. This is a cache-line split.

So the working rule is: use loadu and storeu (the unaligned versions) by default. They accept any address, and they are fast when the address happens to be aligned. Use aligned instructions only when you control the allocation and have measured a benefit. In Rust, a Vec<u8> does not promise 32-byte alignment. To get it you need a type like #[repr(align(32))] wrapped in an array.

Aliasing is two pointers pointing at overlapping memory. Think of dst[i] = src[i + 1]. If dst and src are the same slice, the result depends on the order the loop runs: element i is written before element i + 1's old value would have been read by the next step. A vector version reads 8 src values first, then writes 8 dst values. For a copy-with-shift that changes the answer.

Rust helps a lot here. A function taking dst: &mut [f32] and src: &[f32] can never receive overlapping slices, because the borrow checker forbids a &mut that overlaps any other reference. LLVM knows this (Rust passes the information as noalias), so it can vectorize loops that a C compiler must fear. That is a quiet advantage of Rust for SIMD. The danger returns the moment you use raw pointers or unsafe.

Build it (the "how")

How do you prove you got the tail right? Not by testing one length. By testing every length around the vector width, at every alignment. That is called a differential test: run the SIMD kernel and the scalar oracle on the same input and demand the same answer.

Here is the test from this book's crate. It checks find_byte across all kernels for every length from 0 to 300, at five different start offsets, with needles that do and do not occur:

// tests/differential.rs
/// A buffer whose first byte is 64-byte aligned, so `offset` controls alignment exactly.
fn aligned_bytes(len: usize, seed: u64) -> (Vec<u8>, usize) {
    let mut v = vec![0u8; len + 128];
    let mut r = Rng(seed);
    for b in v.iter_mut() {
        *b = (r.next() >> 24) as u8;
    }
    let off = (64 - (v.as_ptr() as usize % 64)) % 64;
    (v, off)
}

// tests/differential.rs
#[test]
fn find_byte_all_kernels() {
    assert!(is_x86_feature_detected!("avx2") && is_x86_feature_detected!("avx512bw"));
    for len in 0..=300 {
        for offset in [0usize, 1, 7, 31, 63] {
            let (buf, base) = aligned_bytes(len, 1 + len as u64);
            let h = &buf[base + offset..base + offset + len];
            let mut cases: Vec<u8> = vec![0x00, 0xFF, 0x80, 0x7F];
            cases.extend(h.iter().step_by(11).copied()); // needles that DO occur
            for needle in cases {
                let want = bytes::find_byte_scalar(h, needle);
                assert_eq!(bytes::find_byte_swar(h, needle), want, "swar len={len} off={offset}");
                // SAFETY: we asserted avx2 and avx512bw above (sse2 is in every x86-64 CPU).
                unsafe {
                    assert_eq!(bytes::x86::find_byte_sse2(h, needle), want, "sse2 len={len} off={offset}");
                    assert_eq!(bytes::x86::find_byte_avx2(h, needle), want, "avx2 len={len} off={offset}");
                    assert_eq!(bytes::x86::find_byte_avx2_unrolled(h, needle), want, "unrolled len={len}");
                    assert_eq!(bytes::x86::find_byte_avx512(h, needle), want, "avx512 len={len} off={offset}");
                    assert_eq!(bytes::x86::count_byte_avx2(h, needle), bytes::count_byte_scalar(h, needle));
                }
            }
        }
    }
}

And the result:

test classify_every_byte_value ... ok
test csv_scanners ... ok
test filter_all_kernels ... ok
test find_byte_all_kernels ... ok
test floats_and_matmul ... ok
test int8_dot_small_values_agree_big_values_do_not ... ok
test sum_and_clamp ... ok
test result: ok. 7 passed; 0 failed; 0 ignored; 0 measured; 0 filtered out; finished in 0.02s

Why those lengths and offsets? Because the bugs hide at the edges: length 0, length 1, length W − 1, W, W + 1, 2W − 1, and so on. And the offsets 1, 7, 31, 63 make the first vector start at every kind of alignment, including ones that straddle cache lines.

The first instinct when a test fails only at "weird" lengths is to suspect the vector loop. The bug is almost always the tail or the loop bound. For example, while i <= len - W underflows when len < W (an unsigned subtraction wraps to a huge number). The safe form is while i + W <= len. Our code uses len / W * W computed once.

What the test gives us. It turns "I think it is right" into "it agrees with the oracle on 5 × 301 × many inputs." That is evidence, not hope. Chapter 21 builds this into a full verification plan.

Why it is this way (the "why")

Why does Rust care so much about reading one byte past the end?

Because the compiler's optimizations rest on a promise: a slice of length 100 has exactly 100 bytes you may touch. The optimizer uses that promise to reorder, merge and delete memory accesses. A read outside the slice breaks the promise, and then every optimization based on it can go wrong in ways that look random. This is the same reason it is "unsafe" and not just "slow": the problem is not the hardware's reaction, it is the meaning of your program.

A real-world twist: C libraries such as glibc's strlen do read past the end of a string, in a controlled way. They read aligned blocks, and an aligned block can never straddle a page boundary, so the extra bytes are guaranteed to be in mapped memory. That is a deliberate, hardware-aware hack that works in assembly. In Rust, the abstract machine rules still call it undefined behaviour, so Rust SIMD libraries avoid it, or hide it very carefully behind a masked load that does not read the masked-off bytes. We test that claim on real hardware in Chapter 13.

The design lesson. There are two ways to guarantee the tail is safe. You can make the hardware skip the invalid lanes (masks). Or you can make the software only touch valid memory (scalar tail, overlap). Hardware masks are cleaner, which is one of the strongest arguments for AVX-512, SVE and RISC-V's vector design.

What would have to be true for "always use loadu" to be wrong? On old CPUs (before about 2010), unaligned loads were much slower and aligned data mattered a lot. If you must support them, or you have found cache-line splits dominating a hot loop, aligning the start of your data helps.

Common mistakes

• Reading len bytes with a full-width vector at the end. Undefined behaviour, and a crash at the end of a page.
• Underflow in the loop bound. Use i + W <= len.
• Using the overlap trick on a non-idempotent operation. (sum, count, += 1).
• Aligned instructions on unaligned Vecs. Vec<u8> is not 32-byte aligned.
• Calling unsafe code with overlapping slices. Once you leave safe Rust, nobody checks.

Say it simply

The tail is the leftover after the last full vector. Handle it with a scalar loop, an overlapping final vector (only for idempotent operations), a masked operation, or padding, and never read past the end of the slice. Use unaligned loads by default. Rust's &mut rules rule out overlap for you, as long as you stay in safe code. Prove it by testing every length and offset against a scalar oracle.

Check yourself

Quick questions

1. (Warm-up, trace) A slice has 100 bytes and the vector width is 32. How many full vectors? How many tail bytes? Where does the overlapping final vector start?
2. (Core, explain) Which of these can safely use the overlapping-tail trick: x[i] = x[i].min(9), x[i] += 1, x[i] = !x[i] (bitwise not), out[i] = f(in[i]) into a different buffer? Explain each in one sentence.
3. (Core, explain) What is a cache-line split, and why does loadu usually not hurt even though the address is unaligned?
4. (Warm-up, explain) Why can Rust vectorize fn f(dst: &mut [f32], src: &[f32]) more freely than a C compiler can vectorize the pointer version?

Problems

5. (Core, break it and fix it) Here is a buggy tail. Find the bug without running it, then prove it with a test:

   // Fragment: find the sum of bytes. Buggy on purpose.
   let mut i = 0;
   while i <= data.len() - 32 { /* add one 32-byte block */ i += 32; }

   (Hint: what is data.len() - 32 when data.len() is 10?)

6. (Core, implement) Write a differential test for clamp_u8_avx2_overlap over lengths 0 to 100. Then deliberately change clamp to add 1 (non-idempotent) and watch the test catch the overlap bug.

7. (Stretch, implement) Make a 32-byte-aligned buffer using #[repr(align(32))] and write a kernel using _mm256_load_si256. What happens if you pass it buffer + 1? (Do this inside a child process, since it will crash.)

Think harder

8. (Stretch, design and defend) You are designing the public API for a function sum_bytes(&[u8]). Choose how the tail is handled: scalar, masked, or "require the caller to pad". Which do you pick for a library used by strangers, and which for a private function inside one project? Why do they differ?

Part 2 wrap-up

What you can now do

• Look at a loop and name its shape within ten seconds.
• Find loop-carried dependencies, and find independent streams to put in lanes.
• Compare AoS, SoA and AoSoA, and pick one for a given hot loop.
• Handle a tail four ways and say which is legal when.
• Explain why unaligned loads are the default and when alignment matters.
• Write a differential test that sweeps lengths and offsets.

Mini project: the pattern classifier (3 to 4 hours)

Take ten loops from a real code base (yours, or any open-source Rust project). For each, write down: its shape, its loop-carried dependencies, which layout it uses, whether the compiler vectorized it (check with objdump), and your prediction for how much SIMD would help. Then pick the most promising and measure the AoS-versus-SoA difference for it.

Done when: a table with ten rows, each with evidence (the assembly line or the benchmark number), and one measured layout comparison.

Write-up: half a page. How many of your predictions held? Which surprised you?

Reading list

• Intel, Intel 64 and IA-32 Architectures Optimization Reference Manual. The sections on SIMD and alignment are the primary source for the claims about loads in this Part.
• Stonebraker et al., C-Store: A Column-oriented DBMS. VLDB 2005. The same layout argument, in databases.
• The Rust Reference, section on behavior considered undefined, and the standard library docs for core::arch (for the safety contracts of each intrinsic).
• Mike Acton, Data-Oriented Design and C++, CppCon 2014 talk. A strong, practical argument for SoA-style thinking.

Teach it back. In five sentences, explain to a friend why Vec<Particle> can be four to ten times slower than separate arrays for a bulk update, even when the arithmetic is identical.

Answer notes for Part 2

• Ch 6, Q1: (a) map (with a compare), (b) reduction, (c) search / early exit, (d) gather, (e) recurrence. Q3: No loop-carried dependency going upward: iteration i reads a[i+1] which has not yet been written (it holds its original value) and writes a[i]. That is a "read ahead" pattern the compiler can vectorize safely. Going downward it would read a value just written, which is a dependency.
• Ch 7, Q1: AoS: every 24-byte record spans 1 to 2 lines; about 1,000,000 × 24 bytes ≈ 24 MB of lines ≈ 375,000 lines. SoA: 4 MB ≈ 62,500 lines. About 6× fewer. Q3: particle 19 is in block 2 (19 ÷ 8 = 2 remainder 3). Lane 3 of its vy array. The vy array is the 5th array in a block, so its element is at offset (4 × 8 + 3) × 4 bytes if arrays are ordered x, y, z, vx, vy, vz (index 4 for vy).
• Ch 8, Q1: 3 full vectors (96 bytes), 4 tail bytes, final vector starts at 68. Q2: min(9) is idempotent: yes. += 1: no. !x: no (applying it twice gives back the original). out = f(in) into a different buffer: yes, because the input is never modified, so re-reading it gives the same result.

Part 3. Writing it in Rust

Now we write code. Each chapter in this Part adds one layer: first let the compiler try, then use std::arch by hand, then build real byte kernels, then choose between them at runtime, then use AVX-512, then look at other platforms. By the end you can write, test and ship a SIMD kernel.

Chapter 9. Let the compiler try first

In this chapter: you learn what the Rust compiler can vectorize by itself, what stops it, and how to check in one minute with objdump. Often this is the whole job.

You need first: Chapters 4 and 6.

The problem first

You have a hot loop. You are about to spend a day writing intrinsics. Before you do, there is a one-minute check that decides whether that day is needed. The compiler (LLVM, inside rustc) has an auto-vectorizer: a pass that turns plain scalar loops into SIMD loops. It is good. It is also limited in ways you can learn.

The model (the "what")

The auto-vectorizer has two main tools. The loop vectorizer takes a loop and runs several iterations per vector step. The SLP vectorizer (for superword-level parallelism) takes a few similar scalar statements next to each other and merges them into one vector statement.

It will vectorize a loop when it can prove all of these:

1. The iterations are independent (Chapter 6), or the dependency is a reduction it is allowed to reorder.
2. It is allowed to change the order of the operations (this is why integer sums vectorize and float sums do not, Chapter 4).
3. Memory accesses cannot overlap in a harmful way (Rust's &mut helps, Chapter 8).
4. The trip count and bounds checks do not get in the way. A bounds check inside the loop is a branch that can exit early. The compiler can often prove it away, but not always.
5. The CPU features exist. This is the big one. The default Rust target for x86-64 allows only the baseline, and here is what that means in practice:

target_feature="fxsr" target_feature="sse" target_feature="sse2"

That is the whole list: fxsr, sse, sse2. No AVX, no FMA, no popcount. So by default the auto-vectorizer can only use 128-bit registers, even if your CPU has 512-bit ones. This is deliberate: a binary built by default must run on every x86-64 CPU. Here is what a build with -C target-cpu=native allows on the machine used for this book:

target_feature="avx" target_feature="avx2" target_feature="avx512bf16" target_feature="avx512bitalg" target_feature="avx512bw" target_feature="avx512cd" target_feature="avx512dq" target_feature="avx512f" target_feature="avx512fp16" target_feature="avx512ifma" target_feature="avx512vbmi" target_feature="avx512vbmi2" target_feature="avx512vl" target_feature="avx512vnni" target_feature="avx512vpopcntdq" target_feature="avxvnni" target_feature="bmi1" target_feature="bmi2" target_feature="fma" target_feature="pclmulqdq" target_feature="popcnt" target_feature="sse4.1" target_feature="sse4.2" target_

(Cut at 600 characters. There are many more.) So the same source gives very different machine code depending on the flag. Let me show you.

Build it (the "how")

Here are four tiny functions. The file is docs/av.rs.

// docs/av.rs
#[no_mangle] pub fn sum_f32(x: &[f32]) -> f32 { let mut s = 0.0; for &v in x { s += v; } s }
#[no_mangle] pub fn sum_u32(x: &[u32]) -> u32 { x.iter().fold(0u32, |a, &b| a.wrapping_add(b)) }
#[no_mangle] pub fn add_idx(a: &mut [f32], b: &[f32]) { for i in 0..a.len() { a[i] += b[i]; } }
#[no_mangle] pub fn add_zip(a: &mut [f32], b: &[f32]) { for (x, y) in a.iter_mut().zip(b) { *x += *y; } }

I compiled it twice, once with the default target and once with -C target-cpu=x86-64-v3 (this is a named level that means AVX2 + FMA + BMI2, roughly the CPUs since 2013), and then counted instructions with objdump:

rustc --edition 2021 -O --crate-type=cdylib av.rs -o libav_base.so
rustc --edition 2021 -O -C target-cpu=x86-64-v3 --crate-type=cdylib av.rs -o libav_v3.so
objdump -d --no-show-raw-insn -M intel libav_base.so | less     # find <sum_f32>:

The counts of scalar adds (addss), packed adds (addps, paddd, and their v forms), per function and build:

base sum_f32: scalar-add(addss/vaddss)=9 packed-add(addps/paddd/vaddps/vpaddd)=0 bounds-panic-calls=0
base sum_u32: scalar-add(addss/vaddss)=0 packed-add(addps/paddd/vaddps/vpaddd)=5 bounds-panic-calls=0
base add_idx: scalar-add(addss/vaddss)=1 packed-add(addps/paddd/vaddps/vpaddd)=2 bounds-panic-calls=0
v3 sum_f32: scalar-add(addss/vaddss)=9 packed-add(addps/paddd/vaddps/vpaddd)=0 bounds-panic-calls=0
v3 sum_u32: scalar-add(addss/vaddss)=0 packed-add(addps/paddd/vaddps/vpaddd)=13 bounds-panic-calls=0
v3 add_idx: scalar-add(addss/vaddss)=1 packed-add(addps/paddd/vaddps/vpaddd)=4 bounds-panic-calls=0

Three lessons from the table.

Lesson 1: the float sum never vectorizes. sum_f32 uses 9 scalar addss and 0 packed adds, in both builds. Not even x86-64-v3 helps. This is the reorder rule. Here is the start of its machine code in the default build, so you can see what "scalar" looks like:

    8600:	test   rsi,rsi
    8603:	je     8621 <sum_f32+0x21>
    8605:	lea    rcx,[rsi*4+0x0]
    860d:	add    rcx,0xfffffffffffffffc
    8611:	mov    eax,ecx
    8613:	not    eax
    8615:	test   al,0x1c
    8617:	jne    8625 <sum_f32+0x25>
    8619:	xorps  xmm0,xmm0
    861c:	mov    rax,rdi
    861f:	jmp    864d <sum_f32+0x4d>
    8621:	xorps  xmm0,xmm0

Look for xorps xmm0,xmm0 (clearing the single accumulator) and then a loop of addss: one float at a time.

Lesson 2: the integer sum vectorizes. sum_u32: 0 scalar adds, 5 packed adds (SSE2) in the default build, and 13 with v3 (AVX2 unrolls more). Wrapping integer addition can be reordered freely, so the compiler does.

Lesson 3: the same source, wider code with a flag. The v3 build of add_zip uses ymm registers (256-bit) in a 4× unrolled loop:

    85f0:	vmovups ymm0,YMMWORD PTR [rdx+rcx*4]
    85f5:	vmovups ymm1,YMMWORD PTR [rdx+rcx*4+0x20]
    85fb:	vmovups ymm2,YMMWORD PTR [rdx+rcx*4+0x40]
    8601:	vmovups ymm3,YMMWORD PTR [rdx+rcx*4+0x60]
    8607:	vaddps ymm0,ymm0,YMMWORD PTR [rdi+rcx*4]
    860c:	vaddps ymm1,ymm1,YMMWORD PTR [rdi+rcx*4+0x20]

The default build uses only xmm (128-bit). Both are "vectorized", but one is twice as wide.

What about bounds checks? The add_idx function, which indexes a[i] += b[i], also vectorized in both builds on this compiler (1 scalar add and 2 or 4 packed adds, the scalar one being the leftover loop). Rust's iterator style (zip) is the safest way to make this easy for the compiler, but for this simple pattern LLVM handled indexing too. I did not find a case where indexing blocked it in these four examples, and I would not trust that for more complex loops: always look at the output.

The danger of the easy fix. "Just build with -C target-cpu=native and get the wide code for free." Look at what happens to a normal-looking binary. I built the layout benchmark of Chapter 7 with and without the flag and counted wide-register instructions in the binary:

./target/release/e6_layout zmm=0 ymm=0
/tmp/native_t/release/e6_layout zmm=72 ymm=236

The default binary: 0 zmm and 0 ymm instructions. The native binary: 72 zmm and 236 ymm, in code where I wrote no SIMD at all. The compiler sprinkled AVX-512 through the whole program. That binary now crashes with "illegal instruction" on any CPU without AVX-512, and a runtime check cannot save it, because the compiler may have used wide instructions in startup code and ordinary functions. We return to this in Chapter 12.

The first instinct after seeing a slow float sum is to switch on fast-math. Rust deliberately has no global fast-math. The honest move is the one from Chapter 4: you pick the order, write 8 accumulators in plain Rust, and the compiler vectorizes it. You can see this in the Chapter 4 table: 8 accumulators (scalar src) ran at 47.9 GB/s, compared with 5.9 for the naive loop, without a single intrinsic.

What this gives us. A workflow:

1. Write the clear scalar version. It is your oracle.
2. Check the machine code. Is it vectorized? With which width?
3. If it is not, change the source until it is: iterators, chunks_exact, 8 accumulators for floats.
4. Only if that fails (shuffles, compress, movemask, saturating tricks), reach for intrinsics.

Why it is this way (the "why")

Why does the compiler refuse to vectorize a float sum? Why not just do it and hope?

Because the compiler's contract is "the optimized program gives the same result as the source program." Reordering float additions changes the result. If the compiler did it silently, one day a number in your bank statement would be off in the 7th digit and you would have no way to find out why. Languages like C let you opt in with -ffast-math, which tells the compiler "I don't care". It also lets the compiler assume that no value is NaN or infinite, which silently breaks code that checks for them. Rust's designers chose to make you say it explicitly in the source instead.

Why does the default target stop at SSE2? Because a Rust binary is a promise to run on a whole population of CPUs. The people who build an app cannot know whether the user's CPU has AVX-512. The standard way out is: build for the lowest common level, and ship kernels for faster levels inside the same binary, picking at runtime. That is Chapter 12.

The trade-off being made: portability versus peak speed from the compiler alone. You can have both, but only if you pay in explicit code (multiversioning).

What would have to be true for the compiler to be enough? If everyone shipped to a known fleet (your own servers), you could set a baseline with -C target-cpu=... for exactly that CPU, and the auto-vectorizer would use everything. Large cloud companies do this. A public library cannot.

Common mistakes

• Skipping the check. Writing intrinsics for a loop the compiler already vectorized.
• Believing "it's -O3, it must be vectorized." Check with objdump.
• Forgetting the default is 128-bit. Your "vectorized" code might use only a quarter of AVX-512.
• Using target-cpu=native for a binary that you ship.
• Counting xmm instructions as proof. Scalar floats also live in xmm registers (addss). Look for packed instructions (addps, paddd, vmovups).

Say it simply

The compiler auto-vectorizes loops that are independent, reorderable and safe from overlap, but by default only with 128-bit registers. Float sums are not reorderable, so you must pick the order yourself, for example with 8 accumulators. Always read the machine code before writing intrinsics, and never ship a native build to other people's machines.

Check yourself

Quick questions

1. (Warm-up, explain) Why does sum_u32 vectorize and sum_f32 not? Answer in three sentences.
2. (Core, explain) What is the difference between the default build and -C target-cpu=x86-64-v3 for the same source?
3. (Core, predict) Predict: would max over a [f32] slice vectorize in the default build? (Hint: is max order-dependent for floats? What about NaN?) Then check it.
4. (Warm-up, trace) In objdump output, which of these tells you a loop is really vectorized: addss, addps, vaddps ymm, vmovss?

Problems

5. (Core, implement) Take any for loop from your own code. Compile it as a cdylib with #[no_mangle] like av.rs. Count packed versus scalar instructions with objdump. Report.
6. (Core, change one thing) Rewrite sum_f32 with 16 accumulators in plain Rust. Does the compiler vectorize it? With which width in the default build? Use the benchmark from Chapter 4 to check the speed. Is 16 better than 8?
7. (Stretch, break it) Write a loop that looks vectorizable but the compiler cannot vectorize: add an early break, or a function call the compiler cannot see through, or an indexed access through a Vec<usize>. Which one breaks it? Use objdump as proof. How would you fix each?

Think harder

8. (Stretch, design and defend) Your company ships a Rust CLI tool to customers' laptops. Their CPUs are unknown. Choose a build strategy: (a) default, (b) x86-64-v3, (c) native, (d) default plus runtime-dispatched kernels. Defend your choice. What would change it if the tool ran only on your own servers?

Chapter 10. std::arch and the three contracts

In this chapter: you learn how to call real SIMD instructions from stable Rust, and the three separate promises (compile-time, function-level and runtime) that make it safe.

You need first: Chapters 1, 3 and 9.

The problem first

You checked the machine code (Chapter 9) and the compiler did not vectorize your loop. You need the exact instructions: pcmpeqb followed by pmovmskb for a byte search, say. Rust gives them to you through the module std::arch. The catch: if you call one of these instructions on a CPU that does not have it, the program dies with SIGILL (illegal instruction) and no helpful message. So the language wraps them in three layers of protection.

The model (the "what")

The module std::arch::x86_64 has one Rust function for almost every x86 SIMD instruction. These functions are called intrinsics. Their names look cryptic but follow a pattern:

_mm256_add_epi32
 |  |    |    |
 |  |    |    +-- element type: epi32 = packed 32-bit signed ints
 |  |    +------- operation: add
 |  +------------ register width: 256 bits  (nothing = 128, 512 = 512)
 +--------------- prefix: "multimedia" (a historical name; ignore it)

The register types are __m128/__m128d/__m128i (floats, doubles, ints), __m256…, and __m512…. A __m128i holds 128 bits of integers, and which lane width it means is decided by the instruction you apply, not by the type. That is the same lesson as Chapter 4: the instruction carries the meaning.

Now the three contracts.

Contract 1: compile time. #[cfg(target_arch = "x86_64")] and cfg(target_feature = "avx2"). These are decided when you build. target_arch says which family you compile for. If you compile for ARM, the std::arch::x86_64 module does not exist, so code that uses it must sit behind this cfg. cfg(target_feature = "avx2") is true only if you told the compiler (with -C target-feature or -C target-cpu) that every CPU you target has AVX2. It never asks the CPU. It is a statement about the build.

Contract 2: function level. #[target_feature(enable = "avx2")]. This attribute on one function tells the compiler: "inside this function, you may use AVX2 instructions". The rest of the program stays on the baseline. It is not a check. It does not ask the CPU anything. It only gives permission to emit those instructions.

Contract 3: runtime. is_x86_feature_detected!("avx2"). This macro asks the running CPU (through the cpuid instruction, and for AVX also checks that the operating system saves the wide registers when switching tasks). It is the only one of the three that looks at the real machine. It lives in std, not core. After the first call the answer is cached, so later calls cost under a nanosecond (we measure 0.48 ns in Chapter 12).

Put together: you write a function with #[target_feature], and you only call it after the runtime check said yes. Break that order and you get SIGILL.

How Rust makes this safe, and the 2025 change. Before Rust 1.86, a #[target_feature] function had to be unsafe fn: calling it was always unsafe. Since 1.86 you can write a safe fn with #[target_feature]:

• inside a function that has the feature, you can call other functions that need the same feature, safely, with no unsafe block;
• outside, in a function without the feature, calling it is still unsafe.

And since about Rust 1.87, intrinsics that take no pointers (adds, compares, shuffles, set1) are also safe to call inside a function that has the feature. Only intrinsics that take pointers (loads and stores) still need unsafe, because the compiler cannot check that the memory is valid.

Here is the compiler enforcing it. This is real output:

error[E0133]: call to function `add_avx2` with `#[target_feature]` is unsafe and requires unsafe function or block
 --> docs/err.rs:9:5
  |
9 |     add_avx2(a, b) // a plain fn calling it: the compiler stops us
  |     ^^^^^^^^^^^^^^ call to function with `#[target_feature]`
  |
  = help: in order for the call to be safe, the context requires the following additional target feature: avx2

error: aborting due to 1 previous error

For more information about this error, try `rustc --explain E0133`.

And the code that produced it (file docs/err.rs):

// docs/err.rs
use std::arch::x86_64::*;

#[target_feature(enable = "avx2")]
fn add_avx2(a: __m256i, b: __m256i) -> __m256i {
    _mm256_add_epi32(a, b) // no `unsafe` needed here: this fn has avx2
}

fn plain(a: __m256i, b: __m256i) -> __m256i {
    add_avx2(a, b) // a plain fn calling it: the compiler stops us
}

fn main() {}

One more real detail, because it contradicts a common assumption. Even a baseline feature such as SSE2 is not "free": I tried calling _mm_add_epi8 (an SSE2 intrinsic) from a plain function with no attribute and no unsafe. The compiler refused with the same E0133 error, "call to function ... with #[target_feature] is unsafe". The rule is about what the calling function declares, not about what the build happens to allow. In practice this means every function that uses intrinsics either has a matching #[target_feature] or wraps the call in unsafe.

Build it (the "how")

Let us read a real kernel built from these pieces. It finds the first occurrence of a byte, 32 bytes at a time, with AVX2:

// src/bytes.rs
#[target_feature(enable = "avx2")]
pub fn find_byte_avx2(h: &[u8], n: u8) -> Option<usize> {
    const L: usize = 32;
    let end = h.len() / L * L;
    let want = _mm256_set1_epi8(n as i8);
    let mut i = 0;
    while i < end {
        // SAFETY: i + 32 <= end <= len.
        let blk = unsafe { _mm256_loadu_si256(h.as_ptr().add(i).cast()) };
        let m = _mm256_movemask_epi8(_mm256_cmpeq_epi8(blk, want)) as u32;
        if m != 0 {
            return Some(i + m.trailing_zeros() as usize);
        }
        i += L;
    }
    h[end..].iter().position(|&b| b == n).map(|k| end + k)
}

Read it top to bottom.

1. #[target_feature(enable = "avx2")] and a safe fn. Inside we can call _mm256_set1_epi8, _mm256_cmpeq_epi8, _mm256_movemask_epi8 with no unsafe, because they take no pointers.
2. end = len / 32 * 32 is the last index we can process with a full vector. Integer division rounds down, so end is the largest multiple of 32 not above len.
3. _mm256_set1_epi8(n as i8) copies the needle into all 32 lanes. (n as i8 keeps the same 8 bits, and the compare does not care about sign.)
4. The load is unsafe, with a // SAFETY: comment that states why it is in bounds: i + 32 <= end <= len. Writing that comment forces you to check.
5. _mm256_cmpeq_epi8 compares 32 byte pairs at once and writes 0xFF into lanes that match, 0x00 into the others. (That is the mask from Chapter 3.)
6. _mm256_movemask_epi8 takes the top bit of each byte lane and packs the 32 bits into one i32. If it is not zero, at least one byte matched.
7. m.trailing_zeros() counts the zero bits at the bottom, which is the index of the lowest set bit, which is the lane of the first match. Because we scan blocks in increasing address order, the first block with any match holds the first match overall. So this gives the same answer as the scalar position.
8. The tail is a plain scalar loop on h[end..].

Now the caller. It must prove the CPU has AVX2:

if is_x86_feature_detected!("avx2") {
    // SAFETY: we just checked for AVX2.
    let r = unsafe { find_byte_avx2(hay, b'x') };
}

The first instinct for many people is to write #[cfg(target_feature = "avx2")] and think they have checked the CPU. They have not. That cfg is true only if the build enabled AVX2 for the whole program. On a default build it is false, so the kernel is silently left out. And if you build with -C target-feature=+avx2 then it is true for everyone, and the binary dies on non-AVX2 machines. The runtime macro is the only honest check.

A rule that avoids a real bug: detect every feature the function uses. AVX2 does not imply FMA in Rust's feature model. A kernel that calls _mm256_fmadd_ps has #[target_feature(enable = "avx2,fma")], and the check must be avx2 and fma. The AVX-512 kernels in this book enable avx512f,avx512bw and the check must test both.

What this gives us. A three-step shape we will reuse in every kernel:

1. A scalar function (the oracle).
2. One #[target_feature] function per instruction-set level.
3. A small safe wrapper that checks the CPU and calls the best one. (Chapter 12.)

Why it is this way (the "why")

Why three layers and not one? Because they answer three different questions:

Mixing them up is the root of most SIMD bugs. A kernel can be allowed (contract 2) and not run-checked (contract 3): SIGILL. A kernel can be run-checked but not allowed: the compiler refuses to emit the instruction and fails to compile or falls back to slow code.

Why was safe #[target_feature] added so late (1.86)? Because the old rule, "any call is unsafe", is too blunt: it forced unsafe onto code that was perfectly safe, which taught people to stop reading unsafe blocks. The refined rule says exactly what is dangerous: calling it from a place that does not have the feature. That is a narrow, checkable condition, so the compiler now checks it. This is the Rust design habit: make the unsafe part as small and as precise as you can. The result is that unsafe marks only "the CPU might not have this" and "this pointer might be invalid", and everything else is checked for you.

What would have to be true for this to be wrong? If CPUs could change features at runtime (hot-plug a core without AVX into a running process on a heterogeneous chip), a one-time check at startup would not be enough. Real examples are rare on x86, but mixed-core chips and cloud migrations are a reason some libraries re-check in long-lived processes, and a reason to read how Linux asks a process for permission before it uses AMX tile state (the arch_prctl call).

Common mistakes

• Assuming cfg(target_feature) checks the CPU. It does not.
• Missing a feature in the check. Using FMA with only an AVX2 check.
• Calling an unsafe kernel before the check, for example from a constructor or a static initializer.
• Using the wrong module (x86 vs x86_64) on 32-bit builds.
• Unsafe blocks with no explanation. Always write a // SAFETY: comment. If you cannot, the block is wrong.
• Putting #[inline(always)] on a #[target_feature] function. The compiler rejects this combination.

Say it simply

std::arch gives you one Rust function per SIMD instruction. Three layers keep it safe: a build-time cfg for which CPU family you compile for, a per-function #[target_feature] that lets the compiler emit the instructions, and a run-time is_x86_feature_detected! that asks the real CPU. You must run the check before calling the function, and the check must cover every feature the function uses.

Check yourself

Quick questions

1. (Warm-up, trace) Decode _mm512_cmpeq_epi8_mask: what width, what operation, what element type, and what do you think _mask at the end means?

2. (Core, explain) Explain in five sentences the difference between cfg(target_feature = "avx2"), #[target_feature(enable = "avx2")] and is_x86_feature_detected!("avx2"). Which one of them asks the CPU?

3. (Core, predict) Before running it, predict whether this compiles, and why. Then try:

   #[target_feature(enable = "avx2")]
   fn a(x: __m256i) -> __m256i { _mm256_add_epi32(x, x) }
   #[target_feature(enable = "avx2")]
   fn b(x: __m256i) -> __m256i { a(x) }
   fn c(x: __m256i) -> __m256i { a(x) }

4. (Core, explain) Why does _mm256_loadu_si256 require unsafe but _mm256_add_epi32 does not (inside a function with avx2)?

Problems

5. (Core, implement) Write count_byte_sse2 (count matches of a byte using _mm_cmpeq_epi8, _mm_movemask_epi8 and count_ones). Test against iter().filter().count() for all lengths 0 to 100.
6. (Core, break it and fix it) Write a kernel with _mm256_fmadd_ps but #[target_feature(enable = "avx2")] only. What does the compiler say? Fix it. Then write the runtime check that goes with it.
7. (Stretch, implement) Write a function fn best_level() -> &'static str returning "avx512", "avx2", "sse2" using is_x86_feature_detected!. For "avx512" decide which sub-features you require and say why in a comment.

Think harder

8. (Stretch, design and defend) A teammate proposes #![cfg(target_feature = "avx2")] at the top of your SIMD module "so we don't need runtime checks". When is this correct? When does it hurt? What would you do instead for a public crate?

Chapter 11. Byte kernels: find, count, sum, clamp

In this chapter: you build four byte kernels end to end, see them beat the scalar loop by 20× to 40×, compare one of them with the best real-world library (the memchr crate), and learn where and why the advantage disappears.

You need first: Chapters 3, 8 and 10.

The problem first

Bytes are everywhere: text, network packets, file formats, images. Four small operations come up constantly:

• find the first occurrence of a byte (every str::find, every parser looking for a newline),
• count occurrences (line counting, wc -l),
• sum all bytes (checksums, histograms),
• clamp every byte into a range (image processing).

A scalar loop handles one byte per step. A 256-bit register holds 32 bytes. That is the clearest win SIMD offers, and it is where you should start learning.

The model (the "what")

All four kernels follow one shape: load 32 bytes, do a lane-wise operation, and turn the result into something scalar.

 memory:  [ b0 b1 b2 ... b31 ] [ b32 ... b63 ] ...      (blocks of 32 bytes)
              |                      |
          loadu  (1 instr)        loadu
              |
        cmpeq with 32 copies of the needle   -> mask: 0xFF where equal
              |
        movemask                              -> 32-bit integer, bit i = byte i matched
              |
        if mask != 0: trailing_zeros()        -> index of first match

• find uses cmpeq, movemask and trailing_zeros. We just read it in Chapter 10.
• count uses the same cmpeq, but adds up the masks instead of looking at the first one.
• sum uses a single instruction, psadbw (sum of absolute differences), to add 8 bytes at a time into a 64-bit lane.
• clamp uses max then min, both lane-wise, with unsigned byte versions.

Two of them hide a trap that is worth studying.

Count: a narrow counter. The mask for a match is 0xFF, which as a signed byte is −1. So counter = counter - mask adds 1 to each matching lane. We keep 32 per-lane counters in a u8 vector, which is cheap, but a u8 counter can only hold 255. After 255 blocks it could overflow, and the lane would wrap silently. The fix is to run the inner loop for at most 255 blocks, then flush the narrow counters into a wide total, using the psadbw trick, and start again:

// src/bytes.rs
#[target_feature(enable = "avx2")]
pub fn count_byte_avx2(h: &[u8], n: u8) -> usize {
    let end = h.len() / 32 * 32;
    let want = _mm256_set1_epi8(n as i8);
    let mut total = 0usize;
    let mut i = 0;
    while i < end {
        // Each block can add at most 255 per byte lane before a u8 counter would wrap.
        // So we keep a u8 counter for up to 255 blocks, then flush to `total`.
        let mut cnt = _mm256_setzero_si256();
        let stop = (i + 255 * 32).min(end);
        while i < stop {
            // SAFETY: i + 32 <= end <= len.
            let blk = unsafe { _mm256_loadu_si256(h.as_ptr().add(i).cast()) };
            // A match is 0xFF == -1 as a byte. Subtracting -1 adds 1.
            cnt = _mm256_sub_epi8(cnt, _mm256_cmpeq_epi8(blk, want));
            i += 32;
        }
        let sad = _mm256_sad_epu8(cnt, _mm256_setzero_si256()); // 4 x u64 partial sums
        let mut lanes = [0u64; 4];
        // SAFETY: lanes is 32 writable bytes; storeu has no alignment need.
        unsafe { _mm256_storeu_si256(lanes.as_mut_ptr().cast(), sad) };
        total += lanes.iter().sum::<u64>() as usize;
    }
    total + h[end..].iter().filter(|&&b| b == n).count()
}

That is the general recipe for narrow accumulators: prove a bound, flush before it. In the code, 255 * 32 bytes is the largest span before any one lane can see 255 matches. The test suite runs an input of all-matching bytes, 255 * 32 * 3 + 5 long, specifically to hit this boundary.

Sum: use the instruction designed for it. psadbw (_mm256_sad_epu8) computes the absolute difference between two byte vectors and adds each group of 8 differences into one 64-bit lane. Against zero, the absolute difference of a byte is the byte itself, so it is a sum of 8 bytes per lane, into a lane wide enough to never overflow:

// src/bytes.rs
#[target_feature(enable = "avx2")]
pub fn sum_u8_avx2(v: &[u8]) -> u64 {
    let end = v.len() / 32 * 32;
    let zero = _mm256_setzero_si256();
    let mut acc = zero;
    let mut i = 0;
    while i < end {
        // SAFETY: i + 32 <= end <= len.
        let blk = unsafe { _mm256_loadu_si256(v.as_ptr().add(i).cast()) };
        // sad against zero = sum each group of 8 bytes into one u64 lane. No overflow.
        acc = _mm256_add_epi64(acc, _mm256_sad_epu8(blk, zero));
        i += 32;
    }
    let mut lanes = [0u64; 4];
    // SAFETY: lanes is 32 writable bytes.
    unsafe { _mm256_storeu_si256(lanes.as_mut_ptr().cast(), acc) };
    lanes.iter().sum::<u64>() + v[end..].iter().map(|&b| b as u64).sum::<u64>()
}

This is widening from Chapter 4: sum bytes into 64-bit lanes. No overflow handling needed.

Build it (the "how")

The first kernel, SSE2 (16 bytes), is a good warm-up. Its twin for AVX2 (32 bytes) you read in Chapter 10. Then comes the version that gains a lot with one idea.

The first instinct for speed is "wider registers". But look at the AVX2 find_byte: at every 32 bytes it does a compare, a movemask, a test, and a branch. The branch happens once per 32 bytes. We can cut it: test 128 bytes, then branch once. Compare four vectors, OR the four masks together, and only if the combined mask is non-zero do the careful work of finding which block and which byte:

// src/bytes.rs
/// Unrolled: test 128 bytes, then branch once. Fewer branches per byte.
#[target_feature(enable = "avx2")]
pub fn find_byte_avx2_unrolled(h: &[u8], n: u8) -> Option<usize> {
    let want = _mm256_set1_epi8(n as i8);
    let mut i = 0;
    while i + 128 <= h.len() {
        // SAFETY: i + 128 <= len.
        let (e0, e1, e2, e3) = unsafe {
            let p = h.as_ptr().add(i);
            (
                _mm256_cmpeq_epi8(_mm256_loadu_si256(p.cast()), want),
                _mm256_cmpeq_epi8(_mm256_loadu_si256(p.add(32).cast()), want),
                _mm256_cmpeq_epi8(_mm256_loadu_si256(p.add(64).cast()), want),
                _mm256_cmpeq_epi8(_mm256_loadu_si256(p.add(96).cast()), want),
            )
        };
        let any = _mm256_or_si256(_mm256_or_si256(e0, e1), _mm256_or_si256(e2, e3));
        if _mm256_movemask_epi8(any) != 0 {
            for (k, e) in [e0, e1, e2, e3].into_iter().enumerate() {
                let m = _mm256_movemask_epi8(e) as u32;
                if m != 0 {
                    return Some(i + k * 32 + m.trailing_zeros() as usize);
                }
            }
        }
        i += 128;
    }
    // Finish the last < 128 bytes with the plain 32-byte version.
    find_byte_avx2(&h[i..], n).map(|k| i + k)
}

Now the real numbers. Needle absent (worst case: scan everything). GB/s, higher is better:

     size |    scalar    memchr      sse2      avx2    avx2x4    avx512   (GB/s)
     16 B |       2.2       3.1       8.7       2.0       1.2       5.6
     64 B |       2.7       8.0      13.8      17.2      21.5      20.1
    256 B |       2.4      35.7      20.6      31.1      50.3      37.3
    4 KiB |       1.6     100.6      26.2      51.5      85.3     110.2
   64 KiB |       2.8      81.0      26.3      35.1      49.4      52.3
    1 MiB |       2.8      75.1      30.3      28.9      41.5      48.5
   16 MiB |       2.9      28.9      25.2      23.4      27.5      27.3

Columns: scalar is iter().position(). memchr is the memchr crate (the one ripgrep uses), which does its own runtime dispatch. sse2, avx2, avx2x4 (unrolled) and avx512 are the kernels of this book.

Read it like a landscape again.

• Scalar: 1.6 to 2.9 GB/s at every size. A byte per cycle or so.
• At 4 KiB (all in L1): sse2 26, avx2 52, avx2x4 85, avx512 110, memchr 101. The wider and the less branchy, the faster. The unrolled AVX2 gets 1.65× over the plain AVX2, just from fewer branches.
• Our AVX-512 version matches the library at 4 KiB (110 vs 101 GB/s) with 20 lines of code.
• At 64 KiB to 1 MiB (L2): everything drops. The library is still ahead (81 and 75 GB/s) of our best (52 and 48): it has more tricks (it reads in a smarter order and aligns). That is what a serious library adds.
• At 16 MiB (main memory): memchr 29, sse2 25, avx2 23, avx2x4 27.5, avx512 27. Everything converges to about 25 to 29 GB/s. Memory is the limit now, and all the lane width in the world cannot help. This is Chapter 5, live.
• At 16 bytes: scalar 2.2, memchr 3.1, sse2 8.7, avx2 2.0, avx512 5.6. Here the AVX2 kernel is the worst, because 16 bytes is less than one 32-byte vector: it does only its scalar tail plus the setup. Small inputs punish wide kernels. This is the "fixed cost" problem, and it drives Chapter 12.

Even for this "simple" problem, a professional library still beats a clean hand-written kernel at some sizes. That is useful humility: your kernel is a solid baseline, and a mature library is the benchmark to beat before you ship your own.

Finally, clamp and the tail. We built the AVX2 clamp and the overlap version in Chapter 8:

// src/bytes.rs
#[target_feature(enable = "avx2")]
pub fn clamp_u8_avx2(v: &mut [u8], lo: u8, hi: u8) {
    let end = v.len() / 32 * 32;
    // `as i8` keeps the same 8 bits. max_epu8/min_epu8 read them as UNSIGNED.
    let (lo_v, hi_v) = (_mm256_set1_epi8(lo as i8), _mm256_set1_epi8(hi as i8));
    let mut i = 0;
    while i < end {
        // SAFETY: i + 32 <= end <= len; we read then write the same 32 bytes.
        unsafe {
            let p = v.as_mut_ptr().add(i).cast::<__m256i>();
            let b = _mm256_loadu_si256(p);
            _mm256_storeu_si256(p, _mm256_min_epu8(_mm256_max_epu8(b, lo_v), hi_v));
        }
        i += 32;
    }
    for b in &mut v[end..] {
        *b = (*b).clamp(lo, hi);
    }
}

It needs no sign trick because SSE2/AVX2 have unsigned byte min and max (_mm256_min_epu8). The as i8 on the broadcast keeps the same 8 bits, and the unsigned min/max read them back as unsigned.

Why it is this way (the "why")

Why does movemask exist at all? Why not branch on the vector directly?

Because the CPU's branch instructions look at one integer flag, not at 32 byte lanes. movemask is the bridge between the vector world and the scalar world: it takes one bit from each lane and makes a normal integer. Without it, SIMD could only produce vector results, never answer yes/no questions. With it, a whole family of algorithms (search, parse, scan) becomes possible: use SIMD to compute a mask, then use fast scalar bit tricks (trailing_zeros, count_ones, &) on the mask. Remember this pattern. It is the backbone of the fastest text-processing code in the world, including Chapters 18 and 19.

Why does the unrolled version win? Because it changes a throughput-bound-by-branches loop into a throughput-bound-by-loads loop. A branch per 32 bytes limits you to however many branches the front end can predict and issue. Four loads and three ORs per branch lets the load units, which can do two or three loads per cycle, run closer to full speed. It is the same idea as Chapter 2, applied to branches: amortize the fixed cost of a decision over more work.

The trade-off: the unrolled version is worse for short inputs. A 100-byte input has no full 128-byte block at all and falls through to the slower path. Everything in optimization is a trade between the average case and the extremes.

What would have to be true for this design to be wrong? If matches are very common (for example you search for a space in English text, which appears every ~6 bytes), then almost every block has a hit, and the clever "test 128 bytes, branch once" logic does extra work for no benefit. In that regime a different design wins (process the hits, not the misses). The benchmark above is the needle absent case, which is the best case for this design. Always benchmark both ends.

Common mistakes

• trailing_zeros on a zero mask. It returns the bit width (32), which looks like a valid index. Always test mask != 0 first.
• A narrow accumulator that is never flushed. The u8 counter above.
• Using i8 compare for bytes above 127. Chapter 4.
• Benchmarking only the needle-absent case. The needle-first and needle-in-the-middle cases behave differently.
• Beating a library on one size and calling it a win. Check all sizes.

Say it simply

A byte kernel loads 32 bytes, compares or computes lane by lane, and uses movemask to turn the result into an ordinary integer, which trailing_zeros or count_ones then turns into an answer. It is easily 20 to 40 times faster than a byte loop on data in cache, and it converges to memory speed on big data. Mature libraries still win at some sizes. Tiny inputs favour smaller kernels.

Check yourself

Quick questions

1. (Warm-up, predict) For a 64-byte haystack with the needle only at index 40, what would find_byte_avx2 return? How many full 32-byte blocks does it inspect? What is the first block's movemask?
2. (Core, trace) A u8 counter lane sees a match in 300 consecutive blocks. What is its value if never flushed? Why does the code flush at 255 blocks?
3. (Core, explain) In five sentences: why does _mm256_sad_epu8 against zero compute a sum of bytes? What is "absolute difference" doing?
4. (Core, explain) Why does the unrolled kernel lose to the plain one on inputs under 128 bytes?

Problems

5. (Core, implement) Write find_last_byte_avx2 (the last occurrence). Which bit-scan do you use instead of trailing_zeros? In which order do you scan blocks? Test it against iter().rposition().
6. (Core, implement) Write is_ascii_avx2: return true if every byte is below 128. (Hint: movemask reads the top bit of each byte directly. Is there any compare needed at all?) Compare against [u8]::is_ascii from the standard library, which is also optimized.
7. (Stretch, break it and fix it) Change count_byte_avx2 so the inner loop runs for 256 * 32 bytes instead of 255 * 32. Build an all-match input of that length. What is the count it returns, and which lane overflowed? Fix it without changing the constant (hint: change the accumulator).
8. (Stretch, predict then check) Run e4_find and then change the haystack from a..z repeating to every byte equals the needle (a hit in the first block). Predict the GB/s of avx2x4 versus plain avx2. Then measure. Is the gap the same as in the absent case? Explain what you see using the number of branches each kernel runs per byte.

Think harder

9. (Stretch, design and defend) You will write a line counter that counts \n in a 10 GB file. Compare: (a) your AVX2 count_byte on mmaped memory, (b) memchr::memchr_iter, (c) the scalar loop. Which do you pick, and what would you measure to be sure?

Chapter 12. Dispatch: picking the kernel at runtime

In this chapter: you learn how one binary can carry several kernels and choose the right one for the CPU it lands on, what that choice costs, and how to avoid the cost for small inputs.

You need first: Chapters 9, 10 and 11.

The problem first

You wrote three kernels: scalar, AVX2 and AVX-512. A customer runs your program on a laptop with AVX2 and no AVX-512. A second customer runs it on an old server with neither. You have one binary. Who picks the kernel, and when?

If you pick at build time (with -C target-cpu), you must build one binary per CPU level. If you pick at runtime, you ship one binary that adapts. The second is what most libraries do. But runtime dispatch has costs, and small inputs expose them.

The model (the "what")

Dispatch means choosing which implementation to run. There are three styles.

Here is the cached function pointer style, the one this book uses. First the chooser, which runs once:

// src/bin/e10_dispatch.rs
fn choose() -> F {
    if std::is_x86_feature_detected!("avx2") { avx2_thunk } else { bytes::find_byte_scalar }
}

The wrapper that every caller sees. It is a normal safe function with a normal signature:

// src/bin/e10_dispatch.rs
fn dispatched(h: &[u8], n: u8) -> Option<usize> {
    CHOSEN.get_or_init(choose)(h, n)
}

Look at the thunk, the small safe function that is the only place the unsafe call lives:

// src/bin/e10_dispatch.rs
fn avx2_thunk(h: &[u8], n: u8) -> Option<usize> {
    // SAFETY: only returned by `choose` after avx2 was detected.
    unsafe { x86::find_byte_avx2(h, n) }
}

The // SAFETY: comment states the whole proof: this function is only ever returned by choose after AVX2 was detected. Every other function in your program stays safe. This is the pattern: one narrow unsafe seam, with a one-sentence proof.

A small-input cutoff. For tiny inputs a vector kernel is no help (Chapter 11), so we skip dispatch altogether:

// src/bin/e10_dispatch.rs
fn dispatched_with_cutoff(h: &[u8], n: u8) -> Option<usize> {
    if h.len() < 32 { return bytes::find_byte_scalar(h, n); }
    CHOSEN.get_or_init(choose)(h, n)
}

Build it (the "how")

What does dispatch actually cost? Measure it. First, one call to the detection macro:

one is_x86_feature_detected!("avx2") call: 0.48 ns

 bytes |   scalar     avx2     fn-ptr ptr+cutoff detect-each   (ns per call, needle absent)
     1 |     1.61     1.89       3.37       1.59       1.81
     4 |     3.22     3.46       4.10       3.18       3.33
     8 |     4.36     4.79       5.30       4.38       4.56
    16 |     6.87     7.56      10.53       6.91       9.70
    31 |    11.86    12.47      24.23      19.60      22.83
    32 |    21.63     2.66       3.32       2.99       2.36
    64 |    36.41     2.57       3.18       3.50       3.05
   128 |    74.12     3.59       4.06       3.91       3.30
   512 |   160.54     9.63       8.78       8.75       8.82

One is_x86_feature_detected!("avx2") takes 0.48 ns. That is roughly one or two cycles: it reads an atomic that was filled in on the first call. So even "detect every call" is cheap.

Now the table. Columns are nanoseconds per call, needle absent: scalar; avx2 (called directly, no dispatch); fn-ptr (cached pointer); ptr+cutoff (skip below 32 bytes); detect-each (call the macro every time).

Read it row by row.

• 1 to 16 bytes: scalar 1.6 to 6.9 ns, direct AVX2 1.9 to 7.6 ns. The AVX2 kernel does nothing useful below 32 bytes, because everything falls into its scalar tail. It is slightly worse (setup), and fn-ptr adds another ~0.5 to 3.5 ns. ptr+cutoff removes the whole cost: it matches the scalar column.
• 32 bytes: scalar 21.6 ns, AVX2 2.7 ns. Eight times faster, and the dispatch overhead (fn-ptr 3.3, detect-each 2.4) is small by comparison.
• 64 to 512 bytes: the SIMD kernels are 14× to 21× faster than scalar. The dispatch cost, about 0.5 to 1 ns, is invisible.

One row looks odd: at 31 bytes, the function-pointer version took 24.2 ns while the direct AVX2 call took 12.5. Both do the same scalar tail. I do not know why the pointer call came out twice as slow there, and I could not explain it by reading the assembly in the time I had. I would re-run it and look at the generated code before trusting it. I am leaving the number in because it is real output, and a benchmark you can't fully explain is still worth showing.

The conclusion the table supports: dispatch cost does not matter, but the cutoff does. A cached pointer costs about 1 ns. A scalar loop over 16 bytes costs about 7 ns. So the choice is not "dispatch or not" but "use the vector path only when there is enough work to pay for it", and that threshold is the vector width.

The first instinct to avoid the indirect call is to detect features with a static mut flag and a plain if. It works, but static mut is unsound under threads. OnceLock is the safe, standard way (stable since Rust 1.70), and the get_or_init fast path is a single atomic load.

What to do with several kernels. If you have SSE2, AVX2 and AVX-512 versions, order the checks strongest first:

fn choose() -> F {
    if is_x86_feature_detected!("avx512f") && is_x86_feature_detected!("avx512bw") { return avx512_thunk; }
    if is_x86_feature_detected!("avx2") { return avx2_thunk; }
    scalar
}

Note that the AVX-512 check lists both features its kernel uses.

The build-flag danger. All of this assumes the rest of your program was compiled for the baseline. In Chapter 9 you saw a target-cpu=native binary get 72 zmm instructions in code where no one wrote SIMD. If you ship that, no runtime check can save you: the wide instructions are in ordinary functions that run before and after your dispatch. Keep the baseline conservative. Put the wide instructions only in #[target_feature] functions that are reached after the check. If your deployment has a documented CPU floor (your own fleet), set that floor explicitly, record it next to the build configuration, and test on the oldest machine.

Libraries that automate this. The crates multiversion (generates several copies of a function and picks one) and pulp (vectorized closures with runtime dispatch) exist for this job, and wide gives portable SIMD types on stable. I did not benchmark or test them for this book, so I will not describe their performance. Chapter 14 shows how to evaluate a dependency like these.

In no_std code (kernels, firmware), std::is_x86_feature_detected! and OnceLock are not available. You can select with cfg(target_feature = "avx2") when the deployment CPU is fixed and documented, or you can pass a capability value in from the platform at start-up. There is no universal detector, because each platform exposes CPU features differently.

Why it is this way (the "why")

Why must dispatch be a function-level idea instead of something the compiler does for us?

Because CPU feature levels change what is legal to emit, and the compiler needs to know the boundary. Inside a function with AVX2 enabled, LLVM may use AVX2 anywhere, including in its own spills and moves. So the compiler wants whole functions as the unit, with a clear entry where the check has already been done. This also explains why a #[target_feature] function usually cannot be inlined into a caller that lacks the feature: inlining would put AVX2 instructions into a place the check does not cover. A call boundary is the price of that safety.

Another design choice: Rust refuses to let #[inline(always)] coexist with #[target_feature], because "always inline" cannot override that boundary. You have to pay one call per kernel invocation. For a loop that runs thousands of iterations, that is nothing. For a kernel that processes 16 bytes, it matters, which brings us back to the cutoff.

The trade-off being made. Binary size and complexity (several kernels per function) in exchange for one binary that is fast everywhere. The alternative (one binary per CPU level) trades build and release complexity for smaller binaries and zero dispatch cost. Large companies often pick that: they control the fleet. Open-source libraries pick runtime dispatch, because they do not.

What would have to be true for dispatch to be wrong? If every call processed only a few bytes, the indirect call would dominate and you would want the kernel inlined, which means a compile-time choice. Hash functions on short keys are an example, and they are why serious hash libraries pick kernels at build time.

Common mistakes

• Dispatching tiny inputs. Use a cutoff at the vector width.
• Forgetting one feature in the check. AVX-512 kernels that need avx512bw need that checked too.
• target-cpu=native in release builds that others run.
• A static mut flag instead of OnceLock.
• Re-detecting inside the hot loop. Detect once per call at most, never per block.
• Never testing the fallback. The scalar path must be tested too, on a CPU (or a forced code path) that uses it.

Say it simply

A binary can hold several kernels and pick one at runtime with is_x86_feature_detected!. The cheap, safe way is a cached function pointer behind a safe wrapper, with one unsafe call in a tiny thunk whose comment states why it is safe. The dispatch itself costs about a nanosecond, so the real question is the cutoff: below one vector's worth of data the scalar path wins.

Check yourself

Quick questions

1. (Warm-up, predict) A find_byte call on 20 bytes costs about 7 ns scalar. The pointer call adds about 1 ns. If you call it 100 million times, how much time do you lose to dispatch? Is it worth removing?
2. (Core, explain) Why does a #[target_feature] function usually stay a real call instead of being inlined?
3. (Core, trace) Write the order of checks for a kernel family with AVX-512, AVX2 and SSE2 levels. Which features does the AVX-512 level need?
4. (Core, explain) In five sentences: why can target-cpu=native break a program even if you also use runtime dispatch?

Problems

5. (Core, implement) Add a third kernel, find_byte_avx512, to the chooser from this chapter. Make choose return it only when both avx512f and avx512bw are present. Test that it is reached, by adding a fn level_name() you can assert on.
6. (Core, break it and fix it) Make the dispatcher testable: add a function fn find_byte_with(level: Level, ...) that forces a given kernel, and test every level against the scalar oracle on every length from 0 to 300. Why does this matter if your machine only ever chooses AVX-512?
7. (Stretch, predict then check) Run e10_dispatch. Find the crossover size where the fn-ptr version first beats scalar. Predict first. How would the crossover move if the indirect call cost twice as much?

Think harder

8. (Stretch, design and defend) You maintain a public Rust crate with a SIMD byte-search function. Users call it on inputs from 4 bytes to 4 GB. Decide: (a) cutoff policy, (b) which levels, (c) how to test levels you cannot run in CI, (d) what to document. Defend each choice.

Chapter 13. AVX-512 in Rust

In this chapter: you learn what is genuinely new in AVX-512 (mask registers, masked loads, compress), build tails with no scalar code at all, prove with a guard page that masked loads really do not touch masked-off bytes, and learn the honest story of frequency and portability.

You need first: Chapters 3, 8, 10 and 12. A CPU with AVX-512 to run it (many Intel Xeons, and AMD Zen 4 and newer).

The problem first

You have used 256-bit AVX2. Now you meet a CPU with 512-bit registers. The obvious guess is "twice as wide, so about twice as fast". There is some of that. But the bigger change in AVX-512 is not the width. It is a new way to say "only these lanes".

Remember the tail problem (Chapter 8): the last few elements do not fill a vector, and you either wrote a scalar loop or played tricks with overlaps. AVX-512 gives you a better answer in the hardware.

The model (the "what")

AVX-512 is not one extension. It is a family: a base called avx512f (foundation) and many add-ons. You must check each one you use.

The machine for this book reports all of them (see the table in the front matter). Others, especially older Intel desktop chips and some AMD chips, have only a subset. Rust names each one, and stable Rust supports the AVX-512 target features and intrinsics from 1.89 (August 2025). Before that you needed nightly.

What is new, in three parts.

1. Mask registers. AVX-512 has eight 64-bit registers k0 to k7. A mask register holds one bit per lane. A byte compare on 64 bytes gives a 64-bit mask directly in a k register: no movemask needed. In Rust the type is __mmask64 (and __mmask16 for 16 lanes), and it is just an integer, so trailing_zeros() and count_ones() work on it as is.

2. Masked operations. Almost every instruction can take a mask: "do this only in lanes where the mask bit is 1, and in the other lanes either keep the old value (merge) or write zero (zeroing)." The intrinsic names carry it: _mm512_mask_add_epi32 merges, _mm512_maskz_add_epi32 zeroes.

3. Masked loads that do not fault. This one matters most. _mm512_maskz_loadu_epi8(k, ptr) loads only the bytes whose mask bit is 1. The hardware guarantees that masked-off bytes are never accessed, so they cannot cause a page fault even if they lie in unmapped memory. That means the tail of any loop can be one masked load, with no scalar code.

Here is the AVX-512 byte search from Chapter 11, the whole function, so you can read the tail:

// src/bytes.rs
#[target_feature(enable = "avx512f,avx512bw")]
pub fn find_byte_avx512(h: &[u8], n: u8) -> Option<usize> {
    const L: usize = 64;
    let want = _mm512_set1_epi8(n as i8);
    let mut i = 0;
    while i + L <= h.len() {
        // SAFETY: i + 64 <= len.
        let blk = unsafe { _mm512_loadu_si512(h.as_ptr().add(i).cast()) };
        let m: u64 = _mm512_cmpeq_epi8_mask(blk, want);
        if m != 0 {
            return Some(i + m.trailing_zeros() as usize);
        }
        i += L;
    }
    // The tail, with NO scalar loop: a masked load only touches the bytes
    // whose mask bit is 1. Masked-off bytes are never read, so they cannot fault.
    let rem = h.len() - i;
    if rem > 0 {
        let k: u64 = (1u64 << rem) - 1; // rem is 1..=63 here
        // SAFETY: only the low `rem` bytes are read, and they are in bounds.
        let blk = unsafe { _mm512_maskz_loadu_epi8(k, h.as_ptr().add(i).cast()) };
        // Zeroed lanes could match needle 0, so AND the result with k.
        let m = _mm512_cmpeq_epi8_mask(blk, want) & k;
        if m != 0 {
            return Some(i + m.trailing_zeros() as usize);
        }
    }
    None
}

Compare it with the AVX2 version. The main loop is the same shape with _mm512_cmpeq_epi8_mask returning a mask directly. The tail is five lines: build a mask k with rem low bits set (1 << rem) - 1, do a masked load, compare, and AND the result with k. (The AND is needed because masked-off lanes become zero, and if you are searching for the byte 0 they would falsely match.)

Another new instruction: compress. _mm512_maskz_compress_epi32(k, v) takes the lanes selected by k and packs them contiguously at the front. That is the filter primitive that AVX2 lacks, and we use it in Chapter 17.

Build it (the "how")

Claim to test: masked loads never touch masked-off bytes. This is a claim about hardware, and the cleanest test is to try to break it. Put the data at the very end of a memory page and make the next page unreadable, so any touch past the end crashes:

// src/bin/e11_guard.rs
//! Experiment 11: do masked loads really never touch bytes whose mask bit is 0?
//! We put the data at the very end of a page and make the next page unreadable.
use simdbook::bytes::x86;

fn main() {
    let page = 4096usize;
    // SAFETY: plain mmap/mprotect; we own the mapping.
    unsafe {
        let base = libc::mmap(std::ptr::null_mut(), 2 * page, libc::PROT_READ | libc::PROT_WRITE, libc::MAP_PRIVATE | libc::MAP_ANONYMOUS, -1, 0) as *mut u8;
        assert!(base as isize != -1);
        assert_eq!(libc::mprotect(base.add(page).cast(), page, libc::PROT_NONE), 0); // guard page
        for len in 1..=63usize {
            let start = base.add(page - len); // slice ends exactly at the guard page
            for i in 0..len { *start.add(i) = b'a'; }
            let h = std::slice::from_raw_parts(start, len);
            assert_eq!(x86::find_byte_avx512(h, b'#'), None);
        }
        println!("avx512 masked tail: lengths 1..=63 placed against a PROT_NONE page: no fault, correct answers");
        if std::env::args().any(|a| a == "--crash") {
            println!("now doing a plain 64-byte load that crosses into the guard page (this is undefined behaviour in Rust; we do it only to see the hardware react)...");
            let start = base.add(page - 10);
            let v = core::arch::x86_64::_mm512_loadu_si512(start.cast());
            println!("survived?! {:?}", core::hint::black_box(v));
        }
    }
}

The idea, drawn:

   page 1 (readable)                      page 2 (PROT_NONE: any access faults)
 +--------------------------------------+--------------------------------------+
 |  ........ [ the slice: len bytes ]   |  guard                               |
 +--------------------------------------+--------------------------------------+
                                       ^ slice ends exactly at the page boundary

For every length from 1 to 63 we place the slice so that its last byte is the last byte of the page, then call the AVX-512 search. If the masked load touched even one byte past the end, the program would die with SIGSEGV. Result:

avx512 masked tail: lengths 1..=63 placed against a PROT_NONE page: no fault, correct answers

No fault, and correct answers for all 63 lengths. The claim holds on this CPU.

Now the control experiment. With --crash the program also does a plain 64-byte load that crosses into the guard page, which is undefined behaviour in Rust. We do it only to see the hardware react. The program ends with a segmentation fault, and the shell reports exit code 139, which is 128 + 11, and 11 is the number of SIGSEGV:

avx512 masked tail: lengths 1..=63 placed against a PROT_NONE page: no fault, correct answers
now doing a plain 64-byte load that crosses into the guard page (this is undefined behaviour in Rust; we do it only to see the hardware react)...
Segmentation fault
exit code: 139

So the contrast is complete: a masked load of 10 valid bytes is safe at the end of a page, and an unmasked 64-byte load is not.

(A note on Rust's rules: a masked load in the standard library's AVX-512 intrinsics is documented to read only the unmasked elements, so writing it this way is within the language's rules. A plain 64-byte load past the end of a slice is not, even if the page happened to be mapped.)

The first instinct is to use AVX-512 everywhere because it is "the widest". Two cautions come from the data in this book.

First, the width does not always win. In Chapter 17 we filter with AVX-512's compress and with an AVX2 table-lookup trick; the AVX2 version is as fast or faster (0.13 to 0.27 ns per element for AVX2, 0.18 to 0.28 for AVX-512). Work that is bounded by loads or branches does not scale with the register width.

Second, the frequency story. On some Intel chips (notably the first generation of server chips with AVX-512), heavy 512-bit instructions made the core run at a lower clock to stay within power limits, which could slow down the rest of the program too. That is the reason you will read warnings about "AVX-512 downclocking". The honest picture today is: it depends on the chip generation, and newer chips are much less affected. Do not state that "AVX-512 always downclocks" or that "it never matters". On the machine used for this book, 512-bit FMA ran at about 1.8 instructions per cycle, near the 2 per cycle of 256-bit (Chapter 2), and in one of three runs the 512-bit rows were slower, which I could not attribute to anything but noise on a shared host. Measure whole-program time on your target hardware, not just the kernel.

Third, it is not available everywhere. A huge fraction of laptops and desktops has no AVX-512. AMD introduced it in Zen 4, executing 512-bit operations on 256-bit-wide hardware paths, so the gain on those chips comes from the new instructions and masks more than from raw width. Always keep an AVX2 version.

What this gives us. The AVX-512 byte search needs no scalar tail and no overlap trick and it is shorter than the AVX2 version once you count the cleanup code. That alone is a reason to like it.

Why it is this way (the "why")

Why did Intel add a whole new set of mask registers instead of reusing vector registers as in AVX2?

Because masks are used a lot, and used in two different ways. As data, they are the results of comparisons. As control, they say which lanes an instruction may touch. A separate register file lets the hardware read the mask early to decide which memory accesses to perform (which is why masked-off bytes cannot fault) without occupying a whole vector register per mask. It also lets one mask guard many instructions cheaply.

The design has a clear lineage and a clear successor. The idea of predicated vector execution comes from older vector supercomputers. Arm's SVE and the RISC-V Vector extension push it further: the mask (the predicate) controls the loop tail and the vector length is not fixed, so the same program runs on chips with different widths. AVX-512 is the x86 version of the same idea, tied to a fixed width of 512 bits.

The trade-off. More registers, more instruction variants, more silicon, and a fragmented feature list. The fragmentation is the cost: "supports AVX-512" is not one thing, and software has to check each sub-feature. This is also a lesson for any instruction-set design: every optional piece you add multiplies the testing matrix for everyone downstream.

What would have to be true for AVX-512 to be a mistake? If the power cost of wide units pushed clocks down so much that the average program got slower, and few programs could use masks well. History is mixed: it was a real concern on early chips, and it has largely faded as designs improved, but the fragmentation and the weak adoption on consumer CPUs are why many projects still treat it as optional.

Common mistakes

• Checking avx512f only. Byte operations need avx512bw. VNNI needs avx512vnni.
• Forgetting the AND with the tail mask after a zeroing masked load, when the needle could be zero.
• Using the 512-bit path for tiny inputs. The setup cost is the same as AVX2's, and it matters more for small data.
• Assuming wide equals fast without measuring the whole program.
• No fallback. The AVX2 path must exist and be tested.

Say it simply

AVX-512 adds 512-bit registers, but its biggest idea is a set of mask registers with one bit per lane. Masks make a compare result a plain integer, let almost any instruction work on selected lanes only, and let a masked load skip invalid bytes so the tail of a loop needs no scalar code. It is a family of features, not one, so check each one you use, and keep an AVX2 fallback.

Check yourself

Quick questions

1. (Warm-up, trace) A slice has 100 bytes. With 64-byte vectors, how many full vectors? What is the tail mask k for the leftover bytes, written in binary? (Hint: (1 << 36) - 1.)
2. (Core, explain) Why must the tail code in find_byte_avx512 AND the compare result with k?
3. (Core, explain) What is the difference between merge-masking and zero-masking? Which intrinsic prefix (mask_ or maskz_) is which?
4. (Core, explain) In four sentences: why does the guard-page experiment prove more than a test that only checks correct results?

Problems

5. (Core, implement) Write count_byte_avx512 using _mm512_cmpeq_epi8_mask and count_ones(). It needs no flush logic. Why? Test it against the scalar count for all lengths 0 to 300.
6. (Core, break it) Change the tail load in find_byte_avx512 to an unmasked _mm512_loadu_si512, and run it through the guard-page program. What happens? Which lengths crash?
7. (Stretch, predict then check) Run e4_find. At which size does AVX-512 stop beating the unrolled AVX2 kernel? Why do you think so? (Use Chapter 5.)

Think harder

8. (Stretch, design and defend) You are adding an AVX-512 path to a library. Half of your users have laptops without it. Argue for or against shipping it. What would you need to measure first, and what would you document?

Chapter 14. Other worlds: NEON, WebAssembly and std::simd

In this chapter: you learn how the same ideas look on Arm (NEON), in WebAssembly, and in Rust's portable SIMD library, and you learn how to choose between the layers.

You need first: Chapters 3, 10 and 12.

Honesty note. This sandbox has no ARM target and no WebAssembly target installed, so the NEON and WebAssembly code in this chapter was not compiled and not run. It is written from the architecture documentation, and it is marked Not run: check this on your machine. The std::simd code was run.

The problem first

Everything so far was x86. But your code may also need to run on a phone, an Apple laptop, a cloud server on Arm, or in a browser. Do you write each kernel four times? Sometimes yes. But you have options, and the choice is a real design decision.

The model (the "what")

Arm NEON (also called Advanced SIMD) is the SIMD of 64-bit Arm. It uses 128-bit registers, and on 64-bit Arm it is part of the baseline: every chip has it. The types are uint8x16_t (16 bytes), float32x4_t (4 floats), and so on. The Rust module is std::arch::aarch64. The names differ from x86: vaddq_f32 adds four floats (the q means "quad": 128-bit), vld1q_u8 loads 16 bytes, vceqq_u8 compares bytes.

The most important difference for us: NEON has no movemask. To turn a compare result into an integer you need a trick. The common one is a narrowing shift: shift each 16-bit pair right by 4 and keep the low 8 bits, which squashes 16 bytes of 0xFF/0x00 into 64 bits, 4 bits per original byte. The index of the first match is then trailing_zeros() / 4.

// Not run: check this on your machine. Written from the Arm intrinsics documentation.
#[cfg(target_arch = "aarch64")]
pub fn find_byte_neon(h: &[u8], n: u8) -> Option<usize> {
    use core::arch::aarch64::*;
    let end = h.len() / 16 * 16;
    let mut i = 0;
    // SAFETY: NEON is part of the aarch64 baseline; i + 16 <= end <= len for each load.
    unsafe {
        let want = vdupq_n_u8(n);
        while i < end {
            let blk = vld1q_u8(h.as_ptr().add(i));
            let eq = vceqq_u8(blk, want);                         // 0xFF where equal
            let packed = vshrn_n_u16::<4>(vreinterpretq_u16_u8(eq)); // 16 bytes -> 8 bytes (4 bits per input byte)
            let m = vget_lane_u64::<0>(vreinterpret_u64_u8(packed));
            if m != 0 {
                return Some(i + (m.trailing_zeros() / 4) as usize);
            }
            i += 16;
        }
    }
    h[end..].iter().position(|&b| b == n).map(|k| end + k)
}

WebAssembly simd128. WebAssembly is a portable instruction format that runs in browsers and server runtimes. Its SIMD extension is fixed at 128 bits, and the runtime translates each operation into whatever the host CPU has. Unlike NEON, Wasm does have a bitmask instruction.

// Not run: check this on your machine. Build with: RUSTFLAGS="-C target-feature=+simd128" cargo build --target wasm32-unknown-unknown --release
#[cfg(all(target_arch = "wasm32", target_feature = "simd128"))]
pub fn find_byte_wasm(h: &[u8], n: u8) -> Option<usize> {
    use core::arch::wasm32::*;
    let end = h.len() / 16 * 16;
    let want = u8x16_splat(n);
    let mut i = 0;
    while i < end {
        // SAFETY: i + 16 <= end <= len.
        let blk = unsafe { v128_load(h.as_ptr().add(i).cast()) };
        let m = u8x16_bitmask(u8x16_eq(blk, want)); // one bit per byte, as in movemask
        if m != 0 {
            return Some(i + m.trailing_zeros() as usize);
        }
        i += 16;
    }
    h[end..].iter().position(|&b| b == n).map(|k| end + k)
}

Notice that simd128 is a compile-time feature in Wasm: there is no runtime detection from inside the module. A site that must work in all browsers ships two builds and chooses with a feature probe in JavaScript.

Portable SIMD, std::simd. Rust's standard library has a portable SIMD module that compiles for every target and lowers to the best instructions for each. It is, as of the time of writing (October 2026), still unstable: it needs a nightly compiler and #![feature(portable_simd)], tracked in Rust issue #86656. I ran it on a stable compiler by setting RUSTC_BOOTSTRAP=1, a trick only for experiments. Here is the same byte search:

// src/bin/e12_portable.rs
pub fn find_byte_portable(h: &[u8], n: u8) -> Option<usize> {
    const L: usize = 32;
    let want = u8x32::splat(n);
    let mut chunks = h.chunks_exact(L);
    let mut i = 0;
    for c in &mut chunks {
        let m = u8x32::from_slice(c).simd_eq(want).to_bitmask();
        if m != 0 { return Some(i + m.trailing_zeros() as usize); }
        i += L;
    }
    chunks.remainder().iter().position(|&b| b == n).map(|k| i + k)
}

It is almost the same code as the AVX2 one, with no unsafe, no cfg, and no feature checks. u8x32::simd_eq(...) gives a mask type, and .to_bitmask() turns it into an integer. The tail is chunks_exact(...).remainder(). The program asserts that it agrees with the scalar version for every length from 0 to 199. And the speed, against the hand-written AVX2 version:

portable version agrees with scalar for all lengths 0..200

     size |  std::simd  AVX2 arch   (GB/s)
       64 |       34.9       21.4
     4096 |       36.7       61.3
  1048576 |       41.0       52.6

At 64 bytes the portable version is faster (34.9 vs 21.4 GB/s). At 4 KiB it is slower (36.7 vs 61.3). At 1 MiB it is close (41.0 vs 52.6). My portable version is a straight 32-byte loop with one branch per vector, comparable to the plain AVX2 kernel and not the unrolled one, and that is most of the gap at 4 KiB; I did not try to optimize it. The point is that almost no code was needed and no unsafe appeared. When std::simd stabilizes, this will be the default choice for code that only needs the common operations.

A caution on "portable". The module's documentation notes that in some cases the architecture-specific intrinsic is faster or slightly different in behaviour than the portable equivalent, and that you can convert between the two when you need the special instruction. Portable meaning does not mean portable cost. A shuffle that is one instruction on one chip may be several on another.

Choose your layer. Here is the decision table:

Build it (the "how")

How should you evaluate a SIMD dependency before trusting it? This is a checklist, not a recipe:

1. Does it run on your targets? Check x86-64, AArch64 and Wasm, and what it does without SIMD (a fallback).
2. Is it correct? Does it have a scalar oracle and differential tests like the ones in Chapter 8?
3. What is its unsafe surface? Read the // SAFETY: comments in the parts you use.
4. Does it dispatch at runtime, or require a build flag? If a build flag, you may be breaking Chapter 9's native rule.
5. Is it maintained, and does it support a current stable compiler?
6. Is it faster on your data? Benchmark at the sizes and shapes you actually have, against your scalar version and against memchr-class libraries where one exists.

Do this before you write your own kernel. A mature library will often beat the first version you write (we saw memchr ahead of our kernels at 64 KiB in Chapter 11), and it comes with years of bug fixes.

Why it is this way (the "why")

Why does each architecture have its own SIMD, instead of one standard?

Because instruction sets grew from different needs and different decades. x86 grew by adding width again and again (64 → 128 → 256 → 512 bits), keeping every old instruction forever for compatibility. Arm designed a cleaner 128-bit set late, and then Scalable Vector Extension (SVE) in which the width is not fixed: code is written once and runs on chips with 128 to 2048-bit vectors, with predicates handling the loop tail. RISC-V's vector extension makes the same choice. WebAssembly picked a fixed 128 bits as the smallest common denominator that every host can implement.

Each choice trades something. Fixed width is simple to program and to optimize for, and ties the code to a width. Scalable width is portable and future-proof, and harder for compilers and programmers to reason about (you cannot say "8 lanes" at compile time). A portable library has to pick a model, and std::simd picked fixed-size vectors like f32x8, with the compiler lowering them to whatever the target has.

What would have to be true for one standard to work? The hardware would have to agree on semantics that are cheap everywhere. They do not: x86 pshufb zeroes a lane if the index's top bit is set, NEON's equivalent vqtbl1q_u8 also zeroes out-of-range lanes but with a different range rule, and Wasm's i8x16.swizzle too, with its own rule. The exact meaning of "an index out of range" differs, so a portable shuffle must either pick one meaning and sometimes pay extra instructions, or leave it unspecified. Real portability has real cost.

Common mistakes

• Writing x86 intrinsics without cfg(target_arch) guards, which breaks the ARM build.
• Assuming NEON has movemask. It does not.
• Assuming Wasm SIMD can be detected at runtime inside the module. It cannot.
• Calling a "portable" API fast everywhere. Measure on each target.
• Depending on nightly std::simd in a stable project without a pinned toolchain plan.

Say it simply

The same four moves (lane-wise, mask, select, shuffle) exist on every platform, with different names and different costs. NEON is 128-bit and has no movemask. WebAssembly SIMD is fixed 128-bit and chosen at build time. std::simd is a portable layer that is still unstable. Pick the layer by what you must support and what you are willing to maintain, and check a mature library before writing your own.

Check yourself

Quick questions

1. (Warm-up, explain) Why does NEON need the "shift right narrow by 4" trick, and why do we divide trailing_zeros() by 4?
2. (Core, explain) What is the difference between a fixed-width and a scalable-width vector design? Name one platform for each.
3. (Core, explain) Why is "portable SIMD" not the same as "portable performance"?
4. (Warm-up, trace) In the std::simd code, which line plays the role of movemask?

Problems

5. (Core, implement, ARM or Wasm only) If you have an ARM machine or a Wasm toolchain, compile the snippet above, write a differential test (Chapter 8), and report whether it was right the first time. If not, what was wrong? (If you do not have either, do this on paper: trace find_byte_neon for a 16-byte input with the needle at index 5. What is m in hex?)
6. (Core, implement) Write sum_f32 with std::simd (f32x8 accumulators) and compare with sum_f32_avx2_4acc for 4 K floats. Is the portable version close? How many accumulators did you use?
7. (Stretch, break it) In the portable find_byte, replace u8x32 with u8x64. On your CPU does it compile to AVX2, AVX-512 or something else? Check with objdump.

Think harder

8. (Stretch, design and defend) You are writing a library that must run on x86-64 servers, Apple laptops, and in the browser, and needs a fast byte scan. Choose: three hand-written kernels, std::simd on nightly, or a crate. Defend your choice and say what would change it.

Part 3 wrap-up

What you can now do

• Check whether the compiler vectorized a loop, and with what width.
• Explain the three contracts and write a kernel that obeys them.
• Build byte kernels (find, count, sum, clamp) and test them against an oracle.
• Choose a kernel at runtime with a cached pointer and a size cutoff.
• Use AVX-512 masks to write a tail with no scalar code, and prove it is safe.
• Say how the same ideas look on NEON, WebAssembly and std::simd.

Mini project: bytes-kit (4 to 6 hours)

Build a small library with four functions: find_byte, count_byte, sum_u8, clamp_u8. For each: a scalar oracle, an AVX2 kernel, an AVX-512 kernel if you have one, a safe dispatching wrapper with a size cutoff, and differential tests over all lengths 0 to 300 and offsets 0, 1, 7, 31, 63.

Done when: all tests pass; a benchmark table shows scalar, your best kernel and memchr (for find_byte) at 16 B, 256 B, 4 KiB, 1 MiB and 64 MiB; and you can say, for each size, whether you are limited by overhead, instructions or memory.

Write-up: half a page. Where did you beat memchr, where not, and why?

Reading list

• The Rust docs for core::arch (doc.rust-lang.org/core/arch), especially the pages for x86-64 and AArch64, and the Intel Intrinsics Guide (intel.com/content/www/us/en/docs/intrinsics-guide) for the meaning of each x86 intrinsic.
• Rust RFC 2396, target_feature 1.1. The design behind safe #[target_feature] functions.
• Rust tracking issue #86656, for the state of portable SIMD, and the portable-simd repository.
• Andrew Gallant (BurntSushi), the memchr crate source. A professional byte-search implementation to read after this Part.
• Stephens et al., The ARM Scalable Vector Extension, IEEE Micro, 2017, for the scalable-vector idea.

Teach it back. In five sentences, explain to a friend how one binary can run a fast AVX-512 kernel on a new server and a slower plain loop on an old laptop, and why that is safe.

Answer notes for Part 3

• Ch 9, Q3: f32::max is not reorderable the way integer max is, because of NaN handling, so the compiler needs care. Whether it vectorized depends on how you wrote the reduction. Check, don't guess.
• Ch 10, Q3: a compiles. b compiles: it has the same feature as a, so calling a is safe. c fails with E0133 because c has no avx2.
• Ch 11, Q1: returns Some(40). It inspects 2 blocks: the first (bytes 0 to 31) has movemask == 0, the second (bytes 32 to 63) has bit 8 set. Q2: a u8 lane wraps after 255: 300 mod 256 = 44. The flush every 255 blocks keeps each lane at most 255.
• Ch 12, Q1: 100,000,000 × 1 ns = 0.1 s. Probably not worth removing unless the whole job is tiny. If those calls are all under 32 bytes, add a cutoff instead.
• Ch 13, Q1: 1 full vector (64 bytes), 36 tail bytes, k = (1<<36)-1 = 36 ones in binary.

Part 4. Algorithms that real systems use

Parts 1 to 3 gave you the machine, the shapes, and the Rust tools. Part 4 puts them together into five algorithms that real systems depend on: matrix multiply (inference), int8 dot products (quantized inference), filtering (databases), byte classification (parsers), and bit-parallel scanning (text formats). Each chapter starts from a naive version and ends with a measured, tested kernel.

Chapter 15. Dot products and matrix multiply

In this chapter: you build a dot product and a matrix multiply, watch a textbook loop go from 1.8 to 54 GFLOP/s, and learn the one idea that makes every fast matrix multiply work: keep the accumulators in registers.

You need first: Chapters 2, 4 and 5.

The problem first

A neural network layer is a matrix multiply. If you have written an MNIST classifier from scratch, one training step is mostly this loop:

for i in 0..m { for j in 0..n { for p in 0..k { c[i*n + j] += a[i*k + p] * b[p*n + j]; } } }

It does 2·m·n·k floating-point operations. For three 512×512 matrices that is 268 million operations. How fast it runs decides how long your training takes. This one function is where most of the world's AI compute goes.

The model (the "what")

Start with the simpler dot product: sum of a[i] * b[i]. It is a reduction (Chapter 6) with a multiply. On a CPU with FMA, one instruction does the multiply and the add, so the kernel is the 4-accumulator FMA loop from Chapter 2:

// src/linalg.rs
/// Dot product with 4 independent FMA chains.
#[target_feature(enable = "avx2,fma")]
pub fn dot_f32_avx2_fma(a: &[f32], b: &[f32]) -> f32 {
    assert_eq!(a.len(), b.len());
    let end = a.len() / 32 * 32;
    let mut acc = [_mm256_setzero_ps(); 4];
    let mut i = 0;
    while i < end {
        // SAFETY: i + 32 <= end <= len for both slices.
        unsafe {
            for q in 0..4 {
                let x = _mm256_loadu_ps(a.as_ptr().add(i + 8 * q));
                let y = _mm256_loadu_ps(b.as_ptr().add(i + 8 * q));
                acc[q] = _mm256_fmadd_ps(x, y, acc[q]);
            }
        }
        i += 32;
    }
    let s = _mm256_add_ps(_mm256_add_ps(acc[0], acc[1]), _mm256_add_ps(acc[2], acc[3]));
    let mut r = crate::sums::x86::hsum256(s);
    for q in end..a.len() {
        r += a[q] * b[q];
    }
    r
}