# RadixSort implementation design and performance

Credits for the original implementation I based my code on:
Quora - What is the most efficient way to sort a million 32-bit integers?

## Implementation

// RadixSort - works for values up to unsigned_int_max (32-bit)
template<typename Iter>
typedef typename iterator_traits<Iter>::value_type value_type;
vector<value_type> out(__last - __first);

// calculate most-significant-digit (256-base)
value_type __mx = *max_element(__first, __last); //O(n)
int __msb = 0;
do {
__msb += 8;
__mx = __mx >> 8;
} while(__mx);

Iter __i, __j, __s;
bool __swapped = false;
for (int __shift = 0; __shift < __msb; __shift += 8) {

// cycle input/auxiliar vectors
if (__swapped) {
__i = out.begin();
__j = out.end();
__s = __first;
}
else {
__i = __first;
__j = __last;
__s = out.begin();
}

// counting_sort
size_t count[0x100] = {};
for (Iter __p = __i; __p != __j; __p++)
count[(*__p >> __shift) & 0xFF]++;

// prefix-sum
size_t __m, __q = 0;
for (int i = 0; i < 0x100; i++) {
__m = count[i];
count[i] = __q;
__q += __m;
}

// filling result
Iter __v;
for (Iter __p = __i; __p != __j; __p++) {
__v = __s + count[(*__p >> __shift) & 0xFF]++;
*__v = *__p;
}
__swapped = !__swapped;
}

// if ended on auxiliar vector, copy to input vector
if (__swapped) copy(out.begin(), out.end(), __first); //O(n)
}


## Discussion

In order to implement a more general algorithm than the one used as reference, I tried to refactor the code using template and iterators. This is the main difference between my code and the original.

I can't simply swap references when using iterators, so I approached this problem by cycling the iterators used in each counting-sort loop. Doing so prevents the need to copy the output array on every loop.

Another different aspect of my code is that I find max-element and calculate the most-significant-digit on base-256. I use this information to determine how many counting-sort loops are necessary, instead of hardcoding 32 (unsigned_max) and always running the loop 4-times even if all values are less than 255. This actually adds unnecessary overhead if 4-loops are necessary, but it should increase execution time otherwise.

The temporary container used is a vector<value_type> and this is something I'm not quite sure about, I feel like this is a point where I could improve my code. I'd like to hear opinions about it.

## Questions

• What are possible improvements that can be made to my implementation?
• How would you re-factor the iterator cycle I used? (I hate it)
• Which container should I use for the temp-array?
• Should I use std::copy or std::move to move the data from the temp-array to the input container?

### Relevant Information

I decided to run some benchmarks to test if calculating the maximum-digit makes a huge difference. element_values defines the range of values each element can have. The maximum-digit only affects the number of loops if element_values can be represented with 24-bits or less. Here are the results for n = 1E7:

     element_values < 256 (8-bits):
radix_sort_msd - Average time:   51.59 ms
radix_sort_32  - Average time:  109.03 ms

element_values < UINT_MAX(32-bits):
radix_sort_msd - Average time:  107.38 ms
radix_sort_32  - Average time:   89.75 ms


While radix_sort_msd works really well when element_values is small, it really just depends on the dataset. Therefore implementing is a matter of preference.

• Why does radix_sort_32 take longer for element values < 256 than for element values < UINT_MAX? – rcgldr Feb 27 '17 at 3:43
• Assuming an array of near random data, the main overhead in radix sort is all the cache unfriendly random access writes during each radix sort pass, and that's the part that can't be optimized. In the case of 1 million integers, the array would fit in an 8MB L3 cache, but I don't know the impact of cache coherency (for the other cores) when writing to the cache and how writes are committed back to main memory. I don't see much improvement other than multi-threading, for k cores split the array into k parts, sort the k parts in parallel, then merge the k sorted parts. – rcgldr Feb 27 '17 at 3:51
• @rcgldr I noticed that. I double checked results and code. The honest answer is: "I don't know". Regarding the multi-threading solution, I don't fully understand how cache writes works in this scenario either and I'm not sure if I wanna go down this rabbit hole. – Guilherme Barros Avila Feb 27 '17 at 10:08
• I tested element values == 0, element values < 256, and element values <= UINT_MAX. For 1 million elements, times were 0.103 sec, 0.108 sec, 0.128 sec. In this case, the original and working arrays fit within the cache on my system, Intel 3770K 3.5ghz, 32KB L1 and 256KB L2 cache per core, 8MB L3 cache shared by all cores, so the random access (versus sequential) writes aren't creating much overhead. At 16 million integers, the random access writes are outside the range of cache, and times were 0.162 sec, 0.233 sec, 0.444 sec. – rcgldr Feb 28 '17 at 0:06
• Multi-threading: 4 threads on 4 cores at 1 million pseudo random 32 bit unsigned integers provided little gain, 0.128 sec => 0.115 sec. There was > 2x gain at 16 million integers, 0.444 sec => 0.200 sec. Again I assume this is an issue with the cache. – rcgldr Feb 28 '17 at 0:08

One of the things that you were wondering about was whether radix_sort_msd or radix_sort_32 was better. The former called max_element() in order to determine the max width of the values, and the latter always did 4 passes even when they weren't necessary.

I took radix_sort_32 and added this code after generating the counts:

// If this byte is zero for every value, skip the byte entirely.
if (count[0] == (size_t)(__last - __first))
continue;


This is a quick check to find out if the entire input was full of zeroes for the current byte. If so, it moves on to the next byte instead of wasting time making an exact copy of the input to the temp area. It still uses more time than radix_sort_msd in the "values < 256" case, but it is much faster, because for 3 passes it can skip more than 50% of the work.

On my machine, radix_sort_32 got 40% faster in the "values < 256" case compared to before the change. This might make you reconsider which version is better.

Another thing you could do to improve your radix sort is to use a so-called "hybrid" radix sort, which is a hybrid of MSD and LSD radix sorts. The idea is that use you one pass of radix sort on the most significant byte. This breaks the input into 256 "bins". Then you run LSD radix sort (your current sort) on each of the 256 bins.

The advantage of this hybrid sort over a plain LSD radix sort is that the hybrid version is more cache friendly. In the plain LSD sort, if your input doesn't fit in cache, then the second counting pass can't benefit from caching because by the end of the first counting pass, the start of the array will already be evicted from the cache.

With the hybrid sort, the input is first broken up into 256 parts (which hopefully do fit in the cache). So sorting each bin should be faster due to caching. However, there is extra overhead, because each of the 256 smaller sorts needs to do some fixed amount of work to deal with the counting bucket. So we should only do the hybrid sort if the input is larger than our cache size, otherwise we will only have more overhead and no benefit.

I modified your code to do the hybrid sort and it ran faster than the original version on randomized 32-bit input. I did not attempt to optimize it for the "values < 256" case, although there are definitely possibilities for improvements there, similar to what I mentioned in the previous section. Here is the code:

// LSD RadixSort helper function for radix_sort() below.
template<typename Iter>
static void radix_sort_lsd(Iter __first, Iter __last, Iter __out, Iter __outEnd,
int __msb, bool needsSwap)
{
Iter __i, __j, __s;
bool __swapped = false;

for (int __shift = 0; __shift < __msb; __shift += 8) {

// cycle input/auxiliar vectors
if (__swapped) {
__i = __out;
__j = __outEnd;
__s = __first;
}
else {
__i = __first;
__j = __last;
__s = __out;
}

// counting_sort
size_t count[0x100] = {};
for (Iter __p = __i; __p != __j; __p++)
count[(*__p >> __shift) & 0xFF]++;

if (count[0] == (size_t)(__last - __first))
continue;

// prefix-sum
size_t __m, __q = 0;
for (int i = 0; i < 0x100; i++) {
__m = count[i];
count[i] = __q;
__q += __m;
}

// filling result
for (Iter __p = __i; __p != __j; __p++) {
*(__s + count[(*__p >> __shift) & 0xFF]++) = *__p;
}
__swapped = !__swapped;
}

// if ended on auxiliar vector, copy to input vector
if (__swapped != needsSwap) copy(__out, __outEnd, __first); //O(n)
}

// If the input exceeds this threshold (in bytes), we do one pass of MSD
// followed by the rest of the passes done by LSD.  This number should be
// an estimate of the cache size.
#define THRESHOLD 8000000

// RadixSort - works for values up to unsigned_int_max (32-bit)
template<typename Iter>
{
typedef typename iterator_traits<Iter>::value_type value_type;
size_t len = (size_t) (__last - __first);
vector<value_type> out(len);

// First, test if the input exceeds the caching threshold.  If the input
// is smaller than the threshold, just do a straight LSD radix sort.
if (len * sizeof(value_type) < THRESHOLD) {
sizeof(value_type) * 8, false);
return;
}

// Set __shift to the most significant byte.
int  __shift = (sizeof(value_type)-1) * 8;
Iter __s     = out.begin();
Iter __p     = __first;

// counting_sort
size_t count[0x100] = {};

for (size_t i = 0; i < len; i++) {
count[(*__p++ >> __shift) & 0xFF]++;
}

// prefix-sum
size_t __m, __q = 0;
for (int i = 0; i < 0x100; i++) {
__m = count[i];
count[i] = __q;
__q += __m;
}

// filling result
__p = __first;
for (size_t i = 0; i < len; i++) {
*(__s + count[(*__p >> __shift) & 0xFF]++) = *__p++;
}

// For each of the 256 bins, do a LSD radix sort on the bin.  The input
// and auxiliary vectors have been swapped, so we pass needsSwap = true to
// indicate that the LSD sort should end on the auxiliary vector instead.
int startIndex = 0;
for (int i = 0; i < 0x100; i++) {
int endIndex = count[i];
radix_sort_lsd(__s + startIndex, __s + endIndex, __first + startIndex,
__first + endIndex, __shift, true);
startIndex = endIndex;
}
}


### Other things

I think that vector<value_type> is an appropriate container for your auxiliary array. I can't think of anything better than that.

As far as copy vs move, they should be equivalent for numeric types. See this Stackoverflow question and its answers for some good explanations on why that is.

• That is a very elegant approach to handle unnecessary passes on radix_sort_32. Appreciate your input on the matter. – Guilherme Barros Avila Mar 3 '17 at 0:15
• @GuilhermeBarrosAvila Please check out the section I added on the "hybrid" radix sort. – JS1 Mar 3 '17 at 10:06
• hmmm I knew about this approach, but I had no idea how or why it worked. I believe this segment of radix_sort_lsd: (__s + count[(*__p >> __shift) & 0xFF]++) = *__p; should have a * on the left side, right? – Guilherme Barros Avila Mar 3 '17 at 13:05
• Also, if you wanna comment on the use of vector<value_type> instead of another type of container and whether I should use std::copy or std::move to "copy" from the auxiliary vector to the input vector. I could definitely accept your answer as correct, given all the interesting suggestions you gave so far. – Guilherme Barros Avila Mar 3 '17 at 13:10
• @rcgldr The matrix of counts would probably help a little bit. But the counting pass is already somewhat fast because it goes through the input in order. The sorting pass is the big problem, because you are writing values to 256 possible cache lines. If you have a 4-way associative cache and 5+ of the 256 slots happen to occupy the same cache line, then one of those lines will get evicted. And that is what will cause the biggest slowdown. With the hybrid version, if each bin is smaller than the cache, you shouldn't get any cache evictions during the sorting pass. – JS1 Mar 4 '17 at 3:02

I'm wondering what level of optimization was used, as my time for radix_sort_32 averaged 22.8 ms versus 89.75 ms, about 4 times as fast. Converting from using iterators to pointers to first element of vector and then using pointer[index] reduced this to 14.8 ms.

As commented on earlier, in the case of 1 million 32 bit integers, both the original and working array will fit within the 8MB L3/L4 caches found on many processors. This greatly reduces the overhead of the random access writes during each radix sort pass. The vector / array size would need to be larger, beyond the bounds of processor cache, before the random access writes become a much more significant issue.

radix sort times for 1048576 32 bit unsgined integers
Intel 3770K, 3.5ghz, 32 bit mode (Windows XP), Visual Studio 2005

method                                       average time

radix_sort_32 (set msb = 32)                 0.0228 seconds
array                                        0.0132 seconds
array with matrix for counts/indices         0.0127 seconds
array with 4 threads and matrix              0.0115 seconds


As asked for by the original poster Guilherme Barros Avila, here is example code for multi-threaded radix sort. The multi-threaded part is Windows specific. RadixSort() is the same matrix oriented function previously posted (except that the working array b[] is passed as a parameter to avoid multiple allocations for the multiple threads) and not Windows specific. QP define is for Windows performance counter.

//------------------------------------------------------------------//
//------------------------------------------------------------------//
#include <cstdlib>
#include <ctime>
#include <iostream>
#include <windows.h>

// COUNT must be multiple of 4
#define COUNT (1024*1024)

typedef unsigned int uint32_t;

#define QP 0        // if != 0, use queryperformance for timer

#if QP
#include <math.h>
#include <windows.h>
#pragma comment(lib, "winmm.lib")
typedef LARGE_INTEGER LI64;
#endif

#if QP
static LI64     liQPFrequency;  // cpu counter values
static LI64     liStartTime;
static LI64     liStopTime;
static double   dQPFrequency;
static double   dStartTime;
static double   dStopTime;
static double   dElapsedTime;
#else
clock_t ctTimeStart;            // clock values
clock_t ctTimeStop;
#endif

static HANDLE hs0;                      // semaphore handles
static HANDLE hs1;
static HANDLE hs2;
static HANDLE hs3;
static HANDLE ht1;                      // thread handles
static HANDLE ht2;
static HANDLE ht3;

static uint32_t *pa;                    // pointers to buffers
static uint32_t *pb;

void RadixSort(uint32_t a[], uint32_t b[], size_t n);
void Merge(uint32_t a[], uint32_t b[], size_t ll, size_t rr, size_t ee);

uint32_t main()
{
uint32_t *array  = new uint32_t[COUNT];
uint32_t *buffer = new uint32_t[COUNT];
pa = array;
pb = buffer;
for(uint32_t i = 0; i < COUNT; i++){    // generate pseudo random data
uint32_t r;
r  = (((uint32_t)((rand()>>4) & 0xff))<< 0);
r += (((uint32_t)((rand()>>4) & 0xff))<< 8);
r += (((uint32_t)((rand()>>4) & 0xff))<<16);
r += (((uint32_t)((rand()>>4) & 0xff))<<24);
array[i] = r;
}

hs0 = CreateSemaphore(NULL,0,1,NULL);
hs1 = CreateSemaphore(NULL,0,1,NULL);
hs2 = CreateSemaphore(NULL,0,1,NULL);
hs3 = CreateSemaphore(NULL,0,1,NULL);

#if QP
QueryPerformanceFrequency(&liQPFrequency);
timeBeginPeriod(1);                     // set ticker to 1000 hz
Sleep(128);                             // wait for it to settle
QueryPerformanceCounter(&liStartTime);
#else
ctTimeStart = clock();
#endif

ReleaseSemaphore(hs0, 1, NULL);     // start sorts
ReleaseSemaphore(hs1, 1, NULL);
ReleaseSemaphore(hs2, 1, NULL);
ReleaseSemaphore(hs3, 1, NULL);
WaitForSingleObject(ht2, INFINITE); // wait for thead 2
// merge 1st and 2nd halves
Merge(pb, pa, 0, COUNT>>1, COUNT);

#if QP
QueryPerformanceCounter(&liStopTime);
dElapsedTime = (dStopTime - dStartTime) / dQPFrequency;
timeEndPeriod(1);                       // restore ticker to default
std::cout << "# of seconds " << dElapsedTime << std::endl;
#else
ctTimeStop = clock();
std::cout << "# of ticks " << ctTimeStop - ctTimeStart << std::endl;
#endif

for(uint32_t i = 1; i < COUNT; i++){        // check result
if(array[i-1] > array[i]){
std::cout << "failed" << std::endl;
break;
}
}
CloseHandle(ht3);
CloseHandle(ht2);
CloseHandle(ht1);
CloseHandle(hs3);
CloseHandle(hs2);
CloseHandle(hs1);
CloseHandle(hs0);
delete[] buffer;
delete[] array;
return 0;
}

{
WaitForSingleObject(hs0, INFINITE); // wait for semaphore
// sort 1st quarter
RadixSort(pa + 0*(COUNT>>2), pb + 0*(COUNT>>2), COUNT>>2);
WaitForSingleObject(ht1, INFINITE); // wait for thead 1
// merge 1st and 2nd quarter
Merge(pa + 0*(COUNT>>1), pb + 0*(COUNT>>1), 0, COUNT>>2, COUNT>>1);
return 0;
}

{
WaitForSingleObject(hs1, INFINITE); // wait for semaphore
// sort 2nd quarter
RadixSort(pa + 1*(COUNT>>2), pb + 1*(COUNT>>2), COUNT>>2);
return 0;
}

{
WaitForSingleObject(hs2, INFINITE); // wait for semaphore
// sort 3rd quarter
RadixSort(pa + 2*(COUNT>>2), pb + 2*(COUNT>>2), COUNT>>2);
WaitForSingleObject(ht3, INFINITE); // wait for thread 3
// merge 3rd and 4th quarter
Merge(pa + 1*(COUNT>>1), pb + 1*(COUNT>>1), 0, COUNT>>2, COUNT>>1);
return 0;
}

{
WaitForSingleObject(hs3, INFINITE); // wait for semaphore
// sort 4th quarter
RadixSort(pa + 3*(COUNT>>2), pb + 3*(COUNT>>2), COUNT>>2);
return 0;
}

void RadixSort(uint32_t a[], uint32_t b[], size_t count)
{
size_t mIndex[4][256] = {0};            // count / index matrix
size_t i,j,m,n;
uint32_t u;
for(i = 0; i < count; i++){         // generate histograms
u = a[i];
for(j = 0; j < 4; j++){
mIndex[j][(size_t)(u & 0xff)]++;
u >>= 8;
}
}
for(j = 0; j < 4; j++){             // convert to indices
m = 0;
for(i = 0; i < 256; i++){
n = mIndex[j][i];
mIndex[j][i] = m;
m += n;
}
}
for(j = 0; j < 4; j++){             // radix sort
for(i = 0; i < count; i++){     //  sort by current lsb
u = a[i];
m = (size_t)(u>>(j<<3))&0xff;
b[mIndex[j][m]++] = u;
}
std::swap(a, b);                //  swap ptrs
}
}

void Merge(uint32_t a[], uint32_t b[], size_t ll, size_t rr, size_t ee)
{
size_t o = ll;                      // b[]       index
size_t l = ll;                      // a[] left  index
size_t r = rr;                      // a[] right index
while(1){                           // merge data
if(a[l] <= a[r]){               // if a[l] <= a[r]
b[o++] = a[l++];            //   copy a[l]
if(l < rr)                  //   if not end of left run
continue;               //     continue (back to while)
do                          //   else copy rest of right run
b[o++] = a[r++];
while(r < ee);
break;                      //     and return
} else {                        // else a[l] > a[r]
b[o++] = a[r++];            //   copy a[r]
if(r < ee)                  //   if not end of right run
continue;               //     continue (back to while)
do                          //   else copy rest of left run
b[o++] = a[l++];
while(l < rr);
break;                      //     and return
}
}
}