Sorting Networks, in depth: comparators, depth and size, the 0-1 principle, Batcher's constructions, SIMD and GPU use

The model: wires, comparators, size and depth

A network has n wires (positions) and a sequence of comparators. A comparator (i, j) with i < j replaces the values at i and j with their minimum and maximum, in that order. Comparators that touch disjoint wires can run at the same time, so a network is usually drawn in layers. Two numbers describe its cost. Size is the total number of comparators, which is the work. Depth is the number of layers, which is the time on hardware with enough parallel units. A network that sorts every input is a sorting network.

The simplest networks come straight from familiar algorithms. Insertion sort or bubble sort, written as a fixed sequence of compare-exchanges, gives a network of size n(n-1)/2. Once you allow comparators to run in parallel its depth is 2n-3. Odd-even transposition sort alternates comparing pairs (0,1),(2,3),... and (1,2),(3,4),...; it needs n layers and suits hardware where only neighbouring cells can talk. Both are quadratic in size, so they are useful only for very small n or very restricted wiring.

The 0-1 principle and an exhaustive checker

How do you know a network sorts every input? Testing every permutation of n distinct values is n! cases, which is already 3.6 million at n = 10. The 0-1 principle reduces this to 2ⁿ: if a comparator network sorts every input made only of 0s and 1s, it sorts every input.

The proof is short. Let f be any monotone non-decreasing function. A comparator commutes with f, because min(f(a), f(b)) = f(min(a, b)) and the same holds for max. So a whole network commutes with f: applying f to the input and then running the network gives the same result as running the network and then applying f. Now suppose the network fails on some input, leaving an output where a larger value x comes before a smaller value y. Choose f(v) = 1 if v >= x, else 0. Applying f to that input gives a 0-1 input, and its output has a 1 before a 0, so the network also fails on a 0-1 input. By contraposition, sorting all 0-1 inputs is enough.

This turns into a checker that is complete, not just a random test. In pure Python it takes under a second at n = 16 and seconds at n = 20; versions that pack 64 binary inputs into one machine word go further:

from itertools import product

def apply(network, xs):
    xs = list(xs)
    for layer in network:
        for i, j in layer:
            if xs[i] > xs[j]:
                xs[i], xs[j] = xs[j], xs[i]
    return xs

def verify_01(network, n):
    """Return a failing 0-1 input, or None if the network sorts everything."""
    for bits in product((0, 1), repeat=n):
        out = apply(network, bits)
        if any(out[k] > out[k + 1] for k in range(n - 1)):
            return bits
    return None

N4 = [[(0, 1), (2, 3)], [(0, 2), (1, 3)], [(1, 2)]]
assert verify_01(N4, 4) is None
assert verify_01(N4[:2], 4) == (0, 1, 0, 1)   # drop layer 3 and it fails

The failing input (0, 1, 0, 1) shows exactly why the third comparator is needed: after two layers the minimum and maximum are in place, but the middle two values can still be out of order.

Batcher's networks: odd-even merge and bitonic

For n a power of two, Kenneth Batcher's 1968 constructions give networks of size O(n log² n) and depth O(log² n). Both are recursive merge sorts in which the merge step is itself a fixed network.

Odd-even merge sort sorts each half, then merges two sorted sequences by recursively merging their even-indexed and odd-indexed subsequences and finishing with one layer of comparators between neighbours. Bitonic sort sorts one half ascending and the other descending. Together they form a bitonic sequence (one that rises then falls), and a bitonic sequence can be sorted by comparing elements a distance n/2 apart, then n/4, and so on down to 1. Bitonic sort uses more comparators, but every layer has exactly n/2 comparators with a regular stride pattern, which is easier to map onto hardware.

def odd_even_merge_sort(n):
    """Batcher's odd-even merge sort, n a power of two. Returns (i, j) pairs."""
    pairs = []
    def merge(lo, hi, r):
        step = r * 2
        if step < hi - lo:
            merge(lo, hi, step)
            merge(lo + r, hi, step)
            for i in range(lo + r, hi - r, step):
                pairs.append((i, i + r))
        else:
            pairs.append((lo, lo + r))
    def sort(lo, hi):
        if hi - lo >= 1:
            mid = lo + (hi - lo) // 2
            sort(lo, mid)
            sort(mid + 1, hi)
            merge(lo, hi, 1)
    sort(0, n - 1)
    return pairs

def bitonic(n):
    """Bitonic sorter as layers; the pair is ordered (min_wire, max_wire)."""
    layers, k = [], 2
    while k <= n:
        j = k // 2
        while j >= 1:
            layer = []
            for i in range(n):
                l = i ^ j
                if l > i:
                    layer.append((i, l) if i & k == 0 else (l, i))
            layers.append(layer)
            j //= 2
        k *= 2
    return layers

def layerize(pairs, n):
    """Schedule a comparator list into the earliest layer each can run in."""
    ready, layers = [0] * n, []
    for i, j in pairs:
        d = max(ready[i], ready[j])
        if d == len(layers):
            layers.append([])
        layers[d].append((i, j))
        ready[i] = ready[j] = d + 1
    return layers

In the bitonic generator a pair (l, i) with l > i is a descending comparator: the minimum goes to the higher-numbered wire. apply handles both directions because it always puts the minimum on the first wire of the pair. Running both generators and checking them with verify_01 gives:

For n = 2^k both have depth k(k+1)/2. Bitonic has size (n/4)·k·(k+1), and odd-even merge sort has (k² − k + 4)·2^k−2 − 1. Hand-optimised networks for small fixed n save a few comparators, which matters when a network runs billions of times as a base case. In theory, the AKS network achieves O(log n) depth and O(n log n) size, but its constant factors are so large that it is never used in practice.

Worked example: tracing and vectorising a network

Trace the 4-input network from the diagram on 3, 1, 4, 2. Layer 1 compares (w0, w1) and (w2, w3): 3 and 1 swap, 4 and 2 swap, giving 1, 3, 2, 4. Layer 2 compares (w0, w2) and (w1, w3): 1 against 2 stays, 3 against 4 stays. Now w0 holds the minimum of all four and w3 the maximum, because each started layer 2 holding the minimum (or maximum) of its pair. Layer 3 compares (w1, w2): 3 and 2 swap, giving 1, 2, 3, 4. No step looked at a value to decide which comparison to do next. The same five operations run on every input.

Now scale this to SIMD. Load 8 rows of 4 floats into four 8-wide vector registers so that register r holds element r of every row. One comparator is then two vector instructions, min and max, that sort eight rows at once. The whole 4-input network is 10 instructions with no branches. The same idea in NumPy, applying a network column-wise to many rows, makes it easy to see:

import numpy as np

def sort_rows(X, layers):
    Y = X.copy()
    for layer in layers:
        i = np.array([a for a, _ in layer])
        j = np.array([b for _, b in layer])
        lo = np.minimum(Y[:, i], Y[:, j])
        hi = np.maximum(Y[:, i], Y[:, j])
        Y[:, i], Y[:, j] = lo, hi
    return Y

X = np.random.default_rng(0).integers(0, 100, size=(5, 8))
assert (sort_rows(X, layerize(odd_even_merge_sort(8), 8)) == np.sort(X, axis=1)).all()

In C or C++, the comparator is written so the compiler emits conditional moves or vector min/max instead of a branch:

static inline void cswap(int *a, int *b) {
    int x = *a, y = *b;
    *a = x < y ? x : y;     /* typically compiles to cmov, or min with SIMD */
    *b = x < y ? y : x;
}

Where networks win

Small sorts inside larger sorts. Vectorised sort libraries, for example Google's Highway vqsort and Intel's x86-simd-sort, use sorting networks to sort small blocks held in vector registers, and partitioning or merging for the rest. The partition and merge logic is covered in the quicksort and mergesort articles. At this scale a network that does 19 comparisons with no branches beats an insertion sort that does fewer comparisons but mispredicts its branches.

GPUs. A thread block can bitonic-sort a tile in shared memory, with one thread per comparator in each layer and a barrier between layers. Within a warp, register shuffles replace shared memory for the smallest strides. This shows up in per-row top-k selection for sampling and attention sparsity, where many short rows are sorted at once. Large GPU sorts usually use radix sort instead, because O(n log² n) loses to linear passes once n is large. The same blocking and shared-memory thinking appears in the GPU matrix multiply article.

Constant-time code. A sort whose memory accesses depend on secret data can leak through cache timing. Some post-quantum cryptography implementations sort secret arrays with networks built from branchless min/max (djbsort is a well-known example), so timing does not depend on the secret.

Hardware and fixed-size kernels. FPGA and ASIC sorters, median filters over fixed windows in image processing, and sorting the k candidates in a beam all have a fixed n known at design time, which is exactly where a network fits.

Failure modes

• n is not a power of two. Pad with +infinity sentinels and drop them afterwards, or prune the comparators that touch padding wires (a comparator against a +infinity wire never swaps). Check the pruned network with verify_01.
• NaN. Comparisons involving NaN are not a total order, and hardware min/max instructions treat NaN asymmetrically depending on operand order, so a NaN can duplicate or erase values. Map NaN to a sentinel before sorting.
• Stability. Networks are not stable. If equal keys must keep their order, pack the original index into the low bits of the key.
• Records, not keys. Moving a payload requires a key-index pair per comparator, which doubles the register pressure. Sort (key, index) pairs, then gather the payload once at the end.
• The compiler adds a branch. The ternary is a hint, not a guarantee. Check the assembly or use intrinsics in constant-time code.
• Testing with random inputs. Random tests miss rare failing patterns. Use the 0-1 check for every generated or hand-edited network up to about n = 20 (further with a bit-parallel checker), and depend on a proven construction above that.
• Using a network for big n. At a million elements, log² n is about 400 against log n of 20. Switch to merge or radix sort past a few dozen elements (or partial sort if you only need the top k), unless obliviousness is the requirement.

Trade-offs

What to do next

1. Run the generators above for n = 4, 8 and 16, and confirm the sizes and depths with verify_01.
2. Draw the 8-input odd-even merge sort and trace one input by hand.
3. Replace the insertion-sort base case in a quicksort with a branchless 8- or 16-element network and benchmark both on random and sorted data.
4. Use NumPy's column-wise approach to sort many short rows at once and compare it with np.sort(axis=1).
5. On a GPU, write a shared-memory bitonic sort of a 1,024-element tile, with one barrier per layer.
6. For constant-time code, inspect the generated assembly and check there are no data-dependent branches or memory indices.