SWAR: Parallel Bit Operations In A Single Register

In the previous article, I mentioned that finding the n-th set bit within a machine word is the main bottleneck in select queries on some architectures, and promised a follow-up on Sebastiano Vigna's solution (the article - Broadword Implementation of Rank/Select Queries).

The algorithm is built on SWAR (SIMD Within A Register). SWAR is a way to process multiple small integers packed inside a single machine word using ordinary arithmetic. A 64-bit register becomes 8 independent bytes, or 16 four-bit lanes, or 32 two-bit lanes, whatever the algorithm needs.

Vigna's select compiles to ~60 branchless instructions with zero loads, and since modern (and not very modern too) processor cores execute several of these operations in parallel, a call takes ~15 cycles regardless of the input. We'll get to the benchmarks at the end.

In this article I am going to start from the basics required to understand Vigna's algorithm and to build up intuition for similar problems.

SWAR basics

Bitwise operations (AND, OR, XOR, NOT) already work per-lane for free. Every bit is independent, so lane boundaries don't matter.

Arithmetic is the problem, addition and subtraction introduce carries across lane boundaries. There are a few options to solve that

Lane Constants

Most SWAR operations rely on two constants defined by the lane width $k$ :

L_k = \underbrace{00\dots01}_{k}\;\underbrace{00\dots01}_{k}\;\cdots \qquad H_k = \underbrace{10\dots00}_{k}\;\underbrace{10\dots00}_{k}\;\cdots

$L_k$ has the lowest bit of every lane set. $H_k$ has the highest bit of every lane set. For 8-bit lanes in a 64-bit word:

$L_8$ = 0x0101010101010101
$H_8$ = 0x8080808080808080

We'll also use:

${\sim}H_k$ selects everything except the highest bit per lane.
${\sim}L_k$ selects everything except the lowest bit per lane.

fn L(comptime T: type, comptime k: u6) T {
    comptime {
        var acc: T = 0;
        var i = 0;
        while (i < @bitSizeOf(T)) : (i += k) acc |= 1 << i;
        return acc;
    }
}

fn H(comptime T: type, comptime k: u6) T {
    return L(T, k) << (k - 1);
}

pub fn main() !void {
    const print = @import("std").debug.print;
    print("L8 = 0x{X:0>16}\n", .{comptime L(u64, 8)});
    print("H8 = 0x{X:0>16}\n", .{comptime H(u64, 8)});
    print("L4 = 0x{X:0>16}\n", .{comptime L(u64, 4)});
    print("H4 = 0x{X:0>16}\n", .{comptime H(u64, 4)});
}

What multiplication by $L_k$ does

Multiplying a value by $L_k$ creates shifted copies of it, one per lane, and adds them together:

x \times L_8 = x + (x \ll 8) + (x \ll 16) + \cdots + (x \ll 56)

This has two uses depending on what $x$ looks like:

1. Spread a value to all lanes. If $x < 2^k$ (fits in a single lane), the result is $x$ duplicated into every lane. For example, $5 \times L_8$ = 0x0505050505050505. One multiply creates a mask for all 8 bytes at once.

2. Prefix sums across lanes. If each lane of $x$ holds a count that fits in $k$ bits, and the running sums never exceed $2^k - 1$ , byte $i$ of the result accumulates all bytes $0$ through $i$ .

Example

     x = 00 00 00 00 03 01 02 05
	   * 01 01 01 01 01 01 01 01
	     -----------------------
       = 00 00 00 00 03 01 02 05
       + 00 00 00 03 01 02 05 00
       + 00 00 03 01 02 05 00 00
       + 00 03 01 02 05 00 00 00
       + 03 01 02 05 00 00 00 00
       + 01 02 05 00 00 00 00 00
       + 02 05 00 00 00 00 00 00
       + 05 00 00 00 00 00 00 00
x * L8 = 0B 0B 08 07 06 03 07 05

What if the sums do overflow? In such case we can use the even/odd splitting trick to calculate two prefix sums and then add them. A few more operations than a single multiply, but it will still be faster than a loop because we don't have a single dependency (accumulator) for each instruction

What multiplication by $H_k$ does

$H_k = L_k \ll (k-1)$ , so multiplying by $H_k$ is the same as multiplying by $L_k$ then shifting left by $k-1$ . This is less useful than multiplying by $L_k$

Block the carry

How do we add packed values without carries leaking between lanes? Take two packed words with 4-bit lanes:

\begin{align} a &= \underbrace{0101}_{a_0}\;\underbrace{1011}_{a_1} &(5, 11) \\ b &= \underbrace{0100}_{b_0}\;\underbrace{0110}_{b_1} &(4, 6) \end{align}

$a_0 + b_0 = 9$ , fits in 4 bits. But $a_1 + b_1 = 17$ , which needs 5 bits. If we add the whole byte naively, the carry from $a_1 + b_1$ spills into $a_0$ 's lane and corrupts it.

To fix that we'll zero the high bit of each lane before adding.

a' = a \mathbin{\&} {\sim}H_k \qquad b' = b \mathbin{\&} {\sim}H_k

Without the high bit, each value is at most $2^{k-1} - 1$ . The sum of two such values is at most $2^k - 2$ , which fits in $k$ bits and carry can't cross into the next lane.

Example

~H4          =  0111 0111
a            =  0101 1011  (5, 11)
b            =  0100 0110  (4, 6)
a' = a & ~H4 =  0101 0011  (5, 3)
b' = b & ~H4 =  0100 0110  (4, 6)
a' + b'      =  1001 1001  (9, 9)

Without masking, 1011 + 0110 = 10001, a 5-bit result that would have overflowed into $a_0$ . After masking, 0011 + 0110 = 1001, which stays in its lane.

But the high bits we removed still need to be accounted for. How to put them back depends on the operation: wrapping addition uses them in XOR for modulo arithmetic, saturating addition uses them to detect overflow.

Even/Odd Lane Splitting

Some operations produce results wider than the input lanes. Multiplying two 8-bit values can produce a 16-bit result. In these cases we need more room.

The easiest solution is to split lanes with even and odd indexes into two separate words, each with double-width lanes:

\begin{align} mask &= L_{2k} \cdot (2^k - 1)\\ a_{e} &= a \mathbin{\&} mask \\ a_{o} &= (a \gg k) \mathbin{\&} mask \end{align}

$mask$ here selects every other $k$ -bit lane. Even lanes are already sitting at $2k$ -bit intervals, so they need no shifting. Odd lanes shift down by $k$ into the same layout.

For 4-bit lanes in a 16-bit word: $mask$ = 0000 1111 0000 1111

Example

a      = 1001 0101 0011 1100  (9, 5, 3, 12)
mask   = 0000 1111 0000 1111

a_e    = a & mask
       = 0000 0101 0000 1100       (_, 5, _, 12)

a_o    = (a >> 4) & mask
       = 0000 1001 0000 0011       (_, 9, _, 3)

Each 4-bit value now sits in an 8-bit lane with 4 bits of headroom.

After splitting, addition, subtraction and even multiplication all just work without carry blocking, because each result has $k$ extra bits of space.

Operations

Popcount

How do we count set bits in a word using SWAR? The strategy: start with the narrowest possible lanes, count bits within each lane, then merge adjacent lanes into wider ones until we have a single count.

Step 1: 1-bit lanes to 2-bit counts

Each bit is already its own popcount (0 or 1). We need to add adjacent pairs. A 2-bit value $x$ has a high bit and a low bit. The popcount is their sum:

\text{pop}_2(x) = x - ((x \mathbin{\&} H_2) \gg 1)

All four possible 2-bit values:

00 - (00 & 10) >> 1 = 00 - 0 = 00  (0 bits)
01 - (01 & 10) >> 1 = 01 - 0 = 01  (1 bit)
10 - (10 & 10) >> 1 = 10 - 1 = 01  (1 bit)
11 - (11 & 10) >> 1 = 11 - 1 = 10  (2 bits)

We mask with $H_2$ before shifting to prevent bits from leaking across lane boundaries.

Example

x                      =  01 11 10 00
x & H2                 =  00 10 10 00
(x & H2) >> 1          =  00 01 01 00
x - ((x & H2) >> 1)    =  01 10 01 00  (1, 2, 1, 0)

Step 2: 2-bit counts to 4-bit counts

Now each 2-bit lane holds a count (0, 1, or 2). Add adjacent pairs. The mask $L_4 \cdot 3$ = 0011 0011... selects the low 2 bits of every 4-bit group:

\text{pop}_4(x) = (pop_2(x) \mathbin{\&} (L_4 \times 3)) + ((pop_2(x) \gg 2) \mathbin{\&} (L_4 \times 3))

Shift odd 2-bit lanes down by 2, mask both halves, add. Each 4-bit result is at most 4, fits in the lane.

Step 3: 4-bit counts to 8-bit counts

Each 4-bit lane holds at most 4, so the sum of two adjacent lanes is at most 8. That fits in 4 bits so we can add first and mask after without odd-even split:

\text{pop}_8(x) = (pop_4(x) + (pop_4(x) \gg 4)) \mathbin{\&} (L_8 \times 15)

The mask $L_8 \times 15$ = 0x0F0F0F0F0F0F0F0F clears the garbage in the upper 4 bits of each byte.

Step 4: bytes to whole word

Each byte now holds the popcount of its 8 bits (0 to 8). To sum all bytes into a single value we multiply by $L_8$ :

$\text{popcount}(x) = (\text{pop}_8(x) \times L_8) \gg 56$

The multiply computes prefix sums across bytes (as described before in multiplication by $L_k$ ). The total is in the highest byte and we need to shift the word right by 56 to get the result (for 64 bit integers).

Code

const std = @import("std");
const print = std.debug.print;


fn L(comptime T: type, comptime k: u6) T {
    comptime {
        var acc: T = 0;
        var i = 0;
        while (i < @bitSizeOf(T)) : (i += k) acc |= 1 << i;
        return acc;
    }
}

fn H(comptime T: type, comptime k: u6) T {
    return L(T, k) << (k - 1);
}


fn popcount(comptime T: type, x: T) std.math.IntFittingRange(0, @bitSizeOf(T)) {
    comptime std.debug.assert(@bitSizeOf(T) >= 8 and @bitSizeOf(T) <= 128);

    const mask1 = comptime L(T, 4) * 3;
    const mask2 = comptime L(T, 8) * 15;

    var s = x - ((x & comptime H(T, 2)) >> 1);
    s = (s & mask1) + ((s >> 2) & mask1);
    s = (s + (s >> 4)) & mask2;

    return @intCast((s *% comptime L(T, 8)) >> (@bitSizeOf(T) - 8));
}

pub fn main() !void {
    print("popcount(0xFF00FF00FF00FF00) = {}\n", .{popcount(u64, 0xFF00FF00FF00FF00)});
    print("popcount(0x0123456789ABCDEF) = {}\n", .{popcount(u64, 0x0123456789ABCDEF)});
    print("popcount(0xFFFFFFFFFFFFFFFF) = {}\n", .{popcount(u64, 0xFFFFFFFFFFFFFFFF)});
    print("popcount(0xAAAAAAAA)         = {}\n", .{popcount(u32, 0xAAAAAAAA)});
}

The algorithm is only 12 operations with no branches or loops.

Addition

Take two packed words, each holding two 4-bit values:

\begin{align} a &= \underbrace{1001}_{a_0}\;\underbrace{0101}_{a_1} &(9, 5) \\ b &= \underbrace{0011}_{b_0}\;\underbrace{1100}_{b_1} &(3, 12) \end{align}

We want to add $a_0 + b_0$ and $a_1 + b_1$ simultaneously using operations on the whole byte.

$a_0 + b_0 = 12$ , fits in 4 bits, but $a_1 + b_1 = 17$ does not.

What to do in case of overflow?

Use addition with wrapping: $c_i = (a_i + b_i) \bmod 16$
Use addition with saturation: $c_i = \min(a_i + b_i, 15)$
Widen the output using even/odd splitting

Wrapping Addition

\text{wrap\_add}(a, b) = ((a \mathbin{\&} {\sim}H_k) + (b \mathbin{\&} {\sim}H_k)) \oplus ((a \oplus b) \mathbin{\&} H_k)

Derivation

We want $(a + b) \bmod 2^k$ per lane. To achieve that we need to use the carry-blocking pattern: mask out the high bits, add, recombine.

\begin{gather} a' = a \mathbin{\&} {\sim}H_k \qquad b' = b \mathbin{\&} {\sim}H_k \\ c = a' + b' \end{gather}

Now $c$ has the correct lower $k-1$ bits per lane, since the masked addition cannot overflow. We need to put the high bits back.

\begin{gather} h_a = a \mathbin{\&} H_k \qquad h_b = b \mathbin{\&} H_k \qquad h_c = c \mathbin{\&} H_k \\ d = c \mathbin{\&} {\sim}H_k \end{gather}

The correct result per lane is $(a + b) \bmod 2^k$ . We already have the correct lower bits in $d$ . The high bits were removed before addition, so we need to add them back:

\text{result} = (h_a + h_b + h_c + d) \bmod 2^k

$d$ is at most $2^{k-1} - 1$ and doesn't affect the high bit, so it passes through the mod operation.

\text{result} = (h_a + h_b + h_c) \bmod 2^k + d

Now we only need to reduce $h_a + h_b + h_c$ . Each of these numbers is either $2^{k-1}$ or $0$ .

All possible results of $(h_a + h_b + h_c) \bmod 2^k$ :

\begin{align} 0 \cdot 2^{k-1} \bmod 2^k &= 0 \\ 1 \cdot 2^{k-1} \bmod 2^k &= 2^{k-1} \\ 2 \cdot 2^{k-1} \bmod 2^k &= 2^k \bmod 2^k = 0 \\ 3 \cdot 2^{k-1} \bmod 2^k &= (2^k + 2^{k-1}) \bmod 2^k = 2^{k-1} \end{align}

The result is $2^{k-1}$ when an odd number of high bits are set. That's xor: $h_a \oplus h_b \oplus h_c$ .

Combining with $d$ : the result is $d \mid (h_a \oplus h_b \oplus h_c)$ , which can be written as $c \oplus ((a \oplus b) \mathbin{\&} H_k)$ .

const print = @import("std").debug.print;


fn L(comptime T: type, comptime k: u6) T {
    comptime {
        var acc: T = 0;
        var i = 0;
        while (i < @bitSizeOf(T)) : (i += k) acc |= 1 << i;
        return acc;
    }
}

fn H(comptime T: type, comptime k: u6) T {
    return L(T, k) << (k - 1);
}


fn wrapAdd(comptime k: u6, a: u64, b: u64) u64 {
    const hk = comptime H(u64, k);
    return ((a & ~hk) + (b & ~hk)) ^ ((a ^ b) & hk);
}

pub fn main() !void {
    const a: u64 = 0x95; // a = (9, 5)
    const b: u64 = 0x3C; // b = (3, 12)
    const result = wrapAdd(4, a, b);
    // (9 + 3) mod 16 = 12, (5 + 12) mod 16 = 1; expected: 0xC1
    print("wrapAdd(0x{X}, 0x{X}) = 0x{X}\n", .{ a, b, result });
}

Saturating Addition

\begin{align} z &= (a \oplus b) \mathbin{\&} H_k \\ c & = (a \mathbin{\&} {\sim}H_k) + (b \mathbin{\&} {\sim}H_k) \\ m &= ((a \mathbin{\&} b) \mid (z \mathbin{\&} c)) \mathbin{\&} H_k\\ \text{sat\_add}(a, b) &= (c \oplus z) \mid m \mid (m - (m \gg (k-1))) \end{align}

Derivation

Saturating addition clamps the result instead of wrapping: $c_i = \min(a_i + b_i, 2^k - 1)$ .

The difference from wrapping: instead of wrapping around on overflow, we force the lane to all ones ( $2^k - 1$ ).

Let's continue from wrapping add when we computed $h_a$ , $h_b$ , $h_c$ . Overflow happens when 2 or 3 of $h_a$ , $h_b$ , $h_c$ in a lane are set.

m = ((a \mathbin{\&} b) \mid (a \mathbin{\&} c) \mid (b \mathbin{\&} c)) \mathbin{\&} H_k

$m$ has the $k$ -th bit set in every lane that overflowed.

To saturate, we need to fill the entire overflowing lane with ones. Then shift $m$ down by $k-1$ to move the high bit to the lowest position, then subtract:

fill = m - (m \gg (k-1))

If $m$ = 0000 1000, this gives 0000 0111, combine that with $m$ , which is $2^k - 1$ in every overflowing lane, zero elsewhere, and we get the correct result: $m \mid fill$ = 0000 1111.

\text{sat\_add}(a, b) = \text{wrap\_add}(a, b) \mid m \mid (m - (m \gg (k-1)))

We can simplify $m$ by reusing variables from $\text{wrap\_add}$ . Using this fact: $x \mid y = (x \oplus y) \mid (x \mathbin{\&} y)$ :

\begin{align} m &= ((a \mathbin{\&} b) \mid (a \mathbin{\&} c) \mid (b \mathbin{\&} c)) \mathbin{\&} H_k \\ &= ((a \mathbin{\&} b) \mid ((a \mid b) \mathbin{\&} c)) \mathbin{\&} H_k\\ &= ((a \mathbin{\&} b) \mid (((a \oplus b) \mid (a \mathbin{\&} b)) \mathbin{\&} c)) \mathbin{\&} H_k\\ &= ((a \mathbin{\&} b) \mid ((a \oplus b) \mathbin{\&} c) \mid ((a \mathbin{\&} b) \mathbin{\&} c)) \mathbin{\&} H_k\\ &= ((a \mathbin{\&} b) \mid ((a \oplus b) \mathbin{\&} c)) \mathbin{\&} H_k \end{align}

Now let $z = (a \oplus b) \mathbin{\&} H_k$ . Since we only care about bits at $H_k$ positions (the final $\mathbin{\&} H_k$ ) we can replace $(a \oplus b)$ with $z$ :

m = ((a \mathbin{\&} b) \mid (z \mathbin{\&} c)) \mathbin{\&} H_k

$z$ is already computed in $\text{wrap\_add}$ , which means that we saved two operations.

const print = @import("std").debug.print;


fn L(comptime T: type, comptime k: u6) T {
    comptime {
        var acc: T = 0;
        var i = 0;
        while (i < @bitSizeOf(T)) : (i += k) acc |= 1 << i;
        return acc;
    }
}

fn H(comptime T: type, comptime k: u6) T {
    return L(T, k) << (k - 1);
}


fn satAdd(comptime k: u6, a: u64, b: u64) u64 {
    const hk = comptime H(u64, k);
    const z = (a ^ b) & hk;
    const c = (a & ~hk) +% (b & ~hk);
    const m = ((a & b) | (z & c)) & hk;
    return (c ^ z) | m | (m - (m >> (k - 1)));
}

pub fn main() !void {
    const a: u64 = 0x95; // a = (9, 5)
    const b: u64 = 0x3C; // b = (3, 12)
    const result = satAdd(4, a, b);
    // min(9 + 3, 15) = 12, min(5 + 12, 15) = 15; expected: 0xCF
    print("satAdd(0x{X}, 0x{X}) = 0x{X}\n", .{ a, b, result });
}

Widening Addition

Widening addition produces exact results by giving each lane more room. First we split both operands so each $k$ -bit value lands in a $2k$ -bit lane, then we can add normally:

\begin{align} M &= L_{2k} \cdot (2^k - 1)\\ c_{e} &= (a \mathbin{\&} M) + (b \mathbin{\&} M) \\ c_{o} &= ((a \gg k) \mathbin{\&} M) + ((b \gg k) \mathbin{\&} M) \end{align}

After that each value has $k$ extra bits of space, so the addition and multiplication can't overflow. The tradeoff is that we process two half-width words instead of one.

const std = @import("std");


fn L(comptime T: type, comptime k: u6) T {
    comptime {
        var acc: T = 0;
        var i = 0;
        while (i < @bitSizeOf(T)) : (i += k) acc |= 1 << i;
        return acc;
    }
}

fn H(comptime T: type, comptime k: u6) T {
    return L(T, k) << (k - 1);
}


fn widenAdd(comptime T: type, comptime k: u6, a: T, b: T) std.meta.Int(.unsigned, @bitSizeOf(T) * 2) {
    const W = std.meta.Int(.unsigned, @bitSizeOf(T) * 2);
    const mask = comptime L(T, 2 * k) * ((1 << k) - 1);
    const even = (a & mask) + (b & mask);
    const odd = ((a >> k) & mask) + ((b >> k) & mask);
    return @as(W, odd) << @bitSizeOf(T) | @as(W, even);
}

pub fn main() !void {
    const print = std.debug.print;
    const a: u64 = 0x95; // a = (9, 5)
    const b: u64 = 0x3C; // b = (3, 12)
    // 9 + 3 = 12, 5 + 12 = 17; exact results in 8-bit lanes
    print("widenAdd = 0x{X}\n", .{widenAdd(u64, 4, a, b)});
}

Subtraction

The same carry-blocking idea applies for subtraction, but we block borrow instead of carry.

Wrapping Subtraction

\text{wrap\_sub}(a, b) = ((a \mid H_k) - (b \mathbin{\&} {\sim}H_k)) \oplus ((a \oplus {\sim}b) \mathbin{\&} H_k)

We set the high bits of $a$ to 1 (instead of clearing them) so that borrow from the subtraction stays within the lane. Then we have to fix the high bits with xor, same as before.

Derivation

Setting $a \mid H_k$ guarantees that each lane's high bit is 1 before subtracting. Since $b \mathbin{\&} {\sim}H_k$ is at most $2^{k-1} - 1$ , the subtraction cannot borrow from the next lane.

After subtracting, the high bit of each lane holds one of three: the forced 1, the original high bit of $a$ and 0 if we borrowed 1. To clean this up we do the same xor trick as for wrapping add, using $a \oplus {\sim}b$ instead of $a \oplus b$ (negation here accounts for modulo 2 subtraction instead of modulo 2 addition).

const print = @import("std").debug.print;


fn L(comptime T: type, comptime k: u6) T {
    comptime {
        var acc: T = 0;
        var i = 0;
        while (i < @bitSizeOf(T)) : (i += k) acc |= 1 << i;
        return acc;
    }
}

fn H(comptime T: type, comptime k: u6) T {
    return L(T, k) << (k - 1);
}


fn wrapSub(comptime k: u6, a: u64, b: u64) u64 {
    const hk = comptime H(u64, k);
    return ((a | hk) - (b & ~hk)) ^ ((a ^ ~b) & hk);
}

pub fn main() !void {
    const a: u64 = 0x95; // a = (9, 5)
    const b: u64 = 0x3C; // b = (3, 12)
    const result = wrapSub(4, a, b);
    // (9 - 3) mod 16 = 6, (5 - 12) mod 16 = 9; expected: 0x69
    print("wrapSub(0x{X}, 0x{X}) = 0x{X}\n", .{ a, b, result });
}

Saturating Subtraction

\begin{align} z &= (a \oplus {\sim}b) \mathbin{\&} H_k \\ c &= (a \mid H_k) - (b \mathbin{\&} {\sim}H_k) \\ m &= (({\sim}a \mathbin{\&} b) \mid (z \mathbin{\&} {\sim}c)) \mathbin{\&} H_k \\ \text{sat\_sub}(a, b) &= (c \oplus z) \mathbin{\&} {\sim}(m \mid (m - (m \gg (k-1)))) \end{align}

Instead of filling overflowed lanes with all ones (like saturating add), we clear them to zero, since subtracting past zero saturates to 0.

const print = @import("std").debug.print;


fn L(comptime T: type, comptime k: u6) T {
    comptime {
        var acc: T = 0;
        var i = 0;
        while (i < @bitSizeOf(T)) : (i += k) acc |= 1 << i;
        return acc;
    }
}

fn H(comptime T: type, comptime k: u6) T {
    return L(T, k) << (k - 1);
}


fn satSub(comptime k: u6, a: u64, b: u64) u64 {
    const hk = comptime H(u64, k);
    const z = (a ^ ~b) & hk;
    const c = (a | hk) - (b & ~hk);
    const m = ((~a & b) | (z & ~c)) & hk;
    return (c ^ z) & ~(m | (m -% (m >> (k - 1))));
}

pub fn main() !void {
    const a: u64 = 0x95; // a = (9, 5)
    const b: u64 = 0x3C; // b = (3, 12)
    const result = satSub(4, a, b);
    // max(9 - 3, 0) = 6, max(5 - 12, 0) = 0; expected: 0x60
    print("satSub(0x{X}, 0x{X}) = 0x{X}\n", .{ a, b, result });
}

Comparison

Comparing all lanes simultaneously is the core of the "broadword" select. The result is a mask: $H_k$ bit set in each lane where the condition holds.

Equal to Zero

The simplest comparison. A lane is zero when all its bits are zero:

\text{eq0}(x) = {\sim}((x \mathbin{\&} {\sim}H_k) + {\sim}H_k) \mathbin{\&} {\sim}x \mathbin{\&} H_k

Adding ${\sim}H_k$ (all ones except high bits) to the lower bits of $x$ produces a carry into the high bit only if any lower bit was set. Combined with checking the original high bit, we detect all-zero lanes.

Equal / Not Equal

Equality reduces to eq0 with a XOR: two values are equal when their difference is zero.

\text{eq}(x, y) = \text{eq0}(x \oplus y)

Unsigned Less-Than-Or-Equal

This will be required in the final algorithm. Basically, this is just a simplified substraction:

\text{le}_k(x, y) = \left(((y \mid H_k) - (x \mathbin{\&} {\sim}H_k)) \oplus x \oplus y\right) \mathbin{\&} H_k

Setting $y \mid H_k$ prevents borrow from crossing lane boundaries. After subtracting, the high bit of each lane tells us whether $x \le y$ . The XOR with $x \oplus y$ cleans up the forced high bits.

const print = @import("std").debug.print;


fn L(comptime T: type, comptime k: u6) T {
    comptime {
        var acc: T = 0;
        var i = 0;
        while (i < @bitSizeOf(T)) : (i += k) acc |= 1 << i;
        return acc;
    }
}

fn H(comptime T: type, comptime k: u6) T {
    return L(T, k) << (k - 1);
}


fn leX(comptime k: u6, x: u64, y: u64) u64 {
    const hk = comptime H(u64, k);
    return (((y | hk) - (x & ~hk)) ^ x ^ y) & hk;
}

pub fn main() !void {
    // 8-bit lanes: x = (3, 7, 15, 2), y = (5, 7, 10, 8)
    const x: u64 = 0x03_07_0F_02;
    const y: u64 = 0x05_07_0A_08;
    const result: u32 = @truncate(leX(8, x, y));
    // 3<=5? yes, 7<=7? yes, 15<=10? no, 2<=8? yes
    // expected: 0x80800080
    print("leX = 0x{X}\n", .{result});
}

Unsigned Greater-Than (Non-Zero)

We'll also use a greater-than variant that checks $x > 0$ per lane, which is just $\text{ne0}$ for an unsigned case with the result in the high bit:

\text{gt0}(x) = (((x \mid H_k) - L_k) \mid x) \mathbin{\&} H_k

Min / Max

With comparison done, min and max are easy. The comparison gives $H_k$ bits, we expand them to full-lane masks and select between $a$ and $b$ :

\begin{align} m &= \text{le}_k(a, b) \\ M &= m \mid (m - (m \gg (k-1))) \\ \text{min}_k(a, b) &= b \oplus ((a \oplus b) \mathbin{\&} M) \\ \text{max}_k(a, b) &= a \oplus ((a \oplus b) \mathbin{\&} M) \end{align}

$m - (m \gg (k-1))$ is the same lane-filling step from saturating addition. We move the high bit of every "true" lane to the bottom, subtract, 100...0 - 1 = 011...1, or with $m$ gives all ones in the lane.

Then we do $b \oplus ((a \oplus b) \mathbin{\&} M)$ , it picks $a$ where the $M$ is set, $b$ elsewhere.

A shorter option reuses saturating subtraction:

\begin{align} \text{min}_k(a, b) &= a - \text{sat\_sub}(a, b) \\ \text{max}_k(a, b) &= b + \text{sat\_sub}(a, b) \end{align}

const print = @import("std").debug.print;


fn L(comptime T: type, comptime k: u6) T {
    comptime {
        var acc: T = 0;
        var i = 0;
        while (i < @bitSizeOf(T)) : (i += k) acc |= 1 << i;
        return acc;
    }
}

fn H(comptime T: type, comptime k: u6) T {
    return L(T, k) << (k - 1);
}



fn satSub(comptime k: u6, a: u64, b: u64) u64 {
    const hk = comptime H(u64, k);
    const z = (a ^ ~b) & hk;
    const c = (a | hk) - (b & ~hk);
    const m = ((~a & b) | (z & ~c)) & hk;
    return (c ^ z) & ~(m | (m -% (m >> (k - 1))));
}


fn minX(comptime k: u6, a: u64, b: u64) u64 {
    return a -% satSub(k, a, b);
}

fn maxX(comptime k: u6, a: u64, b: u64) u64 {
    return b +% satSub(k, a, b);
}

pub fn main() !void {
    // 4-bit lanes: a = (9, 5), b = (3, 12)
    const a: u64 = 0x95;
    const b: u64 = 0x3C;
    print("min = 0x{X}\n", .{minX(4, a, b)}); // min(9,3)=3, min(5,12)=5 -> 0x35
    print("max = 0x{X}\n", .{maxX(4, a, b)}); // max(9,3)=9, max(5,12)=12 -> 0x9C
}

Average

Unsigned average without overflow. The naive $(a + b) / 2$ overflows if the sum exceeds $2^k - 1$ , we'll use somthing different:

\text{avg}_k(a, b) = (a \mathbin{\&} b) + (((a \oplus b) \mathbin{\&} {\sim}L_k) \gg 1)

$a \mathbin{\&} b$ gives the bits where both are 1 (contributes 1 to each position). $a \oplus b$ gives the bits where they differ (contributes 0.5 each). We shift the XOR result right by 1 to divide it by 2. The mask ${\sim}L_k$ clears the lowest bit before shifting to prevent cross-lane leaks. The result is rounded down.

const print = @import("std").debug.print;


fn L(comptime T: type, comptime k: u6) T {
    comptime {
        var acc: T = 0;
        var i = 0;
        while (i < @bitSizeOf(T)) : (i += k) acc |= 1 << i;
        return acc;
    }
}

fn H(comptime T: type, comptime k: u6) T {
    return L(T, k) << (k - 1);
}


fn avgDown(comptime k: u6, a: u64, b: u64) u64 {
    const lk = comptime L(u64, k);
    return (a & b) + (((a ^ b) & ~lk) >> 1);
}

pub fn main() !void {
    // 8-bit lanes: a = (200, 100), b = (100, 50)
    const a: u64 = 0xC864;
    const b: u64 = 0x6432;
    print("avg = 0x{X}\n", .{avgDown(8, a, b)});
    // avg(200,100)=150, avg(100,50)=75 -> 0x964B
}

Cheatsheet

Everything from above in one table. $w$ is the word width, $k$ is the lane width, $m$ is a comparison result (a mask with bits at $H_k$ positions).

What	How	Notes
Broadcast $s$ to all lanes	$s \times L_k$	$s < 2^k$
Prefix sums across lanes	$x \times L_k$	running sums must fit in $k$ bits, details
Sum of all lanes	$(x \times L_k) \gg (w - k)$	same constraint
Wrapping add	$((a \mathbin{\&} {\sim}H_k) + (b \mathbin{\&} {\sim}H_k)) \oplus ((a \oplus b) \mathbin{\&} H_k)$	derivation
Wrapping sub	$((a \mid H_k) - (b \mathbin{\&} {\sim}H_k)) \oplus ((a \oplus {\sim}b) \mathbin{\&} H_k)$	derivation
Saturating add	wrapping add, then fill overflowed lanes	derivation
Saturating sub	wrapping sub, then clear underflowed lanes	derivation
Exact add / mul	even/odd split into $2k$ -bit lanes	details
Popcount	reduce 1 $\to$ 2 $\to$ 4 $\to$ 8 bit counts, then $\times L_8$	details
$x == 0$	${\sim}((x \mathbin{\&} {\sim}H_k) + {\sim}H_k) \mathbin{\&} {\sim}x \mathbin{\&} H_k$	result at $H_k$ positions
$x == y$	$\text{eq0}(x \oplus y)$
$x \le y$ unsigned	$(((y \mid H_k) - (x \mathbin{\&} {\sim}H_k)) \oplus x \oplus y) \mathbin{\&} H_k$
$x > 0$	$(((x \mid H_k) - L_k) \mid x) \mathbin{\&} H_k$
Min	$a - \text{sat\_sub}(a, b)$
Max	$b + \text{sat\_sub}(a, b)$
Average (floor)	$(a \mathbin{\&} b) + (((a \oplus b) \mathbin{\&} {\sim}L_k) \gg 1)$
$m$ $\to$ full-lane mask	$m \mid (m - (m \gg (k-1)))$	`100..0 - 1 = 011..1`
$m$ $\to$ count of true lanes	$((m \gg (k-1)) \times L_k) \gg (w - k)$	markers to low bits, prefix-sum, take top lane

Vigna's Broadword Select

Now we have all the tools.

The problem: given a word $x$ and a rank $r$ (0-based), find the position of the $r$ -th set bit.

The whole algorithm is two moves applied twice:

Build popcount prefix sums.
Compare all prefix sums against $r$ in parallel, count how many are $\le r$ - the count is the index of the lane holding our bit.

The first pass finds the target byte. The second pass runs on the bits of that byte, spread into byte lanes, so the same steps work again.

Helpers from the Comparison section:

// x <= y per lane, 0x80 marker where true
fn leX(comptime k: u6, x: u64, y: u64) u64 {
    const hk = comptime H(u64, k);
    return (((y | hk) - (x & ~hk)) ^ x ^ y) & hk;
}

// x > 0 per lane, 0x80 marker where true
fn gtX0(comptime k: u6, x: u64) u64 {
    const hk = comptime H(u64, k);
    return (((x | hk) - comptime L(u64, k)) | x) & hk;
}

Step 1: Byte Popcounts to Prefix Sums

const H2 = comptime H(u64, 2);      // 0xAAAAAAAAAAAAAAAA
const M1 = comptime L(u64, 4) * 3;  // 0x3333333333333333
const M2 = comptime L(u64, 8) * 15; // 0x0F0F0F0F0F0F0F0F

var s = x - ((x & H2) >> 1);
s = (s & M1) + ((s >> 2) & M1);
s = ((s + (s >> 4)) & M2) *% L8;

This is popcount stopped at byte level, with the same masks as before. The final *% L8 turns counts into prefix sums: byte $i$ of s holds the number of set bits in bytes $0..i$ of x.

Step 2: Find the Target Byte

var b = ((leX(8, s, r *% L8) >> 7) *% L8 >> 53);

r *% L8 broadcasts the rank to every byte and leX marks every byte where $s_i \le r$ . A marked byte has at most $r$ set bits up to and including itself - too few to reach the $r$ -th bit ( $r$ is 0-based, so it needs $r+1$ ), which must lie further right. Prefix sums never decrease, so the marks form a prefix $0..k-1$ . The first unmarked byte $k$ is the one holding our bit, and counting the markers gives exactly $k$ .

Instead of extracting the top byte with >> 56, we shift it by 53, keeping the count multiplied by 8 so byte index becomes bit offset, one shift instead of two.

Example

Per-byte popcounts of x are (byte 0 to 7): 5, 2, 1, 3, 0, 8, 1, 0, and $r = 9$ .

byte index        7  6  5  4  3  2  1  0

s               = 14 14 13 0B 0B 08 07 05  prefix sums
r *% L8         = 09 09 09 09 09 09 09 09  rank in every byte
leX(8, s, ...)  = 00 00 00 00 00 80 80 80  marker where s_i <= 9
>> 7            = 00 00 00 00 00 01 01 01  markers to low bits
*% L8           = 03 03 03 03 03 03 02 01  prefix-sum the markers
>> 53           = 24 = 3 * 8 = b           3 markers in the top byte,
                                           already scaled to a bit offset

Step 3: Rank Within the Byte

const l = r - ((std.math.shr(u64, s << 8, b)) & 0xFF);

s << 8 moves the prefix sums up one byte: byte $i$ now holds the count of set bits in bytes $0..i-1$ , byte $0$ holds zero. We shift right by b then mask with 0xFF and we get the number of set bits before the target byte. Then we subtract it from $r$ and l is the rank inside the target byte.

Example

Continuing from step 2: $r = 9$ , b = 24 (target byte 3).

byte index      7  6  5  4  3  2  1  0

s             = 14 14 13 0B 0B 08 07 05
s << 8        = 14 13 0B 0B 08 07 05 00  byte i: bits set before byte i
>> b, & 0xFF  = 08                       8 set bits before byte 3

l = r - 8     = 1                        our bit is the 2nd set bit
                                         of the target byte

Note that I use std.math.shr instead of >>. When $r \ge \text{popcount}(x)$ all 8 comparisons in step 2 pass and b becomes 64. Shifting a u64 by 64 is illegal in zig, std.math.shr returns 0 in that case.

Step 4: Find the Bit

s = (gtX0(8, ((std.math.shr(u64, x, b) & 0xFF) *% L8) & 0x8040201008040201) >> 7) *% L8;
b += ((leX(8, s, l *% L8) >> 7) *% L8 >> 56);

Same two moves, now at bit level. std.math.shr(u64, x, b) & 0xFF extracts the target byte, *% L8 copies it into all 8 lanes.

Then we use this constant:

0x8040201008040201 = 10000000 01000000 00100000 ... 00000100 00000010 00000001

Byte $k$ of the constant has a single set bit - bit $k$ . AND it with 8 copies of the target byte and byte $k$ stays non-zero exactly when bit $k$ of the target byte is set. Then gtX0 turns "non-zero" lanes into 0x80, >> 7 and *% L8 prefix-sum the markers, same as step 1.

The second line is step 2 again: count how many prefix sums are $\le l$ , this time shifting by 56 since we want a plain bit index, not times 8. Add it to the byte offset and b is the answer.

Not Found

When $r \ge \text{popcount}(x)$ step 2 gives b = 64, s in step 4 is zero, all bit-level prefix sums are zero and all 8 lanes pass le check and we add 8. So b == 72 is a natural "not found" result.

Complete Implementation

The implementation I use in miara, verified against a naive loop:

const std = @import("std");
const print = std.debug.print;


fn L(comptime T: type, comptime k: u6) T {
    comptime {
        var acc: T = 0;
        var i = 0;
        while (i < @bitSizeOf(T)) : (i += k) acc |= 1 << i;
        return acc;
    }
}

fn H(comptime T: type, comptime k: u6) T {
    return L(T, k) << (k - 1);
}

fn leX(comptime k: u6, x: u64, y: u64) u64 {
    const hk = comptime H(u64, k);
    return (((y | hk) - (x & ~hk)) ^ x ^ y) & hk;
}

fn gtX0(comptime k: u6, x: u64) u64 {
    const hk = comptime H(u64, k);
    return (((x | hk) - comptime L(u64, k)) | x) & hk;
}


pub fn nthSetBitPosU64Broadword(x: u64, r: u6) !u6 {
    const L8 = comptime L(u64, 8);      // 0x0101010101010101
    const H2 = comptime H(u64, 2);      // 0xAAAAAAAAAAAAAAAA
    const M1 = comptime L(u64, 4) * 3;  // 0x3333333333333333
    const M2 = comptime L(u64, 8) * 15; // 0x0F0F0F0F0F0F0F0F
    
    var s = x - ((x & H2) >> 1);
    s = (s & M1) + ((s >> 2) & M1);
    s = ((s + (s >> 4)) & M2) *% L8;

    var b = ((leX(8, s, r *% L8) >> 7) *% L8 >> 53);
    const l = r - ((std.math.shr(u64, s << 8, b)) & 0xFF);

    s = (gtX0(8, ((std.math.shr(u64, x, b) & 0xFF) *% L8) & 0x8040201008040201) >> 7) *% L8;
    b += ((leX(8, s, l *% L8) >> 7) *% L8 >> 56);

    if (b == 72) return error.NotFound;
    return @intCast(b);
}


fn nthSetBitPosNaive(x: u64, r: u6) !u6 {
    var bits = x;
    var remaining: u64 = r;
    for (0..64) |i| {
        if ((bits & 1) != 0) {
            if (remaining == 0) return @intCast(i);
            remaining -= 1;
        }
        bits >>= 1;
    }
    return error.NotFound;
}


pub fn main() !void {
    // bits of 0b10110 are at positions 1, 2, 4
    print("select(0b10110, 0) = {}\n", .{try nthSetBitPosU64Broadword(0b10110, 0)});
    print("select(0b10110, 1) = {}\n", .{try nthSetBitPosU64Broadword(0b10110, 1)});
    print("select(0b10110, 2) = {}\n", .{try nthSetBitPosU64Broadword(0b10110, 2)});

    var prng = std.Random.DefaultPrng.init(0x5eed);
    const rng = prng.random();
    for (0..100_000) |_| {
        const x = rng.int(u64);
        if (x == 0) continue;
        const r: u6 = @intCast(rng.uintLessThan(u64, @popCount(x)));
        const expected = try nthSetBitPosNaive(x, r);
        const actual = try nthSetBitPosU64Broadword(x, r);
        std.debug.assert(expected == actual);
    }
    print("100000 random words match the naive implementation\n", .{});
}

Performance

Setup: 10,000 random u64 words, a random valid rank for each, 1000 passes over the array, -O ReleaseFast. Four implementations:

Naive loop - the bit-by-bit scan from the previous article
Skip loop - clears the lowest set bit $r$ times, then @ctz:

var v = x;
for (0..r) |_| v &= v - 1;
return @ctz(v);

Broadword - the implementation above
BMI2 - PDEP + TZCNT, two instructions, x86-64 only

Benchmark source: bench.zig, zig run -O ReleaseFast bench.zig.

Apple M1:

100000 words, 1000 passes, ns per call
input                       naive     skip   broadword
sparse, ~4 set bits          29.3      6.4         4.9
random, ~32 set bits         89.1     11.4         5.1
dense, ~60 set bits          27.2     17.0         5.1

Intel i7-6700HQ (Skylake, 2.6 GHz):

100000 words, 1000 passes, ns per call
input                       naive     skip   broadword     bmi2
sparse, ~4 set bits          47.4     11.5        10.1      1.0
random, ~32 set bits        143.1     22.1        10.9      1.2
dense, ~60 set bits          54.2     27.7        10.2      1.1

Broadword: flat for every input on both machines. Nothing to mispredict, nothing to miss, the runtime doesn't depend on the data.

The naive loop is 3x slower at 50% density than on sparse or dense inputs, this is due to the impact of branch prediction errors.

The skip loop runs $r$ iterations, so it degrades linearly with rank, it's cheap on sparse words, and worse than broadword or a BMI2 version on dense ones where the average rank is ~30.

On x86 with fast BMI2 PDEP + TZCNT is 3x faster than broadword.

SWAR: Parallel Bit Operations In A Single Register

SWAR basics

Lane Constants

What multiplication by LkL_k does

What multiplication by HkH_k does

Block the carry

Even/Odd Lane Splitting

Operations

Popcount

Step 1: 1-bit lanes to 2-bit counts

Step 2: 2-bit counts to 4-bit counts

Step 3: 4-bit counts to 8-bit counts

Step 4: bytes to whole word

Code

Addition

Wrapping Addition

Saturating Addition

Widening Addition

Subtraction

Wrapping Subtraction

Saturating Subtraction

Comparison

Equal to Zero

Equal / Not Equal

Unsigned Less-Than-Or-Equal

Unsigned Greater-Than (Non-Zero)

Min / Max

Average

Cheatsheet

Vigna's Broadword Select

Step 1: Byte Popcounts to Prefix Sums

Step 2: Find the Target Byte

Step 3: Rank Within the Byte

Step 4: Find the Bit

Not Found

Complete Implementation

Performance

Further Reading

What multiplication by $L_k$ does

What multiplication by $H_k$ does