Programmers often see regular expressions and the respective matching engines as a complex, almost unattainable blackbox. PCRE spans hundreds of thousands of lines of code. It supports a large and evolving feature set. Lookarounds, backreferences, Unicode properties, and multiple matching modes are all present. On top of this, PCRE includes a JIT compiler that translates patterns into machine code for speed.

While all these features are important for a production-grade regex engine, they can obscure the fundamental concepts behind pattern matching. In this blog post, we will start our journey towards the exploration of the process of building a matcher from scratch, all by ourselves.

1. Exact match ⎇

In this section we focus on exact string matching. The goal is to find occurrences of a given pattern, verbatim within a text. The constrained form of this problem opens up many optimisation tricks and algorithmic techniques that we can leverage.

Much ink has been spilled on this topic and many algorithms have been proposed over the years. Nonetheless, we take a stab at deriving some of the folk wisdom techniques ourselves, bottom-up.

1.1. Exact match of a single character ⎇

The simplest case of exact matching is to find a single character within a text. Programmers are already used to the function strchr from the C standard library, which does exactly that:

#include <cstring>
#include <iostream>

int main() {
  const char * text = "Hello, world!";
  char target = 'o';

  const char * result = std::strchr(text, target);
  if (result) {
    std::cout << "Found at position: " << (result - text) << std::endl;
  } else {
    std::cout << "Not found" << std::endl;
  }
}

/* Output:
Found at position: 4
*/

The function strchr scans the input string text character by character, comparing each character to the target. If a match is found, it returns a pointer to the first occurrence; otherwise, it returns nullptr. A naive implementation of this idea could look like this:

const char * simple_strchr(const char * text, char target) {
  while (*text) {
    if (*text == target) {
      return text;
    }
    ++text;
  }
  return nullptr;
}

A quick skim at the generated assembly code (from clang 21.1 with -O3 optimizations) reveals that the compiler does not bother unrolling the loop or using SIMD instructions, preferring to translate the function almost verbatim to x86-64 assembly:

simple_strchr(char const*, char):
        movzx   ecx, byte ptr [rdi]
        test    cl, cl
        je      .LBB0_1
        mov     rax, rdi
.LBB0_4:
        cmp     cl, sil
        je      .LBB0_2
        movzx   ecx, byte ptr [rax + 1]
        inc     rax
        test    cl, cl
        jne     .LBB0_4
.LBB0_1:
        xor     eax, eax
.LBB0_2:
        ret

The use of C-style strings (i.e., with a null terminator) prevents us from knowing the length of the input string in advance. This limitation hinders our ability to apply certain optimizations, such as processing multiple characters simultaneously using SIMD instructions. One alternative is the memchr function, which takes an explicit length parameter:

int main() {
  const char * text = "Hello, world!";
  size_t length = std::strlen(text);
  char target = 'o';

  const char * result = static_cast<const char *>(std::memchr(text, target, length));
  if (result) {
    std::cout << "Found at position: " << (result - text) << std::endl;
  } else {
    std::cout << "Not found" << std::endl;
  }
}

1.1.1. Naive search. ⎇

In this case, the implementation can take advantage of the known length to optimize the search process. Once again we can implement a simple version of memchr ourselves:

#include <cstddef>

const char * simple_memchr(const char * text, char target, size_t length) {
  for (size_t i = 0; i < length; ++i) {
    if (text[i] == target) {
      return &text[i];
    }
  }
  return nullptr;
}

1.1.2. Search via SWAR. ⎇

Disappointingly, the compiler chose to not vectorise it either. Our first idea for optimisation is to process multiple characters in each iteration. We will load 8 bytes at a time, treating them as a single integer, and comparing them via the bit-xor (^) operation. If any byte matches the target character, the corresponding byte in the result of the xor will be zero. We can then check for the presence and location of the zero byte using some bit-twiddling tricks. Here is the implementation:

// __restrict__: a GNU extension that tells the compiler that
// the pointer is the only way to access the underlying data.
void * memchr2(const void * __restrict__ s, int c, size_t n) {
  const unsigned char * __restrict__ p = (const unsigned char * __restrict__) s;
  const unsigned char uc = (unsigned char)c;
  // Align to 8 bytes (optional, helps on some CPUs; also avoids unaligned loads)
  while (n && ((uintptr_t)p & 7u)) {
    if (*p == uc) return (void *)p;
    ++p;  --n;
  }
  // Broadcast target character to all bytes of a 64-bit word
  // E.g. for 'A' (0x41), pattern becomes 0x4141414141414141
  const uint64_t pattern = 0x0101010101010101ULL * (uint64_t)uc;
  while (n >= 8) {
    // Avoids aliasing + unaligned-load issues.
    // __builtin_memcpy is non-standard and optimized away by the compiler.
    uint64_t w;  __builtin_memcpy(&w, p, sizeof(w));
    uint64_t x = w ^ pattern;
    if (((x - 0x0101010101010101ULL) & ~x & 0x8080808080808080ULL) != 0) {
      // Locate first match within this 8-byte block.
      for (int i = 0; i < 8; ++i) {
        if (p[i] == uc) return (void *)(p + i);
      }
    }
    p += 8;  n -= 8;
  }
  // Tail
  while (n--) {
    if (*p == uc) return (void *)p;
    ++p;
  }
  return NULL;
}

Forfeiting some auxiliary explanations and assuming a sufficient level of C-fu in the reader, we can express this more compactly:

const void * memchr2(const void * __restrict__ s, int c, size_t n) {
  const unsigned char * __restrict__ p = (const unsigned char * __restrict__) s;
  const unsigned char uc = (unsigned char)c;
  for (; n && ((uintptr_t)p & 7u); --n, ++p)
    if (*p == uc) return (const void *)p;
  const uint64_t pattern = 0x0101010101010101ULL * (uint64_t)uc;
  for (; n >= 8; n -= 8, p += 8) {
    uint64_t w;  __builtin_memcpy(&w, p, sizeof(w));
    uint64_t x = w ^ pattern;
    if (((x - 0x0101010101010101ULL) & ~x & 0x8080808080808080ULL))
      for (int i = 0; i < 8; ++i)
        if (p[i] == uc) return (const void *)(p + i);
  }
  for (; n--; ++p) if (*p == uc) return (const void *)p;
  return NULL;
}

Using a quickly hacked together benchmark harness, we compare the three implementations of memchr for throughput on my Ryzen 7 PRO 7840U laptop:

buf_size=67108864 bytes, queries=65536, iters=500, seed=1
glibc_memchr    time=4093754278 ns  checksum=0x06d1490b613edffc
simple_memchr   time=13403055073 ns  checksum=0x06d1490b613edffc
memchr2         time=10869952023 ns  checksum=0x06d1490b613edffc

Estimated throughput (GiB/s):
  glibc_memchr : 259.694
  simple_memchr: 79.320         // 3.27x slower than glibc, 1.22x slower than memchr2
  memchr2      : 97.804         // 2.67x slower than glibc

Skewing the test so that the latency becomes more prominent (smaller buffers) yields:

buf_size=1024 bytes, queries=65536, iters=500, seed=1
glibc_memchr    time=401259363 ns  checksum=0x18c93a5a77e26bd1
simple_memchr   time=1493447878 ns  checksum=0x18c93a5a77e26bd1
memchr2         time=977337739 ns  checksum=0x18c93a5a77e26bd1

Estimated throughput (GiB/s):
  glibc_memchr : 36.016
  simple_memchr: 9.677        // 4x slower than glibc, 1.55x slower than memchr2
  memchr2      : 14.787       // 2.57x slower than glibc

Our SWAR (SIMD Within A Register) implementation of memchr is about 2.5-3x slower than the highly optimized glibc version, which is not bad for a first attempt.

1.1.3. Search via SSE4.2 PCMPESTRI. ⎇

There is still room for improvement, and we will reach for an indispensable trick: x86-64 SIMD intrinsics. We will specifically target the PCMPESTRI instruction from the SSE4.2 instruction set, which essentially does the bulk of our work for us.

To quote the Intel® 64 and IA-32 Architectures Software Developer’s Manual, Volume 2B:

PCMPESTRI — Packed Compare Explicit Length Strings, Return Index

Instruction Summary ⎇

Opcode / Instruction	Op/En	64/32-bit Mode	CPUID Flag	Description
66 0F 3A 61 /r imm8 PCMPESTRI xmm1, xmm2/m128, imm8	RMI	V / V	SSE4_2	Perform a packed comparison of string data with explicit lengths, generate an index, and store the result in ECX.
VEX.128.66.0F3A 61 /r ib VPCMPESTRI xmm1, xmm2/m128, imm8	RMI	V / V	AVX	Perform a packed comparison of string data with explicit lengths, generate an index, and store the result in ECX.

---

Instruction Operand Encoding ⎇

Op/En	Operand 1	Operand 2	Operand 3	Operand 4
RMI	ModRM:reg (r)	ModRM:r/m (r)	imm8	N/A

---

Description ⎇

The instruction compares and processes data from two string fragments based on the encoded value in the imm8 control byte (see Section 4.1, “Imm8 Control Byte Operation for PCMPESTRI / PCMPESTRM / PCMPISTRI / PCMPISTRM”).
It generates an index that is written to the count register ECX.

Each string fragment is represented by two values:

1. Data vector
   1. An xmm register for the first operand.
   2. An xmm register or m128 memory operand for the second operand.
   3. Contains byte or word elements of the string.

1. Length
   1. For xmm1: length is taken from EAX/RAX.
   2. For xmm2/m128: length is taken from EDX/RDX.
   3. The length specifies how many bytes or words are valid in the corresponding vector.

---

Length Interpretation ⎇

• The effective length is the absolute value of the length register.
• The absolute-value computation saturates:
  • Bytes: maximum 16
  • Words: maximum 8
• Saturation occurs when the length register value is:
  • Greater than 16 (or 8 for words), or
  • Less than −16 (or −8 for words)
• Whether data elements are bytes or words is determined by imm8[3].

---

Comparison and Result ⎇

• Comparison and aggregation behavior is controlled by fields in imm8 (see Section 4.1).
• The instruction computes an intermediate result mask (IntRes2, see Section 4.1.4).
• The index of the first or last set bit of IntRes2 is returned in ECX, depending on imm8[6]:
  • If imm8[6] = 0: first (least significant) set bit
  • If imm8[6] = 1: last (most significant) set bit
• If no bits are set in IntRes2:
  • ECX = 16 (byte elements)
  • ECX = 8 (word elements)

We proceed with the implementation in the same vein as the previous SWAR memchr2 function:

#include <nmmintrin.h>

void * memchr_sse42(const void * s, int c, size_t n) {
  const unsigned char * p = (const unsigned char *) s;
  const unsigned char uc = (unsigned char) c;
  const __m128i needle = _mm_cvtsi32_si128((int) uc);
  for (; n && ((uintptr_t)p & 15u); --n, ++p) // Align to 16-byte boundary.
    if (*p == uc) return (void *)p;
  // Compare any byte in hay equals any byte in needle-set (size 1).
  const int imm8 = _SIDD_UBYTE_OPS | _SIDD_CMP_EQUAL_ANY | _SIDD_LEAST_SIGNIFICANT;
  for (; n >= 16; n -= 16, p += 16) {
    const int idx = _mm_cmpestri(needle, 1, _mm_load_si128((const __m128i *) p), 16, imm8);
    if (idx != 16) return (void *)(p + (unsigned)idx);
  }
  for (; n; --n, ++p) if (*p == uc) return (void *)p;
  return nullptr;
}

The results are promising. Our implementation is now only about 2x slower than glibc's highly optimized version:

buf_size=67108864 bytes, queries=65536, iters=500, seed=1
glibc_memchr    time=4129417820 ns  checksum=0xcdc64e6fe915f3ed
simple_memchr   time=13248285362 ns  checksum=0xcdc64e6fe915f3ed
memchr2         time=10671500036 ns  checksum=0xcdc64e6fe915f3ed
memchr_sse42    time=8801114604 ns  checksum=0xcdc64e6fe915f3ed

Estimated throughput (GiB/s):
  glibc_memchr : 257.451
  simple_memchr: 80.246
  memchr2      : 99.623
  memchr_sse42 : 120.794

The source of the remaining gap is a disappointing but frequently recurring story in performance engineering on x86-64 platforms. Looking up the instruction on uops.info reveals the following picture for my Zen4 CPU:

The instruction has a high latency (up to 29 cycles when reading from memory!) and consumes multiple micro-operations (uops) that occupy several execution ports. This means that even though we are processing 16 bytes at a time, the instruction's throughput is limited by its latency and port usage. A natural conclusion would be that glibc clearly uses more advanced techniques that certainly don't involve PCMPESTRI. But why? Wasn't the point of CISC CPU architectures to have complex instructions that do more work per instruction?

This is neither the first nor last time we see that complex, purpose-built instructions do not necessarily lead to better performance. It's quite difficult to wrap one's head around such a paradox. Other examples include:

• REP STOSB is not the fastest way to zero out memory (memset/bzero), despite being a single instruction that does exactly that.
• Even more counterintuitively, seemingly unrolled REP STOSQ is typically even slower than REP STOSB for the same purpose. The reason is in fact simple: the microcode that makes REP STOSB relatively efficient is simply not present for this opcode.
• REP MOVSB is not the fastest way to copy memory (memcpy), even though it is designed for that purpose.
• REPNE SCASB is not the fastest way to search for a byte in memory (memchr), despite being a single instruction that performs that operation. And finally, here we are: another purpose-built instruction, PCMPESTRI, that does not deliver the expected performance boost.

Tangentially, on the topic of REPNE SCASB, let's probe the performance of this. The memchr implementation using this instruction is a simple wrapper over inline assembly:

NOINLINE void * memchr_repne_scas(const void * s, int c, size_t n) {
  if (n == 0) return NULL;
  const unsigned char * p = (const unsigned char *) s;
  void * ret;
  __asm__ volatile(
      "cld\n\t"
      "repne scasb\n\t"
      "jne 1f\n\t"              // if ZF==0 (no match), jump to notfound
      "lea -1(%%rdi), %0\n\t"   // match: ret = rdi-1
      "jmp 2f\n\t"
      "1:\n\t"
      "xor %0, %0\n\t"          // not found: ret = NULL
      "2:\n\t"
      : "=&r"(ret), "+D"(p), "+c"(n)
      : "a"((unsigned char)c)
      : "cc", "memory");
  return ret;
}

As expected, the performance is underwhelming, but but the biggest surprise is that it is even worse than our naive simple_memchr implementation:

buf_size=67108864 bytes, queries=65536, iters=500, seed=1
glibc_memchr    time=4364084308 ns  checksum=0x90e66e3000583c48
simple_memchr   time=14458207884 ns  checksum=0x90e66e3000583c48
memchr2         time=12134840048 ns  checksum=0x90e66e3000583c48
memchr_sse42    time=9476436712 ns  checksum=0x90e66e3000583c48
repne_scas      time=17571507777 ns  checksum=0x90e66e3000583c48

Estimated throughput (GiB/s):
  glibc_memchr : 243.608
  simple_memchr: 73.531
  memchr2      : 87.609
  memchr_sse42 : 112.186
  repne_scas   : 60.503

1.1.4. Search via SSE2 PCMPEQ + PMOVMSKB. ⎇

To go one step further, we can come back to the SSE4.2, PCMPESTRI-based implementation and simply use two instructions with lower latency: PCMPEQ (compare), followed by PMOVMSKB (move mask). This combination allows us to achieve similar functionality with reduced latency and better throughput. We will also forfeit the alignment loop: since we have tied to machine specific code, we are free to use unaligned loads. Here is the implementation:

void * memchr_sse2(const void * s, int c, size_t n) {
  const unsigned char * p = (const unsigned char *)s;
  const unsigned char uc = (unsigned char)c;
  const __m128i needle = _mm_set1_epi8((char) uc);
  for (; n && ((uintptr_t)p & 15u); --n, ++p)
    if (*p == uc) return (void *)p;
  for (; n >= 16; n -= 16, p += 16) {
    __m128i hay = _mm_load_si128((const __m128i *)p);
    __m128i eq  = _mm_cmpeq_epi8(hay, needle);
    unsigned mask = (unsigned)_mm_movemask_epi8(eq);
    if (mask)
      return (void *)(p + __builtin_ctz(mask));
  }
  for (; n--; ++p)
    if (*p == uc) return (void *)p;
  return NULL;
}

To justify this, we come back to the instruction table for Zen 4:

It is clear that this combination of opcodes has better throughput guarantees, throughput and requires less Uops. The execution port usage of the CPU is also better, further suggesting better practical performance. Validation of this result is quite easy:

buf_size=67108864 bytes, queries=65536, iters=500, seed=1
glibc_memchr    time=4966514321 ns  checksum=0x42a2457f5bb57a25
simple_memchr   time=12864543169 ns  checksum=0x42a2457f5bb57a25
memchr2         time=10555408424 ns  checksum=0x42a2457f5bb57a25
memchr_sse42    time=8888208502 ns  checksum=0x42a2457f5bb57a25
repne_scas      time=16022187564 ns  checksum=0x42a2457f5bb57a25
memchr_sse2     time=6822718398 ns  checksum=0x42a2457f5bb57a25

Estimated throughput (GiB/s):
  glibc_memchr : 214.058
  simple_memchr: 82.640
  memchr2      : 100.718
  memchr_sse42 : 119.611
  repne_scas   : 66.353
  memchr_sse2  : 155.821

The soundness of this optimisation is further suported by the outputs of e.g. the LLVM machine code analyser; we can compare the movemask version and the the PCMPESTRI version.

The MCA dumps are for the same overall structure (scalar align/peel into a 16B loop, followed by scalar tail), but with different 16B core, as illustrated by the C code. In terms of assembly code, we see:

• Variant 1 core: vpcmpestri xmm0, [rdi], imm8=0 (SSE4.2 string-compare instruction)
• Variant 2 core: vpcmpeqb + vpmovmskb + tzcnt (SIMD compare + mask + bit-scan)

f99 over ~100 iterations. Every indicator favours the second version:

MCA also shows us some auxiliary information about stalls and port usage. In the core, we see vpcmpestri as annotated with "9 uops, lat 23, RThroughput 4.00". It also loads (MayLoad). The reciprocal throughput of 4.00 is particularly damning: in the best case, we can only issue one such instruction every 4 cycles. This is a hard limit on the loop throughput, regardless of other factors.

1.1.5. Search via AVX2 (improved). ⎇

Naturally increasing the vector width to 256 bits (32 bytes) should yield further performance improvements. This is however not the case, or at least it will not be as straightforward as before.

void * memchr_avx2_simple(const void * s, int c, size_t n) {
  const unsigned char * p = (const unsigned char *)s;
  const unsigned char uc = (unsigned char)c;
  const __m256i needle = _mm256_set1_epi8((char) uc);
  for (; n && ((uintptr_t)p & 31u); --n, ++p)
    if (*p == uc) return (void *)p;
  for (; n >= 32; n -= 32, p += 32) {
    __m256i hay = _mm256_load_si256((const __m256i *)p);
    __m256i eq  = _mm256_cmpeq_epi8(hay, needle);
    unsigned mask = (unsigned)_mm256_movemask_epi8(eq);
    if (mask)
      return (void *)(p + __builtin_ctz(mask));
  }
  for (; n--; ++p)
    if (*p == uc) return (void *)p;
  return NULL;
}

The benchmark results are as follows:

Estimated throughput (GiB/s):
  glibc_memchr : 250.816
  simple_memchr: 87.112
  memchr2      : 97.577
  memchr_sse2  : 155.821
  memchr_avx2_simple : 171.643

The AVX2 version uses the same techniques as the glibc memchr, but we are still quite far away from the 250GiB/s ballpark. The way forward is to capitalise on the broader vector width even more and modify the core loop to scan 128 bytes at a time. This can be achieved by unrolling the loop four times and combining the results using bit-wise OR, scanning for the precise location of the match only when the return is guaranteed:

void * memchr_avx2_simple(const void * s, int c, size_t n) {
  const unsigned char * p = (const unsigned char *)s;
  const unsigned char uc = (unsigned char)c;
  const __m256i needle = _mm256_set1_epi8((char) uc);
  for (; n >= 128; n -= 128, p += 128) {
    __m256i hay1 = _mm256_loadu_si256((const __m256i *)p);
    __m256i hay2 = _mm256_loadu_si256((const __m256i *)(p + 32));
    __m256i hay3 = _mm256_loadu_si256((const __m256i *)(p + 64));
    __m256i hay4 = _mm256_loadu_si256((const __m256i *)(p + 96));
    __m256i eq1  = _mm256_cmpeq_epi8(hay1, needle);
    __m256i eq2  = _mm256_cmpeq_epi8(hay2, needle);
    __m256i eq3  = _mm256_cmpeq_epi8(hay3, needle);
    __m256i eq4  = _mm256_cmpeq_epi8(hay4, needle);
    __m256i any  = _mm256_or_si256(_mm256_or_si256(eq1, eq2), _mm256_or_si256(eq3, eq4));
    unsigned mask = (unsigned)_mm256_movemask_epi8(any);
    if (mask) {
      unsigned mmask1 = (unsigned)_mm256_movemask_epi8(eq1);
      if (mmask1) return (void *)(p + __builtin_ctz(mmask1));
      unsigned mmask2 = (unsigned)_mm256_movemask_epi8(eq2);
      if (mmask2) return (void *)(p + 32 + __builtin_ctz(mmask2));
      unsigned mmask3 = (unsigned)_mm256_movemask_epi8(eq3);
      if (mmask3) return (void *)(p + 64 + __builtin_ctz(mmask3));
      unsigned mmask4 = (unsigned)_mm256_movemask_epi8(eq4);
      return (void *)(p + 96 + __builtin_ctz(mmask4));
    }
  }
  for (; n >= 32; n -= 32, p += 32) {
    __m256i hay = _mm256_loadu_si256((const __m256i *)p);
    __m256i eq  = _mm256_cmpeq_epi8(hay, needle);
    unsigned mask = (unsigned)_mm256_movemask_epi8(eq);
    if (mask)
      return (void *)(p + __builtin_ctz(mask));
  }
  for (; n--; ++p)
    if (*p == uc) return (void *)p;
  return NULL;
}

This almost brings us to parity with glibc on large buffers:

buf_size=67108864 bytes, queries=65536, iters=500, seed=1
glibc_memchr    time=3884077381 ns  checksum=0xf3bb7462ab9aa0ae
simple_memchr   time=12991746698 ns  checksum=0xf3bb7462ab9aa0ae
memchr2         time=10150611634 ns  checksum=0xf3bb7462ab9aa0ae
memchr_sse2     time=6738995996 ns  checksum=0xf3bb7462ab9aa0ae
memchr_avx2_simple  time=4027457053 ns  checksum=0xf3bb7462ab9aa0ae

Estimated throughput (GiB/s):
  glibc_memchr : 273.713
  simple_memchr: 81.831
  memchr2      : 104.735
  memchr_sse2  : 157.757
  memchr_avx2_simple : 263.969

1.1.6. Search via SSE2 + BMI (improved). ⎇

The previous AVX2 version uses unaligned loads everywhere (_mm256_loadu_si256), and relies on modern CPUs handling that well enough. This is not a great assumption to make for SSE2 code - if your CPU is old enough to not support AVX2, it is likely to have worse performance for unaligned loads. We can therefore revisit the SSE2 version and improve it by adding proper alignment handling.

The new version will spend a lot of code up front to avoid bad alignment situations and to enable aligned loads (_mm_load_si128) in the hot loops. We will explicitly consider p & 63 (position within a 64-byte cache line), a crosscache path when the initial 16B load would straddle a cache-line boundary (more likely to be slower on some microarchitectures), adjusting p to a 16-byte boundary, then later to a 64-byte boundary (aligned) for the main loop.

The AVX2 code has clean loop tiers: 128B AVX2-unroll, then 32B AVX2, then scalar. Our new version that treats old processors slightly better will do early probing loads that may read past len (still within an aligned block on the same page, which does not trigger a fault) and then checks the found index against len before returning.

For example, if we load 16 bytes even if len < 16, assuming that we don't cross a page boundary (which would cause a segmentation fault), then we can simply do:

idx = ctz(mask); if (idx >= len) return nullptr;

This is a typical speculative read within the same page / valid mapped memory assumption-style trick, often used in high-performance code.

The special cases have crept up and lead us to produce the following code:

static inline __m128i splat_u8_sse2(unsigned char c) {
  __m128i x = _mm_cvtsi32_si128((int)c);
  x = _mm_unpacklo_epi8(x, x);
  x = _mm_unpacklo_epi8(x, x);
  x = _mm_shuffle_epi32(x, 0);
  return x;
}

void * memchr_sse2_bsf(const void * __restrict__  s, int c, size_t n) {
  const unsigned char * __restrict__ p0 = (const unsigned char * __restrict__) s;
  const __m128i needle = splat_u8_sse2((unsigned char)c);
  if (n == 0) return nullptr;
  size_t len = n;
  const unsigned char *p = p0;
  uintptr_t addr = (uintptr_t)p;
  unsigned off64 = (unsigned)(addr & 63u); /* if >48 crosscache */
  if (off64 <= 48u) {
    /* unaligned first 16B */
    __m128i v = _mm_loadu_si128((const __m128i *)(const void *)p);
    __m128i m = _mm_cmpeq_epi8(needle, v);
    unsigned mask = (unsigned)_mm_movemask_epi8(m);
    if (mask) {
      unsigned idx = __builtin_ctz(mask);
      if ((size_t)idx >= len) return nullptr;
      return (void *)(p + idx);
    }
    /* unaligned_no_match_1 */
    if (len <= 16) return nullptr;
    len -= 16;
    const unsigned char *p16 = p + 16;
    unsigned a15 = (unsigned)((uintptr_t)p & 15u);
    p = (const unsigned char *)((uintptr_t)p16 & ~(uintptr_t)15u);
    len += a15;
  } else {
    /* crosscache */
    unsigned off16 = (unsigned)(addr & 15u);
    const unsigned char *pal = (const unsigned char *)(addr & ~(uintptr_t)15u);
    __m128i v = _mm_load_si128((const __m128i *)(const void *)pal);
    __m128i m = _mm_cmpeq_epi8(needle, v);
    unsigned mask = (unsigned)_mm_movemask_epi8(m);
    mask >>= off16;
    if (mask) {
      unsigned idx = __builtin_ctz(mask);
      if ((size_t)idx >= len) return nullptr;
      return (void *)(p0 + idx);
    }
    {
      size_t sum = off16 + len;
      if (sum < off16) sum = (size_t)-1;
      if (sum <= 16) return nullptr;
      len = sum - 16;
    }
    p = pal + 16;
  }
  /* early 2x unroll before 64-byte alignment loop */
  for (;;) {
    if (len <= 64) break;
    len -= 64;
    __m128i v0 = _mm_load_si128((const __m128i *)(const void *)(p + 0));
    __m128i m0 = _mm_cmpeq_epi8(needle, v0);
    unsigned k0 = (unsigned)_mm_movemask_epi8(m0);
    if (k0) return (void *)(p + __builtin_ctz(k0));
    __m128i v1 = _mm_load_si128((const __m128i *)(const void *)(p + 16));
    __m128i m1 = _mm_cmpeq_epi8(needle, v1);
    unsigned k1 = (unsigned)_mm_movemask_epi8(m1);
    if (k1) return (void *)(p + 16 + __builtin_ctz(k1));
    __m128i v2 = _mm_load_si128((const __m128i *)(const void *)(p + 32));
    __m128i m2 = _mm_cmpeq_epi8(needle, v2);
    unsigned k2 = (unsigned)_mm_movemask_epi8(m2);
    if (k2) return (void *)(p + 32 + __builtin_ctz(k2));
    __m128i v3 = _mm_load_si128((const __m128i *)(const void *)(p + 48));
    __m128i m3 = _mm_cmpeq_epi8(needle, v3);
    unsigned k3 = (unsigned)_mm_movemask_epi8(m3);
    if (k3) return (void *)(p + 48 + __builtin_ctz(k3));
    p += 64;
    /* if aligned jump forward, else do another 64 scan */
    if ((((uintptr_t)p) & 63u) == 0) break;
    if (len <= 64) break;
    len -= 64;
    __m128i w0 = _mm_load_si128((const __m128i *)(const void *)(p + 0));
    __m128i n0 = _mm_cmpeq_epi8(needle, w0);
    unsigned q0 = (unsigned)_mm_movemask_epi8(n0);
    if (q0) return (void *)(p + __builtin_ctz(q0));
    __m128i w1 = _mm_load_si128((const __m128i *)(const void *)(p + 16));
    __m128i n1 = _mm_cmpeq_epi8(needle, w1);
    unsigned q1 = (unsigned)_mm_movemask_epi8(n1);
    if (q1) return (void *)(p + 16 + __builtin_ctz(q1));
    __m128i w2 = _mm_load_si128((const __m128i *)(const void *)(p + 32));
    __m128i n2 = _mm_cmpeq_epi8(needle, w2);
    unsigned q2 = (unsigned)_mm_movemask_epi8(n2);
    if (q2) return (void *)(p + 32 + __builtin_ctz(q2));
    __m128i w3 = _mm_load_si128((const __m128i *)(const void *)(p + 48));
    __m128i n3 = _mm_cmpeq_epi8(needle, w3);
    unsigned q3 = (unsigned)_mm_movemask_epi8(n3);
    if (q3) return (void *)(p + 48 + __builtin_ctz(q3));
    p += 64;
    /* align to 64 and carry remainder into len */
    uintptr_t a = (uintptr_t)p;
    unsigned r = (unsigned)(a & 63u);
    p = (const unsigned char *)(a & ~(uintptr_t)63u);
    len += r;
    break;
  }
  /* assumes align of 64 */
  while (len > 64) {
    len -= 64;
    __m128i v0 = _mm_load_si128((const __m128i *)(const void *)(p + 0));
    __m128i v1 = _mm_load_si128((const __m128i *)(const void *)(p + 16));
    __m128i v2 = _mm_load_si128((const __m128i *)(const void *)(p + 32));
    __m128i v3 = _mm_load_si128((const __m128i *)(const void *)(p + 48));
    __m128i m0 = _mm_cmpeq_epi8(needle, v0);
    __m128i m1 = _mm_cmpeq_epi8(needle, v1);
    __m128i m2 = _mm_cmpeq_epi8(needle, v2);
    __m128i m3 = _mm_cmpeq_epi8(needle, v3);
    __m128i t0 = _mm_max_epu8(m0, m2);
    __m128i t1 = _mm_max_epu8(m1, m3);
    __m128i t2 = _mm_max_epu8(t0, t1);
    unsigned any = (unsigned)_mm_movemask_epi8(t2);
    if (any) {
      unsigned k0 = (unsigned)_mm_movemask_epi8(m0);
      if (k0) return (void *)(p + __builtin_ctz(k0));
      unsigned k1 = (unsigned)_mm_movemask_epi8(m1);
      if (k1) return (void *)(p + 16 + __builtin_ctz(k1));
      unsigned k2 = (unsigned)_mm_movemask_epi8(m2);
      if (k2) return (void *)(p + 32 + __builtin_ctz(k2));
      unsigned k3 = (unsigned)_mm_movemask_epi8(m3);
      return (void *)(p + 48 + __builtin_ctz(k3));
    }
    p += 64;
  }
  /* 0 <= len <= 64 */
  {
    if (len == 0) return nullptr;
    if (len > 32) {
      __m128i v0 = _mm_load_si128((const __m128i *)(const void *)(p + 0));
      __m128i m0 = _mm_cmpeq_epi8(needle, v0);
      unsigned k0 = (unsigned)_mm_movemask_epi8(m0);
      if (k0) {
        unsigned idx = __builtin_ctz(k0);
        if ((size_t)idx < len) return (void *)(p + idx);
        return nullptr;
      }
      __m128i v1 = _mm_load_si128((const __m128i *)(const void *)(p + 16));
      __m128i m1 = _mm_cmpeq_epi8(needle, v1);
      unsigned k1 = (unsigned)_mm_movemask_epi8(m1);
      if (k1) {
        unsigned idx = 16 + __builtin_ctz(k1);
        if ((size_t)idx < len) return (void *)(p + idx);
        return nullptr;
      }
      __m128i v2 = _mm_load_si128((const __m128i *)(const void *)(p + 32));
      __m128i m2 = _mm_cmpeq_epi8(needle, v2);
      unsigned k2 = (unsigned)_mm_movemask_epi8(m2);
      if (k2) {
        unsigned idx = 32 + __builtin_ctz(k2);
        if ((size_t)idx < len) return (void *)(p + idx);
        return nullptr;
      }
      if (len <= 48) return nullptr;
      __m128i v3 = _mm_load_si128((const __m128i *)(const void *)(p + 48));
      __m128i m3 = _mm_cmpeq_epi8(needle, v3);
      unsigned k3 = (unsigned)_mm_movemask_epi8(m3);
      if (k3) {
        unsigned idx = 48 + __builtin_ctz(k3);
        if ((size_t)idx < len) return (void *)(p + idx);
        return nullptr;
      }
      return nullptr;
    } else {
      __m128i v0 = _mm_load_si128((const __m128i *)(const void *)(p + 0));
      __m128i m0 = _mm_cmpeq_epi8(needle, v0);
      unsigned k0 = (unsigned)_mm_movemask_epi8(m0);
      if (k0) {
        unsigned idx = __builtin_ctz(k0);
        if ((size_t)idx < len) return (void *)(p + idx);
        return nullptr;
      }
      if (len <= 16) return nullptr;
      __m128i v1 = _mm_load_si128((const __m128i *)(const void *)(p + 16));
      __m128i m1 = _mm_cmpeq_epi8(needle, v1);
      unsigned k1 = (unsigned)_mm_movemask_epi8(m1);
      if (k1) {
        unsigned idx = 16 + __builtin_ctz(k1);
        if ((size_t)idx < len) return (void *)(p + idx);
        return nullptr;
      }
      return nullptr;
    }
  }
}

Benchmarking this version shows that we have attained significantly better throughput:

buf_size=67108864 bytes, queries=65536, iters=500, seed=1
glibc_memchr    time=4525973864 ns  checksum=0x93785bd40474de84
simple_memchr   time=14300729997 ns  checksum=0x93785bd40474de84
memchr2         time=11173968344 ns  checksum=0x93785bd40474de84
memchr_sse2     time=7703221849 ns  checksum=0x93785bd40474de84
memchr_avx2_simple  time=4467790687 ns  checksum=0x93785bd40474de84
memchr_sse2_bsf  time=5780660724 ns  checksum=0x93785bd40474de84

Estimated throughput (GiB/s):
  glibc_memchr : 234.894
  simple_memchr: 74.341
  memchr2      : 95.143
  memchr_sse2  : 138.010
  memchr_avx2_simple : 237.953
  memchr_sse2_bsf : 183.911

1.1.7 Conclusion ⎇

Through a series of optimisations and careful use of SIMD instructions, we have developed a memchr implementation that approaches the performance of glibc's highly optimized version. They also extend to other similar functions like strlen (equivalent to memchr searching for zero) or strchr (equivalent to memchr searching for either a specific character or zero) with analogous techniques.

A straightforward implementation of strlen given the techniques that we have just covered achieves performance that is on-par with glibc's strlen:

size_t strlen_avx2_simple(const void * __restrict__ s) {
  const unsigned char * base = (const unsigned char *) s;
  const unsigned char * p    = base;
  const __m256i zero = _mm256_setzero_si256();
  size_t to_page_end = 4096u - ((uintptr_t)(p) & 4095u);

  for (; to_page_end >= 128u; to_page_end -= 128u, p += 128u) {
    __m256i v1 = _mm256_loadu_si256((const __m256i *)(p + 0));
    __m256i v2 = _mm256_loadu_si256((const __m256i *)(p + 32));
    __m256i v3 = _mm256_loadu_si256((const __m256i *)(p + 64));
    __m256i v4 = _mm256_loadu_si256((const __m256i *)(p + 96));
    __m256i e1 = _mm256_cmpeq_epi8(v1, zero);
    __m256i e2 = _mm256_cmpeq_epi8(v2, zero);
    __m256i e3 = _mm256_cmpeq_epi8(v3, zero);
    __m256i e4 = _mm256_cmpeq_epi8(v4, zero);
    __m256i any = _mm256_or_si256(_mm256_or_si256(e1, e2),
                                  _mm256_or_si256(e3, e4));
    if (_mm256_testz_si256(any, any) == 0) {
      unsigned m1 = (unsigned)_mm256_movemask_epi8(e1);
      if (m1) return (size_t)((p + (size_t)__builtin_ctz(m1)) - base);
      unsigned m2 = (unsigned)_mm256_movemask_epi8(e2);
      if (m2) return (size_t)((p + 32u + (size_t)__builtin_ctz(m2)) - base);
      unsigned m3 = (unsigned)_mm256_movemask_epi8(e3);
      if (m3) return (size_t)((p + 64u + (size_t)__builtin_ctz(m3)) - base);
      unsigned m4 = (unsigned)_mm256_movemask_epi8(e4);
      return (size_t)((p + 96u + (size_t)__builtin_ctz(m4)) - base);
    }
  }

  for (; to_page_end >= 32u; to_page_end -= 32u, p += 32u) {
    __m256i v  = _mm256_loadu_si256((const __m256i *)p);
    __m256i eq = _mm256_cmpeq_epi8(v, zero);
    unsigned mask = (unsigned)_mm256_movemask_epi8(eq);
    if (mask) return (size_t)((p + (size_t)__builtin_ctz(mask)) - base);
  }

  for (; to_page_end; --to_page_end, ++p) {
    if (*p == 0) return (size_t)(p - base);
  }
  
  for (;;) {
    __m256i v1 = _mm256_loadu_si256((const __m256i *)(p + 0));
    __m256i v2 = _mm256_loadu_si256((const __m256i *)(p + 32));
    __m256i v3 = _mm256_loadu_si256((const __m256i *)(p + 64));
    __m256i v4 = _mm256_loadu_si256((const __m256i *)(p + 96));
    __m256i e1 = _mm256_cmpeq_epi8(v1, zero);
    __m256i e2 = _mm256_cmpeq_epi8(v2, zero);
    __m256i e3 = _mm256_cmpeq_epi8(v3, zero);
    __m256i e4 = _mm256_cmpeq_epi8(v4, zero);
    __m256i any = _mm256_or_si256(_mm256_or_si256(e1, e2),
                                  _mm256_or_si256(e3, e4));
    if (_mm256_testz_si256(any, any) == 0) {
      unsigned m1 = (unsigned)_mm256_movemask_epi8(e1);
      if (m1) return (size_t)((p + (size_t)__builtin_ctz(m1)) - base);
      unsigned m2 = (unsigned)_mm256_movemask_epi8(e2);
      if (m2) return (size_t)((p + 32u + (size_t)__builtin_ctz(m2)) - base);
      unsigned m3 = (unsigned)_mm256_movemask_epi8(e3);
      if (m3) return (size_t)((p + 64u + (size_t)__builtin_ctz(m3)) - base);
      unsigned m4 = (unsigned)_mm256_movemask_epi8(e4);
      return (size_t)((p + 96u + (size_t)__builtin_ctz(m4)) - base);
    }
    p += 128;
  }
}

As evidenced by a benchmark in similar vein to the previous ones:

buf_size=67108864 bytes, queries=100000, iters=500, seed=1
strlen_glibc    time=1873349773 ns  gbps=17.911
strlen_avx2_simple  time=2190930772 ns  gbps=15.315

And finally, an implementation of strchr:

char * strchr_avx2_pagesafe(const char * s, int c) {
  const unsigned char * p = (const unsigned char *)s;
  const unsigned char uc = (unsigned char)c;
  const __m256i needle = _mm256_set1_epi8((char)uc);
  const __m256i zero   = _mm256_setzero_si256();
  for (;;) {
    uintptr_t addr = (uintptr_t)p;
    size_t bytes_to_page_end = 4096u - (size_t)(addr & 4095u);
    while (bytes_to_page_end < 32) {
      unsigned char ch = *p;
      if (ch == uc) return (char *)p;
      if (ch == 0)  return NULL;
      ++p;
      ++bytes_to_page_end;
      if (bytes_to_page_end == PAGE_SIZE) break;
    }
    addr = (uintptr_t)p;
    bytes_to_page_end = PAGE_SIZE - (size_t)(addr & PAGE_MASK);
    while (bytes_to_page_end >= 128) {
      __m256i hay1 = _mm256_loadu_si256((const __m256i *)(p +  0));
      __m256i hay2 = _mm256_loadu_si256((const __m256i *)(p + 32));
      __m256i hay3 = _mm256_loadu_si256((const __m256i *)(p + 64));
      __m256i hay4 = _mm256_loadu_si256((const __m256i *)(p + 96));
      __m256i eq1 = _mm256_cmpeq_epi8(hay1, needle);
      __m256i eq2 = _mm256_cmpeq_epi8(hay2, needle);
      __m256i eq3 = _mm256_cmpeq_epi8(hay3, needle);
      __m256i eq4 = _mm256_cmpeq_epi8(hay4, needle);
      __m256i z1  = _mm256_cmpeq_epi8(hay1, zero);
      __m256i z2  = _mm256_cmpeq_epi8(hay2, zero);
      __m256i z3  = _mm256_cmpeq_epi8(hay3, zero);
      __m256i z4  = _mm256_cmpeq_epi8(hay4, zero);
      __m256i m1  = _mm256_or_si256(eq1, z1);
      __m256i m2  = _mm256_or_si256(eq2, z2);
      __m256i m3  = _mm256_or_si256(eq3, z3);
      __m256i m4  = _mm256_or_si256(eq4, z4);
      __m256i any = _mm256_or_si256(_mm256_or_si256(m1, m2), _mm256_or_si256(m3, m4));
      if (_mm256_testz_si256(any, any) == 0) {
        unsigned comb, eqm, bit;
        comb = (unsigned)_mm256_movemask_epi8(m1);
        if (comb) {
          unsigned idx = (unsigned)__builtin_ctz(comb);
          eqm = (unsigned)_mm256_movemask_epi8(eq1);
          bit = 1u << idx;
          return (eqm & bit) ? (char *)(p + idx) : NULL;
        }
        comb = (unsigned)_mm256_movemask_epi8(m2);
        if (comb) {
          unsigned idx = (unsigned)__builtin_ctz(comb);
          eqm = (unsigned)_mm256_movemask_epi8(eq2);
          bit = 1u << idx;
          return (eqm & bit) ? (char *)(p + 32 + idx) : NULL;
        }
        comb = (unsigned)_mm256_movemask_epi8(m3);
        if (comb) {
          unsigned idx = (unsigned)__builtin_ctz(comb);
          eqm = (unsigned)_mm256_movemask_epi8(eq3);
          bit = 1u << idx;
          return (eqm & bit) ? (char *)(p + 64 + idx) : NULL;
        }
        comb = (unsigned)_mm256_movemask_epi8(m4);
        {
          unsigned idx = (unsigned)__builtin_ctz(comb);
          eqm = (unsigned)_mm256_movemask_epi8(eq4);
          bit = 1u << idx;
          return (eqm & bit) ? (char *)(p + 96 + idx) : NULL;
        }
      }
      p += 128;
      bytes_to_page_end -= 128;
    }
    while (bytes_to_page_end >= 32) {
      __m256i hay = _mm256_loadu_si256((const __m256i *)p);
      __m256i eq  = _mm256_cmpeq_epi8(hay, needle);
      __m256i z   = _mm256_cmpeq_epi8(hay, zero);
      __m256i m   = _mm256_or_si256(eq, z);
      unsigned comb = (unsigned)_mm256_movemask_epi8(m);
      if (comb) {
        unsigned idx = (unsigned)__builtin_ctz(comb);
        unsigned eqm = (unsigned)_mm256_movemask_epi8(eq);
        return (eqm & (1u << idx)) ? (char *)(p + idx) : NULL;
      }
      p += 32;
      bytes_to_page_end -= 32;
    }
  }
}

1.2. Exact match of a pair of characters in overlapping regions. ⎇

For the sake of illustrating the techniques involved, we will consider the problem of searching for either of two characters in a memory block (like a variant of wmemchr or such).

1.2.1. Baseline scalar/SWAR implementations. ⎇

We have a couple of the obvious ones. Simple scalar version with 8-bit loads:

size_t memchr2_portable_simple(const void * ptr, size_t n, uint16_t needle) {
  const uint8_t * p = (const uint8_t *) ptr;
  if (n < 2) return n;
  const uint8_t lo = (uint8_t)needle;
  const uint8_t hi = (uint8_t)(needle >> 8);
  for (size_t i = 0; i + 1 < n; ++i)
    if (p[i] == lo && p[i + 1] == hi) return i;
  return n;
}

Simple scalar version with direct, unaligned 16-bit loads:

NOINLINE size_t memchr2_portable_elim2(const void * ptr, size_t n, uint16_t needle) {
  const uint8_t * p = (const uint8_t *) ptr;
  if (n < 2) return n;
  for (size_t i = 0; i + 1 < n; i++)
    if (*(const uint16_t *)(p + i) == needle) return i;
  return n;
}

An attempt to save a load in each iteration:

size_t memchr2_portable_running(const void * ptr, size_t n, uint16_t needle) {
  const uint8_t * p = (const uint8_t *) ptr;
  if (n < 2) return n;
  uint16_t window = *((const uint16_t *)p);
  if (window == needle) return 0;
  for (size_t i = 0; i + 2 < n; i++) {
    window = (uint16_t)((window >> 8) | ((uint16_t)(p[i + 2]) << 8));
    if (window == needle) return i + 1;
  }
  return n;
}

A SWAR-based version that processes 8 bytes at a time, pre-filtering candidates that contain either of the two bytes in any order in succession:

NOINLINE size_t memchr2_portable_qword(const void * data, size_t len, uint16_t needle) {
  if (len < 2) return len;
  uint8_t b1 = (uint8_t)(needle & 0xFF);
  uint8_t b2 = (uint8_t)(needle >> 8);
  const uint64_t pattern1 = 0x0101010101010101ULL * (uint64_t)b1;
  const uint64_t pattern2 = 0x0101010101010101ULL * (uint64_t)b2;
  const uint8_t * p = (const uint8_t *) data;
  const uint8_t * end = (const uint8_t *) data + len - 1;
  while (p + 7 < end) {
    uint64_t chunk = *(uint64_t *)p;
    uint64_t x1 = chunk ^ pattern1;
    uint64_t x2 = chunk ^ pattern2;
    uint64_t hit1 = (x1 - 0x0101010101010101ULL) & ~x1 & 0x8080808080808080ULL;
    uint64_t hit2 = (x2 - 0x0101010101010101ULL) & ~x2 & 0x8080808080808080ULL;
    uint64_t hc1 = (hit1 >> 8) & hit2;
    uint64_t hc2 = (hit2 >> 8) & hit1;
    if (hc1 | hc2) {
      for (int i = 0; i < 7; i++) {
        uint16_t w = *(const uint16_t *)(p + i);
        if (w == needle) return p + i - (const uint8_t *) data;
      }
    }
    p += 7;
  }
  for (; p < end; p++) {
    uint16_t w = *(const uint16_t *)p;
    if (w == needle) return p - (const uint8_t *) data;
  }
  return len;
}

Checking for both hc1 and hc2 helps to convince the compiler to vectorise the inner loop, leading to - paradoxically - better performance than checking for just one of them.

Also unexpectedly, the SWAR version beats all others in benchmarks by a wide margin, and the second runner-up is the naive, 8-bit direct load version:

buf_size=64000000 bytes, queries=65536, iters=100, seed=666
portable_elim2        time=15880241758 ns  checksum=0xc25f073958310200
scalar                time=14393013166 ns  checksum=0xc25f073958310200
portable_qword        time=8435549272 ns  checksum=0xc25f073958310200
portable_running      time=28003921570 ns  checksum=0xc25f073958310200

Estimated throughput (GiB/s):
  portable_elim2   : 13.499
  scalar           : 14.894
  portable_qword   : 25.412
  portable_running : 7.655

1.2.2. SSE2 and SSSE3 implementations. ⎇

The obvious next step is to try to vectorise the search. Without minding unaligned loads, we can do it the obvious way:

NOINLINE size_t memchr2_sse2_bmi1(const void * ptr, size_t n, uint16_t needle) {
  const uint8_t * p = (const uint8_t *) ptr;
  if (n < 2) return n;
  const uint8_t lo = (uint8_t)needle;
  const uint8_t hi = (uint8_t)(needle >> 8);
  const __m128i vlo = _mm_set1_epi8((char)lo);
  const __m128i vhi = _mm_set1_epi8((char)hi);
  size_t i = 0;
  for (; i + 17 <= n; i += 16) {
    __m128i a = _mm_loadu_si128((const __m128i *)(p + i));
    __m128i b = _mm_loadu_si128((const __m128i *)(p + i + 1));
    __m128i eq0 = _mm_cmpeq_epi8(a, vlo);
    __m128i eq1 = _mm_cmpeq_epi8(b, vhi);
    __m128i both = _mm_and_si128(eq0, eq1);
    unsigned mask = (unsigned)_mm_movemask_epi8(both);
    if (mask)
      return i + (size_t)_tzcnt_u32(mask);   // BMI1
  }
  for (; i + 1 < n; ++i)
    if (p[i] == lo && p[i + 1] == hi) return i;
  return n;
}

The result? A twice speedup over the SWAR version!

buf_size=64000000 bytes, queries=65536, iters=100, seed=666
portable_elim2        time=15880241758 ns  checksum=0xc25f073958310200
scalar                time=14393013166 ns  checksum=0xc25f073958310200
portable_qword        time=8435549272 ns  checksum=0xc25f073958310200
portable_running      time=28003921570 ns  checksum=0xc25f073958310200
sse2_bmi1             time=5253894697 ns  checksum=0xc25f073958310200

Estimated throughput (GiB/s):
  portable_elim2   : 13.499
  scalar           : 14.894
  portable_qword   : 25.412
  portable_running : 7.655
  sse2_bmi1        : 40.802

Using SSSE3 we can sidestep the nasty unaligned loads by using PALIGNR, like so:

NOINLINE size_t memchr2_ssse3_bmi1(const void * ptr, size_t n, uint16_t needle) {
  const uint8_t * p = (const uint8_t *) ptr;
  if (n < 2) return n;
  const uint8_t lo = (uint8_t)needle;
  const uint8_t hi = (uint8_t)(needle >> 8);
  const __m128i vlo = _mm_set1_epi8((char)lo);
  const __m128i vhi = _mm_set1_epi8((char)hi);
  size_t i = 0;
  if (n >= 32) {
    __m128i A = _mm_loadu_si128((const __m128i *)(p + 0));
    for (; i + 32 <= n; i += 16) {
      __m128i N = _mm_loadu_si128((const __m128i *)(p + i + 16));
      __m128i B = _mm_alignr_epi8(N, A, 1); // SSSE3
      __m128i eq0 = _mm_cmpeq_epi8(A, vlo);
      __m128i eq1 = _mm_cmpeq_epi8(B, vhi);
      __m128i both = _mm_and_si128(eq0, eq1);
      unsigned mask = (unsigned)_mm_movemask_epi8(both);
      if (mask) return i + (size_t)_tzcnt_u32(mask); // BMI1
      A = N;
    }
  }
  for (; i + 1 < n; ++i)
    if (p[i] == lo && p[i + 1] == hi) return i;
  return n;
}

Disappointingly, that doesn't seem to win us much:

buf_size=64000000 bytes, queries=65536, iters=100, seed=666
portable_elim2        time=15880241758 ns  checksum=0xc25f073958310200
scalar                time=14393013166 ns  checksum=0xc25f073958310200
portable_qword        time=8435549272 ns  checksum=0xc25f073958310200
portable_running      time=28003921570 ns  checksum=0xc25f073958310200
sse2_bmi1             time=5253894697 ns  checksum=0xc25f073958310200
ssse3_bmi1            time=5354999130 ns  checksum=0xc25f073958310200

Estimated throughput (GiB/s):
  portable_elim2   : 13.499
  scalar           : 14.894
  portable_qword   : 25.412
  portable_running : 7.655
  sse2_bmi1        : 40.802
  ssse3_bmi1       : 40.031

Judging by the output of the profiling tool perf, most of the time is spent in the hot loop, suggesting that nothing funny is going on with the tail or tzcnt again. Zen2 analysis via llvm-mca is also somewhat insubtantial:

Iterations:        100
Instructions:      1400
Total Cycles:      364
Total uOps:        1400

Dispatch Width:    4
uOps Per Cycle:    3.85
IPC:               3.85
Block RThroughput: 3.5

Instruction Info:
[1]    [2]    [3]    [4]    [5]    [6]    Instructions:
 1      1     0.25                        vmovdqa       xmm3, xmm2
 1      8     0.33    *                   vmovdqu       xmm2, xmmword ptr [rdi + r8 + 16]
 1      1     0.50                        vpalignr      xmm4, xmm2, xmm3, 1
 1      1     0.33                        vpcmpeqb      xmm3, xmm0, xmm3
 1      1     0.33                        vpcmpeqb      xmm4, xmm4, xmm1
 1      1     0.25                        vpand xmm3, xmm3, xmm4
 1      1     1.00                        vpmovmskb     esi, xmm3
 1      1     0.25                        test  esi, esi
 1      1     0.25                        jne   .over
 1      1     0.25                        lea   rsi, [r8 + 16]
 1      1     0.25                        add   r8, 48
 1      1     0.25                        cmp   r8, rax
 1      1     0.25                        mov   r8, rsi
 1      1     0.25                        jbe   .star

It also appears that we don't have enough execution ports to unroll and interleave the loop further.

1.2.3. AVX2 implementation. ⎇

We start tame, by widening the SSE vectors to AVX2 256-bit ones. Unfortunately, life is never as simple as one would have hoped. AVX2 likes to include opcodes that operate on 128-bit lanes of 256-bit vectors only. Reportedly this was dictated by microarchitectural considerations in Intel's implementation of AVX2 that are now long obsolete, but we have to deal with it nonetheless.

Our implementation will use permute2x128 and valignr to construct the shifted vectors:

NOINLINE size_t memchr2_avx2_bmi1(const void * ptr, size_t n, uint16_t needle) {
  const uint8_t *p = (const uint8_t *)ptr;
  if (n < 2) return n;
  const uint8_t lo = (uint8_t)needle;
  const uint8_t hi = (uint8_t)(needle >> 8);
  const __m256i vlo = _mm256_set1_epi8((char)lo);
  const __m256i vhi = _mm256_set1_epi8((char)hi);
  size_t i = 0;
  if (n >= 64) {
    __m256i A = _mm256_loadu_si256((const __m256i *)(p + 0));
    for (; i + 64 <= n; i += 32) {
      __m256i N = _mm256_loadu_si256((const __m256i *)(p + i + 32));
      __m256i M = _mm256_permute2x128_si256(A, N, 0x21);
      __m256i B = _mm256_alignr_epi8(M, A, 1);
      __m256i eq0  = _mm256_cmpeq_epi8(A, vlo);
      __m256i eq1  = _mm256_cmpeq_epi8(B, vhi);
      __m256i both = _mm256_and_si256(eq0, eq1);
      unsigned mask = (unsigned)_mm256_movemask_epi8(both);
      if (mask) return i + (size_t)_tzcnt_u32(mask);
      A = N;
    }
  }
  for (; i + 1 < n; ++i)
    if (p[i] == lo && p[i + 1] == hi) return i;
  return n;
}

Instinctively, we will also try to shove a couple of _mm256_movemask_epi8 away from the hot path. To do this we will apply the same optimisation idea as in the regular memchr implementation:

static inline __m256i both_from_AN(__m256i A, __m256i N, __m256i vlo, __m256i vhi) {
  __m256i M = _mm256_permute2x128_si256(A, N, 0x21);
  __m256i B = _mm256_alignr_epi8(M, A, 1);
  __m256i eq0  = _mm256_cmpeq_epi8(A, vlo);
  __m256i eq1  = _mm256_cmpeq_epi8(B, vhi);
  return _mm256_and_si256(eq0, eq1);
}
NOINLINE size_t memchr2_avx2_bmi1_alt(const void * ptr, size_t n, uint16_t needle) {
  const uint8_t * p = (const uint8_t *) ptr;
  if (n < 2) return n;
  const uint8_t lo = (uint8_t)needle;
  const uint8_t hi = (uint8_t)(needle >> 8);
  const __m256i vlo = _mm256_set1_epi8((char)lo);
  const __m256i vhi = _mm256_set1_epi8((char)hi);
  size_t i = 0;
  if (n >= 96) {
    __m256i A = _mm256_loadu_si256((const __m256i *)(p + 0));
    for (; i + 96 <= n; i += 64) {
      __m256i N0 = _mm256_loadu_si256((const __m256i *)(p + i + 32));
      __m256i N1 = _mm256_loadu_si256((const __m256i *)(p + i + 64));
      __m256i both0 = both_from_AN(A,  N0, vlo, vhi);
      __m256i both1 = both_from_AN(N0, N1, vlo, vhi);
      __m256i any = _mm256_or_si256(both0, both1);
      if (!_mm256_testz_si256(any, any)) {
        unsigned m0 = (unsigned)_mm256_movemask_epi8(both0);
        if (m0) return i + (size_t)_tzcnt_u32(m0);
        unsigned m1 = (unsigned)_mm256_movemask_epi8(both1);
        return (i + 32) + (size_t)_tzcnt_u32(m1);
      }
      A = N1;
    }
    for (; i + 64 <= n; i += 32) {
      __m256i N = _mm256_loadu_si256((const __m256i *)(p + i + 32));
      __m256i both = both_from_AN(A, N, vlo, vhi);
      if (!_mm256_testz_si256(both, both)) {
        unsigned m = (unsigned)_mm256_movemask_epi8(both);
        return i + (size_t)_tzcnt_u32(m);
      }
      A = N;
    }
  } else if (n >= 64) {
    __m256i A = _mm256_loadu_si256((const __m256i *)(p + 0));
    for (; i + 64 <= n; i += 32) {
      __m256i N = _mm256_loadu_si256((const __m256i *)(p + i + 32));
      __m256i both = both_from_AN(A, N, vlo, vhi);
      if (!_mm256_testz_si256(both, both)) {
        unsigned m = (unsigned)_mm256_movemask_epi8(both);
        return i + (size_t)_tzcnt_u32(m);
      }
      A = N;
    }
  }
  for (; i + 1 < n; ++i)
      if (p[i] == lo && p[i + 1] == hi) return i;
  return n;
}

Benchmarks show that the second AVX2 version is about the same in terms of performance as the first one, except it seems to stall on the vtest instruction on both Zen2 and Zen4 CPUs. A test on a Zen2 CPU (5900X) gives:

buf_size=64000000 bytes, queries=65536, iters=50, seed=16
portable_elim2        time=7122549226 ns  checksum=0x9ea45cc56838615c
scalar                time=3940411071 ns  checksum=0x9ea45cc56838615c
portable_qword        time=1868259066 ns  checksum=0x9ea45cc56838615c
portable_running      time=11036034203 ns  checksum=0x9ea45cc56838615c
sse2_bmi1             time=774990423 ns  checksum=0x9ea45cc56838615c
sse2_bmi1_fast        time=837601040 ns  checksum=0x9ea45cc56838615c
ssse3_bmi1            time=821152775 ns  checksum=0x9ea45cc56838615c
avx2_bmi1             time=710584728 ns  checksum=0x9ea45cc56838615c
avx2_bmi1_alt         time=703469492 ns  checksum=0x9ea45cc56838615c

Estimated throughput (GiB/s):
  portable_elim2   : 14.727
  scalar           : 26.621
  portable_qword   : 56.147
  portable_running : 9.505
  sse2_bmi1        : 135.352
  sse2_bmi1_fast   : 125.234
  ssse3_bmi1       : 127.743
  avx2_bmi1        : 147.620
  avx2_bmi1_alt    : 149.113

The kernel of the improved AVX2 loop appears to be:

   3.81 │ 40:   vmovdqu      ymm4,YMMWORD PTR [rdi+rsi*1+0x20]
   2.90 │       vpcmpeqb     ymm5,ymm0,ymm2
   4.17 │       vperm2i128   ymm3,ymm2,ymm4,0x21
   3.16 │       vpalignr     ymm3,ymm3,ymm2,0x1
   2.78 │       vmovdqu      ymm2,YMMWORD PTR [rdi+rsi*1+0x40]
   3.27 │       vpcmpeqb     ymm3,ymm3,ymm1
   1.95 │       vpand        ymm3,ymm5,ymm3
   2.49 │       vperm2i128   ymm5,ymm4,ymm2,0x21
   1.91 │       vpalignr     ymm5,ymm5,ymm4,0x1
   2.32 │       vpcmpeqb     ymm4,ymm0,ymm4
   3.60 │       vpcmpeqb     ymm5,ymm5,ymm1
   3.15 │       vpand        ymm4,ymm4,ymm5
   5.93 │       vpor         ymm5,ymm3,ymm4
  15.79 │       vptest       ymm5,ymm5
   7.24 │     ↓ jne          154
   7.20 │       lea          r8,[rsi+0x40]
   7.08 │       add          rsi,0xa0
   4.72 │       cmp          rsi,rax
   3.19 │       mov          rsi,r8
   3.13 │     ↑ jbe          40

Unrolling four times gives the following kernel:

   0.79 │ 40:   vmovdqu      ymm4,YMMWORD PTR [rdi+rsi*1+0x20]
   0.53 │       vmovdqu      ymm5,YMMWORD PTR [rdi+rsi*1+0x40]
   0.50 │       vmovdqu      ymm6,YMMWORD PTR [rdi+rsi*1+0x60]
   0.52 │       vpcmpeqb     ymm7,ymm0,ymm2
   1.73 │       vperm2i128   ymm3,ymm2,ymm4,0x21
   1.86 │       vpalignr     ymm3,ymm3,ymm2,0x1
   2.43 │       vmovdqu      ymm2,YMMWORD PTR [rdi+rsi*1+0x80]
   2.29 │       vpcmpeqb     ymm3,ymm3,ymm1
   1.95 │       vpand        ymm3,ymm7,ymm3
   2.29 │       vperm2i128   ymm7,ymm4,ymm5,0x21
   1.10 │       vpalignr     ymm7,ymm7,ymm4,0x1
   0.94 │       vpcmpeqb     ymm4,ymm0,ymm4
   0.84 │       vpcmpeqb     ymm7,ymm7,ymm1
   0.69 │       vpand        ymm4,ymm4,ymm7
   0.71 │       vperm2i128   ymm7,ymm5,ymm6,0x21
   1.04 │       vpor         ymm8,ymm3,ymm4
   0.76 │       vpalignr     ymm7,ymm7,ymm5,0x1
   1.06 │       vpcmpeqb     ymm5,ymm0,ymm5
   0.72 │       vpcmpeqb     ymm7,ymm7,ymm1
   1.04 │       vpand        ymm5,ymm5,ymm7
   0.93 │       vperm2i128   ymm7,ymm6,ymm2,0x21
   0.98 │       vpalignr     ymm7,ymm7,ymm6,0x1
   1.21 │       vpcmpeqb     ymm6,ymm0,ymm6
   1.51 │       vpcmpeqb     ymm7,ymm7,ymm1
   1.23 │       vpand        ymm6,ymm6,ymm7
   3.14 │       vpor         ymm7,ymm5,ymm6
   1.50 │       vpor         ymm7,ymm8,ymm7
  19.06 │       vptest       ymm7,ymm7
   9.21 │     ↓ jne          140
   8.78 │       lea          r8,[rsi+0x80]
   8.29 │       add          rsi,0x120
   8.92 │       cmp          rsi,rax
   0.56 │       mov          rsi,r8
   0.62 │     ↑ jbe          40

Which also stalls on vtest, albeit less frequently.

1.2.4. AVX512 (AVX-512F + AVX-512BW + AVX-512VBMI) implementation. ⎇

Straightforward translation using _mm512_permutex2var_epi8 yields us:

NOINLINE size_t memchr2_avx512vbmi(const void * ptr, size_t n, uint16_t needle) {
  const uint8_t * p = (const uint8_t *) ptr;
  if (n < 2) return n;
  const uint8_t lo = (uint8_t)needle;
  const uint8_t hi = (uint8_t)(needle >> 8);
  const __m512i vlo = _mm512_set1_epi8((char)lo);
  const __m512i vhi = _mm512_set1_epi8((char)hi);
  const __m512i idx = _mm512_setr_epi32(
      0x04030201, 0x08070605, 0x0c0b0a09, 0x100f0e0d,
      0x14131211, 0x18171615, 0x1c1b1a19, 0x201f1e1d,
      0x24232221, 0x28272625, 0x2c2b2a29, 0x302f2e2d,
      0x34333231, 0x38373635, 0x3c3b3a39, 0x403f3e3d
  );
  size_t i = 0;
  if (n >= 128) {
    __m512i A = _mm512_loadu_si512((const void *)(p + 0));
    for (; i + 128 <= n; i += 64) {
      __m512i N = _mm512_loadu_si512((const void *)(p + i + 64));
      __m512i B = _mm512_permutex2var_epi8(A, idx, N);
      __mmask64 k0 = _mm512_cmpeq_epi8_mask(A, vlo);
      __mmask64 k1 = _mm512_cmpeq_epi8_mask(B, vhi);
      __mmask64 k  = k0 & k1;
      if (k) return i + (size_t)_tzcnt_u64((uint64_t)k);
      A = N;  // roll
    }
  }
  for (; i + 1 < n; ++i) if (p[i] == lo && p[i + 1] == hi) return i;
  return n;
}