27 April 2026 · 8 min
We sometimes have to look for a value in a sorted array. The simplest algorithm consists in just going through the values one by one, until we encounter the value, or exhaust the array. We sometimes call this algorithm a linear search. In C++, you can get the desired effect with the std::find function.
For large arrays, you can do better with a binary search. Binary search is a classic algorithm that efficiently locates a target value in a sorted array by repeatedly dividing the search interval in half. Starting with the entire array, it compares the target to the middle element: if the target is smaller, it discards the upper half; if larger, it discards the lower half. This process continues until the target is found or the interval is empty. It is much faster than linear search for large datasets. In C++, this is implemented by the std::binary_search function, which returns a boolean indicating whether the value is present.
The popular Roaring Bitmap format uses arrays of 16-bit integers of size ranging from 1 to 4096. We sometimes have to check whether a value is present. We use a binary search.
I wanted a faster approach. I had two insights.
- Virtually all processors today have data parallel instructions (sometimes called SIMD) that can check several values at once. Both 64-bit ARM and x64 processors (Intel/AMD) always support comparing eight 16-bit integers with a target value using a single instruction. This suggests that you should not bother going down in the binary search to blocks that are smaller than eight elements. And you may also want to cheaply compare sixteen elements or more.
- The binary search checks one value at a time. However, recent processors can load and check more than one value at once. They have excellent memory-level parllelism. This suggest that instead of a binary search, we might want to try a quaternary search: instead of splitting arrays in halves, we might split them in quarters. The net result might generate a few more instructions but the number of instructions is likely not the limiting factor.
Thus, I created something I call the SIMD Quad algorithm. It is an efficient search algorithm for sorted arrays of 16-bit unsigned integers, combining a quaternary interpolation search with SIMD (Single Instruction, Multiple Data). The algorithm divides the array into fixed-size blocks of 16 elements (except maybe for the last block) and uses the last element of each block as interpolation keys to quickly narrow down the search to a single block, then employs SIMD instructions to check all 16 elements in that block simultaneously.
The core idea is to perform a hierarchical search: first, use interpolation search on a coarser level (block boundaries) to find the likely block containing the target value, then switch to SIMD for fine-grained parallel checking within the block. This hybrid approach leverages the strengths of both algorithmic optimization (interpolation search reduces comparisons logarithmically) and hardware acceleration (SIMD checks multiple elements at once).
- Initial Check: If the array has fewer than 16 elements, perform a simple linear search through all elements.
- Block Division: Divide the array into blocks of 16 consecutive elements. For an array of size
cardinality, there arenum_blocks = cardinality / 16full blocks. - Quaternary Interpolation Search: Use the last element of each block (at positions
16-1,32-1, etc.) as keys for interpolation. The search performs a quaternary (base-4) interpolation to find the block where the targetposis likely located. This involves comparing the target against quarter-points of the current search range and adjusting the base accordingly. - Block Selection: After narrowing down, select the appropriate block index
lobased on the interpolation results. - SIMD Check: If a valid block is found, load the 16 elements into SIMD registers (using NEON on ARM or SSE2 on x64) and perform parallel equality comparisons with the target value. If any match is found, return true.
- Remainder Check: For any elements not in full blocks (remainder), perform a linear search.
How does it do? I wrote a benchmark. The benchmark works as follows. For each array size from 2 to 4096 elements, it generates 100,000 sorted arrays of 16-bit unsigned integers. For each size, it performs 10 million membership queries in “cold” mode (each query searches a different array, simulating cache misses) and 10 million queries in “warm” mode (queries are grouped by array, with each array being searched 100 times consecutively, simulating cache hits). The benchmark measures the average time per query for three algorithms: linear search (std::find), binary search (std::binary_search), and the new SIMD Quad algorithm.
I use two systems. An Apple M4 with Apple LLVM and an Intel Emeral Rapids processor with GCC.
Firstly, let us compare the linear search with the binary search.
Intel/GCC:
Apple/LLVM
The result is clear. The binary search beats the linear search as soon as the arrays get large. That is to be expected.
On a cold cache, the linear search is relatively worse. That is to be expected because it accesses more data, causing more cache faults.
We have established that the binary search is the net winner over the linear search. Let us now compare with the SIMD Quad algorithm.
Intel/GCC:
Apple/LLVM
The results differ markedly between the Intel and Apple platform. On the Intel platform the SIMD Quad is more than twice as fast as the binary search on the warm cache. The benefits are lesser on the cold cache. On the Apple platform, the reverse is true, it is with the cold cache that the SIMD Quad is more than twice as fast, whereas the benefits are more marginal on the warm cache.
But the important point is that, in all instances, SIMD Quad is faster than the binary search.
The SIMD component of the algorithm is rather straightforward: we use specialized instructions that save work. So it is easy to see why it might make things faster. There are few instructions, fewer branches.
But what about the ‘quad’ part. Does it matter? So I tried a binary version of the same algorithm. It has the same SIMD optimization, but I am dropping the quaternary interpolation search and replacing it with a standard binary search.
Intel/GCC:
Apple/LLVM
To put it in simple terms, the quad approach has little effect on the Apple platform, but it is a decent optimization on the Intel platform for large arrays in the cold case. The quaternary search better exploits the memory-level parallelism on my Intel server.
Conclusion. What my results suggest is that while a textbook binary search is a decent algorithm, you can do better in ways that matter. Standard algorithms were often not designed for computers that have so much parallelism. The SIMD Quad algorithm tries to leverage both the memory-level and data parallelism. Further, I suspect that we can do even better than my algorithm. Let us get creative!
Further reading: Faster intersections between sorted arrays with shotgun
Appendix (source code)
bool simd_quad(const uint16_t *carr, int32_t cardinality,
uint16_t pos) {
constexpr int32_t gap = 16;
if (cardinality < gap) {
for (int32_t j = 0; j < cardinality; j++) {
if (carr[j] == pos) return true;
}
return false;
}
int32_t num_blocks = cardinality / gap;
int32_t base = 0;
int32_t n = num_blocks;
while (n > 3) {
int32_t quarter = n >> 2;
int32_t k1 = carr[(base + quarter + 1) * gap - 1];
int32_t k2 = carr[(base + 2 * quarter + 1) * gap - 1];
int32_t k3 = carr[(base + 3 * quarter + 1) * gap - 1];
int32_t c1 = (k1 < pos);
int32_t c2 = (k2 < pos);
int32_t c3 = (k3 < pos);
base += (c1 + c2 + c3) * quarter;
n -= 3 * quarter;
}
while (n > 1) {
int32_t half = n >> 1;
base = (carr[(base + half + 1) * gap - 1] < pos)
? base + half : base;
n -= half;
}
int32_t lo = (carr[(base + 1) * gap - 1] < pos)
? base + 1 : base;
if (lo < num_blocks) {
const uint16_t *blk = carr + lo * gap;
#ifdef __ARM_NEON
uint16x8_t needle = vdupq_n_u16(pos);
uint16x8_t v0 = vld1q_u16(blk);
uint16x8_t v1 = vld1q_u16(blk + 8);
uint16x8_t hit = vorrq_u16(vceqq_u16(v0, needle),
vceqq_u16(v1, needle));
return vmaxvq_u16(hit) != 0;
#else
__m128i needle = _mm_set1_epi16((short)pos);
__m128i v0 = _mm_loadu_si128((const __m128i *)blk);
__m128i v1 = _mm_loadu_si128((const __m128i *)(blk + 8));
__m128i hit = _mm_or_si128(_mm_cmpeq_epi16(v0, needle),
_mm_cmpeq_epi16(v1, needle));
return _mm_movemask_epi8(hit) != 0;
#endif
}
for (int32_t j = num_blocks * gap; j < cardinality; j++) {
uint16_t v = carr[j];
if (v >= pos) return (v == pos);
}
return false;
}
Daniel Lemire, "You can beat the binary search," in Daniel Lemire's blog, April 27, 2026, https://lemire.me/blog/2026/04/27/you-can-beat-the-binary-search/.
[BibTeX]
> The results differ markedly between the Intel and AMD platform
Is that supposed to be ARM instead of AMD?
Probably worth adding that you can do better with other data structures if you’re not required to operate on already sorted data.
Something like Eytzinger binary search (see https://algorithmica.org/en/eytzinger) provides a more cache-efficient layout for the same problem compared to binary search. In the SIMD case (as opposed to binary), this layout is similar to B trees, effectively allowing you to read k1/k2/k3 from consecutive locations in memory as opposed to gathering them from different addresses.
I can’t say much about the interpolation search part; it’s not always applicable and has unpredictable latency on non-random data, and for known-random data you might as well use a hash table, so I never used it.
I’d not heard of Eytzinger before, they look a lot like the more well known “heap” data structure (as opposed to heap memory) which tend to be associated with priority queues but variants like B-heaps exhist that seem almost identical.
I wonder if something like a B+tree data structure would be a better fit for SIMD access than either. These are traditionally designed for storage with high latency (e.g. HDDs) but would seem to translate well to fetching a cache line from RAM (instead of a page from disk). It would seem possible to do a compare across the full SIMD register width every time, resulting in an approximately constant factor speedup. I.e. you’d be doing log_16 compares rather than log_2.
I think the demonstration would have been more thorough if the change from branchy to branchless was included (perhaps instead of the linear vs binary search). The remainder of the article uses the branchy std::binary_search as a baseline for comparison to branchless SIMD implementations, which is a poor representation of the performance difference between the scalar and SIMD algorithms.
I added the “lower_bound_overlap” function from Orson Peter’s blog (https://orlp.net/blog/bitwise-binary-search/), into your simple.cpp benchmark (my results below.) Interestingly, it competes with the SIMD functions array size 256 (cold cache) and at array size 4096 (warm cache). While the resulting graph tells a very similar story, the narrower performance difference is also notable.
It would have been interesting to see how a better baseline performs compared to your SIMD functions on Intel.
—– Results on Apple Macbook Air M1, lengths: 64, 256, 1024, 4096 —-
arrays: 100000 size per array: 64 warmth: 100 total queries per mode: 10000000
verification passed: 10000 queries × all algorithms agree with std::binary_search
simd_quad cold : 4.251 ns 0.24 Gv/s
simd_quad warm : 3.730 ns 0.27 Gv/s
simd_binary cold : 4.151 ns 0.24 Gv/s
simd_binary warm : 3.225 ns 0.31 Gv/s
orlp_lower_bound cold : 5.948 ns 0.17 Gv/s
orlp_lower_bound warm : 5.655 ns 0.18 Gv/s
binary_search cold : 8.074 ns 0.12 Gv/s
binary_search warm : 7.843 ns 0.13 Gv/s
arrays: 100000 size per array: 256 warmth: 100 total queries per mode: 10000000
verification passed: 10000 queries × all algorithms agree with std::binary_search
simd_quad cold : 9.251 ns 0.11 Gv/s
simd_quad warm : 5.650 ns 0.18 Gv/s
simd_binary cold : 13.970 ns 0.07 Gv/s
simd_binary warm : 4.780 ns 0.21 Gv/s
orlp_lower_bound cold : 10.724 ns 0.09 Gv/s
orlp_lower_bound warm : 7.064 ns 0.14 Gv/s
binary_search cold : 21.651 ns 0.05 Gv/s
binary_search warm : 10.081 ns 0.10 Gv/s
arrays: 100000 size per array: 1024 warmth: 100 total queries per mode: 10000000
verification passed: 10000 queries × all algorithms agree with std::binary_search
simd_quad cold : 53.188 ns 0.02 Gv/s
simd_quad warm : 8.297 ns 0.12 Gv/s
simd_binary cold : 48.651 ns 0.02 Gv/s
simd_binary warm : 7.063 ns 0.14 Gv/s
orlp_lower_bound cold : 67.511 ns 0.01 Gv/s
orlp_lower_bound warm : 8.851 ns 0.11 Gv/s
binary_search cold : 105.724 ns 0.01 Gv/s
binary_search warm : 13.250 ns 0.08 Gv/s
arrays: 100000 size per array: 4096 warmth: 100 total queries per mode: 10000000
verification passed: 10000 queries × all algorithms agree with std::binary_search
simd_quad cold : 82.593 ns 0.01 Gv/s
simd_quad warm : 13.056 ns 0.08 Gv/s
simd_binary cold : 97.889 ns 0.01 Gv/s
simd_binary warm : 11.552 ns 0.09 Gv/s
orlp_lower_bound cold : 138.318 ns 0.01 Gv/s
orlp_lower_bound warm : 11.357 ns 0.09 Gv/s
binary_search cold : 210.640 ns 0.00 Gv/s
binary_search warm : 17.811 ns 0.06 Gv/s
I’ve re-run the simple benchmark on a i5-6500 CPU with the aforementioned branchless scalar binary search (results below.)
The simd_quad algorithm eventually pulls away for the cold cache case, but the improvements are not as definitive as in the article.
Additionally, it would be trivial to add the optimization to bitwise binary search to use SIMD once the block is the desired power of two. However, I used the generic version in keeping with the spirit of setting a scalar baseline.
——— Results i5-6500, lengths 64, 128, 1024, 2048, 4096 ———
arrays: 100000 size per array: 64 warmth: 100 total queries per mode: 10000000
verification passed: 10000 queries × all algorithms agree with std::binary_search
simd_quad cold : 12.803 ns 0.08 Gv/s 3.06 GHz 39.64 c 75.85 i 0.02 bm 1.91 i/c
simd_quad warm : 7.633 ns 0.13 Gv/s 3.59 GHz 27.45 c 75.85 i 0.02 bm 2.76 i/c
simd_binary cold : 13.847 ns 0.07 Gv/s 2.99 GHz 41.41 c 70.85 i 0.02 bm 1.71 i/c
simd_binary warm : 7.278 ns 0.14 Gv/s 3.59 GHz 26.16 c 70.85 i 0.02 bm 2.71 i/c
orlp_lower_bound cold : 13.032 ns 0.08 Gv/s 3.44 GHz 46.53 c 76.03 i 0.02 bm 1.63 i/c
orlp_lower_bound warm : 8.987 ns 0.11 Gv/s 3.59 GHz 32.31 c 76.03 i 0.02 bm 2.35 i/c
binary_search cold : 35.507 ns 0.03 Gv/s 3.59 GHz 127.47 c 89.28 i 3.32 bm 0.70 i/c
binary_search warm : 29.919 ns 0.03 Gv/s 3.59 GHz 107.52 c 89.28 i 3.17 bm 0.83 i/c
arrays: 100000 size per array: 128 warmth: 100 total queries per mode: 10000000
verification passed: 10000 queries × all algorithms agree with std::binary_search
simd_quad cold : 20.492 ns 0.05 Gv/s 2.79 GHz 57.26 c 93.92 i 0.01 bm 1.64 i/c
simd_quad warm : 9.468 ns 0.11 Gv/s 3.59 GHz 34.04 c 93.92 i 0.01 bm 2.76 i/c
simd_binary cold : 26.221 ns 0.04 Gv/s 2.69 GHz 70.59 c 81.92 i 0.01 bm 1.16 i/c
simd_binary warm : 8.961 ns 0.11 Gv/s 3.59 GHz 32.22 c 81.92 i 0.01 bm 2.54 i/c
orlp_lower_bound cold : 23.303 ns 0.04 Gv/s 2.79 GHz 65.06 c 82.02 i 0.01 bm 1.26 i/c
orlp_lower_bound warm : 9.615 ns 0.10 Gv/s 3.59 GHz 34.56 c 82.02 i 0.01 bm 2.37 i/c
binary_search cold : 46.073 ns 0.02 Gv/s 3.59 GHz 165.39 c 100.64 i 4.18 bm 0.61 i/c
binary_search warm : 36.156 ns 0.03 Gv/s 3.59 GHz 129.91 c 100.64 i 4.07 bm 0.77 i/c
arrays: 100000 size per array: 256 warmth: 100 total queries per mode: 10000000
verification passed: 10000 queries × all algorithms agree with std::binary_search
simd_quad cold : 49.535 ns 0.02 Gv/s 2.29 GHz 113.57 c 101.96 i 0.00 bm 0.90 i/c
simd_quad warm : 11.767 ns 0.08 Gv/s 3.59 GHz 42.30 c 101.96 i 0.00 bm 2.41 i/c
simd_binary cold : 59.820 ns 0.02 Gv/s 2.09 GHz 125.21 c 92.96 i 0.00 bm 0.74 i/c
simd_binary warm : 10.620 ns 0.09 Gv/s 3.59 GHz 38.18 c 92.96 i 0.00 bm 2.43 i/c
orlp_lower_bound cold : 44.852 ns 0.02 Gv/s 2.59 GHz 116.29 c 88.01 i 0.00 bm 0.76 i/c
orlp_lower_bound warm : 10.671 ns 0.09 Gv/s 3.59 GHz 38.35 c 88.01 i 0.00 bm 2.29 i/c
binary_search cold : 64.725 ns 0.02 Gv/s 3.29 GHz 213.25 c 112.07 i 4.95 bm 0.53 i/c
binary_search warm : 42.091 ns 0.02 Gv/s 3.59 GHz 151.27 c 112.07 i 4.90 bm 0.74 i/c
arrays: 100000 size per array: 1024 warmth: 100 total queries per mode: 10000000
verification passed: 10000 queries × all algorithms agree with std::binary_search
simd_quad cold : 146.389 ns 0.01 Gv/s 1.59 GHz 232.96 c 127.99 i 0.00 bm 0.55 i/c
simd_quad warm : 17.109 ns 0.06 Gv/s 3.59 GHz 61.50 c 127.99 i 0.00 bm 2.08 i/c
simd_binary cold : 182.801 ns 0.01 Gv/s 1.44 GHz 272.68 c 114.99 i 0.00 bm 0.42 i/c
simd_binary warm : 17.090 ns 0.06 Gv/s 3.59 GHz 61.38 c 114.99 i 0.00 bm 1.87 i/c
orlp_lower_bound cold : 147.635 ns 0.01 Gv/s 1.50 GHz 222.43 c 100.00 i 0.00 bm 0.45 i/c
orlp_lower_bound warm : 15.366 ns 0.07 Gv/s 3.59 GHz 55.22 c 100.00 i 0.00 bm 1.81 i/c
binary_search cold : 143.199 ns 0.01 Gv/s 2.62 GHz 385.50 c 135.02 i 6.05 bm 0.35 i/c
binary_search warm : 52.658 ns 0.02 Gv/s 3.58 GHz 189.22 c 135.02 i 6.09 bm 0.71 i/c
arrays: 100000 size per array: 2048 warmth: 100 total queries per mode: 10000000
verification passed: 10000 queries × all algorithms agree with std::binary_search
simd_quad cold : 178.018 ns 0.01 Gv/s 1.59 GHz 284.69 c 146.00 i 0.00 bm 0.51 i/c
simd_quad warm : 23.755 ns 0.04 Gv/s 3.30 GHz 78.33 c 146.00 i 0.00 bm 1.86 i/c
simd_binary cold : 261.060 ns 0.00 Gv/s 1.29 GHz 337.12 c 126.00 i 0.00 bm 0.37 i/c
simd_binary warm : 26.227 ns 0.04 Gv/s 3.13 GHz 82.26 c 126.00 i 0.00 bm 1.53 i/c
orlp_lower_bound cold : 220.550 ns 0.00 Gv/s 1.22 GHz 283.33 c 106.00 i 0.00 bm 0.37 i/c
orlp_lower_bound warm : 23.861 ns 0.04 Gv/s 3.00 GHz 71.51 c 106.00 i 0.00 bm 1.48 i/c
binary_search cold : 181.147 ns 0.01 Gv/s 2.59 GHz 470.60 c 146.51 i 6.51 bm 0.31 i/c
binary_search warm : 58.688 ns 0.02 Gv/s 3.59 GHz 210.93 c 146.51 i 6.57 bm 0.69 i/c
arrays: 100000 size per array: 4096 warmth: 100 total queries per mode: 10000000
verification passed: 10000 queries × all algorithms agree with std::binary_search
simd_quad cold : 264.411 ns 0.00 Gv/s 1.39 GHz 370.33 c 154.00 i 0.00 bm 0.42 i/c
simd_quad warm : 35.036 ns 0.03 Gv/s 2.85 GHz 100.08 c 154.00 i 0.00 bm 1.54 i/c
simd_binary cold : 340.375 ns 0.00 Gv/s 1.19 GHz 405.65 c 137.00 i 0.00 bm 0.34 i/c
simd_binary warm : 40.746 ns 0.02 Gv/s 2.70 GHz 109.87 c 137.00 i 0.00 bm 1.25 i/c
orlp_lower_bound cold : 312.105 ns 0.00 Gv/s 1.09 GHz 340.97 c 112.00 i 0.00 bm 0.33 i/c
orlp_lower_bound warm : 35.670 ns 0.03 Gv/s 2.70 GHz 96.17 c 112.00 i 0.00 bm 1.16 i/c
binary_search cold : 214.620 ns 0.00 Gv/s 2.59 GHz 555.95 c 158.00 i 6.97 bm 0.28 i/c
binary_search warm : 66.074 ns 0.02 Gv/s 3.59 GHz 237.40 c 158.00 i 7.05 bm 0.67 i/c
Instead of the quad interpolation, would there be some benefit in loading the last element of each block in SIMD register to compare several keys at the same time ?
Why not use a bitmask of which values are present? 16-bit integers compress to an 8kB bitmask. What am I missing?
The Roaring Bitmap format does adopt a bitmap when there are more than 4096 values. Please see https://github.com/RoaringBitmap/RoaringFormatSpec
I recall reading many years ago that if binary search was made interpolatory, so that given two values that surround the desired value interpolation was used to get the next place to look, the time improved from log(n) to, on average, log(log(n)).
I could be wrong.
Under some conditions, this result is true.
I tried your simd_quad with gap = 64 and two paths for cardinality < 8 and < gap and got faster results on zen3. The ending search loop might need tuning
https://gist.github.com/aconz2/5a1bd36ce34afeeb5479f6ce775bb3d7
Thanks Andrew. I am not surprised!!!
I’ve coded up the following optimization ideas:
1. Eagerly and pessimistically checking the (n % gap) values at the end of the array using SWAR/SIMD. This removes the loop at the end of the function and fixes poor performance with n = 2 * gap, do block-wise binary search on the range 0..(n-gap), then do SIMD on the result of binary search unconditionally.
Results (Apple M1), warm cache:
Optimization #1 results 0.9x to 4x performance (n=29 is atrocious for simd_binary/simd_quad).
Original simd_binary is marginally faster between n=32 to ~1000.
At n ~=1000, my changes begin edging out simd_binary.
At n >= 2000, my changes consistently beat simd_binary, (e.g up to 9.5% faster for troublesome n)
At n >= 3000, up to 12% faster for troublesome n
At n >= 4000, up to 25% faster for troublesome n
Perhaps, this is considered gaming the benchmark since I’ve made all branches solely dependant on the array length, allowing the branch predictor to reach extreme accuracy. Never-the-less, it’s can still provide ideas for optimizations.
Code (made for idea testing, not production):
https://gist.github.com/wald0047/b26f2a1844a0224f96b9b546776cf1c0#file-simple-cpp-L71
The top of my comment was modified by HTML sanitization, fixed version:
1. Eagerly and pessimistically checking the (n % gap) values at the end of the array using SWAR/SIMD. This removes the loop at the end of the function and fixes poor performance with n ≤ 31.
2. When n ≥ 2 * gap, do block-wise binary search on the range 0..(n-gap), then do SIMD on the result of binary search unconditionally.
Do you have a restriction on which SIMD level you use? AVX2 Gather might be useful for the pivot comparison loop, and it goes up to 8-way loads.
It’s fast on recent processors, even if originally it was slow.
In this particular instance, I am using 16-bit unsigned values. I think that gather works with 32-bit or 64-bit integers. In my case, I retrieve three 16-bit values.
I am not hopeful that it will help in this case.
Although, you are right that it is well worth keeping in mind, more generally.
Yes, you get two elements instead of one :).
I’m not sure it’s accurate to describe this as “interpolation search”.
[Interpolation Search](https://en.wikipedia.org/wiki/Interpolation_search) uses the values checked to estimate where next to look, but this algorithm simply checks 3 locations (determined only by size of the current window, not the values), and simply narrows the search to whichever quarter can contain the desired chunk.
I agree with your criticism.
The SIMD Quad algorithm’s performance benefits seem to vary significantly between Intel and Apple platforms, especially concerning cold vs. warm cache scenarios. It’s interesting how the ‘quad’ aspect of the search has a more pronounced effect on Intel, while SIMD alone is more impactful on Apple. Have you explored if there are specific architectural differences in memory access patterns or instruction pipeline efficiencies that contribute to this divergence?
OmniaKey