hybrid_btree.h

#pragma once
 
#include <pixie/bits.h>
#include <pixie/memory_usage.h>
#include <pixie/rmq.h>
 
#include <algorithm>
#include <array>
#include <bit>
#include <cstddef>
#include <cstdint>
#include <functional>
#include <limits>
#include <span>
#include <stdexcept>
#include <type_traits>
#include <utility>
#include <vector>
 
namespace pixie::rmq {

enum class HybridBTreeLeafSelector {
  PrefixSuffix,
  BP,
};

template <class T,
          class Compare = std::less<T>,
          class Index = std::size_t,
          std::size_t LeafSize = 496,
          std::size_t Fanout = 256,
          HybridBTreeLeafSelector LeafSelectorKind =
              HybridBTreeLeafSelector::PrefixSuffix>
class HybridBTree
    : public RmqBase<
          HybridBTree<T, Compare, Index, LeafSize, Fanout, LeafSelectorKind>,
          T> {
 public:
  static_assert(std::is_unsigned_v<Index>,
                "HybridBTree index type must be unsigned");
  static constexpr bool kBpLeafSelector =
      LeafSelectorKind == HybridBTreeLeafSelector::BP;
  static constexpr bool kMaskLeafSelector =
      LeafSelectorKind == HybridBTreeLeafSelector::PrefixSuffix;
  static_assert((kBpLeafSelector && LeafSize == 252) ||
                    (kMaskLeafSelector && (LeafSize == 248 || LeafSize == 496)),
                "HybridBTree requires 252-value BP leaves "
                "or 248/496-value prefix/suffix mask leaves");
  static_assert(Fanout == 256,
                "HybridBTree currently requires a 256 fanout template "
                "argument");
  static_assert(kBpLeafSelector || kMaskLeafSelector,
                "unsupported HybridBTree low-level selector kind");
 
  using Self =
      HybridBTree<T, Compare, Index, LeafSize, Fanout, LeafSelectorKind>;
 
  static constexpr std::size_t npos = RmqBase<Self, T>::npos;
  static constexpr Index invalid_index = std::numeric_limits<Index>::max();
  static constexpr std::size_t kLeafSize = LeafSize;
  static constexpr std::size_t kFanout = Fanout;
  static constexpr std::size_t kMiddleFanout = 192;
  static constexpr std::size_t kMinTopSparseBlockSize = 4096;
  static constexpr std::size_t kMaxTopSparseBlocks = std::size_t{1} << 14;

  HybridBTree() = default;

  HybridBTree(const HybridBTree&) = default;

  HybridBTree(HybridBTree&&) noexcept = default;

  HybridBTree& operator=(const HybridBTree&) = default;

  HybridBTree& operator=(HybridBTree&&) noexcept = default;

  explicit HybridBTree(std::span<const T> values, Compare compare = Compare())
      : values_(values), compare_(compare) {
    build();
  }

  std::size_t size_impl() const { return values_.size(); }

  T value_at_impl(std::size_t position) const { return values_[position]; }

  std::size_t arg_min_impl(std::size_t left, std::size_t right) const {
    if (left >= right || right > values_.size() || level_sizes_.empty()) {
      return npos;
    }
 
    if (right - left <= top_block_size_) {
      return tree_arg_min(left, right);
    }
    const std::size_t top_answer = top_sparse_arg_min(left, right);
    return top_answer != npos ? top_answer : tree_arg_min(left, right);
  }

  static std::size_t top_sparse_block_size_for(std::size_t value_count) {
    if (value_count == 0) {
      return kMinTopSparseBlockSize;
    }
    return std::max(kMinTopSparseBlockSize,
                    1 + (value_count - 1) / kMaxTopSparseBlocks);
  }

  static std::size_t top_sparse_block_count_for(std::size_t value_count) {
    if (value_count == 0) {
      return 0;
    }
    const std::size_t block_size = top_sparse_block_size_for(value_count);
    return 1 + (value_count - 1) / block_size;
  }

  std::size_t top_sparse_block_size() const { return top_block_size_; }

  std::size_t top_sparse_block_count() const { return top_block_count_; }

  std::size_t memory_usage_bytes_impl() const {
    std::size_t bytes = sizeof(*this);
    bytes += pixie::vector_capacity_bytes(leaf_selectors_);
    bytes += pixie::vector_capacity_bytes(medium_selectors_);
    bytes += pixie::vector_capacity_bytes(medium_min_values_);
    bytes += pixie::vector_capacity_bytes(top_sparse_candidates_);
    bytes += pixie::vector_capacity_bytes(medium_level_offsets_);
    bytes += pixie::vector_capacity_bytes(level_sizes_);
    bytes += pixie::vector_capacity_bytes(level_value_spans_);
    bytes += pixie::vector_capacity_bytes(level_fanouts_);
    return bytes;
  }
 
 private:
  static constexpr std::size_t kSelectorEntries = 256;
  static constexpr std::size_t kSelectorBits = 2 * kSelectorEntries;
  static constexpr std::size_t kSelectorWords = kSelectorBits / 64;
  static constexpr std::size_t kEmbeddedOffsetEntries = 252;
  static constexpr std::size_t kEmbeddedPositionEntries = 192;
  static constexpr std::size_t kEmbeddedOffsetBit = 2 * kEmbeddedOffsetEntries;
  static constexpr std::size_t kEmbeddedPositionBit =
      2 * kEmbeddedPositionEntries;
  static constexpr std::size_t kMiddleZeroPrefixWord = 7;
  static constexpr std::size_t kLeafLinearScanThreshold = 64;
  static constexpr std::size_t kLeafAvx2ScanThreshold = 16;
  static constexpr std::uint64_t kEmbeddedOffsetMask =
      std::uint64_t{0xFF} << (kEmbeddedOffsetBit & 63);
  static constexpr bool kInvalidIndexEqualsNpos =
      static_cast<std::size_t>(invalid_index) == npos;
 
  static_assert(kEmbeddedOffsetBit + 8 == kSelectorBits);
  static_assert(kEmbeddedPositionBit + 128 == kSelectorBits);
  static_assert(!kBpLeafSelector || LeafSize <= kEmbeddedOffsetEntries);
  static_assert(kMiddleFanout <= kEmbeddedPositionEntries);
 
  struct MinCandidate {
    std::size_t position = npos;
    const T* value = nullptr;
  };
 
  struct TopCandidate {
    Index position = invalid_index;
  };
  static_assert(sizeof(TopCandidate) == sizeof(Index));
 
  class alignas(64) Bp512Selector {
   public:
    Bp512Selector() = default;

    template <class EntryLess>
    void build(std::size_t entry_count, EntryLess entry_less) {
      if (entry_count > kSelectorEntries) {
        throw std::length_error("HybridBTree local selector too large");
      }
 
      bp_bits_.fill(0);
      if (entry_count == 0) {
        return;
      }
 
      std::array<std::uint16_t, kSelectorEntries> stack{};
      std::size_t stack_size = 0;
      std::size_t write_position = 2 * entry_count;
 
      for (std::size_t i = entry_count; i-- > 0;) {
        while (stack_size != 0 && !entry_less(stack[stack_size - 1], i)) {
          --stack_size;
          prepend_bp_bit(write_position, true);
        }
        stack[stack_size++] = static_cast<std::uint16_t>(i);
        prepend_bp_bit(write_position, false);
      }
 
      while (write_position != 0) {
        prepend_bp_bit(write_position, true);
      }
    }

    std::size_t arg_min(std::size_t slot_left,
                        std::size_t slot_right,
                        std::size_t entry_count) const {
      if (slot_left >= slot_right || slot_right > entry_count ||
          entry_count > kSelectorEntries) {
        return npos;
      }
      if (slot_left + 1 == slot_right) {
        return slot_left;
      }
 
      const std::size_t bit_count = 2 * entry_count;
      const std::size_t first_close = close_position(slot_left);
      const std::size_t last_close = close_position(slot_right - 1);
      if (first_close > last_close || last_close >= bit_count) {
        return npos;
      }
 
      const std::size_t shifted_min =
          depth_arg_min(first_close + 1, last_close + 2, bit_count);
      if (shifted_min == npos || shifted_min == 0) {
        return npos;
      }
 
      const std::size_t zero_rank = rank0_at(shifted_min, bit_count);
      if (zero_rank == 0) {
        return npos;
      }
      const std::size_t entry = zero_rank - 1;
      return entry < entry_count ? entry : npos;
    }

    std::size_t arg_min_with_zero_prefix(std::size_t slot_left,
                                         std::size_t slot_right,
                                         std::size_t entry_count) const {
      if (slot_left >= slot_right || slot_right > entry_count ||
          entry_count > kEmbeddedPositionEntries) {
        return npos;
      }
      if (slot_left + 1 == slot_right) {
        return slot_left;
      }
 
      const std::size_t bit_count = 2 * entry_count;
      const std::size_t first_close =
          close_position_with_zero_prefix(slot_left, bit_count);
      const std::size_t last_close =
          close_position_with_zero_prefix(slot_right - 1, bit_count);
      if (first_close > last_close || last_close >= bit_count) {
        return npos;
      }
 
      const std::size_t shifted_min = depth_arg_min_with_zero_prefix(
          first_close + 1, last_close + 2, bit_count);
      if (shifted_min == npos || shifted_min == 0) {
        return npos;
      }
 
      const std::size_t zero_rank =
          rank0_at_with_zero_prefix(shifted_min, bit_count);
      if (zero_rank == 0) {
        return npos;
      }
      const std::size_t entry = zero_rank - 1;
      return entry < entry_count ? entry : npos;
    }

    std::size_t close_position(std::size_t slot) const {
      return select0_512(bp_bits_.data(), slot);
    }

    std::size_t rank0_at(std::size_t position, std::size_t bit_count) const {
      position = std::min(position, bit_count);
      return position - rank_512(bp_bits_.data(), position);
    }

    void set_embedded_min_offset(std::size_t offset) {
      bp_bits_[kEmbeddedOffsetBit >> 6] =
          (bp_bits_[kEmbeddedOffsetBit >> 6] & ~kEmbeddedOffsetMask) |
          ((static_cast<std::uint64_t>(offset) & 0xFFu)
           << (kEmbeddedOffsetBit & 63));
    }

    std::uint8_t embedded_min_offset() const {
      return static_cast<std::uint8_t>(
          (bp_bits_[kEmbeddedOffsetBit >> 6] & kEmbeddedOffsetMask) >>
          (kEmbeddedOffsetBit & 63));
    }

    void set_embedded_min_position(std::size_t position) {
      bp_bits_[kEmbeddedPositionBit >> 6] =
          static_cast<std::uint64_t>(position);
    }

    std::size_t embedded_min_position() const {
      return static_cast<std::size_t>(bp_bits_[kEmbeddedPositionBit >> 6]);
    }

    void build_zero_prefix_metadata(std::size_t bit_count) {
      const std::size_t word_count = (bit_count + 63) / 64;
      std::uint64_t packed = 0;
      std::size_t zeros = 0;
      for (std::size_t word = 0; word <= word_count; ++word) {
        packed |= (static_cast<std::uint64_t>(zeros) & 0xFFu) << (8 * word);
        if (word == word_count) {
          break;
        }
        const std::size_t word_begin = word * 64;
        const std::size_t word_bits =
            std::min<std::size_t>(64, bit_count - word_begin);
        const std::uint64_t bits = bp_bits_[word] & first_bits_mask(word_bits);
        zeros += word_bits - std::popcount(bits);
      }
      bp_bits_[kMiddleZeroPrefixWord] = packed;
    }
 
   private:
    std::size_t prepend_bp_bit(std::size_t& write_position, bool bit) {
      --write_position;
      if (bit) {
        bp_bits_[write_position >> 6] |= std::uint64_t{1}
                                         << (write_position & 63);
      }
      return write_position;
    }

    int prefix_excess(std::size_t position) const {
      position = std::min(position, kSelectorBits);
      const std::size_t ones = rank_512(bp_bits_.data(), position);
      return static_cast<int>(2 * ones) - static_cast<int>(position);
    }

    std::size_t zero_prefix_at_word(std::size_t word) const {
      return static_cast<std::uint8_t>(bp_bits_[kMiddleZeroPrefixWord] >>
                                       (8 * word));
    }

    std::size_t rank0_at_with_zero_prefix(std::size_t position,
                                          std::size_t bit_count) const {
      position = std::min(position, bit_count);
      const std::size_t full_words = position >> 6;
      std::size_t zeros = zero_prefix_at_word(full_words);
      const std::size_t tail_bits = position & 63;
      if (tail_bits != 0) {
        const std::uint64_t tail =
            bp_bits_[full_words] & first_bits_mask(tail_bits);
        zeros += tail_bits - std::popcount(tail);
      }
      return zeros;
    }

    int prefix_excess_with_zero_prefix(std::size_t position,
                                       std::size_t bit_count) const {
      position = std::min(position, bit_count);
      return static_cast<int>(position) -
             2 * static_cast<int>(
                     rank0_at_with_zero_prefix(position, bit_count));
    }

    std::size_t close_position_with_zero_prefix(std::size_t slot,
                                                std::size_t bit_count) const {
      const std::size_t word_count = (bit_count + 63) / 64;
      for (std::size_t word = 0; word < word_count; ++word) {
        const std::size_t next_zero_prefix = zero_prefix_at_word(word + 1);
        if (next_zero_prefix <= slot) {
          continue;
        }
        const std::size_t local_rank = slot - zero_prefix_at_word(word);
        const std::size_t word_begin = word * 64;
        const std::size_t word_bits =
            std::min<std::size_t>(64, bit_count - word_begin);
        const std::uint64_t zeros =
            (~bp_bits_[word]) & first_bits_mask(word_bits);
        return word_begin + select_64(zeros, local_rank);
      }
      return npos;
    }

    std::size_t depth_arg_min(std::size_t left,
                              std::size_t right,
                              std::size_t bit_count) const {
      const std::size_t depth_count = bit_count + 1;
      if (left >= right || right > depth_count) {
        return npos;
      }
 
      std::size_t position = left;
      int best_depth = prefix_excess(position);
      std::size_t best_position = position;
 
      while (position < right) {
        const std::size_t chunk_begin = (position / 128) * 128;
        const std::size_t local_left = position - chunk_begin;
        const std::size_t local_right =
            std::min<std::size_t>(right - 1, chunk_begin + 128) - chunk_begin;
 
        int candidate_depth;
        std::size_t candidate_position;
        if (chunk_begin >= bit_count) {
          candidate_depth = prefix_excess(chunk_begin);
          candidate_position = chunk_begin;
        } else {
          const std::size_t word = chunk_begin >> 6;
          const ExcessResult result =
              excess_min_128(bp_bits_.data() + word, local_left, local_right);
          candidate_depth = prefix_excess(chunk_begin) + result.min_excess;
          candidate_position = chunk_begin + result.offset;
        }
 
        if (candidate_depth < best_depth) {
          best_depth = candidate_depth;
          best_position = candidate_position;
        }
 
        position = chunk_begin + local_right + 1;
      }
 
      return best_position;
    }

    std::size_t depth_arg_min_with_zero_prefix(std::size_t left,
                                               std::size_t right,
                                               std::size_t bit_count) const {
      const std::size_t depth_count = bit_count + 1;
      if (left >= right || right > depth_count) {
        return npos;
      }
 
      std::size_t position = left;
      int best_depth = prefix_excess_with_zero_prefix(position, bit_count);
      std::size_t best_position = position;
 
      while (position < right) {
        const std::size_t chunk_begin = (position / 128) * 128;
        const std::size_t local_left = position - chunk_begin;
        const std::size_t local_right =
            std::min<std::size_t>(right - 1, chunk_begin + 128) - chunk_begin;
 
        int candidate_depth;
        std::size_t candidate_position;
        if (chunk_begin >= bit_count) {
          candidate_depth =
              prefix_excess_with_zero_prefix(chunk_begin, bit_count);
          candidate_position = chunk_begin;
        } else {
          const std::size_t word = chunk_begin >> 6;
          const ExcessResult result =
              excess_min_128(bp_bits_.data() + word, local_left, local_right);
          candidate_depth =
              prefix_excess_with_zero_prefix(chunk_begin, bit_count) +
              result.min_excess;
          candidate_position = chunk_begin + result.offset;
        }
 
        if (candidate_depth < best_depth) {
          best_depth = candidate_depth;
          best_position = candidate_position;
        }
 
        position = chunk_begin + local_right + 1;
      }
 
      return best_position;
    }
 
    std::array<std::uint64_t, kSelectorWords> bp_bits_{};
  };
 
  static_assert(sizeof(Bp512Selector) == 64);
 
  class alignas(LeafSize == 496 ? 64 : 32) PrefixSuffixMaskLeafSelector {
   public:
    PrefixSuffixMaskLeafSelector() = default;

    template <class EntryLess>
    void build(std::size_t entry_count, EntryLess entry_less) {
      if (entry_count > kMaskEntries) {
        throw std::length_error(
            "HybridBTree prefix/suffix leaf selector too large");
      }
 
      words_.fill(0);
      if (entry_count == 0) {
        set_embedded_min_offset(0);
        return;
      }
 
      std::size_t prefix_best = 0;
      set_mask_bit(0);
      for (std::size_t slot = 1; slot < entry_count; ++slot) {
        if (entry_less(slot, prefix_best)) {
          prefix_best = slot;
          set_mask_bit(slot);
        }
      }
 
      std::size_t suffix_best = entry_count - 1;
      set_mask_bit(suffix_best);
      for (std::size_t slot = entry_count - 1; slot > 0;) {
        --slot;
        if (!entry_less(suffix_best, slot)) {
          suffix_best = slot;
          set_mask_bit(slot);
        }
      }
 
      set_embedded_min_offset(prefix_best);
    }

    std::size_t arg_min(std::size_t slot_left,
                        std::size_t slot_right,
                        std::size_t entry_count) const {
      if (slot_left >= slot_right || slot_right > entry_count ||
          entry_count > kMaskEntries) {
        return npos;
      }
      if (slot_left + 1 == slot_right) {
        return slot_left;
      }
 
      const std::size_t min_offset = embedded_min_offset();
      if (slot_left <= min_offset && min_offset < slot_right) {
        return min_offset;
      }
      if (slot_left == 0 && slot_right <= min_offset) {
        return previous_set_bit_before(slot_right);
      }
      if (slot_right == entry_count && min_offset < slot_left) {
        return next_set_bit_at_or_after(slot_left);
      }
      return npos;
    }

    void set_embedded_min_offset(std::size_t offset) {
      words_[kOffsetWord] =
          (words_[kOffsetWord] & kMaskWordMask) |
          ((static_cast<std::uint64_t>(offset) & kOffsetMask) << kOffsetShift);
    }

    std::size_t embedded_min_offset() const {
      return static_cast<std::size_t>((words_[kOffsetWord] >> kOffsetShift) &
                                      kOffsetMask);
    }
 
   private:
    static constexpr std::size_t kMaskEntries = LeafSize;
    static constexpr std::size_t kOffsetBits = LeafSize == 496 ? 16 : 8;
    static constexpr std::size_t kPackedBits = kMaskEntries + kOffsetBits;
    static constexpr std::size_t kWordCount = (kPackedBits + 63) / 64;
    static constexpr std::size_t kOffsetWord = kMaskEntries / 64;
    static constexpr std::size_t kOffsetShift = kMaskEntries & 63;

    static constexpr std::uint64_t low_bits_mask(std::size_t bits) {
      if (bits == 0) {
        return 0;
      }
      if (bits >= 64) {
        return std::numeric_limits<std::uint64_t>::max();
      }
      return (std::uint64_t{1} << bits) - 1;
    }
 
    static constexpr std::uint64_t kOffsetMask = low_bits_mask(kOffsetBits);
    static constexpr std::uint64_t kMaskWordMask = low_bits_mask(kOffsetShift);
 
    static_assert(!kMaskLeafSelector || kPackedBits % 64 == 0);
    static_assert(!kMaskLeafSelector || kOffsetShift + kOffsetBits == 64);

    void set_mask_bit(std::size_t slot) {
      words_[slot >> 6] |= std::uint64_t{1} << (slot & 63);
    }

    std::uint64_t mask_word(std::size_t word) const {
      return word == kOffsetWord ? words_[word] & kMaskWordMask : words_[word];
    }

    std::size_t previous_set_bit_before(std::size_t limit) const {
      if (limit == 0) {
        return npos;
      }
 
      std::size_t word = (limit - 1) >> 6;
      std::uint64_t bits =
          mask_word(word) & first_bits_mask(((limit - 1) & 63) + 1);
      while (true) {
        if (bits != 0) {
          return word * 64 + 63 - std::countl_zero(bits);
        }
        if (word == 0) {
          break;
        }
        --word;
        bits = mask_word(word);
      }
      return npos;
    }

    std::size_t next_set_bit_at_or_after(std::size_t slot) const {
      std::size_t word = slot >> 6;
      std::uint64_t bits = mask_word(word) & ~first_bits_mask(slot & 63);
      while (word < kWordCount) {
        if (bits != 0) {
          const std::size_t result = word * 64 + std::countr_zero(bits);
          return result < kMaskEntries ? result : npos;
        }
        ++word;
        bits = word < kWordCount ? mask_word(word) : 0;
      }
      return npos;
    }
 
    std::array<std::uint64_t, kWordCount> words_{};
  };
 
  static_assert(!kMaskLeafSelector || sizeof(PrefixSuffixMaskLeafSelector) ==
                                          (LeafSize == 496 ? 64 : 32));
  static_assert(!kMaskLeafSelector || alignof(PrefixSuffixMaskLeafSelector) ==
                                          (LeafSize == 496 ? 64 : 32));
 
  using LeafSelector = std::conditional_t<kMaskLeafSelector,
                                          PrefixSuffixMaskLeafSelector,
                                          Bp512Selector>;

  bool missing_position(std::size_t position) const {
    if constexpr (kInvalidIndexEqualsNpos) {
      return position == npos;
    } else {
      return position == npos ||
             position == static_cast<std::size_t>(invalid_index);
    }
  }

  MinCandidate value_candidate(std::size_t position) const {
    if (missing_position(position) || position >= values_.size()) {
      return {};
    }
    return {position, values_.data() + position};
  }

  MinCandidate subtree_min_candidate(std::size_t level,
                                     std::size_t node) const {
    const std::size_t position = subtree_min_position(level, node);
    if (missing_position(position)) {
      return {};
    }
    return {position, &subtree_min_value(level, node)};
  }

  MinCandidate subtree_child_min_candidate(std::size_t level,
                                           std::size_t node,
                                           std::size_t slot) const {
    const std::size_t child_level = level - 1;
    return subtree_min_candidate(child_level,
                                 node * fanout_at_level(level) + slot);
  }

  bool strictly_better_candidate(MinCandidate left, MinCandidate right) const {
    if (missing_position(left.position)) {
      return false;
    }
    if (missing_position(right.position)) {
      return true;
    }
    return compare_(*left.value, *right.value);
  }

  MinCandidate better_candidate(MinCandidate left, MinCandidate right) const {
    if (missing_position(left.position)) {
      return right;
    }
    if (missing_position(right.position)) {
      return left;
    }
    if (compare_(*right.value, *left.value)) {
      return right;
    }
    if (compare_(*left.value, *right.value)) {
      return left;
    }
    return right.position < left.position ? right : left;
  }

  bool strictly_better_subtree_child_slot(std::size_t level,
                                          std::size_t node,
                                          std::size_t left_slot,
                                          std::size_t right_slot) const {
    return strictly_better_candidate(
        subtree_child_min_candidate(level, node, left_slot),
        subtree_child_min_candidate(level, node, right_slot));
  }

  void build() {
    leaf_selectors_.clear();
    medium_selectors_.clear();
    medium_min_values_.clear();
    top_sparse_candidates_.clear();
    top_block_size_ = kMinTopSparseBlockSize;
    top_block_count_ = 0;
    top_sparse_levels_ = 0;
    medium_level_offsets_.clear();
    level_sizes_.clear();
    level_value_spans_.clear();
    level_fanouts_.clear();
    if (values_.empty()) {
      return;
    }
    if (values_.size() > static_cast<std::size_t>(invalid_index)) {
      throw std::length_error("HybridBTree index type is too small");
    }
 
    initialize_layout((values_.size() + LeafSize - 1) / LeafSize);
    for (std::size_t leaf = 0; leaf < level_sizes_[0]; ++leaf) {
      build_leaf(leaf);
    }
 
    for (std::size_t level = 1; level < level_count(); ++level) {
      for (std::size_t node = 0; node < level_sizes_[level]; ++node) {
        build_internal_node(level, node);
      }
    }
    build_top_sparse_table();
  }

  void initialize_layout(std::size_t leaf_count) {
    level_sizes_.push_back(leaf_count);
    level_value_spans_.push_back(LeafSize);
    level_fanouts_.push_back(0);
 
    std::size_t current_count = leaf_count;
    std::size_t current_span = LeafSize;
    while (current_count > 1) {
      level_fanouts_.push_back(kMiddleFanout);
      current_count = ceil_div(current_count, kMiddleFanout);
      current_span = saturating_product(current_span, kMiddleFanout);
      level_sizes_.push_back(current_count);
      level_value_spans_.push_back(current_span);
    }
 
    leaf_selectors_.resize(level_sizes_[0]);
 
    medium_level_offsets_.assign(level_count(), 0);
    if (level_count() <= 1) {
      return;
    }
 
    std::size_t medium_node_count = 0;
    for (std::size_t level = 1; level < level_count(); ++level) {
      medium_level_offsets_[level] = medium_node_count;
      medium_node_count += level_sizes_[level];
    }
    medium_selectors_.resize(medium_node_count);
    medium_min_values_.reserve(medium_node_count);
  }

  void build_leaf(std::size_t leaf) {
    LeafSelector& selector = leaf_selectors_[leaf];
    const std::size_t begin = node_value_begin(0, leaf);
    const std::size_t count = entry_count(0, leaf);
    selector.build(count, [&](std::size_t left, std::size_t right) {
      return compare_(values_[begin + left], values_[begin + right]);
    });
 
    if constexpr (kBpLeafSelector) {
      const std::size_t slot = selector.arg_min(0, count, count);
      selector.set_embedded_min_offset(slot);
    }
  }

  void build_internal_node(std::size_t level, std::size_t node) {
    Bp512Selector& selector = mutable_selector_at(level, node);
    const std::size_t count = entry_count(level, node);
    const std::size_t first_child = node * fanout_at_level(level);
 
    selector.build(count, [&](std::size_t left, std::size_t right) {
      return strictly_better_subtree_child_slot(level, node, left, right);
    });
 
    const std::size_t slot = selector.arg_min(0, count, count);
    const std::size_t min_position =
        subtree_min_position(level - 1, first_child + slot);
    selector.set_embedded_min_position(min_position);
    medium_min_values_.push_back(values_[min_position]);
    selector.build_zero_prefix_metadata(2 * count);
  }

  static std::size_t saturating_product(std::size_t left, std::size_t right) {
    if (left != 0 && right > std::numeric_limits<std::size_t>::max() / left) {
      return std::numeric_limits<std::size_t>::max();
    }
    return left * right;
  }

  static std::size_t ceil_div(std::size_t value, std::size_t divisor) {
    return value == 0 ? 0 : 1 + (value - 1) / divisor;
  }

  std::size_t tree_arg_min(std::size_t left, std::size_t right) const {
    const std::size_t root_level = level_count() - 1;
    if (left == 0 && right == values_.size()) {
      return subtree_min_position(root_level, 0);
    }
 
    const std::size_t left_leaf = leaf_for_value(left);
    const std::size_t right_leaf = leaf_for_value(right - 1);
    if (left_leaf == right_leaf) {
      return leaf_range_min(left_leaf, left, right);
    }
 
    const auto [level, node_index] = covering_node(left_leaf, right_leaf);
    return query_node(level, node_index, left, right).position;
  }

  std::size_t leaf_range_min(std::size_t leaf,
                             std::size_t left,
                             std::size_t right) const {
    if (left >= right) {
      return npos;
    }
 
    const std::size_t begin = node_value_begin(0, leaf);
    const std::size_t end = node_value_end(0, leaf);
    if (left <= begin && end <= right) {
      return subtree_min_position(0, leaf);
    }
    const std::size_t slot_left = left - begin;
    const std::size_t slot_right = right - begin;
 
    if constexpr (kMaskLeafSelector) {
      const std::size_t slot = leaf_selectors_[leaf].arg_min(
          slot_left, slot_right, entry_count(0, leaf));
      if (slot != npos) {
        return begin + slot;
      }
      return linear_range_min(left, right);
    } else {
      if (right - left <= kLeafLinearScanThreshold) {
        return linear_range_min(left, right);
      }
 
      const std::size_t slot = leaf_selectors_[leaf].arg_min(
          slot_left, slot_right, entry_count(0, leaf));
      if (slot == npos) {
        return npos;
      }
      return begin + slot;
    }
  }

  std::size_t linear_range_min(std::size_t left, std::size_t right) const {
    if (left >= right) {
      return npos;
    }
#ifdef PIXIE_AVX2_SUPPORT
    if constexpr (std::is_same_v<T, std::int64_t> &&
                  std::is_same_v<Compare, std::less<T>>) {
      if (right - left >= kLeafAvx2ScanThreshold) {
        return linear_range_min_i64_avx2(left, right);
      }
    }
#endif
    std::size_t best = left;
    for (std::size_t position = left + 1; position < right; ++position) {
      if (compare_(values_[position], values_[best])) {
        best = position;
      }
    }
    return best;
  }
 
#ifdef PIXIE_AVX2_SUPPORT
  std::size_t linear_range_min_i64_avx2(std::size_t left,
                                        std::size_t right) const {
    const std::int64_t* data = values_.data();
    std::size_t position = left;
 
    __m256i best_values =
        _mm256_loadu_si256(reinterpret_cast<const __m256i*>(data + position));
    __m256i best_positions = _mm256_set_epi64x(
        static_cast<long long>(position + 3),
        static_cast<long long>(position + 2),
        static_cast<long long>(position + 1), static_cast<long long>(position));
    position += 4;
 
    for (; position + 4 <= right; position += 4) {
      const __m256i values =
          _mm256_loadu_si256(reinterpret_cast<const __m256i*>(data + position));
      const __m256i positions =
          _mm256_set_epi64x(static_cast<long long>(position + 3),
                            static_cast<long long>(position + 2),
                            static_cast<long long>(position + 1),
                            static_cast<long long>(position));
      const __m256i take_new = _mm256_cmpgt_epi64(best_values, values);
      best_values = _mm256_blendv_epi8(best_values, values, take_new);
      best_positions = _mm256_blendv_epi8(best_positions, positions, take_new);
    }
 
    alignas(32) std::int64_t value_lanes[4];
    alignas(32) std::uint64_t position_lanes[4];
    _mm256_store_si256(reinterpret_cast<__m256i*>(value_lanes), best_values);
    _mm256_store_si256(reinterpret_cast<__m256i*>(position_lanes),
                       best_positions);
 
    std::int64_t best_value = value_lanes[0];
    std::size_t best_position = static_cast<std::size_t>(position_lanes[0]);
    for (std::size_t lane = 1; lane < 4; ++lane) {
      const std::size_t lane_position =
          static_cast<std::size_t>(position_lanes[lane]);
      if (value_lanes[lane] < best_value ||
          (value_lanes[lane] == best_value && lane_position < best_position)) {
        best_value = value_lanes[lane];
        best_position = lane_position;
      }
    }
 
    for (; position < right; ++position) {
      if (data[position] < best_value) {
        best_value = data[position];
        best_position = position;
      }
    }
    return best_position;
  }
#endif

  std::pair<std::size_t, std::size_t> covering_node(
      std::size_t left_leaf,
      std::size_t right_leaf) const {
    std::size_t level = 0;
    std::size_t left_node = left_leaf;
    std::size_t right_node = right_leaf;
    while (left_node != right_node) {
      ++level;
      const std::size_t fanout = fanout_at_level(level);
      left_node /= fanout;
      right_node /= fanout;
    }
    return {level, left_node};
  }

  std::size_t leaf_for_value(std::size_t position) const {
    return position / LeafSize;
  }

  std::size_t child_for_value(std::size_t child_level,
                              std::size_t position) const {
    return position / level_value_spans_[child_level];
  }

  MinCandidate query_child_slots(std::size_t level,
                                 std::size_t node,
                                 std::size_t slot_left,
                                 std::size_t slot_right,
                                 std::size_t left,
                                 std::size_t right) const {
    if (slot_left >= slot_right) {
      return {};
    }
 
    const std::size_t count = entry_count(level, node);
    const std::size_t slot =
        slot_left + 1 == slot_right
            ? slot_left
            : selector_arg_min(level, node, slot_left, slot_right, count);
    if (slot == npos) {
      return {};
    }
 
    const std::size_t child_level = level - 1;
    const std::size_t first_child = node * fanout_at_level(level);
    const std::size_t child = first_child + slot;
    const MinCandidate child_min = subtree_min_candidate(child_level, child);
    const std::size_t child_begin = node_value_begin(child_level, child);
    const std::size_t child_end = node_value_end(child_level, child);
    if ((left <= child_begin && child_end <= right) ||
        contains_position(left, right, child_min.position)) {
      return child_min;
    }
 
    const std::size_t last_slot = slot_right - 1;
    const std::size_t left_child_begin =
        node_value_begin(child_level, first_child + slot_left);
    const std::size_t left_child_end =
        node_value_end(child_level, first_child + slot_left);
    MinCandidate answer = query_node(child_level, first_child + slot_left,
                                     std::max(left, left_child_begin),
                                     std::min(right, left_child_end));
 
    if (slot_left != last_slot) {
      const std::size_t right_child_begin =
          node_value_begin(child_level, first_child + last_slot);
      const std::size_t right_child_end =
          node_value_end(child_level, first_child + last_slot);
      answer = better_candidate(answer,
                                query_node(child_level, first_child + last_slot,
                                           std::max(left, right_child_begin),
                                           std::min(right, right_child_end)));
    }
 
    if (slot_left + 1 < last_slot) {
      answer = better_candidate(
          answer,
          full_child_slot_range_min(level, node, slot_left + 1, last_slot));
    }
 
    return answer;
  }

  MinCandidate full_child_slot_range_min(std::size_t level,
                                         std::size_t node,
                                         std::size_t slot_left,
                                         std::size_t slot_right) const {
    if (slot_left >= slot_right) {
      return {};
    }
 
    const std::size_t slot =
        slot_left + 1 == slot_right
            ? slot_left
            : selector_arg_min(level, node, slot_left, slot_right,
                               entry_count(level, node));
    if (slot == npos) {
      return {};
    }
 
    return subtree_child_min_candidate(level, node, slot);
  }

  MinCandidate query_node(std::size_t level,
                          std::size_t node,
                          std::size_t left,
                          std::size_t right) const {
    if (left >= right) {
      return {};
    }
    const std::size_t begin = node_value_begin(level, node);
    const std::size_t end = node_value_end(level, node);
    if (left <= begin && end <= right) {
      return subtree_min_candidate(level, node);
    }
    if (level == 0) {
      return value_candidate(leaf_range_min(node, left, right));
    }
 
    const std::size_t child_level = level - 1;
    const std::size_t left_child = child_for_value(child_level, left);
    const std::size_t right_child = child_for_value(child_level, right - 1);
    const std::size_t first_child = node * fanout_at_level(level);
    const std::size_t left_slot = left_child - first_child;
    const std::size_t right_slot = right_child - first_child + 1;
    return query_child_slots(level, node, left_slot, right_slot, left, right);
  }

  bool contains_position(std::size_t left,
                         std::size_t right,
                         std::size_t position) const {
    return !missing_position(position) && left <= position && position < right;
  }

  std::size_t level_count() const { return level_sizes_.size(); }

  std::size_t entry_count(std::size_t level, std::size_t node) const {
    if (level == 0) {
      const std::size_t begin = node_value_begin(0, node);
      return std::min<std::size_t>(LeafSize, values_.size() - begin);
    }
    const std::size_t first_child = node * fanout_at_level(level);
    return std::min<std::size_t>(fanout_at_level(level),
                                 level_sizes_[level - 1] - first_child);
  }

  std::size_t node_value_begin(std::size_t level, std::size_t node) const {
    return node * level_value_spans_[level];
  }

  std::size_t node_value_end(std::size_t level, std::size_t node) const {
    return std::min(values_.size(),
                    node_value_begin(level, node) + level_value_spans_[level]);
  }

  std::size_t subtree_min_position(std::size_t level, std::size_t node) const {
    if (level == 0) {
      return node_value_begin(0, node) +
             leaf_selectors_[node].embedded_min_offset();
    }
    return selector_at(level, node).embedded_min_position();
  }

  const T& subtree_min_value(std::size_t level, std::size_t node) const {
    if (level == 0) {
      return values_[subtree_min_position(0, node)];
    }
    return medium_min_values_[medium_flat_index(level, node)];
  }

  const Bp512Selector& selector_at(std::size_t level, std::size_t node) const {
    return medium_selectors_[medium_flat_index(level, node)];
  }

  Bp512Selector& mutable_selector_at(std::size_t level, std::size_t node) {
    return medium_selectors_[medium_flat_index(level, node)];
  }

  std::size_t selector_arg_min(std::size_t level,
                               std::size_t node,
                               std::size_t slot_left,
                               std::size_t slot_right,
                               std::size_t count) const {
    const Bp512Selector& selector = selector_at(level, node);
    return selector.arg_min_with_zero_prefix(slot_left, slot_right, count);
  }

  std::size_t medium_flat_index(std::size_t level, std::size_t node) const {
    return medium_level_offsets_[level] + node;
  }

  std::size_t fanout_at_level(std::size_t level) const {
    return level_fanouts_[level];
  }

  void build_top_sparse_table() {
    top_sparse_candidates_.clear();
    top_block_size_ = top_sparse_block_size_for(values_.size());
    top_block_count_ = top_sparse_block_count_for(values_.size());
    top_sparse_levels_ =
        top_block_count_ == 0 ? 0 : std::bit_width(top_block_count_);
    if (top_block_count_ == 0) {
      return;
    }
 
    top_sparse_candidates_.assign(top_sparse_levels_ * top_block_count_,
                                  TopCandidate{});
    for (std::size_t block = 0; block < top_block_count_; ++block) {
      const std::size_t begin = block * top_block_size_;
      const std::size_t end = std::min(values_.size(), begin + top_block_size_);
      std::size_t minimum = begin;
      for (std::size_t position = begin + 1; position < end; ++position) {
        if (strictly_better_value_position(position, minimum)) {
          minimum = position;
        }
      }
      top_sparse_candidates_[block] = make_top_candidate(minimum);
    }
 
    for (std::size_t level = 1; level < top_sparse_levels_; ++level) {
      const std::size_t span = std::size_t{1} << level;
      const std::size_t half_span = span >> 1;
      TopCandidate* current =
          top_sparse_candidates_.data() + level * top_block_count_;
      const TopCandidate* previous =
          top_sparse_candidates_.data() + (level - 1) * top_block_count_;
      for (std::size_t block = 0; block + span <= top_block_count_; ++block) {
        current[block] =
            better_top_candidate(previous[block], previous[block + half_span]);
      }
    }
  }

  TopCandidate make_top_candidate(std::size_t position) const {
    if (!valid_value_position(position)) {
      return {};
    }
    return {static_cast<Index>(position)};
  }

  bool valid_value_position(std::size_t position) const {
    return !missing_position(position) && position < values_.size();
  }

  bool strictly_better_value_position(std::size_t left,
                                      std::size_t right) const {
    if (!valid_value_position(left)) {
      return false;
    }
    if (!valid_value_position(right)) {
      return true;
    }
    if (compare_(values_[left], values_[right])) {
      return true;
    }
    if (compare_(values_[right], values_[left])) {
      return false;
    }
    return left < right;
  }

  TopCandidate better_top_candidate(TopCandidate left,
                                    TopCandidate right) const {
    const std::size_t left_position = static_cast<std::size_t>(left.position);
    const std::size_t right_position = static_cast<std::size_t>(right.position);
    return strictly_better_value_position(right_position, left_position) ? right
                                                                         : left;
  }

  TopCandidate top_sparse_block_arg_min(std::size_t block_left,
                                        std::size_t block_right) const {
    if (block_left >= block_right || block_right > top_block_count_ ||
        top_sparse_levels_ == 0) {
      return {};
    }
    const std::size_t length = block_right - block_left;
    const std::size_t level = std::bit_width(length) - 1;
    const std::size_t span = std::size_t{1} << level;
    const TopCandidate* table =
        top_sparse_candidates_.data() + level * top_block_count_;
    return better_top_candidate(table[block_left], table[block_right - span]);
  }

  bool top_candidate_inside(TopCandidate candidate,
                            std::size_t left,
                            std::size_t right) const {
    const std::size_t position = static_cast<std::size_t>(candidate.position);
    return valid_value_position(position) && left <= position &&
           position < right;
  }

  std::size_t top_sparse_arg_min(std::size_t left, std::size_t right) const {
    if (top_block_count_ <= 1) {
      return npos;
    }
 
    const std::size_t padded_block_left = left / top_block_size_;
    const std::size_t padded_block_right = (right - 1) / top_block_size_ + 1;
    if (padded_block_left + 1 >= padded_block_right) {
      return npos;
    }
 
    const TopCandidate padded =
        top_sparse_block_arg_min(padded_block_left, padded_block_right);
    if (top_candidate_inside(padded, left, right)) {
      return static_cast<std::size_t>(padded.position);
    }
 
    const std::size_t first_full_block =
        (left + top_block_size_ - 1) / top_block_size_;
    const std::size_t full_block_right = right / top_block_size_;
    if (first_full_block >= full_block_right) {
      return npos;
    }
 
    TopCandidate answer =
        top_sparse_block_arg_min(first_full_block, full_block_right);
 
    const std::size_t left_border_end = first_full_block * top_block_size_;
    if (left < left_border_end) {
      answer = better_top_candidate(
          answer, make_top_candidate(tree_arg_min(left, left_border_end)));
    }
 
    const std::size_t right_border_begin = full_block_right * top_block_size_;
    if (right_border_begin < right) {
      answer = better_top_candidate(
          answer, make_top_candidate(tree_arg_min(right_border_begin, right)));
    }
 
    return valid_value_position(static_cast<std::size_t>(answer.position))
               ? static_cast<std::size_t>(answer.position)
               : npos;
  }
 
  std::span<const T> values_;
  Compare compare_;
  std::vector<LeafSelector> leaf_selectors_;
  std::vector<Bp512Selector> medium_selectors_;
  std::vector<T> medium_min_values_;
  std::vector<TopCandidate> top_sparse_candidates_;
  std::vector<std::size_t> medium_level_offsets_;
  std::vector<std::size_t> level_sizes_;
  std::vector<std::size_t> level_value_spans_;
  std::vector<std::size_t> level_fanouts_;
  std::size_t top_block_size_ = kMinTopSparseBlockSize;
  std::size_t top_block_count_ = 0;
  std::size_t top_sparse_levels_ = 0;
};
 
}  // namespace pixie::rmq