# Natural merge sort

Note: the next iteration at Natural merge sort - follow-up

I have compiled natural merge sort algorithm from Java to C++.

Natural merge sort sacrifices at most $2N$ amount of work in order to recognize the runs in the input range. A run is a contiguous subsequence which is ascending or strictly descending. Every descending run is reversed. We don't allow non-strictly descending runs, since reversing them would rearrange possible equal elements in it and make the entire algorithm non-stable.

After it has figured the runs, it puts them all into a queue. Then while there is runs to merge, it pops two runs, merges them, and appends the new merged run to the tail of the queue. It continues in this manner until there is only one run in the queue left, which corresponds to the sorted range.

So, all in all, the time complexity of this sort is $\Omega(N) \cap \mathcal{O}(N \log N)$. $\Theta(N)$ space complexity.

Is that a good way of writing (potentially) reusable code? Are there some C++-idioms I should have adhered to?

natural_merge_sort.h:

#ifndef NATURAL_MERGE_SORT_H
#define NATURAL_MERGE_SORT_H

#include <algorithm>
#include <iostream>
#include <iterator>

/*******************************************************************************
* 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_tail;
size_t m_size;
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_tail{0},
m_size{0}
{
capacity = fixCapacity(capacity);
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.                     *
***************************************************************************/
void enqueue(const size_t element)
{
m_tail = (m_tail + 1) & m_mask;
m_size++;
}

/***************************************************************************
* Removes and returns the integer at the head of this queue.               *
***************************************************************************/
size_t dequeue()
{
m_size--;
return ret;
}

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

/*******************************************************************************
* Scans the range [first, lst) 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 lst,
Cmp cmp)
{
const size_t length = std::distance(first, lst);
UnsafeIntQueue* p_q = new UnsafeIntQueue(length / 2 + 1);

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

const RandomIt last = lst - 1;

while (left < last)
{

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

}
else
{
while (left < last && !cmp(*right, *left))
{
++left;
++right;
}

}

++left;
++right;
}

if (left == last)
{
// 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 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);
}

/*******************************************************************************
* Implements the merging routine. Runs in in time O(n + m), where 'n' is the   *
* the length of the left chunk and 'm' is the length of the right chunk.       *
*******************************************************************************/
template<class RandomIt, class Cmp>
void merge(RandomIt source,
RandomIt target,
const size_t offset,
const size_t left_run_length,
const size_t right_run_length,
Cmp cmp)
{
RandomIt left = source + offset;
RandomIt right = left + left_run_length;

const RandomIt left_bound = right;
const RandomIt right_bound = right + right_run_length;

RandomIt target_iter = target + offset;

while (left < left_bound && right < right_bound)
{
*target_iter++ = cmp(*right, *left) ?
*right++:
*left++;
}

std::copy(left, left_bound, target_iter);
std::copy(right, right_bound, target_iter);
}

/*******************************************************************************
* 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;

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

// Request a buffer.
RandomIt buffer = new value_type[length];
std::copy(first, last, buffer);

// 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;
}
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.
const size_t left_run_length = p_queue->dequeue();
const size_t right_run_length = p_queue->dequeue();

// Do the merge.
merge(source,
target,
offset,
left_run_length,
right_run_length,
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;
}
}
}

delete[] buffer;
}

#endif


main.cpp:

#include <algorithm>
#include <chrono>
#include <iostream>
#include <random>

#include "natural_merge_sort.h"

/*******************************************************************************
* Prints an array.                                                             *
*******************************************************************************/
template<class T>
static void print_int_array(const T* begin, const T* last)
{
while (begin < last)
{
std::cout << *begin++ << " ";
}

std::cout << std::endl;
}

/*******************************************************************************
* Checks that the input range is sorted (is in ascending order).               *
*******************************************************************************/
template<class T, class Cmp>
bool is_sorted(T* begin, T* end, Cmp cmp)
{
while (begin < end - 1)
{
if (cmp(*(begin + 1), *begin))
{
return false;
}

++begin;
}

return true;
}

/*******************************************************************************
* 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();
}

/*******************************************************************************
* Checks that the two ranges are of the same length and content.               *
*******************************************************************************/
template <class T>
bool are_equal(const T* begin1,
const T* end1,
const T* begin2,
const T* end2)
{
if (std::distance(begin1, end1) != std::distance(begin2, end2))
{
return false;
}

while (begin1 < end1)
{
if (*begin1++ != *begin2++)
{
return false;
}
}

return true;
}

/*******************************************************************************
* Compares two integer pointers by the values they point to.                   *
*******************************************************************************/
static bool compare_dereference(const int* a, const int* b)
{
return *a < *b;
}

/*******************************************************************************
* Compares two integers.                                                       *
*******************************************************************************/
static bool compare_int(const int a, const int b)
{
return a < b;
}

/*******************************************************************************
* 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: "
<< 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>
static 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: "
<< 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];
std::copy(array1, array1 + sz, array2);

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

profile_stable_sort(array1, array1 + sz, compare_int);
profile_natural_merge_sort(array2, array2 + sz, compare_int);

std::cout << "Same contents: "
<< are_equal(array1, array1 + sz, array2, array2 + 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];
std::copy(array1, array1 + sz, array2);

profile_stable_sort(array1, array1 + sz, compare_dereference);
profile_natural_merge_sort(array2, array2 + sz, compare_dereference);

std::cout << "Same contents: "
<< are_equal(array1, array1 + sz, array2, array2 + 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];
std::copy(array1, array1 + sz, array2);

profile_stable_sort(array1, array1 + sz, compare_int);
profile_natural_merge_sort(array2, array2 + sz, compare_int);

std::cout << "Same contents: "
<< are_equal(array1, array1 + sz, array2, array2 + 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;
}


One thing that I find odd is that you duplicate a great deal of standard library code to implement your sorting algorithm:

• is_sorted can be replaced by std::is_sorted.

• are_equal can be replaced by std::equal. Note that since C++14, there is an overload taking 4 iterators, just like what you did, so that it can use a $O(1)$ std::distance with random-access iterators.

• compare_int can be replaced by std::less<int>, or even the simpler and more generic std::less<> if you have access to the C++14 library.

• Also, you shouldn't bother with a dedicated compare_dereference function when you can use a lambda instead (here, C++14 generic lambda):

profile_stable_sort(array1, array1 + sz, [](auto* a, auto* b){
return *a < *b;
});

• If you use GCC or Clang, you could use __builtin_clz which may be way faster than manually counting leading zeros. A simple macro check could help you to pick the most suitable implementation:

size_t leading_zeros(const size_t num)
{
#if defined(__GNUC__) || defined(__clang__)
return __builtin_clz(num);
#else
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;
#endif
}

• What comes to __builtin_clz, it's ain't an issue since it is called only once per call to sort routine. Apr 10, 2015 at 15:47

In UnsafeIntQueue, you've declared

size_t *m_buffer;


It may be better to instead declare this as

unique_ptr<size_t[]> m_buffer;


and let unique_ptr manage the allocated memory, rather than managing it manually.

Similarly, in build_run_size_queue, you have declared

UnsafeIntQueue* p_q;


however

unique_ptr<UnsafeIntQueue> p_q;


would be better.

In fact, you might consider doing something different: you could pass UnsafeIntQueue around by value rather than by pointer. Passing these objects around with move semantics may turn out to be more efficient.

(you should probably explicitly deleting the copy constructor and copy assignment, to indicate that they shall not be used that way. Although this will automatically happen implicitly if you follow my advice and change m_buffer to be a unique_ptr)

       p_q->enqueue(right - head);