PPC sorting algorithm research attempt.

main.cppp

#include <algorithm>
#include <chrono>
#include <cstring>
#include <iostream>
#include <limits.h>
#include <math.h>
#include <random>
#include <time.h>
#include <limits>
#include <string.h>

#define N (100 * 1000 * 1000)
#define INSERTION_SORT_THRESHOLD 64
#define RADIX 256

using std::cout;
using std::memcpy;

typedef unsigned long long data_t;

static data_t* create_data() {
    data_t* data = new data_t[N];
    std::random_device dev;
    std::mt19937 rng(13);
    std::uniform_int_distribution<int> dist(0, 100000000);

    for (size_t i = 0; i < N; i++) {
        data[i] = dist(rng);
    }

    return data;
}

static void merge_impl(data_t* source,
                       data_t* target,
                       const size_t offset,
                       const size_t left_run_length,
                       const size_t right_run_length)
{
    data_t* left_subarray_pointer = source + offset;
    data_t* right_subarray_pointer = left_subarray_pointer
                                   + left_run_length;

    data_t* left_bound = right_subarray_pointer;
    data_t* right_bound = right_subarray_pointer
                        + right_run_length;

    data_t* target_subarray_pointer = target + offset;

    while (left_subarray_pointer != left_bound &&
           right_subarray_pointer != right_bound) {
        if (*right_subarray_pointer < *left_subarray_pointer) {
            *target_subarray_pointer = *right_subarray_pointer;
            right_subarray_pointer++;
        } else {
            *target_subarray_pointer = *left_subarray_pointer;
            left_subarray_pointer++;
        }

        target_subarray_pointer++;
    }

    memcpy(target_subarray_pointer,
           left_subarray_pointer,
           (left_bound - left_subarray_pointer) * sizeof(data_t));

    memcpy(target_subarray_pointer,
           right_subarray_pointer,
           (right_bound - right_subarray_pointer) * sizeof(data_t));
}

static void insertion_sort(data_t* source,
                           const size_t range_length) {

    for (size_t i = 1; i < range_length;) {
        const data_t x = source[i];
        size_t j = i;

        while (j > 0 && source[j - 1] > x) {
            source[j] = source[j - 1];
            j--;
        }

        source[j] = x;
        i++;
    }
}

static void mergesort_impl(data_t* source,
                           data_t* target,
                           const size_t offset,
                           const size_t range_length) {
    if (range_length < 2) {
        return;
    }

    const size_t half_range_length = range_length / 2;

    mergesort_impl(target,
                   source,
                   offset,
                   half_range_length);

    mergesort_impl(target,
                   source,
                   offset + half_range_length,
                   range_length - half_range_length);

    merge_impl(source,
               target,
               offset,
               half_range_length,
               range_length - half_range_length);
}

void my_mergesort(data_t* base, size_t num)
{
    auto aux = new data_t[num];
    memcpy(aux, base, num * sizeof(data_t));
    mergesort_impl(
        aux,
        base,
        0,
        num);
    delete[] aux;
}

static unsigned long long get_millis() {
    long ms;
    time_t s;
    struct timespec spec;
    clock_gettime(CLOCK_REALTIME, &spec);
    s = spec.tv_sec;
    ms = round(spec.tv_nsec / 1e6);
    if (ms > 999) {
        s++;
        ms -= 1000;
    }

    return s * 1000 + ms;
}

static size_t get_bucket_index(const data_t key,
                               const size_t byte_index) {
    return static_cast<size_t>(key >> (byte_index * CHAR_BIT)) &
           static_cast<size_t>(0xff);
}

static void radix_sort_impl(data_t* source,
                            data_t* target,
                            const size_t num,
                            const size_t byte_index) {

    if (num < INSERTION_SORT_THRESHOLD) {
        insertion_sort(source, num);

        if ((byte_index & 1) == 0) {
            memcpy(target,
                   source,
                   num * sizeof(data_t));
        }

        return;
    }

    size_t bucket_size_map[RADIX];
    size_t start_index_map[RADIX];
    size_t processed_map  [RADIX];

    memset(bucket_size_map, 0, sizeof(size_t) * RADIX);
    memset(processed_map,   0, sizeof(size_t) * RADIX);

    for (data_t* pd = source; pd != source + num; pd++) {
        const data_t current_key = *pd;
        bucket_size_map[get_bucket_index(current_key, byte_index)]++;
    }

    start_index_map[0] = 0;

    for (size_t i = 1; i < RADIX; i++) {
        start_index_map[i] = start_index_map[i - 1]
                           + bucket_size_map[i - 1];
    }

    for (data_t* pd = source; pd != source + num; pd++) {
        const data_t current_key = *pd;
        const size_t bucket_index = get_bucket_index(current_key,
                                                     byte_index);
        *(target + (start_index_map[bucket_index]
                 + processed_map[bucket_index]++)) = current_key;
    }

    if (byte_index) {
        for (size_t i = 0; i != RADIX; i++) {
            if (bucket_size_map[i]) {
                radix_sort_impl(
                    target,
                    source,
                    bucket_size_map[i],
                    byte_index - 1);
            }
        }
    }
}

void radix_sort(data_t* source, size_t num) {
    if (num < 2) {
        return;
    }

    auto aux = new data_t[num];
    memcpy(aux, source, num * sizeof(data_t));
    radix_sort_impl(aux,
                    source,
                    num,
                    3);
    delete[] aux;
}

int main() {
    data_t* data = create_data();
    data_t* copy = new data_t[N];
    data_t* cpy2 = new data_t[N];
    memcpy(copy, data, N * sizeof(data_t));
    memcpy(cpy2, data, N * sizeof(data_t));

    auto st = get_millis();
    std::sort(copy, copy + N);
    auto et = get_millis();
    std::cout << "std::sort in " << et - st << " milliseconds.\n";

    st = get_millis();
    my_mergesort(data, N);
    et = get_millis();
    cout << "my_mergesort in " << et - st << " milliseconds.\n";

    st = get_millis();
    radix_sort(cpy2, N);
    et = get_millis();
    cout << "radix_sort in " << et - st << " milliseconds.\n";

    cout << "Algorithms agree: "
         << std::boolalpha
         << (std::equal(copy, copy + N, data) &&
             std::equal(cpy2, cpy2 + N, data))
         << ".\n";

    delete[] data;
    delete[] copy;
    delete[] cpy2;

    return 0;
}