Rank + Select Structures, in depth: conventions, constant-time select, select-in-word with PDEP, cache-aware layouts and an oracle test harness

Fix the conventions first

Most rank/select bugs are disagreements about conventions, not wrong algorithms. Pick one set, write it down, and enforce it with tests. This page uses the most common set, the one sdsl-lite and most papers use:

• Positions are 0-based. Bit i of the vector is bit i % 64 of word i // 64, least significant bit first.
• rank1(i) counts ones in the half-open range [0, i). So rank1(0) = 0 and rank1(n) is the total number of ones. Valid i runs from 0 to n inclusive.
• select1(j) is 1-based in j: it returns the position of the j-th one, for 1 ≤ j ≤ ones. Calling it with j = 0 or j > ones is an error, not a sentinel.
• rank0(i) = i - rank1(i). select0 needs its own sample table, because zeros are not where the ones are.

Two identities follow, and they make the best tests you will write: rank1(select1(j)) == j - 1 for every valid j, and the bit at select1(j) is always 1. If a library you depend on uses inclusive rank, wrap it once at the boundary instead of adjusting call sites one at a time.

Worked example: 16 bits by hand

Take n = 16 bits, written position 0 first: 1 0 1 1 | 0 0 1 0 | 1 1 1 0 | 0 0 0 1. The ones are at positions 0, 2, 3, 6, 8, 9, 10 and 15, eight in total. For the example, use blocks of 4 bits instead of 64, so each step stays small enough to check by hand. Block popcounts are 3, 1, 3, 1. The cumulative counts before each block are 0, 3, 4, 7.

Check the identity: rank1(select1(5)) = rank1(8) = 4 = 5 - 1. Rank is a lookup plus a masked popcount. Select is a search for the block, then a search within it. The rest of the engineering makes that second search take constant time and touch as few cache lines as possible.

Constant-time select by sampling

Binary search over rank costs O(log n) directory probes. For a billion bits, that is about thirty probes, and most of them miss cache. The standard fix is sampling. Record the position of every S-th one, say S = 4096. To answer select1(j), look up sample (j - 1) // S. That gives a superblock where the search can start. Scan forward over superblock counts until the next one would exceed j. Then scan the words of that superblock, and finish with a select inside one word.

The scan is short when ones are dense. When they are sparse, it can be long: 4,096 ones spread over millions of bits means thousands of superblocks between two samples. Clark's construction, used by sdsl-lite's select_support_mcl, handles this by classifying each sample interval. A long, sparse interval stores the positions of its ones explicitly. That costs little memory, because sparse intervals have few ones. A short, dense interval keeps a smaller second-level directory. Either way, the final scan is bounded by a constant number of words.

A reference implementation

Here is a reference implementation in Python. It follows the conventions above and uses the sampled select. It is meant for a test oracle and for teaching. It is not meant to be fast.

SB = 512          # bits per superblock
SAMPLE = 4096     # sample every SAMPLE-th one

class RankSelect:
    def __init__(self, words, n):
        self.w, self.n = words, n
        self.sb, self.rel, total = [], [], 0
        for k, x in enumerate(words):
            if k % (SB // 64) == 0:
                self.sb.append(total)          # 64-bit absolute count
            self.rel.append(total - self.sb[-1])  # fits in 9 bits: at most 448
            total += x.bit_count()
        self.ones = total
        # samples[t]: superblock holding the (t * SAMPLE + 1)-th one
        self.samples, s = [], 0
        for t in range(0, total, SAMPLE):
            while s + 1 < len(self.sb) and self.sb[s + 1] <= t:
                s += 1
            self.samples.append(s)

    def rank1(self, i):
        assert 0 <= i <= self.n
        k, r = divmod(i, 64)
        if k == len(self.w):
            return self.ones
        return self.sb[k // 8] + self.rel[k] + (self.w[k] & ((1 << r) - 1)).bit_count()

    def select1(self, j):
        if not 1 <= j <= self.ones:
            raise IndexError(j)
        s = self.samples[(j - 1) // SAMPLE]
        while s + 1 < len(self.sb) and self.sb[s + 1] < j:
            s += 1                             # bounded in production by Clark's scheme
        k, last = s * 8, min(len(self.w), s * 8 + 8) - 1
        while k < last and self.sb[s] + self.rel[k + 1] < j:
            k += 1                             # next word still starts before the j-th one
        need = j - self.sb[s] - self.rel[k]    # 1-based rank inside word k
        return k * 64 + select_in_word(self.w[k], need - 1)

def select_in_word(x, k):
    """Position of the (k+1)-th set bit of x, k counted from 0."""
    for _ in range(k):
        x &= x - 1                             # clear lowest set bit
    return (x & -x).bit_length() - 1

The in-superblock loop advances k while the next word still starts before the j-th one. A production version stores the eight relative counts in one 64-bit word and finds the right one with a broadword comparison. The structure is the same either way: a sample lookup, a short scan over counts, and one select inside a word.

Select inside one word

The last step is to find the k-th set bit inside one 64-bit word. On x86 with BMI2, two instructions do it. pdep deposits the bit 1 << k into the k-th set position of the word, and tzcnt reads off where it landed:

#include <stdint.h>
#if defined(__BMI2__)
#include <immintrin.h>
#endif

/* k is 0-based: k = 0 finds the lowest set bit. Requires k < popcount(w). */
static inline unsigned select64(uint64_t w, unsigned k) {
#if defined(__BMI2__)
    return (unsigned)__builtin_ctzll(_pdep_u64(1ULL << k, w));  /* nonzero: k < popcount */
#else
    /* Portable fallback: skip whole bytes by popcount, then clear low bits. */
    unsigned base = 0;
    for (;;) {
        unsigned c = (unsigned)__builtin_popcountll(w & 0xFF);
        if (k < c) break;
        k -= c; w >>= 8; base += 8;
    }
    while (k--) w &= w - 1;
    return base + (unsigned)__builtin_ctzll(w);
#endif
}

Two portability traps hide in this function. First, AMD processors before Zen 3 implement pdep in microcode, with latency that depends on the data and is far slower than on Intel. A binary compiled with -mbmi2 and benchmarked on Intel can be the slowest option on a Zen 2 fleet. Dispatch at runtime on the CPU model, or keep the portable path for those machines. Second, __builtin_popcountll compiled without -mpopcnt or a matching -march becomes a library call or a bit-twiddling sequence. That slows every rank query, not only select. Check the disassembly for a popcnt instruction. On ARM there is no pdep, so the byte-skip fallback is the normal path.

Layouts and cache misses

Popcounts and shifts are cheap. Cache misses are not. A query against a large vector costs roughly one miss per array it touches, so count arrays when you choose a layout.

When the vector fits in last-level cache, simpler layouts win, because instruction count dominates. When it is much larger, the layout that touches the fewest lines wins, and huge pages help because TLB misses grow too. For batch lookups, prefetch the next query's superblock while finishing the current one.

RRR compression stores each small block as its popcount, called its class, plus an offset naming which arrangement of that many ones it is. Skewed vectors shrink well below n bits; random ones do not, and every query pays to decode. sdsl-lite's rrr_vector is the standard implementation. For sorted integer sets, Elias-Fano is often the better compressed choice, and Roaring bitmaps win when you need set operations rather than rank and select.

Testing against an oracle

An oracle-based harness catches nearly every bug in a few seconds. Build the vector from a random bit pattern. Answer every query against a plain prefix-sum list. Make sure the patterns include the shapes that break directories: empty vectors, all zeros, all ones, a single one at the very end, lengths one less and one more than a word and a superblock, and sparse vectors whose gaps exceed the sample interval.

import random

def oracle(bits):
    pre = [0]
    for b in bits:
        pre.append(pre[-1] + b)
    pos = [i for i, b in enumerate(bits) if b]
    return pre, pos

def pack(bits):
    words = [0] * ((len(bits) + 63) // 64)
    for i, b in enumerate(bits):
        words[i // 64] |= b << (i % 64)        # LSB first, matches the convention
    return words

def check(bits):
    rs = RankSelect(pack(bits), len(bits))
    pre, pos = oracle(bits)
    for i in range(len(bits) + 1):
        assert rs.rank1(i) == pre[i], ("rank", i)
    for j, p in enumerate(pos, start=1):
        assert rs.select1(j) == p, ("select", j)
        assert rs.rank1(p) == j - 1
    for bad in (0, len(pos) + 1):
        try:
            rs.select1(bad); raise AssertionError("no error for j=%d" % bad)
        except IndexError:
            pass

for n in (0, 1, 63, 64, 65, 511, 512, 513, 20000):
    for density in (0.0, 0.001, 0.5, 1.0):
        check([1 if random.random() < density else 0 for _ in range(n)])
check([0] * 100000 + [1])                      # gap far beyond one sample

Run the same harness against your fast implementation, and add a test that queries a vector serialized on a different machine to catch byte-order mistakes.

Failure modes

• Inclusive versus exclusive rank. Mixing a library that returns ones in [0, i] with code that expects [0, i) is off by one exactly when bit i is set. That is half the time, so the bug looks random.
• Garbage in the last word. If the bits past n are not zero, rank1(n) and the total count are wrong. Mask the tail when you build the vector.
• Counter overflow. 32-bit superblock counts silently wrap past 2^32 ones. A 16-gigabit vector is not exotic in genomics or web-graph work. Use 64-bit absolute counts, or add a top-level count for every 2^32 bits.
• Relative-count width. A relative count must hold the number of ones before the last word of a superblock: 448 for 512-bit superblocks. That needs 9 bits. Grow the superblock without widening the field, and it wraps.
• Select on an empty or exhausted vector. select1(0) and select1(ones + 1) should raise. Returning n as a sentinel hides bugs in callers.
• Unbounded scans. Plain sampling with no sparse-interval handling is fast in benchmarks on random data and pathologically slow on real, clustered data.
• Mutation. These are static structures. Flipping one bit invalidates every count after it. If you need updates, use a Fenwick tree over word popcounts, or rebuild in batches.

Trade-offs

For a vector that fits in cache, use the simple 25% layout. For one far larger than cache, use an interleaved 3% to 6% layout with huge pages. Build only the select tables you call, and compress only after measuring the vector's entropy.

Rank and select also connect the structures on this site. Wavelet trees call rank once per level, so a 2x faster rank makes every range query 2x faster. Bit tricks such as SWAR popcount and lowest-bit isolation are catalogued in the bit hacks article.

What to do next

1. Write down your conventions: 0-based positions, LSB-first words, half-open rank, 1-based select. Put them in the header of the code.
2. Port the Python class as your oracle and the harness as your test suite before you write any fast code.
3. Implement rank with 64-bit superblock counts and packed relative counts. Check the disassembly for a hardware popcnt.
4. Add sampled select with sparse-interval handling, and test it with a vector whose gaps exceed the sample interval.
5. Dispatch pdep at runtime. Benchmark on every CPU family in your fleet, including any AMD Zen 2 machines.
6. Benchmark with random queries on a vector at least 4x larger than last-level cache. Count cache misses with perf stat before choosing a layout.