# Parallel Integer Sort: Theory and Practice

Xiaojun Dong  UC Riverside<EMAIL_ADDRESS>,  Laxman Dhulipala  University of
Maryland<EMAIL_ADDRESS>,  Yan Gu  UC Riverside<EMAIL_ADDRESS>and  Yihan Sun
UC Riverside<EMAIL_ADDRESS>
(2024)

###### Abstract.

Integer sorting is a fundamental problem in computer science. This paper
studies parallel integer sort both in theory and in practice. In theory, we
show tighter bounds for a class of existing practical integer sort algorithms,
which provides a solid theoretical foundation for their widespread usage in
practice and strong performance. In practice, we design a new integer sorting
algorithm, DovetailSort, that is theoretically-efficient and has good
practical performance.

In particular, DovetailSort overcomes a common challenge in existing parallel
integer sorting algorithms, which is the difficulty of detecting and taking
advantage of duplicate keys. The key insight in DovetailSort is to combine
algorithmic ideas from both integer- and comparison-sorting algorithms. In our
experiments, DovetailSort achieves competitive or better performance than
existing state-of-the-art parallel integer and comparison sorting algorithms
on various synthetic and real-world datasets.

Integer Sort, Radix Sort, Sorting Algorithms, Parallel Algorithms

††journalyear: 2024††copyright: rightsretained††conference: The 29th ACM
SIGPLAN Annual Symposium on Principles and Practice of Parallel Programming;
March 2–6, 2024; Edinburgh, United Kingdom††booktitle: The 29th ACM SIGPLAN
Annual Symposium on Principles and Practice of Parallel Programming (PPoPP
’24), March 2–6, 2024, Edinburgh, United Kingdom††doi:
10.1145/3627535.3638483††isbn: 979-8-4007-0435-2/24/03††ccs: Theory of
computation Parallel algorithms††ccs: Theory of computation Shared memory
algorithms††ccs: Theory of computation Sorting and searching

## 1 Introduction

Sorting is one of the most widely-used primitives in algorithm design, and has
been extensively studied. For many if not most cases, the keys to be sorted
are fixed-length integers. Sorting integer keys is referred to as the _integer
sort_ problem. An integer sorting algorithm takes as input $n$ records with
integer keys in the range $0$ to $r-1$, and outputs the records with keys in
non-decreasing order. Despite decades of effort in studying integer sorting
algorithms, however, obtaining _parallel integer sorting (IS) algorithms that
are efficient both in theory and in practice_ has remained elusive.

_Theoretical Challenges._ _As a special type of sorting, integer sorting
algorithms can outperform comparison sorting algorithms by using the integer
encoding of keys._ This claim is verified in many existing studies (Obeya et
al., 2019; Axtmann et al., 2022) as well as our experiments (see Fig. 1). As a
result, in real-world applications (e.g., (Dhulipala et al., 2021; Sun et al.,
2019; Wang et al., 2023)), integer sort is usually preferred (instead of
comparison sort) when the keys are integers. While it is not surprising that
IS algorithms can outperform comparison sorts, we observe a _significant gap
in connecting the high performance with theory_. Theoretical parallel IS
algorithms with good bounds (Bhatt et al., 1991; Albers and Hagerup, 1997;
Andersson et al., 1998; Hagerup, 1991; Matias and Vishkin, 1991) are quite
complicated, and we are unaware of any implementations of them. On the other
hand, for the practical parallel IS solutions (Blelloch et al., 2020a; Obeya
et al., 2019; Axtmann et al., 2022), we are unaware of “meaningful” analysis
to explain their good performance for general $r$. The best-known bounds for
them, as discussed in (Obeya et al., 2019), are $O(n\log r)$ work (number of
operations) and $\operatorname*{polylog}(rn)$ span (longest dependence chain).
However, note that the main use case for integer sort is when $r=\Omega(n)$,
since otherwise the simpler counting sort (Vishkin, 2010; Rajasekaran and
Reif, 1989) can be used. In this case, the bounds for practical parallel IS
algorithms are no better than comparison sorts ($O(n\log n)$ work and
polylogarithmic span (Blelloch et al., 2020b; Goodrich and Jacob, 2023;
Blelloch et al., 2010, 2020a)). This leads to the following open question in
theory: do the practical parallel IS algorithms indeed have lower asymptotic
costs than comparison sort (and if so, under what circumstances)?

_Practical Challenges._ _As a special type of sorting, integer sorting
algorithms should outperform comparison sorting algorithms by using the
integer encoding of keys._ Since integers are comparable, comparison sort can
be considered as a baseline for sorting integers. Unfortunately, SOTA parallel
IS algorithms do not consistently outperform comparison-based sorting
algorithms. One key reason is the inherent difficulty of dealing with
_duplicate keys_ in integer sorts. In principle, duplicate keys are beneficial
for sorting algorithms. For example, samplesort can skip a recursive
subproblem between two equal pivots; similarly, quicksort can separate keys
equal to the pivot to avoid further processing them. Interestingly, such a
case does not apply to integer sort. Existing parallel IS implementations
follow the most-significant digit (MSD) framework that partitions all keys
into buckets based on the _integer encoding_ (i.e., 8–12 highest bits), and
recurses within each bucket. As such, equal keys cannot be detected until the
last recursion. Although some techniques can be used to detect special
distributions (e.g., all keys are the same), to the best of our knowledge, no
existing parallel IS implementation can benefit from duplicate keys in a
general (provable) and non-trivial manner111 Some techniques in existing IS
implementations (Obeya et al., 2019; Axtmann et al., 2022) have a side-effect
to benefit from duplicates in some cases, at a cost of making the algorithms
_unstable_. We want to overcome this issue without sacrificing stability. .
Fig. 1 and 3 show that on inputs with heavy duplicates, SOTA parallel IS
algorithms can be slower than comparison sorts.

_Our Contributions._ In this paper, we study parallel integer sort to overcome
the aforementioned challenges both in theory and in practice. In theory, _we
answer the open question by proving that a class of parallel IS algorithms,
including DovetailSort proposed in this paper, have $O(n\sqrt{\log r})$ work
for input size $n$ and key range $0$ to $r-1$_. This is asymptotically better
than comparison sorts for a large key range $r=n^{o(\log n)}$. This explains
why the existing parallel IS algorithms outperform comparison sorts in
practice in many cases. We also show that we can achieve
$\operatorname*{polylog}(nr)$ span _whp_ for an unstable integer sort, and
$\tilde{O}(2^{\sqrt{\log r}})$ span for a class of practical stable parallel
IS algorithms, such as our new DovetailSort and the integer sort in ParlayLib
(Blelloch et al., 2020a). We also prove that DovetailSort has $\Theta(n)$ work
for some special input distributions with heavy duplicates. In addition to the
new bounds, our theoretical results (Sec. 4) also explain the choices of
parameters for many (if not all) practical parallel IS implementations. We
believe that our theoretical contribution is important both for analyzing our
own algorithm and for understanding parallel integer sort as a general
algorithmic problem.

Figure 1. Heatmap to compare sorting algorithms on $10^{9}$ records with
32-bit keys and 32-bit values. All numbers are running times relative to the
best for each input. Raw data are in Tab. 3. The baseline algorithms are
described in Tab. 2.

In practice, _we propose DovetailSort (DTSort), a parallel integer sorting
algorithm that is efficient and stable_. DTSort combines the algorithmic ideas
from both integer- and comparison-sorting algorithms. The goal is to detect
and take advantage of heavily duplicate keys (_heavy keys_) as in samplesort
(a commonly-used comparison sort), but preserve the low asymptotic work bound
as in integer sort. The key idea are to 1) use sampling to identify heavy
keys, separate them from the other keys (called _light keys_), and skip them
in the recursive calls, and 2) exploit the bitwise encoding of the keys to
split light keys into subproblems and deal with them recursively. Although the
high-level idea of _separating_ heavy and light keys is natural, an
interesting challenge is that all keys need to be put in increasing order in
the output, and thus the heavy and light keys have to be finally _interleaved_
with each other. Therefore, the heavy-light separation idea requires
_dovetailing_ the keys at a low cost. To do this, DTSort first carefully puts
the heavy keys close to their final destinations when they are identified.
After the light keys have been sorted from the recursive calls, DTSort uses a
dedicated algorithm to combine them into the final sorted order. We present
more details in Sec. 3.

Due to the algorithmic improvements, DTSort achieves consistently strong
performance on various synthetic and real-world input instances. We evaluate
DTSort using SOTA parallel sorting algorithms as baselines. We illustrate the
results on 32-bit keys in Fig. 1, and present more in Sec. 6. Compared to the
best baseline on 32-bit keys (PLIS), our algorithm is competitive
(1.02$\times$ faster) for no or light duplicate keys, and can be up to
2.3$\times$ faster with more duplicates. The advantage of DTSort is larger on
64-bit keys (see Tab. 3). We also evaluate DTSort in several applications with
real-world input data. In tens of the input instances we tested, DTSort is the
fastest in most cases. We note that sorting is one of the most well-studied
problems, and it is almost impossible for one algorithm to be the best for
_all_ input instances. The goal is to design algorithms that are always close
to best in any use case. We believe DTSort achieves very good overall
performance on a variety of input patterns. Our code is publicly available
(Dong et al., 2023a).

## 2 Preliminaries

### 2.1 Notations

We use $O(f(n))$ _with high probability_ (_whp_) in $n$ to mean $O(cf(n))$
with probability at least $1-n^{-c}$ for $c\geq 1$. Throughout the paper, we
consider the high probability bound in $n$ for input size $n$. We omit “in
$n$” with clear context. We use $[r]$ to refer to integers in the range $0$ to
$r-1$. In our analysis, we use the term $2^{\sqrt{\log r}}$, which satisfies
$\omega(\log^{c_{1}}r)=2^{\sqrt{\log r}}=o(r^{c_{2}})$ for any positive
constants $c_{1}$ and $c_{2}$. We use $\tilde{O}(f(nr))$ to denote
${O}(f(nr)\cdot\operatorname*{polylog}(nr))$. We use $\log n$ to represent
$\log_{2}n$. More notations are presented in Tab. 1.

$A[1..n]$ |  Original input array with size $n$
---|---
$[r]$ |  The range of the integer keys $0,\dots,r-1$
$\gamma$ | Number of bits in a “digit” (sorted by each level)
$b=2^{\gamma}$ | The “radix” size
$\theta$ | Base case size threshold
$n^{\prime}$ | Current (recursive) problem size
$d$ | Number of remaining “digits” to be sorted
Table 1. Notations and parameters in this paper.

### 2.2 Computational Models

We use the work-span (or work-depth) model for fork-join parallelism with
binary forking to analyze parallel algorithms (Cormen et al., 2009; Blelloch
et al., 2020b). We assume a set of threads that share a common memory. A
process can fork two child software threads to work in parallel. When both
children complete, the parent process continues. The _work_ of an algorithm is
the total number of instructions and the _span_ (depth) is the length of the
longest sequence of dependent instructions in the computation. We can execute
the computation using a randomized work-stealing scheduler in practice in
$W/P+O(D)$ time _whp_ on $P$ processors (Blumofe and Leiserson, 1998; Arora et
al., 2001; Gu et al., 2022).

### 2.3 Sorting and Integer Sorting

Integer sort takes as input an array $A$ of $n$ records with a key function
that maps each record to an integer key in range $[r]$. The goal is to reorder
the records to be in non-decreasing order. As mentioned in Sec. 1, the main
and most general use scenario of integer sorts is when $r=\Omega(n)$, so we
assume this in our theoretical analysis unless otherwise specified. Many
existing integer sort algorithms consider the integer keys as radix-$b$
(base-$b$) numbers, and work on each “digit” in a round. Usually each digit
refers to every $\mathit{\gamma}$ bits in the integer and
$b=2^{\mathit{\gamma}}$. These algorithms are also referred to as radix sort.
The classic sequential radix sort uses $b=\Theta(n)$, starts from the least
significant digits (LSD), and reorders all the records on every digit using a
stable sort, until all bits are used up (Cormen et al., 2009). It takes
$O(n\log_{n}r)$ work.

Unlike the sequential setting, existing parallel integer sort implementations
start from the most significant digit (MSD). The algorithms (see Alg. 1) have
two steps: 1) _distributing_ all records into $b=2^{\mathit{\gamma}}$ buckets
by their MSD, and 2) _recursing_ to sort all buckets in parallel. This top-
down approach allows all subproblems to be solved independently in parallel,
which is more I/O friendly. The distribution step can be performed by a
counting sort (see Sec. 2.4) with the bucket id as the MSD. Most practical
parallel integer sorting algorithms (Blelloch et al., 2020a; Obeya et al.,
2019; Axtmann et al., 2022; Kokot et al., 2017; Cho et al., 2015; Blelloch et
al., 1991; Lee et al., 2002; Polychroniou and Ross, 2014; Solomonik and Kale,
2010; Wassenberg and Sanders, 2011; Zhang and Wu, 2012) follow the framework
in Alg. 1 and have good performance. However, if a tight bound exists for
these parallel MSD algorithms remains open.

For simplicity, we assume to store the sorted result in the input array (“in-
place”). Here “in-place” refers to the algorithm interface. Our algorithm
still uses $O(n)$ extra space. A sorting algorithm is _stable_ if the output
preserves the relative order among equal keys as in the input, and unstable
otherwise. Stability is required in many sorting applications.

### 2.4 Counting Sort (aka. Distribution)

Counting sort is a commonly-used primitive in parallel sorting algorithms to
distribute records into (a small number of) buckets. It takes an input of
$n^{\prime}$ records with a function to map each record to an integer (the
_bucket id_) in $[r^{\prime}]$ ($r^{\prime}\leq n^{\prime}$). The goal is to
sort the records by the bucket ids. The classic stable counting sort (Vishkin,
2010), which is widely used in practice, has $O(n^{\prime})$ work and
$O(r^{\prime}+\log n^{\prime})$ span. Due to page limits, we put more
discussions about this algorithm in Appendix B, but the audience can treat the
counting sort as a black box.

There exists a randomized algorithm for unstable counting sort on
$n^{\prime}=\Omega(\log^{2}n)$ records with
$r^{\prime}=O(n^{\prime}\log^{O(1)}n)$ buckets in $O(n^{\prime}+r^{\prime})$
work and $O(\log n)$ span _whp_ in $n$ (Rajasekaran and Reif, 1989). However,
this approach is complicated and is less frequently used in practice
(discussions also given in Appendix B).

Input: $A$: input array with keys in $[r]$. $d$: number of digits remained to
be sorted. Each digit refers to $\mathit{\gamma}$ bits in the key. On the
root-level $d=(\log r)/\mathit{\gamma}$. $\theta$: base case threshold.

Notes: We use the $d$-th digit as the $d$-th digit from low to high.

1

2if  _$d=0$_ then return _$A$_

3 // All bits are sorted if  _$|A|\leq\theta$_ then return _ComparisonSort
$(A)$_ // Base case

4

5 _Distributing:_ Apply counting sort on $A$ using the $d$-th digit of the key
as bucket id

6 _Recursing:_ For $i\in[2^{\mathit{\gamma}}]$, let $A_{i}$ be the chunk in
$A$ where the records have the $d$-th digit as $i$

7 parallel_for  _$i\leftarrow 0$ to $2^{\mathit{\gamma}}-1$_ do Parallel-MSD-
Sort($A_{i}$, $d-1$)

8

Algorithm 1 Parallel-MSD-Sort($A[0..n-1],d$)

### 2.5 Comparison Sort, Semisort, and Sampling

Efficient and parallel sorting is widely studied (Inoue and Taura, 2015; Shen
and Ding, 2004; Bhattacharya et al., 2022; Träff, 2018; Katajainen and Träff,
1997; Sato et al., 2016; Shamoto et al., 2016). Comparison sort and semisort
are the other two types of sorting problems. Comparison sort takes an input of
$n$ records and a comparison function “$<$” to order two given records, and
outputs the records in non-decreasing order. Semisort takes an input of $n$
records and an equality-test “$=$” on two records, and reorders the records
such that records with equal keys are adjacent in the output (Dong et al.,
2023b). SOTA parallel algorithms for comparison sort (Axtmann et al., 2022;
Blelloch et al., 2020a) and semisort (Gu et al., 2015b; Dong et al., 2023b)
all use sampling that can detect heavy duplicate keys.

The sampling scheme used in these papers was first proposed in (Rajasekaran
and Reif, 1989). Given a parameter $p$ and input size $n$, the sampling scheme
first selects $p\log n$ samples uniformly at random. It then sorts the
samples, subsampling every $(\log n)$-th element into a list $S^{\prime}$.
Based on Chernoff bounds, any key that appears more than once in $S^{\prime}$
must have $\Omega(n/p)$ occurrences in the input _whp_. In existing
sorting/semisorting algorithms (Axtmann et al., 2022; Blelloch et al., 2020a;
Gu et al., 2015b; Dong et al., 2023b), each such key will be placed in its own
bucket. Since all records in this bucket have the same key, no further sorting
is needed. As such, these algorithms can be significantly faster when the
input contains many duplicate keys. Our DTSort is inspired by this technique
for detecting heavy duplicate keys. In Sec. 3, we will discuss how to apply
this sampling scheme to an integer sort algorithm.

Figure 2. An overview of the approach in the DTSort. Here
$r=16,\mathit{\gamma}=2$. For simplicity and space limit, the sampling scheme
in the figure is not exactly accurate as described in the algorithm. Here we
simply set keys with 2 or more samples as heavy keys.

## 3 The DovetailSort Algorithm

We now present DovetailSort (DTSort), which follows the MSD sort framework in
Alg. 1. By a careful design, DTSort can efficiently handle heavy duplicate
keys, and is _stable_.

To detect and handle duplicate keys, we use the sampling idea in existing
comparison-based sorts. However, we need to carefully combine this approach
with the MSD sort to avoid blowing up the work bound to $O(n\log n)$. As
mentioned, we want our algorithm to separate heavy and light keys, with the
goal of avoiding sorting heavy keys in all recursive levels. The idea of
identifying heavy/light buckets has been used in parallel semisort (Gu et al.,
2015b; Dong et al., 2023b). However, integrating it into a parallel MSD sort
is highly non-trivial. Unlike in semisort where heavy and light keys are
completely independent and can be placed in any order in the output, all
buckets in an integer sort must be finally placed in an increasing order. We
refer to this step as the “ _dovetail merging_ ” step, and discuss an
efficient solution both in theory and in practice in Sec. 3.4.

Input: $A$: input array with keys in $[r]$. $d$: number of digits remained to
be sorted. Each digit refers to $\mathit{\gamma}=\sqrt{\log r}$ bits in the
key. On the root-level $d=(\log r)/\mathit{\gamma}=\sqrt{\log r}$.

Notes: We use “the $d$-th digit” as the $d$-th digit from low to high.

1

2if  _$d=0$_ then return _$A$_

3 // All bits are sorted if  _$|A|\leq 2^{2\mathit{\gamma}}$_ then return
_ComparisonSort $(A)$ _

4 // Base cases

Sampling (find heavy keys & assign buckets):

5 $S\leftarrow\Theta(2^{\gamma}\log n)$ sampled keys from $A$

6 Sort $S$, subsample every ($\log n$)-th key, and store the keys that have
more than one subsamples into $S^{\prime}$

7 for  _$i\leftarrow 0$ to $|S^{\prime}|-1$_ do // each heavy keys’ MSD and
its key

8 $h[i]\leftarrow\langle$the $d$-th digit in $S^{\prime}[i]$,
$S^{\prime}[i]\rangle$

9for  _$i\leftarrow 0$ to $2^{\mathit{\gamma}}-1$_ do // each light bucket’s
MSD and a dummy key

10 $l[i]\leftarrow\langle i,-1\rangle$

11Merge $h$ and $l$ into an array $\mathit{buckets}$

12 Initialize the hash table $H$ and lookup table $L$

13 for  _$i\leftarrow 0$ to $|\mathit{buckets}|$_ do // assign bucket ids in
order

14 $\langle x,y\rangle\leftarrow\mathit{buckets}[i]$

15 if  _$y=-1$_ then $L[x]\leftarrow i$ // light bucket and its id else Insert
key $y$ with value $i$ to $H$ // heavy bucket and its id

16

17

Distributing (Reorder $A$ by bucket ids):

18

19Use counting sort on $A$ with key function $\mathsf{GetBucketID}$ (Alg. 2)

20

21

22

Recursing (sort light buckets):

23 parallel_for_each  _light bucket $B$_ do DTSort($B,d-1$)

Dovetail Merging (interleave light and heavy buckets):

24 parallel_for  _$i\in[2^{\mathit{\gamma}}]$_ do // merge buckets for all MSD
zones

25 Let $B_{0},B_{1},\dots B_{m}$ be all buckets in MSD zone $i$

26 DTMerge($B_{0},B_{1},\dots B_{m}$)

27return _$A$_

28

29 Function  _GetBucketId $(k)$_

30 if  _$k$ is found in $H$_ then return _$H[k]$_ else return _$L[k]$_

Algorithm 2 DTSort($A,d$)

We present DTSort in Alg. 2 and an illustration in Fig. 2. It consists of four
steps, explained in the subsections below. It is an MSD sort using
$\mathit{\gamma}=\sqrt{\log r}$ bits in each digit. We explain how this
parameter enables a good work bound in Sec. 4. We use $n$ as the original
input size, and $n^{\prime}$ as the current recursive subproblem size.
Separating them is needed to rigorously analyze the algorithm and give the
correct high probability bounds. We use _heavy/light bucket_ to mean a bucket
for heavy/light keys, and _heavy/light record_ to mean a record with a
heavy/light key. We use “MSD” to refer to the MSD in the current subproblem
(see the $d$-th digit in Alg. 2).

### 3.1 Step 1: Sampling

As mentioned earlier, the first step is to use sampling to determine heavy
keys, and assign a bucket id to each heavy and light bucket. The theory of
this idea is developed in previous work on samplesort and semisort
(Rajasekaran and Reif, 1989; Gu et al., 2015b).

In this step, $\Theta(2^{\gamma}\log n)$ samples are picked into an array $S$.
Recall that $\gamma$ is the number of bits to be sorted in a level. We then
subsample every $(\log n)$-th key in $S$, and output keys with more than one
subsample to $S^{\prime}$. As discussed in Sec. 2.5, we define all keys in
$S^{\prime}$ as heavy keys, and all other keys are light keys. According to
the MSD (the highest $\mathit{\gamma}$ bits) of the keys, the key range is
divided into $b=2^{\mathit{\gamma}}$ subranges. We call each key range sharing
the same MSD an _MSD zone_. In a regular integer sort, each MSD zone forms a
bucket. In our case, all _light keys_ within the same MSD zone will share a
bucket, but each _heavy key_ has its own bucket.

Unlike samplesort and semisort that do not need to further handle heavy keys
once they are recognized and distributed, in integer sort, we need to finally
combine them with light keys in the same MSD zone. Hence, we wish to keep the
heavy keys close to their final destination in the output from this very first
step, in order to perform an efficient merge at the end of the algorithm. To
do so, we will set the bucket id of a heavy bucket immediately after the light
bucket which it should finally be integrated into. Namely, for each MSD zone,
there is exactly one light bucket followed by some (possible zero) heavy
buckets, which are all the heavy keys with the same MSD as this zone. An
illustration is given in Fig. 2. A unique id is assigned to each bucket in a
serial order. Multiple heavy buckets in the same MSD zone will be totally
ordered by their full keys. In this way, the light buckets and heavy buckets
in the same MSD zone are consecutive, which makes it easier to combine them in
Step 4.

We explain this step in detail on Alg. 2–2 in Alg. 2. To assign the bucket
ids, we first obtain an array $h$ of all heavy keys in a sorted order (Alg.
2), which is a pair $(z,k)$, where $k$ is the key itself and $z$ is the MSD
zone of $k$. We then use an array $l$ (Alg. 2) to maintain a pair $(z,k)$ for
each light bucket (similar to $h$), where $z\in[2^{\mathit{\gamma}}]$ is the
MSD of the light bucket, and $k=-1$ is a dummy key to mark this is a light
bucket. By merging $h$ and $l$, all buckets will be re-organized in the
desired order—all buckets will first be ordered by their MSD; within the same
MSD, the first bucket is always the light one, followed by all heavy buckets
sorted by their key. We can then assign ids to them serially.

After assigning the bucket ids, we build a hash map $H$ to maintain the bucket
id of the heavy keys, and use a lookup array $L$ to map every MSD to its
associated light bucket id. As mentioned, the total number of buckets is
small, so the cost to maintain them is small both in time and space. With $H$
and $L$, we can determine the bucket id for any key $k$ in constant time (the
GetBucketID function in Alg. 2): we first look up whether $k$ is in $H$; if
so, $k$ is a heavy key, and its bucket id can be found in $H$. Otherwise, $k$
is a light key, and we can look up its bucket id in $L$ by using the MSD of
$k$.

### 3.2 Step 2: Distributing

After the first step, we know the bucket id of all keys. In this step, we use
the stable counting sort in Sec. 2.4 with their bucket id as the key, similar
to many existing parallel sorting algorithms (Dong et al., 2023b; Blelloch et
al., 2020a), After this step, all records are stably reordered into the
corresponding buckets.

### 3.3 Step 3: Recursing

After Step 2, keys are sorted across MSD zones, but keys in the same MSD zone
are not necessarily sorted. In this step, we first sort each light bucket
recursively.

As with a regular integer sort, the recursive call deals with the subproblem
using the same integer sort routine on the next “digit”. In DTSort, since all
records in a heavy bucket have the same key, no further sorting is needed.
Hence, we only recursively sort the light buckets (Alg. 2).

Note that by using a recursive approach, we can potentially detect much more
than $2^{\gamma}$ heavy keys if we consider all levels of recursion. Namely,
the heaviest keys will be identified at the topmost level and bypass all
recursive calls, and some medium-level heavy keys will be identified in the
middle levels but still bypass a few recursive calls. This enables a balance
in finding more heavy keys _in total_ to save work, as well as to use a
limited number of heavy keys _in each subproblem_ to avoid a large amount of
extra space in the counting sort in Step 2 (see a discussion in (Dong et al.,
2023b)).

Figure 3. Illustration of the dovetail merging step. The example merges the
buckets in MSD zone 01 in Fig. 2. We use a letter as subscription to
distinguish different records with the same key.

### 3.4 Step 4: Dovetail Merging

After Step 2, all keys are sorted across MSD zones, and all heavy keys are
ordered by their keys. Step 3 further sorts all light keys within each MSD
zone. To compute the final output, we finally must merge the two sorted
sequences (light records and heavy records) in each MSD zone. We refer to this
step as the dovetail merging step, which is the DTMerge function on Alg. 2 in
Alg. 2.

Each DTMerge call deals with an MSD zone with $m+1$ buckets $B_{0..m}$ stored
consecutively in array $A$. The first one $B_{0}$ is the light bucket in this
zone. The rest of the $m$ buckets are heavy buckets. Given these special
properties, the problem is equivalent to viewing the $m$ heavy buckets as
partitioning the light bucket into $m+1$ blocks, and then to interleave them.
Hence, we call this process dovetail merging. A naïve solution is to call a
standard parallel merge to combine $B_{0}$ with the rest $B_{1..m}$, with
linear work and polylogarithmic span. Despite the efficiency in theory, in
practice, it requires two rounds of out-of-place global data movements— _all_
records have to be merged into a temporary array, and then be copied back to
the input array. Such data movements will cause more cache misses and may make
the algorithm inefficient. There exists an in-place parallel merge algorithm
in theory (Gu et al., 2021), but it incurs a large constant hidden in the
asymptotic cost (moving each record at most four times), which may also be
expensive in practice. To enable good performance, our goal is to reduce the
data movements, with the additional consideration of moving less data out of
the input array.

Our idea is based on two observations. First, being strictly in-place is
unnecessary, since $O(n)$ temporary space is required by the counting sort in
Step 2. Therefore, we can reuse the space while minimizing the number of out-
of-place memory access. The second useful observation is that the number of
heavy keys is guaranteed to be small (at most $2^{\mathit{\gamma}}$ of them).
As a result, to achieve a practical solution, our idea is to process buckets
in the same MSD zone one by one _sequentially_ (note that all MSD zones are
still processed in parallel), and aim to achieve _high parallelism within each
bucket_.

With these insights, we present our DTMerge algorithm with an illustration in
Fig. 3 and pseudocode in Alg. 3. The high-level idea is to only copy the
smaller one of the light bucket ($B_{0}$) or the heavy buckets ($B_{1..m}$)
out to a temporary array, so that we copy at most half the records. The larger
half will be handled within the input array. WLOG, assume the current
subproblem has more heavy records than light ones (Alg. 3), and the other case
can be handled symmetrically (Alg. 3). In this case, we first move the light
records to a temporary array $T$ (Alg. 3), and then move all heavy buckets to
their final positions. To obtain the final positions of each heavy bucket, we
first binary search all heavy keys (in parallel) in the light bucket $B_{0}$
(Alg. 3). Combining this information with the sizes of all buckets, we can
infer the starting point of each heavy key. The next step is to (stably) move
the entire bucket $B_{i}$ to this starting point, and do so for all buckets
one by one from $B_{1}$.

Ideally, we wish to copy all records in $B_{i}$ to their final positions in
parallel, but have to be careful to ensure that their final positions are
“safe” to be overwritten. Since all heavy buckets are initially at the end of
this MSD zone, their final destinations must be ahead of their initial
positions. There are three cases for a bucket $B_{i}$’s destination. The first
case is when it contains light records, which is safe since all light records
have been backed up in $T$. The second is when it overlaps with an earlier
heavy bucket $B_{i^{\prime}}$ where $i^{\prime}<i$. It is also safe—since
$i^{\prime}<i$, bucket $i^{\prime}$ must have been moved to its final
destination. In the first two cases, we can safely move the records directly
into the final destination in parallel (Alg. 3). The last and the most
interesting case is when the destination of the bucket overlaps with its own
initial position (Alg. 3). Our idea is inspired by the in-place _circular
shift_ algorithm (Gries and Mills, 1981; Yang et al., 2013) (see Phase 3 in
Fig. 3): we can first flip (reverse the order) the bucket itself (Alg. 3), and
then flip the entire associated region in the array (Alg. 3). Flipping an
array can be performed in parallel and in-place by swapping the corresponding
records. As such, the entire chunk of $B_{i}$ is moved to an earlier position
in a _stable_ way. After dealing with all heavy buckets, we copy all records
in $T$ back to $A$, all in parallel, directly to their final destination (Alg.
3). In this way, our idea incurs out-of-place copies for at most half the
records, and requires moving the remaining records at most twice, but entirely
within $A$. Within each bucket, all operations (direct copy or swapping
records) are highly parallel.

Note that although buckets within the same MSD zone are processed
sequentially, this only increases the span of the entire algorithm by a factor
of $\log n$ because there are only $O(2^{\mathit{\gamma}})$ heavy buckets.
Interestingly, although a parallel merge gives better span bound than our
DTMerge, even in theory it is not necessary, because the counting sort in Step
2 also requires $O(2^{\mathit{\gamma}})$ span. We present more theoretical
analysis in Thm. 4.5.

In Sec. 6.3, we experimentally study DTMerge and show it accelerates Step 4
compared to using a parallel merge. We believe the “dovetail” step is a unique
challenge arising from detecting heavy keys in integer sorts, and our new
technique solves it in an effective yet simple way.

### 3.5 Base Cases

For the recursive structure, we use two base case conditions. The first is as
in regular integer sort when all digits have been sorted (Alg. 2 Alg. 2). The
second is when the problem size is smaller than $2^{2\mathit{\gamma}}$, where
we apply a comparison sort. Using comparison sort as a base case is standard
in implementations for parallel integer sort. Interestingly, instead of
viewing it as an empirical optimization, our theoretical analysis indicates
that this base case, as well as the threshold to use this base case, are very
important to guarantee theoretical efficiency. We present more details in Sec.
4.

Input: $m+1$ buckets consecutive in array $A$. The first one $B_{0}$ is a
sorted light bucket. The other $m$ buckets are heavy buckets sorted by their
keys.

1 Binary search all heavy keys in $B_{0}$, accordingly get $p_{1..m}$, where
$p_{m}$ is the starting index for heavy bucket $B_{i}$ when the array is fully
sorted.

2 if  _$|B_{0}|\leq\sum_{i=1}^{m}|B_{i}|$ _ then // more heavy keys than light

3 copy $B_{0}$ to array $T$

4 for  _$i\leftarrow 1$ to $m$_ do // starting from the first heavy bucket

// Move $B_{i}$ to the final positions starting from $A[p_{i}]$

5 if  _The final positions for $B_{i}$ overlap with the current positions of
$B_{i}$_ then

6 Let $s$ be the current start point of $B_{i}$

7 flip($A,s,s+|B_{i}|$) // Flip the entire region

8 flip($A,p_{i},s+|B_{i}|$) // Flip the destination region

9

10 else // No overlap. Directly copy

11 ParallelForEach  _$j\leftarrow 0$ to $|B_{i}|-1$_ do

12 Copy the $j$-th record in $B_{i}$ to $A[p_{i}+j]$

13 Move light records from $T$ to their final positions in $A$

14else Symmetric case: copy heavy keys out and merge

15

Algorithm 3 DTMerge($A,B_{0},B_{1},\dots B_{m}$)

## 4 The Theory of Parallel MSD Sort

In this section, we provide an analysis for _a class of MSD integer sorting
algorithms_ based on the framework in Alg. 1, including DTSort. Recall that
the algorithm sorts a “digit” of $\mathit{\gamma}$ bits in each round by 1)
distributing all records to the corresponding buckets by the current digit
being sorted, and 2) sorting each bucket recursively. Finally, a base case of
using comparison sort is applied to subproblems with sizes smaller than a
parameter $\theta$. As discussed in Sec. 1 and 2, although this parallel MSD
sort is widely used in practice, we are unaware of any tight analysis of it.
To the best of our knowledge, the best work bound is $O(n\log r)$ (Obeya et
al., 2019), which is no better than a comparison sort with $O(n\log n)$ work
for the realistic setting when $r=\Omega(n)$. This is inconsistent with the
excellent performance of existing parallel IS algorithms in practice, which
are much faster than SOTA comparison-based sorts for many inputs. In this
paper, we bridge this gap between theory and practice by showing Thm. 4.1, 4.4
and 4.5. Recall that the main use case for integer sorts is when $r=\Omega(n)$
(otherwise a counting sort has a better cost bound). We hence make this
assumption on $r$ in our analysis.

### 4.1 The Analysis of the General MSD Sort

As a warmup, and to illustrate some of the key theoretical ideas, we start by
showing that we can achieve an _unstable_ MSD sort with low work and span.

###### Theorem 4.1.

There exists an unstable parallel MSD sorting algorithm with $O(n\sqrt{\log
r})$ work and $O(\log r+\sqrt{\log r}\log n)$ span _whp_.

Under the assumption of the standard RAM model that can handle integers with
$O(\log n)$ bits (i.e., $r=n^{O(1)}$), the theorem implies $O(n\sqrt{\log n})$
work and $O(\log^{1.5}n)$ span _whp_.

The key in our analysis is to properly set up the parameters in Alg. 1. Once
we do so, the analysis becomes relatively clean. The parameters to use are
digit width $\mathit{\gamma}=\sqrt{\log r}$ and base case threshold
$\theta=2^{c\mathit{\gamma}}$ for some constant $c\geq 1$. This indicates
$(\log r)/\mathit{\gamma}=O(\sqrt{\log r})=O(\mathit{\gamma})$ recursive
levels in the entire algorithm. We use the _unstable counting sort_
(Rajasekaran and Reif, 1989) (see Sec. 2.4) for the distribution step, and use
comparison sort with $O(n^{\prime}\log n^{\prime})$ work and $O(\log
n^{\prime})$ span (Goodrich and Jacob, 2023) to handle base cases with size
$n^{\prime}$. The two main parts in Alg. 1 are the base cases (Alg. 1) and the
distribution step (Alg. 1). We first analyze their work with the parameters
shown above.

###### Lemma 4.2.

With the chosen parameters, the total work for base cases (Alg. 1) in Alg. 1
is $O(n\sqrt{\log r})$.

  * Proof.

Let $n_{1},\cdots,n_{k}$ be the base case sizes, where $\sum n_{i}=n$. Due to
the base case condition, $n_{i}\leq 2^{c\mathit{\gamma}}$, so the work to sort
it is $O(n_{i}\log n_{i})=O(n_{i}c\mathit{\gamma})$. The total work for all
base cases is therefore $O(nc\mathit{\gamma})=O(n\sqrt{\log r})$ ($c$ is a
constant). ∎

###### Lemma 4.3.

With the chosen parameters, the distribution step (Alg. 1) Alg. 1 in all
recursive calls has work $O(n\sqrt{\log r})$ _whp_.

  * Proof.

We start with the total work in each recursive level (i.e., those with the
same parameter $d$). Assume all subproblems in this level have sizes
$n_{1},\cdots,n_{k}$, where $\sum n_{i}\leq n$. Using the unstable counting
sort (Rajasekaran and Reif, 1989), the distribution work is $\sum
O(n_{i}+2^{\mathit{\gamma}})$ _whp_ for each subproblem. Due to the base case
condition, $n_{i}>2^{c\mathit{\gamma}}\geq 2^{\mathit{\gamma}}$ since $c\geq
1$. Hence the total work in each recursive level is $\sum O(n_{i})=O(n)$.
Consider all $O(\sqrt{\log r})$ recursive levels, the total work is
$O(n\sqrt{\log r})$. ∎

With the two lemmas, we can now prove Thm. 4.1.

  * Proof of Thm. 4.1.

Combining Lem. 4.2 and 4.3 gives the work bound in Thm. 4.1. We now prove the
span. The span for the base case is $O(\log
2^{c\mathit{\gamma}})=O(\mathit{\gamma})$ (Goodrich and Jacob, 2023). The
total span for one recursive call is the sum of $O(\log n)$ _whp_ (for the
distribution step) (Blelloch et al., 2020b) and $O(\gamma)$ (Alg. 1). Since
there are $O(\sqrt{\log r})$ recursive levels, the total span bound is as
stated. ∎

We note that Thm. 4.1 holds for $\mathit{\gamma}=\Theta(\sqrt{\log r})$, but
for simplicity we omit the $\Theta(\cdot)$ notation here.

_Practical and Stable MSD Sorting Algorithms._ Thm. 4.1 shows the theoretical
efficiency of parallel MSD algorithms. However, the counting sort (Rajasekaran
and Reif, 1989) used in the analysis is not stable, and is unlikely to be
practical due to too many rounds of data movements (Gu et al., 2015b). Many
practical algorithms choose to use the stable counting sort (Vishkin, 2010).
We analyze it as follows.

###### Theorem 4.4.

There exists a stable parallel MSD sorting algorithm with $O(n\sqrt{\log r})$
work and $O(2^{\sqrt{\log r}}\sqrt{\log r})$ span.

  * Proof.

First, note that the work analysis remains the same as Thm. 4.1, since the
stable counting sort has the same work as the unstable one. The only
difference is that the stable version has span related to the number of
buckets, which is $O(\log n+2^{\mathit{\gamma}})$. Assuming $r=\Omega(n)$,
this is $O(2^{\mathit{\gamma}})$. The total span of the entire MSD sorting is
therefore $O(2^{\sqrt{\log r}}\sqrt{\log r})$. ∎

This indicates the same work bound as the unstable version. Although the span
is slightly higher than that in Thm. 4.1, this version enables stability, and
is more practical. In fact, many existing practical parallel IS algorithms,
such as PLIS (Blelloch et al., 2020a), use a similar approach as in Thm. 4.4.

_Theory-Guided Practice._ Thm. 4.4 shows $O(n\sqrt{\log r})$ work for
practical parallel MSD sort. This explains why the integer sort algorithms can
outperform comparison sorts with $O(n\log n)$ work for realistic range of
$r=n^{O(1)}$. Our parameter $\gamma=\sqrt{\log r}$ used in the analysis is
also the choice for most practical MSD algorithms. When $r=2^{64}$, we have
$\gamma=\sqrt{\log r}=8$ and $\theta=2^{c\mathit{\gamma}}=2^{8c}$, which
roughly matches the existing algorithms that empirically sort 8–12 bits in
each recursive level, and use base case sizes from $2^{8}$ to $2^{16}$
(Axtmann et al., 2022; Blelloch et al., 2020a; Obeya et al., 2019).

### 4.2 The Analysis of DTSort

In this section, we show that our DTSort with the additional steps of sampling
and dovetail merging is theoretically efficient. We first prove the following
theorem.

###### Theorem 4.5.

The DTSort algorithm (Alg. 2) is a stable integer sort with $O(n\sqrt{\log
r})$ work and $\tilde{O}(2^{\sqrt{\log r}})$ span.

  * Proof.

DTSort can be viewed as an extension of the generic stable integer sort which
we proved work/depth bounds for in Thm. 4.4. Compared to the algorithm
analyzed in Thm. 4.4, DTSort has two additional steps: the sampling step and
the dovetail merging step. We first show that the asymptotic work bound
remains the same. Note that without sampling and merging, a subproblem with
size $n^{\prime}$ has $O(n^{\prime})$ work. In DTSort, the total number of
samples is $|S|=2^{\mathit{\gamma}}\log n$. We use base case threshold
$\theta=2^{2\mathit{\gamma}}$, and assume $r=\Omega(n)$, so the cost to sort
the samples is $O(|S|\log|S|)=O(2^{\mathit{\gamma}}\log
n\cdot\mathit{\gamma})=o(2^{2\mathit{\gamma}})=o(n^{\prime})$. Other parts in
the sampling step have work proportional to $|S|$, so the overall work in Step
1 is $o(n^{\prime})$.

In the dovetail merging step, we need $O(2^{\mathit{\gamma}})$ binary searches
on $n^{\prime}$ records, so the total cost is $O(2^{\mathit{\gamma}}\log
n^{\prime})$, which is $O(n^{\prime})$ considering $n^{\prime}\geq
2^{2\mathit{\gamma}}$ and $r=\Omega(n)$. When moving a bucket of size $t$, the
total work is $O(t)$ (either a direct copy or two rounds of flipping on at
most $2t$ records). We also copy out and back at most $n^{\prime}/2$ elements.
Combining both parts, the work for the merging step is also $O(n^{\prime})$.
Therefore, compared to the algorithm we used in Thm. 4.4, DTSort uses
asymptotically the same amount of work, $O(n\sqrt{\log r})$.

For the span, note that the algorithm used in Thm. 4.4 has
$O(2^{\mathit{\gamma}})$ span in each recursive call. Again consider the two
additional steps. Taking samples in parallel and sorting them requires
$O(\mathit{\gamma})$ span. We process $L$ and $H$ sequentially since there are
at most $2^{\mathit{\gamma}}$ heavy buckets and exactly $2^{\mathit{\gamma}}$
light buckets. This gives $O(2^{\mathit{\gamma}})$ span in each recursive
call.

For the dovetail merging step, we need to move at most $2^{\mathit{\gamma}}$
heavy buckets. Moving each of them is either a direct copy, or two rounds of
simple flipping. Both can be done by a parallel-for on all affected records,
giving $O(\log n)$ span assuming binary forking. This gives
$O(2^{\mathit{\gamma}}\log n)=\tilde{O}(2^{\mathit{\gamma}})$ span in one
recursive call. Therefore the two additional steps both have
$\tilde{O}(2^{\mathit{\gamma}})$ span in each recursive call. Considering all
$O(\sqrt{\log r})$ levels gives the stated span bound. ∎

We also note that for the dovetail merging step, if we use the simple parallel
merge (i.e., our baseline solution), we can get the same span bound of
$O(2^{\sqrt{\log r}}\sqrt{\log r})$ as in Thm. 4.4. Also, if we assume
arbitrary forking (e.g., PRAM), we can also obtain the same span bound as in
Thm. 4.4. To avoid the minor subtlety here, we hide the logarithmic terms in
the theorem by the $\tilde{O}(\cdot)$ notation. Our main claim is that DTSort
can achieve $O(n\sqrt{\log r})$ work and is highly parallel.

In addition to the worst-case bound, we can also show that DTSort can achieve
asymptotically better work bound on some real-world distributions with heavy
duplicates. We show two special cases that DTSort achieves linear work in Thm.
4.6 and 4.7. Here the exponential distribution means the frequencies of keys
following the probability density function $f(x;\lambda)=\lambda e^{-\lambda
x}$ for $x>0$ and a parameter $\lambda>0$.

###### Theorem 4.6.

DTSort has $O(n)$ work _whp_ if the input key frequency exhibits an
exponential distribution with $\lambda e^{\lambda}\geq\bar{c}/2^{\gamma}$ for
some $\bar{c}>1$. Here $\lambda>0$ is the parameter of the exponential
distribution, which gives probability density function $f(x;\lambda)=\lambda
e^{-\lambda x}$ for $x>0$.

  * Proof.

We first consider the exponential distribution, where the input follows the
probability density function $f(x;\lambda)=\lambda e^{-\lambda x}$ for $x>0$
and a parameter $\lambda>0$. Since we are considering integer sort, the
randomly drawn real number is rounded to an integer, and randomly mapped to an
integer in $[r]$. DTSort can detect a large fraction of the records as heavy
keys when the input is drawn from an exponential distribution. Particularly,
the sampling scheme in Alg. 2 can detect keys with more than
$\bar{c}n/2^{\gamma}$ occurrences as heavy keys _whp_ for some constant
$\bar{c}>1$ (Gu et al., 2015b), or equivalently, keys with frequencies higher
than $\bar{c}/2^{\gamma}$. This indicates that an integer $x$ with $\lambda
e^{-\lambda x}\geq\bar{c}/2^{\gamma}$ will be identified as heavy. The
cumulative distribution function of the exponential distribution is
$F(x;\lambda)=1-e^{-\lambda x}$ for $x\geq 0$. Thus, the total number of light
records is no more than $n\cdot e^{-\lambda
x}=\bar{c}n\gamma/2^{\gamma}=\bar{c}n\sqrt{\log r}/2^{\sqrt{\log r}}$ when
plugging in $x$ from the above threshold for heavy keys. However, since we are
considering integer sort (keys are rounded to integers), we need to make sure
that $x\geq 1$, so we need the assumption that $\lambda
e^{\lambda}\geq\bar{c}/2^{\gamma}$. Even though we assume the worst case that
the rest of the records will be processed for all the rest $\sqrt{\log r}-1$
levels, the cost to process them is $O(\bar{c}n\sqrt{\log r}/2^{\sqrt{\log
r}}\cdot\sqrt{\log r})=O(n\log r/2^{\sqrt{\log r}})=o(n)$. Hence, the total
work for DTSort on exponential distribution under the assumption is $O(n)$. ∎

###### Theorem 4.7.

DTSort has $O(n)$ work _whp_ if there are no more than $c^{\prime}2^{\gamma}$
distinct keys, for some constant $c^{\prime}<1$.

  * Proof.

For simplicity, we use the same notation $\bar{c}$ as the previous proof.
Similarly, keys with frequencies higher than $\bar{c}/2^{\gamma}$ will be
detected as heavy. The maximum number of light records is no more than
$c^{\prime}2^{\gamma}\cdot\bar{c}n^{\prime}/2^{\gamma}<c^{\prime}\bar{c}n^{\prime}$
for a subproblem of size $n^{\prime}$. By setting $c^{\prime}<1/\bar{c}$, the
number of light records will shrink by a constant fraction. The total cost in
a recursive level is linearly proportional to the total size of all
subproblems. Since the number of light records decreases geometrically across
levels, the total work is asymptotically bounded by the root level, which is
$O(n)$. ∎

## 5 Implementation Optimizations

We now discuss some optimizations we find useful when implementing DTSort.

_Overflow Bucket._ Note that integer sort takes the input keys from a range
from $0$ up to $r$, but the actual input keys may be in a smaller range. For
example, although the key can be 32-bit integers, the actual range can be much
smaller than $2^{32}$. To take advantage of this, one possibility is to
compute the maximal key by a parallel reduce, as in many existing
implementations (Kokot et al., 2017; Blelloch et al., 2020a). Here we discuss
a simple yet effective optimization to avoid this computation.

Similar to the heavy/light separation of the keys, here we also use sampling
to detect the key range. In each subproblem, DTSort sets the upper key range
$r^{\prime}$ to be the largest value in the sample set $S$. As such, DTSort
skips leading zeros in $r^{\prime}$ since creating buckets for these bits are
wasteful. Since the upper bound is computed by a sample set, DTSort creates an
_overflow bucket_ to hold records whose keys have at least a “1” for the
skipped bits. The size of the overflow bucket is usually very small. We put
the overflow bucket at the end, and comparison sort records in this bucket. As
such, we do not need to compute the maximal key in the input.

_Minimizing Data Movement._ In Alg. 2, for simplicity, in our description we
assume we use the stable distribution mentioned in Sec. 2.4 and write back the
bucket to $A$ (see Alg. 2). The second copy-back is not necessary and should
be avoided, since the performance of sorting algorithms can suffer due to
excessive data movement. In our implementation, we actually maintain two
arrays $A$ and $T$. After the distribution, we use array $T$ for future
execution, until we see another distribution or a dovetail merging, and we
then distribute or merge the results back to $A$ (or the other array in the
recursion). This approach will reduce the total data movement. Finally, if the
data remains in $T$ at the end of the algorithm, we need to move them back to
$A$ before finishing.

## 6 Experiments

Name | Stable | In-place | Type | Notes
---|---|---|---|---
DTSort | Yes | No | Integer | Our integer sort algorithm
PLIS | Yes | No | Integer |  ParlayLib integer sort (Blelloch et al., 2020a)
IPS2Ra | No | Yes | Integer |  IPS2Ra integer sort (Axtmann et al., 2022)
RS | No | Yes | Integer |  RegionsSort (Obeya et al., 2019)
RD | No | No | Integer |  RADULS (Kokot et al., 2018)
PLSS | Y/N | Y/N | Comparison |  ParlayLib sample sort (Blelloch et al., 2020a)
IPS4o | No | Yes | Comparison |  IPS4o sample sort (Axtmann et al., 2022)
Table 2. Algorithms tested in our experiments. There are two versions of PLSS. Here we use the unstable but faster one. In-place means fully in-place ($o(n)$ extra memory).  Instances | 32-bit key and 32-bit value pairs | 64-bit key and 64-bit value pairs
---|---|---
Integer | Comparison | Integer | Comparison
Ours | PLIS | IPS2Ra | RS | PLSS | IPS4o | Ours | PLIS | IPS2Ra | RS | RD | PLSS | IPS4o
Standard distributions:
Uniform | 109 | .500 | .537 | .671 | .718 | 1.27 | .690 | .994 | 1.14 | 1.09 | 1.43 | 1.86 | 1.65 | 1.11
107 | .501 | .549 | .600 | .705 | 1.14 | .604 | 1.03 | 1.15 | 1.06 | 1.71 | 1.86 | 1.47 | 1.05
105 | .478 | .542 | .595 | .696 | 1.08 | .653 | .859 | 1.22 | 1.01 | 1.36 | 2.66 | 1.35 | 1.10
103 | .506 | .505 | .538 | .613 | .805 | .432 | .795 | 1.41 | 1.66 | 1.54 | 3.23 | 1.01 | .759
10 | .308 | .707 | 1.13 | .438 | .959 | .456 | .581 | 1.93 | 3.78 | 1.26 | 8.25 | 1.12 | .850
Exponential | 1 | .526 | .536 | .574 | .711 | 1.11 | .671 | .976 | 1.16 | 1.04 | 1.39 | 2.28 | 1.33 | 1.16
2 | .502 | .546 | .577 | .711 | 1.12 | .661 | .919 | 1.22 | 1.05 | 1.39 | 2.45 | 1.35 | 1.19
5 | .435 | .567 | .583 | .705 | 1.11 | .612 | .819 | 1.52 | 1.03 | 1.41 | 2.44 | 1.35 | 1.03
7 | .419 | .582 | .554 | .708 | 1.08 | .609 | .782 | 1.69 | 1.01 | 1.48 | 2.53 | 1.32 | .972
10 | .402 | .603 | .560 | .682 | 1.09 | .561 | .763 | 1.87 | 1.05 | 1.53 | 2.61 | 1.28 | .930
Zipfian | 0.6 | .493 | .543 | .630 | .720 | 1.23 | .691 | 1.00 | 1.14 | 1.11 | 1.43 | 1.82 | 1.63 | 1.12
0.8 | .524 | .542 | .619 | .710 | 1.20 | .670 | 1.00 | 1.18 | 1.09 | 1.46 | 1.92 | 1.56 | 1.08
1 | .601 | .631 | .648 | .735 | 1.08 | .590 | 1.04 | 1.44 | 1.71 | 1.53 | 3.10 | 1.30 | 1.08
1.2 | .516 | .832 | 1.07 | .709 | 1.10 | .743 | .918 | 1.95 | 3.29 | 1.72 | 5.85 | 1.45 | 1.22
1.5 | .446 | .946 | 1.90 | .695 | 1.48 | .939 | .883 | 2.56 | 6.57 | 1.74 | 6.78 | 1.99 | 1.65
Avg. | .472 | .601 | .698 | .679 | 1.11 | .629 | .882 | 1.46 | 1.49 | 1.49 | 2.92 | 1.39 | 1.07
Adversarial distributions for integer sort with mixed heavy and light keys:
Bit-Exponential | 10 | 1.11 | .833 | 1.38 | .841 | .857 | .610 | 3.30 | 2.70 | 3.31 | 1.89 | 8.03 | 1.57 | 1.08
30 | .643 | 1.08 | 3.18 | .775 | 1.27 | .908 | 2.75 | 4.04 | 7.90 | 2.30 | 12.8 | 1.36 | 1.26
50 | .550 | 1.20 | 4.29 | .717 | 1.60 | 1.00 | 1.85 | 4.57 | 11.9 | 2.27 | 9.47 | 1.74 | 1.45
100 | .512 | 1.31 | 5.89 | .664 | 1.99 | 1.48 | 1.42 | 4.92 | 17.8 | 2.20 | 4.96 | 2.44 | 2.03
300 | .616 | 1.40 | 8.22 | .606 | 2.32 | 2.02 | 1.32 | 5.25 | 27.9 | 2.12 | 4.15 | 3.26 | 3.31
Avg. | .659 | 1.15 | 3.91 | .716 | 1.52 | 1.11 | 1.99 | 4.19 | 10.9 | 2.15 | 7.25 | 1.97 | 1.68
Table 3. Running time (in seconds) on synthetic data with
$\boldsymbol{n=10^{9}}$. The instances of one distribution is ordered by
increasing number of heavy records from top to bottom. The fastest running
time on each input instance is underlined. “Avg.” $=$ geometric mean. The keys
of RD need to be padded to multiples of 64 bits. We remove it from the 32-bit
experiments because it is too slow after padding to 64 bits.

We implemented DTSort in C++ and tested its performance. DTSort uses ParlayLib
(Blelloch et al., 2020a) to support fork-join parallelism.

_Experimental Setup._ We run our experiments on a 96-core machine (192 hyper-
threads) with 4 $\times$ 2.1 GHz Intel Xeon Gold 6252 CPUs (each with 36MB L3
cache) and 1.5TB of main memory. We always use numactl -i all to interleave
the memory on all CPUs except for sequential tests. We run each test six times
and report the median of the last five runs. All running times are in seconds.

_Parameter Selection._ DTSort mostly follows the theoretical analysis in Sec.
4, with the minor observation that using a variable radix width
$\mathit{\gamma}$ (based on $n$) slightly improved the overall performance.
Hence, DTSort uses $\mathit{\gamma}=\log_{2}\\!\sqrt[3]{n}$ and capped by 8
and 12 (guided by $\Theta(\\!\sqrt{\log r})$, leading to the same asymptotic
work bound), similar to some existing IS algorithms such as PLIS. The base
case threshold is $\theta=2^{14}$.

_Tested Algorithms._ We compare our DTSort with six _baseline_ parallel
algorithms, including four MSD integer sorting algorithms and two comparison-
based algorithms (both are samplesorts). We present their information in Tab.
2.

_Synthetic Distributions._ We used three _standard_ real-world distributions
to test all algorithms: _Unif- $\mu$_ (uniform distribution with $\mu$
distinct keys), _Exp- $\lambda$_ (exponential distribution with parameter
$10^{-5}\lambda$), and _Zipf- $s$_ (Zipfian distribution with parameter $s$),
We also designed a special distribution _Bit-Exponential_ (_BExp- $t$_), which
aims to simulate “worst-case” inputs for DTSort. For _BExp- $t$_, every bit of
a key is a 0 with probability $1/t$, and 1 otherwise. BExp controls the
distribution of the bitwise encoding (rather than frequency) of keys. With a
large $t$, most subproblems will have a mix of light and heavy records, which
is an adversarial case for DTSort due to the necessary dovetail merging for
almost every subproblem. In fact, BExp is generally hard for all MSD
algorithms because it makes subproblem sizes uneven. For simplicity, we say a
distribution is _heavy_ (or _heavier_ than another) if it contains more
duplicates, and _light_ (_lighter_) otherwise.

We choose five parameters for each distribution to give various ratios of
heavy keys (see Tab. 3). For uniform distributions, we test
$\mu=10,10^{3},10^{5},10^{7},10^{9}$. For exponential distributions, we test
$\lambda=1,2,5,7,10$. For Zipfian distributions, we test
$s=0.6,0.8,1,1.2,1.5$. For BExp distribution, we test $t=10,30,50,100,300$. We
round keys to the closest integer for exponential and Zipfian distributions.
For the standard distributions, we map the keys to larger ranges (up to
$2^{32}$ or $2^{64}$) to test the algorithms on the full range of $[r]$. We
use the geometric mean to compare the _average_ performance across multiple
tests.

### 6.1 Overall Performance on Synthetic Data

We present the running time of all tested algorithms on synthetic
distributions in Tab. 3 with different key lengths.

_Standard Distributions._ With 32-bit keys, DTSort has the best overall
performance, which is 1.27–1.47$\times$ faster than the baselines. A heatmap
for 32-bit results is also presented in Fig. 1. As expected, the other integer
sorting algorithms are usually competitive to DTSort on light distributions,
but can be much slower than DTSort on heavy distributions. Take the best
baseline on 32-bit keys, PLIS, as an example. On light distributions, It is
only 2–8% slower than DTSort, and is faster than comparison sorts. However, on
the heaviest several distributions in each category, PLIS is almost always
slower than the best comparison sort IPS4o, making its overall advantage to
IPS4o small (only 4%). On the other hand, DTSort can take advantage of the
heavy keys and outperforms IPS4o in all but two tests (which we study in more
details below). Across all tests, DTSort can be up to 2.3$\times$ faster than
PLIS.

With 64-bit keys, all algorithms perform slower than with 32-bit keys due to
the increase in data size, but the integer sorting algorithms are affected
more significantly. This matches the theory that the work bounds of these
algorithms increase with the value of $r$ (see Sec. 4). In this case, _all_
baseline integer sorting algorithms are slower than the fastest comparison
sort IPS4o on average, and many of them cannot outperform IPS4o even on some
light distributions. On the contrary, DTSort is still faster than both tested
comparison sorts in almost all tests, and is 1.2–1.5$\times$ faster on
average.

Among all tests on the standard distributions, DTSort _performs the best on
all but three tests_ , where samplesort IPS4o is better. Two of them are
_Unif- $10^{3}$_ for both 32- and 64-bits. The reason should be that, based on
our parameter $\gamma$, the average frequency of _Unif- $10^{3}$_ is the
threshold to distinguish heavy and light keys on the root level (about
$10^{6}$). We tested this on purpose to show the bad case for DTSort, since
all “light keys” are reasonably heavy (right around the threshold), but they
have to go through another level of recursion and have to be processed in the
dovetail merging step. This also illustrates the inherent difficulty of
dealing with duplicates for parallel IS algorithms. The last case where DTSort
is slower than IPS4o is _Zipf-1_ , and DTSort is within 5% slower.

In general, DTSort effectively takes advantage of different levels of heavy
duplicates: it exhibits a similar trend as samplesort algorithms where the
running time decreases with more duplicates. Meanwhile, as an integer sort,
DTSort also has the benefit of being asymptotically more efficient than
comparison sorts with realistic input key range $r=n^{O(1)}$.

_The Adversarial Distribution (_BExp_)._ As mentioned above, we design the
_BExp_ distribution to study an adversarial case of DTSort and other MSD
algorithms. The distribution incurs very different sizes of MSD zones for
_all_ subproblems, which incurs load imbalance for all MSD algorithms. It also
generates a mix of light and heavy buckets in almost all subproblems, which
greatly increases the cost of the dovetail merging and offsets the benefit of
identifying heavy keys.

Almost all existing IS algorithms have poor performance on this distribution.
Interestingly, even the comparison-based sorting algorithms perform much worse
when the distribution is very skewed. Among the MSD algorithms, IPS2Ra and RD
are significantly affected and can be up to 10$\times$ slower than the others.
DTSort, PLIS, and RS are lightly affected on 32-bit keys, but in general can
still be within 2$\times$ slower compared to their performance on standard
distributions. A possible explanation is that all the three algorithms use
similar work-stealing scheduler (cilk for RS, ParlayLib scheduler for DTSort
and PLIS), which supports dynamic load balancing. This matches our scalability
tests (reported in in Appendix C), where IPS2Ra and RD do not scale beyond 24
threads on this distribution. RS seems to be particularly resistant to the
BExp distribution. One reason is that RS is _unstable_ , and records with the
same key do not need to be put in a specific order. This sacrifice on
stability has a side-effect to benefit from the highly skewed distributions.

Overall, DTSort still achieves the best performance on 32-bit keys, and is
still faster than comparison sorts in most cases. We believe this advantage
comes from the removal of heavy keys from recursive problems, and thus
alleviates the load-imbalance issue. In addition, the performance of DTSort is
also significantly improved by using our proposed DTMerge algorithm. As we
will show in Sec. 6.3, as expected, on BExp distributions, the merging step is
very expensive, which can take 50% of the entire running time if we simply use
parallel merge. Our new algorithm accelerates this step by 1.7–2.8$\times$,
making the overall running time much lower.

Instances | $\boldsymbol{n}$ | Integer | Comparison
---|---|---|---
Ours | PLIS | IPS2Ra | RS | PLSS | IPS4o
Graph transpose |
| LJ | 69.0M | .043 | .043 | s.g. | .065 | .080 | .159
| TW | 1.47B | .888 | .942 | 3.24 | 1.05 | 1.57 | .891
| CM | 1.61B | .782 | .945 | 1.41 | 1.07 | 1.84 | 1.12
| SD | 2.04B | 1.10 | 1.29 | 2.87 | 1.34 | 2.08 | 1.23
| CW | 42.6B | 28.5 | 37.0 | 32.5 | s.g. | 60.1 | 24.6
| Avg. | .985 | 1.13 | - | - | 1.96 | 1.37
Morton order |
Real World | GL | 24.9M | .026 | .028 | .259 | .024 | .028 | .171
CM | 321M | .184 | .178 | .343 | .209 | .327 | .338
OSM | 2.77B | 2.32 | 2.39 | 3.65 | s.g. | 2.73 | 1.53
Avg. | .223 | .227 | .687 | - | .293 | .445
Varden (Gan and Tao, 2017) | SS2d | 1B | .498 | .557 | .662 | .634 | 1.27 | .775
SS3d | 1B | .512 | .568 | .778 | .611 | 1.12 | .754
SS2d’ | 2B | .973 | 1.16 | 1.27 | 1.17 | 2.44 | 1.39
SS2d’ | 2B | .990 | 1.60 | 2.86 | 1.17 | 2.30 | 1.96
| Avg. | .704 | .875 | 1.17 | .854 | 1.68 | 1.12
Table 4. Running time (in seconds) on multiple applications. We use 32-bit
keys and 32-bit values. See Sec. 6.2 for the detailed explanation on the keys
and values. The fastest running time on each instance is underlined. “$n$” $=$
input sizes. “s.g.” $=$ segmentation fault. “Avg.” $=$ the geometric mean on
instances of the same application. “-” $=$ not applicable.

### 6.2 Applications

Integer sort is widely used in practice. We choose two representative
applications on graph and geometric processing to evaluate the performance of
the tested algorithms.

_Graph Transpose._ Given a directed graph $G=(V,E)$, the graph transpose
problem is to generate $G^{\intercal}=(V,E^{\intercal})$ where
$E^{\intercal}=\\{(v,u):(u,v)\in E\\}$. Here we assume the edges are stored
using the standard format of compressed sparse row. To compute the transposed
graph, we only need to stably integer sort all edges using the key as the
second vertex $v$; here the vertices with larger degrees are the “heavy keys”.
Graph transpose is widely used as a subroutine in many graph algorithms
(Blelloch et al., 2020c; Dhulipala et al., 2021; Ji et al., 2018; Slota et
al., 2014).

We tested five real-world directed graphs, including soc-LiveJournal (LJ)
(Backstrom et al., 2006), twitter (TW) (Kwak et al., 2010), Cosmo50 (CM) (Kwon
et al., 2010; Wang et al., 2021), sd_arc (SD) (Meusel et al., 2014), and
clueweb (CW) (Meusel et al., 2014). The largest graph has 42.6 billion
directed edges. We save each edge as a (32-bit, 32-bit) integer pair since all
vertex ids are within 32 bits. On social networks (LJ and TW) and web graphs
(SD and CW), the degree distributions are more skewed (i.e., more duplicate
keys). CM is a $k$-NN graph, and the degrees are more evenly-distributed. We
present our results in Tab. 4.

Compared to the best baseline (PLIS), DTSort is always faster by up to
1.3$\times$. The samplesort IPS4o performs the best on CW, but performs much
slower on small graphs. Overall, DTSort has the best performance across all
graphs.

_Morton Sort._ For a $d$-dimensional point, a z-value is calculated by
interleaving the binary representations of each coordinate. The process to
sort the z-values of a set of points is called Morton sort, which orders
multidimensional data in one dimension while preserving the locality of data
points. It is widely used in geometry, graphics, and machine learning
applications (Blelloch et al., 2016; Wang et al., 2020; Blelloch and Dobson,
2022; Gu et al., 2013, 2015a; Li et al., 2012).

We tested seven datasets for Morton sort and list the result in Tab. 4. Three
are real-world datasets, including GeoLife (GL) (Zheng et al., 2008; Wang et
al., 2021), Cosmo50 (CM) (Kwon et al., 2010; Wang et al., 2021) and
OpenStreetMap (OSM) (roa, 2010). The other four are synthetic datasets by the
widely-used generator Varden (Gan and Tao, 2017) that creates points with
varying densities. On the synthetic distribution with varying densities,
DTSort outperforms all other baselines by 1.2–3.4$\times$. On the three
datasets with real-world points, existing integer sorts (RS or PLIS) achieved
the best performance on GL and CM, while DTSort is competitive (3–17% slower).
On OSM, none of the integer sorts outperform the samplesort IPS4o.
Interestingly, DTSort still achieved the best average performance on the three
datasets (PLIS is also close). This may be due to the reason that DTSort
_combines the benefit for both comparison and integer sorts_ , and hence _does
not have a very bad worst-case performance_.

(a)

(c)

(e)

(b)

(d)

(f)

Figure 4.  (a) and (b): Analysis for the performance of heavy-key detection.
Numbers are running time (lower is better) with or without heavy-key
detection. (a) is for 32-bit keys and (b) is for 64-bit keys. (c) and (d):
Analysis for the performance of dovetail merging. Numbers are running time
(lower is better) using our dovetail merging algorithm or a baseline merging
algorithm. (c) is for 32-bit keys and (d) is for 64-bit keys. (e) and (f):
Scalability (higher is better) with varying number of threads and running time
(lower is better) with varying input sizes on 32-bit key and 32-bit value
pairs on one instance: _Zipf-0.8_. Full analysis is given in Appendix C.
Discussions are in Sec. 6.3.

### 6.3 Performance Study

_Heavy Keys Detection Evaluation._ To understand the performance gain by
detecting heavy keys and processing them separately, we compare our
implementation with (labeled as DTSort) and without (labeled as Plain) heavy
keys detection in Fig. 4 (a) and (b). We select eight representative
distributions, which are the lightest and heaviest cases for each distribution
respectively.

As discussed in Sec. 3, the overhead of DTSort compared to Plain includes
sampling to detect heavy keys, and merging heavy keys with light keys. On
32-bit inputs, when there are no or fewer heavy keys, DTSort is comparable
with Plain, with at most 12% performance loss on _BExp-10_. It indicates that
the overhead of sampling is negligible compared to other steps. With heavy
duplicate keys, since heavy keys are gathered in separate buckets and do not
need to be sorted in later recursions, DTSort outperforms Plain by up to
2.12$\times$ (on _BExp-300_). Overall, DTSort outperforms Plain by 25% on
average. On 64-bit input, the performance improvement by using heavy keys
detection becomes larger. DTSort consistently outperforms Plain, with a
1.50$\times$ speedup on average.

_Dovetail Merge Evaluation._ In this section, we present an in-depth study of
the dovetail merging step that implements the heavy-light separation idea. As
mentioned in Sec. 3.4, our baseline solution (parallel merge), is correct and
even asymptotically better, but it requires to move all records between the
input array and a temporary array. To achieve the best performance, we
utilized the special property of this problem (i.e., very few heavy buckets)
and proposed our DTMerge algorithm in Sec. 3.4. We select seven representative
distributions to evaluate the merge performance. For uniform distribution, we
only use _Unif- $10^{3}$_, since no merge is needed in the lightest and
heaviest cases. For other distributions, we select the lightest and heaviest
cases.

To illustrate the effectiveness of our DTMerge approach, we compare it against
a standard parallel merge in ParlayLib (Blelloch et al., 2020a) (referred to
as PLMerge) as the baseline. When reporting the running time, we also run the
algorithms _without_ the merging step in all recursive calls. This will not
give a correct sorting result, but indicates the running time of the other
steps in DTSort. We present the result for both 32 bits and 64 bits in Fig. 4
(c) and (d). Since they generally have similar trends, we discuss the
experiments based on 32 bits.

On lighter distributions (e.g., _Exp-1_ and _Zipf-0.6_), the overall
improvement is marginal given that the dovetail merging step is fast anyway.
However, on heavier (more skewed) distributions, this step can take a
significant portion of the overall time if we use PLMerge, which is 17% on
_Exp-10_ and 42% on _Zipf-1.5_. With our DTMerge, this step only takes 6% in
_Exp-10_ and 28% for _Zipf-1.5_ , improving the overall performance for both
by about 20%.

For BExp, since almost all subproblems have a mix of heavy and light keys, the
merging time is significant. Using PLMerge, the merging step takes 27–43%
running time. Our DTMerge improved the performance of the merging step by up
to 2.8$\times$, leading to 1.1–1.38$\times$ faster sorting time.

_Scalability Tests._ DTSort has good scalability to both the number of
processors and input sizes. Due to page limit, we provide results on one
distribution (_Zipf-0.8_) as an example in Fig. 4 (e) and (f), and give the
full discussions in Appendix C.

## 7 Conclusion

We study the theory and practice of parallel integer sort. In theory, we show
that a class of widely-used practical MSD integer sorting algorithms have
$O(n\sqrt{\log r})$ work and low span. The bounds are asymptotically lower
than comparison sort for a large key range $r=n^{o(\log n)}$, and provide a
solid theoretical foundation for explaining the high performance of parallel
integer sort. In practice, we propose DTSort to overcome a common challenge in
existing MSD algorithms: the difficulty to detect and take advantage of
duplicate keys. DTSort combines algorithmic insights from both integer and
comparison sorting algorithms, and uses a special dovetail merging step to
make it more efficient. DTSort achieves the same $O(n\sqrt{\log r})$ work
bound as other MSD integer sorting algorithms, and has the same advantage of
processing heavy duplicates efficiently as samplesort algorithms. In our
experiments, DTSort achieves competitive or better performance than state-of-
the-art parallel integer and comparison sorting algorithms on various
synthetic and real-world datasets. As a future work, we consider further
improving the performance of parallel sorting algorithm by accelerating the
distributing step, inspired by some existing work (Obeya et al., 2019; Axtmann
et al., 2022; Gu et al., 2021).

## References

  * (1)
  * roa (2010) 2010\.  OpenStreetMap © OpenStreetMap contributors.  https://www.openstreetmap.org/.
  * Adve and Hill (1990) Sarita V Adve and Mark D Hill. 1990.  Weak ordering—a new definition. In _ACM International Symposium on Computer Architecture (ISCA)_ , Vol. 18. 2–14.
  * Albers and Hagerup (1997) Susanne Albers and Torben Hagerup. 1997.  Improved parallel integer sorting without concurrent writing.  _Information and Computation_ 136, 1 (1997), 25–51.
  * Andersson et al. (1998) Arne Andersson, Torben Hagerup, Stefan Nilsson, and Rajeev Raman. 1998.  Sorting in Linear Time?  _J. Computer and System Sciences_ 57, 1 (1998), 74–93.
  * Arora et al. (2001) N. S. Arora, R. D. Blumofe, and C. G. Plaxton. 2001.  Thread Scheduling for Multiprogrammed Multiprocessors.  _Theory of Computing Systems (TOCS)_ 34, 2 (01 Apr 2001).
  * Axtmann et al. (2022) Michael Axtmann, Sascha Witt, Daniel Ferizovic, and Peter Sanders. 2022.  Engineering in-place (shared-memory) sorting algorithms.  _ACM Transactions on Parallel Computing (TOPC)_ 9, 1 (2022), 1–62.
  * Backstrom et al. (2006) Lars Backstrom, Dan Huttenlocher, Jon Kleinberg, and Xiangyang Lan. 2006.  Group formation in large social networks: membership, growth, and evolution. In _ACM International Conference on Knowledge Discovery and Data Mining (SIGKDD)_. 44–54.
  * Bhatt et al. (1991) Pramod Chandra P Bhatt, Krzysztof Diks, Torben Hagerup, VC Prasad, Tomasz Radzik, and Sanjeev Saxena. 1991.  Improved deterministic parallel integer sorting.  _Information and Computation_ 94, 1 (1991), 29–47.
  * Bhattacharya et al. (2022) Arghya Bhattacharya, Abiyaz Chowdhury, Helen Xu, Rathish Das, Rezaul A Chowdhury, Rob Johnson, Rishab Nithyanand, and Michael A Bender. 2022.  When Are Cache-Oblivious Algorithms Cache Adaptive? A Case Study of Matrix Multiplication and Sorting. In _European Symposium on Algorithms (ESA)_.
  * Blelloch et al. (2020a) Guy E. Blelloch, Daniel Anderson, and Laxman Dhulipala. 2020a.  ParlayLib — a toolkit for parallel algorithms on shared-memory multicore machines. In _ACM Symposium on Parallelism in Algorithms and Architectures (SPAA)_. 507–509.
  * Blelloch and Dobson (2022) Guy E Blelloch and Magdalen Dobson. 2022.  Parallel Nearest Neighbors in Low Dimensions with Batch Updates. In _Algorithm Engineering and Experiments (ALENEX)_. SIAM, 195–208.
  * Blelloch et al. (2016) Guy E. Blelloch, Daniel Ferizovic, and Yihan Sun. 2016.  Just Join for Parallel Ordered Sets. In _ACM Symposium on Parallelism in Algorithms and Architectures (SPAA)_.
  * Blelloch et al. (2012) Guy E. Blelloch, Jeremy T. Fineman, Phillip B. Gibbons, and Julian Shun. 2012.  Internally deterministic parallel algorithms can be fast. In _ACM Symposium on Principles and Practice of Parallel Programming (PPOPP)_. 181–192.
  * Blelloch et al. (2020b) Guy E. Blelloch, Jeremy T. Fineman, Yan Gu, and Yihan Sun. 2020b.  Optimal parallel algorithms in the binary-forking model. In _ACM Symposium on Parallelism in Algorithms and Architectures (SPAA)_. 89–102.
  * Blelloch et al. (2010) Guy E. Blelloch, Phillip B. Gibbons, and Harsha Vardhan Simhadri. 2010.  Low depth cache-oblivious algorithms. In _ACM Symposium on Parallelism in Algorithms and Architectures (SPAA)_.
  * Blelloch et al. (2020c) Guy E. Blelloch, Yan Gu, Julian Shun, and Yihan Sun. 2020c.  Parallelism in Randomized Incremental Algorithms.  _J. ACM_ 67, 5 (2020), 1–27.
  * Blelloch et al. (1991) Guy E. Blelloch, Charles E. Leiserson, Bruce M. Maggs, C. Greg Plaxton, Stephen J. Smith, and Marco Zagha. 1991.  A Comparison of Sorting Algorithms for the Connection Machine CM-2. In _ACM Symposium on Parallelism in Algorithms and Architectures (SPAA)_.
  * Blumofe and Leiserson (1998) Robert D. Blumofe and Charles E. Leiserson. 1998.  Space-Efficient Scheduling of Multithreaded Computations.  _SIAM J. on Computing_ 27, 1 (1998).
  * Cho et al. (2015) Minsik Cho, Daniel Brand, Rajesh Bordawekar, Ulrich Finkler, Vincent Kulandaisamy, and Ruchir Puri. 2015.  PARADIS: an efficient parallel algorithm for in-place radix sort.  _Proceedings of the VLDB Endowment (PVLDB)_ 8, 12 (2015), 1518–1529.
  * Cormen et al. (2009) Thomas H. Cormen, Charles E. Leiserson, Ronald L. Rivest, and Clifford Stein. 2009.  _Introduction to Algorithms (3rd edition)_.  MIT Press.
  * Dhulipala et al. (2021) Laxman Dhulipala, Guy E. Blelloch, and Julian Shun. 2021.  Theoretically efficient parallel graph algorithms can be fast and scalable.  _ACM Transactions on Parallel Computing (TOPC)_ 8, 1 (2021), 1–70.
  * Dong et al. (2023a) Xiaojun Dong, Laxman Dhulipala, Yan Gu, and Yihan Sun. 2023a.  DovetailSort: A Parallel Integer Sort Algorithm.  https://github.com/ucrparlay/DovetailSort.
  * Dong et al. (2023b) Xiaojun Dong, Yunshu Wu, Zhongqi Wang, Laxman Dhulipala, Yan Gu, and Yihan Sun. 2023b.  High-Performance and Flexible Parallel Algorithms for Semisort and Related Problems. In _ACM Symposium on Parallelism in Algorithms and Architectures (SPAA)_.
  * Gan and Tao (2017) Junhao Gan and Yufei Tao. 2017.  On the hardness and approximation of Euclidean DBSCAN.  _ACM Transactions on Database Systems (TODS)_ 42, 3 (2017), 1–45.
  * Gharachorloo et al. (1990) Kourosh Gharachorloo, Daniel Lenoski, James Laudon, Phillip Gibbons, Anoop Gupta, and John Hennessy. 1990.  Memory consistency and event ordering in scalable shared-memory multiprocessors. In _ACM International Symposium on Computer Architecture (ISCA)_ , Vol. 18. ACM New York, NY, USA, 15–26.
  * Goodrich and Jacob (2023) Michael T Goodrich and Riko Jacob. 2023.  Optimal Parallel Sorting with Comparison Errors. In _ACM Symposium on Parallelism in Algorithms and Architectures (SPAA)_. 355–365.
  * Gries and Mills (1981) David Gries and Harlan Mills. 1981.  _Swapping sections_.  Technical Report. Cornell University.
  * Gu et al. (2015a) Yan Gu, Yong He, and Guy E Blelloch. 2015a.  Ray specialized contraction on bounding volume hierarchies. In _Computer Graphics Forum_ , Vol. 34. 309–318.
  * Gu et al. (2013) Yan Gu, Yong He, Kayvon Fatahalian, and Guy Blelloch. 2013.  Efficient BVH construction via approximate agglomerative clustering. In _High-Performance Graphics (HPG)_.
  * Gu et al. (2022) Yan Gu, Zachary Napier, and Yihan Sun. 2022.  Analysis of Work-Stealing and Parallel Cache Complexity. In _SIAM Symposium on Algorithmic Principles of Computer Systems (APOCS)_. SIAM, 46–60.
  * Gu et al. (2021) Yan Gu, Omar Obeya, and Julian Shun. 2021.  Parallel In-Place Algorithms: Theory and Practice. In _SIAM Symposium on Algorithmic Principles of Computer Systems (APOCS)_. 114–128.
  * Gu et al. (2015b) Yan Gu, Julian Shun, Yihan Sun, and Guy E. Blelloch. 2015b.  A Top-Down Parallel Semisort. In _ACM Symposium on Parallelism in Algorithms and Architectures (SPAA)_. 24–34.
  * Hagerup (1991) Torben Hagerup. 1991.  Constant-time parallel integer sorting. In _ACM Symposium on Theory of Computing (STOC)_. 299–306.
  * Inoue and Taura (2015) Hiroshi Inoue and Kenjiro Taura. 2015.  SIMD-and cache-friendly algorithm for sorting an array of structures.  _Proceedings of the VLDB Endowment_ 8, 11 (2015), 1274–1285.
  * Ji et al. (2018) Yuede Ji, Hang Liu, and H Howie Huang. 2018.  ispan: Parallel identification of strongly connected components with spanning trees. In _International Conference for High Performance Computing, Networking, Storage, and Analysis (SC)_. IEEE, 731–742.
  * Katajainen and Träff (1997) Jyrki Katajainen and Jesper Larsson Träff. 1997.  A meticulous analysis of mergesort programs. In _Italian Conference on Algorithms and Complexity (ICAC)_. Springer, 217–228.
  * Kokot et al. (2017) Marek Kokot, Sebastian Deorowicz, and Agnieszka Debudaj-Grabysz. 2017.  Sorting data on ultra-large scale with RADULS: New incarnation of radix sort. In _Beyond Databases, Architectures and Structures (BDAS)_. Springer, 235–245.
  * Kokot et al. (2018) Marek Kokot, Sebastian Deorowicz, and Maciej Długosz. 2018.  Even faster sorting of (not only) integers. In _Man-Machine Interactions 5: 5th International Conference on Man-Machine Interactions, ICMMI 2017 Held at Kraków, Poland, October 3-6, 2017_. Springer, 481–491.
  * Kwak et al. (2010) Haewoon Kwak, Changhyun Lee, Hosung Park, and Sue Moon. 2010.  What is Twitter, a social network or a news media?. In _International World Wide Web Conference (WWW)_. 591–600.
  * Kwon et al. (2010) YongChul Kwon, Dylan Nunley, Jeffrey P Gardner, Magdalena Balazinska, Bill Howe, and Sarah Loebman. 2010.  Scalable clustering algorithm for N-body simulations in a shared-nothing cluster. In _International Conference on Scientific and Statistical Database Management_. Springer, 132–150.
  * Lee et al. (2002) Shin-Jae Lee, Minsoo Jeon, Dongseung Kim, and Andrew Sohn. 2002.  Partitioned parallel radix sort.  _J. Parallel Distrib. Comput._ 62, 4 (2002), 656–668.
  * Li et al. (2012) Shengren Li, Lance C Simons, Jagaseesh Bhaskar Pakaravoor, Fatemeh Abbasinejad, John D Owens, and Nina Amenta. 2012.  $k$ANN on the GPU with shifted sorting. In _High-Performance Graphics (HPG)_.
  * Matias and Vishkin (1991) Yossi Matias and Uzi Vishkin. 1991.  On parallel hashing and integer sorting.  _J. Algorithms_ 12, 4 (1991), 573–606.
  * Meusel et al. (2014) Robert Meusel, Oliver Lehmberg, Christian Bizer, and Sebastiano Vigna. 2014.  Web Data Commons — Hyperlink Graphs.  http://webdatacommons.org/hyperlinkgraph.
  * Obeya et al. (2019) Omar Obeya, Endrias Kahssay, Edward Fan, and Julian Shun. 2019.  Theoretically-Efficient and Practical Parallel In-Place Radix Sorting. In _ACM Symposium on Parallelism in Algorithms and Architectures (SPAA)_. 213–224.
  * Polychroniou and Ross (2014) Orestis Polychroniou and Kenneth A Ross. 2014.  A comprehensive study of main-memory partitioning and its application to large-scale comparison-and radix-sort. In _ACM SIGMOD International Conference on Management of Data (SIGMOD)_. 755–766.
  * Rajasekaran and Reif (1989) Sanguthevar Rajasekaran and John H. Reif. 1989.  Optimal and sublogarithmic time randomized parallel sorting algorithms.  _SIAM J. on Computing_ 18, 3 (1989), 594–607.
  * Sato et al. (2016) Hitoshi Sato, Ryo Mizote, Satoshi Matsuoka, and Hirotaka Ogawa. 2016.  I/O chunking and latency hiding approach for out-of-core sorting acceleration using GPU and flash NVM. In _IEEE International Conference on Big Data (Big Data)_. 398–403.
  * Shamoto et al. (2016) Hideyuki Shamoto, Koichi Shirahata, Aleksandr Drozd, Hitoshi Sato, and Satoshi Matsuoka. 2016.  GPU-accelerated large-scale distributed sorting coping with device memory capacity.  _IEEE Trans. on Big Data_ 2, 1 (2016), 57–69.
  * Shen and Ding (2004) Xipeng Shen and Chen Ding. 2004.  Adaptive data partition for sorting using probability distribution. In _International Conference on Parallel Processing (ICPP)_. IEEE, 250–257.
  * Slota et al. (2014) George M. Slota, Sivasankaran Rajamanickam, and Kamesh Madduri. 2014.  BFS and coloring-based parallel algorithms for strongly connected components and related problems. In _IEEE International Parallel and Distributed Processing Symposium (IPDPS)_. IEEE, 550–559.
  * Solomonik and Kale (2010) Edgar Solomonik and Laxmikant V Kale. 2010.  Highly scalable parallel sorting. In _IEEE International Parallel and Distributed Processing Symposium (IPDPS)_. IEEE, 1–12.
  * Sun et al. (2019) Yihan Sun, Guy E Blelloch, Wan Shen Lim, and Andrew Pavlo. 2019.  On supporting efficient snapshot isolation for hybrid workloads with multi-versioned indexes.  _Proceedings of the VLDB Endowment (PVLDB)_ 13, 2 (2019), 211–225.
  * Träff (2018) Jesper Larsson Träff. 2018.  Parallel quicksort without pairwise element exchange.  _arXiv preprint:1804.07494_ (2018).
  * Vishkin (2010) Uzi Vishkin. 2010.  Thinking in parallel: Some basic data-parallel algorithms and techniques.  (2010).
  * Wang et al. (2023) Letong Wang, Xiaojun Dong, Yan Gu, and Yihan Sun. 2023.  Parallel Strong Connectivity Based on Faster Reachability. In _ACM SIGMOD International Conference on Management of Data (SIGMOD)_.
  * Wang et al. (2020) Yiqiu Wang, Yan Gu, and Julian Shun. 2020.  Theoretically-Efficient and Practical Parallel DBSCAN. In _ACM SIGMOD International Conference on Management of Data (SIGMOD)_. 2555–2571.
  * Wang et al. (2021) Yiqiu Wang, Shangdi Yu, Laxman Dhulipala, Yan Gu, and Julian Shun. 2021.  GeoGraph: A Framework for Graph Processing on Geometric Data.  _ACM SIGOPS Operating Systems Review_ 55, 1 (2021), 38–46.
  * Wassenberg and Sanders (2011) Jan Wassenberg and Peter Sanders. 2011.  Engineering a multi-core radix sort. In _European Conference on Parallel Processing (Euro-Par)_. Springer, 160–169.
  * Yang et al. (2013) Qingxuan Yang, John Ellis, Khalegh Mamakani, and Frank Ruskey. 2013.  In-place permuting and perfect shuffling using involutions.  _Inform. Process. Lett._ 113, 10-11 (2013), 386–391.
  * Zhang and Wu (2012) Keliang Zhang and Baifeng Wu. 2012.  A novel parallel approach of radix sort with bucket partition preprocess. In _International Conference on High Performance Computing (HPCC)_. IEEE, 989–994.
  * Zheng et al. (2008) Yu Zheng, Like Liu, Longhao Wang, and Xing Xie. 2008.  Learning transportation mode from raw gps data for geographic applications on the web. In _International World Wide Web Conference (WWW)_. 247–256.

##  Appendix A Determinacy Race

A _determinacy race_ is when two logically parallel operations access the same
memory location and at least one of them is a write (Cormen et al., 2009). We
say an algorithm is _race-free_ if it has no determinacy race. Note that this
requirement is stronger than data race (Gharachorloo et al., 1990; Adve and
Hill, 1990). A race-free algorithm is (internally) deterministic (Blelloch et
al., 2012), and has many advantages, including ease of reasoning about the
code, verifying correctness, debugging, and analyzing the performance. For
DTSort, all operations in the algorithm are deterministic once we fix the
random seed, which is a useful feature in practice.

##  Appendix B More Details about Counting Sort

Counting sort is a commonly-used primitive in parallel sorting algorithms to
distribute records into (a small number of) buckets. It takes an input of
$n^{\prime}$ records with a function to map each record to an integer (the
_bucket id_) in $[r^{\prime}]$ ($r^{\prime}\leq n^{\prime}$). The goal is to
sort the records by the bucket ids. The classic stable counting sort (Vishkin,
2010), which is widely used in practice, has $O(n^{\prime})$ work and
$O(r^{\prime}+\log n^{\prime})$ span, and can be made I/O-efficient (Blelloch
et al., 2010). Here we do not use $n$ since counting sort is used in the
recursive structure of many sorting algorithms and thus the subproblem size
$n^{\prime}$ can be smaller than $n$. The algorithm partitions the input array
into $l=\Theta(n^{\prime}/r^{\prime})$ subarrays each with size
$\Theta(r^{\prime})$. For each subarray, we sequentially count the appearances
of each key, and in total we store $l\cdot r^{\prime}=\Theta(n)$ such counts
in a matrix (the counting matrix). The total count of each key can be computed
by a reduce on the counts for all subarrays. We then apply a prefix-sum on the
counting matrix to get the offset of each bucket in each subarray, and finally
distribute the records to their destinations in the buckets. In practice, to
make it faster, we pick the number of subarrays smaller than $m/r^{\prime}$ so
the counting matrix is smaller (and all operations on it can be faster). Under
this premise, a smaller number of subarrays is preferred to improve the
performance, since it makes the counting matrix small to fit in cache (Dong et
al., 2023b).

There exists a randomized algorithm for unstable counting sort on
$n^{\prime}=\Omega(\log^{2}n)$ records with
$r^{\prime}=O(n^{\prime}\log^{O(1)}n)$ buckets in $O(n^{\prime}+r^{\prime})$
work and $O(\log n)$ span _whp_ in $n$ (Rajasekaran and Reif, 1989). However,
this approach is complicated and is less frequently used in practice due to
some overhead. The reason is that this algorithm requires allocating an
$O(n)$-size array but with a large constant hidden in the Big-O, and uses many
random accesses to distribute the records to the associated buckets which is
I/O-unfriendly. Finally, it also needs to pack the large array. Putting all
pieces together makes the algorithm less practical than the stable counting
sort mentioned above.

##  Appendix C More Experimental Results

_Parallel Scalability._ We test all algorithms using various numbers of
threads (from 1 to 192) to show how they scale with the number of threads. For
each distribution, we choose two representative tests (a heavy one and a light
one), and present our results in Fig. 12, 12, 12, 12, 12, 12, 12, 12, 20, 20,
20, 20, 20, 20, 20 and 20. On the three standard distributions, most
algorithms present favorably good speedups. Our algorithm achieves the top-
tier (almost linear) speedup, with up to 72$\times$ on 192 hyper-threads.

On the BExp distribution, the self-relative speedups are slightly worst for
all algorithms. In particular, IPS2Ra, IPS4o and RD do not scale well (or
become even worse) with increasing number of threads. As discussed in the main
body, this is mainly because BExp creates subproblems with very different
sizes, and thus is challenging to be handled by integer sorts due to load
imbalance. In most cases, DTSort, RS, PLSS and PLIS still scale to at least 96
cores (RS’s performance drops slightly with hyperthreading). This may be due
to the use of a work-stealing scheduler in these implementations, which
enables dynamic load balancing. Our algorithm still achieves almost-linear
scalability, and has the best speedups in most cases.

_Input Size Scalability._ We test all algorithms with various input sizes for
both 32-bit and 64-bit keys. For each distribution, we choose two
representative tests (a heavy one and a light one), and present our results in
Fig. 28, 28, 28, 28, 28, 28, 28, 28, 36, 36, 36, 36, 36, 36, 36 and 36. We
remove RD from these figures since it has a serious scale-down issue on some
of the distributions (needs more than two seconds with $n=10^{7}$). We measure
the input size scalability starting from $n=10^{7}$, because on an even
smaller size, most algorithms only take less than ten milliseconds. Overall,
our algorithm DTSort, PLIS, and PLSS scale well in a large range ($n=10^{7}$
to $2\times 10^{9}$). and our algorithm achieves the best performance in most
cases. RS, IPS2Ra, and IPS4o do not perform well when the input sizes are
relatively small ($n\leq 10^{8}$), but they also achieve good performance when
the input sizes are sufficiently large.

Self-speedup with Varying Thread Counts ($n=10^{9}$, 32-bit keys and 32-bit
values). Higher is better.
Uniform Distribution
Figure 5. Self-speedup with varying thread counts of all tested
implementations on _Unif–_ $\boldsymbol{10^{7}}$. Figure 6. Self-speedup with
varying thread counts of all tested implementations on _Unif–_
$\boldsymbol{10^{3}}$.

Exponential Distribution
Figure 7. Self-speedup with varying thread counts of all tested
implementations on _Exp–_ $\boldsymbol{2}$. Figure 8. Self-speedup with
varying thread counts of all tested implementations on _Exp–_
$\boldsymbol{7}$.

Zipfian Distribution
Figure 9. Self-speedup with varying thread counts of all tested
implementations on _Zipf–_ $\boldsymbol{0.8}$. Figure 10. Self-speedup with
varying thread counts of all tested implementations on _Zipf–_
$\boldsymbol{1.2}$.

BExp Distribution
Figure 11. Self-speedup with varying thread counts of all tested
implementations on _BExp–_ $\boldsymbol{30}$. Figure 12. Self-speedup with
varying thread counts of all tested implementations on _BExp–_
$\boldsymbol{100}$.

Self-speedup with Varying Thread Counts ($n=10^{9}$, 64-bit keys and 64-bit
values). Higher is better.
Uniform Distribution
Figure 13. Self-speedup with varying thread counts of all tested
implementations on _Unif–_ $\boldsymbol{10^{7}}$. Figure 14. Self-speedup with
varying thread counts of all tested implementations on _Unif–_
$\boldsymbol{10^{3}}$.

Exponential Distribution
Figure 15. Self-speedup with varying thread counts of all tested
implementations on _Exp–_ $\boldsymbol{2}$. Figure 16. Self-speedup with
varying thread counts of all tested implementations on _Exp–_
$\boldsymbol{7}$.

Zipfian Distribution
Figure 17. Self-speedup with varying thread counts of all tested
implementations on _Zipf–_ $\boldsymbol{0.8}$. Figure 18. Self-speedup with
varying thread counts of all tested implementations on _Zipf–_
$\boldsymbol{1.2}$.

BExp Distribution
Figure 19. Self-speedup with varying thread counts of all tested
implementations on _BExp–_ $\boldsymbol{30}$. Figure 20. Self-speedup with
varying thread counts of all tested implementations on _BExp–_
$\boldsymbol{100}$.

Scalability with increasing input size (32-bit keys). Lower is better.

Uniform Distribution
Figure 21. Scalability with increasing input size ($n$) of all tested
implementations on _Unif–_ $\boldsymbol{10^{7}}$. Figure 22. Scalability with
increasing input size ($n$) of all tested implementations on _Unif–_
$\boldsymbol{10^{3}}$.

Exponential Distribution
Figure 23. Scalability with increasing input size ($n$) of all tested
implementations on _Exp–_ $\boldsymbol{2}$. Figure 24. Scalability with
increasing input size ($n$) of all tested implementations on _Exp–_
$\boldsymbol{7}$.

Zipfian Distribution
Figure 25. Scalability with increasing input size ($n$) of all tested
implementations on _Zipf–_ $\boldsymbol{0.8}$. Figure 26. Scalability with
increasing input size ($n$) of all tested implementations on _Zipf–_
$\boldsymbol{1.2}$.

BExp Distribution
Figure 27. Scalability with increasing input size ($n$) of all tested
implementations on _BExp–_ $\boldsymbol{30}$. Figure 28. Scalability with
increasing input size ($n$) of all tested implementations on _BExp–_
$\boldsymbol{100}$.

Scalability with increasing input size (64-bit keys). Lower is better.

Uniform Distribution
Figure 29. Scalability with increasing input size ($n$) of all tested
implementations on _Unif–_ $\boldsymbol{10^{7}}$. Figure 30. Scalability with
increasing input size ($n$) of all tested implementations on _Unif–_
$\boldsymbol{10^{3}}$.

Exponential Distribution
Figure 31. Scalability with increasing input size ($n$) of all tested
implementations on _Exp–_ $\boldsymbol{2}$. Figure 32. Scalability with
increasing input size ($n$) of all tested implementations on _Exp–_
$\boldsymbol{7}$.

Zipfian Distribution
Figure 33. Scalability with increasing input size ($n$) of all tested
implementations on _Zipf–_ $\boldsymbol{0.8}$. Figure 34. Scalability with
increasing input size ($n$) of all tested implementations on _Zipf–_
$\boldsymbol{1.2}$.

BExp Distribution
Figure 35. Scalability with increasing input size ($n$) of all tested
implementations on BExp-$\boldsymbol{30}$. Figure 36. Scalability with
increasing input size ($n$) of all tested implementations on
BExp-$\boldsymbol{100}$.