cartesian_hybrid_btree.h

#pragma once
 
#include <pixie/bits.h>
#include <pixie/detail/serialization.h>
#include <pixie/memory_usage.h>
#include <pixie/rank_select/support.h>
#include <pixie/rmq.h>
#include <pixie/rmq/utils/succinct_monotone_stack.h>
#include <pixie/storage/aligned.h>
 
#include <algorithm>
#include <array>
#include <bit>
#include <cstddef>
#include <cstdint>
#include <functional>
#include <limits>
#include <optional>
#include <span>
#include <stdexcept>
#include <type_traits>
#include <utility>
#include <vector>
 
namespace pixie::rmq {
 
template <class T,
          class Compare,
          class Index,
          std::size_t LeafSize,
          bool UseTopSparseOverlay>
class CartesianHybridBTree;
 
namespace detail {

template <class Index = std::size_t,
          std::size_t LeafSize = 512,
          bool UseHighSparseLayout = true,
          std::size_t HighSparseLayoutLevels = 2>
class HybridBTreePlusMinusOne {
  template <class, class, class, std::size_t, bool>
  friend class ::pixie::rmq::CartesianHybridBTree;
 
 public:
  static_assert(std::is_unsigned_v<Index>,
                "HybridBTreePlusMinusOne index type must be unsigned");
  static_assert(LeafSize != 0 && LeafSize % 512 == 0,
                "HybridBTreePlusMinusOne leaf size must be a positive "
                "multiple of 512");
  static_assert(!UseHighSparseLayout || HighSparseLayoutLevels > 0,
                "HybridBTreePlusMinusOne high sparse layout must cover "
                "at least one level when enabled");
 
  static constexpr std::size_t npos = std::numeric_limits<std::size_t>::max();
  static constexpr Index invalid_index = std::numeric_limits<Index>::max();
  static constexpr std::size_t kLeafSize = LeafSize;
  static constexpr std::size_t kHighLevelFanout = 256;
  static constexpr std::size_t kMiddleFanout = 192;
  static constexpr bool kUseHighSparseLayout = UseHighSparseLayout;
  static constexpr std::size_t kHighSparseLayoutLevels = HighSparseLayoutLevels;

  HybridBTreePlusMinusOne() = default;

  HybridBTreePlusMinusOne(std::span<const std::uint64_t> bits,
                          std::size_t depth_count) {
    build(bits, depth_count);
  }

  HybridBTreePlusMinusOne(std::span<const std::uint64_t> bits,
                          std::size_t depth_count,
                          const RankSelectSupport<>& rank_index) {
    build(bits, depth_count, rank_index);
  }

  void build(std::span<const std::uint64_t> bits, std::size_t depth_count) {
    input_bits_ = bits;
    depth_count_ = depth_count;
    external_rank_index_ = nullptr;
    build();
  }

  void build(std::span<const std::uint64_t> bits,
             std::size_t depth_count,
             const RankSelectSupport<>& rank_index) {
    input_bits_ = bits;
    depth_count_ = depth_count;
    external_rank_index_ = &rank_index;
    build();
  }

  std::size_t size() const { return depth_count_; }

  bool empty() const { return depth_count_ == 0; }

  std::size_t memory_usage_bytes() const {
    std::size_t bytes = sizeof(*this);
    if (external_rank_index_ == nullptr) {
      bytes += pixie::optional_nested_owned_memory_bytes(owned_rank_index_);
    }
    bytes += pixie::vector_capacity_bytes(internal_selectors_);
    bytes += pixie::vector_capacity_bytes(internal_min_positions_);
    bytes += pixie::vector_capacity_bytes(internal_min_depths_);
    bytes += pixie::vector_capacity_bytes(high_child_metadata_);
    bytes += pixie::vector_capacity_bytes(high_sparse_min_slots_);
    bytes += pixie::vector_capacity_bytes(internal_level_offsets_);
    bytes += pixie::vector_capacity_bytes(min_summary_level_offsets_);
    bytes += pixie::vector_capacity_bytes(high_level_offsets_);
    bytes += pixie::vector_capacity_bytes(level_sizes_);
    bytes += pixie::vector_capacity_bytes(level_position_spans_);
    bytes += pixie::vector_capacity_bytes(level_fanouts_);
    return bytes;
  }

  std::size_t arg_min(std::size_t left, std::size_t right) const {
    if (left >= right || right > depth_count_ || level_sizes_.empty()) {
      return npos;
    }
    if (left + 1 == right) {
      return left;
    }
 
    const std::size_t root_level = level_count() - 1;
    if (left == 0 && right == depth_count_) {
      return subtree_min_candidate(root_level, 0).position;
    }
 
    const std::size_t left_leaf = leaf_for_position(left);
    const std::size_t right_leaf = leaf_for_position(right - 1);
    if (left_leaf == right_leaf) {
      return leaf_range_arg_min_relative(left_leaf, left, right);
    }
 
    const auto [level, node] = covering_node(left_leaf, right_leaf);
    return query_node(level, node, left, right).position;
  }

  std::size_t select0(std::size_t rank) const {
    if (rank == 0 || depth_count_ == 0 || rank_index_or_null() == nullptr) {
      return npos;
    }
    const std::size_t delta_count = depth_count_ - 1;
    const std::size_t position = rank_index().select0(rank);
    return position < delta_count ? position : npos;
  }
 
 private:
  static constexpr std::size_t kLeafWords = LeafSize / 64;
  static constexpr std::size_t kLeafChunks = LeafSize / 128;
  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 kEmbeddedSummaryWords = 2;
  static constexpr std::size_t kEmbeddedSummaryBits =
      64 * kEmbeddedSummaryWords;
  static constexpr std::size_t kEmbeddedSummaryMaxEntries =
      (kSelectorBits - kEmbeddedSummaryBits) / 2;
  static constexpr std::size_t kEmbeddedSummaryPositionWord =
      kSelectorWords - 2;
  static constexpr std::size_t kEmbeddedSummaryDepthWord = kSelectorWords - 1;
  static constexpr std::size_t kHighSparseTableLevels =
      static_cast<std::size_t>(std::bit_width(kHighLevelFanout));
  static constexpr std::size_t kHighSparseSlotsPerNode =
      kHighSparseTableLevels * kHighLevelFanout;
  static constexpr bool kInvalidIndexEqualsNpos =
      static_cast<std::size_t>(invalid_index) == npos;
  static_assert(kEmbeddedSummaryMaxEntries == 192);
  static_assert(kMiddleFanout == kEmbeddedSummaryMaxEntries);
  static_assert(sizeof(std::size_t) <= sizeof(std::uint64_t));
 
  struct DepthCandidate {
    std::size_t position = npos;
    std::int64_t depth = std::numeric_limits<std::int64_t>::max();
  };
 
  struct HighChildMetadata {
    std::size_t position_begin = 0;
    std::size_t position_end = 0;
    Index min_position = invalid_index;
    std::int64_t min_depth = std::numeric_limits<std::int64_t>::max();
  };
 
  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(
            "HybridBTreePlusMinusOne 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;
    }

    void set_embedded_min_summary(std::size_t position, std::int64_t depth) {
      bp_bits_[kEmbeddedSummaryPositionWord] =
          static_cast<std::uint64_t>(position);
      bp_bits_[kEmbeddedSummaryDepthWord] = std::bit_cast<std::uint64_t>(depth);
    }

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

    std::int64_t embedded_min_depth() const {
      return std::bit_cast<std::int64_t>(bp_bits_[kEmbeddedSummaryDepthWord]);
    }
 
   private:
    friend class HybridBTreePlusMinusOne;

    void 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);
      }
    }

    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);
    }

    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 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(bit_count);
          candidate_position = bit_count;
        } 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::array<std::uint64_t, kSelectorWords> bp_bits_{};
  };
 
  static_assert(sizeof(Bp512Selector) == 64);
 
  void serialize_metadata(BinaryWriter& writer) const {
    if (depth_count_ != 0 &&
        (owned_rank_index_.has_value() || external_rank_index_ == nullptr)) {
      throw std::invalid_argument(
          "Depth RMQ serialization requires external rank support");
    }
    if (depth_count_ != 0) {
      validate_serialized_state(*external_rank_index_,
                                DeserializationValidation::kQuick);
    } else {
      validate_serialized_state(RankSelectSupport<>(),
                                DeserializationValidation::kQuick);
    }
 
    writer.write_size(depth_count_);
    writer.write_size(internal_selectors_.size());
    for (const Bp512Selector& selector : internal_selectors_) {
      for (const std::uint64_t word : selector.bp_bits_) {
        writer.write_u64(word);
      }
    }
    pixie::detail::write_vector(
        writer, std::span<const Index>(internal_min_positions_));
    pixie::detail::write_vector(
        writer, std::span<const std::int64_t>(internal_min_depths_));
 
    writer.write_size(high_child_metadata_.size());
    for (const HighChildMetadata& metadata : high_child_metadata_) {
      writer.write_size(metadata.position_begin);
      writer.write_size(metadata.position_end);
      pixie::detail::write_integral(writer, metadata.min_position);
      writer.write_i64(metadata.min_depth);
    }
    pixie::detail::write_vector(
        writer, std::span<const std::uint8_t>(high_sparse_min_slots_));
    pixie::detail::write_vector(
        writer, std::span<const std::size_t>(internal_level_offsets_));
    pixie::detail::write_vector(
        writer, std::span<const std::size_t>(min_summary_level_offsets_));
    pixie::detail::write_vector(
        writer, std::span<const std::size_t>(high_level_offsets_));
    pixie::detail::write_vector(writer,
                                std::span<const std::size_t>(level_sizes_));
    pixie::detail::write_vector(
        writer, std::span<const std::size_t>(level_position_spans_));
    pixie::detail::write_vector(writer,
                                std::span<const std::size_t>(level_fanouts_));
    writer.write_size(high_level_begin_);
  }
 
  static HybridBTreePlusMinusOne deserialize_metadata(
      std::span<const std::uint64_t> bits,
      const RankSelectSupport<>& rank_index,
      BinaryReader& reader) {
    HybridBTreePlusMinusOne result;
    result.input_bits_ = bits;
    result.depth_count_ = reader.read_size();
    result.external_rank_index_ = &rank_index;
 
    const std::size_t selector_count = reader.read_size();
    const std::vector<Bp512Selector> empty_selectors;
    if (selector_count > empty_selectors.max_size()) {
      throw std::length_error(
          "Serialized depth RMQ selector count is too large");
    }
    if (selector_count >
        reader.remaining() / (kSelectorWords * sizeof(std::uint64_t))) {
      throw std::invalid_argument("Truncated serialized depth RMQ selectors");
    }
    result.internal_selectors_.resize(selector_count);
    for (Bp512Selector& selector : result.internal_selectors_) {
      for (std::uint64_t& word : selector.bp_bits_) {
        word = reader.read_u64();
      }
    }
    result.internal_min_positions_ = pixie::detail::read_vector<Index>(reader);
    result.internal_min_depths_ =
        pixie::detail::read_vector<std::int64_t>(reader);
 
    const std::size_t metadata_count = reader.read_size();
    const std::vector<HighChildMetadata> empty_metadata;
    if (metadata_count > empty_metadata.max_size()) {
      throw std::length_error(
          "Serialized depth RMQ child metadata is too large");
    }
    constexpr std::size_t kSerializedMetadataBytes =
        2 * sizeof(std::uint64_t) + sizeof(Index) + sizeof(std::int64_t);
    if (metadata_count > reader.remaining() / kSerializedMetadataBytes) {
      throw std::invalid_argument(
          "Truncated serialized depth RMQ child metadata");
    }
    result.high_child_metadata_.resize(metadata_count);
    for (HighChildMetadata& metadata : result.high_child_metadata_) {
      metadata.position_begin = reader.read_size();
      metadata.position_end = reader.read_size();
      metadata.min_position = pixie::detail::read_integral<Index>(reader);
      metadata.min_depth = reader.read_i64();
    }
    result.high_sparse_min_slots_ =
        pixie::detail::read_vector<std::uint8_t>(reader);
    result.internal_level_offsets_ =
        pixie::detail::read_vector<std::size_t>(reader);
    result.min_summary_level_offsets_ =
        pixie::detail::read_vector<std::size_t>(reader);
    result.high_level_offsets_ =
        pixie::detail::read_vector<std::size_t>(reader);
    result.level_sizes_ = pixie::detail::read_vector<std::size_t>(reader);
    result.level_position_spans_ =
        pixie::detail::read_vector<std::size_t>(reader);
    result.level_fanouts_ = pixie::detail::read_vector<std::size_t>(reader);
    result.high_level_begin_ = reader.read_size();
    return result;
  }
 
  void restore_external_sources(
      std::span<const std::uint64_t> bits,
      const RankSelectSupport<>& rank_index) noexcept {
    input_bits_ = bits;
    owned_rank_index_.reset();
    external_rank_index_ = &rank_index;
  }
 
  void restore_empty_sources() noexcept {
    input_bits_ = {};
    owned_rank_index_.reset();
    external_rank_index_ = nullptr;
  }
 
  void validate_serialized_state(const RankSelectSupport<>& rank_index,
                                 DeserializationValidation validation) const {
    if (depth_count_ == 0) {
      if (!input_bits_.empty() || owned_rank_index_.has_value() ||
          !internal_selectors_.empty() || !internal_min_positions_.empty() ||
          !internal_min_depths_.empty() || !high_child_metadata_.empty() ||
          !high_sparse_min_slots_.empty() || !internal_level_offsets_.empty() ||
          !min_summary_level_offsets_.empty() || !high_level_offsets_.empty() ||
          !level_sizes_.empty() || !level_position_spans_.empty() ||
          !level_fanouts_.empty() ||
          high_level_begin_ != std::numeric_limits<std::size_t>::max()) {
        throw std::invalid_argument(
            "Invalid serialized empty depth RMQ metadata");
      }
      return;
    }
 
    const std::size_t delta_count = depth_count_ - 1;
    const std::size_t required_words =
        delta_count == 0 ? 0 : 1 + (delta_count - 1) / 64;
    if (required_words > input_bits_.size() ||
        rank_index.size() < delta_count || owned_rank_index_.has_value()) {
      throw std::invalid_argument(
          "Invalid serialized depth RMQ source metadata");
    }
 
    const std::size_t level_count = level_sizes_.size();
    if (level_count == 0 || internal_level_offsets_.size() != level_count ||
        min_summary_level_offsets_.size() != level_count ||
        high_level_offsets_.size() != level_count ||
        level_position_spans_.size() != level_count ||
        level_fanouts_.size() != level_count) {
      throw std::invalid_argument(
          "Invalid serialized depth RMQ topology sizes");
    }
 
    const std::size_t expected_leaf_count = 1 + (depth_count_ - 1) / LeafSize;
    if (level_sizes_[0] != expected_leaf_count ||
        level_position_spans_[0] != LeafSize || level_fanouts_[0] != 0 ||
        internal_level_offsets_[0] != 0 ||
        min_summary_level_offsets_[0] != npos || high_level_offsets_[0] != 0) {
      throw std::invalid_argument("Invalid serialized depth RMQ leaf topology");
    }
 
    std::size_t current_count = expected_leaf_count;
    std::size_t current_span = LeafSize;
    std::size_t expected_internal_count = 0;
    for (std::size_t level = 1; level < level_count; ++level) {
      const std::size_t expected_fanout =
          current_count > kHighLevelFanout * kHighLevelFanout
              ? kMiddleFanout
              : kHighLevelFanout;
      current_count = ceil_div(current_count, expected_fanout);
      current_span = saturating_product(current_span, expected_fanout);
      if (level_fanouts_[level] != expected_fanout ||
          level_sizes_[level] != current_count ||
          level_position_spans_[level] != current_span ||
          internal_level_offsets_[level] != expected_internal_count) {
        throw std::invalid_argument(
            "Invalid serialized depth RMQ level topology");
      }
      if (current_count >
          std::numeric_limits<std::size_t>::max() - expected_internal_count) {
        throw std::length_error(
            "Serialized depth RMQ internal count is too large");
      }
      expected_internal_count += current_count;
    }
    if (current_count != 1 ||
        internal_selectors_.size() != expected_internal_count) {
      throw std::invalid_argument(
          "Invalid serialized depth RMQ selector topology");
    }
 
    const std::size_t root_level = level_count - 1;
    const std::size_t expected_high_begin =
        level_count == 1
            ? std::numeric_limits<std::size_t>::max()
            : (level_fanouts_[root_level] == kHighLevelFanout ? root_level
                                                              : level_count);
    if (high_level_begin_ != expected_high_begin) {
      throw std::invalid_argument(
          "Invalid serialized depth RMQ high-level topology");
    }
 
    std::size_t high_node_count = 0;
    std::size_t side_summary_count = 0;
    for (std::size_t level = 1; level < level_count; ++level) {
      const bool high_level = level >= expected_high_begin;
      if (high_level_offsets_[level] != high_node_count) {
        throw std::invalid_argument(
            "Invalid serialized depth RMQ high-level offsets");
      }
      if (high_level) {
        high_node_count += level_sizes_[level];
      }
      const bool embeds =
          !high_level && level_fanouts_[level] <= kEmbeddedSummaryMaxEntries;
      const std::size_t expected_summary_offset =
          embeds ? npos : side_summary_count;
      if (min_summary_level_offsets_[level] != expected_summary_offset) {
        throw std::invalid_argument(
            "Invalid serialized depth RMQ summary offsets");
      }
      if (!embeds) {
        side_summary_count += level_sizes_[level];
      }
    }
    if (internal_min_positions_.size() != side_summary_count ||
        internal_min_depths_.size() != side_summary_count ||
        high_node_count >
            std::numeric_limits<std::size_t>::max() / kHighLevelFanout ||
        high_child_metadata_.size() != high_node_count * kHighLevelFanout ||
        high_node_count >
            std::numeric_limits<std::size_t>::max() / kHighSparseSlotsPerNode ||
        high_sparse_min_slots_.size() !=
            high_node_count * kHighSparseSlotsPerNode) {
      throw std::invalid_argument(
          "Invalid serialized depth RMQ metadata counts");
    }
 
    for (const Index position : internal_min_positions_) {
      if (position == invalid_index ||
          static_cast<std::size_t>(position) >= depth_count_) {
        throw std::invalid_argument(
            "Invalid serialized depth RMQ minimum position");
      }
    }
    for (std::size_t level = expected_high_begin; level < level_count;
         ++level) {
      for (std::size_t node = 0; node < level_sizes_[level]; ++node) {
        const std::size_t count = entry_count(level, node);
        const std::size_t first_child = node * level_fanouts_[level];
        const std::size_t flat = high_level_offsets_[level] + node;
        for (std::size_t slot = 0; slot < count; ++slot) {
          const HighChildMetadata& metadata =
              high_child_metadata_[flat * kHighLevelFanout + slot];
          const std::size_t child = first_child + slot;
          const std::size_t expected_begin =
              child * level_position_spans_[level - 1];
          const std::size_t expected_end = std::min(
              depth_count_, expected_begin + level_position_spans_[level - 1]);
          if (metadata.position_begin != expected_begin ||
              metadata.position_end != expected_end ||
              metadata.min_position == invalid_index ||
              metadata.min_position < expected_begin ||
              metadata.min_position >= expected_end) {
            throw std::invalid_argument(
                "Invalid serialized depth RMQ child metadata");
          }
        }
      }
    }
    if (validation == DeserializationValidation::kFull) {
      validate_exact_metadata();
    }
  }
 
  void validate_exact_metadata() const {
    for (std::size_t level = 1; level < level_count(); ++level) {
      for (std::size_t node = 0; node < level_sizes_[level]; ++node) {
        const std::size_t count = entry_count(level, node);
        const std::size_t first_child = node * level_fanouts_[level];
        std::array<DepthCandidate, kSelectorEntries> child_minima{};
        for (std::size_t slot = 0; slot < count; ++slot) {
          child_minima[slot] =
              subtree_min_candidate(level - 1, first_child + slot);
        }
 
        Bp512Selector expected_selector;
        expected_selector.build(count,
                                [&](std::size_t left, std::size_t right) {
                                  return strictly_better_candidate(
                                      child_minima[left], child_minima[right]);
                                });
        std::size_t best_slot = 0;
        for (std::size_t slot = 1; slot < count; ++slot) {
          if (strictly_better_candidate(child_minima[slot],
                                        child_minima[best_slot])) {
            best_slot = slot;
          }
        }
        const DepthCandidate expected_minimum = child_minima[best_slot];
        if (level_embeds_min_summary(level)) {
          expected_selector.set_embedded_min_summary(expected_minimum.position,
                                                     expected_minimum.depth);
        } else {
          const std::size_t flat = min_summary_flat_index(level, node);
          if (internal_min_positions_[flat] != expected_minimum.position ||
              internal_min_depths_[flat] != expected_minimum.depth) {
            throw std::invalid_argument(
                "Serialized depth RMQ minimum disagrees with source");
          }
        }
        if (selector_at(level, node).bp_bits_ != expected_selector.bp_bits_) {
          throw std::invalid_argument(
              "Serialized depth RMQ selector disagrees with source");
        }
 
        if (!is_high_level(level)) {
          continue;
        }
        const std::size_t high_flat = high_flat_index(level, node);
        for (std::size_t slot = 0; slot < kHighLevelFanout; ++slot) {
          HighChildMetadata expected;
          if (slot < count) {
            const std::size_t child = first_child + slot;
            expected.position_begin = node_position_begin(level - 1, child);
            expected.position_end = node_position_end(level - 1, child);
            expected.min_position =
                static_cast<Index>(child_minima[slot].position);
            expected.min_depth = child_minima[slot].depth;
          }
          const HighChildMetadata& actual =
              high_child_metadata_at(high_flat, slot);
          if (actual.position_begin != expected.position_begin ||
              actual.position_end != expected.position_end ||
              actual.min_position != expected.min_position ||
              actual.min_depth != expected.min_depth) {
            throw std::invalid_argument(
                "Serialized depth RMQ child metadata disagrees with source");
          }
        }
 
        std::array<std::uint8_t, kHighSparseSlotsPerNode> expected_sparse{};
        for (std::size_t slot = 0; slot < count; ++slot) {
          expected_sparse[slot] = static_cast<std::uint8_t>(slot);
        }
        for (std::size_t sparse_level = 1;
             sparse_level < kHighSparseTableLevels; ++sparse_level) {
          const std::size_t span = std::size_t{1} << sparse_level;
          if (span > count) {
            break;
          }
          const std::size_t half_span = span >> 1;
          const std::size_t previous = (sparse_level - 1) * kHighLevelFanout;
          const std::size_t current = sparse_level * kHighLevelFanout;
          for (std::size_t slot = 0; slot + span <= count; ++slot) {
            const std::size_t left = expected_sparse[previous + slot];
            const std::size_t right =
                expected_sparse[previous + slot + half_span];
            expected_sparse[current + slot] = static_cast<std::uint8_t>(
                strictly_better_candidate(child_minima[right],
                                          child_minima[left])
                    ? right
                    : left);
          }
        }
        const std::uint8_t* actual_sparse =
            high_sparse_min_slots_begin(high_flat);
        if (!std::ranges::equal(
                expected_sparse,
                std::span(actual_sparse, kHighSparseSlotsPerNode))) {
          throw std::invalid_argument(
              "Serialized depth RMQ sparse metadata disagrees with source");
        }
      }
    }
  }

  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);
    }
  }

  DepthCandidate better_candidate(DepthCandidate left,
                                  DepthCandidate right) const {
    if (missing_position(left.position)) {
      return right;
    }
    if (missing_position(right.position)) {
      return left;
    }
    if (right.depth < left.depth) {
      return right;
    }
    if (left.depth < right.depth) {
      return left;
    }
    return right.position < left.position ? right : left;
  }

  bool strictly_better_candidate(DepthCandidate left,
                                 DepthCandidate right) const {
    if (missing_position(left.position)) {
      return false;
    }
    if (missing_position(right.position)) {
      return true;
    }
    if (left.depth != right.depth) {
      return left.depth < right.depth;
    }
    return left.position < right.position;
  }

  void build() {
    owned_rank_index_.reset();
    internal_selectors_.clear();
    internal_min_positions_.clear();
    internal_min_depths_.clear();
    high_child_metadata_.clear();
    high_sparse_min_slots_.clear();
    internal_level_offsets_.clear();
    min_summary_level_offsets_.clear();
    high_level_offsets_.clear();
    level_sizes_.clear();
    level_position_spans_.clear();
    level_fanouts_.clear();
    high_level_begin_ = std::numeric_limits<std::size_t>::max();
 
    if (depth_count_ == 0) {
      return;
    }
    if (depth_count_ > static_cast<std::size_t>(invalid_index)) {
      throw std::length_error(
          "HybridBTreePlusMinusOne index type is too small");
    }
    if (depth_count_ >
        static_cast<std::size_t>(std::numeric_limits<std::int64_t>::max())) {
      throw std::length_error(
          "HybridBTreePlusMinusOne depth range is too large");
    }
    if (input_bits_.size() < (depth_count_ - 1 + 63) / 64) {
      throw std::invalid_argument(
          "HybridBTreePlusMinusOne bit span is too small");
    }
 
    const std::size_t delta_count = depth_count_ - 1;
    if (external_rank_index_ == nullptr) {
      owned_rank_index_.emplace(input_bits_, delta_count,
                                RankSelectSupport<>::SelectSupport::kSelect0);
    } else if (external_rank_index_->size() < delta_count) {
      throw std::invalid_argument(
          "HybridBTreePlusMinusOne external rank index is too small");
    }
 
    initialize_layout((depth_count_ + LeafSize - 1) / LeafSize);
    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);
      }
    }
  }

  void initialize_layout(std::size_t leaf_count) {
    level_sizes_.push_back(leaf_count);
    level_position_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 > kHighLevelFanout * kHighLevelFanout) {
      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_position_spans_.push_back(current_span);
    }
    while (current_count > 1) {
      level_fanouts_.push_back(kHighLevelFanout);
      current_count = ceil_div(current_count, kHighLevelFanout);
      current_span = saturating_product(current_span, kHighLevelFanout);
      level_sizes_.push_back(current_count);
      level_position_spans_.push_back(current_span);
    }
    internal_level_offsets_.assign(level_count(), 0);
    min_summary_level_offsets_.assign(level_count(), npos);
    high_level_offsets_.assign(level_count(), 0);
    if (level_count() <= 1) {
      return;
    }
 
    std::size_t internal_count = 0;
    for (std::size_t level = 1; level < level_count(); ++level) {
      internal_level_offsets_[level] = internal_count;
      internal_count += level_sizes_[level];
    }
    internal_selectors_.resize(internal_count);
 
    if constexpr (UseHighSparseLayout) {
      const std::size_t root_level = level_count() - 1;
      std::size_t high_layout_levels = 0;
      for (std::size_t level = root_level;
           level > 0 && high_layout_levels < HighSparseLayoutLevels &&
           fanout_at_level(level) == kHighLevelFanout;
           --level) {
        ++high_layout_levels;
      }
      high_level_begin_ = high_layout_levels == 0
                              ? level_count()
                              : root_level + 1 - high_layout_levels;
 
      std::size_t high_node_count = 0;
      for (std::size_t level = high_level_begin_; level < level_count();
           ++level) {
        high_level_offsets_[level] = high_node_count;
        high_node_count += level_sizes_[level];
      }
      high_child_metadata_.resize(high_node_count * kHighLevelFanout);
      high_sparse_min_slots_.resize(high_node_count * kHighSparseSlotsPerNode);
    } else {
      high_level_begin_ = level_count();
    }
 
    std::size_t side_summary_count = 0;
    for (std::size_t level = 1; level < level_count(); ++level) {
      if (!level_embeds_min_summary(level)) {
        min_summary_level_offsets_[level] = side_summary_count;
        side_summary_count += level_sizes_[level];
      }
    }
    internal_min_positions_.resize(side_summary_count, invalid_index);
    internal_min_depths_.resize(side_summary_count,
                                std::numeric_limits<std::int64_t>::max());
  }

  void build_internal_node(std::size_t level, std::size_t node) {
    const std::size_t count = entry_count(level, node);
    const std::size_t first_child = node * fanout_at_level(level);
    const bool high_level = is_high_level(level);
    const std::size_t high_flat = high_level ? high_flat_index(level, node) : 0;
 
    std::array<DepthCandidate, kSelectorEntries> child_minima{};
    for (std::size_t slot = 0; slot < count; ++slot) {
      child_minima[slot] = subtree_min_candidate(level - 1, first_child + slot);
    }
 
    if (high_level) {
      for (std::size_t slot = 0; slot < count; ++slot) {
        const std::size_t child = first_child + slot;
        const DepthCandidate child_min = child_minima[slot];
        HighChildMetadata& metadata =
            mutable_high_child_metadata_at(high_flat, slot);
        metadata.position_begin = node_position_begin(level - 1, child);
        metadata.position_end = node_position_end(level - 1, child);
        metadata.min_position = static_cast<Index>(child_min.position);
        metadata.min_depth = child_min.depth;
      }
      build_high_sparse_min_slots(level, node, count);
    }
 
    Bp512Selector& selector = mutable_selector_at(level, node);
    selector.build(count, [&](std::size_t left, std::size_t right) {
      return strictly_better_candidate(child_minima[left], child_minima[right]);
    });
 
    const std::size_t slot = selector.arg_min(0, count, count);
    const DepthCandidate minimum = child_minima[slot];
    if (level_embeds_min_summary(level)) {
      selector.set_embedded_min_summary(minimum.position, minimum.depth);
    } else {
      const std::size_t flat = min_summary_flat_index(level, node);
      internal_min_positions_[flat] = static_cast<Index>(minimum.position);
      internal_min_depths_[flat] = minimum.depth;
    }
  }

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

  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;
  }

  std::uint64_t word_or_zero(std::size_t word) const {
    return word < input_bits_.size() ? input_bits_[word] : 0;
  }

  const std::uint64_t* chunk_words_or_copy(
      std::size_t first_word,
      std::array<std::uint64_t, 2>& storage) const {
    if (first_word + 1 < input_bits_.size()) {
      return input_bits_.data() + first_word;
    }
    storage[0] = word_or_zero(first_word);
    storage[1] = word_or_zero(first_word + 1);
    return storage.data();
  }

  const RankSelectSupport<>* rank_index_or_null() const {
    return external_rank_index_ != nullptr
               ? external_rank_index_
               : (owned_rank_index_ ? &*owned_rank_index_ : nullptr);
  }

  const RankSelectSupport<>& rank_index() const {
    return *rank_index_or_null();
  }

  std::int64_t depth_at_position(std::size_t position) const {
    const std::size_t delta_count = depth_count_ == 0 ? 0 : depth_count_ - 1;
    position = std::min(position, delta_count);
    const std::uint64_t ones = rank_index().rank(position);
    return static_cast<std::int64_t>(ones) -
           static_cast<std::int64_t>(position - ones);
  }

  DepthCandidate scan_leaf_range_with_base(std::size_t leaf,
                                           std::size_t left_offset,
                                           std::size_t right_offset,
                                           std::int64_t base_depth) const {
    const std::size_t begin = node_position_begin(0, leaf);
    const std::size_t count = entry_count(0, leaf);
    if (count == 0 || left_offset > right_offset || left_offset >= count) {
      return {};
    }
    right_offset = std::min(right_offset, count - 1);
 
    DepthCandidate answer;
    std::int64_t chunk_base_excess = 0;
    const std::size_t first_word = leaf * kLeafWords;
    for (std::size_t chunk = 0; chunk < kLeafChunks; ++chunk) {
      const std::size_t chunk_begin = chunk * 128;
      if (chunk_begin >= count || chunk_begin > right_offset) {
        break;
      }
 
      std::array<std::uint64_t, 2> chunk_storage{};
      const std::uint64_t* chunk_words =
          chunk_words_or_copy(first_word + 2 * chunk, chunk_storage);
 
      const std::size_t chunk_end =
          std::min<std::size_t>(count - 1, chunk_begin + 127);
      if (left_offset > chunk_end) {
        chunk_base_excess += prefix_excess_128(chunk_words, 128);
        continue;
      }
 
      const std::size_t local_left =
          std::max(left_offset, chunk_begin) - chunk_begin;
      const std::size_t local_right =
          std::min(right_offset, chunk_end) - chunk_begin;
      const ExcessResult result =
          excess_min_128(chunk_words, local_left, local_right);
      const std::size_t offset = chunk_begin + result.offset;
      if (result.offset != npos && offset < count) {
        answer = better_candidate(
            answer, {begin + offset,
                     base_depth + chunk_base_excess + result.min_excess});
      }
 
      chunk_base_excess += prefix_excess_128(chunk_words, 128);
    }
    return answer;
  }

  std::size_t leaf_range_arg_min_relative(std::size_t leaf,
                                          std::size_t left,
                                          std::size_t right) const {
    if (left >= right) {
      return npos;
    }
 
    const std::size_t begin = node_position_begin(0, leaf);
    const std::size_t count = entry_count(0, leaf);
    std::size_t left_offset = left - begin;
    std::size_t right_offset = right - begin - 1;
    if (count == 0 || left_offset >= count || left_offset > right_offset) {
      return npos;
    }
    right_offset = std::min(right_offset, count - 1);
 
    std::size_t best_position = npos;
    std::int64_t best_depth = std::numeric_limits<std::int64_t>::max();
    std::int64_t chunk_base_excess = 0;
    bool first_chunk = true;
 
    const std::size_t first_word = leaf * kLeafWords;
    std::size_t chunk = left_offset / 128;
    for (; chunk < kLeafChunks; ++chunk) {
      const std::size_t chunk_begin = chunk * 128;
      if (chunk_begin >= count || chunk_begin > right_offset) {
        break;
      }
 
      std::array<std::uint64_t, 2> chunk_storage{};
      const std::uint64_t* chunk_words =
          chunk_words_or_copy(first_word + 2 * chunk, chunk_storage);
 
      const std::size_t chunk_end =
          std::min<std::size_t>(count - 1, chunk_begin + 127);
      const std::size_t local_left =
          std::max(left_offset, chunk_begin) - chunk_begin;
      const std::size_t local_right =
          std::min(right_offset, chunk_end) - chunk_begin;
      const int left_prefix =
          first_chunk ? prefix_excess_128(chunk_words, local_left) : 0;
      const std::int64_t local_base =
          first_chunk ? -static_cast<std::int64_t>(left_prefix)
                      : chunk_base_excess;
      const ExcessResult result =
          excess_min_128(chunk_words, local_left, local_right);
      const std::size_t offset = chunk_begin + result.offset;
      if (result.offset != npos && offset < count) {
        const std::int64_t candidate_depth = local_base + result.min_excess;
        if (best_position == npos || candidate_depth < best_depth) {
          best_position = begin + offset;
          best_depth = candidate_depth;
        }
      }
 
      const int chunk_excess = prefix_excess_128(chunk_words, 128);
      if (first_chunk) {
        chunk_base_excess =
            static_cast<std::int64_t>(chunk_excess) - left_prefix;
        first_chunk = false;
      } else {
        chunk_base_excess += chunk_excess;
      }
    }
 
    return best_position;
  }

  DepthCandidate leaf_range_min(std::size_t leaf,
                                std::size_t left,
                                std::size_t right) const {
    if (left >= right) {
      return {};
    }
 
    const std::size_t begin = node_position_begin(0, leaf);
    const std::size_t slot_left = left - begin;
    const std::size_t slot_right = right - begin;
    return scan_leaf_range_with_base(leaf, slot_left, slot_right - 1,
                                     depth_at_position(begin));
  }

  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_position(std::size_t position) const {
    return position / LeafSize;
  }

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

  DepthCandidate 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 HighChildMetadata* high_children =
        is_high_level(level) ? high_child_metadata_begin(level, node) : nullptr;
    const DepthCandidate child_min =
        high_children != nullptr ? high_child_min_candidate(level, node, slot)
                                 : subtree_min_candidate(child_level, child);
    const std::size_t child_begin =
        high_children != nullptr ? high_children[slot].position_begin
                                 : node_position_begin(child_level, child);
    const std::size_t child_end = high_children != nullptr
                                      ? high_children[slot].position_end
                                      : node_position_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 =
        high_children != nullptr
            ? high_children[slot_left].position_begin
            : node_position_begin(child_level, first_child + slot_left);
    const std::size_t left_child_end =
        high_children != nullptr
            ? high_children[slot_left].position_end
            : node_position_end(child_level, first_child + slot_left);
    DepthCandidate 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 =
          high_children != nullptr
              ? high_children[last_slot].position_begin
              : node_position_begin(child_level, first_child + last_slot);
      const std::size_t right_child_end =
          high_children != nullptr
              ? high_children[last_slot].position_end
              : node_position_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;
  }

  DepthCandidate 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 child_min_candidate(level, node, slot);
  }

  DepthCandidate 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_position_begin(level, node);
    const std::size_t end = node_position_end(level, node);
    if (left <= begin && end <= right) {
      return subtree_min_candidate(level, node);
    }
    if (level == 0) {
      return leaf_range_min(node, left, right);
    }
 
    const std::size_t child_level = level - 1;
    const std::size_t left_child = child_for_position(child_level, left);
    const std::size_t right_child = child_for_position(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_position_begin(0, node);
      return std::min<std::size_t>(LeafSize, depth_count_ - 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_position_begin(std::size_t level, std::size_t node) const {
    return node * level_position_spans_[level];
  }

  std::size_t node_position_end(std::size_t level, std::size_t node) const {
    return std::min(depth_count_, node_position_begin(level, node) +
                                      level_position_spans_[level]);
  }

  DepthCandidate subtree_min_candidate(std::size_t level,
                                       std::size_t node) const {
    if (level == 0) {
      return leaf_range_min(node, node_position_begin(0, node),
                            node_position_end(0, node));
    }
    if (level_embeds_min_summary(level)) {
      const Bp512Selector& selector = selector_at(level, node);
      return {selector.embedded_min_position(), selector.embedded_min_depth()};
    }
    const std::size_t flat = min_summary_flat_index(level, node);
    return {static_cast<std::size_t>(internal_min_positions_[flat]),
            internal_min_depths_[flat]};
  }

  DepthCandidate child_min_candidate(std::size_t level,
                                     std::size_t node,
                                     std::size_t slot) const {
    if (is_high_level(level)) {
      return high_child_min_candidate(level, node, slot);
    }
    return subtree_min_candidate(level - 1,
                                 node * fanout_at_level(level) + slot);
  }

  DepthCandidate high_child_min_candidate(std::size_t level,
                                          std::size_t node,
                                          std::size_t slot) const {
    const HighChildMetadata& metadata =
        high_child_metadata_at(high_flat_index(level, node), slot);
    return {static_cast<std::size_t>(metadata.min_position),
            metadata.min_depth};
  }

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

  Bp512Selector& mutable_selector_at(std::size_t level, std::size_t node) {
    return internal_selectors_[internal_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 {
    if constexpr (UseHighSparseLayout) {
      if (is_high_level(level)) {
        return high_sparse_arg_min(level, node, slot_left, slot_right, count);
      }
    }
    return selector_at(level, node).arg_min(slot_left, slot_right, count);
  }

  bool is_high_level(std::size_t level) const {
    if constexpr (!UseHighSparseLayout) {
      (void)level;
      return false;
    }
    return level > 0 && level >= high_level_begin_ && level < level_count();
  }

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

  bool level_embeds_min_summary(std::size_t level) const {
    return level > 0 && !is_high_level(level) &&
           fanout_at_level(level) <= kEmbeddedSummaryMaxEntries;
  }

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

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

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

  const HighChildMetadata* high_child_metadata_begin(std::size_t level,
                                                     std::size_t node) const {
    return high_child_metadata_.data() +
           high_flat_index(level, node) * kHighLevelFanout;
  }

  const HighChildMetadata& high_child_metadata_at(std::size_t high_flat,
                                                  std::size_t slot) const {
    return high_child_metadata_[high_flat * kHighLevelFanout + slot];
  }

  HighChildMetadata& mutable_high_child_metadata_at(std::size_t high_flat,
                                                    std::size_t slot) {
    return high_child_metadata_[high_flat * kHighLevelFanout + slot];
  }

  std::uint8_t* mutable_high_sparse_min_slots_begin(std::size_t high_flat) {
    return high_sparse_min_slots_.data() + high_flat * kHighSparseSlotsPerNode;
  }

  const std::uint8_t* high_sparse_min_slots_begin(std::size_t high_flat) const {
    return high_sparse_min_slots_.data() + high_flat * kHighSparseSlotsPerNode;
  }

  std::size_t better_high_child_slot(std::size_t level,
                                     std::size_t node,
                                     std::size_t left_slot,
                                     std::size_t right_slot) const {
    const DepthCandidate left =
        high_child_min_candidate(level, node, left_slot);
    const DepthCandidate right =
        high_child_min_candidate(level, node, right_slot);
    return better_candidate(left, right).position == right.position ? right_slot
                                                                    : left_slot;
  }

  void build_high_sparse_min_slots(std::size_t level,
                                   std::size_t node,
                                   std::size_t count) {
    const std::size_t high_flat = high_flat_index(level, node);
    std::uint8_t* table = mutable_high_sparse_min_slots_begin(high_flat);
    for (std::size_t slot = 0; slot < count; ++slot) {
      table[slot] = static_cast<std::uint8_t>(slot);
    }
 
    for (std::size_t table_level = 1; table_level < kHighSparseTableLevels;
         ++table_level) {
      const std::size_t span = std::size_t{1} << table_level;
      if (span > count) {
        break;
      }
      const std::size_t half_span = span >> 1;
      const std::uint8_t* previous =
          table + (table_level - 1) * kHighLevelFanout;
      std::uint8_t* current = table + table_level * kHighLevelFanout;
      for (std::size_t slot = 0; slot + span <= count; ++slot) {
        current[slot] = static_cast<std::uint8_t>(better_high_child_slot(
            level, node, previous[slot], previous[slot + half_span]));
      }
    }
  }

  std::size_t high_sparse_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 {
    if (slot_left >= slot_right || slot_right > count) {
      return npos;
    }
    const std::size_t length = slot_right - slot_left;
    if (length == 1) {
      return slot_left;
    }
 
    const std::size_t high_flat = high_flat_index(level, node);
    const std::size_t table_level = std::bit_width(length) - 1;
    const std::size_t span = std::size_t{1} << table_level;
    const std::uint8_t* table =
        high_sparse_min_slots_begin(high_flat) + table_level * kHighLevelFanout;
    return better_high_child_slot(level, node, table[slot_left],
                                  table[slot_right - span]);
  }
 
  std::span<const std::uint64_t> input_bits_;
  std::size_t depth_count_ = 0;
  std::optional<RankSelectSupport<>> owned_rank_index_;
  const RankSelectSupport<>* external_rank_index_ = nullptr;
  std::vector<Bp512Selector> internal_selectors_;
  std::vector<Index> internal_min_positions_;
  std::vector<std::int64_t> internal_min_depths_;
  std::vector<HighChildMetadata> high_child_metadata_;
  std::vector<std::uint8_t> high_sparse_min_slots_;
  std::vector<std::size_t> internal_level_offsets_;
  std::vector<std::size_t> min_summary_level_offsets_;
  std::vector<std::size_t> high_level_offsets_;
  std::vector<std::size_t> level_sizes_;
  std::vector<std::size_t> level_position_spans_;
  std::vector<std::size_t> level_fanouts_;
  std::size_t high_level_begin_ = std::numeric_limits<std::size_t>::max();
};
 
}  // namespace detail

template <class T,
          class Compare = std::less<T>,
          class Index = std::size_t,
          std::size_t LeafSize = 512,
          bool UseTopSparseOverlay = true>
class CartesianHybridBTree
    : public RmqBase<CartesianHybridBTree<T,
                                          Compare,
                                          Index,
                                          LeafSize,
                                          UseTopSparseOverlay>,
                     T>,
      public SerializationBase<CartesianHybridBTree<T,
                                                    Compare,
                                                    Index,
                                                    LeafSize,
                                                    UseTopSparseOverlay>> {
 private:
  using BpDepthRmq = detail::HybridBTreePlusMinusOne<Index, LeafSize, true, 1>;
 
  struct TopCandidate {
    Index position = std::numeric_limits<Index>::max();
  };
  static_assert(sizeof(TopCandidate) == sizeof(Index));
 
  struct LoadedState {
    pixie::AlignedStorage bp_bits_;
    std::size_t bp_bit_count_ = 0;
    std::vector<TopCandidate> top_sparse_candidates_;
    std::size_t top_block_size_ = 4096;
    std::size_t top_block_count_ = 0;
    std::size_t top_sparse_levels_ = 0;
    std::optional<RankSelectSupport<>> bp_index_;
    BpDepthRmq bp_depth_rmq_;
  };
 
  struct LoadTag {};
 
  static constexpr std::array<std::uint8_t, 8> kSerializationMagic = {
      'P', 'I', 'X', 'I', 'E', 'R', 'M', 'Q'};
  static constexpr std::uint32_t kSerializationVersion = 4;
  static constexpr std::size_t kSerializationHeaderBytes = 48;
 
 public:
  static_assert(std::is_unsigned_v<Index>,
                "CartesianHybridBTree index type must be unsigned");
  static_assert(LeafSize != 0 && LeafSize % 512 == 0,
                "CartesianHybridBTree leaf size must be a positive "
                "multiple of 512");
 
  using Self =
      CartesianHybridBTree<T, Compare, Index, LeafSize, UseTopSparseOverlay>;
 
  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 kMinTopSparseBlockSize = 4096;
  static constexpr std::size_t kMaxTopSparseBlocks = std::size_t{1} << 14;
  static constexpr bool kUseTopSparseOverlay = UseTopSparseOverlay;
  static constexpr bool kSerializationSupported =
      std::same_as<T, std::int64_t> &&
      std::same_as<Compare, std::less<std::int64_t>> &&
      std::same_as<Index, std::size_t> && LeafSize == 512 &&
      UseTopSparseOverlay;

  CartesianHybridBTree() = default;

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

  void serialize_impl(BinaryWriter& writer) const
    requires(kSerializationSupported)
  {
    validate_serialized_state(DeserializationValidation::kQuick);
 
    const std::size_t artifact_begin = writer.size_bytes();
    pixie::detail::write_magic(writer, kSerializationMagic);
    writer.write_u32(kSerializationVersion);
    writer.write_u8(pixie::detail::kLittleEndianMarker);
    writer.write_u8(sizeof(std::uint64_t));
    writer.write_u8(sizeof(T));
    writer.write_u8(sizeof(Index));
    const std::size_t artifact_size_position = writer.write_u64_placeholder();
    writer.write_size(values_.size());
    writer.write_size(LeafSize);
    writer.write_u32(static_cast<std::uint32_t>(UseTopSparseOverlay));
    writer.write_u32(0);
 
    writer.write_size(bp_bit_count_);
    writer.write_size(top_block_size_);
    writer.write_size(top_block_count_);
    writer.write_size(top_sparse_levels_);
    bp_bits_.serialize(writer);
    writer.write_size(top_sparse_candidates_.size());
    for (const TopCandidate candidate : top_sparse_candidates_) {
      pixie::detail::write_integral(writer, candidate.position);
    }
    writer.write_u8(static_cast<std::uint8_t>(bp_index_.has_value()));
    if (bp_index_) {
      bp_index_->serialize(writer);
    }
    bp_depth_rmq_.serialize_metadata(writer);
 
    const std::size_t unpadded_size = writer.size_bytes() - artifact_begin;
    writer.write_zeros(
        (sizeof(std::uint64_t) - unpadded_size % sizeof(std::uint64_t)) %
        sizeof(std::uint64_t));
    const std::size_t artifact_size = writer.size_bytes() - artifact_begin;
    writer.patch_u64(artifact_size_position,
                     static_cast<std::uint64_t>(artifact_size));
  }

  static Self deserialize_impl(
      BinaryReader& reader,
      std::span<const std::int64_t> values,
      DeserializationValidation validation = DeserializationValidation::kQuick)
    requires(kSerializationSupported)
  {
    BinaryReader candidate = reader;
    const std::size_t available_size = candidate.remaining();
    pixie::detail::require_magic(candidate, kSerializationMagic);
    if (candidate.read_u32() != kSerializationVersion ||
        candidate.read_u8() != pixie::detail::kLittleEndianMarker ||
        candidate.read_u8() != sizeof(std::uint64_t) ||
        candidate.read_u8() != sizeof(T) ||
        candidate.read_u8() != sizeof(Index)) {
      throw std::invalid_argument("Incompatible serialized RMQ artifact");
    }
    const std::size_t artifact_size = pixie::detail::checked_artifact_size(
        candidate.read_u64(), kSerializationHeaderBytes, available_size);
    const std::size_t source_value_count = candidate.read_size();
    if (candidate.read_size() != LeafSize ||
        candidate.read_u32() !=
            static_cast<std::uint32_t>(UseTopSparseOverlay) ||
        candidate.read_u32() != 0) {
      throw std::invalid_argument("Incompatible serialized RMQ configuration");
    }
 
    BinaryReader payload =
        candidate.read_subreader(artifact_size - kSerializationHeaderBytes);
    LoadedState state;
    state.bp_bit_count_ = payload.read_size();
    state.top_block_size_ = payload.read_size();
    state.top_block_count_ = payload.read_size();
    state.top_sparse_levels_ = payload.read_size();
    state.bp_bits_ = deserialize_aligned_storage(payload);
 
    const std::size_t candidate_count = payload.read_size();
    const std::vector<TopCandidate> empty_candidates;
    if (candidate_count > empty_candidates.max_size()) {
      throw std::length_error(
          "Serialized RMQ top candidate count is too large");
    }
    if (candidate_count > payload.remaining() / sizeof(Index)) {
      throw std::invalid_argument("Truncated serialized RMQ top candidates");
    }
    state.top_sparse_candidates_.resize(candidate_count);
    for (TopCandidate& candidate : state.top_sparse_candidates_) {
      candidate.position = pixie::detail::read_integral<Index>(payload);
    }
 
    const std::uint8_t has_rank_index = payload.read_u8();
    if (has_rank_index > 1) {
      throw std::invalid_argument("Invalid serialized RMQ rank-index marker");
    }
    const std::size_t bp_word_count =
        state.bp_bit_count_ == 0 ? 0 : 1 + (state.bp_bit_count_ - 1) / 64;
    if (bp_word_count > state.bp_bits_.as_words64().size()) {
      throw std::invalid_argument("Serialized RMQ BP storage is too small");
    }
    if (has_rank_index != 0) {
      state.bp_index_ = RankSelectSupport<>::deserialize(
          payload, state.bp_bits_.as_words64().first(bp_word_count),
          validation);
    }
 
    const RankSelectSupport<> empty_rank_index;
    const RankSelectSupport<>& rank_index =
        state.bp_index_ ? *state.bp_index_ : empty_rank_index;
    state.bp_depth_rmq_ = BpDepthRmq::deserialize_metadata(
        state.bp_bits_.as_words64(), rank_index, payload);
    payload.require_zero_padding(sizeof(std::uint64_t) - 1);
 
    if (source_value_count != values.size()) {
      throw std::invalid_argument(
          "Serialized RMQ source value count is inconsistent");
    }
    validate_loaded_state(values, state, validation);
    reader = candidate;
    return Self(LoadTag{}, values, std::move(state));
  }

  CartesianHybridBTree(const CartesianHybridBTree& other)
      : values_(other.values_),
        compare_(other.compare_),
        bp_bits_(other.bp_bits_),
        bp_bit_count_(other.bp_bit_count_),
        top_sparse_candidates_(other.top_sparse_candidates_),
        top_block_size_(other.top_block_size_),
        top_block_count_(other.top_block_count_),
        top_sparse_levels_(other.top_sparse_levels_) {
    reset_bp_indexes();
  }

  CartesianHybridBTree& operator=(const CartesianHybridBTree& other) {
    if (this == &other) {
      return *this;
    }
    values_ = other.values_;
    compare_ = other.compare_;
    bp_bits_ = other.bp_bits_;
    bp_bit_count_ = other.bp_bit_count_;
    top_sparse_candidates_ = other.top_sparse_candidates_;
    top_block_size_ = other.top_block_size_;
    top_block_count_ = other.top_block_count_;
    top_sparse_levels_ = other.top_sparse_levels_;
    reset_bp_indexes();
    return *this;
  }

  CartesianHybridBTree(CartesianHybridBTree&& other) noexcept
      : values_(other.values_),
        compare_(std::move(other.compare_)),
        bp_bits_(std::move(other.bp_bits_)),
        bp_bit_count_(other.bp_bit_count_),
        top_sparse_candidates_(std::move(other.top_sparse_candidates_)),
        top_block_size_(other.top_block_size_),
        top_block_count_(other.top_block_count_),
        top_sparse_levels_(other.top_sparse_levels_) {
    other.values_ = std::span<const T>();
    other.bp_bit_count_ = 0;
    other.top_block_size_ = kMinTopSparseBlockSize;
    other.top_block_count_ = 0;
    other.top_sparse_levels_ = 0;
    reset_bp_indexes();
  }

  CartesianHybridBTree& operator=(CartesianHybridBTree&& other) noexcept {
    if (this == &other) {
      return *this;
    }
    values_ = other.values_;
    compare_ = std::move(other.compare_);
    bp_bits_ = std::move(other.bp_bits_);
    bp_bit_count_ = other.bp_bit_count_;
    top_sparse_candidates_ = std::move(other.top_sparse_candidates_);
    top_block_size_ = other.top_block_size_;
    top_block_count_ = other.top_block_count_;
    top_sparse_levels_ = other.top_sparse_levels_;
    other.values_ = std::span<const T>();
    other.bp_bit_count_ = 0;
    other.top_block_size_ = kMinTopSparseBlockSize;
    other.top_block_count_ = 0;
    other.top_sparse_levels_ = 0;
    reset_bp_indexes();
    return *this;
  }

  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()) {
      return npos;
    }
    if constexpr (!UseTopSparseOverlay) {
      return cartesian_arg_min(left, right);
    }
    if (right - left <= top_block_size_) {
      return cartesian_arg_min(left, right);
    }
    const std::size_t top_answer = top_sparse_arg_min(left, right);
    if (top_answer != npos) {
      return top_answer;
    }
    return cartesian_arg_min(left, right);
  }

  std::size_t bp_bit_count() const { return bp_bit_count_; }

  std::span<const std::uint64_t> bp_words() const {
    return bp_storage_words().first(bp_word_count());
  }

  static std::size_t top_sparse_block_size_for(std::size_t value_count) {
    if (value_count == 0) {
      return kMinTopSparseBlockSize;
    }
    return std::max(kMinTopSparseBlockSize,
                    ceil_div(value_count, kMaxTopSparseBlocks));
  }

  static std::size_t top_sparse_block_count_for(std::size_t value_count) {
    if (value_count == 0) {
      return 0;
    }
    return ceil_div(value_count, top_sparse_block_size_for(value_count));
  }

  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 {
    return sizeof(*this) + bp_bits_.allocated_bytes() +
           pixie::vector_capacity_bytes(top_sparse_candidates_) +
           pixie::optional_nested_owned_memory_bytes(bp_index_) +
           pixie::nested_owned_memory_bytes(bp_depth_rmq_);
  }
 
 private:
  static pixie::AlignedStorage deserialize_aligned_storage(
      BinaryReader& reader) {
    const std::size_t size = reader.read_size();
    const std::span<const std::byte> bytes = reader.read_bytes(size);
    if (size > std::numeric_limits<std::size_t>::max() / 8) {
      throw std::length_error("Serialized RMQ storage is too large");
    }
    pixie::AlignedStorage result(size * 8);
    std::ranges::copy(bytes, result.writable_bytes().begin());
    return result;
  }
 
  template <class State>
  static void validate_loaded_state(std::span<const T> values,
                                    const State& state,
                                    DeserializationValidation validation) {
    if (values.size() > (static_cast<std::size_t>(invalid_index) - 1) / 2) {
      throw std::length_error("Serialized RMQ value count is too large");
    }
    const std::size_t expected_bp_bit_count = 2 * values.size();
    if (state.bp_bit_count_ != expected_bp_bit_count) {
      throw std::invalid_argument("Invalid serialized RMQ BP bit count");
    }
 
    std::size_t expected_padded_bits = 0;
    if (expected_bp_bit_count != 0) {
      const std::size_t depth_count = expected_bp_bit_count + 1;
      if (depth_count > std::numeric_limits<std::int64_t>::max()) {
        throw std::length_error("Serialized RMQ BP depth count is too large");
      }
      const std::size_t leaf_count = 1 + (depth_count - 1) / LeafSize;
      if (leaf_count > std::numeric_limits<std::size_t>::max() / LeafSize) {
        throw std::length_error("Serialized RMQ padded BP size is too large");
      }
      expected_padded_bits = leaf_count * LeafSize;
    }
    if (state.bp_bits_.size_bits() != expected_padded_bits) {
      throw std::invalid_argument("Invalid serialized RMQ BP storage size");
    }
 
    if (values.empty()) {
      if (state.top_block_size_ != kMinTopSparseBlockSize ||
          state.top_block_count_ != 0 || state.top_sparse_levels_ != 0 ||
          !state.top_sparse_candidates_.empty() ||
          state.bp_index_.has_value() || !state.bp_depth_rmq_.empty()) {
        throw std::invalid_argument("Invalid serialized empty RMQ metadata");
      }
      const RankSelectSupport<> empty_rank;
      state.bp_depth_rmq_.validate_serialized_state(empty_rank, validation);
      return;
    }
 
    const std::size_t expected_block_size =
        top_sparse_block_size_for(values.size());
    const std::size_t expected_block_count =
        top_sparse_block_count_for(values.size());
    const std::size_t expected_levels = std::bit_width(expected_block_count);
    if (state.top_block_size_ != expected_block_size ||
        state.top_block_count_ != expected_block_count ||
        state.top_sparse_levels_ != expected_levels ||
        expected_block_count >
            std::numeric_limits<std::size_t>::max() / expected_levels ||
        state.top_sparse_candidates_.size() !=
            expected_block_count * expected_levels ||
        !state.bp_index_.has_value() ||
        state.bp_index_->size() != expected_bp_bit_count ||
        state.bp_index_->supports_select1() ||
        !state.bp_index_->supports_select0() ||
        state.bp_depth_rmq_.size() != expected_bp_bit_count + 1) {
      throw std::invalid_argument("Invalid serialized RMQ index metadata");
    }
 
    for (std::size_t level = 0;
         validation == DeserializationValidation::kQuick &&
         level < expected_levels;
         ++level) {
      const std::size_t span = std::size_t{1} << level;
      for (std::size_t block = 0; block < expected_block_count; ++block) {
        const std::size_t position = static_cast<std::size_t>(
            state.top_sparse_candidates_[level * expected_block_count + block]
                .position);
        const bool populated = block + span <= expected_block_count;
        if (!populated) {
          if (position != static_cast<std::size_t>(invalid_index)) {
            throw std::invalid_argument(
                "Invalid serialized RMQ sparse-table padding");
          }
          continue;
        }
        const std::size_t begin = block * expected_block_size;
        const std::size_t end =
            std::min(values.size(), (block + span) * expected_block_size);
        if (position < begin || position >= end) {
          throw std::invalid_argument(
              "Invalid serialized RMQ sparse-table candidate");
        }
      }
    }
    state.bp_depth_rmq_.validate_serialized_state(*state.bp_index_, validation);
    if (validation == DeserializationValidation::kFull) {
      validate_exact_source(values, state, expected_padded_bits);
    }
  }
 
  template <class State>
  static void validate_exact_source(std::span<const T> values,
                                    const State& state,
                                    std::size_t expected_padded_bits) {
    const std::span<const std::uint64_t> bp_words = state.bp_bits_.as_words64();
    const auto bp_bit = [bp_words](std::size_t position) {
      return ((bp_words[position >> 6] >> (position & 63)) & 1u) != 0;
    };
    const auto require_bp_bit = [&](std::size_t position, bool expected) {
      if (bp_bit(position) != expected) {
        throw std::invalid_argument(
            "Serialized RMQ Cartesian BP disagrees with source values");
      }
    };
    const auto better_position = [&](std::size_t left, std::size_t right) {
      if (right == static_cast<std::size_t>(invalid_index)) {
        return left;
      }
      if (left == static_cast<std::size_t>(invalid_index)) {
        return right;
      }
      if (std::less<T>{}(values[right], values[left])) {
        return right;
      }
      if (std::less<T>{}(values[left], values[right])) {
        return left;
      }
      return std::min(left, right);
    };
 
    utils::SuccinctIncreasingStack stack(values.size());
    std::size_t write_position = 2 * values.size();
    std::size_t block_minimum = static_cast<std::size_t>(invalid_index);
    for (std::size_t position = values.size(); position-- > 0;) {
      while (!stack.empty()) {
        const std::size_t top_position = values.size() - 1 - stack.top();
        if (std::less<T>{}(values[top_position], values[position])) {
          break;
        }
        stack.pop();
        require_bp_bit(--write_position, true);
      }
      stack.push(values.size() - 1 - position);
      require_bp_bit(--write_position, false);
 
      block_minimum = better_position(block_minimum, position);
      if (position % state.top_block_size_ == 0) {
        const std::size_t block = position / state.top_block_size_;
        if (static_cast<std::size_t>(
                state.top_sparse_candidates_[block].position) !=
            block_minimum) {
          throw std::invalid_argument(
              "Serialized RMQ block minimum disagrees with source values");
        }
        block_minimum = static_cast<std::size_t>(invalid_index);
      }
    }
    while (write_position != 0) {
      require_bp_bit(--write_position, true);
    }
 
    for (std::size_t position = 2 * values.size();
         position < expected_padded_bits; ++position) {
      if (bp_bit(position)) {
        throw std::invalid_argument("Serialized RMQ BP padding is non-zero");
      }
    }
 
    for (std::size_t level = 1; level < state.top_sparse_levels_; ++level) {
      const std::size_t span = std::size_t{1} << level;
      const std::size_t half_span = span >> 1;
      const std::size_t current_offset = level * state.top_block_count_;
      const std::size_t previous_offset = (level - 1) * state.top_block_count_;
      for (std::size_t block = 0; block < state.top_block_count_; ++block) {
        const std::size_t actual = static_cast<std::size_t>(
            state.top_sparse_candidates_[current_offset + block].position);
        if (block + span > state.top_block_count_) {
          if (actual != static_cast<std::size_t>(invalid_index)) {
            throw std::invalid_argument(
                "Invalid serialized RMQ sparse-table padding");
          }
          continue;
        }
        const std::size_t left = static_cast<std::size_t>(
            state.top_sparse_candidates_[previous_offset + block].position);
        const std::size_t right = static_cast<std::size_t>(
            state.top_sparse_candidates_[previous_offset + block + half_span]
                .position);
        if (actual != better_position(left, right)) {
          throw std::invalid_argument(
              "Serialized RMQ sparse-table candidate disagrees with source");
        }
      }
    }
  }
 
  void validate_serialized_state(DeserializationValidation validation) const {
    validate_loaded_state(values_, *this, validation);
  }
 
  CartesianHybridBTree(LoadTag, std::span<const T> values, LoadedState&& state)
      : values_(values),
        compare_(),
        bp_bits_(std::move(state.bp_bits_)),
        bp_bit_count_(state.bp_bit_count_),
        top_sparse_candidates_(std::move(state.top_sparse_candidates_)),
        top_block_size_(state.top_block_size_),
        top_block_count_(state.top_block_count_),
        top_sparse_levels_(state.top_sparse_levels_),
        bp_index_(std::move(state.bp_index_)),
        bp_depth_rmq_(std::move(state.bp_depth_rmq_)) {
    if (bp_index_) {
      bp_depth_rmq_.restore_external_sources(bp_bits_.as_words64(), *bp_index_);
    } else {
      bp_depth_rmq_.restore_empty_sources();
    }
  }

  std::size_t cartesian_arg_min(std::size_t left, std::size_t right) const {
    const std::size_t first_close = select_close_position(left + 1);
    const std::size_t last_close = select_close_position(right);
    if (first_close == npos || last_close == npos || first_close > last_close) {
      return npos;
    }
 
    const std::size_t shifted_min =
        bp_depth_rmq_.arg_min(first_close + 1, last_close + 2);
    if (shifted_min == npos || shifted_min == 0) {
      return npos;
    }
    const RankSelectSupport<>& bp_index = *bp_index_;
    const std::size_t answer = bp_index.rank0(shifted_min) - 1;
    return answer < values_.size() ? answer : npos;
  }

  void build() {
    bp_bits_.resize(0);
    bp_bit_count_ = 0;
    top_sparse_candidates_.clear();
    top_block_size_ = kMinTopSparseBlockSize;
    top_block_count_ = 0;
    top_sparse_levels_ = 0;
    reset_bp_indexes();
 
    if (values_.empty()) {
      return;
    }
    if (values_.size() > (static_cast<std::size_t>(invalid_index) - 1) / 2) {
      throw std::length_error(
          "CartesianHybridBTree RMQ index type is too small");
    }
 
    bp_bit_count_ = 2 * values_.size();
    bp_bits_.resize(padded_bp_bit_capacity());
    std::ranges::fill(bp_storage_words(), std::uint64_t{0});
    build_bp_bits();
    if constexpr (UseTopSparseOverlay) {
      build_top_sparse_table();
    }
    reset_bp_indexes();
  }

  void build_bp_bits() {
    utils::SuccinctIncreasingStack stack(values_.size());
    std::size_t write_position = bp_bit_count_;
    std::span<std::uint64_t> words = bp_storage_words();
    const auto prepend_open = [&]() {
      --write_position;
      words[write_position >> 6] |= std::uint64_t{1} << (write_position & 63);
    };
 
    for (std::size_t i = values_.size(); i-- > 0;) {
      while (!stack.empty() &&
             !compare_(values_[stack_index(stack.top())], values_[i])) {
        stack.pop();
        prepend_open();
      }
      stack.push(stack_key(i));
      --write_position;
    }
 
    while (write_position != 0) {
      prepend_open();
    }
  }

  std::size_t stack_key(std::size_t value_index) const {
    return values_.size() - 1 - value_index;
  }

  std::size_t stack_index(std::size_t key) const {
    return values_.size() - 1 - key;
  }

  void reset_bp_indexes() {
    bp_index_.reset();
    bp_depth_rmq_ = BpDepthRmq();
    if (bp_bit_count_ == 0) {
      return;
    }
    const std::span<const std::uint64_t> words = bp_words();
    const std::span<const std::uint64_t> padded_words = bp_storage_words();
    // TODO: try incorporating rank/select information into the tree.
    bp_index_.emplace(words, bp_bit_count_,
                      RankSelectSupport<>::SelectSupport::kSelect0,
                      values_.size());
    bp_depth_rmq_ = BpDepthRmq(padded_words, bp_bit_count_ + 1, *bp_index_);
  }

  std::size_t bp_word_count() const { return ceil_div(bp_bit_count_, 64); }

  std::size_t padded_bp_bit_capacity() const {
    if (bp_bit_count_ == 0) {
      return 0;
    }
    const std::size_t depth_count = bp_bit_count_ + 1;
    return ceil_div(depth_count, LeafSize) * LeafSize;
  }

  std::span<std::uint64_t> bp_storage_words() {
    return bp_bits_.writable_words64();
  }

  std::span<const std::uint64_t> bp_storage_words() const {
    return bp_bits_.as_words64();
  }

  void build_top_sparse_table() {
    top_sparse_candidates_.clear();
    top_block_size_ = top_sparse_block_size_for(values_.size());
    top_block_count_ = ceil_div(values_.size(), top_block_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]);
      }
    }
  }

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

  TopCandidate make_top_candidate(std::size_t position) const {
    if (position >= values_.size()) {
      return {};
    }
    return {static_cast<Index>(position)};
  }

  bool valid_value_position(std::size_t position) const {
    return position != npos &&
           position != static_cast<std::size_t>(invalid_index) &&
           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(cartesian_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(cartesian_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::size_t select_close_position(std::size_t rank) const {
    if (rank == 0 || rank > values_.size() || !bp_index_) {
      return npos;
    }
    const std::size_t position = bp_index_->select0(rank);
    return position < bp_bit_count_ ? position : npos;
  }
 
  std::span<const T> values_;
  Compare compare_;
  pixie::AlignedStorage bp_bits_;
  std::size_t bp_bit_count_ = 0;
  std::vector<TopCandidate> top_sparse_candidates_;
  std::size_t top_block_size_ = kMinTopSparseBlockSize;
  std::size_t top_block_count_ = 0;
  std::size_t top_sparse_levels_ = 0;
  std::optional<RankSelectSupport<>> bp_index_;
  BpDepthRmq bp_depth_rmq_;
};

template <class T,
          class Compare = std::less<T>,
          class Index = std::size_t,
          std::size_t LeafSize = 512>
using CartesianBTree = CartesianHybridBTree<T, Compare, Index, LeafSize, false>;
 
}  // namespace pixie::rmq