Parallel O(sqrt n) Overhead LSD Radix Sort

# Parallel 𝒪(√𝑛) Overhead LSD Radix Sort
Source: [https://arxiv.org/html/2607.05302](https://arxiv.org/html/2607.05302)
Zuse Institute Berlin, Germanyclausecker@zib\.deNHR\-Verein Graduate School programme Zuse Institute Berlin, Germanyschintke@zib\.dehttps://orcid\.org/0000\-0003\-4548\-788X\{CCSXML\}¡ccs2012¿ ¡concept¿ ¡concept\_id¿10003752\.10003809\.10010031\.10010033¡/concept\_id¿ ¡concept\_desc¿Theory of computation Sorting and searching¡/concept\_desc¿ ¡concept\_significance¿500¡/concept\_significance¿ ¡/concept¿ ¡/ccs2012¿

This work was produced during a secondment of the first author at the Sanders workgroup at KIT\. The author would like to thank Marvin Williams for his invaluable input\.

## Parallel𝒪⁡\(n\)\\mathcal\{O\}\(\\sqrt\{n\}\)Overhead LSD Radix Sort

###### Abstract

We present Radsort, a variant of LSD radix sort, sorting data with𝒪⁡\(n\)\\mathcal\{O\}\(\\sqrt\{n\}\)additional space\. Radsort is stable, admits a simple implementation and is easy to parallelise\. For arrays exceeding a size of around2​MiB2\\,\\mathrm\{MiB\}it outperforms a conventional out\-of\-place LSD radix sort\.

###### ccs

Theory of computation Sorting and searching

###### keywords

cache locality, radix sort, sorting

††runningauthor:R\. Clausecker and F\. Schintke††copyright:Robert Clausecker and Florian Schintke††supplement:https://github\.com/clausecker/radsort††editors:John Q\. Open and Joan R\. Access††event\-title:42nd Conference on Very Important Topics \(CVIT 2016\)††event\-shorttitle:CVIT 2016††event\-acronym:CVIT††year:2016††event\-date:December 24–27, 2016††event\-location:Little Whinging, United Kingdom††series\-volume:42††articleno:23## 1Introduction

An LSD \(least significant digit\) radix sort algorithm sorts an arrayAAof lengthnnholdingntn\_\{t\}\-tuples of keys\(key​\(A⁡\[i\],0\),key​\(A⁡\[i\],1\),…,key​\(A⁡\[i\],nt−1\)\)\\bigl\(\\textit\{key\}\(A\[i\],0\),\\textit\{key\}\(A\[i\],1\),\\ldots,\\textit\{key\}\(A\[i\],n\_\{t\}\-1\)\\bigr\)of alphabetΣ\\Sigmawithσ\\sigmasymbols into lexicographic order by repeatedly sorting according tokey​\(A​\[i\],nt−1\)\\textit\{key\}\(A\[i\],n\_\{t\}\-1\), thenkey​\(A​\[i\],nt−2\)\\textit\{key\}\(A\[i\],n\_\{t\}\-2\), and so on untilkey​\(A​\[i\],0\)\\textit\{key\}\(A\[i\],0\)\. With a stable sort, identical keys end up sorted by their suffixes, giving a lexicographic order\. A conventional LSD radix sort sorts bykey​\(A​\[i\],t\)\\textit\{key\}\(A\[i\],t\)by first taking a histogramHHof how often eachc∈Σc\\in\\Sigmaoccurs\. A prefix sum overHHobtains the bucket startsSS, with which we sortAAinto a new arrayA′A^\{\\prime\}\(Alg\.[1](https://arxiv.org/html/2607.05302#alg1)\)\.

Algorithm 1out\-of\-place LSD radix sort round1procedureradixSortStep\(

A′,A,tA^\{\\prime\},A,t\)⊳\\trianglerightsort arrayAAby keyttintoA′A^\{\\prime\}

2

H←0​…​0H\\leftarrow 0\\ldots 0
3for

i←0​…​n−1i\\leftarrow 0\\ldots n\-1do⊳\\trianglerighttake a histogram ofkey​\(A​\[i\],t\)\\textit\{key\}\(A\[i\],t\)

4

H⁡\[key​\(A⁡\[i\],t\)\]←H⁡\[key​\(A⁡\[i\],t\)\]\+1H\[\\textit\{key\}\(A\[i\],t\)\]\\leftarrow H\[\\textit\{key\}\(A\[i\],t\)\]\+1
5endfor

6

S⁡\[0\]←0S\[0\]\\leftarrow 0
7for

i←1​…​σ−1i\\leftarrow 1\\ldots\\sigma\-1do

8

S⁡\[i\]←S⁡\[i−1\]\+H⁡\[i−1\]S\[i\]\\leftarrow S\[i\-1\]\+H\[i\-1\]⊳\\trianglerightprefix sum overHH

9endfor

10for

i←0​…​n−1i\\leftarrow 0\\ldots n\-1do

11

c←key​\(A​\[i\],t\)c\\leftarrow\\textit\{key\}\(A\[i\],t\)⊳\\trianglerightdetermine bucket

12

A′​\[S⁡\[c\]\]←A⁡\[i\]A^\{\\prime\}\\bigl\[S\[c\]\\bigr\]\\leftarrow A\[i\]⊳\\trianglerightsortA⁡\[i\]A\[i\]intoA′A^\{\\prime\}

13

S⁡\[c\]←S⁡\[c\]\+1S\[c\]\\leftarrow S\[c\]\+1⊳\\trianglerightadvance bucketccpast new element

14endfor

15endprocedure

With the input arrayAAbeing consumed, output is produced in a separate arrayA′A^\{\\prime\}\. While onlynnelements are occupied in both arrays together at any point in the algorithm, nevertheless2​n2nelements of storage, i\. e\.nnelements of overhead in addition to the input, need to be provided\. We would like to address this shortcoming with a novel LSD radix sort algorithm that does the job with only𝒪⁡\(n\)\\mathcal\{O\}\(\\sqrt\{n\}\)extra space at a similar performance\. We call this algorithm*Radsort*due to it being aradical overheadradix sort\.

Algorithm 2the Radsort algorithm1procedureradSort\(

A,n,ntA,n,n\_\{t\}\)

2

setup​\(A,n\)\\textsc\{setup\}\(A,n\)⊳\\trianglerightAlg\.[4](https://arxiv.org/html/2607.05302#alg4)

3for

t←nt−1​…​0t\\leftarrow n\_\{t\}\-1\\ldots 0do

4

sortPhase​\(t\)\\textsc\{sortPhase\}\(t\)⊳\\trianglerightAlg\.[5](https://arxiv.org/html/2607.05302#alg5)

5fixupPhase⊳\\trianglerightAlg\.[6](https://arxiv.org/html/2607.05302#alg6)

6endfor

7Finalize⊳\\trianglerightAlg\.[7](https://arxiv.org/html/2607.05302#alg7)

8endprocedure

## 2Radsort

Radsort \(Alg\.[2](https://arxiv.org/html/2607.05302#alg2)\) treats the input arrayAAas a sequence of*blocks*of some block sizebb111see §[3](https://arxiv.org/html/2607.05302#S3)for discussion on the selection ofbb\. After a block of input is consumed, its storage is reused for subsequently produced output, reducing the storage overhead to a constant number of*scratch blocks*, as well as some bookkeeping\.

Each round of sorting comprises two*phases*\. During the*sort phase*, one output block is initially assigned to each of theσ\\sigma*buckets*we sort into\. Once such a block is full, a new block is drawn from previously consumed input blocks, overwriting the input as it is processed\. The result of the sort phase areσ\\sigmarandomly interleaved sequences of output blocks, each forming one of the output buckets\. In the*fixup phase*following each*sort phase*, these sequences are deinterleaved, giving the data block in order of buckets\. The deinterleaving does not move data blocks around—instead, a permutationπ\\pitracks the order in which the data blocks need to be processed for the data to appear in order\.

Like in a conventional LSD radix sort, these two phases are repeated for each key position in turn\. After allntn\_\{t\}rounds of sorting, a*finalisation step*shuffles the data blocks around, giving a sorted result inAA\. While Radsort allocates space in a more complicated way than a traditional LSD radix sort, it works according to the same principles, giving the same stability guarantees at lower memory overhead and better cache locality\.

### 2\.1Data Structures

In addition to the⌊n/b⌋\\lfloor n/b\\rfloorblocks ofAA,2​σ2\\sigma*scratch blocks*stored in the temporary arrayTTare required:σ\\sigmablocks of*head start*to cover the initial output blocks for the buckets of the current round, andσ\\sigmaadditional covering the partial blocks of the previous round of sorting\. The algorithm treats the concatenation ofTTandAAas one long array ofnπn\_\{\\pi\}blocks ofbbelements each, indexed through theblockatfunction\. The first2​σ2\\sigmaindices refer to blocks inTT, the other indices to blocks overlayingAA\. Theblockat​\(i\)\\textit\{blockat\}\(i\)function returns a pointer to the beginning of the block at indexii\.

blockat​\(i\)=\{address ofT⁡\[i​b\]if0≤i<2​σaddress ofA⁡\[\(i−2​σ\)​b\]if2​σ≤i<nπ\\textit\{blockat\}\(i\)=\\begin\{cases\}\\mbox\{address of $T\[ib\]$\}&\\mbox\{if $0\\leq i<2\\sigma$\}\\\\ \\mbox\{address of $A\[\(i\-2\\sigma\)b\]$\}&\\mbox\{if $2\\sigma\\leq i<n\_\{\\pi\}$\}\\end\{cases\}\(1\)Thenmodbn\\bmod belements at the end ofAAare copied into a scratch block at the start of the algorithm\. They are unused until the contents ofTTare copied back toAAat the end of the algorithm\. The arrayπ\\pidescribes a permutation of thesenπn\_\{\\pi\}blocks, such thatπ⁡\[i\]\\pi\[i\]holds the index of the block at*logical index*ii\. Fig\.[1](https://arxiv.org/html/2607.05302#S2.F1)shows an example of this data structure for strings ofnt=4n\_\{t\}=4characters on initialisation\.

Each block is either*allocated*or*unallocated*\. An allocated block is either*full*, meaning that all of its elements hold data, or*partial*, meaning that some of its elements hold data\. Unallocated blocks do not hold any data and are ready for reuse\. At the firstσ\\sigmalogical indices, a*head start*ofσ\\sigmaunallocated blocks is found\. They are followed by allocated blocks at logical indicesσ≤i<f\\sigma\\leq i<f, and finally unallocated blocks are logical indicesf≤i<nπf\\leq i<n\_\{\\pi\}\.

ArrayPPtracks up toσ\\sigmapartial blocks\. It holds in order of logical index partial blocks\(i,l\)\(i,l\)described by logical indexiiand block length0≤l<b0\\leq l<b\. Of each such partial block, the firstllelements are used, and the otherb−lb\-lelements are unused\. To simplify the implementation, the block at logical indexf−1f\-1is always a partial block\.

Algorithm 3finding the size of the next block1functionsizeOfNextBlock\(

i,iPi,i\_\{P\}\)

2

\(iiP,liP\)←P⁡\[iP\]\(i\_\{i\_\{P\}\},l\_\{i\_\{P\}\}\)\\leftarrow P\[i\_\{P\}\]
3if

π⁡\[i\]=iiP\\pi\[i\]=i\_\{i\_\{P\}\}then⊳\\trianglerightisπ⁡\[i\]\\pi\[i\]a partial block?

4

iP←iP\+1i\_\{P\}\\leftarrow i\_\{P\}\+1
5return

liPl\_\{i\_\{P\}\}
6elsereturn

bb
7endif

8endfunction

Concatenated in the order given byπ\\pi, the contents of full and partial blocks hold the elements to be sorted in the order induced by the rounds of sorting performed so far\. We can iterate over the allocated blocks in logical order, discovering the length of each block by traversingPPat the same time\. This idea is wrapped in the functionsizeOfNextBlock​\(i,iP\)\\textsc\{sizeOfNextBlock\}\(i,i\_\{P\}\)which gives the length of the block at logical indexii, assuming the next partial block is tracked inP⁡\[iP\]P\[i\_\{P\}\]\. If the block is partial,iPi\_\{P\}is updated to track the next partial block\.

In summary, the state of the Radsort algorithm is tracked in the variables given below\. Additional variables are needed within individual algorithm steps and are explained there\.

AAinput arrayarray of stringsnninput lengthn=\|A\|n=\|A\|σ\\sigmaalphabet sizeσ=\|Σ\|\\sigma=\|\\Sigma\|ntn\_\{t\}key lengthintegerbbblock sizeintegernπn\_\{\\pi\}block countnπ=⌊n/b⌋\+2​σn\_\{\\pi\}=\\lfloor n/b\\rfloor\+2\\sigmaTTscratch bufferarray of2​σ​b2\\sigma bstringsfffill levelintegerπ\\pipermutationarray ofnπn\_\{\\pi\}integersPPpartial blocksarray ofσ\\sigmapairs of integers
TT2​σ2\\sigma⋯\\cdotsAA⌊n/b⌋\\lfloor n/b\\rfloor⋯\\cdotsnmodbn\\bmod btail elementsπ\\piσ\\vphantom\{\\leq\}\\sigmahead start≤σ\\leq\\sigmaunallocateddata entries01σ\\sigmaffnπ−1n\_\{\\pi\}\-1logical indexa block withbbelementsofnt=4n\_\{t\}=4characters eachfrob quux⋯\\cdotsbuzz⋯\\cdots⋯\\cdots⋯\\cdotsFigure 1:Radsort’s data structures at initialisationacbaccabacbacabbcabbcabcbc0011223344556677889910101111TTAA*physical**logical*bacbccacabbacbaabbcabcbcca6611111133227700559910108844π⁡\[0\]\\pi\[0\]=iout=0,i\_\{\\mathrm\{out\}\}=0,iin=3i\_\{\\mathrm\{in\}\}=3P⁡\[0\]=\(2,3\)P\[0\]=\(2,3\)P⁡\[1\]=\(5,1\)P\[1\]=\(5,1\)P⁡\[2\]=\(8,2\)P\[2\]=\(8,2\)f=11f=11before sort phasebbacacba*exemplary elements:*ccccbcaaaaabbbbacabbcabcbc0011223344556677889910101111TTAAaaaabbbbcccccbaabbcabcbcca6611111133227700559910108844iouti\_\{\\mathrm\{out\}\}iin=7i\_\{\\mathrm\{in\}\}=7B⁡\[c\]B\[c\]=\[ , \]B⁡\[b\]B\[b\]=\[ , \]B⁡\[a\]B\[a\]=\[ , \]f=11f=11middle ofsort phasebbbbccccaaaccccbcaaaaabbbb0011223344556677889910101111TTAAcccbbbaaabcaaaabbbbcccccba6611111133227700559910108844iouti\_\{\\mathrm\{out\}\}iin=fi\_\{\\mathrm\{in\}\}=\{\\color\[rgb\]\{0,0,0\}f\}B⁡\[c\]B\[c\]=\[ , \]B⁡\[b\]B\[b\]=\[ , \]B⁡\[a\]B\[a\]=\[ , \]aftersort phasebccccbbbbcccccaaaaaaaabbbb0011223344556677889910101111TTAAaaaaaaaabbbbbbbbbccccccccc1010884466779911112200113355π⁡\[0\]\\pi\[0\]=P⁡\[0\]=\(9,0\)P\[0\]=\(9,0\)P⁡\[1\]=\(0,1\)P\[1\]=\(0,1\)P⁡\[2\]=\(5,1\)P\[2\]=\(5,1\)f=12f=12*block withb=4b=\\mathit\{4\}elements**tail*afterfixup phaseFigure 2:Algorithm state during the second round of sorting an array of strings withn=26n=26,nt=4n\_\{t\}=4,b=4b=4, andΣ=\{a,b,c\}\\Sigma=\\\{\\texttt\{a\},\\texttt\{b\},\\texttt\{c\}\\\}\. Changes highlighted in red\. Unallocated cells are blanked out\.
### 2\.2Initialisation

The variables are initially set up \(cf\. Fig\.[1](https://arxiv.org/html/2607.05302#S2.F1)\) such that the contents ofAAare represented by the blocks at logical indicesσ\\sigmatoσ\+⌈n/b⌉\\sigma\+\\lceil n/b\\rceil\. The finalnmodbn\\bmod belements ofAAare copied to blockσ\\sigma\. This establishes the*data invariant*of the state variables, shown below\. Each round of sorting assumes that this invariant holds initially and reestablishes it at conclusion\.

1. 1\.π\\piis a permutation of0​…​nπ−10\\ldots n\_\{\\pi\}\-1
2. 2\.PPtracks the partial blocks in ascending order of logical index
3. 3\.the block at logical indexf−1f\-1is a partial block
4. 4\.logical indices00toσ−1\\sigma\-1andfftonπ−1~n\_\{\\pi\}\-1refer to unallocated blocks
5. 5\.logical indicesσ\\sigmatof−1f\-1refer to allocated blocks, the concatenation of whose used elements holds the input permuted according to the current progress of the algorithm

Algorithm 4initialisation of variables1proceduresetup\(

A,nA,n\)

2

f←σ\+⌊n/b⌋\+1f\\leftarrow\\sigma\+\\lfloor n/b\\rfloor\+1
3

π⁡\[0​…​σ−1\]←0,…,σ−1\\pi\[0\\ldots\\sigma\-1\]\\leftarrow 0,\\ldots,\\sigma\-1⊳\\trianglerightassign head start

4

π⁡\[σ​…​f−2\]←2​σ,…,σ\+f−2\\pi\[\\sigma\\ldots f\-2\]\\leftarrow 2\\sigma,\\ldots,\\sigma\+f\-2⊳\\trianglerightassign contents ofAA

5

π⁡\[f−1​…​nπ−1\]←σ,…,2​σ−1\\pi\[f\-1\\ldots n\_\{\\pi\}\-1\]\\leftarrow\\sigma,\\ldots,2\\sigma\-1⊳\\trianglerightassign remaining blocks

6

T⁡\[0​…​2​σ​b−1\]←0T\[0\\dots 2\\sigma b\-1\]\\leftarrow 0⊳\\trianglerightfillyTTwith dummy values \(optional\)

7

blockat​\(σ\)​\[0​…​\(nmodb\)−1\]←A⁡\[b⁡⌊n/b⌋​…​n−1\]\\textit\{blockat\}\(\\sigma\)\[0\\dots\(n\\bmod b\)\-1\]\\leftarrow A\[b\\lfloor n/b\\rfloor\\ldots n\-1\]⊳\\trianglerightcopy tail to scratch block

8

P⁡\[0\]←\(f−1,nmodb\)P\[0\]\\leftarrow\(f\-1,n\\bmod b\)⊳\\trianglerighttrack input tail as partial block

9

P⁡\[1​…​σ−1\]←\(nπ,0\)P\[1\\ldots\\sigma\-1\]\\leftarrow\(n\_\{\\pi\},0\)⊳\\trianglerightfillPPwith dummy values \(optional\)

10endprocedure

### 2\.3Sorting

Each round of sorting \(cf\. Fig\.[2](https://arxiv.org/html/2607.05302#S2.F2)\) sorts the elements ofAAstably according to some key positiontt\. We writekey​\(A​\[i\],t\)∈Σ\\textit\{key\}\(A\[i\],t\)\\in\\Sigmafor key positionttof elementA⁡\[i\]A\[i\]\. First, the*sort phase*sorts elements into buckets according to the value of their key\. Then, the*fixup phase*computes newπ\\piandPPbased on the results of the sort phase, restoring the data invariant\. The sort and fixup phases require several additional variables that can be discarded afterwards:

BBbucket allocationarray ofσ\\sigmapairs of pointersCCblock countsarray ofσ\\sigmaintegersSSbucket startsarray ofσ\\sigmaintegersUUblock usagearray ofnπn\_\{\\pi\}charactersπ′\\pi^\{\\prime\}new permutationarray ofnπn\_\{\\pi\}integersiini\_\{\\mathrm\{in\}\}next input blockintegeriouti\_\{\\mathrm\{out\}\}next output blockintegeriPi\_\{P\}next partial blockinteger
#### 2\.3\.1Sort phase

The sort phase allocates one block for each bucket in the bucket arrayBB\. This array holds for each bucketcca pair of pointers\(pnext,pend\)\(p\_\{\\mathrm\{next\}\},p\_\{\\mathrm\{end\}\}\)pointing to the next free element, and the end of the block\. It then traverses the allocated blocks byiini\_\{\\mathrm\{in\}\}in logical order, and for each block sorts its elements into the correct output block\. If an output block is full, a new block is drawn from the next logical indexiouti\_\{\\mathrm\{out\}\}using thenewBucketBlock​\(c\)\\textsc\{newBucketBlock\}\(c\)procedure\. The bucket each output block is used for is tracked in the usage arrayUU, and the number of blocks assigned to each of the buckets is tracked inCC\.

As the sort proceeds through the input array, consumed blocks of input are reassigned into output blocks\. It is critical that blocks of input are fully consumed before they are reassigned into output blocks, lest elements are overwritten before they can be sorted\. This is ensured throughiini\_\{\\mathrm\{in\}\}having a*head start*ofσ\\sigmaunallocated blocks:

###### Lemma 1\.

When a new output block needs to be assigned,iout<iini\_\{\\mathrm\{out\}\}<i\_\{\\mathrm\{in\}\}\.

###### Proof\.

Letjjbe the number of elements processed so far\. The input comprises of blocks holding at mostbbelements, so at least⌊j/b⌋\\lfloor j/b\\rfloorblocks of input have been processed so far, henceσ\+⌊j/b⌋≤iin\\sigma\+\\lfloor j/b\\rfloor\\leq i\_\{\\mathrm\{in\}\}\. As for output blocks, a new one is assigned whenever the current block of some bucket is full\. Thus there are nowσ−1\\sigma\-1partial output blocks \(one for each bucket, except for the bucket we need to assign a new output block to\) as well as up to⌊j/b⌋\\lfloor j/b\\rfloorfull output blocks\. Henceiout≤σ−1\+⌊j/b⌋i\_\{\\mathrm\{out\}\}\\leq\\sigma\-1\+\\lfloor j/b\\rfloor\. Joining these two together gives

iout≤σ−1\+⌊j/b⌋<σ\+⌊j/b⌋≤iini\_\{\\mathrm\{out\}\}\\leq\\sigma\-1\+\\lfloor j/b\\rfloor<\\sigma\+\\lfloor j/b\\rfloor\\leq i\_\{\\mathrm\{in\}\}\(2\)as required\. ∎

Algorithm 5sort phase1procedurenewBucketBlock\(

cc\)⊳\\trianglerightallocate a new bucket block for keycc

2

pnext←blockat​\(π⁡\[iout\]\)p\_\{\\mathrm\{next\}\}\\leftarrow\\textit\{blockat\}\(\\pi\[i\_\{\\mathrm\{out\}\}\]\)
3

B⁡\[c\]←\(pnext,pnext\+b\)B\[c\]\\leftarrow\(p\_\{\\mathrm\{next\}\},p\_\{\\mathrm\{next\}\}\+b\)
4

C⁡\[c\]←C⁡\[c\]\+1C\[c\]\\leftarrow C\[c\]\+1
5

U⁡\[iout\]←cU\[i\_\{\\mathrm\{out\}\}\]\\leftarrow c
6

iout←iout\+1i\_\{\\mathrm\{out\}\}\\leftarrow i\_\{\\mathrm\{out\}\}\+1
7endprocedure

8proceduresortPhase\(t\)

9

C⁡\[0​…​σ−1\]←iP←iout←0C\[0\\ldots\\sigma\-1\]\\leftarrow i\_\{P\}\\leftarrow i\_\{\\mathrm\{out\}\}\\leftarrow 0⊳\\trianglerightinitialise variables

10forc∈Σc\\in\\SigmadonewBucketBlock\(c\)⊳\\trianglerightassign initial output blocks

11for

iin←σ​…​f−1i\_\{\\mathrm\{in\}\}\\leftarrow\\sigma\\ldots f\-1do

12

pin←blockat​\(π⁡\[iin\]\)p\_\{\\mathrm\{in\}\}\\leftarrow\\textit\{blockat\}\(\\pi\[i\_\{\\mathrm\{in\}\}\]\)⊳\\trianglerightget next block of input

13

l←sizeOfNextBlock​\(iin,iP\)l\\leftarrow\\textsc\{sizeOfNextBlock\}\(i\_\{\\mathrm\{in\}\},i\_\{P\}\)
14for

j←0​…​l−1j\\leftarrow 0\\ldots l\-1do

15

c←key​\(pin​\[j\],t\)c\\leftarrow\\textit\{key\}\(p\_\{\\mathrm\{in\}\}\[j\],t\)⊳\\trianglerightfind bucket of current element

16

\(pnext,pend\)←B⁡\[c\]\(p\_\{\\mathrm\{next\}\},p\_\{\\mathrm\{end\}\}\)\\leftarrow B\[c\]
17copy

pin​\[j\]p\_\{\\mathrm\{in\}\}\[j\]to memory at

pnextp\_\{\\mathrm\{next\}\}
18if

pnext\+1=pendp\_\{\\mathrm\{next\}\}\+1=p\_\{\\mathrm\{end\}\}then⊳\\trianglerightbucket full?

19

newBucketBlock​\(c\)\\textsc\{newBucketBlock\}\(c\)
20else

21

B⁡\[c\]←\(pnext\+1,pend\)B\[c\]\\leftarrow\(p\_\{\\mathrm\{next\}\}\+1,p\_\{\\mathrm\{end\}\}\)
22endif

23endfor

24endfor

25endprocedure

#### 2\.3\.2Fixup phase

After the sort phase, the blocks at logical index00toiouti\_\{\\mathrm\{out\}\}hold the sorted buckets in some interleaving, with up tobbpartial blocks described byBB\. The fixup phase computes a new permutationπ′\\pi^\{\\prime\}that starts with a head start ofσ\\sigmaunallocated blocks, followed by the deinterleaved buckets in order of keys and finally any remaining unallocated blocks\. We are guaranteed to always findσ\\sigmaunallocated blocks for the head start:

###### Lemma 2\.

After the sort phase, there are at leastσ\\sigmaunallocated blocks\.

###### Proof\.

A full block holdsbbelements\. As there arennelements in total, there are at most⌊n/b⌋\\lfloor n/b\\rfloorfull output blocks\. Furthermore, there areσ\\sigmapartial blocks tracked inBB\. As each block after sorting is either used as an output block or unallocated, and as there arenπn\_\{\\pi\}blocks in total, we find for the number of unallocated blocksnunallocatedn\_\{\\mathrm\{unallocated\}\}

nunallocated≥nπ−σ−⌊n/b⌋=σn\_\{\\mathrm\{unallocated\}\}\\geq n\_\{\\pi\}\-\\sigma\-\\lfloor n/b\\rfloor=\\sigma\(3\)as required\. ∎

The buckets are deinterleaved by first computing the starting indicesS⁡\[0​…​σ−1\]S\[0\\ldots\\sigma\-1\]of the buckets and then by shuffling the blocks according to the usagesUUintoπ′\\pi^\{\\prime\}\. Meanwhile, each allocated block is compared with the pointers inBBto populatePPwith the whereabouts of the new partial blocks\. As there is one partial block for each bucket, each of theσ\\sigmaelements ofPPends up being initialised\. Finally,π\\piis set toπ′\\pi^\{\\prime\}to restore the data invariant\.

Algorithm 6fixup phase1procedurefixupPhase

2

S⁡\[0\]←σS\[0\]\\leftarrow\\sigma
3fori←1​…​σ−1i\\leftarrow 1\\ldots\\sigma\-1do

S⁡\[i\]←S⁡\[i−1\]\+C⁡\[i−1\]S\[i\]\\leftarrow S\[i\-1\]\+C\[i\-1\]⊳\\trianglerightprefix sum overCCplusσ\\sigma

4

π′​\[0​…​σ−1\]←π⁡\[iout​…​iout\+σ−1\]\\pi^\{\\prime\}\[0\\ldots\\sigma\-1\]\\leftarrow\\pi\[i\_\{\\mathrm\{out\}\}\\ldots i\_\{\\mathrm\{out\}\}\+\\sigma\-1\]⊳\\trianglerightassign head start

5for

i←0​…​iout−1i\\leftarrow 0\\ldots i\_\{\\mathrm\{out\}\}\-1do⊳\\trianglerightshuffle used blocks into order of buckets

6

c←U⁡\[i\]c\\leftarrow U\[i\]
7

π′​\[S⁡\[c\]\]←π⁡\[i\]\\pi^\{\\prime\}\[S\[c\]\]\\leftarrow\\pi\[i\]
8

\(pnext,pend\)←B⁡\[c\]\(p\_\{\\mathrm\{next\}\},p\_\{\\mathrm\{end\}\}\)\\leftarrow B\[c\]
9if

pend=blockat​\(π⁡\[i\]\)\+bp\_\{\\mathrm\{end\}\}=\\textit\{blockat\}\(\\pi\[i\]\)\+bthen⊳\\trianglerightisπ⁡\[i\]\\pi\[i\]a partial block?

10

P⁡\[c\]←\(S⁡\[c\],pnext−blockat​\(π⁡\[i\]\)\)P\[c\]\\leftarrow\\bigl\(S\[c\],p\_\{\\mathrm\{next\}\}\-\\textit\{blockat\}\(\\pi\[i\]\)\\bigr\)
11endif

12

S⁡\[c\]←S⁡\[c\]\+1S\[c\]\\leftarrow S\[c\]\+1
13endfor

14

π′​\[iout\+σ​…​nπ−1\]←π⁡\[iout\+σ​…​nπ−1\]\\pi^\{\\prime\}\[i\_\{\\mathrm\{out\}\}\+\\sigma\\ldots n\_\{\\pi\}\-1\]\\leftarrow\\pi\[i\_\{\\mathrm\{out\}\}\+\\sigma\\ldots n\_\{\\pi\}\-1\]⊳\\trianglerightassign remaining blocks

15

f←iout\+σf\\leftarrow i\_\{\\mathrm\{out\}\}\+\\sigma
16

π←π′\\pi\\leftarrow\\pi^\{\\prime\}
17endprocedure

Algorithm 7restoration of the contents ofAA1procedureswap\(

i1,j1,i2,j2i\_\{1\},j\_\{1\},i\_\{2\},j\_\{2\}\)

2

π⁡\[i1\],π⁡\[i2\]←j2,j1\\pi\[i\_\{1\}\],\\pi\[i\_\{2\}\]\\leftarrow j\_\{2\},j\_\{1\}
3

π−1​\[j1\],π−1​\[j2\]←i2,i1\\pi^\{\-1\}\[j\_\{1\}\],\\pi^\{\-1\}\[j\_\{2\}\]\\leftarrow i\_\{2\},i\_\{1\}
4

i1,i2←i2,i1i\_\{1\},i\_\{2\}\\leftarrow i\_\{2\},i\_\{1\}
5endprocedure

6procedurefinalize

7

k←iP←0k\\leftarrow i\_\{P\}\\leftarrow 0
8

ifree,jfree←0,π⁡\[0\]i\_\{\\mathrm\{free\}\},j\_\{\\mathrm\{free\}\}\\leftarrow 0,\\pi\[0\]⊳\\trianglerightguaranteed unallocated as per invariant

9fori←0​…​nπ−1i\\leftarrow 0\\dots n\_\{\\pi\}\-1doπ−1​\[π​\[i\]\]←i\\pi^\{\-1\}\[\\pi\[i\]\]\\leftarrow i⊳\\trianglerightcompute inverse permutation

10for

iin←σ​…​σ\+⌊n/b⌋−1i\_\{\\mathrm\{in\}\}\\leftarrow\\sigma\\ldots\\sigma\+\\lfloor n/b\\rfloor\-1do⊳\\trianglerightpull first⌊n/b⌋\\lfloor n/b\\rfloorblocks

11

jin←π⁡\[iin\]j\_\{\\mathrm\{in\}\}\\leftarrow\\pi\[i\_\{\\mathrm\{in\}\}\]
12

jout,iout←iin−σ,π−1​\[iin−σ\]j\_\{\\mathrm\{out\}\},i\_\{\\mathrm\{out\}\}\\leftarrow i\_\{\\mathrm\{in\}\}\-\\sigma,\\pi^\{\-1\}\[i\_\{\\mathrm\{in\}\}\-\\sigma\]
13if

iin<iout<fi\_\{\\mathrm\{in\}\}<i\_\{\\mathrm\{out\}\}<fthen⊳\\trianglerightjoutj\_\{\\mathrm\{out\}\}not the right block, but allocated?

14copy block at

joutj\_\{\\mathrm\{out\}\}to block at

jfreej\_\{\\mathrm\{free\}\}⊳\\trianglerightpush block atjoutj\_\{\\mathrm\{out\}\}away

15

swap​\(iout,jout,ifree,jfree\)\\textsc\{swap\}\(i\_\{\\mathrm\{out\}\},j\_\{\\mathrm\{out\}\},i\_\{\\mathrm\{free\}\},j\_\{\\mathrm\{free\}\}\)⊳\\trianglerighttrack swapping ofifreei\_\{\\mathrm\{free\}\}andiouti\_\{\\mathrm\{out\}\}

16endif

17

l←sizeOfNextBlock​\(iin,iP\)l\\leftarrow\\textsc\{sizeOfNextBlock\}\(i\_\{\\mathrm\{in\}\},i\_\{P\}\)
18

A⁡\[k​…​k\+l−1\]←blockat​\(jin\)​\[0​…​l−1\]A\[k\\ldots k\+l\-1\]\\leftarrow\\textit\{blockat\}\(j\_\{\\mathrm\{in\}\}\)\[0\\ldots l\-1\]⊳\\trianglerightpull block fromjinj\_\{\\mathrm\{in\}\}

19

k←k\+lk\\leftarrow k\+l
20

swap​\(iout,jout,iin,jin\)\\textsc\{swap\}\(i\_\{\\mathrm\{out\}\},j\_\{\\mathrm\{out\}\},i\_\{\\mathrm\{in\}\},j\_\{\\mathrm\{in\}\}\)⊳\\trianglerighttrack swapping ofiini\_\{\\mathrm\{in\}\}andiouti\_\{\\mathrm\{out\}\}

21if

iin≠iouti\_\{\\mathrm\{in\}\}\\neq i\_\{\\mathrm\{out\}\}then⊳\\trianglerightdid we pull from elsewhere?

22

ifree,jfree←iin,jini\_\{\\mathrm\{free\}\},j\_\{\\mathrm\{free\}\}\\leftarrow i\_\{\\mathrm\{in\}\},j\_\{\\mathrm\{in\}\}⊳\\trianglerightif yes, there is now a new free spot

23endif

24endfor

25for

iin←σ\+⌊n/b⌋​…​f−1i\_\{\\mathrm\{in\}\}\\leftarrow\\sigma\+\\lfloor n/b\\rfloor\\ldots f\-1do⊳\\trianglerightpull any remaining blocks \(all inTT\)

26

jin←π⁡\[iin\]j\_\{\\mathrm\{in\}\}\\leftarrow\\pi\[i\_\{\\mathrm\{in\}\}\]
27

l←sizeOfNextBlock​\(iin,iP\)l\\leftarrow\\textsc\{sizeOfNextBlock\}\(i\_\{\\mathrm\{in\}\},i\_\{P\}\)
28

A⁡\[k​…​k\+l−1\]←blockat​\(jin\)​\[0​…​l−1\]A\[k\\ldots k\+l\-1\]\\leftarrow\\textit\{blockat\}\(j\_\{\\mathrm\{in\}\}\)\[0\\ldots l\-1\]⊳\\trianglerightpull block fromjinj\_\{\\mathrm\{in\}\}

29

k←k\+lk\\leftarrow k\+l
30endfor

31endprocedure

### 2\.4Finalisation

Following some rounds of sorting, the used elements of the blocks at logical indicesσ\\sigmatof−1f\-1hold the elements of the input in sorted order\. To wrap things up, we need to shuffle this collection of logically ordered blocks back intoAA, preserving the ordering\.

This is achieved by shuffling the used blocks around such that they are found at indices2​σ​…​2​σ\+f−12\\sigma\\ldots 2\\sigma\+f\-1, When the blocks are moved into place, their elements are actually moved right behind the elements of the previous block, ignoring gaps of unused elements at the end of partial blocks, and leaving the elements in their final locations\. After⌊n/b⌋\\lfloor n/b\\rfloorblocks have been shuffled to indices2​σ2\\sigmatonπ−1n\_\{\\pi\}\-1, all remaining blocks must be located inTT\. As an optimisation, we can directly copy them to their final positions without moving other data away, as there is no risk of other blocks being in their way\.

In this step, pairs of indicesi∘,j∘i\_\{\\circ\},j\_\{\\circ\}refer to the logical and physical indices of the same block\. These are related through the identity

j∘=π⁡\[i∘\],i∘=π−1​\[j∘\]\.j\_\{\\circ\}=\\pi\[i\_\{\\circ\}\],\\quad i\_\{\\circ\}=\\pi^\{\-1\}\[j\_\{\\circ\}\]\.\(4\)The following extra variables are used:

π−1\\pi^\{\-1\}inverse permutationarray ofnπn\_\{\\pi\}integersiin,jini\_\{\\mathrm\{in\}\},j\_\{\\mathrm\{in\}\}block holding next data chunkintegersiout,jouti\_\{\\mathrm\{out\}\},j\_\{\\mathrm\{out\}\}where next block goesintegersifree,jfreei\_\{\\mathrm\{free\}\},j\_\{\\mathrm\{free\}\}some unallocated blockintegersiPi\_\{P\}next partial blockintegerkkno\. of elements written intoAAinteger
To help us carry out these shuffles, we first compute an inverse permutationπ−1\\pi^\{\-1\}such thatπ∘π−1=\(0,…,nπ−1\)\\pi\\circ\\pi^\{\-1\}=\(0,\\ldots,n\_\{\\pi\}\-1\)\. The assignment of two blocksj1j\_\{1\}andj2j\_\{2\}can be swapped with theswap​\(i1,j1,i2,j2\)\\textsc\{swap\}\(i\_\{1\},j\_\{1\},i\_\{2\},j\_\{2\}\)function, while moving their contents in a separate step\.

For each block index2​σ≤j<nπ2\\sigma\\leq j<n\_\{\\pi\}, we distinguish three cases: \(a\) the right block is already at indexjj, in which case we just move it ahead to fill the gap, \(b\) the block is unallocated, in which case we find the right block atπ−1​\[j\]\\pi^\{\-1\}\[j\]and*pull*it to its final position, or \(c\) there is some other allocated block at indexjj, in which case we*push*it to some free block and reduce to case \(b\)\. The variablejfreej\_\{\\mathrm\{free\}\}tracks a spot to push an empty block to\. whenever this spot is occupied, we are going to subsequently pull a block from somewhere else, restoring the unallocated spot\.

## 3Complexity Analysis

For analysis, we assume array element size and alphabet sizeσ\\sigmato be constants andb∈𝒪⁡\(n\)b\\in\\mathcal\{O\}\(n\)\.

##### Runtime

Thesetup​\(A,n\)\\textsc\{setup\}\(A,n\)procedure from Alg\.[4](https://arxiv.org/html/2607.05302#alg4)initialises arrays ofnπ∈𝒪⁡\(n/b\)n\_\{\\pi\}\\in\\mathcal\{O\}\(n/b\)andbbelements respectively for a total runtime of𝒪⁡\(n/b\+b\)⊂𝒪⁡\(n\)\\mathcal\{O\}\(n/b\+b\)\\subset\\mathcal\{O\}\(n\)\. ThesortPhase​\(t\)\\textsc\{sortPhase\}\(t\)procedure from Alg\.[5](https://arxiv.org/html/2607.05302#alg5)executes one inner loop iteration for each of thennarray elements for a total runtime of𝒪⁡\(n\)\\mathcal\{O\}\(n\)\. ThefixupPhaseprocedure from Alg\.[6](https://arxiv.org/html/2607.05302#alg6)writes to each element ofπ′\\pi^\{\\prime\}once, for a total runtime of𝒪⁡\(n/b\)\\mathcal\{O\}\(n/b\)\. Lastly, thefinalizeprocedure from Alg\.[7](https://arxiv.org/html/2607.05302#alg7)traversesπ\\pi, copying⌊n/b⌋\\lfloor n/b\\rfloorblocks at most twice, while copying the remainingf−⌊n/b⌋−σ≤σf\-\\lfloor n/b\\rfloor\-\\sigma\\leq\\sigmablocks at most once for a total runtime of𝒪⁡\(b⁡\(n/b\+σ\)\)⊂𝒪⁡\(n\)\\mathcal\{O\}\(b\(n/b\+\\sigma\)\)\\subset\\mathcal\{O\}\(n\)\. For the full sort given in Alg\.[2](https://arxiv.org/html/2607.05302#alg2), thesortPhase​\(t\)\\textsc\{sortPhase\}\(t\)andfixupPhaseprocedures are called once for each key position, raising the total runtime to𝒪⁡\(nt​n\)\\mathcal\{O\}\(n\_\{t\}n\), which matches the conventional LSD radix sort\.

##### Space

In addition to the input stored inAA, arrayTToccupies𝒪⁡\(b\)\\mathcal\{O\}\(b\)bytes of storage, arraysUU,π\\pi,π′\\pi^\{\\prime\}, andπ−1\\pi^\{\-1\}occupy𝒪⁡\(n/b\)\\mathcal\{O\}\(n/b\)bytes of storage, and the other variables are of constant size, for a total extra storage requirement of𝒪⁡\(b\+n/b\)\\mathcal\{O\}\(b\+n/b\)\. This requirement can be minimised by choosing someb∈Θ⁡\(n\)b\\in\\Theta\(\\sqrt\{n\}\), giving just𝒪⁡\(n\)\\mathcal\{O\}\(\\sqrt\{n\}\)bytes of extra storage\.

## 4Implementation Notes and Optimisations

We have produced a variety of Radsort implementations in the C language\. We make our code freely available\.222See[https://github\.com/clausecker/radsort](https://github.com/clausecker/radsort)\.Implementations in other languages should be straightforward; in languages that do not provide pointers or pointer arithmetic, theBBarray can be altered to use the same data structure as thePParray, although at a loss of performance\.

Using the pointer\-based data structure of theBBarray for thePParray is not advisable: WhilefixupPhase\(Alg\.[6](https://arxiv.org/html/2607.05302#alg6)\) is simplified by not needing to translateBBintoPP, performance ofsortPhase​\(t\)\\textsc\{sortPhase\}\(t\)\(Alg\.[5](https://arxiv.org/html/2607.05302#alg5)\) is reduced due to the longer dependency chain when checking if the next input block is partial\. Additionally, thefinalizeprocedure \(Alg\.[7](https://arxiv.org/html/2607.05302#alg7)\) needs partial blocks to be tracked by their logical indices, mandating that the data structure be translated at least once\.

A radix ofσ=256\\sigma=256appears to be optimal for sorting key\-value pairs of integers, though the best choice depends on cache and element size, as well as key type and distribution and must be determined empirically\. The basic tradeoff is that larger radices allow for less rounds of sorting, but require more cache to track the buckets, eventually escalating into higher levels of the cache hierarchy at great performance cost\.

For a type\-generic implementation, thesortPhase​\(t\)\\textsc\{sortPhase\}\(t\)procedure should be monomorphised, ideally for each key position\. In the common case of integer keys split into byte\-sized digits \(σ=256\\sigma=256\), retrieving the byte at positionttof the key by means of shifts and masks performs better than direct loading of the key byte from memory, as the latter generally incurs an extra memory operation, while the former shifts the load to arithmetic execution units that are otherwise underutilised\.

While a block size of orderb∈Θ⁡\(n\)b\\in\\Theta\(\\sqrt\{n\}\)minimises space overhead, choosing a fixed block size simplifies the implementation\. The space overhead then becomes some fixed fraction of the input size, which is often acceptable\. For example, we usedb=512b=512in our experiments\. At an element size of8​B8\\,\\textrm\{B\}, this gives a fixed overhead of2​MiB2\\,\\textrm\{MiB\}for theTTarray,8​KiB8\\,\\textrm\{KiB\}for arraysBB,CC,PP, andSStogether, and9​B9\\,\\textrm\{B\}per block of input \(i\. e\. 0\.22 % the input data set\) to storeπ\\pi,π′\\pi^\{\\prime\},π−1\\pi^\{\-1\}, andUU, reusingπ′\\pi^\{\\prime\}forπ−1\\pi^\{\-1\}\.

Like with a conventional LSD radix sort, performance of Radsort is memory bound and eliminating just a single load or store from the inner loop ofSortPhase​\(t\)\\textsc\{SortPhase\}\(t\)can improve performance dramatically, as those memory accesses compete for resources with the write in Alg\.[5](https://arxiv.org/html/2607.05302#alg5), L\.[17](https://arxiv.org/html/2607.05302#alg5.l17), the bottleneck of the algorithm\. One technique to eliminate such a load is discussed in §[4\.2](https://arxiv.org/html/2607.05302#S4.SS2), and gives a 5–50 % speedup depending on array size and microarchitecture\.

### 4\.1Avoiding Finalisation

In some use cases, it suffices to iterate over the sorted data\. In such cases, thefinalizeprocedure can be skipped, employing instead thesizeOfNextBlock​\(i,iP\)\\textsc\{sizeOfNextBlock\}\(i,i\_\{P\}\)function from Alg\.[3](https://arxiv.org/html/2607.05302#alg3)to iterate over the array elements in permuted representation as in Alg\.[5](https://arxiv.org/html/2607.05302#alg5)\(see also §[2\.1](https://arxiv.org/html/2607.05302#S2.SS1)\)\.

Using such an iterator pattern, the dataset can be modified in place, and even be truncated \(cf\. Alg\.[8](https://arxiv.org/html/2607.05302#alg8)\) without disturbing the Radsort algorithm state\. This enables algorithm designs where the same dataset is repeatedly sorted by different keys, and then modified, without having to finalise the sort or recreate the algorithm state for every sort\.

Algorithm 8truncate the dataset such that logical blockii, elementjjis the final datumproceduretruncate\(

i,ji,j\)

f←i\+1f\\leftarrow i\+1⊳\\trianglerightmark logical blockiias the final block

for

iP←0​…​σ−1i\_\{P\}\\leftarrow 0\\ldots\\sigma\-1do⊳\\trianglerightskip partial blocks before logical blockii

\(iiP,liP\)←P⁡\[iP\]\(i\_\{i\_\{P\}\},l\_\{i\_\{P\}\}\)\\leftarrow P\[i\_\{P\}\]

if

i≤iiPi\\leq i\_\{i\_\{P\}\}then

P⁡\[iP\]←\(i,j\+1\)P\[i\_\{P\}\]\\leftarrow\(i,j\+1\)⊳\\trianglerightlogical blockiiends at elementjj

return

endif

endfor

endprocedure

### 4\.2Faster End\-of\-Block Checking

ThesortPhase​\(t\)\\textsc\{sortPhase\}\(t\)procedure of Alg\.[5](https://arxiv.org/html/2607.05302#alg5)tracks output blocks through pairs\(pnext,pend\)\(p\_\{\\mathrm\{next\}\},p\_\{\\mathrm\{end\}\}\)pointing to the next free element and the end of the block, respectively\. It is beneficial to reduce this pair to just one pointer, as to eliminate a load of the second pointer in the hot loop, and to reduce cache pressure\.

In a simplified model, where each element is one byte in size and the allocation of all data structures includingAAcan be controlled, this can be achieved by choosing some block sizeb=2qb=2^\{q\}and ensuring bothAAandTTare aligned to a multiple ofbb\. Then, nopendp\_\{\\mathrm\{end\}\}is needed, as a pointer into some block ofAAorTTpoints to a block boundary iff it is aligned to a multiple ofbbbytes\-\-\-which can easily be checked with a bitwise\-and on the pointer\.333Thepend=blockat​\(π⁡\[i\]\)\+bp\_\{\\mathrm\{end\}\}=\\textit\{blockat\}\(\\pi\[i\]\)\+bcheck of Alg\.[6](https://arxiv.org/html/2607.05302#alg6), L\.[9](https://arxiv.org/html/2607.05302#alg6.l9)must then be realised with a range check likeblockat​\(π⁡\[i\]\)≤pnext<blockat​\(π⁡\[i\]\)\+b\\textit\{blockat\}\(\\pi\[i\]\)\\leq p\_\{\\mathrm\{next\}\}<\\textit\{blockat\}\(\\pi\[i\]\)\+b\. This is not strictly legal in the C23 language\[[6](https://arxiv.org/html/2607.05302#bib.bib8), § 6\.5\.8 ¶ 6\], \(pointers may only be compared for ordering if they point into the same array, butpnextp\_\{\\mathrm\{next\}\}may point into eitherAAorTT\) but unproblematic in practice\.

In real applications, we usually cannot control the alignment ofAAand neither is the element size guaranteed to be some nice number\. However, with some modifications, the same idea still works: let the element size be2r​e2^\{r\}e, witheeodd and the block size again beb=2qb=2^\{q\}\. Then in an array of such elements of no particular alignment, everybbelements there is an element with addressp0p\_\{0\}such that

p0≫r≡0modb,p\_\{0\}\\gg r\\equiv 0\\mod b,\(5\)where≫\\ggdenotes a bitwise logical right\-shift\. This condition is just as easy to check using a bitwise\-and on the pointer and allows us to find block boundaries if the first element of each block satisfies Eq\.[5](https://arxiv.org/html/2607.05302#S4.E5)\. Furthermore, for the address of an arbitrary array elementpp, we can find the offset0≤o<b0\\leq o<bof such an element fromppas

o=−\(p≫r\)​e−1modbo=\-\(p\\gg r\)e^\{\-1\}\\bmod b\(6\)wheree−1e^\{\-1\}is the modular inverse ofeemodulobb\. Due tobbbeing a power of two,e−1e^\{\-1\}is easy to compute\[[4](https://arxiv.org/html/2607.05302#bib.bib5)\]\. This gives the following modifications: allocate2​σ\+12\\sigma\+1blocks444An extra scratch block is needed as the total number of elements shunted to head and tail may be up to2​b−22b\-2, possibly exceeding a whole block of elements\.forTTaligned to a multiple of2q​r2^\{qr\}, ensuring that its block address satisfy Eq\.[5](https://arxiv.org/html/2607.05302#S4.E5)\. Compute the offsetooof the first element ofAAsatisfying Eq\.[5](https://arxiv.org/html/2607.05302#S4.E5)using Eq\.[6](https://arxiv.org/html/2607.05302#S4.E6)\. Theblockat​\(i\)\\textit\{blockat\}\(i\)function is adjusted to

blockat​\(i\)=\{address ofT⁡\[i\]if0≤i<2​σ\+1address ofA⁡\[\(i−2​σ−1\)​b\+o\]if2​σ\+1≤i<nπ\\textit\{blockat\}\(i\)=\\begin\{cases\}\\mbox\{address of $T\[i\]$\}&\\mbox\{if $0\\leq i<2\\sigma\+1$\}\\\\ \\mbox\{address of $A\[\(i\-2\\sigma\-1\)b\+o\]$\}&\\mbox\{if $2\\sigma\+1\\leq i<n\_\{\\pi\}$\}\\end\{cases\}\(7\)andnπ=⌊\(n−o\)/b⌋\+2​σ\+1n\_\{\\pi\}=\\lfloor\(n\-o\)/b\\rfloor\+2\\sigma\+1\. Transfer the firstooelements to a scratch block and adjust Alg\.[4](https://arxiv.org/html/2607.05302#alg4)such that this scratch block \(the head\) ends up at logical blockσ\\sigma, followed by the blocks overlayingAA, following the tail\. Various other steps of the algorithm must also be adjusted to account for the extra scratch block\. The end\-of\-block check then becomes a check for the condition of Eq\.[5](https://arxiv.org/html/2607.05302#S4.E5), which is realised by masking the pointer with a bitmask of the form\(b−1\)≪r\(b\-1\)\\ll rand checking if the result is zero\. Thependp\_\{\\mathrm\{end\}\}pointers no longer need to be read and can be eliminated entirely\.

This approach yielded a 5–50 % speedup depending on array size and microarchitecture\. It should be considered if the programming language and environment permit it\.

### 4\.3Parallel Operation

The algorithm can be parallelised with some changes to the data structures\. Before each sort phase, we distribute the blocks of the input array intontn\_\{t\}roughly evenly sized chunks, for sorting withntn\_\{t\}threads\. Each chunk gets its own head start ofσ\\sigmablocks, requiring aTTarray of2​σ​nt2\\sigma n\_\{t\}scratch blocks in total\. The sort phase then proceeds in parallel with each thread sorting its chunk with its ownBBarray\. The fixup phase is sequential and must interleave the individual thread’sBBarrays into one globalPParray ofnt​σn\_\{t\}\\sigmaentries each\.

The finalisation procedure is more troublesome to parallelise\. It can be left as a sequential operation, as it only contributes a small amount of the total run time, or it can be implemented by first shuffling all the blocks into order using a parallel permutation algorithm\[[3](https://arxiv.org/html/2607.05302#bib.bib1)\], and then moving the array elements to eliminate gaps caused by partial blocks\.

As Radsort is memory bound and exhibits only moderate cache locality, a performance ceiling quickly is reached as the available memory channels are saturated\. See §[5\.2](https://arxiv.org/html/2607.05302#S5.SS2)for more discussion\. It may be of interest to use an initial round of sorting on the most significant key position to split the input into buckets that can be parcelled into threads, each of which processes the remaining key positions of its bucket\(s\) from the least\-significant digit sequentially\. This reduces the working set of each thread fromnnton/ntn/n\_\{t\}elements on average, improving cache efficacy and reducing NUMA effects\.

## 5Evaluation

We evaluated the performance of Radsort in comparison to a classic out\-of\-place LSD radix sort on a variety of systems\. Five algorithm variants are compared in total:

genericA generic out\-of\-place LSD radix sort \(Alg\.[1](https://arxiv.org/html/2607.05302#alg1)\) with output prefetching\[[2](https://arxiv.org/html/2607.05302#bib.bib6)\]\. An initial pass takes histograms of the 4 key bytes, followed by 4 rounds of sorting\.

swcAn out\-of\-place LSD radix sort implemented with Wassenberg’s software\-defined write\-combining\[[8](https://arxiv.org/html/2607.05302#bib.bib4)\]with a block size of512​B512\\,\\mathrm\{B\}\.

radsortSingle\-threaded Radsort \(Alg\.[2](https://arxiv.org/html/2607.05302#alg2)\) as described in §[2](https://arxiv.org/html/2607.05302#S2)with none of the improvements mentioned in §[4](https://arxiv.org/html/2607.05302#S4)using a fixed block size ofb=512b=512\(i\. e\.4​KiB4\\,\\mathrm\{KiB\}\)\.

bitmanipSingle\-threaded Radsort implemented using the bit\-manipulation based end\-of\-block checking described in §[4\.2](https://arxiv.org/html/2607.05302#S4.SS2), otherwise the same as*radsort*\.

parallelMulti\-threaded Radsort using 4 threads by default\. Otherwise the same as*radsort*\.

### 5\.1Setup

All benchmarks sort pairs of 32\-bit keys and 32\-bit values innt=4n\_\{t\}=4rounds by each key byte in turn \(σ=256\\sigma=256\)\. The key values are uniformly distributed over the range0​…​232−10\\ldots 2^\{32\}\-1, as generated by a xorshift RNG\[[5](https://arxiv.org/html/2607.05302#bib.bib7)\]\. While we have evaluated radsort on a variety of platforms, we have selected benchmark results as measured on the following machines for this paper:

Two sets of benchmarks were performed: in the*size*sets \(Fig\.[4](https://arxiv.org/html/2607.05302#S5.F4)\), all algorithms were measured on arrays with2n2^\{n\}and3×2n−13\\times 2^\{n\-1\}elements with array sizes from22to as many elements as the memory fits\. In the*threads*sets \(Fig\.[3](https://arxiv.org/html/2607.05302#S5.F3)\), the*parallel*implementation was measured with an array size of128​GiB128\\,\\mathrm\{GiB\}\(i\. e\.2342^\{34\}elements\) and thread counts from 1 to 128\.

20 runs of each benchmark were performed\. In the plots, dots show the individual runs, with a line drawn through the average\. For the size sets, execution was pinned to a single socket to avoid NUMA effects as much as possible\.

Figure 3:Sorting speed compared on a128​GiB128\\,\\mathrm\{GiB\}array with thread count from 1 to 128\.Figure 4:Sorting speeds on three machines with arrays of8​B8\\,\\mathrm\{B\}to96​GiB96\\,\\mathrm\{GiB\}\(192​GiB192\\,\\mathrm\{GiB\}bottom\)\.
### 5\.2Results

For small arrays up to around twice the L2 cache size, the*generic*implementation is the clear winner\. As cache misses are mostly absent at these array sizes, the overhead of Radsort cannot be outweighed by its cache\-locality benefits\. As input size exceeds this limit, Radsort and its variants quickly outperform the generic implementation\. Implementing bit\-manipulation based end\-of\-block checking yields a speedup of 5–50 % on all systems tested \(including others not shown here\), except for the*power*system where it is slightly slower\.

As expected, the*swc*variant performs very consistently over the whole range of array lengths, as it minimises cache effects through manual write combining\. This consistency is maintained even as the array size approaches a significant chunk of the main memory, while the generic implementation drops to around half of its performance, matching the results observed by Wassenberg et al\.\[[8](https://arxiv.org/html/2607.05302#bib.bib4)\]\. The*bitmanip*variant exceeds this performance on all systems for array sizes above the L2 cache size, showing that its timely reuse of input blocks for output blocks avoid read\-to\-own transfers just as effectively as*swc*\.

Whenever it is worth using Radsort over an out\-of\-place radix sort, the parallel variant outperforms the scalar variants, though the advantages quickly diminish \(cf\. Fig\.[3](https://arxiv.org/html/2607.05302#S5.F3)\), topping out around 24 threads on the*icelake*machine and around 48 threads on the*power*machine\. NUMA heavily impacts the*grace*machine, causing high measurement variance based on thread placement\. Performance peaks here at around 72 threads\.

In summary, it seems advantageous to use*generic*for short arrays, and to switch to*bitmanip*or*parallel*as input exceeds some empirically determined threshold\. This threshold could be set to the size ofTT, allowing for the reuse ofTTas the output arrayA′A^\{\\prime\}of Alg\.[1](https://arxiv.org/html/2607.05302#alg1)\.

## 6Related Work

Radsort is an attempt to adapt the ideas of IPS2RA\[[1](https://arxiv.org/html/2607.05302#bib.bib2)\]to LSD radix sort\. While IPS2RA is designed to use𝒪⁡\(1\)\\mathcal\{O\}\(1\)extra memory, the author found this difficult to achieve under the LSD approach’s stability requirements\. By raising the overhead to𝒪⁡\(n\)\\mathcal\{O\}\(\\sqrt\{n\}\), the permutation of blocks can be explicitly tracked inπ\\pi, addressing these challenges\. As a consequence, blocks only have to actually be moved around during finalisation, giving a significant speedup\.

The result is similar to the “out\-of\-cache, in\-place, list of blocks” partitioning described by Polychroniou et al\.\[[7](https://arxiv.org/html/2607.05302#bib.bib3), § 3\.2\.3\], but uses arrays over linked lists to trackπ\\pi, improving cache utilisation, and permitting a straightforward finalisation procedure\.

A major performance limitation of radix sorts is the need for output buckets to be read into cache so that they can be written to \(read for ownership\), effectively halving the write bandwidth and causing cache misses\. While CPUs provide write\-combining buffers to avoid this problem, radix sorts have a fan\-out that exceeds their number, rendering them ineffective\. The performance impact can be reduced by prefetching the output buckets speculatively\[[2](https://arxiv.org/html/2607.05302#bib.bib6)\], or more thoroughly by manual software\-defined write\-combining\[[8](https://arxiv.org/html/2607.05302#bib.bib4)\]\.

Radsort solves the problem more elegantly: as input blocks are reused for output blocks after only a short delay, the output blocks are still hot in cache555As long as the cache fits up toσ​b\\sigma belements that have been consumed but not yet reused as output blocks\.when written to, avoiding an extra read for ownership\. Consequently, we have found prefetching to have little to no performance impact and software write\-combining to not be needed\.

## 7Future Work

It is promising to adapt Radsort to the sorting of dissimilarly sized elements, such as lines of text or JSON records, directly, i\. e\. without resorting to sorting an array of pointers to variable\-length records\. Likewise, adapting the general approach to external sorting should be investigated\.

During the research for this paper, the authors investigated whether Radsort’s data structures would work for an MSD \(i\. e\. recursive\) radix sort procedure\. Keeping the data structures identical, we quickly run into problems: while each recursive iteration consumes one partial block from the bucket it recurses over, it produces up toσ\\sigmanew partial blocks, leading to an uncompetitive memory overhead\. The authors found that this problem can be addressed by compacting the partial blocks into a sequence of dense blocks after each iteration, tracking the beginning and end of partial data within this sequence using an extra data structure\. This reduces the extra space for partial blocks to one block per recursion level plus2​σ2\\sigmawords to track the number of elements in the partial tail of each bucket, requiring𝒪⁡\(\(σ\+b\)​nt\)\\mathcal\{O\}\\bigl\(\(\\sigma\+b\)n\_\{t\}\\bigr\)extra space\. Due to time constraints, this idea was not explored further\.

## 8Conclusion

In Radsort, we provide a straightforward way to implement LSD radix sorting with just𝒪⁡\(n\)\\mathcal\{O\}\(\\sqrt\{n\}\)space overhead\. Radsort readily parallelises to a moderate number of threads\. The performance can be further improved using bit manipulation techniques, at the cost of a more complex implementation\.

As a consequence of more effective cache utilisation, the performance of Radsort is competitive with a standard out\-of\-place LSD radix sort for arrays as small as2​MiB2\\,\\mathrm\{MiB\}, even when using only one thread\. The read\-to\-own bottleneck on large datasets is avoided without the requirement of platform\-specific techniques like software\-defined write\-combining\.

## References

- \[1\]M\. Axtmann, S\. Witt, D\. Ferizovic, and P\. Sanders\(2022\)Engineering in\-place \(shared\-memory\) sorting algorithms\.ACM Trans\. Parallel Comput\.9\(1\),pp\. 1–62\.External Links:[Document](https://dx.doi.org/10.1145/3505286)Cited by:[§6](https://arxiv.org/html/2607.05302#S6.p1.1)\.
- \[2\]T\. Downs\(2019\)Beating up on qsort\.External Links:[Link](https://travisdowns.github.io/blog/2019/05/22/sorting.html)Cited by:[item generic](https://arxiv.org/html/2607.05302#S5.I1.ix1.p1.1),[§6](https://arxiv.org/html/2607.05302#S6.p3.1)\.
- \[3\]T\. Hagerup and J\. Keller\(1995\)Fast parallel permutation algorithms\.Parallel Process\. Lett\.5\(2\),pp\. 139–148\.External Links:[Document](https://dx.doi.org/10.1142/S0129626495000126)Cited by:[§4\.3](https://arxiv.org/html/2607.05302#S4.SS3.p2.1)\.
- \[4\]J\. Hurchalla\(2022\)An improved integer modular multiplicative inverse \(modulo2w2^\{w\}\)\.External Links:2204\.04342Cited by:[§4\.2](https://arxiv.org/html/2607.05302#S4.SS2.p3.3)\.
- \[5\]G\. Marsaglia\(2003\)Random number generators\.J\. Mod\. Appl\. Stat\. Methods\.2\(1\),pp\. 2–13\.External Links:[Document](https://dx.doi.org/10.22237/jmasm/1051747320)Cited by:[§5\.1](https://arxiv.org/html/2607.05302#S5.SS1.p1.1)\.
- \[6\]J\. Meneide and F\. Wiedijk\(2024\)Programming languages – C\.StandardTechnical ReportISO/IEC 9899:2023,International Organization for Standardization\.Cited by:[footnote 3](https://arxiv.org/html/2607.05302#footnote3)\.
- \[7\]O\. Polychroniou and K\. A\. Ross\(2014\)A comprehensive study of main\-memory partitioning and its application to large\-scale comparison\- and radix\-sort\.InProceedings of the 2014 ACM SIGMOD International Conference on Management of Data,SIGMOD ’14,pp\. 755–766\.External Links:ISBN 9781450323765,[Link](https://doi.org/10.1145/2588555.2610522),[Document](https://dx.doi.org/10.1145/2588555.2610522)Cited by:[§6](https://arxiv.org/html/2607.05302#S6.p2.1)\.
- \[8\]J\. Wassenberg and P\. Sanders\(2010\)Faster radix sort via virtual memory and write\-combining\.External Links:1008\.2849Cited by:[item swc](https://arxiv.org/html/2607.05302#S5.I1.ix2.p1.1),[§5\.2](https://arxiv.org/html/2607.05302#S5.SS2.p2.1),[§6](https://arxiv.org/html/2607.05302#S6.p3.1)\.