poplar.h

/*
 * The MIT No Attribution License (MIT-0)
 *
 * Copyright (c) 2017-2020 Morwenn
 *
 * Permission is hereby granted, free of charge, to any person obtaining a copy
 * of this software and associated documentation files (the "Software"), to deal
 * in the Software without restriction, including without limitation the rights
 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
 * copies of the Software, and to permit persons to whom the Software is
 * furnished to do so, subject to the following conditions:
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
 * THE SOFTWARE.
 */
#ifndef POPLAR_HEAP_H_
#define POPLAR_HEAP_H_

////////////////////////////////////////////////////////////
// Headers
////////////////////////////////////////////////////////////
#include <cstddef>
#include <functional>
#include <iterator>
#include <limits>
#include <type_traits>
#include <utility>

namespace poplar
{
    namespace detail
    {
        ////////////////////////////////////////////////////////////
        // Generic helper functions
        ////////////////////////////////////////////////////////////

        // Returns 0 if n == 0 else 2^floor(log2(n))

        template<typename Unsigned>
        constexpr auto bit_floor(Unsigned n) noexcept
            -> Unsigned
        {
            constexpr auto bound = std::numeric_limits<Unsigned>::digits / 2;
            for (std::size_t i = 1 ; i <= bound ; i <<= 1) {
                n |= (n >> i);
            }
            return n & ~(n >> 1);
        }

        // Returns 2^floor(log2(n)), assumes n > 0

        template<typename Unsigned>
        constexpr auto unguarded_bit_floor(Unsigned n) noexcept
            -> Unsigned
        {
            return bit_floor(n);
        }

#if defined(__GNUC__) || defined(__clang__)
        constexpr auto unguarded_bit_floor(unsigned int n) noexcept
            -> unsigned int
        {
            constexpr auto k = std::numeric_limits<unsigned>::digits;
            return 1u << (k - 1 - __builtin_clz(n));
        }

        constexpr auto unguarded_bit_floor(unsigned long n) noexcept
            -> unsigned long
        {
            constexpr auto k = std::numeric_limits<unsigned long>::digits;
            return 1ul << (k - 1 - __builtin_clzl(n));
        }

        constexpr auto unguarded_bit_floor(unsigned long long n) noexcept
            -> unsigned long long
        {
            constexpr auto k = std::numeric_limits<unsigned long long>::digits;
            return 1ull << (k - 1 - __builtin_clzll(n));
        }
#endif

        ////////////////////////////////////////////////////////////
        // Insertion sorts
        ////////////////////////////////////////////////////////////

        // Insertion sort which doesn't check for empty sequences
        template<typename BidirectionalIterator, typename Compare>
        auto unchecked_insertion_sort(BidirectionalIterator first, BidirectionalIterator last,
                            Compare compare)
            -> void
        {
            for (auto cur = std::next(first) ; cur != last ; ++cur) {
                auto sift = cur;
                auto sift_1 = std::prev(cur);

                // Compare first so we can avoid 2 moves for
                // an element already positioned correctly
                if (compare(*sift, *sift_1)) {
                    auto tmp = std::move(*sift);
                    do {
                        *sift = std::move(*sift_1);
                    } while (--sift != first && compare(tmp, *--sift_1));
                    *sift = std::move(tmp);
                }
            }
        }

        template<typename BidirectionalIterator, typename Compare>
        auto insertion_sort(BidirectionalIterator first, BidirectionalIterator last,
                            Compare compare)
            -> void
        {
            if (first == last) return;
            unchecked_insertion_sort(std::move(first), std::move(last), std::move(compare));
        }

        ////////////////////////////////////////////////////////////
        // Poplar heap specific helper functions
        ////////////////////////////////////////////////////////////

        template<typename RandomAccessIterator, typename Size, typename Compare>
        auto sift(RandomAccessIterator first, Size size, Compare compare)
            -> void
        {
            if (size < 2) return;

            auto root = first + (size - 1);
            auto child_root1 = root - 1;
            auto child_root2 = first + (size / 2 - 1);

            while (true) {
                auto max_root = root;
                if (compare(*max_root, *child_root1)) {
                    max_root = child_root1;
                }
                if (compare(*max_root, *child_root2)) {
                    max_root = child_root2;
                }
                if (max_root == root) return;

                using std::swap;
                swap(*root, *max_root);

                size /= 2;
                if (size < 2) return;

                root = max_root;
                child_root1 = root - 1;
                child_root2 = max_root - (size - size / 2);
            }
        }

        template<typename RandomAccessIterator, typename Size, typename Compare>
        auto pop_heap_with_size(RandomAccessIterator first, RandomAccessIterator last,
                                Size size, Size poplar_size, Compare compare)
            -> void
        {
            auto last_root = std::prev(last);
            auto bigger = last_root;
            auto bigger_size = poplar_size;

            // Look for the bigger poplar root
            auto it = first;
            while (true) {
                auto root = std::next(it, poplar_size - 1);
                if (root == last_root) break;
                if (compare(*bigger, *root)) {
                    bigger = root;
                    bigger_size = poplar_size;
                }
                it = std::next(root);

                size -= poplar_size;
                poplar_size = unguarded_bit_floor(size + 1u) - 1u;
            }

            // If a poplar root was bigger than the last one, exchange
            // them and sift
            if (bigger != last_root) {
                std::iter_swap(bigger, last_root);
                sift(bigger - (bigger_size - 1), bigger_size, std::move(compare));
            }
        }
    }

    ////////////////////////////////////////////////////////////
    // Standard-library-style make_heap and sort_heap
    ////////////////////////////////////////////////////////////

    template<typename RandomAccessIterator, typename Compare=std::less<>>
    auto push_heap(RandomAccessIterator first, RandomAccessIterator last, Compare compare={})
        -> void
    {
        using poplar_size_t = std::make_unsigned_t<
            typename std::iterator_traits<RandomAccessIterator>::difference_type
        >;
        poplar_size_t size = std::distance(first, last);

        // Find the size of the poplar that will contain the new element in O(log n)
        poplar_size_t last_poplar_size = detail::bit_floor(size + 1u) - 1u;
        while (size - last_poplar_size != 0) {
            size -= last_poplar_size;
            last_poplar_size = detail::unguarded_bit_floor(size + 1u) - 1u;
        }

        // Sift the new element in its poplar in O(log n)
        detail::sift(std::prev(last, last_poplar_size), last_poplar_size, std::move(compare));
    }

    template<typename RandomAccessIterator, typename Compare=std::less<>>
    auto pop_heap(RandomAccessIterator first, RandomAccessIterator last, Compare compare={})
        -> void
    {
        using poplar_size_t = std::make_unsigned_t<
            typename std::iterator_traits<RandomAccessIterator>::difference_type
        >;
        poplar_size_t size = std::distance(first, last);
        auto poplar_size = detail::bit_floor(size + 1u) - 1u;
        detail::pop_heap_with_size(std::move(first), std::move(last),
                                   size, poplar_size, std::move(compare));
    }

    template<typename RandomAccessIterator, typename Compare=std::less<>>
    auto make_heap(RandomAccessIterator first, RandomAccessIterator last, Compare compare={})
        -> void
    {
        using poplar_diff_t = std::make_unsigned_t<
            typename std::iterator_traits<RandomAccessIterator>::difference_type
        >;

        poplar_diff_t size = std::distance(first, last);
        if (size < 2) return;

        // A sorted collection is a valid poplar heap; whenever the heap
        // is small, using insertion sort should be faster, which is why
        // we start by constructing 15-element poplars instead of 1-element
        // ones as the base case
        constexpr poplar_diff_t small_poplar_size = 15;
        if (size <= small_poplar_size) {
            detail::unchecked_insertion_sort(std::move(first), std::move(last),
                                             std::move(compare));
            return;
        }

        // Determines the "level" of the poplars seen so far; the log2 of this
        // variable will be used to make the binary carry sequence
        poplar_diff_t poplar_level = 1;

        auto it = first;
        auto next = std::next(it, small_poplar_size);
        while (true) {
            // Make a 15 element poplar
            detail::unchecked_insertion_sort(it, next, compare);

            poplar_diff_t poplar_size = small_poplar_size;
            // Bit trick iterate without actually having to compute log2(poplar_level)
            for (auto i = (poplar_level & -poplar_level) >> 1 ; i != 0 ; i >>= 1) {
                it -= poplar_size;
                poplar_size = 2 * poplar_size + 1;
                detail::sift(it, poplar_size, compare);
                ++next;
            }

            if (poplar_diff_t(std::distance(next, last)) <= small_poplar_size) {
                detail::insertion_sort(std::move(next), std::move(last),
                                       std::move(compare));
                return;
            }

            it = next;
            std::advance(next, small_poplar_size);
            ++poplar_level;
        }
    }

    template<typename RandomAccessIterator, typename Compare=std::less<>>
    auto sort_heap(RandomAccessIterator first, RandomAccessIterator last, Compare compare={})
        -> void
    {
        using poplar_size_t = std::make_unsigned_t<
            typename std::iterator_traits<RandomAccessIterator>::difference_type
        >;
        poplar_size_t size = std::distance(first, last);
        if (size < 2) return;

        auto poplar_size = detail::bit_floor(size + 1u) - 1u;
        do {
            detail::pop_heap_with_size(first, last, size, poplar_size, compare);
            --last;
            --size;
            poplar_size = detail::unguarded_bit_floor(size + 1u) - 1u;
        } while (size > 1);
    }

    template<typename RandomAccessIterator, typename Compare=std::less<>>
    auto is_heap_until(RandomAccessIterator first, RandomAccessIterator last, Compare compare={})
        -> RandomAccessIterator
    {
        if (std::distance(first, last) < 2) {
            return last;
        }

        using poplar_diff_t = std::make_unsigned_t<
            typename std::iterator_traits<RandomAccessIterator>::difference_type
        >;
        // Determines the "level" of the poplars seen so far; the log2 of this
        // variable will be used to make the binary carry sequence
        poplar_diff_t poplar_level = 1;

        auto it = first;
        auto next = std::next(it);
        while (true) {
            poplar_diff_t poplar_size = 1;

            // Bit trick iterate without actually having to compute log2(poplar_level)
            for (auto i = (poplar_level & -poplar_level) >> 1 ; i != 0 ; i >>= 1) {
                // Beginning and size of the poplar to track
                it -= poplar_size;
                poplar_size = 2 * poplar_size + 1;

                // Check poplar property against child roots
                auto root = it + (poplar_size - 1);
                auto child_root1 = root - 1;
                if (compare(*root, *child_root1)) {
                    return next;
                }
                auto child_root2 = it + (poplar_size / 2 - 1);
                if (compare(*root, *child_root2)) {
                    return next;
                }

                if (next == last) return last;
                ++next;
            }

            if (next == last) return last;
            it = next;
            ++next;
            ++poplar_level;
        }
    }

    template<typename RandomAccessIterator, typename Compare=std::less<>>
    auto is_heap(RandomAccessIterator first, RandomAccessIterator last, Compare compare={})
        -> bool
    {
        return poplar::is_heap_until(first, last, compare) == last;
    }
}

#endif // POPLAR_HEAP_H_