Parallel natural merge sort

Question

Continuing working on Natural merge sort, I have parallelized it. Requires \$\Theta(N)\$ space and runs in $$\Omega(N + \frac{N}{P}) \cap \mathcal{O}(N + \frac{N}{P}\log_2\frac{N}{P})$$ time, where \$P\$ is the amount of cores available.

parallel_natural_merge_sort.h:

#ifndef NATURAL_MERGE_SORT_H
#define NATURAL_MERGE_SORT_H

#include <algorithm>
#include <iterator>
#include <thread>

/*******************************************************************************
* Implements a simple, array-based queue of integers. All three operations run *
* in constant time. This queue, however, does not check for under-/overflow of *
* underlying buffer because of performance considerations.                     *
*******************************************************************************/
class UnsafeIntQueue {
private:
    const size_t MINIMUM_CAPACITY = 256;

    size_t m_head;
    size_t m_tail;
    size_t m_size;
    size_t m_mask;
    size_t* m_buffer;

    /***************************************************************************
    * Makes sure a capacity is at least 'MINIMUM_CAPACITY' and is a power of   *
    * two.                                                                     *
    ***************************************************************************/
    size_t fixCapacity(size_t capacity)
    {
        capacity = std::max(capacity, MINIMUM_CAPACITY);
        size_t s = 1;

        while (s < capacity)
        {
            s <<= 1;
        }

        return s;
    }

public:

    /***************************************************************************
    * Constructs a new integer queue, which can accommodate 'capacit' amount   *
    * integers.                                                                *
    ***************************************************************************/
    UnsafeIntQueue(size_t capacity) :
    m_head{0},
    m_tail{0},
    m_size{0}
    {
        capacity = fixCapacity(capacity);
        m_mask = capacity - 1;
        m_buffer = new size_t[capacity];
    }

    /***************************************************************************
    * Destroys this queue, which releases the underlying buffer.               *
    ***************************************************************************/
    ~UnsafeIntQueue()
    {
        delete[] m_buffer;
    }

    /***************************************************************************
    * Appends the input integer to the tail of this queue.                     *
    ***************************************************************************/
    inline void enqueue(const size_t element)
    {
        m_buffer[m_tail & m_mask] = element;
        m_tail = (m_tail + 1) & m_mask;
        m_size++;
    }

    /***************************************************************************
    * Removes and returns the integer at the head of this queue.               *
    ***************************************************************************/
    inline size_t dequeue()
    {
        const size_t ret = m_buffer[m_head];
        m_head = (m_head + 1) & m_mask;
        m_size--;
        return ret;
    }

    /***************************************************************************
    * Returns the amount of integers in this queue.                            *
    ***************************************************************************/
    inline size_t size() const
    {
        return m_size;
    }
};

/*******************************************************************************
* Scans the range [first, last) and returns the queue containing sizes of each *
* run in the order they appear while scanning from left to right.              *
*******************************************************************************/
template<class RandomIt, class Cmp>
std::unique_ptr<UnsafeIntQueue> build_run_size_queue(RandomIt first,
                                                     RandomIt last,
                                                     Cmp cmp)
{
    const size_t length = std::distance(first, last);
    UnsafeIntQueue* p_q = new UnsafeIntQueue(length / 2 + 1);

    RandomIt head;
    RandomIt left = first;
    RandomIt right = left + 1;

    const RandomIt lst = last - 1;

    while (left < lst)
    {
        head = left;

        if (cmp(*right++, *left++))
        {
            // Reading a strictly descending run.
            while (left < lst && cmp(*right, *left))
            {
                ++left;
                ++right;
            }

            p_q->enqueue(right - head);
            std::reverse(head, right);
        }
        else
        {
            // Reading a ascending run.
            while (left < lst && !cmp(*right, *left))
            {
                ++left;
                ++right;
            }

            p_q->enqueue(left - head + 1);
        }

        ++left;
        ++right;
    }

    if (left == lst)
    {
        // Handle the case of an orphan element at the end of the range.
        p_q->enqueue(1);
    }

    return std::unique_ptr<UnsafeIntQueue>(p_q);
}

/*******************************************************************************
* Returns the amount of leading zeros in 'num'.                                *
*******************************************************************************/
size_t leading_zeros(const size_t num)
{
    size_t count = 0;

    for (size_t t = (size_t) 1 << (8 * sizeof(t) - 1); t; t >>= 1, ++count)
    {
        if ((t & num))
        {
            return count;
        }
    }

    return count;
}

/*******************************************************************************
* Returns the amount of merge passes needed to sort a range with 'run_amount'  *
* runs.                                                                        *
*******************************************************************************/
size_t get_pass_amount(size_t run_amount)
{
    return 8 * sizeof(run_amount) - leading_zeros(run_amount - 1);
}

/*******************************************************************************
* The actual implementation of natural merge sort.                             *
*******************************************************************************/
template<class RandomIt, class Cmp>
void natural_merge_sort_impl(RandomIt first, 
                             RandomIt last, 
                             RandomIt buffer,
                             Cmp cmp)
{
    const size_t length = std::distance(first, last);

    if (length < 2)
    {
        // Trivially sorted.
        return;
    }

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

    // Scan the runs.
    std::unique_ptr<UnsafeIntQueue> p_queue = build_run_size_queue(first, last, cmp);

    // Count the amount of merge passes over the array required to bring order.
    const size_t merge_passes = get_pass_amount(p_queue->size());

    RandomIt source;
    RandomIt target;

    // Make sure that after the last merge pass, all data ends up in the input
    // container.
    if ((merge_passes & 1) == 1)
    {
        source = buffer;
        target = first;
        std::copy(first, last, buffer);
    }
    else
    {
        source = first;
        target = buffer;
    }

    size_t runs_left = p_queue->size();
    size_t offset = 0;

    // While there is runs to merge, do...
    while (p_queue->size() > 1)
    {
        // Remove two runs from the head of the run queue.
        size_t left_run_length = p_queue->dequeue();
        size_t right_run_length = p_queue->dequeue();

        std::merge(source + offset,
                   source + offset + left_run_length,
                   source + offset + left_run_length,
                   source + offset + left_run_length + right_run_length,
                   target + offset,
                   cmp);

        // Append the merged run to the tail of the queue.
        p_queue->enqueue(left_run_length + right_run_length);
        runs_left -= 2;
        offset += left_run_length + right_run_length;

        // The current pass over the array is almost complete.
        switch (runs_left)
        {
            case 1:
            {
                const size_t single_length = p_queue->dequeue();

                std::copy(source + offset,
                          source + offset + single_length,
                          target + offset);

                p_queue->enqueue(single_length);
            }

            // FALL THROUGH!

            case 0:
            {
                runs_left = p_queue->size();
                offset = 0;
                RandomIt tmp = source;
                source = target;
                target = tmp;
                break;
            }
        }
    }
}

/*******************************************************************************
* Implements the natural merge sort, which sacrifices one pass over the input  *
* range in order to establish an implicit queue of runs. A run is the longest  *
* consecutive subsequence, in which all elements are ascending or strictly     *
* descending. Every descending run is reversed to ascending run. We cannot     *
* consider non-strictly descending runs, since that would sacrifice the stabi- *
* lity of the algorithm. After the run queue is establish, the algorithm re-   *
* moves two runs from the head of the queue, merges them into one run, which   *
* is then appended to the tail of the run queue. Merging continues until the   *
* queue contains only one run, which denotes that the entire input range is    *
* sorted.                                                                      *
*                                                                              *
* The best-case complexity is O(N), the average and worst-case complexity is   *
* O(N log N). Space complexity is O(N).                                        *
*******************************************************************************/
template<class RandomIt, class Cmp>
void natural_merge_sort(RandomIt first, RandomIt last, Cmp cmp)
{
    const size_t length = std::distance(first, last);

    if (length < 2)
    {
        // Trivially sorted.
        return;
    }

    typedef typename std::iterator_traits<RandomIt>::value_type value_type;
    RandomIt buffer = new value_type[length];
    natural_merge_sort_impl(first, last, buffer, cmp);
    delete[] buffer;
}

/*******************************************************************************
* Implements parallel merge sort.                                              *
*******************************************************************************/
template<class RandomIt, class Cmp>
void parallel_natural_merge_sort_impl(RandomIt source, 
                                      RandomIt target, 
                                      const size_t length, 
                                      const size_t thread_quota,
                                      Cmp cmp)
{
    if (thread_quota == 1)
    {
        natural_merge_sort_impl(target, target + length, source, cmp);
        return;
    }

    const size_t left_quota = thread_quota / 2;
    const size_t right_quota = thread_quota - left_quota;
    const size_t left_length = length / 2;

    if (thread_quota == 2)
    {
        std::thread thread_(natural_merge_sort_impl<RandomIt, Cmp>,
                            source,
                            source + left_length,
                            target,
                            cmp);

        natural_merge_sort_impl(source + left_length, 
                                source + length, 
                                target + left_length, 
                                cmp);

        thread_.join();

        std::merge(source, 
                   source + left_length, 
                   source + left_length, 
                   source + length,
                   target,
                   cmp);
        return;
    }

    std::thread left_thread(parallel_natural_merge_sort_impl<RandomIt, Cmp>,
                            target, 
                            source, 
                            left_length, 
                            left_quota, 
                            cmp);

    parallel_natural_merge_sort_impl(target + left_length, 
                                     source + left_length, 
                                     length - left_length, 
                                     right_quota, 
                                     cmp);
    // Wait for the left thread.
    left_thread.join();

    // Merge the two chunks.
    std::merge(source, 
               source + left_length, 
               source + left_length, 
               source + length,
               target,
               cmp);
}

/*******************************************************************************
* The actual parallel merge sort. If the system has N CPU cores, this sort     *     
* will split the range into N chunks of equal length assuming that N is a      *
* power of two, sort them concurrently and merge.                              *
*******************************************************************************/
template<class RandomIt, class Cmp>
void parallel_natural_merge_sort(RandomIt begin, RandomIt end, Cmp cmp)
{
    // At least 16384 elements per thread.
    constexpr size_t MINIMUM_THREAD_LOAD = 1 << 14;
    const size_t cores = std::thread::hardware_concurrency();
    const size_t length = std::distance(begin, end);
    const size_t spawn = std::min(cores, length / MINIMUM_THREAD_LOAD);

    if (spawn < 2)
    {
        natural_merge_sort(begin, end, cmp);
        return;
    }

    typedef typename std::iterator_traits<RandomIt>::value_type value_type;
    RandomIt buffer = new value_type[length];
    std::copy(begin, end, buffer);
    parallel_natural_merge_sort_impl(buffer, begin, length, spawn, cmp);
}

#endif

main.cpp:

#include <chrono>
#include <functional>
#include <iostream>
#include <random>

#include "parallel_natural_merge_sort.h"

/*******************************************************************************
* Creates a random integer array of length 'length', minimum integer           *
* 'minimum', maximum integer 'maximum', using seed 'seed'.                     *
*******************************************************************************/
static int* get_random_int_array(const size_t length,
                                 const int minimum,
                                 const int maximum,
                                 const unsigned int seed)
{
    std::default_random_engine generator(seed);
    std::uniform_int_distribution<int> distribution(minimum, maximum);

    int* array = new int[length];

    for (size_t i = 0; i < length; ++i)
    {
        array[i] = distribution(generator);
    }

    return array;
}

/*******************************************************************************
* Create an array of pointers to integers.                                     *
*******************************************************************************/
static int** get_random_int_pointer_array(const size_t length,
                                          const int minimum,
                                          const int maximum,
                                          const unsigned seed)
{
    std::default_random_engine generator(seed);
    std::uniform_int_distribution<int> distribution(minimum, maximum);

    int** array = new int*[length];

    for (size_t i = 0; i < length; ++i)
    {
        array[i] = new int(distribution(generator));
    }

    return array;
}

/*******************************************************************************
* Returns a strongly presorted array of integers.                              *
*******************************************************************************/
static int* get_presorted_int_array(const size_t length)
{
    int* array = new int[length];
    int num = 0;

    for (size_t i = 0; i < length / 2; ++i)
    {
        array[i] = num++;
    }

    for (size_t i = length / 2; i < length; ++i)
    {
        array[i] = num--;
    }

    return array;
}

/*******************************************************************************
* Returns the milliseconds since the Unix epoch.                               *
*******************************************************************************/
static unsigned long long get_milliseconds()
{
    return std::chrono::duration_cast<std::chrono::milliseconds>(
           std::chrono::system_clock::now().time_since_epoch()).count();
}

/*******************************************************************************
* Profiles the 'std::stable_sort' agains the range ['begin', 'end') using the  *
* comparator 'cmp'.                                                            *
*******************************************************************************/
template<class T, class Cmp>
static void profile_stable_sort(T begin, T end, Cmp cmp)
{
    unsigned long long ta = get_milliseconds();
    std::stable_sort(begin, end, cmp);
    unsigned long long tb = get_milliseconds();

    std::cout << "std::stable_sort in "
              << (tb - ta)
              << " milliseconds. "
              << "Sorted: "
              << std::is_sorted(begin, end, cmp)
              << std::endl;
}

/*******************************************************************************
* Profiles the 'natural_merge_sort' agains the range ['begin', 'end') using    *
* the comparator 'cmp'.                                                        *
*******************************************************************************/
template<class T, class Cmp>
void profile_natural_merge_sort(T begin, T end, Cmp cmp)
{
    unsigned long long ta = get_milliseconds();
    natural_merge_sort(begin, end, cmp);
    unsigned long long tb = get_milliseconds();

    std::cout << "natural_merge_sort in "
              << (tb - ta)
              << " milliseconds. "
              << "Sorted: "
              << std::is_sorted(begin, end, cmp)
              << std::endl;
}

/*******************************************************************************
* Profiles the 'natural_merge_sort' agains the range ['begin', 'end') using    *
* the comparator 'cmp'.                                                        *
*******************************************************************************/
template<class T, class Cmp>
void profile_parallel_natural_merge_sort(T begin, T end, Cmp cmp)
{
    unsigned long long ta = get_milliseconds();
    parallel_natural_merge_sort(begin, end, cmp);
    unsigned long long tb = get_milliseconds();

    std::cout << "parallel_natural_merge_sort in "
              << (tb - ta)
              << " milliseconds. "
              << "Sorted: "
              << std::is_sorted(begin, end, cmp)
              << std::endl;
}

/*******************************************************************************
* Profiles the sorting algorithms on a random integer array.                   *
*******************************************************************************/
static void profile_on_random_array(const size_t sz,
                                    const int minimum,
                                    const int maximum,
                                    const unsigned seed)
{
    int* array1 = get_random_int_array(sz, minimum, maximum, seed);
    int* array2 = new int[sz];
    int* array3 = new int[sz];

    std::copy(array1, array1 + sz, array2);
    std::copy(array1, array1 + sz, array3);

    std::cout << "--- PROFILING ON RANDOM ARRAY OF LENGTH "
              << sz
              << " ---"
              << std::endl;

    profile_stable_sort(array1, 
                        array1 + sz, 
                        std::less<>());

    profile_natural_merge_sort(array2, 
                               array2 + sz, 
                               std::less<>());

    profile_parallel_natural_merge_sort(array3,
                                        array3 + sz,
                                        std::less<>());

    std::cout << "Same contents: "
              << (std::equal(array1, array1 + sz, array2, array2 + sz)
                 && std::equal(array1, array1 + sz, array3, array3 + sz))
              << std::endl
              << std::endl;
}

/*******************************************************************************
* Profiles the sorting algorithms on an array of pointers to random integers.  *
*******************************************************************************/
static void profile_on_integer_pointer_array(const size_t sz,
                                             const int minimum,
                                             const int maximum,
                                             const unsigned seed)
{
    std::cout << "--- PROFILING ON RANDOM POINTER ARRAY OF LENGTH "
              << sz
              << " ---"
              << std::endl;

    int** array1 = get_random_int_pointer_array(sz,
                                                minimum,
                                                maximum,
                                                seed);
    int** array2 = new int*[sz];
    int** array3 = new int*[sz];

    std::copy(array1, array1 + sz, array2);
    std::copy(array1, array1 + sz, array3);

    auto lambda = [](int* a, int* b){
        return *a < *b;
    };

    profile_stable_sort(array1, 
                        array1 + sz, 
                        lambda);

    profile_natural_merge_sort(array2, 
                               array2 + sz, 
                               lambda);

    profile_parallel_natural_merge_sort(array3,
                                        array3 + sz,
                                        lambda);
    std::cout << "Same contents: "
              << (std::equal(array1, array1 + sz, array2, array2 + sz)
                 && std::equal(array1, array1 + sz, array3, array3 + sz))
              << std::endl
              << std::endl;
}

/*******************************************************************************
* Profiles the sorting algorithms on a presorted array.                        *
*******************************************************************************/
static void profile_on_presorted_array(const size_t sz)
{
    std::cout << "--- PROFILING ON PRESORTED ARRAY OF LENGTH "
              << sz
              << " ---"
              << std::endl;

    int* array1 = get_presorted_int_array(sz);
    int* array2 = new int[sz];
    int* array3 = new int[sz];

    std::copy(array1, array1 + sz, array2);
    std::copy(array1, array1 + sz, array3);

    profile_stable_sort(array1, 
                        array1 + sz, 
                        std::less<>());

    profile_natural_merge_sort(array2, 
                               array2 + sz, 
                               std::less<>());

    profile_parallel_natural_merge_sort(array3,
                                        array3 + sz,
                                        std::less<>());

    std::cout << "Same contents: "
              << (std::equal(array1, array1 + sz, array2, array2 + sz) 
                 && std::equal(array1, array1 + sz, array3, array3 + sz))
              << std::endl
              << std::endl;
}

/*******************************************************************************
* The entry point to a demo program.                                           *
*******************************************************************************/
int main(int argc, const char * argv[]) {
    unsigned long long seed = get_milliseconds();

    std::cout << "Seed: "
              << seed
              << std::endl
              << std::endl;

    const size_t length = 5000000;
    const int min_int = -100;
    const int max_int = 300;

    std::cout << std::boolalpha;

    profile_on_random_array(length, min_int, max_int, seed);
    profile_on_integer_pointer_array(length, min_int, max_int, seed);
    profile_on_presorted_array(length);

    return 0;
}

Is there room for improvement? Efficiency? Style? Naming?

Answer 1

The style looks pretty idiomatic to me. Comments:

---

UnsafeIntQueue* p_q = new UnsafeIntQueue(length / 2 + 1);
... cmp(*right, *left) ...
return std::unique_ptr<UnsafeIntQueue>(p_q);

If the user-provided comparator cmp throws an exception, you leak memory. Prefer to wrap the heap-allocated object in a unique_ptr right away, so that if an exception is thrown it will be destroyed properly.

auto p_q = std::make_unique<UnsafeIntQueue>(length / 2 + 1);
... cmp(*right, *left) ...
return p_q;

---

In fact, build_run_size_queue should simply return a value of type UnsafeIntQueue — you don't need to muck around with heap-allocation here at all. Get rid of the unique_ptr and just return a simple UnsafeIntQueue directly.

---

Which brings up the fact that UnsafeIntQueue is not safe for copying or moving. (That's probably why you wrapped it in a unique_ptr: because you couldn't get it to work as a value type.) The solution to that is to get rid of that manual memory management; replace size_t *m_buffer with std::vector<size_t> m_buffer, and replace all uses of m_mask with m_buffer.size(). Don't worry about a loss of efficiency: std::vector works exactly like a new'ed array under the hood, as long as you stick to indexing with m_buffer[i] instead of m_buffer.at(i).

---

size_t leading_zeros(const size_t num)
{
    size_t count = 0;

    for (size_t t = (size_t) 1 << (8 * sizeof(t) - 1); t; t >>= 1, ++count)
    {
        if ((t & num))
        {
            return count;
        }
    }

    return count;
}

I would write this as

// Use __builtin_clzll(x) if available; otherwise...
size_t clz(size_t x)
{
    size_t result = CHAR_BIT * sizeof(x);
    while (x != 0) {
        x >>= 1;
        result -= 1;
    }
    return result;
}

Four lines shorter, including the comment (although it will get longer if you actually add the #ifdefs to use the non-standard compiler builtin). The main savings is due to brace style. Putting in the explicit != 0 is a good habit because it will silence warnings from some compilers. (Some compilers these days have blacklisted "assignment-used-as-an-expression" to such an extent that they actually warn about if ((expr)) — you put two pairs of parentheses where only one was needed, so clearly you've made some kind of mistake!)

Parameters and return values in C++11-or-later should never be const-qualified; it inhibits move semantics in general, and is thus a bad habit that you should break ASAP.

Note that std::numeric_limits<unsigned char>::digits is a valid synonym for CHAR_BIT, and in fact so is 8 in my opinion; but then in my opinion 64 is a valid synonym for 8 * sizeof(x), and you're clearly not willing to make that assumption yet.

---

In parallel_natural_merge_sort:

RandomIt buffer = new value_type[length];

This doesn't even compile, unless RandomIt happens to be value_type*. Try it with std::vector<value_type>::iterator or the like, and see what happens.

Furthermore, if it does compile, you get a memory leak because nobody ever frees the new'ed buffer. Replace your manual memory management with std::vector<value_type> buffer(length);, and change parallel_natural_merge_sort_impl to take two template type parameters RandomIt1 and RandomIt2, which are not necessarily synonymous.

---

Being able to call new value_type[length] (or std::vector<value_type>(length)) depends on value_type being default-constructible, which might not be the case: it might be moveable but not default-constructible nor even copyable. Does the "natural merge sort" algorithm really require that the type be copyable, or could you refactor this code to work with moves alone?

---

/* Scans the range [first, last) ... */

The usual names in C++ for the endpoints of a half-open range are begin and end. I don't see a good reason not to use those names here. This would free up the name last (or, idiomatically, back) for the pointer-to-last-item variable which is currently confusingly named lst.

---

const size_t length = std::distance(first, last);

The actual return type here is std::iterator_traits<RandomIt>::difference_type, which is most likely signed ptrdiff_t but could be anything under the sun (and is not necessarily convertible to size_t). This is a great place to use auto to make your code more generic. I haven't traced all the way down the rabbit hole, but I'm pretty sure that ultimately you'll want to make UnsafeIntQueue into a class template so that it can store difference_types instead of merely size_ts.

---

Inside natural_merge_sort_impl:

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

This is scrupulously correct, but (A) in C++11-and-later it's idiomatic to use the new syntax

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

and (B) you don't actually use this typedef anywhere in the function.

---

In parallel_natural_merge_sort_impl:

if (thread_quota == 2)

You can kill this whole block, can't you? The code immediately below the if block appears to do exactly the same thing (spawn two threads, join them, merge the resulting arrays).

---

In build_run_size_queue, if it's permitted by the rules of the game, consider simply invoking std::sort on the buffer (or maybe on medium-sized chunks of it). Right now you're wasting a lot of time comparing elements in order to come up with run-lengths that will generally be only 2 or 3, and then you're going to waste time again comparing those elements to each other in order to merge them. You could instead do something like "Count up the run-length starting at p; if you find that the run-length is less than 100, then std::sort the 100 elements starting at p and then continue counting from p+100" (with suitable bounds-checking, of course).

---

As usual, no comment on the correctness of the algorithm. However, I find it hard to believe that any sort algorithm with so many std::copys on large ranges could outperform an algorithm that uses only move and swap.