code:

#include <time.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>

#include <algorithm>
#include <thread>
#include <vector>

template <class T>
void find_mids(T first1, T last1, T first2, T last2,
               T *mid1_out, T *mid2_out) {
    size_t n = ((last1 - first1) + (last2 - first2)) / 2;

    T w1 = first1;
    T w2 = first2;

    while (n != 0) {
        size_t m = n - n / 2;
        if (m > last1 - w1) {
            w2 += m;
        } else if (m > last2 - w2) {
            w1 += m;
        } else {
            if (*(w1 + (m - 1)) < *(w2 + (m - 1))) {
                w1 += m;
            } else {
                w2 += m;
            }
        }
        n -= m;
    }
    *mid1_out = w1;
    *mid2_out = w2;
}

template <class T, class U>
void recombine(T first1, T last1, T first2, T last2,
               U dest_first, int levels) {
    T mid1, mid2;
    find_mids(first1, last1, first2, last2, &mid1, &mid2);

    auto first_part_size = (mid1 - first1) + (mid2 - first2);

    if (levels == 1) {
        std::thread t(std::merge<T, T, U>,
                      first1, mid1, first2, mid2,
                      dest_first);
        std::merge(mid1, last1, mid2, last2,
                   dest_first + first_part_size);
        t.join();
    } else {
        std::thread t(recombine<T, U>, first1, mid1,
                      first2, mid2,
                      dest_first, levels - 1);
        recombine(mid1, last1, mid2, last2,
                  dest_first + first_part_size,
                  levels - 1);
        t.join();
    }
}

template <class T, class U>
void do_fastsort(T first, T last,
                 U aux_first, U aux_last,
                 int levels) {
    if (levels == 0) {
        std::sort(first, last);
    } else {
        T mid = first + (last - first) / 2;
        U aux_mid = aux_first
            + (aux_last - aux_first) / 2;
        std::thread t(&do_fastsort<T, U>,
                      first, mid, aux_first, aux_mid,
                      levels - 1);
        do_fastsort(mid, last, aux_mid, aux_last,
                    levels - 1);
        t.join();

        if (levels % 2 == 1) {
            recombine(first, mid, mid, last,
                      aux_first, levels);
        } else {
            recombine(aux_first, aux_mid, aux_mid,
                      aux_last, first, levels);
        }
    }
}

template <class T>
void fastsort(T first, T last) {
    using V = typename T::value_type;
    std::unique_ptr<V[]> aux(new V[last - first]);

    do_fastsort(first, last,
                aux.get(), aux.get() + (last - first), 6);
}

double clock_secs() {
    struct timespec ts;
    clock_gettime(CLOCK_REALTIME, &ts);
    return ts.tv_sec
        + static_cast<double>(ts.tv_nsec) / 1000000000;
}

int main() {
    srand(time(NULL));
    puts("Generating randoms...");
    std::vector<int> vec(pow(2, 26));

    for (int& x : vec) {
        x = rand();
    }

    puts("Sorting...");
    double start = clock_secs();
    fastsort(vec.begin(), vec.end());
    double duration = clock_secs() - start;

    if (!std::is_sorted(vec.begin(), vec.end())) {
        printf("Oops, not sorted!\n");
    } else {
        printf("Sort runtime: %f\n", duration);
    }
}