Finding and Filtering with xsimd

In the previous lesson, we used SIMD to perform some basic math. Vectorizing arithmetic is useful, but for many data structures and algorithms - like searching an array or probing a hash table - we aren't doing math. We are searching and branching.

How do you write an if statement when a single register holds 32 distinct numbers simultaneously? In scalar code, the == operator evaluates two values and returns a single bool (true or false), which the CPU uses to jump to a different instruction address.

In SIMD code, the CPU cannot jump in 32 different directions at once. Instead, a SIMD comparison is a branchless mask operation, reusing the bitwise techniques and branchless logic we introduced previously.

For example, if we compare two xsimd::batch<int8_t> registers using ==, the CPU compares all lanes in parallel and outputs a batch_bool<int8_t> register:

If a lane matches, that lane's corresponding bits in the output register are all set to 1. Given our input type is 8 bits in this example (int8_t), we'd have a sequence of 8 bits set to 1s in the output register - an 0xFF byte.

If a lane does not match, its corresponding bits in the output register are all set to 0 - an 0x00 byte in our int8_t example.

We're showing this with 128-bit registers in these examples. These let us check 16 int8_t values concurrently, but the same process applies across any number of lanes:

The following example shows this mechanism in action:

1#include <iostream>
2#include <vector>
3#include <cstdint>
4#include <xsimd/xsimd.hpp>
5
6int main() {
7  using Batch = xsimd::batch<int8_t>;
8  using MaskBatch = xsimd::batch_bool<int8_t>;
9
10  std::vector<int8_t> Buffer(Batch::size, 0);
11  Buffer[3] = 42;
12  Buffer[11] = 42;
13
14  Batch DataBatch{Batch::load_unaligned(Buffer.data())};
15  Batch TargetBatch{Batch::broadcast(42)};
16
17  // Compare all lanes in parallel without branching
18  MaskBatch Matches{DataBatch == TargetBatch}; 
19
20  for (size_t i{0}; i < 16; ++i) {
21    std::cout << Matches.get(i) << " ";
22  }
23}

10 0 0 1 0 0 0 0 0 0 0 1 0 0 0 0

Compressing to a Bitmask

Having a large vector register full of 0xFF and 0x00 bytes is great, but to work with that output using standard C++ logic, we need to pull that information back into a simple scalar register. We do this by calling mask() on the batch_bool result.

Under the hood, mask() invokes a dedicated hardware instruction (like vpmovmskb on x86) that tightly packs our result into a standard integer bitmask. This output has predictable properties - for example:

• If Bit 0 of our integer mask is 1, then Lane 0 matched.
• If Bit 5 of our integer mask is 1, then Lane 5 matched.
• If the entire integer mask is 0, then none of the elements matched.

Note that binary sequences and array indicies tend to be visualized in opposite orders. Arrays are usually represented from left-to-right, where index 0 is the leftmost element, whilst bits are usually represented right-to-left with bit 0 being the rightmost bit.

The mask() function handles this remapping for us, so it's not something we need to worry about:

Putting this all together, this sequence allows us to check many array slots (usually 16, 32, or 64 depending on the hardware) in a single clock cycle, without a single branch.

Then, by checking if (mask != 0) once, we replace 16, 32, or 64 unpredictable loop branches with one consolidated check to determine whether our value was found:

1#include <iostream>
2#include <vector>
3#include <cstdint>
4#include <bitset> 
5#include <xsimd/xsimd.hpp>
6
7int main() {
8  using Batch = xsimd::batch<int8_t>;
9  using MaskBatch = xsimd::batch_bool<int8_t>;
10
11  std::vector<int8_t> Buffer(Batch::size, 0);
12  Buffer[3] = 42;
13  Buffer[11] = 42;
14
15  Batch DataBatch{Batch::load_unaligned(Buffer.data())};
16  Batch TargetBatch{Batch::broadcast(42)};
17
18  MaskBatch Matches{DataBatch == TargetBatch};
19
20  // Compress the vector result into a scalar integer bitmask
21  uint64_t mask{Matches.mask()}; 
22
23  std::cout << "Bitmask: " << std::bitset<16>{mask} << "\n";
24
25  if (mask != 0) { 
26    std::cout << "Match found in this batch!";
27  }
28}

1Bitmask: 0000100000001000
2Match found in this batch!

Finding the Exact Index

Once we have a non-zero uint64_t bitmask - say, 00001000 in binary, where Bit 3 is set to 1 - how do we translate that bit into the array index 3 without writing a loop to check each bit?

We reach back into C++20's <bit> header, which we introduced earlier in the course and used again when implementing HyperLogLog .

While HyperLogLog used std::countl_zero() to count leading zeros from the most significant bit, here we use std::countr_zero() to count trailing zeros starting from the least significant (rightmost) bit. For example:

• If Lane 0 matched, Bit 0 is 1 (...00000001). There are 0 trailing zeros to the right of Bit 0 so std::countr_zero(mask) returns 0.
• If Lane 3 matched, Bit 3 is 1 (...00001000). There are 3 trailing zeros (Bits 0, 1, and 2) to the right of Bit 3 so std::countr_zero(mask) returns 3.

On x86 hardware, std::countr_zero() compiles to a single-cycle tzcnt (Trailing Zero Count) instruction. There is zero looping - the silicon hands us the exact matching lane index in one clock cycle.

Let's combine xsimd::batch<int8_t>, mask(), and std::countr_zero() to build a SimdFind() function in SimdAlgorithms.h.

Our function will accept a pointer to an array of int8_t bytes, the length of the array, and a target byte to search for. It will return the index of the first match, or count if the target does not exist in the array:

1#pragma once
2#include <xsimd/xsimd.hpp>
3#include <cstddef>
4#include <cstdint>
5#include <bit>
6
7namespace dsa {
8
9inline size_t SimdFind(
10  const int8_t* data,
11  size_t count,
12  int8_t target
13) {
14  using Batch = xsimd::batch<int8_t>;
15  constexpr size_t Lanes = Batch::size;
16
17  // Broadcast our target byte into all 16 or 32 lanes
18  Batch targetBatch(target);
19
20  size_t i{0};
21
22  // Step through the array 32 bytes at a time
23  for (; i + Lanes <= count; i += Lanes) {
24    Batch chunk = Batch::load_unaligned(data + i);
25
26    // Compare all 32 bytes in a single cycle and
27    // extract the 32-bit integer bitmask
28    uint64_t mask = (chunk == targetBatch).mask();
29
30    // Perform a single branch for all 32 elements
31    if (mask != 0) {
32      // Find the exact lane that matched using tzcnt
33      size_t laneOffset = std::countr_zero(mask);
34      return i + laneOffset;
35    }
36  }
37
38  // Scalar Tail Loop for any remaining bytes
39  for (; i < count; ++i) {
40    if (data[i] == target) return i;
41  }
42
43  return count;
44}
45
46}

When we use this to search an array of 1,024 bytes on a processor with 256-bit (32-byte) registers, instead of executing 1,024 loop iterations with 1,024 conditional branches, the CPU executes at most 32 iterations.

A single 64-byte L1 cache line holds two entire 32-byte batches. The CPU pulls one cache line from RAM and searches the entire 64-byte block in two instructions.

As soon as mask != 0 triggers, std::countr_zero(mask) pinpoints the exact byte offset of the first match, adds it to our base index i, and returns immediately.

Using the function would look like this:

1#include <iostream>
2#include <vector>
3#include <cstdint>
4#include <dsa/SimdAlgorithms.h> 
5
6int main() {
7  std::vector<int8_t> Buffer(16, 0);
8  Buffer[3] = 42;
9  Buffer[11] = 42;
10
11  size_t Index{
12    dsa::SimdFind(Buffer.data(), Buffer.size(), 42
13  )}; 
14
15  if (Index != Buffer.size()) {
16    std::cout << "Match found at index " << Index << "\n";
17  } else {
18    std::cout << "No match found\n";
19  }
20}

Finding the first match with std::countr_zero(mask) is enough for a basic searching algorithm. However, what if multiple slots inside the same 32-byte batch match our target, and we want to find all of them?

When we use these techniques to build a hash table later in the chapter, we still have to deal with the classic problem where different inputs hash to the same value.

When we search for that hash, our comparison will produce a bitmask where multiple lanes match, for example:

1// Both lanes 2 and 5 match
2// Binary: 00100100 (Bits 2 and 5 are set)
3uint64_t mask = 0b00100100;

When we check Lane 2 in our hash table, we might discover it was just a hash collision (a false positive), meaning the actual full key we are looking might be sitting in Lane 5.

Therefore, we need a fast way to iterate through every 1 bit inside mask one by one, and we don't want to waste cycles scanning the 0 bits in between them.

Kernighan's Bit-Clearing

We can do this using a classic bitwise algorithm known as Kernighan's Bit-Clearing Trick:

1mask &= (mask - 1);

We explain why this bitwise logic works below, but let's first see an example of it in action.

In this case, we'll imagine mask() returns 00100100, representing matches on lanes 2 (the 1 on the right) and 5 (the 1 on the left). We execute std::countr_zero() and mask &= (mask - 1) in a loop to retrieve all the matches:

1#include <iostream>
2#include <cstdint>
3#include <bit>
4#include <bitset>
5
6int main() {
7  // Binary: 00100100 (Bits 2 and 5 are set)
8  uint64_t mask{0b00100100};
9
10  while (mask != 0) {
11    size_t lane{
12      static_cast<size_t>(std::countr_zero(mask))
13    };
14
15    std::cout << "Found match at lane: " << lane
16      << " | Mask: " << std::bitset<8>{mask} << "\n";
17
18    // Erase the lowest set 1-bit
19    mask &= (mask - 1);
20  }
21}

1Found match at lane: 2 | Mask: 00100100
2Found match at lane: 5 | Mask: 00100000

This trick works because, when we subtract 1 from any binary number, it flips the lowest 1 bit to 0, and flips every 0 bit to the right of it to 1.

When we bitwise AND (&) that result back with the original mask, it erases the lowest set 1 bit while leaving every higher 1 bit untouched.

Here is a step-by-step walkthrough of our code:

Because CPUs execute this pattern so frequently, x86 processors even have a single-cycle hardware instruction - blsr - that the compiler will use whenever it sees mask &= (mask - 1).

Implementing FindAll()

Let's add a SimdFindAll() function, combining the familiar two-loop strategy of SIMD algorithms with the mask &= (mask - 1) technique to identify all of the indices where our matches were found, not just the first one:

1// ...
2#include <vector> 
3
4namespace dsa {
5
6// ...
7
8// Appends the index of every matching element into 'outIndices'
9inline void SimdFindAll(
10  const int8_t* data,
11  size_t count,
12  int8_t target,
13  std::vector<size_t>& outIndices
14) {
15  using Batch = xsimd::batch<int8_t>;
16  constexpr size_t Lanes = Batch::size;
17
18  Batch targetBatch(target);
19  size_t i{0};
20
21  for (; i + Lanes <= count; i += Lanes) {
22    Batch chunk = Batch::load_unaligned(data + i);
23    uint64_t mask = (chunk == targetBatch).mask();
24
25    // Hop from match to match, skipping all non-matches
26    while (mask != 0) {
27      size_t lane = std::countr_zero(mask);
28      outIndices.push_back(i + lane);
29
30      // Clear the lowest set bit (compiles to single-cycle BLSR)
31      mask &= (mask - 1);
32    }
33  }
34
35  for (; i < count; ++i) {
36    if (data[i] == target) outIndices.push_back(i);
37  }
38}
39
40}

If we use SimdFindAll() to scan a 32-byte chunk that contains 2 matches, the outer loop runs 1 SIMD comparison, and the inner while (mask != 0) loop executes exactly 2 times - hopping directly to the 2 matching indices while ignoring the 30 non-matching bytes.

Finally, let's implement one more common use case - counting how many elements in the collection match our target.

The core structure of this algorithm matches the common two-loop pattern. The only changes from SimdFindAll() to SimdCount() is that we don't need to identify which specific bits in each batch matched - we just need to count how many matched.

This is a simpler problem - we're can just use std::popcount() on each mask to count how many bits are 1. We use that value to update a totalMatches counter that keeps track of the number of matches across all batches:

1// ...
2
3namespace dsa {
4
5// ...
6
7// Counts how many times 'target' appears in 'data'
8inline size_t SimdCount(
9  const int8_t* data,
10  size_t count, 
11  int8_t target
12) {
13  using Batch = xsimd::batch<int8_t>;
14  constexpr size_t Lanes = Batch::size;
15
16  Batch targetBatch(target);
17  size_t totalMatches{0};
18  size_t i{0};
19
20  for (; i + Lanes <= count; i += Lanes) {
21    Batch chunk = Batch::load_unaligned(data + i);
22    uint64_t mask = (chunk == targetBatch).mask();
23
24    // Count all set 1-bits in a single CPU cycle using popcnt
25    totalMatches += std::popcount(mask);
26  }
27
28  for (; i < count; ++i) {
29    if (data[i] == target) ++totalMatches;
30  }
31
32  return totalMatches;
33}
34
35}

Using SimdFind(), SimdCount(), and SimdFindAll() would look something like this:

1#include <iostream>
2#include <vector>
3#include <cstdint>
4#include <dsa/SimdAlgorithms.h>
5
6int main() {
7  // Create a 64-byte buffer filled with 0s
8  std::vector<int8_t> Buffer(64, 0);
9
10  // Plant our target value (42) at specific indices
11  Buffer[5] = 42;
12  Buffer[19] = 42;
13  Buffer[31] = 42;
14  Buffer[50] = 42;
15
16  size_t firstIndex = dsa::SimdFind(
17    Buffer.data(), Buffer.size(), 42
18  );
19  std::cout << "First match at index: "
20    << firstIndex << "\n";
21
22  size_t total = dsa::SimdCount(
23    Buffer.data(), Buffer.size(), 42
24  );
25  std::cout << "Total matches found:  "
26    << total << "\n";
27
28  std::vector<size_t> allMatches;
29  dsa::SimdFindAll(
30    Buffer.data(), Buffer.size(), 42, allMatches
31  );
32
33  std::cout << "All matching indices: ";
34  for (size_t idx : allMatches) {
35    std::cout << idx << " ";
36  }
37}

1First match at index: 5
2Total matches found:  4
3All matching indices: 5 19 31 50

Let's benchmark our SIMD functions against basic scalar implementations. We'll search a 64KB buffer for a target byte near the end of the buffer, and count occurrences across the buffer:

1#include <benchmark/benchmark.h>
2#include <vector>
3#include <algorithm>
4#include <cstdint>
5#include <random>
6#include <dsa/SimdAlgorithms.h>
7
8constexpr size_t kBufferSize = 65'536; // 64 KB
9
10std::vector<int8_t> GenerateSearchBuffer() {
11  std::vector<int8_t> data(kBufferSize, 1);
12  // Place our target byte (42) near the end of the array
13  data[kBufferSize - 13] = 42;
14  return data;
15}
16
17std::vector<int8_t> GenerateCountBuffer() {
18  std::vector<int8_t> data(kBufferSize);
19  std::mt19937 rng(42);
20  std::uniform_int_distribution<int> dist(0, 15);
21  for (size_t i = 0; i < kBufferSize; ++i) {
22    data[i] = static_cast<int8_t>(dist(rng));
23  }
24  return data;
25}
26
27const std::vector<int8_t> SearchData = GenerateSearchBuffer();
28const std::vector<int8_t> CountData = GenerateCountBuffer();
29
30// Scalar Find
31size_t ScalarFind(
32  const int8_t* data,
33  size_t count,
34  int8_t target
35) {
36  for (size_t i = 0; i < count; ++i) {
37    if (data[i] == target) return i;
38  }
39  return count;
40}
41
42// Scalar Count
43size_t ScalarCount(
44  const int8_t* data,
45  size_t count,
46  int8_t target
47) {
48  size_t result = 0;
49  for (size_t i = 0; i < count; ++i) {
50    if (data[i] == target) ++result;
51  }
52  return result;
53}
54
55static void BM_ScalarFind(benchmark::State& state) {
56  for (auto _ : state) {
57    size_t idx = ScalarFind(
58      SearchData.data(), SearchData.size(), 42
59    );
60    benchmark::DoNotOptimize(idx);
61  }
62}
63
64static void BM_SimdFind(benchmark::State& state) {
65  for (auto _ : state) {
66    size_t idx = dsa::SimdFind(
67      SearchData.data(), SearchData.size(), 42
68    );
69    benchmark::DoNotOptimize(idx);
70  }
71}
72
73static void BM_ScalarCount(benchmark::State& state) {
74  for (auto _ : state) {
75    size_t c = ScalarCount(
76      CountData.data(), CountData.size(), 7
77    );
78    benchmark::DoNotOptimize(c);
79  }
80}
81
82static void BM_SimdCount(benchmark::State& state) {
83  for (auto _ : state) {
84    size_t c = dsa::SimdCount(
85      CountData.data(), CountData.size(), 7
86    );
87    benchmark::DoNotOptimize(c);
88  }
89}
90
91BENCHMARK(BM_ScalarFind)->Unit(benchmark::kMicrosecond);
92BENCHMARK(BM_SimdFind)->Unit(benchmark::kMicrosecond);
93BENCHMARK(BM_ScalarCount)->Unit(benchmark::kMicrosecond);
94BENCHMARK(BM_SimdCount)->Unit(benchmark::kMicrosecond);

1------------------------
2Benchmark            CPU
3------------------------
4BM_ScalarFind    16.4 us
5BM_SimdFind     0.900 us
6
7BM_ScalarCount   33.0 us
8BM_SimdCount     1.34 us

In these experiments, our SimdFind() loads 32 bytes at a time into a register, compares all 32 bytes in a single clock cycle, and sweeps through all 65,536 bytes in less than a single microsecond - over 16x faster than the scalar loop.

Similarly, SimdCount() crushes the scalar implementation by combining 32-lane parallel comparisons with the hardware popcnt instruction via std::popcount(), avoiding per-element increments entirely.

Some compilers may auto-vectorize our ScalarCount implementation, which will reduce the performance gap considerably. However, our explicit vectorization approach is still likely to be more efficient, and allows our code to perform more consistently across different toolchains regardless of their auto-vectorization strategy.

Here are the complete files we built in this lesson:

In this lesson, we replaced slow element-by-element branching with parallel SIMD comparisons and bitwise hardware intrinsics:

1. Branchless SIMD Comparisons: Comparing two xsimd::batch registers with == checks all lanes in a single cycle, producing a batch_bool mask filled with 0xFF (1 bits) for matches and 0x00 (0 bits) for misses.
2. Bitmask Compression: Calling mask() collapses the wide vector register down into a compact scalar integer, mapping each vector lane to a single bit.
3. Trailing Zero Scans: Passing a non-zero bitmask to std::countr_zero() invokes the single-cycle tzcnt instruction, translating the lowest matching bit directly into an array index offset.
4. Kernighan's Bit-Clearing: Using mask &= (mask - 1) strips the lowest set 1 bit in a single cycle (blsr), allowing us to iterate exclusively over matching lanes while skipping all non-matches.
5. Hardware Population Counts: Pairing mask() with std::popcount() lets us tally every matching byte across an entire register in one instruction.

Our SimdFindAll() function gives us the ability to inspect 16 or 32 values in a single clock cycle and instantly extract the indices of every match. This allows us to rethink the hash-searching logic we were implementing in the previous chapter.

Rather than trying to reduce lookup times using clever techniques like Robin-Hood hashing, we can instead just use the wide registers and vector instructions of modern hardware to solve the problem using brute force.

This is the idea behind Swiss Tables - the modern standard for hash-based containers - which we'll explore in the next lesson.