Add up ten million copies of 0.1 in float32, one after another, and you get 1,087,937, not 1,000,000. Add the same numbers in a different order, pairing them up in a tree, and you get 1,000,000 to every digit float32 can show. Neither number is a bug. Floating-point addition rounds after every step, so the order of additions changes the answer, and the tree order, which happens to be the natural order on a GPU, is more accurate.

Every kernel of §11.3 gave each thread its own output. A vector add, a matrix scale, and a transpose never ask two threads to agree on anything. A sum does. A billion inputs must become one number, so threads have to combine their work, first through the registers of a warp, then through the shared memory and barriers of a block, and finally across blocks that cannot wait for each other. This is where the tools of §11.3 start working together, and where the order of operations, which never mattered for a vector add, starts to change the answer.

Summing is the simplest example of a reduction, combining many values into one with an associative operation, and its sibling the scan (prefix sum) produces every partial result along the way. Reductions are inside every softmax, LayerNorm, loss, and gradient norm in a model, and inside every all-reduce across GPUs. Scans are less visible but just as basic. They compute the offsets when tokens are routed to experts, the cumulative probabilities of top-p sampling, the positions in packed sequences, and, most strikingly, they turn a linear recurrence, a computation that looks hopelessly sequential, into one that runs in logarithmic depth. That last trick is what makes modern recurrent architectures trainable on GPUs.

This section builds both from the ground up. We derive tree reductions and their cost, see why the way a block pairs up its elements decides whether it diverges and conflicts (§11.1 and §11.3 again), implement warp shuffles, quantify why parallel sums are more accurate and how to make them deterministic, merge means, variances, and softmax statistics, handle segmented reductions, and then derive the Blelloch work-efficient scan with a proof of correctness. It ends with scans over operators other than addition, including the one that parallelizes $h_t = a_t h_{t-1} + b_t$.

Tree reductions

To sum $n$ numbers sequentially takes $n - 1$ additions, each depending on the previous one, so even with infinitely many processors it takes $n - 1$ steps. Pair them up instead. In the first step, $n/2$ processors each add one pair. In the next, $n/4$ add pairs of those results, and so on, until one value remains after $\lceil\log_2 n\rceil$ steps.

Two quantities describe a parallel algorithm. Its work $W(n)$ is the total number of operations, and its depth $D(n)$ is the length of the longest chain of operations that must happen one after another. The tree reduction has $W(n) = n - 1$, the same as the sequential sum, and $D(n) = \lceil \log_2 n\rceil$. Brent's theorem says that with $p$ processors an algorithm can run in time at most $W/p + D$, so the tree finishes in about $n/p + \log_2 n$ steps, essentially perfect speedup until $p$ approaches $n/\log_2 n$.

In plain words, a tree reduction does no more work than a loop, and its depth is only logarithmic, so it uses as many processors as you have.

The cell runs the tree the way a GPU block does, with sequential addressing. At each step the first half of the active elements adds the second half. The next subsection shows why this particular pairing is the right one.

Pairing elements in shared memory

Inside a block, the tree runs in shared memory, with one thread per pair and a barrier between levels, since each level reads sums the previous level wrote. Which thread adds which pair looks like a detail. It is not, and the three obvious choices show why (Harris, 2007).

• Interleaved addressing. At the level with stride $s = 1, 2, 4, \dots$, thread $t$ is active if $t$ is a multiple of $2s$, and adds x[t + s] into x[t]. Active threads are scattered, so every warp contains both active and idle threads. Every warp diverges (§11.1), and most lanes do nothing while still occupying the warp.
• Strided indexing. Make the active threads contiguous instead, thread $t < n/(2s)$ adding x[2st + s] into x[2st]. Now whole warps go idle together, but lane $t$ reads word $2st$, a stride of $2s$ words, and by the bank rule of §11.3 that is a $2s$-way conflict, up to 32-way.
• Sequential addressing. Thread $t < h$ adds x[t + h] into x[t], with $h = n/2, n/4, \dots, 1$. Active threads are contiguous, so warps never diverge until fewer than 32 threads remain, and consecutive lanes read consecutive words, so there are no conflicts.

__global__ void block_sum_shared(const float* x, float* block_sums, int64_t n) {
    extern __shared__ float sdata[];                   // blockDim.x floats, a power of two
    int t = threadIdx.x;
    int64_t i = (int64_t)blockIdx.x * blockDim.x + t;
    sdata[t] = (i < n) ? x[i] : 0.0f;                  // pad with the identity
    __syncthreads();
    for (int h = blockDim.x / 2; h > 0; h >>= 1) {     // sequential addressing
        if (t < h) sdata[t] += sdata[t + h];           // threads 0..h-1, consecutive words
        __syncthreads();                               // outside the if, so every thread arrives
    }
    if (t == 0) block_sums[blockIdx.x] = sdata[0];
}
// block_sum_shared<<<blocks, 256, 256 * sizeof(float)>>>(x, block_sums, n);

The barrier sits outside the if, as §11.3 required, so idle threads still reach it. The three pairings are Harris's kernels 1 to 3, and his sequence continues. Kernel 4 adds two elements as each thread loads, so no thread is idle at the first level. Kernel 5 unrolls the last six levels, where only one warp is active, and drops their barriers. (That relied on the lanes of a warp running in lockstep, which is no longer guaranteed since Volta, so modern code finishes with the warp shuffles below, or at least __syncwarp().) Kernel 6 unrolls the whole loop for a block size known at compile time. Kernel 7 has each thread sum many elements before the tree starts. On a G80 those steps took a sum of 4 million integers from about 2 to about 63 GB/s of the chip's 86 GB/s, thirty times faster with the same arithmetic.

import numpy as np

def pairing_costs(n, scheme):
    levels = int(np.log2(n))
    return [(n // 32, 0, 1)] * levels                # placeholder: every level looks perfect

From a block to the whole GPU

A real kernel adds three refinements. Harris (2007) walked a reduction kernel through a sequence of such optimizations, each a large speedup over the last, and the modern version looks like this.

Each thread first sums many elements by itself. A grid-stride loop (§11.3) has each thread accumulate a private partial sum over its share of the input, sequentially, in a register. This keeps the work at $n$ while cutting the tree to the number of threads, and it makes the loads perfectly coalesced. It is Brent's theorem applied on purpose.

Warps reduce with shuffles. Threads of a warp can read each other's registers directly with __shfl_down_sync, without shared memory or barriers. Five shuffle-and-add steps (offsets 16, 8, 4, 2, 1) reduce 32 values to one.

Blocks combine through a second pass or atomics. Lane 0 of each warp writes its result to shared memory, one warp reduces those, and each block then either writes its partial sum to an array for a second, tiny kernel, or adds it to the output with atomicAdd. Blocks cannot wait for each other inside a kernel (§11.3), which is why the last level needs one of these two.

__inline__ __device__ float warp_reduce_sum(float v) {
    for (int offset = 16; offset > 0; offset >>= 1)
        v += __shfl_down_sync(0xffffffff, v, offset);  // lane i adds the value held by lane i + offset
    return v;                                          // lane 0 now holds the sum of all 32 lanes
}

__global__ void sum_kernel(const float* x, float* block_sums, int64_t n) {
    __shared__ float warp_sums[32];                    // one slot per warp, up to 1,024 threads
    float v = 0.0f;
    for (int64_t i = (int64_t)blockIdx.x * blockDim.x + threadIdx.x; i < n;
         i += (int64_t)gridDim.x * blockDim.x)
        v += x[i];                                     // private, sequential, coalesced
    v = warp_reduce_sum(v);
    int lane = threadIdx.x % 32, warp = threadIdx.x / 32;
    if (lane == 0) warp_sums[warp] = v;
    __syncthreads();
    if (warp == 0) {
        int n_warps = (blockDim.x + 31) / 32;
        v = (lane < n_warps) ? warp_sums[lane] : 0.0f;
        v = warp_reduce_sum(v);
        if (lane == 0) block_sums[blockIdx.x] = v;     // deterministic; a second launch sums these
    }
}

The shuffle's semantics are precise. __shfl_down_sync(mask, v, delta) returns, to lane $i$, the value of v in lane $i + \delta$, and lanes with $i + \delta \ge 32$ get their own value back. After the step with offset 16, lane $i < 16$ holds $v_i + v_{i+16}$, after offset 8 lane $i < 8$ holds the sum of four values, and after offset 1 lane 0 holds all 32. The other lanes hold partial garbage, which is fine. The variant __shfl_xor_sync(mask, v, m) exchanges values between lanes $i$ and $i \oplus m$ (bitwise exclusive or), and five steps with $m = 16, 8, 4, 2, 1$ leave the full sum in every lane, a butterfly all-reduce, useful when every thread needs the result, as in the normalization of a softmax.

How fast can it be? A float32 sum does one addition per 4-byte element, an arithmetic intensity of $1/4$ FLOP per byte, about a thousand times below the H100's ridge point (§11.2). So a reduction is memory-bound, and the only figure of merit is achieved bandwidth. Summing $10^9$ float32 values, 4 GB, takes at least $4\times10^{9}/3.35\times10^{12} \approx 1.2$ ms on an H100. Every refinement above, coalesced grid-stride loads, shuffles instead of barriers, one write per block, exists to bring a reduction close to that bound, which is about the speed of a plain memory copy. In production you rarely write one. CUB's BlockReduce and DeviceReduce, and torch.sum, already do all of this.

import numpy as np

def shfl_down(v, delta):
    return np.roll(v, -delta)                         # placeholder: wraps around, unlike CUDA

def shfl_xor(v, m):
    return v.copy()

Numerically safe parallel sums

Back to the puzzle. Every floating-point addition rounds its result to the nearest representable number, introducing a relative error of at most $u$, the unit roundoff ($u = 2^{-24} \approx 6\times10^{-8}$ in float32, half the machine epsilon of §0.9). In sequential summation, the first element passes through $n - 1$ roundings, and the classical bound (Higham, 1993) is

where $S$ is the exact sum, $\hat S$ the computed one, and the bounds hold to first order in $u$. In a tree each element passes through only $\log_2 n$ roundings.

In plain words, the error of a sum grows with the number of additions each input passes through, which is $n$ in a loop and $\log_2 n$ in a tree.

The sequential bound is also the typical behavior when the partial sum grows large. Once the running total of our 0.1s passes $2^{19} \approx 524{,}000$, neighboring float32 numbers are $2^{-4} = 0.0625$ apart, so each added 0.1, which is 1.6 gaps, is rounded to two gaps, 0.125. That 25 percent excess, repeated over millions of additions, is what pushes the total so far past a million. Kahan summation fixes the sequential loop by carrying a compensation term that recovers the low-order bits each addition loses, with an error bound of about $2u\sum|x_i|$ independent of $n$, at four times the operations. On GPUs the tree, plus float32 accumulation for low-precision inputs, is usually enough.

The tree is several hundred times more accurate than the loop here, as accurate as Kahan's compensated loop. For the low-precision formats of training the lesson is sharper. A bfloat16 number has only 8 bits of significand, so a bfloat16 accumulator stops changing once the total is about 256 times the typical addend. Every reduction in a mixed-precision model (the softmax denominator, the LayerNorm statistics, the loss, gradient norms, the accumulators of matmuls) is therefore done in float32, even when inputs and outputs are 16-bit (§8.5).

Determinism

Because the order matters, a reduction whose order varies gives results that vary. atomicAdd of per-block partials adds them in whatever order blocks happen to finish, so the last bits change from run to run. Changing the number of blocks or threads, as a library might for a different input size, changes the tree. This is why some PyTorch operations are nondeterministic unless you ask for deterministic algorithms (§0.9). A deterministic reduction fixes the tree. Each block writes its partial sum to its own slot, and a second pass combines the slots in a fixed order, at the cost of a second launch.

A subtler version matters for inference. Some kernels choose how to split a reduction based on the batch size, so the same request can produce slightly different logits depending on what else is in the batch. Thinking Machines Lab (He, 2025) traced much of the run-to-run variation of LLM inference servers to exactly this lack of batch invariance, and showed that kernels with fixed reduction orders make sampling at temperature 0 reproducible.

Merging statistics, not just sums

A LayerNorm needs a mean and a variance, and a numerically safe one cannot be computed as $\operatorname{mean}(x^2) - \mu^2$, which cancels catastrophically when the mean $\mu$ is much larger than the standard deviation $\sigma$ (§5.2 showed this and fixed it with Welford's sequential update). To parallelize it, each thread or warp summarizes its chunk as a triple $(n_a, \mu_a, M_a)$, the count, the mean, and $M_a = \sum (x_i - \mu_a)^2$, and two summaries are merged by Chan, Golub & LeVeque's formula. With $n = n_a + n_b$ and $\delta = \mu_b - \mu_a$,

In plain words, the combined mean moves toward the second chunk's mean in proportion to its size, and the combined spread is the two spreads plus a term for how far apart the two means are.

The merge is associative, so it can be applied in a tree like any sum. That is how fused LayerNorm kernels split a row across the threads of a block.

The softmax denominator is a reduction of the same kind. To sum $e^{x_j}$ over a row without overflow, each chunk keeps its maximum $m$ and its sum $\ell = \sum_j e^{x_j - m}$ relative to that maximum, and two chunk summaries combine by

with identity $(-\infty, 0)$. In plain words, both partial sums are rescaled to the larger of the two maxima before they are added, so every exponent stays at most zero. The operator is associative (and commutative), so a row's softmax statistics reduce in a tree, across the lanes of a warp, or as a stream of tiles, like any sum. This is the online softmax. §12.2 derives it and proves the merge exact, §11.6 uses it to stream rows too long for one block, and FlashAttention is built on it.

import numpy as np

def merge(a, b):
    n = a[0] + b[0]
    return n, (a[0] * a[1] + b[0] * b[1]) / max(n, 1), a[2] + b[2]   # placeholder: forgets the delta term

def tree_mean_var(x, leaf=32):
    x = np.asarray(x, dtype=float)

Segmented reductions

Many reductions in ML are not over a whole tensor but over groups of varying size laid end to end. Tokens routed to the same expert in a mixture-of-experts layer (§13.2), documents packed into one sequence (§4.4), the bags of an embedding-bag lookup, the neighbors of each node in a graph network. A standard representation is offsets, as in a sparse CSR matrix. Segment $j$ is values[offsets[j]:offsets[j + 1]], so offsets has one more entry than there are segments, and an empty segment has equal consecutive offsets.

On a GPU there are three strategies. One block per segment is simple but balances badly when segment sizes vary a lot. Sorting by segment and running a reduce-by-key over the flattened array balances perfectly. And a segmented scan uses the scan machinery below with an operator that restarts at each segment boundary. Two tempting NumPy shortcuts are wrong in instructive ways. np.add.reduceat returns a wrong value, not zero, for an empty segment, and differencing a global prefix sum, $c_{o_{j+1}} - c_{o_j}$, cancels catastrophically when an earlier segment is large, the same disease as $\operatorname{mean}(x^2) - \mu^2$.

import numpy as np

def segment_sum(values, offsets):
    values = np.asarray(values, dtype=float)
    return np.add.reduceat(values, np.asarray(offsets)[:-1], axis=0)   # placeholder: wrong for empty segments

Prefix sums

A scan applies an associative operator $\oplus$ with identity $e$ cumulatively. The inclusive scan of $(x_0, \dots, x_{n-1})$ is $y_i = x_0 \oplus \cdots \oplus x_i$, and the exclusive scan is $y_i = x_0 \oplus \cdots \oplus x_{i-1}$, with $y_0 = e$. For addition, these are np.cumsum and np.cumsum shifted right by one with a leading 0. The exclusive form is the one that computes offsets. If $x_i$ is the number of items segment $i$ will write, the exclusive scan says where each segment starts, which is how a GPU compacts an array, sorts by radix, or lays out the tokens bound for each expert.

A sequential scan does $n - 1$ operations in $n - 1$ dependent steps. The question is how to do it in parallel without doing much more work.

Hillis and Steele, fast but wasteful

The first idea (Hillis & Steele, 1986, also known as the Kogge-Stone pattern) is to double the reach at each step. For $d = 1, 2, 4, \dots$, every element adds the element $d$ positions to its left, $y_i \leftarrow y_{i-d} \oplus y_i$ for $i \ge d$. After the step with reach $d$, $y_i$ holds the sum of the $2d$ elements ending at $i$ (or all of them, if there are fewer), so after $\lceil\log_2 n\rceil$ steps it is the inclusive scan.

The depth is $\log_2 n$, as low as it gets, but each step does nearly $n$ additions, so the work is about $n\log_2 n$, a factor of $\log_2 n$ more than the loop. With $10^9$ elements that is 30 times the work, and on a GPU work costs energy and, for large arrays, memory traffic. Hillis-Steele is the right choice when processors are plentiful relative to elements, as within a single warp, where shuffles do each step in one instruction.

__inline__ __device__ float warp_inclusive_scan(float v) {
    int lane = threadIdx.x % 32;
    for (int d = 1; d < 32; d <<= 1) {
        float up = __shfl_up_sync(0xffffffff, v, d);   // lane i receives v from lane i - d
        if (lane >= d) v += up;                        // the first d lanes keep their value
    }
    return v;                                          // lane i now holds x_0 + ... + x_i
}

A block scan is built from this in three steps, the same shape as the block reduction. Each warp scans its 32 values, one warp scans the warp totals, and every thread adds the total of the warps before its own.

Blelloch's work-efficient scan

Blelloch (1990) showed how to scan with $O(n)$ work and $O(\log n)$ depth, using two passes over a balanced binary tree whose leaves are the input. Assume $n$ is a power of two (pad with the identity otherwise), and work in place on an array x.

The up-sweep is the tree reduction of the first part of this section, done in place. For levels $l = 0, 1, \dots, \log_2 n - 1$, with stride $2^{l+1}$, every block of $2^{l+1}$ elements adds its left half's last element into its right half's last element,

Afterward, the last position of every block of size $2^{l+1}$ holds the total of that block, and x[n-1] holds the total of everything.

The down-sweep walks back down. Set x[n-1] to the identity. Then for levels $l = \log_2 n - 1$ down to 0, in every block of $2^{l+1}$ elements, let $L = i + 2^l - 1$ (the end of the left half) and $R = i + 2^{l+1} - 1$ (the end of the right half), and do

When it finishes, x holds the exclusive scan. Step through it on 8 or 16 elements in the widget, and compare with Hillis-Steele's work.

Count the cost. The up-sweep does $n/2 + n/4 + \dots + 1 = n - 1$ operations, and so does the down-sweep, so $W(n) = 2(n - 1)$, about twice the loop and far less than $n\log_2 n$. Each sweep has $\log_2 n$ levels, so $D(n) = 2\log_2 n$.

In plain words, the up-sweep computes the sum of every subtree, and the down-sweep hands each subtree the sum of everything to its left.

import numpy as np

def exclusive_scan(x):
    x = np.asarray(x)
    y = np.cumsum(x) - x                             # placeholder: right answer, sequential, no op count
    return y, 0

Scanning arrays larger than a block

A block can scan a few thousand elements in shared memory. For a large array, the standard scan-then-propagate scheme takes three steps. Each block scans its own chunk and records the chunk's total. A second, small scan computes an exclusive scan of the chunk totals. A third pass adds each chunk's offset to all of its elements. The work stays linear, but the array crosses HBM about four times (read and write in the first pass, read and write in the third), and a scan, like a sum, is memory-bound. Merrill & Garland (2016) reduced this to the minimum, one read and one write, with a single-pass scan in which each block, after computing its local total, looks back at the published totals of its predecessors (decoupled look-back) rather than waiting for a separate kernel. This is how NVIDIA's CUB library scans at close to memory-copy speed.

Scans of other operators, and linear recurrences

Nothing in the Blelloch scan used the fact that $\oplus$ is addition, only that it is associative. Running maximum, running product, and logical or are scans. So is the segmented sum, with elements that are pairs $(f, v)$ of a flag (1 at the start of a segment) and a value, and the operator

which restarts the sum wherever a new segment begins and is associative (check it on three elements). And so are operators on richer states, like the $(m, \ell)$ pairs of the online softmax.

The most consequential example is the linear recurrence

with $a_t$ and $b_t$ given for $t = 1, \dots, T$. It looks inherently sequential, since every state needs the previous one. But each step is an affine map $h \mapsto a_t h + b_t$, and affine maps compose into affine maps. Applying $(a_1, b_1)$ and then $(a_2, b_2)$ gives $h \mapsto a_2(a_1 h + b_1) + b_2$, so define

This operator is associative because function composition is associative, and its identity is $(1, 0)$. The inclusive scan of the pairs $(a_t, b_t)$ under $\bullet$ gives, at position $t$, the composed map from $h_0$ to $h_t$, whose second component is $h_t$ itself (since $h_0 = 0$). So all $T$ states can be computed with depth $O(\log T)$ and work $O(T)$, instead of $T$ dependent steps. The same holds when $h_t$ is a vector and $a_t$ acts elementwise, as in a diagonal state-space model, and, at higher cost, when $a_t$ is a matrix.

In plain words, each step of a linear recurrence is a "multiply then add" rule, two such rules combine into one rule of the same kind, and a scan combines them all in a logarithmic number of rounds.

This is why linear recurrences have come back. An LSTM or GRU (§3.3) cannot be parallelized over time, because its state passes through a nonlinearity at every step. If the recurrence is kept linear in the state, with the gates computed from the input alone, training becomes a scan and parallelizes over the sequence (Martin & Cundy, 2018). Linear attention can be written as such a recurrence with a matrix-valued state, and selective state-space models such as Mamba make $a_t$ and $b_t$ depend on the input token and evaluate the recurrence with a hardware-aware parallel scan that keeps the expanded state in on-chip memory instead of writing it to HBM, the keep-it-on-chip argument behind kernel fusion (§11.7) applied to a scan. Many of these models split the sequence into chunks, process each chunk with dense matrix operations that suit tensor cores, and carry the state from chunk to chunk with exactly the operator above, which is how parallel scans for recurrences are made fast on real GPUs.

In PyTorch

PyTorch exposes the common scans and segmented reductions directly. The block checks a segmented sum done with index_add_, a running maximum with cummax, the linear recurrence by scan against a loop, and shows that the grouping of float32 additions changes the result and that a bfloat16 accumulator stalls.

Read more

• Blelloch (1990). Prefix Sums and Their Applications. The work-efficient scan, its tree derivation, and a catalog of algorithms built from scans.
• Harris, Sengupta & Owens (2007). Parallel Prefix Sum (Scan) with CUDA. GPU Gems 3, chapter 39, the Blelloch scan in CUDA with bank-conflict avoidance and large arrays.
• Harris (2007). Optimizing Parallel Reduction in CUDA. The classic sequence of reduction kernels, each a measured improvement over the last.
• Merrill & Garland (2016). Single-pass Parallel Prefix Scan with Decoupled Look-back. How CUB scans at the speed of a memory copy.
• Higham (1993). The Accuracy of Floating Point Summation. Error bounds for recursive, pairwise, and compensated summation, and when each matters.
• Chan, Golub & LeVeque (1983). Algorithms for Computing the Sample Variance: Analysis and Recommendations. The pairwise merge of means and variances, and why the textbook formula fails.
• Martin & Cundy (2018). Parallelizing Linear Recurrent Neural Nets Over Sequence Length. Linear recurrences as parallel scans, the idea behind the scans of modern state-space models.
• He (2025). Defeating Nondeterminism in LLM Inference. How reduction order and batch-dependent kernels cause nondeterministic inference, and batch-invariant kernels that fix it.

Watch next

▶

▶

▶