Parallel Prefix Scan: Hillis-Steele Versus Blelloch

A scan takes an array and an associative operator and returns every partial fold: out[i] = x₀ ⊕ x₁ ⊕ … ⊕ xᵢ. It is the one parallel primitive that is genuinely not a reduction, and it is the substrate under compaction, radix sorting, run-length encoding, sparse-matrix construction, and parser combinator libraries. The naive sequential scan does n−1 operations and is trivially correct. The parallel versions trade differently, and the choice between them is one of the few places in parallel algorithm design where a constant factor in work is a real design decision rather than a detail.

Why the tree does not suffice

A reduction tree computes the total but discards the partial sums that a scan must emit. Keeping them all in a naive way requires every node of the tree to know the sum of everything to its left, which is not available locally. The two standard algorithms solve that in opposite ways: Hillis-Steele recomputes more, Blelloch saves and reuses.

Both need the operator to be associative, and neither needs it commutative if the implementation is careful about order. Both are expressed most cleanly on a power-of-two length and extended to arbitrary lengths by padding with the identity.

Hillis-Steele: doubling the distance

The algorithm runs ⌈log₂ n⌉ rounds, one for each stride d = 1, 2, 4, …, up to the largest power of two below n. In the round with stride d, every element at index i ≥ d replaces its value with x[i − d] ⊕ x[i]. All updates in a round read the array as it was at the start of the round, so they are independent and can run concurrently.

start:    x0  x1  x2  x3  x4  x5  x6  x7
d = 1:    x0  x0x1  x1x2  x2x3  x3x4  x4x5  x5x6  x6x7
d = 2:    x0  x0x1  x0..2 x1..3 x2..4 x3..5 x4..6 x5..7
d = 4:    x0  x0..1 x0..2 x0..3 x0..4 x0..5 x0..6 x0..7

After k rounds each position holds the fold of the 2^k elements ending at it (or all of them, near the start). Three rounds suffice for eight elements, and ⌈log₂ n⌉ in general. Work is Θ(n log n): each of the log n rounds touches almost every element. Span is Θ(log n). Parallelism is therefore Θ(n) — the full width of the array — which is why the algorithm is attractive on wide SIMD hardware.

The practical cost is that it does Θ(n log n) additions where Θ(n) are needed. On a GPU where the rounds are single wide instructions, that is often the right trade: the extra arithmetic is cheaper than the extra synchronisation a more complex algorithm would require. On a CPU with narrow vectors, the log n factor is real work and it shows.

#include <stddef.h>
#include <string.h>

/* Hillis-Steele inclusive scan, in place, n a power of two.
   Requires a scratch buffer so a round reads pre-round values. */
void scan_hs(int *x, int *tmp, size_t n)
{
    for (size_t d = 1; d < n; d <<= 1) {
        for (size_t i = 0; i < n; ++i)
            tmp[i] = (i >= d) ? x[i - d] + x[i] : x[i];
        memcpy(x, tmp, n * sizeof *x);
    }
}

The memcpy is the sequential stand-in for a barrier: a real parallel version must ensure no thread writes x[i] before every thread has read x[i−d] for that round. That barrier, once per round, is log n barriers total.

Blelloch: up-sweep and down-sweep

Blelloch’s scan is the reduction tree plus a second pass that distributes the prefixes back down. It runs in two phases over an array stored as a complete binary tree in level order, leaves at positions n/2…n−1.

The up-sweep is exactly a reduction: for d from 1 to log₂ n, every node at level d adds its two children. At the end the root holds the total, and every internal node holds the sum of its subtree.

The down-sweep reverses direction. Set the root to the identity, then for d from log₂ n down to 1, for each internal node at that level: save its value to a temporary, write the parent’s value into the node (the exclusive prefix of its subtree), and write parent ⊕ saved into its right child. After the last round, every leaf holds the exclusive prefix of everything before it.

up-sweep, array held in level order, internal nodes overwritten:
  start   [3  1  7  0  4  1  6  3]
  d = 1   [3  4  7  7  4  5  6  9]
  d = 2   [3  4  7 11  4  5  6 14]
  d = 4   [3  4  7 11  4  5  6 25]

down-sweep, root := identity, then the swap-and-add pass:
  start   [3  4  7 11  4  5  6  0]
  d = 4   [3  4  7  0  4  5  6 11]
  d = 2   [3  0  7  4  4 11  6 16]
  d = 1   [0  3  4 11 11 15 16 22]   <- exclusive prefixes

The final row is the exclusive prefix array: each entry is the sum of everything strictly before it, so index 5 holds 3+1+7+0+4 = 15, and the last entry holds 22, the sum of the first seven. The total, 25, is what the root held before it was overwritten with the identity.

Work is Θ(n): the up-sweep performs n−1 additions and the down-sweep performs n−1 more. Span is Θ(log n) — 2 log n rounds. So the algorithm is work-optimal and has the same span class as Hillis-Steele while doing Θ(log n) times less work. Parallelism is Θ(n/log n) rather than Θ(n), the price of the extra pass.

#include <stddef.h>

/* Blelloch exclusive scan on n a power of two, in place. */
void scan_blelloch(int *x, size_t n)
{
    /* up-sweep */
    for (size_t d = 1; d < n; d <<= 1)
        for (size_t i = 0; i < n; i += 2 * d)
            x[i + 2 * d - 1] += x[i + d - 1];
    x[n - 1] = 0;                          /* identity at the root */
    /* down-sweep */
    for (size_t d = n >> 1; d > 0; d >>= 1)
        for (size_t i = 0; i < n; i += 2 * d) {
            int t = x[i + d - 1];
            x[i + d - 1] = x[i + 2 * d - 1];
            x[i + 2 * d - 1] += t;
        }
}

Note the direction of the swap in the down-sweep: the parent’s value moves into the left child’s slot and the saved value is added into the right. Getting the two assignments in the wrong order destroys the saved value, which is the most common way this algorithm is written incorrectly. The result is exclusive; a shift-and-add converts it to inclusive.

Choosing between them, and the hybrid

The decision rests on three numbers: work, span, and the cost ratio between an arithmetic operation and a barrier.

	work	span	barriers	parallelism
Sequential	n	n	0	1
Hillis-Steele	n log n	log n	log n	Θ(n)
Blelloch	2n	2 log n	2 log n	Θ(n/log n)

At n = 2²⁰, Hillis-Steele does about 21 million additions and Blelloch about 2 million, a factor of ten. Blelloch pays 40 barrier rounds against Hillis-Steele’s 20. Put the arithmetic at one cycle per addition and eight cores: Hillis-Steele’s share is about 2.5 million cycles against Blelloch’s 0.25 million, and the extra 20 barriers, at a thousand cycles each, add roughly 20 thousand cycles. On those assumptions the work term dominates and Blelloch comes out close to nine times faster. Make the arithmetic cheaper instead, as happens when the scan is bound by streaming memory rather than by the ALU, and the barrier count becomes the term that decides. On a GPU with thousands of lanes, where a round is a few wide instructions and intra-warp synchronisation is cheap, the terms swap places: the extra rounds Blelloch needs cost more than the extra arithmetic Hillis-Steele does.

The usual hybrid avoids the choice: process the array in blocks that fit in shared memory or in a few cache lines, run a work-optimal scan inside each block, then scan the block totals with a small Blelloch tree, then add each block’s exclusive offset back into its elements. Total work is Θ(n), the number of barrier rounds is O(log(n/B)) for n/B blocks of size B, and the inner scan runs entirely in fast memory. This is the standard shape for large scans on both GPUs and CPUs, and it is the same decomposition that makes merge-path partitioning work for merging.

Where scans go wrong

Padding is the first trap. Padding with the identity gives correct inclusive results for the real elements, but if the exclusive offsets are then used to index into the original array, a padded slot may produce an offset past the end. The fix is to compute the scan over exactly n elements by treating the padding as an identity and clamping indices, not by trusting the padded result.

Non-associative operators are the second. A scan of “last element satisfying a predicate” is expressible with the operator that returns the right-hand value when it is valid and the left otherwise — associative, and it needs a distinguished invalid value. A scan of running maximum is associative; a scan of running average is not, and needs the (sum, count) pair trick described in the reduction article.

The third is overflow and saturation. A byte-wide scan of a 256-element array of ones overflows; the offsets are wrong and the error is silent. Widening the accumulator, or scanning in blocks small enough that the block total fits, is the usual remedy.

The fourth is the difference between inclusive and exclusive, which libraries disagree about in naming. A scan that returns “the sum of everything before me” and one that returns “the sum of everything up to and including me” differ by one element, and code that uses one where the other is intended produces an off-by-one that looks like an indexing bug. The block-level decomposition of a merge uses the same idea: split into cache-sized pieces, solve the piece boundaries with a small search, then fix up.