MSD radix sort in C++

Question

I have this MSD radix sort in C++:

#ifndef RADIXSORT_H
#define RADIXSORT_H

#include <algorithm>
#include <type_traits>

static constexpr size_t BITS_PER_BYTE = 8;
static constexpr size_t BUCKET_AMOUNT = 256;
static constexpr size_t QUICKSORT_THRESHOLD = 256;

template<class T>
static size_t getBucketIndexUnsigned(T element, size_t byteIndex)
{
    return static_cast<size_t>((element >> (byteIndex * BITS_PER_BYTE)) & 0xff);
}

template<class RandomIt>
void unsigned_radix_sort(RandomIt firstSource, RandomIt lastSource,
                         RandomIt firstTarget, RandomIt lastTarget,
                         size_t byteIndex)
{
    auto rangeLength = std::distance(firstSource, lastSource);

    if (rangeLength < QUICKSORT_THRESHOLD)
    {
        std::sort(firstSource, lastSource);

        if ((byteIndex & 1) == 0)
        {
            std::copy(firstSource, lastSource, firstTarget);
        }

        return;
    }

    typedef typename std::iterator_traits<RandomIt>::value_type value_type;

    size_t bucketSizeMap[BUCKET_AMOUNT];
    size_t startIndexMap[BUCKET_AMOUNT];
    size_t processedMap [BUCKET_AMOUNT];

    for (size_t i = 0; i < BUCKET_AMOUNT; ++i)
    {
        bucketSizeMap[i] = 0;
        startIndexMap[i] = 0;
        processedMap[i]  = 0;
    }

    for (RandomIt it = firstSource; it != lastSource; ++it)
    {
        bucketSizeMap[getBucketIndexUnsigned(*it, byteIndex)]++;
    }

    for (size_t i = 1; i < BUCKET_AMOUNT; ++i)
    {
        startIndexMap[i] = startIndexMap[i - 1] + bucketSizeMap[i - 1];
    }

    for (RandomIt it = firstSource; it != lastSource; ++it)
    {
        value_type element = *it;
        size_t bucket = getBucketIndexUnsigned(element, byteIndex);
        *(firstTarget + startIndexMap[bucket]
          + processedMap[bucket]++) = element;
    }

    if (byteIndex > 0)
    {
        for (size_t i = 0; i < BUCKET_AMOUNT; ++i)
        {
            if (bucketSizeMap[i] != 0)
            {
                unsigned_radix_sort(firstTarget + startIndexMap[i],
                                    firstTarget + startIndexMap[i]
                                    + bucketSizeMap[i],
                                    firstSource + startIndexMap[i],
                                    firstSource + startIndexMap[i]
                                    + bucketSizeMap[i],
                                    byteIndex - 1);
            }
        }
    }
}

template<class RandomIt>
void signed_radix_sort(RandomIt firstSource, RandomIt lastSource,
                       RandomIt firstTarget, RandomIt lastTarget,
                       size_t byteIndex)
{
    auto rangeLength = std::distance(firstSource, lastSource);

    if (rangeLength < QUICKSORT_THRESHOLD)
    {
        std::sort(firstSource, lastSource);

        if ((byteIndex & 1) == 0)
        {
            std::copy(firstSource, lastSource, firstTarget);
        }

        return;
    }

    typedef typename std::iterator_traits<RandomIt>::value_type value_type;

    size_t bucketSizeMap[BUCKET_AMOUNT];
    size_t startIndexMap[BUCKET_AMOUNT];
    size_t processedMap [BUCKET_AMOUNT];

    for (size_t i = 0; i < BUCKET_AMOUNT; ++i)
    {
        bucketSizeMap[i] = 0;
        startIndexMap[i] = 0;
        processedMap[i]  = 0;
    }

    for (RandomIt it = firstSource; it != lastSource; ++it)
    {
        bucketSizeMap[getBucketIndexUnsigned(*it, byteIndex)]++;
    }

    startIndexMap[BUCKET_AMOUNT >> 1] = 0;

    for (size_t i = (BUCKET_AMOUNT >> 1) + 1; i < BUCKET_AMOUNT; ++i)
    {
        startIndexMap[i] = startIndexMap[i - 1] + bucketSizeMap[i - 1];
    }

    startIndexMap[0] = startIndexMap[BUCKET_AMOUNT - 1] +
                       bucketSizeMap[BUCKET_AMOUNT - 1];

    for (size_t i = 1; i < BUCKET_AMOUNT >> 1; ++i)
    {
        startIndexMap[i] = startIndexMap[i - 1] + bucketSizeMap[i - 1];
    }

    for (RandomIt it = firstSource; it != lastSource; ++it)
    {
        value_type element = *it;
        size_t bucket = getBucketIndexUnsigned(element, byteIndex);
        *(firstTarget + startIndexMap[bucket]
          + processedMap[bucket]++) = element;
    }

    if (byteIndex > 0)
    {
        for (size_t i = 0; i < BUCKET_AMOUNT; ++i)
        {
            if (bucketSizeMap[i] != 0)
            {
                unsigned_radix_sort(firstTarget + startIndexMap[i],
                                    firstTarget + startIndexMap[i]
                                    + bucketSizeMap[i],
                                    firstSource + startIndexMap[i],
                                    firstSource + startIndexMap[i]
                                    + bucketSizeMap[i],
                                    byteIndex - 1);
            }
        }
    }
}

template<class RandomIt>
void unsigned_radix_sort(RandomIt first, RandomIt last)
{
    auto rangeLength = std::distance(first, last);
    typedef typename std::iterator_traits<RandomIt>::value_type value_type;
    value_type* aux = new value_type[rangeLength];

    unsigned_radix_sort(first,
                        last,
                        aux,
                        aux + rangeLength,
                        sizeof(value_type) - 1);
    delete[] aux;
}

template<class RandomIt>
void signed_radix_sort(RandomIt first, RandomIt last)
{
    auto rangeLength = std::distance(first, last);
    typedef typename std::iterator_traits<RandomIt>::value_type value_type;
    value_type* aux = new value_type[rangeLength];

    signed_radix_sort(first,
                      last,
                      aux,
                      aux + rangeLength,
                      sizeof(value_type) - 1);
    delete[] aux;
}

template<class RandomIt>
void sort(RandomIt first, RandomIt last)
{
    if (std::distance(first, last) < 2)
    {
        // Nothing to sort.
        return;
    }

    typedef typename std::iterator_traits<RandomIt>::value_type value_type;

    if (std::is_unsigned<value_type>::value)
    {
        // The input integers are unsigned.
        unsigned_radix_sort(first, last);
    }
    else if (std::is_signed<value_type>::value)
    {
        // The input integers are signed.
        signed_radix_sort(first, last);
    }
    else
    {
        // The input objects not integers at all; delegate to std::sort.
        std::sort(first, last);
    }
}

#endif // RADIXSORT_H

Performance figures

Looks like this:

--- Unsigned array ---
std::sort() in 6992 milliseconds.
Radix sort in 3494 milliseconds.
Equal: true
--- Signed array ---
std::sort() in 6852 milliseconds.
Radix sort in 3445 milliseconds.
Equal: true
--- Unsigned vector ---
std::sort() in 7028 milliseconds.
Radix sort in 7580 milliseconds.
Equal: true
--- Signed vector ---
std::sort() in 6975 milliseconds.
Radix sort in 7050 milliseconds.
Equal: true
Bye!

Critique request

I have only two questions:

1. How can I tell the compiler that my radix sorts accepts only integer types?
2. Is there any possibility of improving performance on std::vectors?

Please, tell me anything that comes to mind.

Answer 1

Use System costants as appropriate

This is already defined by the Library:

static constexpr size_t BITS_PER_BYTE = 8;

The number of bits in a byte is: CHAR_BITS.

Static

Static is one of those overloaded keywords. Using it on a function in a header file (that has a definition). Means you get one definition in every compilation unit (that includes the header file). This is probably not what you want.

template<class T>
static size_t getBucketIndexUnsigned(T element, size_t byteIndex)
{
    return static_cast<size_t>((element >> (byteIndex * BITS_PER_BYTE)) & 0xff);
}

DRY code

The implementation of unsigned_radix_sort and signed_radix_sort is nearly identical. Why not make the code identical? I am sure with only minor tweaks you can combine the code bases of these two functions.

Iterator types.

There are to ranges here:

template<class RandomIt>
void unsigned_radix_sort(RandomIt firstSource, RandomIt lastSource,
                         RandomIt firstTarget, RandomIt lastTarget,
                         size_t byteIndex)

Each controlled by a different container type. So I would guess that you probably want to use two different kids of iterator.

template<typename RandomIt1, typename RandomIt2>
void unsigned_radix_sort(RandomIt1 firstSource, RandomIt1 lastSource,
                         RandomIt2 firstTarget, RandomIt2 lastTarget,
                         size_t byteIndex)

Comment

        if ((byteIndex & 1) == 0)
        {
            std::copy(firstSource, lastSource, firstTarget);
        }

This needs a comment. It took me a while to figure out that you are only doing this on odd recursive calls because you swap data from one set of container to the other.

Using.

    typedef typename std::iterator_traits<RandomIt>::value_type value_type;

I have started using the new using style to declare type synonyms.

    using value_type = typename std::iterator_traits<RandomIt>::value_type;

Initialize to zero automatically

You can declare and initialize thes buckets using one line (no need for the following loop).

    size_t bucketSizeMap[BUCKET_AMOUNT];
    size_t startIndexMap[BUCKET_AMOUNT];
    size_t processedMap [BUCKET_AMOUNT];

    for (size_t i = 0; i < BUCKET_AMOUNT; ++i)
    {
        bucketSizeMap[i] = 0;
        startIndexMap[i] = 0;
        processedMap[i]  = 0;
    }

Just do:

    size_t bucketSizeMap[BUCKET_AMOUNT]();    // Zero initialize values.
    size_t startIndexMap[BUCKET_AMOUNT] = {0};// or zero init C style.
    size_t processedMap [BUCKET_AMOUNT]();

Algorithms:

You should be looking at using a standard algorithm with a lambda (or at least the new range based for).

    for (RandomIt it = firstSource; it != lastSource; ++it)
    {
        bucketSizeMap[getBucketIndexUnsigned(*it, byteIndex)]++;
    }

    for (size_t i = 1; i < BUCKET_AMOUNT; ++i)
    {
        startIndexMap[i] = startIndexMap[i - 1] + bucketSizeMap[i - 1];
    }

Using standard algorithms.

    std::for_each(firstSource, lastSource, [](value_type const& v)
                  {++bucketSizeMap[getBucketIndexUnsigned(v, byteIndex)];});

    std::transform(startIndexMap, startIndexMap + BUCKET_AMOUNT - 1, bucketSizeMap,
                   startIndexMap + 1, [](value_type const& v1, value_type const& v2) { return v1 + v2;});

Again. You could probably do the next loop using a standard algorithm. But this time I am not 100% convinced that is a good idea as it is slightly longer. But you may want to look into it.

    for (RandomIt it = firstSource; it != lastSource; ++it)
    {
        value_type element = *it;
        size_t bucket = getBucketIndexUnsigned(element, byteIndex);
        *(firstTarget + startIndexMap[bucket]
          + processedMap[bucket]++) = element;
    }

RAII

This is a classic anti pattern for resource allocation.

    value_type* aux = new value_type[rangeLength];

    // STUFF

    delete[] aux;

Technically not exception safe. This may become more important if you extend this past integer types. But it is just as easy to do the following instead.

    std::vector<value_type>   aux(rangeLength);

Compile Time Vs Runtime Time

This test is a classic example of when to perform compile time test and plant the correct code using templates rather than a run-time test.

    if (std::is_unsigned<value_type>::value)
    {
        // The input integers are unsigned.
        unsigned_radix_sort(first, last);
    }
    else if (std::is_signed<value_type>::value)
    {
        // The input integers are signed.
        signed_radix_sort(first, last);
    }
    else
    {
        // The input objects not integers at all; delegate to std::sort.
        std::sort(first, last);
    }

I would do:

    template<typename I, 
             typename T = typename std::iterator_traits<RandomIt>::value_type,
             bool isSigned = std::is_signed<T>::value,
             bool isUnsigned = std::is_unsigned<T>::value>
    void doSort(I begin, I end)
    {
        std::sort(begin, end);
    }

    template<typename I, typename T>
    void doSort<I, T, true, false>(first, last)
    {
        signed_radix_sort(first, last);
    }

    template<typename I, typename T>
    void doSort<I, T, false, true>(first, last)
    {
        unsigned_radix_sort(first, last);
    }

template<class RandomIt>
void sort(RandomIt first, RandomIt last)
{
    if (std::distance(first, last) < 2)
    {
        // Nothing to sort.
        return;
    }
    std::doSort(first, last);
}