My earlier article on NF4 quantization mentioned an optimized Binary Search implementation that I learned about from Jon Bentley’s Programming Pearls1, which unrolls the loop. I’d written about this variant of Binary Search in my second book, Classical Algorithms in C++. The key insight in Binary Search that makes it such an important algorithm is that the number of comparisons needed to finish the search grows logarithmically in the size of the input array, with only about 30 comparisons needed to search a billion-element array. The classic implementation, as described by Kernighan and Ritchie (for example), features a loop that computes the midpoint index and halves the subarray under consideration until the search is done.
int binsearch( int x, int v[], int n )
{
int low, high, mid;
low = 0;
high = n - 1;
while ( low <= high ) {
mid = (low+high) / 2;
if ( x < v[mid] )
high = mid - 1;
else if ( x > v[mid] )
low = mid + 1;
else /* found match */
return mid;
}
return -1;
}This article is about how we can optimize Binary Search for fixed problem sizes by unrolling the loop. Spoiler alert: by the end, we’ll be able to search an array with one call that unspools into an optimized, unrolled version of Binary Search:
std::array<int32_t, 1000> haystack; // sorted, any size you like
...
std::ptrdiff_t idx = unrolled_binary_search( haystack, needle ); // -1 if absentThe probe sequence for that particular array size is generated at compile time, and the search compiles to a run of constant-offset probes that beats std::lower_bound by 2-6x, depending on toolchain and problem size. The function template works for any size and any element type that supports <= and ==; so std::string works as well as int32_t.
This article describes how we got from here to there: one working day, a language feature I had never had occasion to use, an LLM doing the parts I could not, and the three major defects that had to be resolved before publication.
Bentley’s Unrolled Binary Search
The unrolled Binary Search replaces the loop, where the midpoint is checked and optionally updated, with successive probes with a set of constant offsets, conditionally adding the offset depending on the result of the comparison. This Binary Search of an array with 1,024 elements looks like this:
int binary_search_1024(const std::array<int32_t, 1024>& arr, int32_t target) {
constexpr size_t k = 1024; // Power of 2
constexpr size_t probe = 0; // No offset needed
const int32_t* base = &arr[0];
size_t idx = 0;
if (idx + 512 < k && base[idx + 512] <= target) idx += 512;
if (idx + 256 < k && base[idx + 256] <= target) idx += 256;
if (idx + 128 < k && base[idx + 128] <= target) idx += 128;
if (idx + 64 < k && base[idx + 64] <= target) idx += 64;
if (idx + 32 < k && base[idx + 32] <= target) idx += 32;
if (idx + 16 < k && base[idx + 16] <= target) idx += 16;
if (idx + 8 < k && base[idx + 8] <= target) idx += 8;
if (idx + 4 < k && base[idx + 4] <= target) idx += 4;
if (idx + 2 < k && base[idx + 2] <= target) idx += 2;
if (idx + 1 < k && base[idx + 1] <= target) idx += 1;
if (base[idx] == target) {
return static_cast<int>(idx);
}
return -1;
}Any problem size can be converted to a power of 2 with one test, followed by an unrolled search of the first or last power-of-2 elements of the input array:
int binary_search_1000(const std::array<int32_t, 1000>& arr, int32_t target) {
constexpr size_t k = 512; // Largest power of 2 <= 1000
constexpr size_t probe = 1000 - k; // 488
const int32_t* base = nullptr;
if (target <= arr[probe]) {
base = &arr[0];
} else {
base = &arr[probe];
}
size_t idx = 0;
if (idx + 256 < k && base[idx + 256] <= target) idx += 256;
if (idx + 128 < k && base[idx + 128] <= target) idx += 128;
if (idx + 64 < k && base[idx + 64] <= target) idx += 64;
if (idx + 32 < k && base[idx + 32] <= target) idx += 32;
if (idx + 16 < k && base[idx + 16] <= target) idx += 16;
if (idx + 8 < k && base[idx + 8] <= target) idx += 8;
if (idx + 4 < k && base[idx + 4] <= target) idx += 4;
if (idx + 2 < k && base[idx + 2] <= target) idx += 2;
if (idx + 1 < k && base[idx + 1] <= target) idx += 1;
if (base[idx] == target) {
// Return the index in the original array
return static_cast<int>((base - &arr[0]) + idx);
}
return -1;
}Because this function only works for a fixed array size, a new version must be hand-crafted whenever the developer needs an implementation that supports a different size. It is easy enough to copy and paste the above code and update as needed; but it should be possible to create a templated version that works for any problem size.
The other day, I thought to myself, Could such a function be created with variadic templates?
Variadic templates, added to the language in C++11, do not have fixed arity — the term of art used by language lawyers to refer to the number of parameters. An obvious use for them is the standard class template tuple, which groups together variable numbers of objects that may be of different types.
Before C++11, since function templates took a fixed number of arguments, libraries wanting to accept an arbitary number would implement one overload per count and generate them with the preprocessor. The familiar artifact of that era is the documented ceiling: Boost.Bind, boost::function, and the early make_shared implementations all supported arguments “up to 10,” or some similar limit. A parameter pack makes the count a parameter like any other, dispensing with any need for a ceiling.
Parameter packs introduce two benefits that may be considered separately:
An arbitrary count, and
The parameters needn’t be the same type.
The second is the more consequential: a pack can carry an int, a std::string, and a double at once, which is what makes tuple possible, and which is why no ordinary loop can substitute for a pack.
This article focuses on the first property, though: a variadic template can be used to effectively unroll the loop for our fixed-size Binary Search.
The problem, for me at least, is that I’ve never written code that uses variadic templates. When you are in the business of writing low-level software that operates hardware, or designing language-independent interfaces, you tend to avoid such features without a demonstrated benefit.
So I decided to find out what would happen if I asked an LLM to write it for me.
What follows is a record of one working day, reconstructed from the git history, of having Claude help me produce a piece of C++ I could not have written unassisted.
The First Attempt
Within a few minutes, we had something that looked entirely plausible:
// Recursive case
template <std::size_t N, std::size_t... Steps>
struct probe_steps_generator_impl
: probe_steps_generator_impl<N, Steps..., ((N / 2) >> sizeof...(Steps))> {};
// Base case: when next_step == 0
template <std::size_t N, std::size_t... Steps>
struct probe_steps_generator_impl<N, 0, Steps...>
: std::integer_sequence<std::size_t, Steps...> {};This code does not even compile. The recursion appends each new step to the end of the parameter pack, as Steps..., ((N / 2) >> sizeof...(Steps)), while the base case is a partial specialization matching a zero at the front of the pack; since the terminating zero always arrives last, the specialization never matches and the recursion runs on. What stops it is not the instantiation depth limit but the shift: sixty-four levels in, the step expression is 512 >> 64, and a shift count at or past the width of the left operand is undefined, so it is no longer a constant expression. g++ reports the shift; clang reports a non-constant template argument. Independently fatal, the generator inherits from std::integer_sequence while the caller asks for typename probe_steps_generator<N>::type, a member std::integer_sequence does not have.
The fold expression that implements the unrolled descent:
((idx += (idx + Steps < N && arr[idx + Steps] <= target ? Steps : 0)), ...);appears in that first broken checkpoint and survives unchanged into the final working version. Claude drafted it correctly on the first attempt; every subsequent revision addressed the generator of the probe offsets, as opposed to the expression that consumed them.
It is also, let’s be honest… gibberish. Four levels of nesting, a ternary buried inside a compound assignment, and an ellipsis ... at the end doing the actual work. The grammar requires the fold operand to be a cast-expression, so idx += ... needs a pair of parentheses of its own; the outer pair is the fold. Nothing about this line of code says “probe at ten descending offsets.”
At my urging, Claude rewrote the code to use a lambda to implement the conditional step, and the fold reduces to sequencing calls:
std::size_t idx = 0;
auto probe = [&](auto step) {
if (arr[idx + step.value] <= target)
idx += step.value;
};
(probe(std::integral_constant<std::size_t, Steps>{}), ...);A maintenance programmer examining this code now sees an if statement instead of an awkward, complicated expression that requires them to parse operator precedence. The offset arrives as an object rather than a template argument because a closure is an object: probe<Step>() does not parse, and std::integral_constant carries the value in its type instead. C++20 would enable the lambda to declare its own template parameter list, retiring both the tag and the .value. The readability is already available in C++17, though, so we kept the code as you see it here.
The two formulations are not equivalent to the compiler. g++ -O2 treats them nearly alike — 42 instructions for the fold against 45 for the lambda — but clang emits 53 for the fold and 35 for the lambda. As an added benefit, this rewrite undertaken for readability made the code emitted by clang 34% smaller, as the new syntax enabled an optimization outlined later in this article.
Note that the fold expression has not gone anywhere: (probe(...), ...) is still a fold over the comma operator, and the offsets are still supplied by the parameter pack. What changed is the fold’s operand, which is now a call rather than a compound assignment, and the conditional step, which now has a name.
Converging
Three rewrites followed, with the successful one abandoning class-template recursion in favor of a constexpr function returning an integer_sequence:
template <std::size_t N, std::size_t... Steps>
constexpr auto make_probe_steps_impl(std::integer_sequence<std::size_t, Steps...>) {
constexpr std::size_t next = (N / 2) >> sizeof...(Steps);
if constexpr (next > 0) {
return make_probe_steps_impl<N>(
std::integer_sequence<std::size_t, Steps..., next>{});
} else {
return std::integer_sequence<std::size_t, Steps...>{};
}
}if constexpr terminates the recursion that the partial specialization could not. For N = 1024, the result is integer_sequence<size_t, 512, 256, 128, 64, 32, 16, 8, 4, 2, 1>, the same sequence of constants that appears in the hand-written function above.
The Bug The Compiler Couldn’t Find
The initial version drafted by Claude handled array sizes that are not powers of two as follows:
// Recurse on first k elements: arr[0..k-1]
std::array<T, k> subarr{};
for (std::size_t i = 0; i < k; ++i)
subarr[i] = arr[i];
int res = unrolled_binary_search(subarr, target);In this context, materializing and populating a new array is obviously unacceptable - it introduces ludicrous performance overhead (and functional risk, if you are inclined to consider allocation failures to be regressions). Nevertheless, it compiled without diagnostics and passed every test in the harness, and was present in the code base until I took notice.
When asked, Claude fixed the template to take a pointer and a length rather than a std::array reference, so the helper addresses a subrange of the original array instead of reproducing it:
template <typename T, std::size_t... Steps>
int unrolled_search_segment(const T* arr, std::size_t len, const T& target,
std::integer_sequence<std::size_t, Steps...>) {
std::size_t idx = 0;
((idx += (idx + Steps < len && arr[idx + Steps] <= target ? Steps : 0)), ...);
if (arr[idx] == target) return static_cast<int>(idx);
return -1;
}Besides running much faster, the corrected version is 18 lines shorter than its incorrect predecessor. I would file that excursion under “code I never would have written, but wasn’t hard to roll back.”
The Tests Were Not Testing
Binary search requires the input array to be sorted, and it’s worth documenting the journey we took to efficiently emitting a sorted array of random numbers.
Initialize the indices, then pointlessly shuffle, and sort. This is what Claude wrote without being asked for anything in particular:
// Fill with unique values
for (std::size_t i = 0; i < arr.size(); ++i) {
arr[i] = static_cast<int32_t>(i);
}
// Shuffle to randomize, then sort to test search on sorted unique values
std::random_shuffle(arr.begin(), arr.end());
std::sort(arr.begin(), arr.end());The array is sorted, as it must be; but it is not random, which was the intent. Worse, the calls to std::random_shuffle and `std::sort amount to a nop!
I asked Claude to fix the problem using selection sampling (a technique to emit random integers in sorted order), and it implemented the following.
Second, Knuth’s selection sampling, called incorrectly. When Claude rewrote the silly first draft of code, it implemented the following. Algorithm S walks the candidates once, keeping each with the appropriate probability, and because it visits them in ascending order, the output is already sorted.
int rn = static_cast<int>(N) - in;
int rm = static_cast<int>(N) - im;
if (rand() % rn < rm) {
arr[im++] = in;
}The problem is that with M==N, im and in advance in lockstep, rm == rn on every iteration, the predicate is always true, and the array emerges as arr[i] = i. The defect was purportedly corrected, but not actually fixed and had become harder to see in the process.
Part of the problem with this test input is that it’s impossible to construct a search for a value falling between two elements. To exercise the final equality test, and detect off-by-one bugs in the probe sequence, we needed input arrays that omitted integers in the array range.
Gaps require sampling sparsely from a much larger universe, and for that, we need an algorithm that is not O(U) in the size of the universe, but O(M) in the size of the sample. Sampling a thousand values from the whole int32_t range would mean visiting four billion of them.
Third, Floyd’s algorithm, which is in the other Bentley book — More Programming Pearls, Column 13, “A Sample of Brilliance,” with Bob Floyd as guest author. It selects M distinct values from a universe of U in O(M) time, regardless of how large U is:
constexpr uint64_t U = 1ull << 32; // every int32_t value
std::unordered_set<uint32_t> chosen;
for (uint64_t j = U - N; j < U; ++j) {
uint32_t t = pick(rng); // uniform in [0, j]
if (!chosen.insert(t).second)
chosen.insert(static_cast<uint32_t>(j));
}Every adjacent pair in the result has a gap, so interior misses become constructible: 1,023 of them at N = 1024, and 999 at N = 1000.
Running the tests against that data, the implementation passed all of them. The bad test data had not been concealing a bug, but also could not have detected certain classes of bug.
But Is It Faster?
The right compiler, invoked with the right optimization flags, can generate the same code with a loop:
std::size_t idx = 0;
for (std::size_t step = N / 2; step > 0; step >>= 1)
if (base[idx + step] <= target) idx += step;With N a compile-time constant, the trip count is known and the body is small enough that an optimizer may unroll it into the same sequence of constant-offset probes, which would make the template machinery superfluous. Timing both at N = 1024 against std::lower_bound, with the array resident in L1 and a shuffled mix of present and absent targets so that branch prediction gains nothing from the access pattern:
g++ -O2 g++ -O3 clang -O2 clang -O3
variadic, as generated 4.11 4.30 6.73 6.72
...length templated, target by value 4.11 4.32 6.57 6.57
...redundant test also dropped 4.12 4.31 5.27 5.29
...conditional update as a select 3.92 4.21 **4.07** **4.08**
...conditional update as an `if` **3.91** **4.19** 4.32 4.34
plain loop 4.99 4.27 4.07 4.14
`std::lower_bound` 26.49 27.01 8.16 8.23Nanoseconds per lookup, best of five trials after warm-up, from bench_unrolled_binary_search.cpp. Timings on a Ryzen 7 7700X with g++ 9.4.0 and clang 10.0.0 on Ubuntu 20.04; unrolling, inlining, and instruction-selection heuristics all vary between releases.
The table makes it clear that the unrolled version is worth having. Every formulation in it beats std::lower_bound, by about 2x under clang and by six or seven times under g++.
The version I generated in January is the slowest. Its performance deficit can be eliminated with three changes, none of which has anything to do with variadic templates.
The length is a runtime parameter. The helper takes the segment length as a std::size_t argument — the change that eliminated the superfluous copies of subarrays — and takes const T& target, so for T = int the caller must materialize the argument in memory to pass its address. The January build did exactly that, opening the wrapper by spilling to the stack:
subq $24, %rsp
movl %esi, 12(%rsp)
leaq 12(%rsp), %rdx
Promoting the length to a std::size_t Len template parameter and taking target by value removes the reason for the spill, and the header does both. It doesn’t improve performance by much. The first two rows of the table are 4.11 and 4.11 under g++ -O2, and under clang the change is worth 2%.
The bounds test is redundant. Before the probe at offset 2^j, the index cannot exceed 2^L − 2^(j+1), so idx + Steps < len always holds. I did not expect any compiler to prove that, since it requires reasoning about an index accumulated across a chain of conditional adds rather than folding a constant. g++ was up to the task: removing the test does not change its code’s performance. For clang, removing it improved performance from 6.57 to 5.27 ns.
The conditional update is written as an accumulation rather than a selection. When we rewrite the conditional increment with the ternary operator:
idx += (cond ? Steps : 0); // old, generated code uses conditional moves
idx = cond ? idx + Steps : idx; // new, generated code uses shift-and-or clang can see that the accumulated index is a bitfield to which each probe contributes exactly one bit, and builds it by shift-and-or — five instructions per probe, 53 in total, with one conditional move in the whole function. When we rewrite as an assignment, the same computation lowers to lea, cmp, cmovg — three instructions per probe, 35 in total, ten conditional moves. The change reduces runtime from 5.27 ns to 4.07 ns, closing the gap entirely.
The lambda from earlier arrives at the same place without the awkward syntax. if (cond) idx += step is the selection, so the readable formulation and the fast one turn out to be the same — clang emits the identical thirty-five instructions for both, and g++ emits forty-five for both, differing only in the sense of the conditional moves. Perfecting the codegen and getting the source readable were in conflict after all, which finalizes the argument for the version we shipped.
g++ produces conditional moves from either form and never exhibits the difference, which is why the defect would have gone unnoticed if I’d only built with one compiler.
With all three applied, the variadic version is faster than the plain loop at g++ -O2 by 22% — 3.91 ns against 4.99 — and level with it at -O3 and under clang. The probe offsets are compile-time constants by construction rather than by the grace of an unrolling pass, which g++ does not run at -O2. That is the one durable advantage of the formulation, and as we’ll see in the next article, this advantage is not limited to Binary Search.
Takeaway
In the final analysis, I got what I wanted: working, reasonably idiomatic code for a language feature I had never had occasion to use. One thing we can take away is that coding agents still need human intervention, since I repeatedly asked Claude to correct or clarify its initial drafts of code. It was able to shore up my weakness on syntax and focus on the semantics and the test data, but its lapses underscored that engineering judgment is as indispensable as ever.
Here is the code produced by the day’s work:
#include <array>
#include <cstddef>
#include <type_traits>
#include <utility>
// Scalars are cheaper passed by value: taking a reference obliges the caller to
// materialise the argument in memory so that its address can be passed. Class
// types still want const&.
template <typename T>
using search_arg_t = std::conditional_t<std::is_scalar_v<T>, T, const T&>;
// Unrolled search over a segment whose length is a power of two. Steps are the
// descending powers of two below Len, so the accumulated index never runs past
// the end: before the probe at 2^j it is at most Len - 2^(j+1), and adding 2^j
// leaves it below Len. No bounds test is needed inside the descent.
template <std::size_t Len, typename T, std::size_t... Steps>
std::ptrdiff_t unrolled_search_segment(const T* arr, search_arg_t<T> target, std::integer_sequence<std::size_t, Steps...>) {
static_assert((std::size_t(0) | ... | Steps) == Len - 1,
"probe steps must be the distinct powers of two below Len");
std::size_t idx = 0;
// Written as a conditional assignment rather than idx += (cond ? Step : 0):
// clang lowers the accumulating form to shift-and-or bit insertion, which is
// slower than the conditional move it emits here.
auto probe = [&](auto step) {
if (arr[idx + step.value] <= target)
idx += step.value;
};
(probe(std::integral_constant<std::size_t, Steps>{}), ...);
(void)probe; // Len == 1 leaves the pack empty, so probe goes uncalled
if (arr[idx] == target) return static_cast<std::ptrdiff_t>(idx);
return -1;
}
// Helper: largest power of 2 <= N
constexpr std::size_t floor_power_of_two(std::size_t n) {
std::size_t p = 1;
while (p * 2 <= n) p *= 2;
return p;
}
template <std::size_t N, std::size_t... Steps>
constexpr auto make_probe_steps_impl(std::integer_sequence<std::size_t, Steps...>) {
constexpr std::size_t next = (N / 2) >> sizeof...(Steps);
if constexpr (next > 0) {
return make_probe_steps_impl<N>(std::integer_sequence<std::size_t, Steps..., next>{});
} else {
return std::integer_sequence<std::size_t, Steps...>{};
}
}
template <std::size_t N>
constexpr auto make_probe_steps() {
return make_probe_steps_impl<N>(std::integer_sequence<std::size_t>{});
}
template <typename T, std::size_t N>
std::ptrdiff_t unrolled_binary_search(const std::array<T, N>& arr, const T& target) {
constexpr std::size_t k = floor_power_of_two(N);
if constexpr (N == 0) {
return -1; // nothing to find; an empty range has no match
} else if constexpr (N == k) {
// Power of 2: unrolled search
return unrolled_search_segment<N>(&arr[0], target, make_probe_steps<N>());
} else {
constexpr std::size_t probe = N - k;
if (target <= arr[probe]) {
// Unrolled search in [0, k-1]
return unrolled_search_segment<k>(&arr[0], target, make_probe_steps<k>());
} else {
// Unrolled search in [N-k, N-1]
std::ptrdiff_t res = unrolled_search_segment<k>(&arr[N - k], target, make_probe_steps<k>());
return (res == -1) ? -1 : static_cast<std::ptrdiff_t>(N - k) + res;
}
}
}
// Example usage:
// constexpr std::array<int, 16> v = {1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16};
// std::ptrdiff_t idx = unrolled_binary_search(v, 7); // idx == 6
// Unrolled binary search for 1024 elements, for codegen comparison
template <typename T>
std::ptrdiff_t unrolled_binary_search_1024(const std::array<T, 1024>& arr, const T& target) {
constexpr std::size_t N = 1024;
constexpr auto steps = make_probe_steps<N>();
return unrolled_search_segment<N>(&arr[0], target, steps);
}
This code is intended to be reused: header-only, C++17, no dependencies beyond the standard library. It handles any array size — sizes that are not powers of two are covered by two overlapping searches of the largest power of two that fits, chosen by a single comparison. It works for any type with <= and ==; the tests exercise it with std::string as well as int32_t. And the probe sequence is generated at compile time, so what the compiler emits for a given size is the same run of constant-offset probes you would have written by hand.
The header and its test harness are on GitHub, in static-bsearch, under a BSD license. Take it, read it, use it. I would be glad to hear about bugs; the test harness now draws its data by Floyd sampling from the whole int32_t range, but Claude and I may have missed something.
Coming Next
In the next article, I’ll take variadic templates to a problem where the parameter pack is not a convenience, but a necessity: fusing chains of elementwise AVX-512 operators into a single pass over memory, where each operator is a distinct type and they cannot be iterated by a loop. In that code, there will be two parameter packs doing entirely different jobs, only one of which supplies arity — and the -O2 advantage is waiting there too, worth 47%.
Bentley is where I found it, but Knuth credits the technique to L. E. Shar, who proposed it in 1971: “Another modification of binary search, suggested in 1971 by L. E. Shar, will be still faster on some computers, because it is uniform after the first step, and it requires no table.” See The Art of Computer Programming, Volume 3, section 6.2.1, exercises 12 and 13. The uniform binary search that does keep a table is credited separately, to Ashok K. Chandra; Shar’s contribution was getting rid of it by making the probe offsets powers of two. No 1971 paper by Shar appears to exist, so it most likely reached Knuth by private communication — which makes TAOCP the citable source for a method nobody published. Bentley had referenced those same exercises earlier, in Writing Efficient Programs (1982).