poplar_sort implemenation using a priority queue

poplar_sort.hpp

#include <cstddef>
#include <iterator>
#include <type_traits>
#include <utility>
#include <vector>

//////////////////////////////////////////////////////////////
// NOTE: if you want to try this algorithm, you will need to
// borrow bits from the internals of cpp-sort, the version in
// this gist only includes the poplar_sort-specific parts

//////////////////////////////////////////////////////////////
// Heap algorithms: sift_up is reimplemented because I needed
// to tweak it for a micro-optimization but the other heap
// algorithms used are otherwise equivalent to the ones in
// libc++

template<int N, typename Compare, typename RandomAccessIterator, typename Projection>
auto sift_up(RandomAccessIterator first, RandomAccessIterator last,
             Compare compare, Projection projection,
             difference_type_t<RandomAccessIterator> len)
    -> void
{
    using utility::iter_move;
    auto&& comp = utility::as_function(compare);
    auto&& proj = utility::as_function(projection);
    using difference_type = difference_type_t<RandomAccessIterator>;

    // We don't always want to sift a value up to the top of the heap:
    // sometimes we know that the top N values are already bigger than
    // the one we are about to sift up, so we can skip those comparisons
    constexpr difference_type nb_elems_to_skip = N;

    // WARNING: we only use it with len > 1, but if we ever change that the
    //          algorithm is supposed to be guarded by a size check
    len = (len - 2) / 2;
    RandomAccessIterator it = first + len;
    if (comp(proj(*it), proj(*--last))) {
        auto value = iter_move(last);
        auto&& value_proj = proj(value);
        do {
            *last = iter_move(it);
            last = it;
            if (len <= nb_elems_to_skip) {
                break;
            }
            len = (len - 1) / 2;
            it = first + len;
        } while (comp(proj(*it), value_proj));
        *last = std::move(value);
    }
}

//////////////////////////////////////////////////////////////
// Supporting types

template<typename RandomAccessIterator>
struct poplar;
template<typename RandomAccessIterator>
struct poplar_ptr;

template<typename RandomAccessIterator>
struct poplar
{
    RandomAccessIterator begin, end;
    std::make_unsigned_t<difference_type_t<RandomAccessIterator>> size;

    // Leave it uninitialized, the algorithm doesn't check against
    // nullptr being dereferencing it anyway
    poplar_ptr<RandomAccessIterator>* back_ref;

    auto root() const
        -> RandomAccessIterator
    {
        auto res = end;
        return --res;
    }
};

template<typename RandomAccessIterator>
struct poplar_ptr
{
    using poplar_iterator = typename std::vector<poplar<RandomAccessIterator>>::iterator;

    poplar_iterator it = nullptr;

    poplar_ptr() = default;

    poplar_ptr(const poplar_ptr&) = delete;
    poplar_ptr operator=(const poplar_ptr&) = delete;

    poplar_ptr(poplar_ptr&& other) noexcept:
        it(std::move(other).it)
    {
        it->back_ref = this;
    }

    poplar_ptr(const poplar_iterator& pop_it) noexcept:
        it(pop_it)
    {
        it->back_ref = this;
    }

    auto operator=(poplar_ptr&& other) noexcept
        -> poplar_ptr&
    {
        it = other.it;
        it->back_ref = this;
        return *this;
    }

    auto operator=(const poplar_iterator& pop_it) noexcept
        -> poplar_ptr&
    {
        it = pop_it;
        it->back_ref = this;
        return *this;
    }
};

template<typename RandomAccessIterator>
auto swap(poplar_ptr<RandomAccessIterator>& lhs, poplar_ptr<RandomAccessIterator>& rhs)
    -> void
{
    // Improves codegen a bit on Clang
    auto tmp = lhs.it;
    lhs = rhs.it;
    rhs = tmp;
}

//////////////////////////////////////////////////////////////
// Poplar sort

template<typename RandomAccessIterator, typename Compare, typename Projection>
auto poplar_sort(RandomAccessIterator first, RandomAccessIterator last,
                 Compare compare, Projection projection)
    -> void
{
    using poplar_t = poplar<RandomAccessIterator>;
    using poplar_size_t = std::make_unsigned_t<difference_type_t<RandomAccessIterator>>;

    // Size of the unsorted subsequence
    poplar_size_t size = std::distance(first, last);
    if (size < 2) return;

    std::vector<poplar<RandomAccessIterator>> poplars;
    poplars.reserve(log2(size + 1) + 1);  // TODO: take care of the overflow?

    //
    // Size of the biggest poplar in the array, which always is a number
    // of the form 2^n - 1
    //
    // It's worth noting that the +1 never causes problems: we're only
    // using unsigned integers, so when size is the biggest representable
    // value for its type, size + 1 == 0 thanks to the behaviour of
    // unsigned overflow; hyperfloor(0) == 0, and subtracting 1 to that
    // gives back the biggest representable value, which happens to be
    // a number of the form 2^n - 1
    //
    poplar_size_t poplar_size = hyperfloor(size + 1u) - 1u;

    // Make the poplar heap
    auto it = first;
    do {
        if (poplar_size_t(std::distance(it, last)) >= poplar_size) {
            auto begin = it;
            auto end = it + poplar_size;
            make_poplar(begin, end, compare, projection);
            poplars.push_back({begin, end, poplar_size});
            it = end;
        } else {
            poplar_size = (poplar_size + 1) / 2 - 1;
        }
    } while (poplar_size > 0);

    // The next section sorts the poplar heap as follows: every iteration
    // of a loop checks the root of every poplar to find the biggest
    // element and swaps it with the root of the last poplar, turning the
    // poplar with the biggest root into a semipoplar, which is then
    // adjusted again to become a proper poplar

    // An array is used to store the poplar's size and position, and a
    // binary heap is used as a priority queue to avoid having to check
    // the root of every poplar at every iteration; the algorithm therefore
    // proceeds as follows: extract the biggest value from the heap, swap
    // it with the last root, adjust the semipoplar, adjust the heap, then
    // compute the new poplars and add their roots to the heap

    // Function to compare roots of poplars
    auto&& comp = utility::as_function(compare);
    auto&& proj = utility::as_function(projection);
    auto ind_comp = [&comp, &proj](const auto& lhs, const auto& rhs) {
        return comp(proj(*lhs.it->root()), proj(*rhs.it->root()));
    };

    // Create a priority queue for poplar roots
    std::vector<poplar_ptr<RandomAccessIterator>> poplars_by_max_root;
    poplars_by_max_root.reserve(log2(size + 1) + 1); // TODO: take care of the overflow
    for (auto it = poplars.begin() ; it != poplars.end() ; ++it) {
        poplars_by_max_root.emplace_back(it);
    }

    detail::make_heap(
        poplars_by_max_root.begin(), poplars_by_max_root.end(),
        ind_comp, utility::identity{}
    );

    std::size_t one_element = 0;
    std::size_t several_elements = 0;

    // Sort the poplar heap
    do {
        // Retrieve the max element from the priority queue
        auto last = std::prev(poplars.end());
        auto bigger = poplars_by_max_root.front().it;

        if (bigger != last) {
            // Swap and sift if needed
            using utility::iter_swap;
            iter_swap(bigger->root(), last->root());
            sift(bigger->begin, bigger->size, compare, projection);

            // Adjust pointers to make sure that the top of the heap
            // points to the last poplar
            auto save_ptr = std::next(last->back_ref);
            using std::swap;
            swap(*last->back_ref, poplars_by_max_root.front());

            // The value that was sift up from the poplar is likely not in its
            // correct place in the heap, so we need to sift it up: however we
            // know that the top of the heap is bigger than this value, so we
            // only sift if it is after the  third element, and we also skip
            // comparisons with the top of the heap
            auto len = save_ptr - poplars_by_max_root.data();
            if (len > 3) {
                detail::sift_up<1>(
                    poplars_by_max_root.data(), save_ptr,
                    ind_comp, utility::identity{},
                    len
                );
            }
        }

        // If the last poplar had one element, destroy it
        if (poplars.back().size == 1) {
            // This is the common case, the other ones (see later conditions)
            // occur so infrequently that we can just optimize this case

            // pop_heap, except that we avoid swapping back the front value
            // to the back of the array
            poplars_by_max_root.front() = poplars_by_max_root.back().it;
            poplars_by_max_root.pop_back();
            detail::sift_down(
                poplars_by_max_root.begin(), poplars_by_max_root.end(),
                ind_comp, utility::identity{},
                poplars_by_max_root.size(), poplars_by_max_root.begin()
            );
            // Remove the last poplar
            poplars.pop_back();

            // The following cases occur only a few times during the sort
            if (poplars.size() < 2) {
                if (poplars.size() == 0) {
                    return;
                }
                // assume poplars.size() == 1
                if (poplars.back().size == 1) {
                    return;
                }
                auto& back = poplars.back();
                auto old_end = back.end;
                auto new_size = (back.size - 1) / 2;
                auto middle = back.begin + new_size;
                back.end = middle;
                back.size = new_size;
                poplars.push_back({middle, --old_end, new_size});

                // Add the new poplar roots to the priority queue
                // poplars.size() == 2, just put two elements in reverse order
                if (comp(proj(*poplars[0].root()), proj(*poplars[1].root()))) {
                    poplars_by_max_root[0] = std::next(poplars.begin());
                    poplars_by_max_root.emplace_back(poplars.begin());
                } else {
                    poplars_by_max_root[0] = poplars.begin();
                    poplars_by_max_root.emplace_back(std::next(poplars.begin()));
                }
            }
        } else {
            // Split the remains of the last poplar in two
            auto& back = poplars.back();
            auto old_end = back.end;
            auto new_size = (back.size - 1) / 2;
            auto middle = back.begin + new_size;
            back.end = middle;
            back.size = new_size;
            poplars.push_back({middle, --old_end, new_size});

            // Add the new poplar roots to the priority queue: the beginning of the
            // second-to-last poplar is already in place so the top of the heap is
            // already pointing to the new element, which means that we can directly
            // sift down without having to explicitly assign anything (since we didn't
            // remove said poplar but replaced its field, the back_ref is also already
            // pointing to the right memory location)
            detail::sift_down(
                poplars_by_max_root.begin(), poplars_by_max_root.end(),
                ind_comp, utility::identity{},
                poplars_by_max_root.size(), poplars_by_max_root.begin()
            );
            // Add a the new last poplar to the heap
            poplars_by_max_root.emplace_back(std::prev(poplars.end()));
            detail::sift_up<0>(
                poplars_by_max_root.begin(), poplars_by_max_root.end(),
                ind_comp, utility::identity{},
                poplars_by_max_root.size()
            );
        }
    } while (poplars.size() > 1);
}