Wavelet Trees and Range Queries on Sequences

A static sequence of n symbols over an alphabet of size σ supports a range of queries: how many times does symbol c occur in positions [l, r); which symbol is the k-th smallest in that range; which positions hold c in that range. A per-symbol index of positions answers the third cheaply and the others badly. A wavelet tree answers all three in O(log σ) time while storing the sequence in n log σ bits plus a rank index, which for byte-sized alphabets is the same order as the raw bytes and for small alphabets is substantially less.

The structure is a binary trie over the alphabet’s bit patterns, with each node storing a bitvector that says which way each symbol in its subsequence went. The recursion is on the alphabet, not on the positions, which is what keeps the depth at log σ rather than log n.

The partition, level by level

Take the alphabet Σ and split it into two halves by the top bit of each symbol’s code: symbols whose leading bit is 0 go left, those with 1 go right. Store a bitvector over the whole sequence with a 1 at each position whose symbol went right. The left child then receives the subsequence of symbols with leading bit 0, in their original relative order, and the right child receives the complement. Recurse within each child on the next bit, until the alphabet has one symbol.

Two properties fall out of this construction and are the whole reason the queries work.

The bitvector’s rank tells you where a position went. If rank₁(i) counts the ones before i, then a position i whose bit is 0 maps to position rank₀(i) in the left child, and a position whose bit is 1 maps to position rank₁(i) in the right child. The mapping is monotone and preserves order.

A contiguous range stays contiguous. A range [l, r) in the parent maps to [rank₀(l), rank₀(r)) in the left child and [rank₁(l), rank₁(r)) in the right. This is the invariant that lets every query descend by narrowing a range rather than by looking up individual positions.

#include <stdint.h>
#include <stddef.h>

/* One level of the tree: the bitvector over the current subsequence plus
   its rank index. Levels are stored in a flat array, one per depth. */
struct level {
    const uint64_t *bits;
    const uint64_t *super;   /* rank at each 512-bit boundary */
    const uint16_t *sub;     /* bits before this word within its superblock */
};

/* Count of symbol c in the positions [l, r) of the whole sequence,
   over a byte alphabet: the descent has eight levels. */
static size_t range_count(const struct level *lv, size_t l, size_t r,
                          unsigned c) {
    for (unsigned bit = 0; bit < 8; bit++) {
        size_t l1 = rank1(lv[bit].bits, lv[bit].super, lv[bit].sub, l);
        size_t r1 = rank1(lv[bit].bits, lv[bit].super, lv[bit].sub, r);
        if (((c >> (7 - bit)) & 1) == 0) { l -= l1; r -= r1; }  /* zeros map left */
        else                             { l = l1; r = r1; }    /* ones map right */
    }
    return r - l;                     /* the leaf range holds only symbol c */
}

At the leaf the range contains only occurrences of c, so its length is the count. The loop performs two rank queries per level, eight levels for a byte alphabet, for sixteen rank operations per range count. That is the cost model to keep in mind: O(log σ) levels, each with a constant number of rank queries, each rank query a handful of instructions.

Quantile and the selection descent

The k-th smallest symbol in [l, r) is found by descending and counting. At each level, the number of zeros in the mapped range is z = (r0 - l0). If k < z, the answer lies in the left child and the range becomes [l0, r0) with the same k. Otherwise the answer lies in the right child, the range becomes [l1, r1), and k decreases by z.

/* k-th smallest (0-based) symbol in [l, r), byte alphabet. */
static unsigned quantile(const struct level *lv, size_t l, size_t r, size_t k) {
    unsigned sym = 0;
    for (unsigned bit = 0; bit < 8; bit++) {
        size_t l1 = rank1(lv[bit].bits, lv[bit].super, lv[bit].sub, l);
        size_t r1 = rank1(lv[bit].bits, lv[bit].super, lv[bit].sub, r);
        size_t zeros = (r - l) - (r1 - l1);       /* symbols that went left */
        if (k < zeros) {
            l -= l1;  r -= r1;                    /* descend left */
        } else {
            k -= zeros;
            l = l1;  r = r1;                      /* descend right */
            sym |= 1u << (7 - bit);
        }
    }
    return sym;
}

The same descent answers “how many symbols are less than c in [l, r)” by accumulating the zero counts whenever the descent goes right, which is the rank query on the value axis rather than the position axis.

The reason both work is that the structure is simultaneously a trie on values and a sequence on positions, and the two views are linked by the order-preserving rank mapping. Nothing else in the structure needs to be stored.

Space: n log σ bits plus the rank indexes

The bitvectors total exactly n bits per level, summed over all nodes at that level, because each symbol appears exactly once at each level. With ⌈log₂ σ⌉ levels the total is n ⌈log₂ σ⌉ bits. For bytes that is 8n bits, or n bytes, exactly the size of the raw input. For a four-symbol alphabet such as nucleotide codes it is 2n bits, which matches the obvious two-bits-per-symbol packing rather than beating it.

The rank indexes add roughly 0.3 bits per stored bit, as computed in Bitset Representation and Word-Level Rank and Select, so the total is about 1.3 n ⌈log₂ σ⌉ bits. At σ = 256 that is 10.4 bits per symbol — larger than the input. The structure is therefore not a compressor for byte alphabets; its value there is the query capability, not the space. It becomes a space win only when the alphabet’s effective size is well below the code width, as with text over a restricted character set, or when the sequence itself is compressible and the bitvectors are stored compressed, which is a different structure with O(log σ) unchanged but query costs that rise with the compression rate.

The σ in the space bound is the alphabet size, not the number of distinct symbols that actually occur. A sequence of 8-bit codes with only 20 distinct values stored at depth 8 wastes three-quarters of the space. Two standard remedies exist: remap the alphabet to its distinct symbols first, which requires a dictionary and is only worth it when the distinct count is known to be small; or store the tree at the depth implied by the code width and accept the waste.

Range reporting and the position list

The third query — which positions hold c in [l, r) — is answered by descending to the leaf for c and then walking back up. At the leaf, the occurrences of c occupy a contiguous range [a, b) of the leaf’s subsequence, and the number of them in the query range is known. The positions in the original sequence are recovered by inverting the rank mapping at each level on the way up: a position p in a child’s subsequence corresponds to the p-th occurrence of its bit in the parent, which is select of p on the parent’s bitvector.

That is why select is needed as well as rank, and why the two operations are treated as a pair. The inverse walk costs O(log σ) select operations per reported position, so reporting k positions costs O(k log σ).

A structure that supports reporting is necessarily slower than one that only counts, because each reported position needs the full inverse path. Applications that need only counts — approximate membership, selectivity estimation — use the count-only descent and never pay for select indexes at all, which removes a third of the space overhead.

Where a wavelet tree is the wrong tool

Three cases where the structure loses to something simpler.

Dynamic sequences. Inserting a symbol shifts every subsequent position at every level, and the bitvectors are not built for that. Mutable alternatives are the wavelet matrix over a dynamic bitvector (with the same log σ query cost but a much larger constant), or a balanced BST keyed by position with subtree counts. Neither is a small change.

Point queries only. If the only question is “what is the symbol at position i”, the raw array is ⌈log₂ σ⌉ bits per symbol and one load, which beats the tree on both counts. The wavelet tree pays for range structure it does not use.

Very large alphabets with sparse occurrence. For σ = 2³² and n small, the depth is 32 and the tree stores 32 bits per symbol plus indexes — worse than a hash map from symbol to position list for every query except range counting.

The structure’s real home is the static sequence with a mixed query load: range counts, quantiles, and range reporting over the same data. That combination appears in compressed text indexes, in columnar stores answering predicate selectivity, and in sequence aligners indexing a reference. The memory layout of the bitvectors, and the choice between storing them level by level and interleaving them, is the same trade discussed in Cache Behaviour and Node Layout in Disjoint-Set Arrays: the descent visits one bitvector per level, so levels stored contiguously give one cache line per level while interleaved levels give one line per query.

A variant worth knowing about is the wavelet matrix, which stores all levels in a single bitvector array with a per-level offset instead of a tree of nodes. It has the same query algorithms and the same space, but the levels are contiguous and the node boundaries disappear, which removes a pointer indirection per level. For alphabets of size 2^k the two are equivalent in query cost; the matrix is simpler to implement and friendlier to the cache.

The next article closes the cluster by applying the same rank/select machinery to trees rather than sequences: Balanced Parentheses and Succinct Tree Navigation.