bp_vector.cpp

#include "bp_vector.hpp"
#include "util.hpp"

namespace succinct {

    namespace {

        // XXX(ot): remove useless tables

        class excess_tables
        {
        public:
            excess_tables() {
                for (int c = 0; c < 256; ++c) {
                    for (uint8_t i = 0; i < 9; ++i) {
                        m_fwd_pos[c][i] = 0;
                        m_bwd_pos[c][i] = 0;
                    }
                    // populate m_fwd_pos, m_fwd_min, and m_fwd_min_idx
                    int excess = 0;
                    m_fwd_min[c] = 0;
                    m_fwd_min_idx[c] = 0;

                    for (char i = 0; i < 8; ++i) {
                        if ((c >> i) & 1) { // opening
                            ++excess;
                        } else { // closing
                            --excess;
                            if (excess < 0 &&
                                m_fwd_pos[c][-excess] == 0) { // not already found
                                m_fwd_pos[c][-excess] = uint8_t(i + 1);
                            }
                        }

                        if (-excess > m_fwd_min[c]) {
                            m_fwd_min[c] = uint8_t(-excess);
                            m_fwd_min_idx[c] = uint8_t(i + 1);
                        }
                    }
                    m_fwd_exc[c] = (char)excess;

                    // populate m_bwd_pos and m_bwd_min
                    excess = 0;
                    m_bwd_min[c] = 0;

                    for (uint8_t i = 0; i < 8; ++i) {
                        if ((c << i) & 128) { // opening
                            ++excess;
                            if (excess > 0 &&
                                m_bwd_pos[c][(uint8_t)excess] == 0) { // not already found
                                m_bwd_pos[c][(uint8_t)excess] = uint8_t(i + 1);
                            }
                        } else { // closing
                            --excess;
                        }

                        m_bwd_min[c] = uint8_t(std::max(excess, (int)m_bwd_min[c]));
                    }
                }
            }

            char m_fwd_exc[256];

            uint8_t m_fwd_pos[256][9];
            uint8_t m_bwd_pos[256][9];

            uint8_t m_bwd_min[256];
            uint8_t m_fwd_min[256];

            uint8_t m_fwd_min_idx[256];
        };

        const static excess_tables tables;

        inline bool find_close_in_word(uint64_t word, uint64_t byte_counts, bp_vector::excess_t cur_exc, uint64_t& ret)
        {
            assert(cur_exc > 0 && cur_exc <= 64);
            const uint64_t cum_exc_step_8 = (uint64_t(cur_exc) + ((2 * byte_counts - 8 * broadword::ones_step_8) << 8)) * broadword::ones_step_8;

            uint64_t min_exc_step_8 = 0;
            for (size_t i = 0; i < 8; ++i) {
                size_t shift = i * 8;
                min_exc_step_8 |= ((uint64_t)(tables.m_fwd_min[(word >> shift) & 0xFF])) << shift;
            }

            const uint64_t has_result = broadword::leq_step_8(cum_exc_step_8, min_exc_step_8);

            unsigned long shift;
            if (broadword::lsb(has_result, shift)) {
                uint8_t bit_pos = tables.m_fwd_pos[(word >> shift) & 0xFF][(cum_exc_step_8 >> shift) & 0xFF];
                assert(bit_pos > 0);
                ret = shift + bit_pos - 1;
                return true;
            }
            return false;
        }

        inline bool find_open_in_word(uint64_t word, uint64_t byte_counts, bp_vector::excess_t cur_exc, uint64_t& ret) {
            assert(cur_exc > 0 && cur_exc <= 64);
            const uint64_t rev_byte_counts = broadword::reverse_bytes(byte_counts);
            const uint64_t cum_exc_step_8 = (uint64_t(cur_exc) - ((2 * rev_byte_counts - 8 * broadword::ones_step_8) << 8)) * broadword::ones_step_8;

            uint64_t max_exc_step_8 = 0;
            for (size_t i = 0; i < 8; ++i) {
                size_t shift = i * 8;
                max_exc_step_8 |= ((uint64_t)(tables.m_bwd_min[(word >> (64 - shift - 8)) & 0xFF])) << shift;
            }

            const uint64_t has_result = broadword::leq_step_8(cum_exc_step_8, max_exc_step_8);

            unsigned long shift;
            if (broadword::lsb(has_result, shift)) {
                uint8_t bit_pos = tables.m_bwd_pos[(word >> (64 - shift - 8)) & 0xFF][(cum_exc_step_8 >> shift) & 0xFF];
                assert(bit_pos > 0);
                ret = 64 - (shift + bit_pos);
                return true;
            }
            return false;
        }

        inline void
        excess_rmq_in_word(uint64_t word, bp_vector::excess_t& exc, uint64_t word_start,
                           bp_vector::excess_t& min_exc, uint64_t& min_exc_idx)
        {
            bp_vector::excess_t min_byte_exc = min_exc;
            uint64_t min_byte_idx = 0;

            for (size_t i = 0; i < 8; ++i) {
                size_t shift = i * 8;
                size_t byte = (word >> shift) & 0xFF;
                // m_fwd_min is negated
                bp_vector::excess_t cur_min = exc - tables.m_fwd_min[byte];

                min_byte_idx = (cur_min < min_byte_exc) ? i : min_byte_idx;
                min_byte_exc = (cur_min < min_byte_exc) ? cur_min : min_byte_exc;

                exc += tables.m_fwd_exc[byte];
            }

            if (min_byte_exc < min_exc) {
                min_exc = min_byte_exc;
                uint64_t shift = min_byte_idx * 8;
                min_exc_idx = word_start + shift + tables.m_fwd_min_idx[(word >> shift) & 0xFF];
            }
        }
    }

    inline bool bp_vector::find_close_in_block(uint64_t block_offset, bp_vector::excess_t excess, uint64_t start, uint64_t& ret) const {
        if (excess > excess_t((bp_block_size - start) * 64)) {
            return false;
        }
        assert(excess > 0);
        for (uint64_t sub_block_offset = start; sub_block_offset < bp_block_size; ++sub_block_offset) {
            uint64_t sub_block = block_offset + sub_block_offset;
            uint64_t word = m_bits[sub_block];
            uint64_t byte_counts = broadword::byte_counts(word);
            assert(excess > 0);
            if (excess <= 64) {
                if (find_close_in_word(word, byte_counts, excess, ret)) {
                    ret += sub_block * 64;
                    return true;
                }
            }
            excess += static_cast<excess_t>(2 * broadword::bytes_sum(byte_counts) - 64);
        }
        return false;
    }

    uint64_t bp_vector::find_close(uint64_t pos) const
    {
        assert((*this)[pos]); // check there is an opening parenthesis in pos
        uint64_t ret = -1U;
        // Search in current word
        uint64_t word_pos = (pos + 1) / 64;
        uint64_t shift = (pos + 1) % 64;
        uint64_t shifted_word = m_bits[word_pos] >> shift;
        // Pad with "open"
        uint64_t padded_word = shifted_word | (-!!shift & (~0ULL << (64 - shift)));
        uint64_t byte_counts = broadword::byte_counts(padded_word);

        excess_t word_exc = 1;
        if (find_close_in_word(padded_word, byte_counts, word_exc, ret)) {
            ret += pos + 1;
            return ret;
        }

        // Otherwise search in the local block
        uint64_t block = word_pos / bp_block_size;
        uint64_t block_offset = block * bp_block_size;
        uint64_t sub_block = word_pos % bp_block_size;
        uint64_t local_rank = broadword::bytes_sum(byte_counts) - shift; // subtract back the padding
        excess_t local_excess = static_cast<excess_t>((2 * local_rank) - (64 - shift));
        if (find_close_in_block(block_offset, local_excess + 1, sub_block + 1, ret)) {
            return ret;
        }

        // Otherwise, find the first appropriate block
        excess_t pos_excess = excess(pos);
        uint64_t found_block = search_min_tree<1>(block + 1, pos_excess);
        uint64_t found_block_offset = found_block * bp_block_size;
        excess_t found_block_excess = get_block_excess(found_block);

        // Search in the found block
        bool found = find_close_in_block(found_block_offset, found_block_excess - pos_excess, 0, ret);
        assert(found); (void)found;
        return ret;
    }

    inline bool bp_vector::find_open_in_block(uint64_t block_offset, bp_vector::excess_t excess, uint64_t start, uint64_t& ret) const {
        if (excess > excess_t(start * 64)) {
            return false;
        }
        assert(excess >= 0);

        for (uint64_t sub_block_offset = start - 1; sub_block_offset + 1 > 0; --sub_block_offset) {
            assert(excess > 0);
            uint64_t sub_block = block_offset + sub_block_offset;
            uint64_t word = m_bits[sub_block];
            uint64_t byte_counts = broadword::byte_counts(word);
            if (excess <= 64) {
                if (find_open_in_word(word, byte_counts, excess, ret)) {
                    ret += sub_block * 64;
                    return true;
                }
            }
            excess -= static_cast<excess_t>(2 * broadword::bytes_sum(byte_counts) - 64);
        }
        return false;
    }

    uint64_t bp_vector::find_open(uint64_t pos) const
    {
        assert(pos);
        uint64_t ret = -1U;
        // Search in current word
        uint64_t word_pos = (pos / 64);
        uint64_t len = pos % 64;
        // Rest is padded with "close"
        uint64_t shifted_word = -!!len & (m_bits[word_pos] << (64 - len));
        uint64_t byte_counts = broadword::byte_counts(shifted_word);

        excess_t word_exc = 1;
        if (find_open_in_word(shifted_word, byte_counts, word_exc, ret)) {
            ret += pos - 64;
            return ret;
        }

        // Otherwise search in the local block
        uint64_t block = word_pos / bp_block_size;
        uint64_t block_offset = block * bp_block_size;
        uint64_t sub_block = word_pos % bp_block_size;
        uint64_t local_rank = broadword::bytes_sum(byte_counts); // no need to subtract the padding
        excess_t local_excess = -static_cast<excess_t>((2 * local_rank) - len);
        if (find_open_in_block(block_offset, local_excess + 1, sub_block, ret)) {
            return ret;
        }

        // Otherwise, find the first appropriate block
        excess_t pos_excess = excess(pos) - 1;
        uint64_t found_block = search_min_tree<0>(block - 1, pos_excess);
        uint64_t found_block_offset = found_block * bp_block_size;
        // Since search is backwards, have to add the current block
        excess_t found_block_excess = get_block_excess(found_block + 1);

        // Search in the found block
        bool found = find_open_in_block(found_block_offset, found_block_excess - pos_excess, bp_block_size, ret);
        assert(found); (void)found;
        return ret;
    }

    template <int direction>
    inline bool bp_vector::search_block_in_superblock(uint64_t block, excess_t excess, size_t& found_block) const
    {
        size_t superblock = block / superblock_size;
        excess_t superblock_excess = get_block_excess(superblock * superblock_size);
        if (direction) {
            for (size_t cur_block = block;
                 cur_block < std::min((superblock + 1) * superblock_size, (size_t)m_block_excess_min.size());
                 ++cur_block) {
                if (excess >= superblock_excess + m_block_excess_min[cur_block]) {
                    found_block = cur_block;
                    return true;
                }
            }
        } else {
            for (size_t cur_block = block;
                 cur_block + 1 >= (superblock * superblock_size) + 1;
                 --cur_block) {
                if (excess >= superblock_excess + m_block_excess_min[cur_block]) {
                    found_block = cur_block;
                    return true;
                }
            }
        }

        return false;
    }

    inline bp_vector::excess_t bp_vector::get_block_excess(uint64_t block) const {
        uint64_t sub_block_idx = block * bp_block_size;
        uint64_t block_pos = sub_block_idx * 64;
        excess_t excess = static_cast<excess_t>(2 * sub_block_rank(sub_block_idx) - block_pos);
        assert(excess >= 0);
        return excess;
    }

    inline bool bp_vector::in_node_range(uint64_t node, excess_t excess) const {
        assert(m_superblock_excess_min[node] != excess_t(size()));
        return excess >= m_superblock_excess_min[node];
    }

    template <int direction>
    inline uint64_t bp_vector::search_min_tree(uint64_t block, excess_t excess) const
    {
        size_t found_block = -1U;
        if (search_block_in_superblock<direction>(block, excess, found_block)) {
            return found_block;
        }

        size_t cur_superblock = block / superblock_size;
        size_t cur_node = m_internal_nodes + cur_superblock;
        while (true) {
            assert(cur_node);
            bool going_back = (cur_node & 1) == direction;
            if (!going_back) {
                size_t next_node = direction ? (cur_node + 1) : (cur_node - 1);
                if (in_node_range(next_node, excess)) {
                    cur_node = next_node;
                    break;
                }
            }
            cur_node /= 2;
        }

        assert(cur_node);

        while (cur_node < m_internal_nodes) {
            uint64_t next_node = cur_node * 2 + (1 - direction);
            if (in_node_range(next_node, excess)) {
                cur_node = next_node;
                continue;
            }

            next_node = direction ? (next_node + 1) : (next_node - 1);
            // if it is not one child, it must be the other
            assert(in_node_range(next_node, excess));
            cur_node = next_node;
        }

        size_t next_superblock = cur_node - m_internal_nodes;
        bool ret = search_block_in_superblock<direction>(next_superblock * superblock_size + (1 - direction) * (superblock_size - 1),
                                                         excess, found_block);
        assert(ret); (void)ret;

        return found_block;
    }


    bp_vector::excess_t
    bp_vector::excess(uint64_t pos) const
    {
        return static_cast<excess_t>(2 * rank(pos) - pos);
    }

    void
    bp_vector::excess_rmq_in_block(uint64_t start, uint64_t end,
                                   bp_vector::excess_t& exc,
                                   bp_vector::excess_t& min_exc,
                                   uint64_t& min_exc_idx) const
    {
        assert(start <= end);
        if (start == end) return;

        assert((start / bp_block_size) == ((end - 1) / bp_block_size));
        for (size_t w = start; w < end; ++w) {
            excess_rmq_in_word(m_bits[w], exc, w * 64,
                               min_exc, min_exc_idx);
        }
    }

    void
    bp_vector::excess_rmq_in_superblock(uint64_t block_start, uint64_t block_end,
                                        bp_vector::excess_t& block_min_exc,
                                        uint64_t& block_min_idx) const
    {
        assert(block_start <= block_end);
        if (block_start == block_end) return;

        uint64_t superblock = block_start / superblock_size;

        assert(superblock == ((block_end - 1) / superblock_size));
        excess_t superblock_excess = get_block_excess(superblock * superblock_size);

        for (uint64_t block = block_start; block < block_end; ++block) {
            if (superblock_excess + m_block_excess_min[block] < block_min_exc) {
                block_min_exc = superblock_excess + m_block_excess_min[block];
                block_min_idx = block;
            }
        }
    }


    void
    bp_vector::find_min_superblock(uint64_t superblock_start, uint64_t superblock_end,
                                   bp_vector::excess_t& superblock_min_exc,
                                   uint64_t& superblock_min_idx) const {

        if (superblock_start == superblock_end) return;

        uint64_t cur_node = m_internal_nodes + superblock_start;
        uint64_t rightmost_span = superblock_start;

        excess_t node_min_exc = m_superblock_excess_min[cur_node];
        uint64_t node_min_idx = cur_node;

        // code below assumes that there is at least one right-turn in
        // the node-root-node path, so we must handle this case
        // separately
        if (superblock_end - superblock_start == 1) {
            superblock_min_exc = node_min_exc;
            superblock_min_idx = superblock_start;
            return;
        }

        // go up the tree until we find the lowest node that spans the
        // whole superblock range
        size_t h = 0;
        while (true) {
            assert(cur_node);

            if ((cur_node & 1) == 0) { // is a left child
                // add right subtree to candidate superblocks
                uint64_t right_sibling = cur_node + 1;
                rightmost_span += uint64_t(1) << h;

                if (rightmost_span < superblock_end &&
                    m_superblock_excess_min[right_sibling] < node_min_exc) {
                    node_min_exc = m_superblock_excess_min[right_sibling];
                    node_min_idx = right_sibling;
                }

                if (rightmost_span >= superblock_end - 1) {
                    cur_node += 1;
                    break;
                }
            }

            cur_node /= 2; // parent
            h += 1;
        }

        assert(cur_node);

        // go down until we reach superblock_end
        while (rightmost_span > superblock_end - 1) {
            assert(cur_node < m_superblock_excess_min.size());
            assert(h > 0);

            h -= 1;
            uint64_t left_child = cur_node * 2;
            uint64_t right_child_span = uint64_t(1) << h;
            if ((rightmost_span - right_child_span) >= (superblock_end - 1)) {
                // go to left child
                rightmost_span -= right_child_span;
                cur_node = left_child;
            } else {
                // go to right child and add left subtree to candidate
                // subblocks
                if (m_superblock_excess_min[left_child] < node_min_exc) {
                    node_min_exc = m_superblock_excess_min[left_child];
                    node_min_idx = left_child;
                }
                cur_node = left_child + 1;
            }
        }

        // check last left-turn
        if (rightmost_span < superblock_end &&
            m_superblock_excess_min[cur_node] < node_min_exc) {
            node_min_exc = m_superblock_excess_min[cur_node];
            node_min_idx = cur_node;
        }

        assert(rightmost_span == superblock_end - 1);

        // now reach the minimum leaf in the found subtree (cur_node),
        // which is entirely contained in the range
        if (node_min_exc < superblock_min_exc) {
            cur_node = node_min_idx;
            while (cur_node < m_internal_nodes) {
                cur_node *= 2;
                // remember that past-the-end nodes are filled with size()
                if (m_superblock_excess_min[cur_node + 1] <
                    m_superblock_excess_min[cur_node]) {
                    cur_node += 1;
                }
            }

            assert(m_superblock_excess_min[cur_node] == node_min_exc);
            superblock_min_exc = node_min_exc;
            superblock_min_idx = cur_node - m_internal_nodes;

            assert(superblock_min_idx >= superblock_start);
            assert(superblock_min_idx < superblock_end);
        }
    }

    uint64_t bp_vector::excess_rmq(uint64_t a, uint64_t b, excess_t& min_exc) const
    {
        assert(a <= b);

        excess_t cur_exc = excess(a);
        min_exc = cur_exc;
        uint64_t min_exc_idx = a;

        if (a == b) {
            return min_exc_idx;
        }

        uint64_t range_len = b - a;

        uint64_t word_a_idx = a / 64;
        uint64_t word_b_idx = (b - 1) / 64;

        // search in word_a
        uint64_t shift_a = a % 64;
        uint64_t shifted_word_a = m_bits[word_a_idx] >> shift_a;
        uint64_t subword_len_a = std::min(64 - shift_a, range_len);

        uint64_t padded_word_a =
            (subword_len_a == 64)
            ? shifted_word_a
            : (shifted_word_a | (~0ULL << subword_len_a));

        excess_rmq_in_word(padded_word_a, cur_exc, a,
                           min_exc, min_exc_idx);

        if (word_a_idx == word_b_idx) {
            // single word
            return min_exc_idx;
        }

        uint64_t block_a = word_a_idx / bp_block_size;
        uint64_t block_b = word_b_idx / bp_block_size;

        cur_exc -= 64 - excess_t(subword_len_a); // remove padding

        if (block_a == block_b) {
            // same block
            excess_rmq_in_block(word_a_idx + 1, word_b_idx,
                                cur_exc, min_exc, min_exc_idx);

        } else {
            // search in partial block of word_a
            excess_rmq_in_block(word_a_idx + 1, (block_a + 1) * bp_block_size,
                                cur_exc, min_exc, min_exc_idx);

            // search in blocks
            excess_t block_min_exc = min_exc;
            uint64_t block_min_idx = -1U;

            uint64_t superblock_a = (block_a + 1) / superblock_size;
            uint64_t superblock_b = block_b / superblock_size;

            if (superblock_a == superblock_b) {
                // same superblock
                excess_rmq_in_superblock(block_a + 1, block_b,
                                         block_min_exc, block_min_idx);
            } else {
                // partial superblock of a
                excess_rmq_in_superblock(block_a + 1,
                                         (superblock_a + 1) * superblock_size,
                                         block_min_exc,
                                         block_min_idx);

                // search min superblock in the min tree
                excess_t superblock_min_exc = min_exc;
                uint64_t superblock_min_idx = -1U;
                find_min_superblock(superblock_a + 1, superblock_b,
                                    superblock_min_exc, superblock_min_idx);

                if (superblock_min_exc < min_exc) {
                    excess_rmq_in_superblock(superblock_min_idx * superblock_size,
                                             (superblock_min_idx + 1) * superblock_size,
                                             block_min_exc,
                                             block_min_idx);
                }

                // partial superblock of b
                excess_rmq_in_superblock(superblock_b * superblock_size,
                                         block_b,
                                         block_min_exc,
                                         block_min_idx);
            }

            if (block_min_exc < min_exc) {
                cur_exc = get_block_excess(block_min_idx);
                excess_rmq_in_block(block_min_idx * bp_block_size,
                                    (block_min_idx + 1) * bp_block_size,
                                    cur_exc, min_exc, min_exc_idx);
                assert(min_exc == block_min_exc);
            }

            // search in partial block of word_b
            cur_exc = get_block_excess(block_b);
            excess_rmq_in_block(block_b * bp_block_size, word_b_idx,
                                cur_exc, min_exc, min_exc_idx);
        }

        // search in word_b
        uint64_t word_b = m_bits[word_b_idx];
        uint64_t offset_b = b % 64;
        uint64_t padded_word_b =
            (offset_b == 0)
            ? word_b
            : (word_b | (~0ULL << offset_b));

        excess_rmq_in_word(padded_word_b, cur_exc, word_b_idx * 64,
                           min_exc, min_exc_idx);

        assert(min_exc_idx >= a);
        assert(min_exc == excess(min_exc_idx));

        return min_exc_idx;
    }


    void bp_vector::build_min_tree()
    {
        if (!size()) return;

        std::vector<block_min_excess_t> block_excess_min;
        excess_t cur_block_min = 0, cur_superblock_excess = 0;
        for (uint64_t sub_block = 0; sub_block < m_bits.size(); ++sub_block) {
            if (sub_block % bp_block_size == 0) {
                if (sub_block % (bp_block_size * superblock_size) == 0) {
                    cur_superblock_excess = 0;
                }
                if (sub_block) {
                    assert(cur_block_min >= std::numeric_limits<block_min_excess_t>::min());
                    assert(cur_block_min <= std::numeric_limits<block_min_excess_t>::max());
                    block_excess_min.push_back((block_min_excess_t)cur_block_min);
                    cur_block_min = cur_superblock_excess;
                }
            }
            uint64_t word = m_bits[sub_block];
            uint64_t mask = 1ULL;
            // for last block stop at bit boundary
            uint64_t n_bits =
                (sub_block == m_bits.size() - 1 && size() % 64)
                ? size() % 64
                : 64;
            // XXX(ot) use tables.m_fwd_{min,max}
            for (uint64_t i = 0; i < n_bits; ++i) {
                cur_superblock_excess += (word & mask) ? 1 : -1;
                cur_block_min = std::min(cur_block_min, cur_superblock_excess);
                mask <<= 1;
            }
        }
        // Flush last block mins
        assert(cur_block_min >= std::numeric_limits<block_min_excess_t>::min());
        assert(cur_block_min <= std::numeric_limits<block_min_excess_t>::max());
        block_excess_min.push_back((block_min_excess_t)cur_block_min);

        size_t n_blocks = util::ceil_div(data().size(), bp_block_size);
        assert(n_blocks == block_excess_min.size());

        size_t n_superblocks = (n_blocks + superblock_size - 1) / superblock_size;

        size_t n_complete_leaves = 1;
        while (n_complete_leaves < n_superblocks) n_complete_leaves <<= 1; // XXX(ot): I'm sure this can be done with broadword::msb...
        // n_complete_leaves is the smallest power of 2 >= n_superblocks
        m_internal_nodes = n_complete_leaves;
        size_t treesize = m_internal_nodes + n_superblocks;

        std::vector<excess_t> superblock_excess_min(treesize);

        // Fill in the leaves of the tree
        for (size_t superblock = 0; superblock < n_superblocks; ++superblock) {
            excess_t cur_super_min = static_cast<excess_t>(size());
            excess_t superblock_excess = get_block_excess(superblock * superblock_size);

            for (size_t block = superblock * superblock_size;
                 block < std::min((superblock + 1) * superblock_size, n_blocks);
                 ++block) {
                cur_super_min = std::min(cur_super_min, superblock_excess + block_excess_min[block]);
            }
            assert(cur_super_min >= 0 && cur_super_min < excess_t(size()));

            superblock_excess_min[m_internal_nodes + superblock] = cur_super_min;
        }

        // fill in the internal nodes with past-the-boundary values
        // (they will also serve as sentinels in debug)
        for (size_t node = 0; node < m_internal_nodes; ++node) {
            superblock_excess_min[node] = static_cast<excess_t>(size());
        }

        // Fill bottom-up the other layers: each node updates the parent
        for (size_t node = treesize - 1; node > 1; --node) {
            size_t parent = node / 2;
            superblock_excess_min[parent] = std::min(superblock_excess_min[parent], // same node
                                                     superblock_excess_min[node]);
        }

        m_block_excess_min.steal(block_excess_min);
        m_superblock_excess_min.steal(superblock_excess_min);
    }
}