RadixSort implementation design and performance

Question

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>
void radix_sort(Iter __first, Iter __last){
    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.

Answer 1

Improvement to radix_sort_32

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.

Hybrid radix sort

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>
void radix_sort(Iter __first, Iter __last)
{
    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) {
        radix_sort_lsd(__first, __last, out.begin(), out.end(),
                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.

Answer 2

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
radix_sort_32 using pointer[index]           0.0148 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.

//------------------------------------------------------------------//
//  mtrsrt.cpp      mult-threaded radix sort                        //
//------------------------------------------------------------------//
#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 DWORD WINAPI Thread0(LPVOID);    // thread functions
static DWORD WINAPI Thread1(LPVOID);
static DWORD WINAPI Thread2(LPVOID);
static DWORD WINAPI Thread3(LPVOID);

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);
    ht1 = CreateThread(NULL, 0, Thread1, 0, 0, 0);
    ht2 = CreateThread(NULL, 0, Thread2, 0, 0, 0);
    ht3 = CreateThread(NULL, 0, Thread3, 0, 0, 0);

#if QP
    QueryPerformanceFrequency(&liQPFrequency);
    dQPFrequency = (double)liQPFrequency.QuadPart;
    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);
    Thread0((LPVOID)NULL);              // run thread 0
    WaitForSingleObject(ht2, INFINITE); // wait for thead 2
    // merge 1st and 2nd halves
    Merge(pb, pa, 0, COUNT>>1, COUNT);

#if QP
    QueryPerformanceCounter(&liStopTime);
    dStartTime = (double)liStartTime.QuadPart;
    dStopTime  = (double)liStopTime.QuadPart;
    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;
}

static DWORD WINAPI Thread0(LPVOID lpvoid)
{
    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;
}

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

static DWORD WINAPI Thread2(LPVOID lpvoid)
{
    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;
}

static DWORD WINAPI Thread3(LPVOID lpvoid)
{
    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
        }
    }
}