# A Laser-based Time Calibration System
for the MEG II Timing Counter

M. Nakao, G. Boca, P. W. Cattaneo, M. De Gerone, F. Gatti, M. Nishimura, W.
Ootani, M. Rossella,
Y. Uchiyama, M. Usami, and K. Yoshida Manuscript received December 10, 2018.
This work was supported by INFN, Italy and JSPS, Japan KAKENHI Grant Number
JP26000004 and JP17J04114.M. Nakao, M. Nishimura, and M. Usami are with the
Department of Physics, The University of Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo
113-0033, Japan (e-mail: nakao@icepp.s.u-tokyo.ac.jp)K. Yoshida was with the
Department of Physics, The University of Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo
113-0033, Japan (e-mail: nakao@icepp.s.u-tokyo.ac.jp)G. Boca, P. W. Cattaneo,
and M. Rossella are with INFN Pavia, 6-27100, Pavia, Italy, and G. Boca is
also with the Department of Physics, The University of Pavia, 6-27100, Pavia,
Italy.M. De Gerone and F. Gatti are with INFN Genova, 33-16146 Genova, Italy,
and also with the University of Genova, 33-16146 Genova, Italy.W. Ootani and
Y. Uchiyama are with ICEPP, The University of Tokyo, 7-3-1 Hongo, Bunkyo-ku,
Tokyo 113-0033, Japan

###### Abstract

We have developed a new laser-based time calibration system for the MEG II
timing counter dedicated to timing measurement of positrons. The detector
requires precise timing alignment between $\sim\,$500 scintillation counters.
In this study, we present the calibration system which can directly measure
the time offset of each counter relative to the laser-synchronized pulse. We
thoroughly tested all the optical components and the uncertainty of this
method is estimated to be 24 ps. In 2017, we installed the full system into
the MEG II environment and performed a commissioning run. This method shows
excellent stability and consistency with another method. The proposed system
provides a precise timing alignment for SiPM-based timing detectors. It also
has potential in areas such as TOF-PET.

##  I Introduction

Lepton flavor violating decays are powerful tools to search for new physics
beyond the standard model. The MEG II experiment at Paul Scherrer Institut
(PSI) in Switzerland will search for one of these decays, $\mu^{+}\to
e^{+}\gamma$, with a sensitivity of $6\times 10^{-14}$ improving the existing
limit of almost an order of magnitude[1]. The sensitivity of MEG, the first
phase of the experiment, was limited by the detector performance. The
$\mu^{+}$ beam intensity was reduced to satisfy the requirements of stable
operation of the detectors. Therefore, in order to make full use of the
available beam intensity of $7\times 10^{7}\mu^{+}/\mathrm{s}$, we developed
new detectors with resolutions in energy, position, and timing improved by a
factor of two. In particular, good time resolution plays an important role to
reduce the dominant background that depends quadratically on the beam
intensity.

The pixelated Timing Counter (pTC), one of the key detectors in MEG II, is
dedicated to the measurement of the positron timing. Its design is based on a
new approach to improve the time resolution: the multiple hit scheme[2]. The
pTC consists of 512 scintillation counters (Fig. 1) and is divided into two
equal sectors: one placed upstream the MEG II target and the other downstream.
Each counter, with a size of 40 (50) $\times$ 120 $\times$ 5 mm3, consists of
a fast plastic scintillator (SAINT-GOBAIN, BC422) read out by six series-
connected silicon photomultipliers (SiPM, AdvanSiD, ASD-NUV3S-P High-Gain, 3
$\times$ 3 mm2, 50 $\times$ 50 $\mu$ m 2, V${}_{\mathrm{breakdown}}\sim\,$24
V) at both short ends. Signal positrons hit $\sim\,$9 counters on overage and
the total time resolution is expected to improve with the number of hit
counters. This approach enables us to achieve a time resolution of 35 ps, two
times better than in MEG.

Exploiting the multiple hit scheme requires a precise relative time
calibration between counters, which is a challenging task. We have developed a
new time calibration method using a pulsed laser and commercially available
optical components. This method measures directly the time offsets for the
timing alignment.

Figure 1: A single scintillation counter with height 50 mm (or 40 mm depending
on its position in the pTC).

##  II Overview of the laser-based time calibration system

(a)

(b)

Figure 2: (a)Schematic and (b)picture of the laser-based time calibration
system for the MEG II timing counter. A detailed description can be found in
Sec.II.

The laser-based time calibration system is shown in Fig. 2. The laser pulse
(“Laser controller” and “Laser head”) is split into each scintillation counter
simultaneously via several optical components: a mode scrambler, optical
splitters, and optical fibers. The laser device is produced by Hamamatsu with
a part number of PLP-10 and has the following features: wavelength: 401 nm,
pulse duration: 60(typical)–100(max) ps, peak power: 200 mW. The mode
scrambler produced by OPTO SCIENCE, INC. stabilizes transmission mode in a
fiber. The optical switch produced by LEONI switches light paths. Optical
splitters divide the input light into several channels. 1 to 2 splitter is
produced by OZ OPTICS and the others are produced by Lightel Technologies Inc.
We use multimode and graded index fibers with a core diameter of 50 $\mu$m and
a cladding diameter of 125 $\mu$m. There are produced by OZ OPTICS (for 2.5 m)
and Lightel Technologies Inc. (for 10 m). The selection of these items was
based on a detailed study presented in [3]. To get enough laser power for each
counter, one-eighth of counters ($\sim\,64$) are illuminated at one data
taking. The time offset of each counter is measured relative to laser-
synchronized pulse. We have developed a new method to fix an optical fiber
into the scintillator: a small hole with a diameter of 2.5 mm and with a depth
of 1 mm is drilled into the bottom long side of the scintillator and a
polycarbonate screw (at the bottom in Fig. 1), which has a fiber-shaped hole
inside, fixes the optical fiber together with the black support bar made of
ABS resin.

Slow control tasks such as monitoring temperatures, controlling the laser, and
switching the optical switch are performed based on the Midas Slow Control
System[4] developed at PSI and TRIUMF.

##  III Research and development

We focused on the following items in research and development of the laser
calibration system.

###  III-A Fiber insertion method

Firstly, we confirmed that the hole in the scintillator does not degrade the
time resolution. Then we considered whether we should put optical grease
inside the hole expecting better optical contact. However, it turned out that
the time offset depends on the amount of the grease. The amount of grease
during the assembly could not be reproduced satisfactorily. Therefore, we
decided not to use optical grease. We checked the reproducibility: the fiber
was inserted and the time offset was measured repeatedly. The reproducibility
of the time offset is highly satisfying (3.8 ps in standard deviation). The
stability was also confirmed during the commissioning phase (Sect. IV).

###  III-B Temperature dependence

Time offset of all optical components depends significantly on temperature. We
measured the temperature coefficient of the whole system to be 7.6$\pm$1.5
ps/∘C. The main contribution was from the “SYNC signal” cable (RG174/U) in
Fig. 2a; we will replace it with another cable with a smaller temperature
coefficient (FSJ1-50A) before the engineering run in 2018. In addition, the
temperature in the detector area is expected to be stabilized within 1 ∘C,
thus the temperature effect on the time offset will be negligible.

###  III-C Mass test

All optical components including $\sim\,$600 fibers and $\sim\,$80 splitters,
were tested before installation. Particularly the optical length in the
optical splitter has a large difference among different production batches.
However, it does not affect the performance because this difference of optical
length up to counter are subtracted from the time offset.

###  III-D Uncertainty of the system

Finally, the uncertainty of the laser calibration is calculated to be 24 ps by
summing up all the contributions which can affect the time offset such as
transit time variation inside the SiPMs, reproducibility of fiber insertion
method, and statistical uncertainty in the mass test.

A detailed study of the mass test and an estimation of the uncertainty are
presented in [5][6].

##  IV Commissioning and performance

In 2016 one-fourth of the system was operated and we reported the result in
[7]. In 2017 we finished the construction and the assembly of the pTC and then
we installed it into the MEG II experimental area at PSI to operate it under
the MEG II $\mu^{+}$ beam (7$\times 10^{7}\mu/\mathrm{s}$).

Data with positrons from Michel decay ($\mu^{+}\to e^{+}\nu\nu$) were
collected for two weeks and the detector was operated for two months. The time
offset calculated with this method was stable enough over time (2.5 ps in
standard deviation for one month). As shown in Fig. 3, this method was
compared with another calibration method, which uses positron tracks, and they
are in good agreement within $\sim\,$50 ps.

Figure 3: The difference of time offsets between the laser-based method and
the track-based method. Each error bar includes the systematic and statistical
uncertainties of the laser-based method ($\sim\,$24 ps). The standard
deviation is 54 ps downstream (DS) and 46 ps upstream (US).

##  V Conclusion

We have developed the laser-based time calibration system for the MEG II pTC
with an uncertainty of 24 ps. The time offset calculated with this method is
stable enough for the detector operation. In addition, the consistency of the
time offset was checked with another method.

Laser-based methods for TOF detectors in the past were commisioned only for
plastic scintillator detectors with photomultiplier tube read-out; however,
SiPM read-out enables us to design more compact and flexible timing detectors
and the rapid development of SiPM technology suggests that scintillation
detector with SiPM read-out will be extensively used in a wide variety of
applications.

Therefore, laser-based time calibration methods are expected to become more
and more widespread. Furthermore, this method is precise enough to be
compatible with the multiple hit scheme. It also has potential in areas such
as Time-Of-Flight positron emission tomography (TOF-PET).

The full engineering run of the MEG II experiment will start in 2019 followed
by the physics run, searching for the lepton flavor violating muon decay.

## References

  * [1] A.M. Baldini, et al., “The design of the MEG II experiment”, Eur. Phys. J. C (2018) 78: 380, arXiv:1801.04688
  * [2] Y. Uchiyama, et al., “30-ps Time Resolution with Segmented Scintillation Counter for MEG II”, Nucl. Instrum. Methods A, 845:507–510, 2017.
  * [3] R. Bertoni, et al., “A laser diode based system for calibration of fast time-of-flight detectors”. J. Instrum., 11(05):P05024, 2016.
  * [4] MIDAS website: midas.triumf.ca
  * [5] M. Rossella, et al., “Design and test of the calibration system of the MEGII Pixelated timing Counter”, Nucl. Instrum. Methods A, 2018 (in press)
  * [6] M. Nakao, et al., “A laser–based time calibration system for highly segmented scintillation counters” (in preparation)
  * [7] M. Nakao, et al., “Results from Pilot Run for MEG II Positron Timing Counter”, Springer Proc. Phys. 213 (2018) 237-241, arXiv:1808.07279


# ThielSort: Implementing the Diverting Fast Radix Algorithm

Stuart Thiel<EMAIL_ADDRESS>0000-0002-2172-3916 Concordia
University1455 De Maisonneuve Blvd. W.MontrealQuebecCanadaH3G 1M8 ,  Larry
Thiel<EMAIL_ADDRESS>Concordia University1455 De Maisonneuve Blvd.
W.MontrealQuebecCanadaH3G 1M8 and  Gregory Butler<EMAIL_ADDRESS>0000-0002-6938-0879 Concordia University1455 De Maisonneuve Blvd.
W.MontrealQuebecCanadaH3G 1M8

(2022)

###### Abstract.

This paper presents ThielSort, a practical implementation of the Diverting
Fast Radix (DFR) Algorithm. The theoretical improvements over classical radix
sorts are outlined and implementation details are specified to demonstrate
that the algorithm is competitive with the state of the art. The effectiveness
of this implementation of the DFR algorithm is shown by considering a variety
of standard distributions of data and input sizes.

sorting, radix sort, diverting

††copyright: acmcopyright††journalyear: 2022††doi:
XXXXXXX.XXXXXXX††conference: ACM Journal of Experimental Algorithmics; Feb.
20, 2022; Woodstock, NY††price: 15.00††isbn: 978-1-4503-XXXX-X/18/06††ccs:
Theory of computation Sorting and searching††ccs: Theory of computation Divide
and conquer

##  1\. Introduction

Exactly how Radix Sort came to be is unclear. Knuth cite’s Hollerith’s sorting
and tabulating machines that were used in the US Census Office of the late
1800s and into the early 1900s(Knuth, 1998, p.384) as a solid early example of
sorting, but Cormen et al. summarizes it well with the following(Cormen et
al., 2009, p.211):

> Knuth credits H.H. Seward with inventing counting sort in 1954, as well as
> with the idea of combining counting sort with radix sort. Radix sorting
> starting with the least significant digit appears to be a folk algorithm
> widely used by operators of mechanical card-sorting machines.

Its origins as a folk algorithm, often referenced to the 1880s in terms of
mechanical use, the fact that there are both most and least significant
flavors of the algorithm and even changes in use of language over time lead to
some confusion. In (Friend, 1956), care is taken to distinguish the term digit
used in the “Internal Digital Sorting” as separate from the radix, which need
not be 10-based, but they note that it is “…important to recognize that the
positions in the control field must be considered in order from ”least
significant” to ”most significant” with each one requring[sic] a separate
pass”(Friend, 1956, p.18). In (Rahman and Raman, 2002) the standard 8-bit
radix is used and the terms least/most-significant-bit (LSB/MSB) first radix
sort are considered, which appears to be easily interpreted as least or most-
significant-byte, so that form occurs regularly. After 2000, least/most-
significant-digit (LSD/MSD) radix sort, the terms we use in this paper,
becomes more common in spite of Friend’s careful admonishment on the use of
language; in particular Cormen et al., p.211’s analysis of the linear time of
LSD radix sort (though they do not call it exactly that) refers to “$d$-digit
numbers in which each digit can take on up to $k$ possible values.”

There is occasional conflation with counting sort, though when only one pass
is involved, they are the same. Bucket sort is occasionally conflated with MSD
radix sort. In some works, for example (Lee et al., 2002), one sees the terms
“left-to-right” or “right-to-left” sort that appear to be based on a big-
endian interpretation of the actual encoding. (Maus, 2002) refers to MSD radix
sort as also being “top down radix sort”, with the “left radix” (again, big-
endian) variety being a “bottom up radix sort.” (Lee et al., 2002) and (Al-
Badarneh and El-Aker, 2004) both refer to LSD radix sort as a “straight radix
sort”. (Lorin, 1975) and (Andersson and Nilsson, 1994) use the term
“lexographic sort,” though this is implied as being for variable key sized
inputs (e.g. variable length strings), and it is generally of the MSD variety,
as (Andersson and Nilsson, 1994) also refer to it as “forward radix sort” and
note that it is often best to divert to comparison-based-sorts when buckets
become small.

The terminology used for the components of a radix sort is more consistent,
with most variance being early on or based on the appearance of terms when
they are considered important, and their absence otherwise. Pass is used
consistently going as far back as (Friend, 1956). (Friend, 1956) refers to the
“recipient area” for a particular position, but this is synonymous to what is
called the buffer in most latter papers; they do note that the “sending area”
(what is now referred to as input, in general) and the “recipient area”
alternate with each pass in out-of-place implementations, which are the
dominant variety of radix sort. (Cormen et al., 2009) refers to bins. The
majority of other papers refer to the specific target area for elements with
like digits as buckets.

Counting, or frequency counting is quite consistent. Where terminology around
counting varies is in consideration of performance impacts of those counts, so
papers that concern themselves with these issues include mention of more
details related to the counting, such as the initialization of the counts and
the “prefix sum” of the counts(Rahman and Raman, 2002), and not just the count
itself, which is their primary focus.

The term for dealing varies significantly across papers. (Friend, 1956, p.18)
uses the term disbursement. (Rahman and Raman, 2002) calls this a permutation
phase, what appears to be a nod to in-place MSD radix sort’s unstable approach
to dealing. There are also some papers that refer to the dealing phase as the
distribution phase, as radix sorts fall into the class of distribution sorts;
(Cormen et al., 2009) describes how mechanical card-sorting machines would “…
examine a given column of each card in a deck and distribute the card…”(Cormen
et al., 2009, p.197) though they do not use the term deal that one might
readily associate with deck and card. It is difficult to say if the term deal
is the most common term in the literature, but it is in the running.

Table 1. This table identifies literature where improvements have been made or discussed, organized by parts of the radix sort where improvements have been applied and specific areas of research that affect the radix sort as a whole. Whether the improvement was specific to LSD or MSD radix sorts was specified for clarity. Topic | Source
---|---
Count |  (Friend, 1956), (Rahman and Raman, 2002), (Thiel et al., 2016), (Thiel, 2019), (Kumar, 2019), (Hanel et al., 2020)
Deal |  LSD: (LaMarca and Ladner, 1999), (Rahman and Raman, 2002), (Satish et al., 2009), (Adinets and Merrill, 2022) MSD: (Jiménez-González et al., 2003), (Kokot et al., 2017), (Kokot et al., 2018), (Hanel et al., 2020)
Diversion |  LSD: (Sedgewick, 1998), (Al-Badarneh and El-Aker, 2004), (Thiel, 2019) MSD: (Maus, 2002), (Jiménez-González et al., 2003), (Kokot et al., 2017), (Kokot et al., 2018)
Parallelization |  (Lee et al., 2002), (Satish et al., 2009), (Delorme et al., 2013), (Adinets and Merrill, 2022)
Variable Radix Size |  (LaMarca and Ladner, 1999), (Maus, 2002)
Variable Key Size |  (Paige and Tarjan, 1987), (Mcllroy et al., 1993), (Andersson and Nilsson, 1994), (Bentley and Sedgewick, 1997), (Nilsson, 2000), (Kärkkäinen and Rantala, 2009)

##  2\. Preliminaries

A key subset of sorting algorithms are those that sort based on the radix of
data instead of by using comparisons. Each element in an input has a key that
is considered in terms of a radix to break it into digits. These sorts achieve
linear asymptotic complexity in all cases, and are thus generally more
favorable when dealing with data that can be broken conveniently into digits.

The four common terms used by this approach are pass, bucket, deal and count.
We pass through our input. During a pass we process each element, and this
process is either a count or a deal, with both the dealing and the counting
based on a specific digit of the key. During a single pass we only deal based
on one digit. The count gathers information about the occurrence of a digit in
the data in order to build buckets that we will deal into during a subsequent
pass. We deal records into buckets, which means somehow moving them into a
different location. Our initial input can be considered a single bucket, but
after the first pass where we deal, the input is organized into buckets, and
we will process the records in subsequent passes based on which buckets they
are in.

A simple form of radix sort is a counting sort, sometimes referred to as a
bucket sort, though Knuth does not make the distinction between radix and
bucket (or digital) sorts (Knuth, 1998, p.169). In this algorithm, we consider
that there is a single digit and we can pass through the input twice, counting
the occurrence of each key in the first pass, and dealing into buckets during
the second pass. In between the two passes, we convert the count into indices
where records should be placed. With the indices established, a second pass
through the input places everything in order. The areas where records will be
placed, where the indices point, are referred to as buckets.

Consider an input whose contents were $N$ records with random keys between $0$
and $9$. To put this in order one could create an array $count[10]$,
initialize its contents to zero, then examine each element in $N$, increasing
the count $count[N[i].key]++$ for each element in $N$. This set of counts
could then be converted into a set of indices where records with the same keys
would be directed to the same index. The $0^{th}$ index would point to enough
space to hold all the records that have a $0$ for the given digit, the next
index would point to enough space to hold all the records that have a $1$ and
so forth. Passing through the input once more, each record would be placed
into the space allocated based on the counts, incrementing the corresponding
index so it is pointing to the next free space $buff[index[N[i].key]++]=N$.
When complete, the records would be sorted based on their keys, and each index
would point to the space immediately following the records that had been
placed.

###  2.1. Least Significant Digit (LSD) Radix Sorts

Algorithm 1 Least Significant Digit (LSD) Radix Sort

1: initialize space, including a buffer the same size as the input

2: for each record do

3: count the occurrence of the least significant digit

4: end for

5: convert counts of first digit into bucket indices in the empty buffer space

6: for each digit save the last, starting with the least significant do

7: for each record do

8: count the occurrence of the next least significant digit

9: deal into buckets specified for the current digit

10: end for

11: convert counts for the next digit into bucket indices in available space

12: end for

13: for each record do

14: deal into buckets specified for the last digit

15: end for

LSD sorts stably process data from the least significant digit to the most
significant digit, performing fewer additional operations per input value,
relative to a Most Significant Digit (MSD) approach that does not divert. This
is how most people sort a deck of playing cards, dealing into piles with like
face values before stacking those up in order and dealing into piles of like
suit.

Historically, improvements to LSD radix sorts focus on improving cache-
locality and skipping unneeded passes. Cache-locality improvements are
generally architecture dependent, whereas pass-skipping is dependent on data
distribution. Modern work on LSD radix sorts, particularly out of NVIDA,
focuses on the application of LSD radix sorts parralleized across GPUs (Satish
et al., 2009; Adinets and Merrill, 2022).

The pseudocode in Algorithm 1 shows a traditional approach to LSD radix sort.
An initial counting pass is performed (lines 2 \- 5). For each digit but the
last, each record is dealt into its identified space in the buffer, counting
the subsequent digit and prepare the next round of counts (lines 6 \- 11). The
final digit is dealt back into the original space (line 14).

###  2.2. Most Significant Digit (MSD) Radix Sorts

MSD radix sorts process each digit, starting at the most significant digit.
For each digit processed, each record will be passed over twice. In the first
pass, the digits being considered are counted in order to create buckets of
the appropriate size for dealing. In the second pass, the digit being
considered is used to identify the properly sized bucket to deal the record
into. After all digits are dealt, each bucket is sorted recursively. After
each digit is processed anything in a bucket will be greater than anything in
all buckets to its left and less than anything in all buckets to its right. By
the time the least significant digit is processed, each bucket for that digit
will either be empty or will have one or more records with the same key, and
thus is entirely sorted.

A common variant on MSD radix sort is to apply diversion, what we describe as
a diverting MSD radix sort. Unless all records are dealt into the same bucket,
when each bucket is processed recursively it will be smaller than the original
input. Processing smaller buckets is more cache efficient as buckets fit into
ever more local levels of cache. Further, if buckets become small enough,
diverting to another algorithm can be efficient, with anywhere from 10 to 50%
of total time of the sort involving the sorting of these tiny buckets Kokot et
al. (2018). Adinets and Merrill also identify that diversion allows parallel
variants of MSD radix sort to perform significantly fewer passes as opposed to
their proposed LSD algorithm(Adinets and Merrill, 2022).

MSD radix sort can be implemented in-place, but will not be stable due to the
way records are swapped around. The in-place process swaps records from the
free spots in identified buckets, from left to right, into the first free
spots in their destination buckets. The next record from the first free spot
in the left-most bucket with free spots is then processed. The savings in
space come at a potential increase in cost to access records and a guaranteed
increase in swaps. Ska Sort, implemented by Skarupke (2017), is an example of
a diverting MSD radix sort that makes use of many technical optimizations and
ease-of-use measures to demonstrate that this type of sort is still relevant
when memory must be conserved. As the performance times of in-place
implementations are generally slower than corresponding ${\Theta}(n)$ space
implementations, and we will not consider in-place algorithms any further.

Algorithm 2 Diverting Most Significant Digit (MSD) Radix Sort

1: if this is the initial call then

2: initialize buffer the same size as the input

3: use the most significant digit as the current digit and the input as the
current bucket

4: end if

5: initialize space for counters/indices

6: if the current bucket is small enough to divert then

7: run a diversion algorithm on the current bucket

8: return

9: end if

10: for each record in the current bucket do

11: count the occurrence of the current digit

12: end for

13: convert counts for the current digit into bucket indices in available
space

14: for each record in the current bucket do

15: deal into buckets specified for the current digit

16: end for

17: if the current digit is the least significant digit then

18: return

19: end if

20: for each bucket do

21: recurse using the current bucket, the next most significant digit and the
available space

22: end for {Ensuring the details of where the sorted values end up is an
implementation detail.}

MSD radix sorts can be stable when implemented using an extra buffer of the
same size as the input. However, diverting using a non-stable algorithm
removes the stability guarantee. Kokot et al. (2018) uses sorting networks to
divert in their RADULS2 implementation, and thus give up stability for the
performance gains of their diversion algorithm. Using a Merge Sort as a
diversion algorithm would be stable, and would take advantage of the already
existing buffer space. Using Insertion Sort would be stable and would perform
very well, but only on the smallest buckets.

Algorithm 2 shows diverting MSD radix sort pseudocode given records in
$input$, a suitably sized extra buffer in $buf$ and a set of digits ordered
least-to-most significant in $digits$. The pseudocode starts with a check for
the input being smaller than the diversion threshold (line 6), in which case
it diverts to Insertion Sort (line 7) or some other diversion algorithm. For
the current most-significant digit, the algorithm counts the different values
for that digit in the input and then deals into buckets located in $buf$
(lines 10–12). The algorithm is then recursively called on each bucket, using
the now available space in the original input and no longer considering the
now-counted digit (lines 20–22). The recursion ends all digits have been
processed (lines 17–19).

For simplicity sake, the pseudocode in Algorithm 2 hides some details that are
less immediately relevant to the general flow of the algorithm. The pseudocode
leaves out the creation and management of buckets between the count and deal
functions (lines 10–12). The pseudocode hides any fiddling with offsets which
may need to be done when reusing the buffer during recursion (lines 20–22).
Lastly, the pseudocode of Algorithm 2 hides the assumption that there are an
even number of digits, though an actual implementation would account for
either possibility.

Historically, improvements to MSD radix sorts focus on parallelization,
improving cache-locality and effective diversion algorithms, with examples of
all these improvements in Kokot et al. (2018), and an example of diversion
with an in-place and broadly usable implementation in Skarupke (2017).
Parallelization and cache-locality improvements are generally architecture
dependent, whereas diversion is more dependent on data distribution.

##  3\. Diverting Fast Radix

In previous work, Thiel et al. (2016), we have demonstrated the improvements
of Simple Fast Radix (referred to as “Fast Radix” in that work) over radix
sorts, showing that omitting the first counting pass leads to faster
performance. Here we will show how this approach can be extended and how
diversion can reliably be applied to any LSD radix sort in order to create a
competitively performant algorithm.

Fast Radix performs its dealing passes using estimations instead of counting
for all but the last pass, including an overflow processing step in between
passes. Each pass estimates bucket sizes, deals into buckets or the overflow
and then processes overflow. Subsequent passes deal from each predicted bucket
and then from its overflow counterpart if necessary, such that the algorithm
still performs the next deal in stable, bucket order. The algorithm counts the
last digit in the penultimate pass to ensure exact buckets in the last pass,
avoiding an extraneous copy pass and putting everything in order.

Algorithm 3 Fast Radix Sort

1: initialize space, including a buffer the same size as the input

2: estimate bucket sizes for each digit but the most significant

3: for each digit but the last, starting at the least significant do

4: for each record, in bucket order do

5: if this is the next to last digit then

6: count the occurrence of the most significant digit

7: end if

8: if the bucket this record would be dealt to is full then

9: deal into overflow space

10: count the occurrence for the overflown digit

11: else

12: deal into the estimated bucket

13: end if

14: end for

15: if the existing overflow space is not large enough then

16: create more overflow space

17: end if

18: convert counts for overflown digits into bucket indices in available
overflow space

19: for for each overflown record do

20: deal into buckets specified for the overflown digits {virtually extending
overflown buckets}

21: end for

22: end for

23: convert counts for the most significant digit into bucket indices in
available space

24: for each record in bucket order do

25: deal into buckets specified for the last digit

26: end for

Fast Radix differs from Simple Fast Radix in that only the last pass gets
exact counts, with all other passes using estimates. While Simple Fast Radix
could deal overflow directly to the beginning of input during its initial
estimation pass, as there was no opportunity for conflict, Fast Radix cannot
make that assumption on subsequent passes. When overflow found on a pass
exceeds the overflow space already available in overflow, the Fast Radix
algorithm creates a new overflow buffer that can hold twice the amount of
overflow found. The algorithm can then assign overflow buckets to the first
half of overflow and overflow from the subsequent pass can be dealt into the
second half, with the algorithm dealing any excess safely into the beginning
of input. The algorithm avoids a conflict where dealing overflow would
overwrite parts of input that it had not yet dealt out by first dealing
overflow into an overflow space that is at least as big as the existing
overflow buckets.

Diverting LSD radix sort estimates the number of most significant digits over
which an LSD radix sort must run to achieve optimal diversion for the
remaining partially ordered records. After an LSD radix sort runs on those
high-order digits, the diverting LSD radix sort will determine if any large
buckets remain. Small buckets are sorted in bulk using a simple, stable sort,
most readily Insertion Sort. Any large buckets are sorted recursively with
diverting LSD radix sort.

Algorithm 4 Diverting Fast Radix (DFR) Sort

1: if this is the initial call then

2: initialize buffer the same size as the input

3: end if

4: estimate the top passes to yield advantageous diversion

5: initialize space for counters/indices

6: perform LSD radix sort only on top passes using the existing buffer/counter
space

7: while there are more records in the input do

8: scan for large buckets

9: if we find a large bucket then

10: perform a simple diversion sort on previous small buckets

11: recurse on the large bucket

12: end if

13: end while

14: perform a simple diversion sort on previous small buckets {Ensuring the
details of where the sorted values end up is an implementation detail.}

Diverting LSD radix sort differs from diverting MSD radix sort. It estimates
the number of passes required to reach optimal diversion and then performs an
LSD radix sort on those passes, recursing on remaining large buckets or
diverting on runs of small buckets. Conversely, diverting MSD Radix sort runs
recursively on each level at an ever-increasing cost before diverting. Suppose
the algorithms use the same diversion threshold. In that case, we have shown
that diverting LSD radix sort can accurately estimate the average number of
passes for input with random records from a uniform distribution.

###  3.1. Analysis

In considering the average number of passes performed by diverting radix
sorts, we will consider diverting MSD Radix Sort and our proposed diverting
LSD Radix Sort. Current research on MSD Radix Sorts consider diversion as a
means to reduce cost(Kokot et al., 2018). There is limited occurrence of any
diverting LSD Radix Sorts, with the summary consideration of diversion after a
constant number of passes(Sedgewick, 1998, Ch.10) and a more recent approach
involving LSD Radix Sorting floating point numbers that proposed diverting
after “$Math.min(\log{2}{n+9},32)$”(Maus, 2019) bits. A consideration of the
model for passes performed by both diverting MSD Radix Sort and our diverting
LSD Radix Sort, where diversion is based on bucket size, shows that the
average number of passes performed is the same in either case.

The concept of diversion is relevant to hybrid sort algorithms which utilise
multiple sort algorithms at the different passes during recursion (Kokot et
al., 2018). Diversion strategies in radix sort are predicated on a threshold
of a bucket size $N_{d}$. If the size is below the threshold then the
recursive call to radix sort on that bucket defaults to a different sort
algorithm. For our analysis, we will assume the diverting sort algorithm is
non-recursive so does not involve further passes. We deal with the non-
diverting case by regarding it as adopting a null diversion strategy which
never diverts from the radix sort.

Figure 1. A comparison of ThielSort and RADULS2 on uniformly distributed
integers with random inputs ranging from $N==256^{3}\sim=16$ million to
$N\sim=256^{3.75}==1$ billion. The y-axis is time in seconds and the x-axis is
$\log_{256}(N)$. At $N==1$ billion, RADULS takes $\sim 50\%$ longer.

With respect to the cost of counting, we note that traditionally it must be
done for each $N$ for every pass, with the cost of managing counting in MSD
radix sorts growing exponentially with $M$, a significant cost noted in
(Adinets and Merrill, 2022). Conversely, estimating bucket sizes and
considering overflow as a fraction of $N$, which avoids the cost of counting,
can be modeled recursively as $g_{o}(N,M)$ for random inputs with uniform
distributions. $N$ is the size of the input and $M$ is the number of different
buckets one might deal into (the digit size).

(1)
$g_{o}(N,M)=g_{o}(N-1,M)-\frac{1}{M}g_{o}(N-1,M,\left\lceil\frac{N}{M}\right\rceil-1)$

Where $g_{o}(N,M,k)$ is defined as:

(2)
$g_{o}(N,M,k)=(\frac{N}{M}-k)\frac{M}{N}\binom{N}{k}(\frac{1}{M})^{k}(1-\frac{1}{M})^{N-k}$

The overflow fraction $g_{o}(N,M)$ is clearly monotonically decreasing and can
experimentally be shown to be a very close to zero even while $N$ is on the
order of $10^{3}$. Those passes that use estimation, on average, save
effectively the entire cost of counting, and the final pass that does count,
allowing records to be placed in their sorted positions, costs no more than it
would in traditional LSD or MSD radix sorts. We note also that in the event
that overflow estimation is entirely wrong, it still only performs a glorified
counting pass, costing comparably to traditional LSD and MSD counting passes.

While MSD radix sorts do not readily allow counting to be replaced, these
results demonstrate that counting is superfluous for all passes but the last
for LSD radix sorts with inputs whose least significant digits are effectively
randomized according to a uniform distribution.

Considering the Marsaglia polar method (or the Knuth polar implementation) for
generating values from a normal distribution, we make the conjecture that each
digit of subsequently lower significance comes closer to being uniformly
distributed, making estimation naturally more effective on the least
significant digits even with normal distributions so long as they do not have
very small standard deviations.

Figure 2. A comparison of ThielSort and RADULS2 on normally distributed
integers with random inputs ranging from $N==256^{3}\sim=16$ million to
$N\sim=256^{3.75}==1$ billion. The average values for the random inputs is
$\sim 2^{63}$ with a standard deviation of $\sim 2^{61}$, what we consider a
wide normal distribution. The y-axis is time in seconds and the x-axis is
$\log_{256}(N)$. At $N==1$ billion, RADULS takes $\sim 25\%$ longer.

##  4\. Implementation

In implementing ThielSort, optimizations were made to render it competitive
with RADULS2(Kokot et al., 2017), the current fastest radix sort. Cache-
sensitive random writes, block writes and a reduction on overflow checks were
the three main technical optimizations applied.

As in most current radix sort research, tackling the cost of random writes
during the dealing pass of when inputs are large was required. While ThielSort
does have approaches that make use of SIMD ”non-temporal” writes, similarly to
RADULS2(Kokot et al., 2017), in the primary range of data presented, we have
found that simple pre-fetching (for writes) was sufficient to yield the
desired performance gains. There are still ranges of very large ($n$ greater
than 1 billion on many machines tested) where the non-temporal write strategy
is more effective than pre-fetching. As well, with smaller input sizes, the
cost of pre-fetching renders it slower than a vanilla deal. ThielSort’s
current implementation focuses on the performance in the 10s of millions to
one billion range to demonstrate that with such technical optimizations DFR
can outclass diverting MSD radix sorts.

In managing overflow, some dealing is performed on blocks of elements. While
in uniform distributions of data and even many other distributions, very
little time is spent on this process, in the cases were overflow is a more
significant issue, taking advantage of faster block-writing significantly
reduces the cost of these writes and modern processors have operations that
support this. ThielSort makes use of this within its ”ladle” system for
dynamic destinations to further reduce the overhead cost of managing overflow.

A significant concern noted with DFR was that it replaced counting with an
overflow check for each record on each pass. This cost can be significantly
amortized across a large number of attempts by noting that if each buck starts
with $\frac{N}{256}$ empty space, one can deal that many records before
performing 256 bucket checks. Each subsequent round of dealing would be able
to deal without checks so long as the dealing amount did not exceed the
smallest bucket noted in the last round of bucket checks. When the smallest
amount of space is sufficiently reduced that this trade-off is no longer
advantageous, either counting every time or using some other technique can
finish off. ThielSort uses this approach, using its dynamic ”ladle” system and
block writes when overflow becomes common. As overflow is rare, the vast
majority of data is processed with a constant number of checks, rendering this
cost at least constant for most data and vanishingly small for the rest for
buckets with good overflow estimates. On the order of only
$\frac{2N}{256^{2}}<1\%$ records would remain after 256 such amortized checks
if bucket estimates were accurate.

##  5\. Experiments

The experiment examines the average case performance of ThielSort, our current
implementation of the Diverting Fast Radix. Results contrast ThielSort with
the RADULS2 implementation of their corresponding algorithm, as it is
currently the fastest published radix sort available for large n on a single
processor. std::sort is orders of magnitude worse and not shown.

Tests were run on inputs of varying sizes from $N\sim=10^{3}$ to
$N\sim=10^{9}$. Uniform distributions and a variety of normal distributions
were considered, reflecting common testing approaches. Zipf distributions were
also considered as Kokot et al. have demonstrated their implementation’s
performance against those distributions, and the results reflect the narrow
normal distributions shown in Figure 3 while generating their input sequences
much more rapidly.

While the performance reported by the RADULS2 implementation, as provided by
their authors, does not include setup steps that are traditionally included in
timing tests, such as allocating the extra buffer, results presented herein
have included these steps in RADULS2 timings by adjusting where timings start
to reflect all this work.

Figure 3. A comparison of ThielSort and RADULS2 on normally distributed
integers with random inputs ranging from $N==256^{3}\sim=16$ million to
$N\sim=256^{3.75}==1$ billion. The average values for the random inputs is
$\sim 2^{32}$ with a standard deviation of $\sim 2^{30}$, what we consider a
narrow normal distribution. The y-axis is time in seconds and the x-axis is
$\log_{256}(N)$. At $N==1$ billion, RADULS takes $\sim 23\%$ longer.

Tests were run on a variety of available platforms, including HPC computing
systems, with Concordia’s speed HPC system being used for the test results
shown here. This implementation of the Diverting Fast Radix is single-threaded
and thus results are compared against RADULS2 running using a single thread.
In the experimental results presented here, for each input size, three random
inputs were generated and the timing results were averaged for each sort. Jobs
were queues on the HPC system with sufficient ram to avoid paging to disk and
each job ran alone on a single thread.

For larger input sizes, ThielSort performed significantly better in all
distributions. While a wide normal distributions adversely affect ThielSort,
it similarly adversely affects any other radix sort, as shown here with
RADULS2. On uniform distributions, RADULS2 took $\sim 50\%$ longer with inputs
of size $1$ billion. On wide normal distributions, RADULS2 took $\sim 25\%$
longer with inputs of size $1$ billion, which we interpret as a benefit from
more early diversion, but still not making up for the extra work of an MSD
radix sort. On narrow normal distributions RADULS2 took $\sim 25\%$ longer
with inputs of size $1$ billion, with both implementations taking longer as
the distribution effectively forcing both sorts to perform more passes.

While RADULS2 still performed well on smaller input sizes, this is a function
of the difference in technical optimizations used at those input-sizes;
ThielSort presents a practical demonstration that the implementation of the
Diverting Fast Radix can yield the predicted benefits over diverting MSD radix
sorts.

##  6\. Conclusion and Future Research

We presented ThielSort, an implementation of the Diverting FastRadix
algorithm. The experimental results comparing ThielSort with RADULS2, the
current fastest radix sort, show that the Diverting Fast Radix algorithm can
be effectively implemented with technical optimizations in the same class as
existing diverting MSD radix sorts to outperform them. The algorithmic
improvements offered by Diverting Fast Radix are independent of hardware and
exploring the affordances of the Diverting Fast Radix algorithm should allow
for even more effective implementations.

Current literature highlights that current radix sort optimization focuses on
parallel implementations and writing in cache efficient manners. We have
already demonstrated some cache efficient approaches with this implementation,
but determining optimal approaches for a given system at compile-time is an
interesting avenue of future work. Parallel implementations would be
advantageous, and applying previous advances in the area of LSD Radix Sorts to
the Diverting Fast Radix algorithm is an obvious future target. Adinets and
Merrill’s recent pre-print(Adinets and Merrill, 2022) highlights how actively
pursued parallel LSD radix sorts are and their work represents an optimal
avenue to bring ThielSort and the Diverting Fast Radix algorithm into
mainstream usage.

The Diverting Fast Radix algorithm performs estimations to achieve some of its
gains. Even the trivial estimations currently used offers very tangible
improvements. However, approaches like sampling should be provably and
efficiently able to increase performance further with non-uniform
distributions. Strategies, including sampling, for dealing with smaller input
sizes would make implementations of this algorithm even more broadly
applicable.

Even while considering future improvements relating to implementations of our
algorithm, we argue that these results already conclusively identify that LSD
radix sorts are superior to MSD radix sorts when sorting data that is
randomized according to a uniform distribution. Future research may support
the conjecture that Diverting LSD radix sorts are superior to Diverting MSD
radix sorts in all cases, i.e. there exists a corresponding LSD radix sort
algorithm that provides a faster result, on average, for any MSD radix sort
algorithm. We hope that this returns the attention of those researchers
developing technical optimizations in this area to LSD radix sorts.

## References

  * (1)
  * Adinets and Merrill (2022) Andy Adinets and Duane Merrill. 2022.  Onesweep: A Faster Least Significant Digit Radix Sort for GPUs.  _arXiv preprint arXiv:2206.01784_ (2022).
  * Al-Badarneh and El-Aker (2004) Amer Al-Badarneh and Fouad El-Aker. 2004.  Efficient Adaptive In‐Place Radix Sorting.  _Informatica_ 15, 3 (2004), 295 – 302.  https://lib-ezproxy.concordia.ca/login?url=https://search.ebscohost.com/login.aspx?direct=true&db=a9h&AN=18106146&site=ehost-live&scope=site
  * Andersson and Nilsson (1994) Arne Andersson and Stefan Nilsson. 1994.  A new efficient radix sort. In _Proceedings, 35th Annual Symposium on Foundations of Computer Science, 1994_. IEEE, New York, NY, USA, 714–721.  https://doi.org/10.1109/SFCS.1994.365721
  * Bentley and Sedgewick (1997) Jon L. Bentley and Robert Sedgewick. 1997.  Fast Algorithms for Sorting and Searching Strings. In _Proceedings of the Eighth Annual ACM-SIAM Symposium on Discrete Algorithms_ (New Orleans, Louisiana, USA). ACM, Philadelphia, PA, USA, 360–369.
  * Cormen et al. (2009) Thomas H. Cormen, Charles E. Leiserson, Ronald L. Rivest, and Clifford Stein. 2009.  _Introduction to Algorithms, 3rd Edition_.  MIT Press, Cambridge, MA, USA.  http://mitpress.mit.edu/books/introduction-algorithms
  * Delorme et al. (2013) Michael C. Delorme, Tarek S. Abdelrahman, and Chengyan Zhao. 2013.  Parallel Radix Sort on the AMD Fusion Accelerated Processing Unit. In _2013 42nd International Conference on Parallel Processing_. IEEE, New York, NY, USA, 339–348.  https://doi.org/10.1109/ICPP.2013.43
  * Friend (1956) Edward H. Friend. 1956\.  Sorting on Electronic Computer Systems.  _J. ACM_ 3, 3 (July 1956), 134–168.  https://doi.org/10.1145/320831.320833
  * Hanel et al. (2020) Carson Hanel, Arif Arman, Di Xiao, John Keech, and Dmitri Loguinov. 2020.  _Vortex: Extreme-Performance Memory Abstractions for Data-Intensive Streaming Applications_.  Association for Computing Machinery, New York, NY, USA, 623–638.  https://doi.org/10.1145/3373376.3378527
  * Jiménez-González et al. (2003) Daniel Jiménez-González, E Guinovart, J-L Larriba-Pey, and Juan J Navarro. 2003.  CC-Radix: a cache conscious sorting based on Radix sort. In _Parallel, Distributed and Network-Based Processing, 2003. Proceedings. Eleventh Euromicro Conference on_. IEEE, New York, NY, USA, 101–108.  https://doi.org/10.1109/EMPDP.2003.1183573
  * Kärkkäinen and Rantala (2009) Juha Kärkkäinen and Tommi Rantala. 2009.  Engineering Radix Sort for Strings. In _String Processing and Information Retrieval_ , Amihood Amir, Andrew Turpin, and Alistair Moffat (Eds.). Springer Berlin Heidelberg, Berlin, Heidelberg, 3–14.
  * Knuth (1998) Donald E. Knuth. 1998\.  _The Art of Computer Programming_.  Number v. 3 in Addison-Wesley series in computer science and information processing. Addison-Wesley, Boston, MA, USA.  https://books.google.ca/books?id=CH1GAAAAYAAJ
  * Kokot et al. (2017) Marek Kokot, Sebastian Deorowicz, and Agnieszka Debudaj-Grabysz. 2017\.  Sorting Data on Ultra-Large Scale with RADULS. In _Beyond Databases, Architectures and Structures. Towards Efficient Solutions for Data Analysis and Knowledge Representation_ , Stanisław Kozielski, Dariusz Mrozek, Paweł Kasprowski, Bożena Małysiak-Mrozek, and Daniel Kostrzewa (Eds.). Springer International Publishing, New York, NY, USA, 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_ , Aleksandra Gruca, Tadeusz Czachórski, Katarzyna Harezlak, Stanisław Kozielski, and Agnieszka Piotrowska (Eds.). Springer, New York, NY, USA, 481–491.  https://doi.org/10.1007/978-3-319-67792-7_47
  * Kumar (2019) Ravin Kumar. 2019\.  Modified Counting Sort.  In _System Performance and Management Analytics_ , P. K. Kapur, Yury Klochkov, Ajit Kumar Verma, and Gurinder Singh (Eds.). Springer Singapore, Singapore, 251–258.  https://doi.org/10.1007/978-981-10-7323-6_21
  * LaMarca and Ladner (1999) Anthony LaMarca and Richard E Ladner. 1999.  The Influence of Caches on the Performance of Sorting.  _Journal of Algorithms_ 31, 1 (1999), 66 – 104.  https://doi.org/10.1006/jagm.1998.0985
  * Lee et al. (2002) Shin-Jae Lee, Minsoo Jeon, Dongseung Kim, and Andrew Sohn. 2002\.  Partitioned Parallel Radix Sort.  _J. Parallel and Distrib. Comput._ 62, 4 (2002), 656 – 668.  https://doi.org/10.1006/jpdc.2001.1808
  * Lorin (1975) Harold Lorin. 1975\.  _Sorting and Sort Systems (The Systems Programming Series)_.  Addison-Wesley Longman Publishing Co., Inc., Boston, MA, USA.
  * Maus (2002) Arne Maus. 2002\.  ARL, a faster in-place, cache friendly sorting algorithm. In _Norsk Informatik Konferranse NIK’2002_. NIKT, Norway, 85–95.
  * Maus (2019) Arne Maus. 2019\.  RadixInsert, a much faster stable algorithm for sorting floating-point numbers. In _Norsk IKT-Konferanse for Forskning og Utdanning_. NIKT, Norway.
  * Mcllroy et al. (1993) Peter M Mcllroy, Keith Bostic, and Malcolm Douglas Mcllroy. 1993\.  Engineering Radix Sort.  _Computing Systems_ 6, 1 (1993), 5–27.
  * Nilsson (2000) Stefan Nilsson. 2000\.  The Fastest Sorting Algorithm?  _Dr. Dobb’s journal (1989)_ 25, 4 (2000), 38–+.
  * Paige and Tarjan (1987) Robert Paige and Robert E Tarjan. 1987.  Three partition refinement algorithms.  _SIAM J. Comput._ 16, 6 (1987), 973–989.
  * Rahman and Raman (2002) Naila Rahman and Rajeev Raman. 2002.  Adapting Radix Sort to the Memory Hierarchy.  _ACM J. Exp. Algorithmics_ 6 (dec 2002), 7–es.  https://doi.org/10.1145/945394.945401
  * Satish et al. (2009) Nadathur Satish, Mark Harris, and Michael Garland. 2009\.  Designing efficient sorting algorithms for manycore GPUs. In _2009 IEEE International Symposium on Parallel Distributed Processing_. IEEE, New York, NY, USA, 1–10.  https://doi.org/10.1109/IPDPS.2009.5161005
  * Sedgewick (1998) Robert Sedgewick. 1998\.  _Algorithms in C++_.  Number v. 1-2 in Algorithms in C++. Addison-Wesley, Boston, MA, USA.  https://books.google.ca/books?id=331YAAAAYAAJ
  * Skarupke (2017) Malte Skarupke. 2017\.  Ska Sort.  https://github.com/skarupke/ska_sort.
  * Thiel (2019) Stuart Thiel. 2019\.  Diverting Fast Radix.  https://github.com/ramou/dfr
  * Thiel et al. (2016) Stuart Thiel, Greg Butler, and Larry Thiel. 2016\.  Improving GraphChi for Large Graph Processing: Fast Radix Sort in Pre-Processing. In _Proceedings of the 20th International Database Engineering & Applications Symposium_ (Montreal, QC, Canada) _(IDEAS ’16)_. Association for Computing Machinery, New York, NY, USA, 135–141.  https://doi.org/10.1145/2938503.2938554