Parallel Prefix / Scan

Why Parallel Prefix / Scan Matters

Suppose five records need output regions of lengths $[3,1,7,0,4]$. Their starting offsets are $[0,3,4,11,11]$: each record begins after the lengths of all preceding records. That is an exclusive prefix sum. Including each record's own length instead gives cumulative ends $[3,4,11,11,15]$, an inclusive prefix sum.

A running accumulator computes prefixes sequentially, but its dependency chain appears to force every result to wait for the previous one. Parallel prefix, or scan, reorganizes the calculation with associative combinations so many partial reductions can proceed at the same time.

This lesson covers two classic approaches. Hillis–Steele doubles the covered interval at each level, giving a short dependency depth with extra repeated work. Blelloch first reduces a balanced tree upward, then propagates prefixes downward, giving linear total combine work. Its conventional result is exclusive; the supplied programs also convert it to inclusive for an exact comparison.

Python, C++, and Java simulate deterministic worker assignments, phase barriers, and complete intermediate arrays. The sample has five elements, so Blelloch visibly pads to eight identities. Both algorithms preserve order, handle empty/singleton inputs, and support a generic associative operator internally; additional tests use string concatenation to catch operand-order bugs hidden by addition.

Scans support output offsets, stable compaction, radix sorting, cumulative distributions, and many parallel primitives. The work/span ideas from the preceding lesson explain their tradeoff: a short dependency path is useful, while avoiding unnecessary total work becomes essential as input grows.

Core Concepts and Terminology

Let the ordered input be $[a_0,\ldots,a_{n-1}]$, with a binary operator $\oplus$ and identity e. Require associativity, $$ (x\oplus y)\oplus z=x\oplus(y\oplus z), $$ and a two-sided identity $e\oplus x=x\oplus e=x$. Together these form a monoid. Commutativity is not required: input order must stay unchanged. Addition with zero, multiplication with one, and string concatenation with the empty string are examples; subtraction does not meet the associativity requirement.

The inclusive output is $$ y_i=a_0\oplus a_1\oplus\cdots\oplus a_i. $$ The exclusive output is $x_0=e$ and $$ x_i=a_0\oplus\cdots\oplus a_{i-1}\quad(i\ge1). $$ Both have length n. A reduction returns only the total, whereas scan returns all prefixes. For nonempty input, $y_i=x_i\oplus a_i$ and the total is $y_{n-1}$. The empty input has empty inclusive/exclusive outputs; a separate total-reduction API could return e, but these programs do not add an extra output element.

Hillis–Steele levels. At stride $d=1,2,4,\ldots\lt n$, each active target $i\ge d$ computes $$ \mathrm{next}[i]=\mathrm{old}[i-d]\oplus\mathrm{old}[i], $$ while indices below d keep their old values. All reads come from the same previous-level snapshot. After this level, index i reduces the ordered interval from $\max(0,i-2d+1)$ through i. Doubling until the interval reaches the beginning produces inclusive scan.

Blelloch tree storage. For $n\gt0$, let $m=2^{\lceil\log_2n\rceil}$ and append $m-n$ identities. Padding does not change any valid prefix or the total. The up-sweep stores the reduction of each power-of-two interval at its rightmost index. Reset the root entry to e, then use the down-sweep to give every interval the reduction of all original elements strictly before it. Leaf entries become exclusive prefixes; trim to n.

Guy Blelloch's Prefix Sums and Their Applications develops scan as a general associative primitive and describes many uses. NVIDIA's scan chapter illustrates Hillis–Steele and the work-efficient tree approach. The states and counts below come from the supplied programs, with explicit padding and conversion conventions.

Worker model. P logical workers execute each phase's jobs in ascending target order. Job ordinal j uses worker $j\bmod P$ in zero-based wave $\lfloor j/P\rfloor$. All jobs in a wave are conceptually simultaneous, and the next algorithmic phase waits for the entire current phase. A phase with q jobs takes $\lceil q/P\rceil$ modeled rounds.

A Hillis–Steele job performs one ordered combine. An up-sweep job combines one endpoint pair; a down-sweep job reads a pair and writes its two child prefixes using one combine. Root reset is one assignment job with zero combines. Inclusive conversion has one combine job per valid element. Each such job costs one abstract round unit when assigned a worker; this is a defined constant-cost model, not an instruction-by-instruction hardware model. Input copies, untouched carries, snapshot recording, and barrier overhead are outside the printed round count.

The simulator executes these jobs sequentially in the host language and prints phase-end arrays plus each job's lane/wave. It creates no threads. Results do not depend on P, while the modeled number of waves does.

How Parallel Prefix / Scan Works, Step by Step

Hillis–Steele, producing inclusive scan. Start from a copy of the input. For each doubling stride, preserve a stable old array and create a separate destination. Targets from d through $n-1$ combine old[i-d] before old[i]; lower positions carry unchanged. Assign jobs to lane/wave by their ordinal, finish the level, record the destination, and double d. Even if a lane receives several targets in different waves, every target in that level reads old data.

Blelloch up-sweep. Pad to m identities, then visit strides $d=1,2,4,\ldots,m/2$. For each interval of length $2d$, let r be its right endpoint and $l=r-d$ its left-half endpoint. Update $$ \mathrm{tree}[r]\leftarrow\mathrm{tree}[l]\oplus\mathrm{tree}[r]. $$ Right endpoints are $2d-1,4d-1,\ldots,m-1$. Each endpoint pair is disjoint from the other pairs at that level, so jobs can update those slots in place once the previous level has finished. At the top, tree[m-1] holds the complete reduction.

Root reset. Replace the total at the root with e. This means no elements precede the whole input interval. Forgetting this reset propagates prefixes that already include the total, corrupting exclusive scan. The demonstration records the reset as one assignment phase and one modeled round, with zero combine operations.

Blelloch down-sweep. Visit strides $d=m/2,m/4,\ldots,1$ over the same endpoint pairs. Before writing, save the left interval reduction $L=\mathrm{tree}[l]$ and incoming prefix $P_0=\mathrm{tree}[r]$. Then write $$ \mathrm{tree}[l]\leftarrow P_0,\qquad \mathrm{tree}[r]\leftarrow P_0\oplus L. $$ The left child receives the prefix before the parent interval. The right child receives that prefix followed by the entire left interval. Saving L avoids reading a value already overwritten by the first assignment. Keeping $P_0$ before L is essential for noncommutative operators.

After the final down level, the first n entries are exclusive prefixes. The program builds inclusive output with $x_i\oplus a_i$ for every valid i. Original input values remain available for this conversion. The diagram follows the first down-sweep pair of the worked example, showing both old values and both new prefixes.

Empty input skips both algorithms and every phase. For a singleton, Hillis–Steele returns its input with no combine level. Blelloch has no up/down level, resets its single root to e, then optionally combines it with the original input to produce inclusive output. These explicit cases avoid invalid root access or ambiguous padding.

flowchart TD
    A["Validate input and worker count"] --> B{"Input empty?"}
    B -->|Yes| C["Return empty outputs and zero counts"]
    B -->|No| D["Hillis-Steele: copy input; stride one"]
    D --> E{"Stride below input length?"}
    E -->|Yes| F["Read old snapshot; combine active targets; barrier; double stride"]
    F --> E
    E -->|No| G["Save inclusive result; pad original input with identity"]
    G --> H["Blelloch: up-sweep levels with a barrier after each"]
    H --> I["Reset root to identity"]
    I --> J["Down-sweep levels: save left total; propagate ordered prefixes"]
    J --> K["Trim exclusive output to original length"]
    K --> L["Combine exclusive prefix with original element for inclusive output"]
    L --> M["Return both results, phase states and worker counts"]

Worked Example with a Full Trace

Use input $[3,1,7,0,4]$ and two workers. Hillis–Steele operates on five elements directly. Blelloch pads to $m=8$ using $[3,1,7,0,4,0,0,0]$. The sum identity is zero; the last three entries are padding only.

Every Hillis–Steele combine. Stride one performs $3+1=4$, $1+7=8$, $7+0=7$, and $0+4=4$ at targets one through four, giving $[3,4,8,7,4]$. Stride two performs $3+8=11$, $4+7=11$, and $8+4=12$ at targets two through four, giving $[3,4,11,11,12]$. Stride four performs $3+12=15$ at target four, giving $[3,4,11,11,15]$. Each arithmetic expression uses the preceding level, never a value just updated in its own level.

Every Blelloch pair. Up-sweep uses left reduction followed by right reduction. At stride one, endpoint pairs $(0,1),(2,3),(4,5),(6,7)$ compute $3+1=4$, $7+0=7$, $4+0=4$, and $0+0=0$. Stride two combines endpoints $(1,3)$ as $4+7=11$ and $(5,7)$ as $4+0=4$. Stride four combines $(3,7)$ as $11+4=15$.

Reset root slot seven from fifteen to zero. Down-sweep stride four maps old $(L,P_0)=(11,0)$ at slots $(3,7)$ to $(0,11)$. Stride two maps slots $(1,3)$ from $(4,0)$ to $(0,4)$ and slots $(5,7)$ from $(4,11)$ to $(11,15)$. Stride one maps $(0,1)$ from $(3,0)$ to $(0,3)$; $(2,3)$ from $(7,4)$ to $(4,11)$; $(4,5)$ from $(4,11)$ to $(11,15)$; and $(6,7)$ from $(0,15)$ to $(15,15)$.

The full phase-end trace is below. Targets name the right endpoint for tree pair jobs; the companion left endpoint is target minus stride. For root reset, the one target is an assignment, not a combine. Conversion targets are ordinary output indices.

Every worker assignment. Within a phase, jobs are ordered by the target list. With P equal to two, the first two jobs use workers zero and one in wave zero, the next two workers zero and one in wave one, and so on. Thus up stride one's targets $[1,3,5,7]$ map to lane/wave pairs $[(0,0),(1,0),(0,1),(1,1)]$. Up stride two's targets $[3,7]$ run together in wave zero. The root reset and first down level each use only worker zero. The program prints all lane and wave arrays for every phase.

The final padded exclusive state is $[0,3,4,11,11,15,15,15]$. Trim to five values, then convert each exclusive prefix using the original element:

Both inclusive algorithms return $[3,4,11,11,15]$; Blelloch exclusive returns $[0,3,4,11,11]$. Hillis–Steele performs eight combines in five modeled rounds. Blelloch performs seven up-sweep and seven down-sweep combines, plus a root assignment, taking nine exclusive rounds; conversion adds five combines and three rounds, making twelve inclusive rounds.

This small padded example costs more operations and rounds under Blelloch. Linear asymptotic work does not promise fewer operations at every small n, and the inclusive conversion has a separate cost. Compare the same output convention and include padding/setup assumptions when making a performance claim.

Implementation in Python, C++, and Java

Each panel is a complete program with both scans, deterministic lane/wave assignments, full phase states, and the same standard output. Run Python directly, compile C++ as C++17, or save the Java panel as Main.java and run Main. These are worker simulations, not threaded implementations or GPU timing experiments.

Numeric entry point. Call sum_scan(values,processors) in Python/C++, or sumScan(values,processors) in Java. Require $0\le n\le256$, integer inputs in $[-10^6,10^6]$, and $1\le P\le64$. Invalid worker counts give PROCESSORS, excessive length SIZE, and invalid numeric values VALUE, checked in that order. Python requires exact integers and rejects booleans/floats/strings. C++ uses long long values; Java uses long. Maximum sum magnitude is at most $256\cdot10^6$, safely within their 64-bit arithmetic. Non-null flat containers are expected; inputs remain unchanged.

Generic core. scan_core / scanCore takes values, an identity, a pure binary combine function, and P. It validates the worker/count limits and assumes the operator is associative with that two-sided identity. Stored values are treated as immutable and the combine function must not mutate them. The generic contract cannot prove these algebraic requirements at runtime. Tests use strings with concatenation as well as integer sums. String concatenation has variable cost, so its correctness tests do not imply constant-time combine performance.

Result fields. The result contains padded length, Hillis–Steele inclusive output, Blelloch exclusive output, converted Blelloch inclusive output, all stages, combine counts, and modeled round totals. Each stage stores its phase, stride, target list, worker lanes, zero-based wave numbers, number of rounds, combine count, and phase-end array. A root phase has one job and zero combines; its stride field is zero. The final inclusive phase also has stride zero, with one combine per valid index.

Hillis–Steele keeps old values stable while creating a new level. Blelloch updates tree endpoint pairs in place because different jobs in the same level use disjoint slots. Each down job saves both old values before writing. Phase snapshots are copies, so later changes cannot rewrite earlier trace records.

Empty input has padded length zero, empty outputs/stages, and all counts zero. A singleton has no Hillis–Steele combine and zero Hillis arithmetic rounds; its output still requires ordinary host validation/copy work. Blelloch's singleton path resets its root and converts if inclusive output is requested, giving one exclusive and two inclusive modeled rounds. These round counts intentionally exclude input-copy and reporting costs.

from dataclasses import dataclass


@dataclass
class Stage:
    phase: str
    stride: int
    targets: list
    lanes: list
    waves: list
    rounds: int
    combines: int
    state: list


@dataclass
class Result:
    padded: int
    hillis: list
    exclusive: list
    inclusive: list
    stages: list
    hs_ops: int = 0
    tree_ops: int = 0
    convert_ops: int = 0
    hs_rounds: int = 0
    exclusive_rounds: int = 0
    inclusive_rounds: int = 0


def scan_core(values, identity, combine, processors):
    if type(processors) is not int or not 1 <= processors <= 64:
        raise ValueError("PROCESSORS")
    n = len(values)
    if n > 256:
        raise ValueError("SIZE")
    r = Result(0, [], [], [], [])
    if n == 0:
        return r

    def record(phase, stride, targets, state, combines):
        rounds = (len(targets) + processors - 1) // processors
        r.stages.append(Stage(phase, stride, targets,
                              [j % processors for j in range(len(targets))],
                              [j // processors for j in range(len(targets))], rounds, combines, state.copy()))
        return rounds

    current = values.copy()
    stride = 1
    while stride < n:
        old = current
        current = old.copy()
        targets = list(range(stride, n))
        for i in targets:
            current[i] = combine(old[i - stride], old[i])
        r.hs_ops += len(targets)
        r.hs_rounds += record("hs", stride, targets, current, len(targets))
        stride *= 2
    r.hillis = current
    size = 1
    while size < n:
        size *= 2
    r.padded = size
    tree = values.copy() + [identity] * (size - n)
    stride = 1
    while stride < size:
        targets = list(range(2 * stride - 1, size, 2 * stride))
        for right in targets:
            tree[right] = combine(tree[right - stride], tree[right])
        r.tree_ops += len(targets)
        r.exclusive_rounds += record("up", stride, targets, tree, len(targets))
        stride *= 2
    tree[-1] = identity
    r.exclusive_rounds += record("root", 0, [size - 1], tree, 0)
    stride = size // 2
    while stride:
        targets = list(range(2 * stride - 1, size, 2 * stride))
        for right in targets:
            left = right - stride
            left_total, prefix = tree[left], tree[right]
            tree[left] = prefix
            tree[right] = combine(prefix, left_total)
        r.tree_ops += len(targets)
        r.exclusive_rounds += record("down", stride, targets, tree, len(targets))
        stride //= 2
    r.exclusive = tree[:n]
    r.inclusive = [combine(r.exclusive[i], values[i]) for i in range(n)]
    r.convert_ops = n
    r.inclusive_rounds = r.exclusive_rounds + record("inclusive", 0, list(range(n)), r.inclusive, n)
    return r


def sum_scan(values, processors):
    if type(processors) is not int or not 1 <= processors <= 64:
        raise ValueError("PROCESSORS")
    if len(values) > 256:
        raise ValueError("SIZE")
    if any(type(v) is not int or not -1000000 <= v <= 1000000 for v in values):
        raise ValueError("VALUE")
    return scan_core(values, 0, lambda left, right: left + right, processors)


def compact(values):
    return '[' + ','.join(map(str, values)) + ']'


def main():
    values, processors = [3, 1, 7, 0, 4], 2
    r = sum_scan(values, processors)
    assert r.hillis == r.inclusive
    print(f"input={compact(values)} processors={processors} padded={r.padded}")
    for s in r.stages:
        print(f"{s.phase} stride={s.stride} targets={compact(s.targets)} lanes={compact(s.lanes)} waves={compact(s.waves)} "
              f"rounds={s.rounds} combines={s.combines} state={compact(s.state)}")
    print("hillis_inclusive=" + compact(r.hillis))
    print("blelloch_exclusive=" + compact(r.exclusive))
    print("blelloch_inclusive=" + compact(r.inclusive))
    print(f"combines hillis={r.hs_ops} blelloch_tree={r.tree_ops} inclusive_conversion={r.convert_ops}")
    print(f"model_rounds hillis={r.hs_rounds} blelloch_exclusive={r.exclusive_rounds} blelloch_inclusive={r.inclusive_rounds}")


if __name__ == "__main__":
    main()

#include <algorithm>
#include <cassert>
#include <iostream>
#include <stdexcept>
#include <string>
#include <utility>
#include <vector>
using namespace std;
using Word=long long;
template<class T>struct Stage {
    string phase;
    int stride;
    vector<int> targets,lanes,waves;
    int rounds,combines;
    vector<T> state;
};
template<class T>struct Result {
    int padded=0;
    vector<T> hillis,exclusive,inclusive;
    vector<Stage<T>> stages;
    int hs_ops=0,tree_ops=0,convert_ops=0,hs_rounds=0,exclusive_rounds=0,inclusive_rounds=0;
};
template<class T,class Combine>Result<T> scan_core(const vector<T>& values,const T& identity,Combine combine,int processors) {
    if(processors<1 || processors>64)throw invalid_argument("PROCESSORS");
    if(values.size()>256)throw invalid_argument("SIZE");
    int n=static_cast<int>(values.size());Result<T> r;if(n==0)return r;
    auto record=[&](const string& phase,int stride,const vector<int>& targets,const vector<T>& state,int combines) {
        vector<int> lanes,waves;
        for(int j=0;j<static_cast<int>(targets.size());++j) { lanes.push_back(j%processors);waves.push_back(j/processors); }
        int rounds=(static_cast<int>(targets.size())+processors-1)/processors;
        r.stages.push_back(Stage<T>{phase,stride,targets,move(lanes),move(waves),rounds,combines,state});return rounds;
    };
    vector<T> current=values;
    for(int stride=1;stride<n;stride*=2) {
        vector<T> old=current;vector<int> targets;
        for(int i=stride;i<n;++i) { targets.push_back(i);current[i]=combine(old[i-stride],old[i]); }
        int count=static_cast<int>(targets.size());r.hs_ops+=count;r.hs_rounds+=record("hs",stride,targets,current,count);
    }
    r.hillis=current;int size=1;while(size<n)size*=2;r.padded=size;
    vector<T> tree=values;tree.resize(size,identity);
    for(int stride=1;stride<size;stride*=2) {
        vector<int> targets;
        for(int right=2*stride-1;right<size;right+=2*stride) { targets.push_back(right);tree[right]=combine(tree[right-stride],tree[right]); }
        int count=static_cast<int>(targets.size());r.tree_ops+=count;r.exclusive_rounds+=record("up",stride,targets,tree,count);
    }
    tree.back()=identity;r.exclusive_rounds+=record("root",0,vector<int>{size-1},tree,0);
    for(int stride=size/2;stride>0;stride/=2) {
        vector<int> targets;
        for(int right=2*stride-1;right<size;right+=2*stride) {
            targets.push_back(right);int left=right-stride;T left_total=tree[left],prefix=tree[right];
            tree[left]=prefix;tree[right]=combine(prefix,left_total);
        }
        int count=static_cast<int>(targets.size());r.tree_ops+=count;r.exclusive_rounds+=record("down",stride,targets,tree,count);
    }
    r.exclusive.assign(tree.begin(),tree.begin()+n);vector<int> targets;
    for(int i=0;i<n;++i) { targets.push_back(i);r.inclusive.push_back(combine(r.exclusive[i],values[i])); }
    r.convert_ops=n;r.inclusive_rounds=r.exclusive_rounds+record("inclusive",0,targets,r.inclusive,n);return r;
}
Result<Word> sum_scan(const vector<Word>& values,int processors) {
    if(processors<1 || processors>64)throw invalid_argument("PROCESSORS");
    if(values.size()>256)throw invalid_argument("SIZE");
    for(Word v:values)if(v < -1000000 || v > 1000000)throw invalid_argument("VALUE");
    return scan_core(values,Word{0},[](Word left,Word right){return left+right;},processors);
}
template<class T>string compact(const vector<T>& values) {
    string s="[";for(size_t i=0;i<values.size();++i) { if(i)s+=',';s+=to_string(values[i]); }return s+"]";
}
int main() {
    vector<Word> values{3,1,7,0,4};int processors=2;Result<Word> r=sum_scan(values,processors);assert(r.hillis==r.inclusive);
    cout<<"input="<<compact(values)<<" processors="<<processors<<" padded="<<r.padded<<'\n';
    for(const auto& s:r.stages)cout<<s.phase<<" stride="<<s.stride<<" targets="<<compact(s.targets)<<" lanes="<<compact(s.lanes)<<" waves="<<compact(s.waves)
        <<" rounds="<<s.rounds<<" combines="<<s.combines<<" state="<<compact(s.state)<<'\n';
    cout<<"hillis_inclusive="<<compact(r.hillis)<<"\nblelloch_exclusive="<<compact(r.exclusive)<<"\nblelloch_inclusive="<<compact(r.inclusive)
        <<"\ncombines hillis="<<r.hs_ops<<" blelloch_tree="<<r.tree_ops<<" inclusive_conversion="<<r.convert_ops
        <<"\nmodel_rounds hillis="<<r.hs_rounds<<" blelloch_exclusive="<<r.exclusive_rounds<<" blelloch_inclusive="<<r.inclusive_rounds<<'\n';
    return 0;
}

import java.util.ArrayList;
import java.util.Arrays;
import java.util.Collections;
import java.util.List;
import java.util.function.BinaryOperator;

public class Main {
    static final class Stage<T> {
        final String phase;
        final int stride,rounds,combines;
        final List<Integer> targets,lanes=new ArrayList<>(),waves=new ArrayList<>();
        final List<T> state;
        Stage(String phase,int stride,List<Integer> targets,List<T> state,int combines,int processors) {
            this.phase=phase;this.stride=stride;this.targets=new ArrayList<>(targets);this.state=new ArrayList<>(state);this.combines=combines;
            rounds=(targets.size()+processors-1)/processors;
            for(int j=0;j<targets.size();++j) { lanes.add(j%processors);waves.add(j/processors); }
        }
    }
    static final class Result<T> {
        int padded,hsOps,treeOps,convertOps,hsRounds,exclusiveRounds,inclusiveRounds;
        List<T> hillis=new ArrayList<>(),exclusive=new ArrayList<>(),inclusive=new ArrayList<>();
        final List<Stage<T>> stages=new ArrayList<>();
        int record(String phase,int stride,List<Integer> targets,List<T> state,int combines,int processors) {
            Stage<T> s=new Stage<>(phase,stride,targets,state,combines,processors);stages.add(s);return s.rounds;
        }
    }
    static <T> Result<T> scanCore(List<T> values,T identity,BinaryOperator<T> combine,int processors) {
        if(processors<1 || processors>64)throw new IllegalArgumentException("PROCESSORS");
        if(values.size()>256)throw new IllegalArgumentException("SIZE");
        int n=values.size();Result<T> r=new Result<>();if(n==0)return r;
        List<T> current=new ArrayList<>(values);
        for(int stride=1;stride<n;stride*=2) {
            List<T> old=new ArrayList<>(current);List<Integer> targets=new ArrayList<>();
            for(int i=stride;i<n;++i) { targets.add(i);current.set(i,combine.apply(old.get(i-stride),old.get(i))); }
            r.hsOps+=targets.size();r.hsRounds+=r.record("hs",stride,targets,current,targets.size(),processors);
        }
        r.hillis=current;int size=1;while(size<n)size*=2;r.padded=size;
        List<T> tree=new ArrayList<>(values);tree.addAll(Collections.nCopies(size-n,identity));
        for(int stride=1;stride<size;stride*=2) {
            List<Integer> targets=new ArrayList<>();
            for(int right=2*stride-1;right<size;right+=2*stride) { targets.add(right);tree.set(right,combine.apply(tree.get(right-stride),tree.get(right))); }
            r.treeOps+=targets.size();r.exclusiveRounds+=r.record("up",stride,targets,tree,targets.size(),processors);
        }
        tree.set(size-1,identity);r.exclusiveRounds+=r.record("root",0,Arrays.asList(size-1),tree,0,processors);
        for(int stride=size/2;stride>0;stride/=2) {
            List<Integer> targets=new ArrayList<>();
            for(int right=2*stride-1;right<size;right+=2*stride) {
                targets.add(right);int left=right-stride;T leftTotal=tree.get(left),prefix=tree.get(right);
                tree.set(left,prefix);tree.set(right,combine.apply(prefix,leftTotal));
            }
            r.treeOps+=targets.size();r.exclusiveRounds+=r.record("down",stride,targets,tree,targets.size(),processors);
        }
        r.exclusive=new ArrayList<>(tree.subList(0,n));List<Integer> targets=new ArrayList<>();
        for(int i=0;i<n;++i) { targets.add(i);r.inclusive.add(combine.apply(r.exclusive.get(i),values.get(i))); }
        r.convertOps=n;r.inclusiveRounds=r.exclusiveRounds+r.record("inclusive",0,targets,r.inclusive,n,processors);return r;
    }
    static Result<Long> sumScan(long[] values,int processors) {
        if(processors<1 || processors>64)throw new IllegalArgumentException("PROCESSORS");
        if(values.length>256)throw new IllegalArgumentException("SIZE");
        List<Long> list=new ArrayList<>();
        for(long v:values) { if(v < -1000000 || v > 1000000)throw new IllegalArgumentException("VALUE");list.add(v); }
        return scanCore(list,0L,(left,right)->left+right,processors);
    }
    static String compact(List<?> values) { return values.toString().replace(" ",""); }
    public static void main(String[]args) {
        long[] values={3,1,7,0,4};int processors=2;Result<Long> r=sumScan(values,processors);
        if(!r.hillis.equals(r.inclusive))throw new AssertionError("scan");
        System.out.println("input="+Arrays.toString(values).replace(" ","")+" processors="+processors+" padded="+r.padded);
        for(Stage<Long> s:r.stages)System.out.println(s.phase+" stride="+s.stride+" targets="+compact(s.targets)+" lanes="+compact(s.lanes)+" waves="+compact(s.waves)
            +" rounds="+s.rounds+" combines="+s.combines+" state="+compact(s.state));
        System.out.println("hillis_inclusive="+compact(r.hillis));System.out.println("blelloch_exclusive="+compact(r.exclusive));System.out.println("blelloch_inclusive="+compact(r.inclusive));
        System.out.println("combines hillis="+r.hsOps+" blelloch_tree="+r.treeOps+" inclusive_conversion="+r.convertOps);
        System.out.println("model_rounds hillis="+r.hsRounds+" blelloch_exclusive="+r.exclusiveRounds+" blelloch_inclusive="+r.inclusiveRounds);
    }
}

Expected standard output, identical in all three languages:

• input=[3,1,7,0,4] processors=2 padded=8
• hs stride=1 targets=[1,2,3,4] lanes=[0,1,0,1] waves=[0,0,1,1] rounds=2 combines=4 state=[3,4,8,7,4]
• hs stride=2 targets=[2,3,4] lanes=[0,1,0] waves=[0,0,1] rounds=2 combines=3 state=[3,4,11,11,12]
• hs stride=4 targets=[4] lanes=[0] waves=[0] rounds=1 combines=1 state=[3,4,11,11,15]
• up stride=1 targets=[1,3,5,7] lanes=[0,1,0,1] waves=[0,0,1,1] rounds=2 combines=4 state=[3,4,7,7,4,4,0,0]
• up stride=2 targets=[3,7] lanes=[0,1] waves=[0,0] rounds=1 combines=2 state=[3,4,7,11,4,4,0,4]
• up stride=4 targets=[7] lanes=[0] waves=[0] rounds=1 combines=1 state=[3,4,7,11,4,4,0,15]
• root stride=0 targets=[7] lanes=[0] waves=[0] rounds=1 combines=0 state=[3,4,7,11,4,4,0,0]
• down stride=4 targets=[7] lanes=[0] waves=[0] rounds=1 combines=1 state=[3,4,7,0,4,4,0,11]
• down stride=2 targets=[3,7] lanes=[0,1] waves=[0,0] rounds=1 combines=2 state=[3,0,7,4,4,11,0,15]
• down stride=1 targets=[1,3,5,7] lanes=[0,1,0,1] waves=[0,0,1,1] rounds=2 combines=4 state=[0,3,4,11,11,15,15,15]
• inclusive stride=0 targets=[0,1,2,3,4] lanes=[0,1,0,1,0] waves=[0,0,1,1,2] rounds=3 combines=5 state=[3,4,11,11,15]
• hillis_inclusive=[3,4,11,11,15]
• blelloch_exclusive=[0,3,4,11,11]
• blelloch_inclusive=[3,4,11,11,15]
• combines hillis=8 blelloch_tree=14 inclusive_conversion=5
• model_rounds hillis=5 blelloch_exclusive=9 blelloch_inclusive=12

Complexity, Correctness, and Invariants

Hillis–Steele work and span. For positive n, let $k=\lceil\log_2n\rceil$ and $m=2^k$. There are k combine levels. Level $r$ uses stride $2^r$ and $n-2^r$ combines, so the exact combine count is $$ H(n)=\sum_{r=0}^{k-1}(n-2^r)=nk-(2^k-1). $$ For n one, k and H are zero. For growing n, work is $\Theta(n\log n)$ and dependency span $\Theta(\log n)$ under constant-time operations. A sequential scan uses linear work, so Hillis–Steele is not work-efficient asymptotically.

Blelloch work and span. A padded tree of m leaves has $m-1$ internal nodes. Up-sweep and down-sweep each use one combine per internal node: $$ B_{\mathrm{tree}}(n)=2(m-1). $$ Root reset adds one assignment, not a combine. The demonstrated inclusive conversion adds n combines, giving $2(m-1)+n$ binary operations for converted inclusive output. Padding and output copies are linear because $n\le m\lt2n$ for every positive n, including powers of two. The algorithm's total work is $\Theta(n)$ and span $\Theta(1+\log n)$ when ordinary initialization/output work is included. There are $2k+1$ tree/root phases for exclusive output, then one parallel conversion phase for inclusive output.

Exact finite-worker model. The printed round totals are $$ R_{\mathrm{HS}}(n,P)=\sum_{r=0}^{k-1}\left\lceil\frac{n-2^r}{P}\right\rceil, $$ $$ R_{\mathrm{exclusive}}(n,P)=1+2\sum_{r=1}^{k}\left\lceil\frac{m/2^r}{P}\right\rceil, \qquad R_{\mathrm{inclusive}}=R_{\mathrm{exclusive}}+\left\lceil\frac nP\right\rceil. $$ These formulas apply to positive n; the empty case has zero rounds. For the sample, k three and m eight give five Hillis rounds, nine Blelloch exclusive rounds, and twelve converted inclusive rounds. With enough workers, the sums reduce to k, $2k+1$, and $2k+2$ respectively. These are defined job-stage costs with phase barriers, not measured hardware cycles.

For $n\ge2$, the phase model scales as $\Theta(n\log n/P+\log n)$ for Hillis–Steele and $\Theta(n/P+\log n)$ for Blelloch. This connects to work/span analysis: enough workers expose logarithmic depth, while finite processors must also perform the total work. Both algorithms are data-oblivious within the declared shape: input values do not change targets, phase counts, or scheduling. Thus best/average/worst combine counts for fixed n are identical; there is no probability assumption about data values.

Host implementation overhead. The demonstration deliberately copies every phase-end array for its trace. Hillis–Steele also maintains a stable source snapshot at each level. Recording the tree's logarithmically many full m-entry states adds $O(m(1+\log m))$ host time and retained storage, even though its arithmetic tree algorithm has linear combine work. The complete simulator therefore uses $O(n(1+\log n))$ host work and trace storage for constant-size values, independent of the modeled P-worker speed. Removing saved history permits $O(n)$ working array space; performing only the disjoint Blelloch endpoint updates/padding/conversion gives linear host arithmetic/data work. Returning all prefixes already requires $\Omega(n)$ output storage or writes.

Hillis–Steele invariant. Before stride d, index i stores the ordered reduction of the last at most d input elements ending at i. The old interval for $i-d$ ends immediately before the old interval for i. Combining earlier interval first yields the ordered last at most $2d$ elements. Lower indices retain the complete prefix they already cover. Induction over doubled strides proves every final inclusive prefix. Stable snapshots preserve this induction even when several target jobs are assigned to one worker over different waves.

Blelloch up-sweep invariant. At each interval's right endpoint, the stored value becomes the ordered reduction of that interval. Its left and right child reductions are adjacent and disjoint; associativity permits combining them without changing order. Identity padding preserves original reductions. The root therefore holds the total.

Down-sweep invariant. Each active parent endpoint holds the reduction strictly before its interval, while its left endpoint still holds the left child's total. The left child receives the same preceding reduction; the right child receives that reduction followed by the left total. Saving the old left total before assignment and preserving operand order maintain this property at the next level. Resetting the root supplies the empty preceding reduction e. At leaves, the invariant is exactly exclusive scan, and $x_i\oplus a_i$ yields inclusive scan.

The algorithms never rely on commutativity or inverse operations. Incorrect identities, nonassociative arithmetic, missing barriers, or mutated old values invalidate the proof assumptions rather than merely changing performance.

Variants, Pitfalls, and Optimizations

Stable Hillis–Steele reads are essential. In a left-to-right in-place loop, an index can read its neighbor's newly updated value and double-count input. The diagram's four-element example reaches wrong result $[3,4,14,15]$. A descending sequential loop can sometimes preserve the old-read property for this particular indexing, but that is not a substitute for race-free simultaneous parallel levels. Use distinct level buffers or an implementation that explicitly guarantees the same data dependencies.

Blelloch's in-place rule is different. Within a single tree level, endpoint pairs are disjoint, so those jobs can safely read/write their own pair. Across levels, parents depend on completed child totals or prefixes. Save the left total before replacing it during down-sweep, reset the root, and synchronize the phase transitions. A barrier among only one subset of GPU workers cannot act as a global barrier for unrelated blocks.

Preserve noncommutative order. Correct down-sweep writes prefix ⊕ left_total to the right child. Addition hides operand reversal; concatenation exposes it. Likewise, Hillis–Steele puts the earlier interval before the later one, and hierarchical scans prepend a chunk's preceding prefix to its local prefix. General scan does not grant permission to sort, reverse, or swap operands.

Use the correct identity and output convention. Addition pads with zero, multiplication with one, maximum with an appropriate negative-infinity identity, and concatenation with the empty string. Padding with an ordinary value that is not the identity changes the tree totals. Exclusive output starts with e and excludes the current element; inclusive includes it. Shifting with the wrong direction or forgetting to trim padding produces a different result shape.

Inclusive conversion can use a saved total. The demonstration combines $x_i\oplus a_i$ at every position because the rule is easy to verify. Another valid conversion shifts exclusive results left and appends the original total saved before root reset. This can avoid n conversion combines, though moving/materializing n outputs still has a cost. Retain that saved total explicitly; the reset root no longer contains it. No subtraction or inverse operation is needed.

Small inputs and padding can favor simpler code. The five-element example has eight Hillis combines versus fourteen Blelloch tree combines plus conversion. Near a power-of-two boundary, m can approach twice n. Work-efficient asymptotics remain valuable for large n, but do not determine small-case constants, memory traffic, or barrier costs. Compare the same requested output with initialization and conversion included consistently.

Floating point is not exactly associative. Parallel regrouping may produce a different rounded result from a sequential accumulator. A fixed tree can give repeatable results for a fixed implementation, but repeatability is not equality to every other tree. Choose a numerical error/tolerance model or a reproducible summation method when exact equality is required. The supplied integer bounds avoid overflow and rounding so exact comparisons are appropriate.

Generic operations have their own cost. The work/span table assumes fixed-size values and constant-time combines. Concatenating growing strings can copy increasingly large payloads; matrix products or symbolic expressions have different costs. The noncommutative string suite checks operand order, not the unit-cost performance assumption. Integer addition also needs a sufficiently wide type when input bounds are relaxed.

Hierarchical scan handles large arrays. Scan each chunk locally and retain each chunk total. Exclusive-scan those totals, then combine each preceding-chunk prefix before the local prefix. This composes ordered scans without requiring every element to share one physical execution group. Tail chunks need identities or valid-length handling, and inter-chunk stages need a real global ordering mechanism.

Segmented scans add boundaries. To reset prefixes at marked starts, lift values with head flags and use an associative operator that either starts a new aggregate at a right-hand head or extends the preceding aggregate. An ordinary scan ignoring those flags crosses segment boundaries. Segmentation is an extension of the algebra, not an arbitrary reset inserted into simultaneous jobs.

Application preconditions matter. Signed prefix sums are supported, but offset allocation and cumulative distributions generally require nonnegative lengths/weights. Negative inputs can make prefixes decrease. Stable compaction also needs each selected element to use its exclusive flag count as its destination; inclusive counts would shift every selected destination by one.

Practice Problems and Real-World Use

Verify both outputs against a direct sequential ordered fold of each prefix. Do not generate the oracle by repeating the same tree updates. For stage checks, derive Hillis intervals directly from original input; derive up-sweep endpoint totals from original interval folds and down-sweep entries from the original prefix before each child interval. This independently checks both arithmetic and untouched entries.

Useful properties include unchanged input, equal lengths, inclusive agreement between algorithms, exclusive first element e, and exclusive output shifted by one matching earlier inclusive outputs. Check all intermediate arrays, target lists, lane/wave mappings, combine counts, and round totals. Include empty/singleton input, negative/zero values, powers of two and adjacent lengths, more workers than jobs, and a single worker. Noncommutative operators reveal order mistakes that sums cannot.

All three implementations passed $6313$ shared inputs per language: $5021$ integer-sum cases and $1292$ concatenation cases. There are $6306$ valid cases and $7$ rejected numeric inputs. The suite includes exhaustive short integer arrays, every numeric length through $256$ across several worker counts, randomized arrays, empty strings/multi-character concatenations, padding boundaries, invalid limits/values, and $8$ additional Python type rejections. Direct interval/prefix-fold oracles match every output and phase state, all worker lane/wave assignments, operation counts, and model rounds. An instrumented generic combine checks actual call counts; a deliberate left-to-right snapshot bug reproduces the diagram's wrong result. C++ ran with address and undefined-behavior sanitizers. All standalone samples print byte-identical output.

For stable compaction, use flags $[0,1,1,0,1]$. Their exclusive sums are $[0,0,1,2,2]$, so selected elements at original indices one, two, four go to destinations zero, one, two. The final count is the last exclusive prefix plus the last flag, equal to three. This is a practical reason the exclusive convention is useful.

For allocation, the worked input's offsets define regions $[0,3)$, $[3,4)$, $[4,11)$, $[11,11)$, and $[11,15)$. The zero-length region is empty and legitimately shares an offset with the next region. The same pattern assigns positions in sparse-array construction, variable-length output, and radix-sort scatter steps when their own ordering rules are respected.