Radix sort on sixteen cores, in .NET and Elixir

August 29, 2026 · personal · algorithms · dotnet · elixir · bake-off

the short version

Use the library sort. Introsort, pdqsort, timsort: whichever yours is, it takes any comparator, has had decades of tuning, and is right almost every time. Radix sort is the special tool for one situation, and when that situation shows up nothing else is close.

Reach for it when all of these hold:

• The keys are fixed-width and there are a lot of them. Integers, floats, timestamps, ids, hashes, fixed-length strings. At a few thousand keys everything is in cache and the library sort wins on constant factors; by a million the gap in this post is already sixfold.
• The keys sit in a contiguous array: a column rather than a list of objects. Sorting objects by a property pays for the indirection and the key extraction on every access, and that swamps the inner loop.
• The order can be written as bytes. Flip the sign bit of a signed int, flip every bit of a negative float, store big-endian, and (country, year, id) becomes one byte-comparable key. If the order needs a comparator, a locale or a lambda, radix sort isn't a candidate.
• It's on a hot path: a query operator, an index build, a per-frame pass, a pipeline stage. Something that sorts all day.

Leave it alone when the input is small or nearly sorted (timsort and pdqsort go linear on existing runs; radix does its passes regardless), when keys are variable-length strings (MSD radix handles them, but it's a different algorithm with a higher crossover), when memory is tight (radix needs a second buffer of n keys; introsort needs none), when keys are wide relative to n (128-bit keys are many passes), or when you'd be maintaining it in application code. The version that wins is a few hundred lines that need a benchmark harness to keep honest, and the places that use it own it.

Those places: DuckDB radix-sorts on binary-comparable keys, one thread per block, then merges. ClickHouse keeps an LSD and an MSD radix sort for numeric columns, the MSD one for LIMIT queries that only need the top few percent. NumPy's stable sort is a radix sort for integers of 16 bits or less. On a GPU it's the default: CUB's DeviceRadixSort is an LSD radix sort over every primitive type, and the Onesweep paper made it another 1.5× faster in 2022.

where it earns its keep: ORDER BY in a columnar database

The case I'd point at is an analytical database sorting a hundred million rows, because every condition above holds by construction. DuckDB's sort operator, as its authors described it in 2021, works like this.

First, every ORDER BY key is encoded into a fixed-width byte string per row. Integers are swapped to big-endian so byte order is numeric order, the sign bit is flipped so negatives sort before positives, a descending column has every bit inverted, NULL costs one extra byte up front, and a long string contributes a prefix, with the full string consulted only when two prefixes tie. The point of all that work is one property: memcmp order on the encoded keys is the query's order, and a byte-by-byte radix sort produces exactly memcmp order.

Second, each thread radix-sorts the rows it scanned, stably, which matters because ties in ORDER BY and the window functions that sit on top of it are defined by the previous order. Third, the sorted runs are merged in parallel, with Merge Path finding the cut points so that the merge itself splits across cores.

That's the pattern from this post in production shape: the keys are fixed-width because the encoder made them so, the sort is the private-then-meet-once shape, and the reason to pay for a custom sort is that it sits under every ORDER BY, merge join and top-N in the engine. Their numbers, from 2021: a hundred million integers in just under five seconds on one thread, and roughly three seconds in parallel, level with ClickHouse. That's far slower than the 74 ms below because an engine sorts rows, not bare ints: the payload columns travel with the keys, and the key encoding is a pass of its own. The algorithm is the same.

why I came back to it

In March 2015 I quoted a line from Wikipedia, "radix sorts are often, in practice, the fastest and most useful sorts on parallel machines". Eleven years and a sixteen-core laptop later, here is what it means with numbers attached. Ten million random integers: .NET's Array.Sort takes 525 ms, radix sort across sixteen threads takes 15, and with two tricks from a 2010 paper, 10. Elixir's Enum.sort takes 1.2 s, and radix sort across sixteen processes takes 101 ms. Same algorithm in both, and the interesting part is the one moment per pass where the workers have to meet.

why a comparison sort is hard to cut up

A comparison sort is a chain of decisions: which pair you compare next depends on how the last one went. Merging two sorted lists is one front-to-back scan, and you cannot place the 500th item before the 499th. Quicksort's partition is the same shape, and its recursion tree is only balanced if the pivots are lucky, so hand two halves to two cores and one of them may get most of the work. Even on one core, the branch at every comparison is unpredictable on random data, and a mispredicted branch stalls the pipeline.

Parallel comparison sorts exist (sample sort, bitonic networks), but they either do more than n log n work or need a clever divide step to keep every core busy.

the pass that parallelises itself

The 2015 post has the sequential mechanism: sort on the lowest digit with a stable counting sort, then the next digit up, until you run out. That's LSD, least significant digit first, as opposed to MSD, which partitions on the top digit and recurses into each run; LSD is the one that parallelises cleanly because every pass is the same fixed amount of work over the whole input. A counting sort is three steps, and each one falls apart into independent pieces. Cut the input into one slice per core:

1. Histogram. Each core counts the digits in its own slice into its own private array. Nothing shared, nothing to branch on.
2. Prefix sum. Add the counts up in a fixed order: digit 0 from core 0, digit 0 from core 1, ... then digit 1 from core 0, and so on. Each (core, digit) pair now knows where its run starts in the output. This is the only moment the cores meet, and it's over sixteen small arrays, not n keys.
3. Scatter. Each core writes its slice into the slots it was given in step 2. No locks, because no two keys want the same slot.

Stability survives the split because of the order in step 2: every 3 from slice 0 lands before every 3 from slice 1, and inside a slice keys go out in the order they came in. Sixteen cores produce exactly the array one core would have.

Nothing about the shape depends on the data. Every core gets n/16 keys on every pass, there's no pivot to get unlucky with, and a 32-bit key is done in three passes of 11 bits (2048 buckets, so the scan is over 16 × 2048 counts). Linear work, and parallel apart from one tiny meeting per pass. That is the claim.

.NET: threads and one array

Threads share memory, so the picture maps straight onto Parallel.For and one shared output array. The whole pass is the three steps above and nothing else:

System.Threading.Tasks.Parallel.For(0, workers, w =>
{
    var local = counts[w];
    Array.Clear(local);
    var hi = Math.Min((w + 1) * chunk, n);
    for (var i = w * chunk; i < hi; i++) local[(s[i] >> sh) & mask]++;
});

var running = 0;
for (var digit = 0; digit < radix; digit++)
    for (var w = 0; w < workers; w++)
    {
        var c = counts[w][digit];
        counts[w][digit] = running;
        running += c;
    }

System.Threading.Tasks.Parallel.For(0, workers, w =>
{
    var offsets = counts[w];
    var hi = Math.Min((w + 1) * chunk, n);
    for (var i = w * chunk; i < hi; i++)
    {
        var key = s[i];
        d[offsets[(key >> sh) & mask]++] = key;
    }
});

(src, dst) = (dst, src);

counts[w] is worker w's private histogram, and the scan overwrites it in place with the same worker's start offsets, so the scatter loop needs nothing it didn't already own. src and dst swap at the end of each pass.

Ten million random non-negative 32-bit ints, medians of seven runs, .NET 8 in a Linux arm64 container on an M3 Max with sixteen cores:

Two things stand out. Radix on one thread already beats the library sort four to one, from the branch-free inner loop and the fixed pass count. And the scaling flattens after eight threads: each pass streams the whole array through memory twice, and sixteen cores run out of memory bandwidth long before they run out of work. That's why 11 bits beats 8: three passes instead of four is a quarter less traffic, and traffic is the ceiling. Hoisting the bounds checks with spans and refs buys another 5% (15.4 ms), which tells you the loop was already memory-bound.

and then the literature

The version above is the textbook one. Wassenberg and Sanders got a CPU radix sort to within 12% of the machine's memory bandwidth in 2010 with three ideas, and all three port to C# without leaving managed code.

1. Reverse sorting. One MSD pass on the top ten bits partitions the input into 1,024 buckets using the shared scatter above. After that every bucket is sorted on its low 22 bits by the worker that owns it, privately, so the workers meet once instead of once per pass and the remaining passes run inside one core's cache.
2. Software write-combining. A scatter over 2,048 digits writes to 2,048 different places, which is the worst thing you can do to a cache. Instead, each digit's keys are staged in a 64-byte buffer and flushed a whole cache line at a time, so the random writes become sequential bursts. The paper flushes with non-temporal stores; .NET has no intrinsic for those on arm64, so this is the plain-store version.
3. One read for all histograms. A bucket's histograms for every remaining pass are counted in a single sweep.

The write-combining staging is one function:

static void Put(ref int buf, ref int fill, ref int offsets, ref int dst, int digit, int key)
{
    ref var f = ref Unsafe.Add(ref fill, digit);
    Unsafe.Add(ref buf, digit * Wc + f) = key;
    if (++f == Wc)
    {
        ref var off = ref Unsafe.Add(ref offsets, digit);
        Unsafe.CopyBlockUnaligned(ref Unsafe.As<int, byte>(ref Unsafe.Add(ref dst, off)),
                                  ref Unsafe.As<int, byte>(ref Unsafe.Add(ref buf, digit * Wc)), Wc * sizeof(int));
        off += Wc;
        f = 0;
    }
}

Wc is 16 ints, one cache line. CopyBlockUnaligned with a constant size is inlined by the JIT into four vector stores. And the private per-bucket work, once the MSD pass has run and each worker has been handed a contiguous run of buckets balanced by key count:

for (var b = firstBucket; b < lastBucket; b++)
{
    int lo = bucketStart[b], len = bucketStart[b + 1] - lo;
    if (len <= 1) continue;
    var bucket = dst.AsSpan(lo, len);
    var tmp = scratch.AsSpan(0, len);

    // all this bucket's histograms from one read
    for (var p = 0; p < passes; p++) Array.Clear(hists[p]);
    foreach (var key in bucket)
        for (var p = 0; p < passes; p++) hists[p][(key >> (p * localBits)) & mask]++;

    Span<int> src = bucket, dstSpan = tmp;
    for (var p = 0; p < passes; p++)
    {
        LocalPass(src, dstSpan, hists[p], p * localBits, mask, radix, buf, fill);
        var swap = src; src = dstSpan; dstSpan = swap;
    }
    if (passes % 2 == 1) tmp.CopyTo(bucket);   // odd pass count ends in scratch
}

The scratch buffer, histograms and write-combine buffers are allocated once per worker and reused across its buckets. Same container, medians of seven runs at 10M and three at 100M:

Write-combining is worth 1.5× on a single thread, which is the paper's number, and the bigger the array the more it matters on sixteen: 1.2× at 10M, 1.6× at 100M, because the cost it removes is the scatter's random writes missing cache, and a 400 MB destination misses more. A hundred million ints in 74 ms is 82× Array.Sort.

Elixir: processes and :atomics

There's no shared array to start from on the BEAM: process heaps are private and a message is a copy into the receiver's heap. Lists and maps in messages would put a copy of every key on every pass, and the copying would cost more than the bucketing. What the runtime does have is :atomics (OTP 21.2): an off-heap array of 64-bit integers, referenced like a binary and shared between processes by reference, with atomic get, put, add and add_get. It's the C array this algorithm was written for. Two n-slot arrays as source and destination, swapped each pass; each worker's histogram and offsets in a small :atomics of its own, so the inner loops allocate nothing. The pass:

hists =
  ranges
  |> Enum.map(fn {lo, hi} -> Task.async(fn -> h = :atomics.new(radix, []); histogram(src, lo, hi, shift, mask, h); h end) end)
  |> Enum.map(&Task.await(&1, :infinity))

Enum.reduce(0..(radix - 1), 1, fn d, running ->
  Enum.reduce(hists, running, fn h, running ->
    c = :atomics.get(h, d + 1)
    :atomics.put(h, d + 1, running)
    running + c
  end)
end)

ranges
|> Enum.zip(hists)
|> Enum.map(fn {{lo, hi}, off} -> Task.async(fn -> scatter(src, dst, lo, hi, shift, mask, off) end) end)
|> Enum.each(&Task.await(&1, :infinity))

Same three steps as the C#, same in-place scan turning each worker's histogram into its offsets. The scatter loop is five BIF calls per key and nothing else:

defp scatter(_src, _dst, lo, hi, _shift, _mask, _off) when lo > hi, do: :ok
defp scatter(src, dst, lo, hi, shift, mask, off) do
  x = :atomics.get(src, lo)
  pos = :atomics.add_get(off, ((x >>> shift) &&& mask) + 1, 1) - 1
  :atomics.put(dst, pos, x)
  scatter(src, dst, lo + 1, hi, shift, mask, off)
end

add_get bumps the offset for that digit and returns the new value, so pos is this key's slot and the next key with the same digit gets the one after. The offsets array is the worker's own and the slots it writes are disjoint from every other worker's by construction, so there is nothing to lock and no compare-and-swap anywhere.

Ten million keys, medians of five, Elixir 1.19 on OTP 27, same machine:

Twelve times faster than Enum.sort, and the same shape as the .NET curve: 11 bits beats 8, and sixteen workers give 8.3× over one. Loading the keys from a binary into the array is 13 ms with sixteen workers each filling their own range, and reading the sorted array back out to a binary is 12, so a process that receives ten million ints as bytes can hand back sorted bytes in about 130 ms.

The one-core gap to .NET (843 ms against 85) is the price of the runtime: every :atomics call is a BIF call at tens of nanoseconds, where the C# loop does a memory access. Parallelism buys most of that back. The only way to buy the rest is a NIF, and at that point you are writing the C# version in Rust.

Nx, for the record

If the keys are already a tensor, there's a one-liner. A tensor is a contiguous binary, and Nx.sort on the EXLA backend hands that buffer to compiled XLA code:

t = Nx.from_binary(bin, :s32, backend: EXLA.Backend)
sorted = Nx.sort(t)

192 ms for ten million keys and 17 for a million, without the BEAM touching an element. It's a comparison sort on one core (I sampled the OS process while it looped: 104% CPU), so it's six times faster than Enum.sort and still behind the sixteen-process :atomics version, and getting a list into a tensor and back (336 and 382 ms at 10M) costs more than the sort. On a CUDA GPU the same call dispatches to a radix sort, which is the parallel machine the Wikipedia line was written about.

the two runtimes, side by side

Ten million random 32-bit ints on the same sixteen cores:

what to take from it

The parallelism is in the algorithm's shape. A radix pass is private histograms, one tiny meeting, and a scatter into slots nobody else wants, and that shape survives any runtime that can give sixteen workers one mutable array. Past eight cores the ceiling is memory bandwidth, so fewer passes beats more threads, and the literature's tricks are all ways of touching memory less: meet once, then stay in cache; stage writes into cache lines.

On the BEAM the array is the design decision. Anything that crosses a process boundary by value pays a copy per key per pass; :atomics is the shared mutable array that doesn't, and it carries this algorithm to within a runtime constant of the C#.

where to look

• The 2015 post this follows on from: MIT: Introduction to Algorithms, with the sequential mechanism and the decimal-digit versions in C# and Elixir.
• The next post: Sorting a billion integers from a stream, the same sort at the scale where the network is the bottleneck.
• The lecture: 6.006 Counting Sort, Radix Sort, Lower Bounds for Sorting on MIT OpenCourseWare.
• Zagha and Blelloch, "Radix sort for vector multiprocessors", Supercomputing '91: the paper behind the Wikipedia line.
• Laurens Kuiper, Fastest table sort in the West: redesigning DuckDB's sort, 2021: binary-comparable keys, per-thread radix sort, Merge Path.
• ClickHouse's RadixSort.h: LSD stable and MSD partial, with the sign and float bit transforms in the header comment.
• NumPy's 1.17 release notes: radix sort for integer types of 16 bits or less under kind="stable".
• CUB's DeviceRadixSort and Adinets and Merrill, Onesweep: a faster least significant digit radix sort for GPUs, 2022.
• Jan Wassenberg and Peter Sanders, Faster Radix Sort via Virtual Memory and Write-Combining, 2010: reverse sorting and software write-combining, with the bandwidth measurements.
• Satish et al., "Fast sort on CPUs and GPUs: a case for bandwidth oblivious SIMD sort", SIGMOD 2010: the Intel radix sort the paper above beats by 1.5×.
• Marc Gravell, Sorting myself out, extreme edition, 2018: radix sort in C# on Span<T>.
• HPCsharp: a C# library with parallel radix sorts, if you'd rather not maintain your own.
• Erlang's atomics module.
• Nx and the EXLA backend.