diff --git a/README.md b/README.md
index 564058f..516fa64 100644
--- a/README.md
+++ b/README.md
@@ -58,7 +58,7 @@ After building:
 
 ### BitVector
 
-Benchmarks are random 50/50 0-1 bitvectors up to $2^34$ bits.
+Benchmarks are random 50/50 0-1 bitvectors up to $2^{34}$ bits.
 
 ```bash
 ./benchmarks
@@ -70,7 +70,7 @@ Benchmarks are random 50/50 0-1 bitvectors up to $2^34$ bits.
 ./bench_rmm
 ```
 
-For visualization, write the CSV output to a file using `--benchmark_out=<file>` (e.g. `./bench_rmm --benchmark_out=rmm_bench.csv`) and plot it with `misc/plot_rmm.py`.
+For visualization, write the JSON output to a file using `--benchmark_out=<file>` (e.g. `./bench_rmm --benchmark_out=rmm_bench.json`) and plot it with `misc/plot_rmm.py`.
 
 ---
 
diff --git a/include/bits.h b/include/bits.h
index 4f0bb25..15f0e6c 100644
--- a/include/bits.h
+++ b/include/bits.h
@@ -405,3 +405,108 @@ void rank_32x8(const uint8_t* x, uint8_t* result) {
   }
 #endif
 }
+
+/**
+ * @brief Efficiently searches for the first occurrence of a 16-bit value in
+ * the range [@p begin, @p end_excl) using AVX2 when available.
+ * @details Loads 16 consecutive int16_t elements (256 bits) per iteration.
+ * Compares them against the @p target value using vectorized equality.
+ * If any match is found, extracts the index of the first matching lane from
+ * the comparison mask. Falls back to a scalar tail loop for leftover
+ * elements, or to a fully scalar search if AVX2 is not supported.
+ * @returns The index of the first match, or @p npos if the value is not found.
+ */
+static inline size_t find_forward_equal_i16_avx2(const int16_t* arr,
+                                                 const size_t& begin,
+                                                 const size_t& end_excl,
+                                                 const int16_t& target,
+                                                 const size_t& npos) noexcept {
+#ifdef PIXIE_AVX2_SUPPORT
+  static constexpr size_t STEP = 16;
+  __m256i vtarget = _mm256_set1_epi16(target);
+  size_t i = begin;
+  size_t n = end_excl;
+  for (; i + STEP <= n; i += STEP) {
+    unsigned mask = _mm256_movemask_epi8(_mm256_cmpeq_epi16(
+        _mm256_loadu_si256(reinterpret_cast<const __m256i*>(arr + i)),
+        vtarget));
+    if (mask) {
+      return i + (std::countr_zero(mask) >> 1);
+    }
+  }
+  for (; i < n; ++i) {
+    if (arr[i] == target) {
+      return i;
+    }
+  }
+#else
+  for (size_t i = begin; i < end_excl; ++i) {
+    if (arr[i] == target) {
+      return i;
+    }
+  }
+#endif
+  return npos;
+}
+
+/**
+ * @brief Performs a backward search for a 16-bit value in a given range.
+ * @details Scans the array segment [@p begin .. @p end_incl] from right to
+ * left.
+ * If AVX2 is available, processes data in 256-bit blocks (16 × int16_t) using
+ * vectorized equality comparison for higher throughput. Falls back to a
+ * scalar backward scan when AVX2 is not supported. Returns the index of the
+ * rightmost occurrence of @p target, or @p npos if no match is found.
+ */
+static inline size_t find_backward_equal_i16_avx2(const int16_t* arr,
+                                                  const size_t& begin,
+                                                  const size_t& end_incl,
+                                                  const int16_t& target,
+                                                  const size_t& npos) noexcept {
+  if (begin > end_incl) {
+    return npos;
+  }
+#ifdef PIXIE_AVX2_SUPPORT
+  static constexpr size_t STEP = 16;
+  size_t len = end_incl + 1 - begin;
+  size_t nblocks = len / STEP;
+  __m256i vtarget = _mm256_set1_epi16(target);
+  if (nblocks > 0) {
+    size_t first_block = begin + (len % STEP);
+    for (size_t p = first_block + (nblocks - 1) * STEP;;) {
+      unsigned mask = _mm256_movemask_epi8(_mm256_cmpeq_epi16(
+          _mm256_loadu_si256(reinterpret_cast<const __m256i*>(arr + p)),
+          vtarget));
+      if (mask) {
+        return p + ((31u - std::countl_zero(mask)) >> 1);
+      }
+      if (p == first_block) {
+        break;
+      }
+      p -= STEP;
+    }
+
+    for (size_t i = first_block; i > begin;) {
+      --i;
+      if (arr[i] == target) {
+        return i;
+      }
+    }
+  } else {
+    for (size_t i = end_incl + 1; i > begin;) {
+      --i;
+      if (arr[i] == target) {
+        return i;
+      }
+    }
+  }
+#else
+  for (size_t i = end_incl + 1; i > begin;) {
+    --i;
+    if (arr[i] == target) {
+      return i;
+    }
+  }
+#endif
+  return npos;
+}
diff --git a/include/rmm_tree.h b/include/rmm_tree.h
index 263b668..30b2ad9 100644
--- a/include/rmm_tree.h
+++ b/include/rmm_tree.h
@@ -1,4 +1,6 @@
 #pragma once
+#include <immintrin.h>
+
 #include <algorithm>
 #include <array>
 #include <bit>
@@ -9,6 +11,8 @@
 #include <string>
 #include <vector>
 
+#include "bits.h"
+
 namespace pixie {
 /**
  * @brief Range min–max tree over a bitvector (LSB-first) tailored for
@@ -71,6 +75,13 @@ class RmMTree {
   // segment (to handle "10" crossing) both needed for: rank10, select10.
   std::vector<uint8_t> node_first_bit, node_last_bit;
 
+  // leaf_prefix[k * leaf_prefix_stride + j] =
+  // excess(block_begin + j) - excess(block_begin), j in [0..leaf_size],
+  // block_begin = k * block_bits. needed for:
+  // fwdsearch/bwdsearch/close/open/enclose
+  std::vector<int16_t> leaf_prefix;
+  size_t leaf_prefix_stride = 0;
+
  public:
   /**
    * @brief Sentinel for "not found".
@@ -90,23 +101,23 @@ class RmMTree {
 
   /**
    * @brief Build from a '0'/'1' string.
-   * @param bp Bitstring (characters '0' and '1').
+   * @param bit_string Bitstring (characters '0' and '1').
    * @param leaf_block_bits Desired leaf size (power of two, 0 = auto).
    * @param max_overhead Max allowed overhead fraction (<0 to disable
    * constraint).
    * @details Block size priority: (1) respect @p max_overhead, (2) explicit @p
    * leaf_block_bits, (3) set to ceil_pow2(log2(num_bits)).
    */
-  explicit RmMTree(const std::string& bp,
+  explicit RmMTree(const std::string& bit_string,
                    const size_t& leaf_block_bits /*0=auto*/ = 0,
                    const float& max_overhead /*<0=off*/ = -1.0) {
-    build_from_string(bp, leaf_block_bits, max_overhead);
+    build_from_string(bit_string, leaf_block_bits, max_overhead);
   }
 
   /**
    * @brief Build from 64-bit words (LSB-first).
    * @param words Array of words holding bits LSB-first.
-   * @param Nbits Number of valid bits.
+   * @param bit_count Number of valid bits.
    * @param leaf_block_bits Desired leaf size (power of two, 0 = auto).
    * @param max_overhead Max allowed overhead fraction (<0 to disable
    * constraint).
@@ -114,720 +125,775 @@ class RmMTree {
    * leaf_block_bits, (3) set to ceil_pow2(log2(num_bits)).
    */
   explicit RmMTree(const std::vector<std::uint64_t>& words,
-                   size_t Nbits,
+                   size_t bit_count,
                    const size_t& leaf_block_bits /*0=auto*/ = 0,
                    const float& max_overhead /*<0=off*/ = -1.0) {
-    build_from_words(words, Nbits, leaf_block_bits, max_overhead);
+    build_from_words(words, bit_count, leaf_block_bits, max_overhead);
   }
 
   // --------- queries: rank/select/excess ----------
 
   /**
-   * @brief Number of ones in prefix [0, i).
-   * @details Returns 0 for i==0. Complexity: O(log n) with small constants.
+   * @brief Number of ones in prefix [0, @p end_position).
+   * @details Returns 0 for @p end_position == 0.
    */
-  size_t rank1(const size_t& i) const {
-    if (i == 0) {
+  size_t rank1(const size_t& end_position) const {
+    if (end_position == 0) {
       return 0;
     }
-    const size_t blk = block_of(i - 1);
-    size_t ans = 0;
-    if (blk > 0) {
-      const auto nodes = cover_blocks(0, blk - 1);
-      for (const size_t& v : nodes) {
-        ans += ones_in_node(v);
-      }
-    }
-    const size_t Lb = blk * block_bits;
-    const size_t Rb = std::min(num_bits, Lb + block_bits);
-    ans += rank1_in_block(Lb, std::min(i, Rb));
-    return ans;
+    const size_t block_index = block_of(end_position - 1);
+    size_t ones_count = 0;
+    if (block_index > 0) {
+      size_t nodes_buffer[64];
+      const size_t node_count =
+          cover_blocks_collect(0, block_index - 1, nodes_buffer);
+      for (size_t j = 0; j < node_count; ++j) {
+        ones_count += ones_in_node(nodes_buffer[j]);
+      }
+    }
+    const size_t block_begin = block_index * block_bits;
+    const size_t block_end = std::min(num_bits, block_begin + block_bits);
+    ones_count +=
+        rank1_in_block(block_begin, std::min(end_position, block_end));
+    return ones_count;
   }
 
   /**
-   * @brief Number of zeros in prefix [0, i).
-   * @details Computed as i - rank1(i).
+   * @brief Number of zeros in prefix [0, @p end_position).
+   * @details Computed as @p end_position - rank1(@p end_position).
    */
-  size_t rank0(const size_t& i) const { return i - rank1(i); }
+  size_t rank0(const size_t& end_position) const {
+    return end_position - rank1(end_position);
+  }
 
   /**
-   * @brief 1-based select of the k-th one.
-   * @return Position of k-th '1' or npos if not found.
+   * @brief 1-based select of the @p target_one_rank-th one.
+   * @return Position of @p target_one_rank-th '1' or npos if not found.
    */
-  size_t select1(size_t k) const {
-    if (k == 0 || num_bits == 0) {
+  size_t select1(size_t target_one_rank) const {
+    if (target_one_rank == 0 || num_bits == 0) {
       return npos;
     }
-    size_t v = 1;
-    if (ones_in_node(v) < k) {
+    size_t node_index = 1;
+    if (ones_in_node(node_index) < target_one_rank) {
       return npos;
     }
-    size_t base = 0;
-    const size_t leaf0 = first_leaf_index();
-    while (v < leaf0) {
-      const size_t Lc = v << 1, Rc = Lc | 1;
-      const uint32_t o_l = ones_in_node(Lc);
-      if (o_l >= k) {
-        v = Lc;
+    size_t segment_base = 0;
+    while (node_index < first_leaf_index) {
+      const size_t left_child = node_index << 1;
+      const size_t right_child = left_child | 1;
+      const uint32_t ones_in_left_child = ones_in_node(left_child);
+      if (ones_in_left_child >= target_one_rank) {
+        node_index = left_child;
       } else {
-        k -= o_l;
-        base += segment_size_bits[Lc];
-        v = Rc;
+        target_one_rank -= ones_in_left_child;
+        segment_base += segment_size_bits[left_child];
+        node_index = right_child;
       }
     }
-    return select1_in_block(base,
-                            std::min(base + segment_size_bits[v], num_bits), k);
+    return select1_in_block(
+        segment_base,
+        std::min(segment_base + segment_size_bits[node_index], num_bits),
+        target_one_rank);
   }
 
   /**
-   * @brief 1-based select of the k-th zero.
-   * @return Position of k-th '0' or npos if not found.
+   * @brief 1-based select of the @p target_zero_rank-th zero.
+   * @return Position of @p target_zero_rank-th '0' or npos if not found.
    */
-  size_t select0(size_t k) const {
-    if (k == 0 || num_bits == 0) {
+  size_t select0(size_t target_zero_rank) const {
+    if (target_zero_rank == 0 || num_bits == 0) {
       return npos;
     }
-    size_t v = 1;
-    const auto zeros = [&](const size_t& x) noexcept {
-      return segment_size_bits[x] - ones_in_node(x);
+    size_t node_index = 1;
+    const auto zeros_in_node = [&](const size_t& node) noexcept {
+      return segment_size_bits[node] - ones_in_node(node);
     };
-    if (zeros(v) < k) {
+    if (zeros_in_node(node_index) < target_zero_rank) {
       return npos;
     }
-    size_t base = 0;
-    const size_t leaf0 = first_leaf_index();
-    while (v < leaf0) {
-      const size_t Lc = v << 1, Rc = Lc | 1;
-      const size_t z_l = zeros(Lc);
-      if (z_l >= k) {
-        v = Lc;
+    size_t segment_base = 0;
+    while (node_index < first_leaf_index) {
+      const size_t left_child = node_index << 1;
+      const size_t right_child = left_child | 1;
+      const size_t zeros_in_left_child = zeros_in_node(left_child);
+      if (zeros_in_left_child >= target_zero_rank) {
+        node_index = left_child;
       } else {
-        k -= z_l;
-        base += segment_size_bits[Lc];
-        v = Rc;
+        target_zero_rank -= zeros_in_left_child;
+        segment_base += segment_size_bits[left_child];
+        node_index = right_child;
       }
     }
-    return select0_in_block(base,
-                            std::min(base + segment_size_bits[v], num_bits), k);
+    return select0_in_block(
+        segment_base,
+        std::min(segment_base + segment_size_bits[node_index], num_bits),
+        target_zero_rank);
   }
 
   /**
-   * @brief Rank of the pattern "10" (starts) within [0, i).
-   * @details Counts p where bit[p]==1 and bit[p+1]==0 with p+1<i.
+   * @brief Rank of the pattern "10" (starts) within [0, @p end_position).
+   * @details Counts p where bit[p]==1 and bit[p+1]==0 with p+1<@p
+   * end_position.
    */
-  size_t rank10(const size_t& i) const {
-    if (i <= 1) {
+  size_t rank10(const size_t& end_position) const {
+    if (end_position <= 1) {
       return 0;
     }
-    const size_t blk = block_of(i - 1);
-    size_t ans = 0;
-    int prev_last = -1;
-
-    if (blk > 0) {
-      const auto nodes = cover_blocks(0, blk - 1);
-      for (const size_t& v : nodes) {
-        ans += node_pattern10_count[v];
-        if (prev_last != -1 && prev_last == 1 && node_first_bit[v] == 0) {
-          ++ans;
+    const size_t block_index = block_of(end_position - 1);
+    size_t pattern_count = 0;
+    int previous_last_bit = -1;
+
+    if (block_index > 0) {
+      const auto covered_nodes = cover_blocks(0, block_index - 1);
+      for (const size_t& node_index : covered_nodes) {
+        pattern_count += node_pattern10_count[node_index];
+        if (previous_last_bit != -1 && previous_last_bit == 1 &&
+            node_first_bit[node_index] == 0) {
+          ++pattern_count;
         }
-        prev_last = node_last_bit[v];
+        previous_last_bit = node_last_bit[node_index];
       }
     }
-    const size_t Lb = blk * block_bits;
-    ans += rr_in_block(Lb, i);
+    const size_t block_begin = block_index * block_bits;
+    pattern_count += rr_in_block(block_begin, end_position);
     // boundary between the last full node and the leaf tail
-    if (blk > 0 && i > Lb && prev_last == 1 && bit(Lb) == 0) {
-      ++ans;
+    if (block_index > 0 && end_position > block_begin &&
+        previous_last_bit == 1 && bit(block_begin) == 0) {
+      ++pattern_count;
     }
-    return ans;
+    return pattern_count;
   }
 
   /**
-   * @brief 1-based select of the k-th "10" start.
+   * @brief 1-based select of the @p target_pattern_rank-th "10" start.
    * @return Position p such that bits[p..p+1]=="10", or npos if not found.
    */
-  size_t select10(size_t k) const {
-    if (k == 0 || num_bits == 0) {
+  size_t select10(size_t target_pattern_rank) const {
+    if (target_pattern_rank == 0 || num_bits == 0) {
       return npos;
     }
-    size_t ind = 1;
-    if (node_pattern10_count[ind] < k) {
+    size_t node_index = 1;
+    if (node_pattern10_count[node_index] < target_pattern_rank) {
       return npos;
     }
-    size_t base = 0;
-    const size_t leaf0 = first_leaf_index();
-    while (ind < leaf0) {
-      const size_t Lc = ind << 1, Rc = Lc | 1;
-      const size_t cross =
-          (node_last_bit[Lc] == 1 && node_first_bit[Rc] == 0) ? 1u : 0u;
-      if (node_pattern10_count[Lc] >= k) {
-        ind = Lc;
+    size_t segment_base = 0;
+    while (node_index < first_leaf_index) {
+      const size_t left_child = node_index << 1, right_child = left_child | 1;
+      const size_t crossing_pattern =
+          (node_last_bit[left_child] == 1 && node_first_bit[right_child] == 0)
+              ? 1u
+              : 0u;
+      if (node_pattern10_count[left_child] >= target_pattern_rank) {
+        node_index = left_child;
         continue;
       }
-      size_t rem = k - node_pattern10_count[Lc];
-      if (cross) {
-        if (rem == 1) {
-          return base + segment_size_bits[Lc] - 1;
+      size_t remaining_rank =
+          target_pattern_rank - node_pattern10_count[left_child];
+      if (crossing_pattern) {
+        if (remaining_rank == 1) {
+          return segment_base + segment_size_bits[left_child] - 1;
         }
-        --rem;
+        --remaining_rank;
       }
-      base += segment_size_bits[Lc];
-      ind = Rc;
-      k = rem;
+      segment_base += segment_size_bits[left_child];
+      node_index = right_child;
+      target_pattern_rank = remaining_rank;
     }
     return select10_in_block(
-        base, std::min(base + segment_size_bits[ind], num_bits), k);
+        segment_base,
+        std::min(segment_base + segment_size_bits[node_index], num_bits),
+        target_pattern_rank);
   }
 
   /**
-   * @brief Prefix excess on [0, i): +1 for '1', −1 for '0'.
+   * @brief Prefix excess on [0, @p end_position): +1 for '1', −1 for '0'.
    */
-  inline int excess(const size_t& i) const {
-    return int64_t(rank1(i)) * 2 - int64_t(i);
+  inline int excess(const size_t& end_position) const {
+    return int64_t(rank1(end_position)) * 2 - int64_t(end_position);
   }
 
   /**
-   * @brief Forward search: first position p ≥ i where excess(p) = excess(i) +
-   * d.
+   * @brief Forward search: first position p ≥ @p start_position where
+   * excess(p) = excess(@p start_position) + @p delta.
    * @details Scans remainder of current leaf, then descends using precomputed
    * bounds. Returns npos if no such position exists.
    */
-  size_t fwdsearch(const size_t& i, const int& d) const {
-    if (i >= num_bits) {
+  size_t fwdsearch(const size_t& start_position, const int& delta) const {
+    if (start_position >= num_bits) {
       return npos;
     }
 
-    const int start_excess = excess(i);
-    const int target = start_excess + d;
-
     // 1) scan the remainder of the current leaf
-    const size_t blk = block_of(i);
-    const size_t Lb = blk * block_bits;
-    const size_t Rb = std::min(num_bits, Lb + block_bits);
-    int cur = start_excess;
-    for (size_t p = i; p < Rb; ++p) {
-      cur += bit(p) ? +1 : -1;
-      if (cur == target) {
-        return p;
-      }
-    }
-
-    // 2) suffix after the leaf: cover full blocks [blk+1 .. leaf_count-1]
-    const int excess_at_Rb = excess(Rb);
-    int need = target - excess_at_Rb;  // target expressed in coordinates of the
-                                       // current node start
-    if (blk + 1 <= (leaf_count ? leaf_count - 1 : 0)) {
-      const auto nodes = cover_blocks(blk + 1, leaf_count - 1);
-      size_t base = (blk + 1) * block_bits;
-      for (const size_t& v : nodes) {
-        if (need == 0) {
-          return base;
+    const size_t leaf_block_index = block_of(start_position);
+    const size_t block_begin = leaf_block_index * block_bits;
+    const size_t block_end = std::min(num_bits, block_begin + block_bits);
+    int leaf_delta = 0;
+    const size_t leaf_result = leaf_fwd_bp_simd(
+        leaf_block_index, block_begin, start_position, delta, leaf_delta);
+    if (leaf_result != npos) {
+      return leaf_result;
+    }
+
+    int remaining_delta = delta - leaf_delta;
+    if (leaf_block_index + 1 < leaf_count) {
+      size_t nodes_buffer[64];
+      const size_t node_count = cover_blocks_collect(
+          leaf_block_index + 1, leaf_count - 1, nodes_buffer);
+      size_t segment_base = (leaf_block_index + 1) * block_bits;
+      for (size_t j = 0; j < node_count; ++j) {
+        const size_t node_index = nodes_buffer[j];
+        if (remaining_delta == 0) {
+          return segment_base;
         }
-        if (node_min_prefix_excess[v] <= need &&
-            need <= node_max_prefix_excess[v]) {
-          return descend_fwd(v, need, base);
+        if (node_min_prefix_excess[node_index] <= remaining_delta &&
+            remaining_delta <= node_max_prefix_excess[node_index]) {
+          return descend_fwd(node_index, remaining_delta, segment_base);
         }
-        need -= node_total_excess[v];
-        base += segment_size_bits[v];
+        remaining_delta -= node_total_excess[node_index];
+        segment_base += segment_size_bits[node_index];
       }
     }
     return npos;
   }
 
   /**
-   * @brief Backward search: last position p ≤ i where excess(p) = excess(i) +
-   * d.
+   * @brief Backward search: last position p ≤ @p start_position where
+   * excess(p) = excess(@p start_position) + @p delta.
    * @details Scans inside the leaf to the left, then climbs to examine left
    * siblings. Returns npos if no such position exists.
    */
-  size_t bwdsearch(const size_t& i, const int& d) const {
-    if (i > num_bits || i == 0) {
+  size_t bwdsearch(const size_t& start_position, const int& delta) const {
+    if (start_position > num_bits || start_position == 0) {
       return npos;
     }
-    const int start_excess = excess(i);
-    const int target = start_excess + d;
 
     // 1) scan inside the block
-    const size_t blk = block_of(i - 1);
-    const size_t Lb = blk * block_bits;
-    int cur = start_excess;
-    for (size_t p = i; p > Lb;) {
-      --p;
-      cur += bit(p) ? -1 : +1;
-      if (cur == target) {
-        return p;
-      }
-      if (p == 0) {
-        break;
-      }
-    }
-    const int excess_at_Lb = excess(Lb);
-    if (Lb < i && excess_at_Lb == target) {
-      return Lb;
-    }
-
-    // 2) climb up
-    int need = target - excess_at_Lb;
-    size_t v = leaf_index_of(Lb);
-    size_t base = Lb;
-    while (v > 1) {
-      if (v & 1) {                 // v is the right child
-        const size_t sib = v ^ 1;  // left sibling
-        const size_t border =
-            base;  // right border of the sibling (== start(v))
-        const int need_node =
-            need + node_total_excess[sib];  // target in coordinates relative to
-                                            // the start of sib
-        const bool allow_rb = (border != i);  // j must be < i
+    const size_t leaf_block_index = block_of(start_position - 1);
+    const size_t block_begin = leaf_block_index * block_bits;
+    int leaf_delta = 0;  // excess(start_position) - excess(block_begin)
+    const size_t leaf_result = leaf_bwd_bp_simd(
+        leaf_block_index, block_begin, start_position, delta, leaf_delta);
+    if (leaf_result != npos) {
+      return leaf_result;
+    }
+
+    // need = target - excess(block_begin) = excess(start_position) + delta -
+    // excess(block_begin) = leaf_delta + delta
+    int remaining_delta = leaf_delta + delta;
+    size_t node_index = leaf_index_of(block_begin);
+    size_t segment_base = block_begin;
+    while (node_index > 1) {
+      if (node_index & 1) {  // node_index is the right child
+        const size_t sibling_index = node_index ^ 1;  // left sibling
+        const size_t sibling_border =
+            segment_base;  // right border of the sibling (== start(node_index))
+        const int needed_inside_sibling =
+            remaining_delta +
+            node_total_excess[sibling_index];  // target in coordinates relative
+                                               // to the start of sibling
+        const bool allow_right_border =
+            (sibling_border != start_position);  // j must be < start_position
 
         // try inside the sibling, but return only if a position is found
-        if (need_node == 0 || (node_min_prefix_excess[sib] <= need_node &&
-                               need_node <= node_max_prefix_excess[sib])) {
-          const size_t ans = descend_bwd(sib, border - segment_size_bits[sib],
-                                         need_node, border, allow_rb);
-          if (ans != npos) {
-            return ans;
+        if (needed_inside_sibling == 0 ||
+            (node_min_prefix_excess[sibling_index] <= needed_inside_sibling &&
+             needed_inside_sibling <= node_max_prefix_excess[sibling_index])) {
+          const size_t result = descend_bwd(
+              sibling_index, sibling_border - segment_size_bits[sibling_index],
+              needed_inside_sibling, sibling_border, allow_right_border);
+          if (result != npos) {
+            return result;
           }
         }
         // junction between children is a separate branch (allowed only if < i)
-        if (need_node == node_total_excess[sib] && border < i) {
-          return border;
+        if (needed_inside_sibling == node_total_excess[sibling_index] &&
+            sibling_border < start_position) {
+          return sibling_border;
         }
 
         // stepped over the sibling, shifted the zero point of the coordinates
-        need += node_total_excess[sib];
-        base -= segment_size_bits[sib];
+        remaining_delta += node_total_excess[sibling_index];
+        segment_base -= segment_size_bits[sibling_index];
       }
-      v >>= 1;
+      node_index >>= 1;
     }
     return npos;
   }
 
   /**
-   * @brief Position of the first minimum of excess on [i, j] (inclusive).
+   * @brief Position of the first minimum of excess on [@p range_begin, @p
+   * range_end]
+   * (inclusive).
    * @return Position of first occurrence of minimum, or npos on invalid range.
    */
-  size_t range_min_query_pos(const size_t& i, const size_t& j) const {
-    if (i > j || j >= num_bits) {
+  size_t range_min_query_pos(const size_t& range_begin,
+                             const size_t& range_end) const {
+    if (range_begin > range_end || range_end >= num_bits) {
       return npos;
     }
 
-    const size_t blk_i = block_of(i);
-    const size_t Lbi = blk_i * block_bits;
-    const size_t Rbi = std::min(num_bits, Lbi + block_bits);
-    const size_t blk_j = block_of(j);
-    const size_t Lbj = blk_j * block_bits;
+    const size_t begin_block_index = block_of(range_begin);
+    const size_t begin_block_start = begin_block_index * block_bits;
+    const size_t begin_block_end =
+        std::min(num_bits, begin_block_start + block_bits);
+    const size_t end_block_index = block_of(range_end);
+    const size_t end_block_start = end_block_index * block_bits;
 
-    int best_val = INT_MAX;
-    size_t best_pos = npos;
-    size_t chosen_node = 0;
-    int pref = 0, pref_at_choice = 0;
+    int best_value = INT_MAX;
+    size_t best_position = npos;
+    size_t chosen_node_index = 0;
+    int prefix_excess = 0, prefix_at_choice = 0;
 
     // prefix
-    int mn_first = INT_MAX;
-    size_t first_pos = npos;
-    const size_t end_first = std::min(j, (size_t)(Rbi ? Rbi - 1 : 0));
-    if (i <= end_first) {
-      first_min_value_pos8(i, end_first, mn_first, first_pos);
-      pref = (int64_t)rank1_in_block(i, end_first + 1) * 2 -
-             int64_t(end_first + 1 - i);
-      best_val = mn_first;
-      best_pos = first_pos;
-      chosen_node = 0;
+    int min_prefix_first_chunk = INT_MAX;
+    size_t first_chunk_position = npos;
+    const size_t end_of_first_chunk = std::min(
+        range_end, (size_t)(begin_block_end ? begin_block_end - 1 : 0));
+    if (range_begin <= end_of_first_chunk) {
+      first_min_value_pos8(range_begin, end_of_first_chunk,
+                           min_prefix_first_chunk, first_chunk_position);
+      prefix_excess =
+          (int64_t)rank1_in_block(range_begin, end_of_first_chunk + 1) * 2 -
+          int64_t(end_of_first_chunk + 1 - range_begin);
+      best_value = min_prefix_first_chunk;
+      best_position = first_chunk_position;
+      chosen_node_index = 0;
     }
 
     // middle
-    const size_t leaf0 = first_leaf_index();
-    if (blk_i + 1 <= blk_j - 1) {
-      size_t l = leaf0 + (blk_i + 1);
-      size_t r = leaf0 + (blk_j - 1);
-      size_t Rnodes[64];
-      int rn = 0;
-
-      while (l <= r) {
-        if (l & 1) {
-          const size_t v = l++;
-          const int cand = pref + node_min_prefix_excess[v];
-          if (cand < best_val) {
-            best_val = cand;
-            best_pos = npos;
-            chosen_node = v;
-            pref_at_choice = pref;
+    if (begin_block_index + 1 <= end_block_index - 1) {
+      size_t left_index = first_leaf_index + (begin_block_index + 1);
+      size_t right_index = first_leaf_index + (end_block_index - 1);
+      size_t right_nodes[64];
+      int right_nodes_count = 0;
+
+      while (left_index <= right_index) {
+        if (left_index & 1) {
+          const size_t node_index = left_index++;
+          const int candidate =
+              prefix_excess + node_min_prefix_excess[node_index];
+          if (candidate < best_value) {
+            best_value = candidate;
+            best_position = npos;
+            chosen_node_index = node_index;
+            prefix_at_choice = prefix_excess;
           }
-          pref += node_total_excess[v];
+          prefix_excess += node_total_excess[node_index];
         }
-        if ((r & 1) == 0) {
-          Rnodes[rn++] = r--;
+        if ((right_index & 1) == 0) {
+          right_nodes[right_nodes_count++] = right_index--;
         }
-        l >>= 1;
-        r >>= 1;
-      }
-      while (rn--) {
-        const size_t v = Rnodes[rn];
-        const int cand = pref + node_min_prefix_excess[v];
-        if (cand < best_val) {
-          best_val = cand;
-          best_pos = npos;
-          chosen_node = v;
-          pref_at_choice = pref;
+        left_index >>= 1;
+        right_index >>= 1;
+      }
+      while (right_nodes_count--) {
+        const size_t node_index = right_nodes[right_nodes_count];
+        const int candidate =
+            prefix_excess + node_min_prefix_excess[node_index];
+        if (candidate < best_value) {
+          best_value = candidate;
+          best_position = npos;
+          chosen_node_index = node_index;
+          prefix_at_choice = prefix_excess;
         }
-        pref += node_total_excess[v];
+        prefix_excess += node_total_excess[node_index];
       }
     }
 
     // tail
-    if (blk_j != blk_i) {
-      int mn_last;
-      size_t last_pos;
-      first_min_value_pos8(Lbj, j, mn_last, last_pos);
-      const int cand = pref + mn_last;
-      if (cand < best_val) {
-        best_val = cand;
-        best_pos = last_pos;
-        chosen_node = 0;
+    if (end_block_index != begin_block_index) {
+      int min_prefix_last_chunk;
+      size_t last_chunk_position;
+      first_min_value_pos8(end_block_start, range_end, min_prefix_last_chunk,
+                           last_chunk_position);
+      const int candidate = prefix_excess + min_prefix_last_chunk;
+      if (candidate < best_value) {
+        best_value = candidate;
+        best_position = last_chunk_position;
+        chosen_node_index = 0;
       }
     }
 
-    if (best_pos != npos) {
-      return best_pos;
+    if (best_position != npos) {
+      return best_position;
     }
 
-    return descend_first_min(chosen_node, best_val - pref_at_choice,
-                             node_base(chosen_node));
+    return descend_first_min(chosen_node_index, best_value - prefix_at_choice,
+                             node_base(chosen_node_index));
   }
 
   /**
-   * @brief Value of the minimum prefix excess on [i, j] relative to i.
-   * @details Equivalent to min_{t in [i..j]} (excess(t+1) - excess(i)).
+   * @brief Value of the minimum prefix excess on [@p range_begin, @p range_end]
+   * relative to @p range_begin.
+   * @details Equivalent to min_{t in [@p range_begin..@p range_end]}
+   * (excess(t+1) - excess(@p range_begin)).
    */
-  int range_min_query_val(const size_t& i, const size_t& j) const {
-    if (i > j || j >= num_bits) {
+  int range_min_query_val(const size_t& range_begin,
+                          const size_t& range_end) const {
+    if (range_begin > range_end || range_end >= num_bits) {
       return 0;
     }
-    size_t p = range_min_query_pos(i, j);
-    if (p == npos) {
+    size_t min_position = range_min_query_pos(range_begin, range_end);
+    if (min_position == npos) {
       return 0;
     }
-    return excess(p + 1) - excess(i);
+    return excess(min_position + 1) - excess(range_begin);
   }
 
   /**
-   * @brief Position of the first maximum of excess on [i, j] (inclusive).
+   * @brief Position of the first maximum of excess on [@p range_begin, @p
+   * range_end]
+   * (inclusive).
    * @return Position of first occurrence of maximum, or npos on invalid range.
    */
-  size_t range_max_query_pos(const size_t& i, const size_t& j) const {
-    if (i > j || j >= num_bits) {
+  size_t range_max_query_pos(const size_t& range_begin,
+                             const size_t& range_end) const {
+    if (range_begin > range_end || range_end >= num_bits) {
       return npos;
     }
 
-    const size_t blk_i = block_of(i);
-    const size_t Lbi = blk_i * block_bits;
-    const size_t Rbi = std::min(num_bits, Lbi + block_bits);
-    const size_t blk_j = block_of(j);
-    const size_t Lbj = blk_j * block_bits;
+    const size_t begin_block_index = block_of(range_begin);
+    const size_t begin_block_start = begin_block_index * block_bits;
+    const size_t begin_block_end =
+        std::min(num_bits, begin_block_start + block_bits);
+    const size_t end_block_index = block_of(range_end);
+    const size_t end_block_start = end_block_index * block_bits;
 
-    int best_val = INT_MIN;
-    size_t best_pos = npos;
-    size_t chosen_node = 0;
-    int pref = 0, pref_at_choice = 0;
+    int best_value = INT_MIN;
+    size_t best_position = npos;
+    size_t chosen_node_index = 0;
+    int prefix_excess = 0, prefix_at_choice = 0;
 
     // prefix
-    int mx_first = INT_MIN;
-    size_t first_pos = npos;
-    const size_t end_first = std::min(j, (size_t)(Rbi ? Rbi - 1 : 0));
-    if (i <= end_first) {
-      first_max_value_pos8(i, end_first, mx_first, first_pos);
-      pref = (int64_t)rank1_in_block(i, end_first + 1) * 2 -
-             int64_t(end_first + 1 - i);
-      best_val = mx_first;
-      best_pos = first_pos;
-      chosen_node = 0;
+    int max_prefix_first_chunk = INT_MIN;
+    size_t first_chunk_position = npos;
+    const size_t end_of_first_chunk = std::min(
+        range_end, (size_t)(begin_block_end ? begin_block_end - 1 : 0));
+    if (range_begin <= end_of_first_chunk) {
+      first_max_value_pos8(range_begin, end_of_first_chunk,
+                           max_prefix_first_chunk, first_chunk_position);
+      prefix_excess =
+          (int64_t)rank1_in_block(range_begin, end_of_first_chunk + 1) * 2 -
+          int64_t(end_of_first_chunk + 1 - range_begin);
+      best_value = max_prefix_first_chunk;
+      best_position = first_chunk_position;
+      chosen_node_index = 0;
     }
 
     // middle
-    const size_t leaf0 = first_leaf_index();
-    if (blk_i + 1 <= blk_j - 1) {
-      size_t l = leaf0 + (blk_i + 1);
-      size_t r = leaf0 + (blk_j - 1);
-      size_t Rnodes[64];
-      int rn = 0;
-
-      while (l <= r) {
-        if (l & 1) {
-          const size_t v = l++;
-          const int cand = pref + node_max_prefix_excess[v];
-          if (cand > best_val) {
-            best_val = cand;
-            best_pos = npos;
-            chosen_node = v;
-            pref_at_choice = pref;
+    if (begin_block_index + 1 <= end_block_index - 1) {
+      size_t left_index = first_leaf_index + (begin_block_index + 1);
+      size_t right_index = first_leaf_index + (end_block_index - 1);
+      size_t right_nodes[64];
+      int right_nodes_count = 0;
+
+      while (left_index <= right_index) {
+        if (left_index & 1) {
+          const size_t node_index = left_index++;
+          const int candidate =
+              prefix_excess + node_max_prefix_excess[node_index];
+          if (candidate > best_value) {
+            best_value = candidate;
+            best_position = npos;
+            chosen_node_index = node_index;
+            prefix_at_choice = prefix_excess;
           }
-          pref += node_total_excess[v];
+          prefix_excess += node_total_excess[node_index];
         }
-        if ((r & 1) == 0) {
-          Rnodes[rn++] = r--;
+        if ((right_index & 1) == 0) {
+          right_nodes[right_nodes_count++] = right_index--;
         }
-        l >>= 1;
-        r >>= 1;
-      }
-      while (rn--) {
-        const size_t v = Rnodes[rn];
-        const int cand = pref + node_max_prefix_excess[v];
-        if (cand > best_val) {
-          best_val = cand;
-          best_pos = npos;
-          chosen_node = v;
-          pref_at_choice = pref;
+        left_index >>= 1;
+        right_index >>= 1;
+      }
+      while (right_nodes_count--) {
+        const size_t node_index = right_nodes[right_nodes_count];
+        const int candidate =
+            prefix_excess + node_max_prefix_excess[node_index];
+        if (candidate > best_value) {
+          best_value = candidate;
+          best_position = npos;
+          chosen_node_index = node_index;
+          prefix_at_choice = prefix_excess;
         }
-        pref += node_total_excess[v];
+        prefix_excess += node_total_excess[node_index];
       }
     }
 
     // tail
-    if (blk_j != blk_i) {
-      int mx_last;
-      size_t last_pos;
-      first_max_value_pos8(Lbj, j, mx_last, last_pos);
-      const int cand = pref + mx_last;
-      if (cand > best_val) {
-        best_val = cand;
-        best_pos = last_pos;
-        chosen_node = 0;
+    if (end_block_index != begin_block_index) {
+      int max_prefix_last_chunk;
+      size_t last_chunk_position;
+      first_max_value_pos8(end_block_start, range_end, max_prefix_last_chunk,
+                           last_chunk_position);
+      const int candidate = prefix_excess + max_prefix_last_chunk;
+      if (candidate > best_value) {
+        best_value = candidate;
+        best_position = last_chunk_position;
+        chosen_node_index = 0;
       }
     }
 
-    if (best_pos != npos) {
-      return best_pos;
+    if (best_position != npos) {
+      return best_position;
     }
 
-    return descend_first_max(chosen_node, best_val - pref_at_choice,
-                             node_base(chosen_node));
+    return descend_first_max(chosen_node_index, best_value - prefix_at_choice,
+                             node_base(chosen_node_index));
   }
 
   /**
-   * @brief Value of the maximum prefix excess on [i, j] relative to i.
+   * @brief Value of the maximum prefix excess on [@p range_begin, @p
+   * range_end] relative to @p range_begin.
    */
-  int range_max_query_val(const size_t& i, const size_t& j) const {
-    if (i > j || j >= num_bits) {
+  int range_max_query_val(const size_t& range_begin,
+                          const size_t& range_end) const {
+    if (range_begin > range_end || range_end >= num_bits) {
       return 0;
     }
-    size_t p = range_max_query_pos(i, j);
-    if (p == npos) {
+    size_t max_position = range_max_query_pos(range_begin, range_end);
+    if (max_position == npos) {
       return 0;
     }
-    return excess(p + 1) - excess(i);
+    return excess(max_position + 1) - excess(range_begin);
   }
 
   /**
-   * @brief How many times the minimum prefix excess occurs on [i, j].
+   * @brief How many times the minimum prefix excess occurs on [@p range_begin,
+   * @p range_end].
    */
-  size_t mincount(const size_t& i, const size_t& j) const {
-    if (i > j || j >= num_bits) {
+  size_t mincount(const size_t& range_begin, const size_t& range_end) const {
+    if (range_begin > range_end || range_end >= num_bits) {
       return 0;
     }
 
-    const size_t blk_i = block_of(i);
-    const size_t Lbi = blk_i * block_bits;
-    const size_t Rbi = std::min(num_bits, Lbi + block_bits);
-    const size_t blk_j = block_of(j);
-    const size_t Lbj = blk_j * block_bits;
+    const size_t begin_block_index = block_of(range_begin);
+    const size_t begin_block_start = begin_block_index * block_bits;
+    const size_t begin_block_end =
+        std::min(num_bits, begin_block_start + block_bits);
+    const size_t end_block_index = block_of(range_end);
+    const size_t end_block_start = end_block_index * block_bits;
 
-    int best_val = INT_MAX;
-    size_t cnt = 0;
-    int pref = 0;
+    int best_value = INT_MAX;
+    size_t min_count = 0;
+    int prefix_excess = 0;
 
     // first chunk
     {
-      int cur = 0, mn = INT_MAX, c = 0;
-      const size_t end = std::min(j, Rbi - 1);
-      for (size_t p = i; p <= end; ++p) {
-        cur += bit(p) ? +1 : -1;
-        if (cur < mn) {
-          mn = cur;
-          c = 1;
-        } else if (cur == mn) {
-          ++c;
+      int current_excess = 0, min_value = INT_MAX, local_count = 0;
+      const size_t end_of_first_chunk =
+          std::min(range_end, begin_block_end - 1);
+      for (size_t position = range_begin; position <= end_of_first_chunk;
+           ++position) {
+        current_excess += bit(position) ? +1 : -1;
+        if (current_excess < min_value) {
+          min_value = current_excess;
+          local_count = 1;
+        } else if (current_excess == min_value) {
+          ++local_count;
         }
       }
-      best_val = mn;
-      cnt = c;
-      pref = cur;  // offset toward the middle
+      best_value = min_value;
+      min_count = local_count;
+      prefix_excess = current_excess;  // offset toward the middle
     }
 
     // middle
-    if (blk_i + 1 <= blk_j - 1) {
-      const auto mids = cover_blocks(blk_i + 1, blk_j - 1);
-      for (const size_t& v : mids) {
-        const int cand = pref + node_min_prefix_excess[v];
-        if (cand < best_val) {
-          best_val = cand;
-          cnt = node_min_count[v];
-        } else if (cand == best_val) {
-          cnt += node_min_count[v];
+    if (begin_block_index + 1 <= end_block_index - 1) {
+      const auto middle_nodes =
+          cover_blocks(begin_block_index + 1, end_block_index - 1);
+      for (const size_t& node_index : middle_nodes) {
+        const int candidate =
+            prefix_excess + node_min_prefix_excess[node_index];
+        if (candidate < best_value) {
+          best_value = candidate;
+          min_count = node_min_count[node_index];
+        } else if (candidate == best_value) {
+          min_count += node_min_count[node_index];
         }
-        pref += node_total_excess[v];
+        prefix_excess += node_total_excess[node_index];
       }
     }
 
     // last chunk
-    if (blk_j != blk_i) {
-      int cur = 0, mn = INT_MAX, c = 0;
-      for (size_t p = Lbj; p <= j; ++p) {
-        cur += bit(p) ? +1 : -1;
-        if (cur < mn) {
-          mn = cur;
-          c = 1;
-        } else if (cur == mn) {
-          ++c;
+    if (end_block_index != begin_block_index) {
+      int current_excess = 0, min_value = INT_MAX, local_count = 0;
+      for (size_t position = end_block_start; position <= range_end;
+           ++position) {
+        current_excess += bit(position) ? +1 : -1;
+        if (current_excess < min_value) {
+          min_value = current_excess;
+          local_count = 1;
+        } else if (current_excess == min_value) {
+          ++local_count;
         }
       }
-      const int cand = pref + mn;
-      if (cand < best_val) {
-        best_val = cand;
-        cnt = c;
-      } else if (cand == best_val) {
-        cnt += c;
+      const int candidate = prefix_excess + min_value;
+      if (candidate < best_value) {
+        best_value = candidate;
+        min_count = local_count;
+      } else if (candidate == best_value) {
+        min_count += local_count;
       }
     }
-    return cnt;
+    return min_count;
   }
 
   /**
-   * @brief Position of the q-th (1-based) occurrence of the minimum on [i, j].
-   * @return Position or npos if q exceeds the number of minima.
+   * @brief Position of the @p target_min_rank-th (1-based) occurrence of the
+   * minimum on [@p range_begin, @p range_end].
+   * @return Position or npos if @p target_min_rank exceeds the number of
+   * minima.
    */
-  size_t minselect(const size_t& i, const size_t& j, size_t q) const {
-    if (i > j || j >= num_bits || q == 0) {
+  size_t minselect(const size_t& range_begin,
+                   const size_t& range_end,
+                   size_t target_min_rank) const {
+    if (range_begin > range_end || range_end >= num_bits ||
+        target_min_rank == 0) {
       return npos;
     }
 
-    const size_t blk_i = block_of(i);
-    const size_t Lbi = blk_i * block_bits;
-    const size_t Rbi = std::min(num_bits, Lbi + block_bits);
-    const size_t blk_j = block_of(j);
-    const size_t Lbj = blk_j * block_bits;
+    const size_t begin_block_index = block_of(range_begin);
+    const size_t begin_block_start = begin_block_index * block_bits;
+    const size_t begin_block_end =
+        std::min(num_bits, begin_block_start + block_bits);
+    const size_t end_block_index = block_of(range_end);
+    const size_t end_block_start = end_block_index * block_bits;
 
     // prefix
-    const size_t end_first = std::min(j, Rbi - 1);
-    int cur_first = 0, mn_first = 0;
-    uint32_t c_first = 0;
-
-    if (i <= end_first) {
-      scan_range_min_count8(i, end_first, cur_first, mn_first, c_first);
+    const size_t end_of_first_chunk = std::min(range_end, begin_block_end - 1);
+    int current_first_chunk_excess = 0, min_first_chunk = 0;
+    uint32_t count_first_chunk = 0;
+
+    if (range_begin <= end_of_first_chunk) {
+      scan_range_min_count8(range_begin, end_of_first_chunk,
+                            current_first_chunk_excess, min_first_chunk,
+                            count_first_chunk);
     } else {
-      cur_first = 0;
-      mn_first = INT_MAX;
-      c_first = 0;
+      current_first_chunk_excess = 0;
+      min_first_chunk = INT_MAX;
+      count_first_chunk = 0;
     }
 
-    int best_val = (mn_first == INT_MAX ? INT_MAX : mn_first);
-    size_t total_cnt = (mn_first == INT_MAX ? 0u : (size_t)c_first);
-    int pref = cur_first;  // offset for middle
+    int best_value = (min_first_chunk == INT_MAX ? INT_MAX : min_first_chunk);
+    size_t total_count =
+        (min_first_chunk == INT_MAX ? 0u : (size_t)count_first_chunk);
+    int prefix_excess = current_first_chunk_excess;  // offset for middle
 
-    const size_t leaf0 = first_leaf_index();
-    size_t l = leaf0 + blk_i + 1;
-    size_t r = leaf0 + blk_j - 1;
-    size_t Rnodes[64];
-    int rn = 0;
+    size_t left_index = first_leaf_index + begin_block_index + 1;
+    size_t right_index = first_leaf_index + end_block_index - 1;
+    size_t right_nodes[64];
+    int right_nodes_count = 0;
 
     // middle
-    if (blk_i + 1 <= blk_j - 1) {
-      while (l <= r) {
-        if (l & 1) {
-          const int cand = pref + node_min_prefix_excess[l];
-          if (cand < best_val) {
-            best_val = cand;
-            total_cnt = node_min_count[l];
-          } else if (cand == best_val) {
-            total_cnt += node_min_count[l];
+    if (begin_block_index + 1 <= end_block_index - 1) {
+      while (left_index <= right_index) {
+        if (left_index & 1) {
+          const int candidate =
+              prefix_excess + node_min_prefix_excess[left_index];
+          if (candidate < best_value) {
+            best_value = candidate;
+            total_count = node_min_count[left_index];
+          } else if (candidate == best_value) {
+            total_count += node_min_count[left_index];
           }
-          pref += node_total_excess[l++];
+          prefix_excess += node_total_excess[left_index++];
         }
-        if ((r & 1) == 0) {
-          Rnodes[rn++] = r--;
+        if ((right_index & 1) == 0) {
+          right_nodes[right_nodes_count++] = right_index--;
         }
-        l >>= 1;
-        r >>= 1;
-      }
-      while (rn--) {
-        const size_t v = Rnodes[rn];
-        const int cand = pref + node_min_prefix_excess[v];
-        if (cand < best_val) {
-          best_val = cand;
-          total_cnt = node_min_count[v];
-        } else if (cand == best_val) {
-          total_cnt += node_min_count[v];
+        left_index >>= 1;
+        right_index >>= 1;
+      }
+      while (right_nodes_count--) {
+        const size_t node_index = right_nodes[right_nodes_count];
+        const int candidate =
+            prefix_excess + node_min_prefix_excess[node_index];
+        if (candidate < best_value) {
+          best_value = candidate;
+          total_count = node_min_count[node_index];
+        } else if (candidate == best_value) {
+          total_count += node_min_count[node_index];
         }
-        pref += node_total_excess[v];
+        prefix_excess += node_total_excess[node_index];
       }
     }
 
     // tail
-    int cur_last = 0, mn_last = INT_MAX;
-    uint32_t c_last = 0;
-    if (blk_j != blk_i) {
-      scan_range_min_count8(Lbj, j, cur_last, mn_last, c_last);
-      const int cand = pref + mn_last;
-      if (cand < best_val) {
-        best_val = cand;
-        total_cnt = c_last;
-      } else if (cand == best_val) {
-        total_cnt += c_last;
-      }
-    }
-
-    if (q > total_cnt) {
+    int current_last_chunk_excess = 0, min_last_chunk = INT_MAX;
+    uint32_t count_last_chunk = 0;
+    if (end_block_index != begin_block_index) {
+      scan_range_min_count8(end_block_start, range_end,
+                            current_last_chunk_excess, min_last_chunk,
+                            count_last_chunk);
+      const int candidate = prefix_excess + min_last_chunk;
+      if (candidate < best_value) {
+        best_value = candidate;
+        total_count = count_last_chunk;
+      } else if (candidate == best_value) {
+        total_count += count_last_chunk;
+      }
+    }
+
+    if (target_min_rank > total_count) {
       return npos;
     }
 
     // prefix
-    if (mn_first == best_val && c_first) {
-      if (q <= c_first) {
-        return qth_min_in_block(i, end_first, q);
+    if (min_first_chunk == best_value && count_first_chunk) {
+      if (target_min_rank <= count_first_chunk) {
+        return qth_min_in_block(range_begin, end_of_first_chunk,
+                                target_min_rank);
       }
-      q -= c_first;
+      target_min_rank -= count_first_chunk;
     }
 
     // middle
-    pref = cur_first;
-    if (blk_i + 1 <= blk_j - 1) {
-      l = leaf0 + (blk_i + 1);
-      r = leaf0 + (blk_j - 1);
-      rn = 0;
-      while (l <= r) {
-        if (l & 1) {
-          const size_t v = l++;
-          const int cand = pref + node_min_prefix_excess[v];
-          if (cand == best_val) {
-            if (q <= node_min_count[v]) {
-              return descend_qth_min(v, best_val - pref, q, node_base(v));
+    prefix_excess = current_first_chunk_excess;
+    if (begin_block_index + 1 <= end_block_index - 1) {
+      left_index = first_leaf_index + (begin_block_index + 1);
+      right_index = first_leaf_index + (end_block_index - 1);
+      right_nodes_count = 0;
+      while (left_index <= right_index) {
+        if (left_index & 1) {
+          const size_t node_index = left_index++;
+          const int candidate =
+              prefix_excess + node_min_prefix_excess[node_index];
+          if (candidate == best_value) {
+            if (target_min_rank <= node_min_count[node_index]) {
+              return descend_qth_min(node_index, best_value - prefix_excess,
+                                     target_min_rank, node_base(node_index));
             }
-            q -= node_min_count[v];
+            target_min_rank -= node_min_count[node_index];
           }
-          pref += node_total_excess[v];
+          prefix_excess += node_total_excess[node_index];
         }
-        if (!(r & 1)) {
-          Rnodes[rn++] = r--;
+        if (!(right_index & 1)) {
+          right_nodes[right_nodes_count++] = right_index--;
         }
-        l >>= 1;
-        r >>= 1;
-      }
-      while (rn--) {
-        const size_t v = Rnodes[rn];
-        const int cand = pref + node_min_prefix_excess[v];
-        if (cand == best_val) {
-          if (q <= node_min_count[v]) {
-            return descend_qth_min(v, best_val - pref, q, node_base(v));
+        left_index >>= 1;
+        right_index >>= 1;
+      }
+      while (right_nodes_count--) {
+        const size_t node_index = right_nodes[right_nodes_count];
+        const int candidate =
+            prefix_excess + node_min_prefix_excess[node_index];
+        if (candidate == best_value) {
+          if (target_min_rank <= node_min_count[node_index]) {
+            return descend_qth_min(node_index, best_value - prefix_excess,
+                                   target_min_rank, node_base(node_index));
           }
-          q -= node_min_count[v];
+          target_min_rank -= node_min_count[node_index];
         }
-        pref += node_total_excess[v];
+        prefix_excess += node_total_excess[node_index];
       }
     }
 
     // tail
-    if (blk_j != blk_i && (pref + mn_last) == best_val) {
-      return qth_min_in_block(Lbj, j, q);
+    if (end_block_index != begin_block_index &&
+        (prefix_excess + min_last_chunk) == best_value) {
+      return qth_min_in_block(end_block_start, range_end, target_min_rank);
     }
 
     return npos;
@@ -836,232 +902,249 @@ class RmMTree {
   // ----- parentheses navigation (BP) -----
 
   /**
-   * @brief close(i): matching ')' for '(' at i.
+   * @brief close(@p open_position): matching ')' for '(' at @p open_position.
    * @return Position of matching ')', or npos.
    */
-  inline size_t close(const size_t& i) const {
-    if (i >= num_bits) {
+  inline size_t close(const size_t& open_position) const {
+    if (open_position >= num_bits) {
       return npos;
     }
-    return fwdsearch(i, -1);
+    return fwdsearch(open_position, -1);
   }
 
   /**
-   * @brief open(i): matching '(' for ')' at i.
+   * @brief open(@p close_position): matching '(' for ')' at @p close_position.
    * @return Position of matching '(', or npos.
    */
-  inline size_t open(const size_t& i) const {
+  inline size_t open(const size_t& close_position) const {
     // bwdsearch allows i in [1..num_bits]
-    if (i == 0 || i > num_bits) {
+    if (close_position == 0 || close_position > num_bits) {
       return npos;
     }
-    const size_t r = bwdsearch(i, 0);
-    return (r == npos ? npos : r + 1);
+    const size_t result = bwdsearch(close_position, 0);
+    return (result == npos ? npos : result + 1);
   }
 
   /**
-   * @brief enclose(i): opening '(' that strictly encloses position i.
+   * @brief enclose(@p position): opening '(' that strictly encloses @p
+   * position.
    * @return Position of enclosing '(', or npos.
    */
-  inline size_t enclose(const size_t& i) const {
-    if (i == 0 || i > num_bits) {
+  inline size_t enclose(const size_t& position) const {
+    if (position == 0 || position > num_bits) {
       return npos;
     }
-    const size_t r = bwdsearch(i, -2);
-    return (r == npos ? npos : r + 1);
+    const size_t result = bwdsearch(position, -2);
+    return (result == npos ? npos : result + 1);
   }
 
  private:
   /**
    * @brief Count "10" occurrences inside a 64-bit slice of given logical
-   * length.
-   * @details Only positions fully inside the slice are counted.
+   * length @p length.
+   * @details Only positions fully inside the slice @p slice are counted.
    */
   static inline size_t pop10_in_slice64(const std::uint64_t& slice,
-                                        const int& len) noexcept {
-    if (len <= 1) {
+                                        const int& length) noexcept {
+    if (length <= 1) {
       return 0;
     }
-    std::uint64_t P = slice & ~(slice >> 1);  // candidates for "10"
-    if (len < 64) {
-      P &= ((std::uint64_t(1) << (len - 1)) - 1);
+    std::uint64_t pattern_mask = slice & ~(slice >> 1);  // candidates for "10"
+    if (length < 64) {
+      pattern_mask &= ((std::uint64_t(1) << (length - 1)) - 1);
     } else {
-      P &= 0x7FFFFFFFFFFFFFFFull;
+      pattern_mask &= 0x7FFFFFFFFFFFFFFFull;
     }
-    return (size_t)std::popcount(P);
+    return (size_t)std::popcount(pattern_mask);
   }
 
   /**
-   * @brief Rank of ones within [Lb, Rb).
-   * @details Works on word boundaries; Rb may equal Lb.
+   * @brief Rank of ones within [@p block_begin, @p block_end).
+   * @details Works on word boundaries; @p block_end may equal @p block_begin.
    */
-  size_t rank1_in_block(const size_t& Lb, const size_t& Rb) const noexcept {
-    if (Rb <= Lb) {
+  size_t rank1_in_block(const size_t& block_begin,
+                        const size_t& block_end) const noexcept {
+    if (block_end <= block_begin) {
       return 0;
     }
-    size_t w_l = Lb >> 6;
-    const size_t w_r = Rb >> 6;
-    size_t off_l = Lb & 63;
-    const size_t off_r = Rb & 63;
-    size_t cnt = 0;
-    if (w_l == w_r) {
+    size_t left_word_index = block_begin >> 6;
+    const size_t right_word_index = block_end >> 6;
+    size_t left_offset = block_begin & 63;
+    const size_t right_offset = block_end & 63;
+    size_t count = 0;
+    if (left_word_index == right_word_index) {
       const std::uint64_t mask =
-          ((off_r == 0) ? 0 : ((std::uint64_t(1) << off_r) - 1)) &
-          (~std::uint64_t(0) << off_l);
-      return (size_t)std::popcount(bits[w_l] & mask);
+          ((right_offset == 0) ? 0 : ((std::uint64_t(1) << right_offset) - 1)) &
+          (~std::uint64_t(0) << left_offset);
+      return (size_t)std::popcount(bits[left_word_index] & mask);
     }
-    if (off_l) {
-      cnt += (size_t)std::popcount(bits[w_l] & (~std::uint64_t(0) << off_l));
-      ++w_l;
+    if (left_offset) {
+      count += (size_t)std::popcount(bits[left_word_index] &
+                                     (~std::uint64_t(0) << left_offset));
+      ++left_word_index;
     }
-    while (w_l < w_r) {
-      cnt += (size_t)std::popcount(bits[w_l]);
-      ++w_l;
+    while (left_word_index < right_word_index) {
+      count += (size_t)std::popcount(bits[left_word_index]);
+      ++left_word_index;
     }
-    if (off_r) {
-      cnt +=
-          (size_t)std::popcount(bits[w_r] & ((std::uint64_t(1) << off_r) - 1));
+    if (right_offset) {
+      count += (size_t)std::popcount(bits[right_word_index] &
+                                     ((std::uint64_t(1) << right_offset) - 1));
     }
-    return cnt;
+    return count;
   }
 
   /**
-   * @brief Count "10" starts within [Lb, Rb).
+   * @brief Count "10" starts within [@p block_begin, @p block_end).
    * @details Accounts for cross-word boundaries.
    */
-  size_t rr_in_block(const size_t& Lb, const size_t& Rb) const noexcept {
-    if (Rb <= Lb + 1) {
+  size_t rr_in_block(const size_t& block_begin,
+                     const size_t& block_end) const noexcept {
+    if (block_end <= block_begin + 1) {
       return 0;
     }
-    size_t w_l = Lb >> 6;
-    const size_t w_r = (Rb - 1) >> 6;
-    const int off_l = Lb & 63;
-    const int off_r = (Rb - 1) & 63;
-    size_t cnt = 0;
+    size_t left_word_index = block_begin >> 6;
+    const size_t right_word_index = (block_end - 1) >> 6;
+    const int left_offset = block_begin & 63;
+    const int right_offset = (block_end - 1) & 63;
+    size_t count = 0;
 
-    if (w_l == w_r) {
-      const int len = off_r - off_l + 1;
-      const std::uint64_t slice = bits[w_l] >> off_l;
-      return pop10_in_slice64(slice, len);
+    if (left_word_index == right_word_index) {
+      const int length = right_offset - left_offset + 1;
+      const std::uint64_t slice = bits[left_word_index] >> left_offset;
+      return pop10_in_slice64(slice, length);
     }
 
     // prefix word
     {
-      const int len = 64 - off_l;
-      const std::uint64_t slice = bits[w_l] >> off_l;
-      cnt += pop10_in_slice64(slice, len);
+      const int length = 64 - left_offset;
+      const std::uint64_t slice = bits[left_word_index] >> left_offset;
+      count += pop10_in_slice64(slice, length);
     }
     // full interior words
-    for (size_t w = w_l + 1; w < w_r; ++w) {
-      const std::uint64_t x = bits[w];
-      cnt += pop10_in_slice64(x, 64);
+    for (size_t word_index = left_word_index + 1; word_index < right_word_index;
+         ++word_index) {
+      const std::uint64_t word = bits[word_index];
+      count += pop10_in_slice64(word, 64);
     }
     // suffix word
     {
-      const int len = off_r + 1;
-      const std::uint64_t mask =
-          (len == 64) ? ~std::uint64_t(0) : ((std::uint64_t(1) << len) - 1);
-      const std::uint64_t slice = bits[w_r] & mask;
-      cnt += pop10_in_slice64(slice, len);
+      const int length = right_offset + 1;
+      const std::uint64_t mask = (length == 64)
+                                     ? ~std::uint64_t(0)
+                                     : ((std::uint64_t(1) << length) - 1);
+      const std::uint64_t slice = bits[right_word_index] & mask;
+      count += pop10_in_slice64(slice, length);
     }
     // cross-word boundaries (bit 63 of w and bit 0 of w+1)
-    for (size_t w = w_l; w < w_r; ++w) {
-      if (((bits[w] >> 63) & 1u) && ((bits[w + 1] & 1u) == 0)) {
-        ++cnt;
+    for (size_t word_index = left_word_index; word_index < right_word_index;
+         ++word_index) {
+      if (((bits[word_index] >> 63) & 1u) &&
+          ((bits[word_index + 1] & 1u) == 0)) {
+        ++count;
       }
     }
-    return cnt;
+    return count;
   }
 
   /**
-   * @brief 1-based select of k-th "10" within [Lb, Rb).
+   * @brief 1-based select of @p target_pattern_rank-th "10" within [@p
+   * block_begin, @p block_end).
    * @return Position or npos.
    */
-  size_t select10_in_block(const size_t& Lb,
-                           const size_t& Rb,
-                           size_t k) const noexcept {
-    if (Rb <= Lb + 1) {
+  size_t select10_in_block(const size_t& block_begin,
+                           const size_t& block_end,
+                           size_t target_pattern_rank) const noexcept {
+    if (block_end <= block_begin + 1) {
       return npos;
     }
-    size_t w_l = Lb >> 6;
-    const size_t w_r = (Rb - 1) >> 6;
-    const int off_l = Lb & 63;
-    const int off_r = (Rb - 1) & 63;
+    size_t left_word_index = block_begin >> 6;
+    const size_t right_word_index = (block_end - 1) >> 6;
+    const int left_offset = block_begin & 63;
+    const int right_offset = (block_end - 1) & 63;
 
-    const auto select_in_masked_slice = [&](const std::uint64_t& slice,
-                                            const int& len,
-                                            const size_t& kk) noexcept -> int {
-      if (len <= 1) {
+    const auto select_in_masked_slice =
+        [&](const std::uint64_t& slice, const int& length,
+            const size_t& target_index) noexcept -> int {
+      if (length <= 1) {
         return -1;
       }
-      std::uint64_t P = slice & ~(slice >> 1);
-      if (len < 64) {
-        P &= ((std::uint64_t(1) << (len - 1)) - 1);
+      std::uint64_t pattern_mask = slice & ~(slice >> 1);
+      if (length < 64) {
+        pattern_mask &= ((std::uint64_t(1) << (length - 1)) - 1);
       } else {
-        P &= 0x7FFFFFFFFFFFFFFFull;
+        pattern_mask &= 0x7FFFFFFFFFFFFFFFull;
       }
-      return select_in_word(P, kk);
+      return select_in_word(pattern_mask, target_index);
     };
 
-    if (w_l == w_r) {
-      const int len = off_r - off_l + 1;
-      const std::uint64_t slice = bits[w_l] >> off_l;
-      const int off = select_in_masked_slice(slice, len, k);
-      return off >= 0 ? (Lb + (size_t)off) : npos;
+    if (left_word_index == right_word_index) {
+      const int length = right_offset - left_offset + 1;
+      const std::uint64_t slice = bits[left_word_index] >> left_offset;
+      const int offset =
+          select_in_masked_slice(slice, length, target_pattern_rank);
+      return offset >= 0 ? (block_begin + (size_t)offset) : npos;
     }
 
     // prefix word
     {
-      const int len = 64 - off_l;
-      const std::uint64_t slice = bits[w_l] >> off_l;
-      std::uint64_t P = slice & ~(slice >> 1);
-      P &= ((std::uint64_t(1) << (len - 1)) - 1);
-      const int c = std::popcount(P);
-      if (k <= (size_t)c) {
-        const int off = select_in_masked_slice(slice, len, k);
-        return Lb + (size_t)off;
+      const int length = 64 - left_offset;
+      const std::uint64_t slice = bits[left_word_index] >> left_offset;
+      std::uint64_t pattern_mask = slice & ~(slice >> 1);
+      pattern_mask &= ((std::uint64_t(1) << (length - 1)) - 1);
+      const int count = std::popcount(pattern_mask);
+      if (target_pattern_rank <= (size_t)count) {
+        const int offset =
+            select_in_masked_slice(slice, length, target_pattern_rank);
+        return block_begin + (size_t)offset;
       }
-      k -= c;
+      target_pattern_rank -= count;
     }
 
     // walk interior boundaries and words
-    for (size_t w = w_l; w + 1 < w_r; ++w) {
+    for (size_t word_index = left_word_index; word_index + 1 < right_word_index;
+         ++word_index) {
       // boundary between w and w+1
-      if (((bits[w] >> 63) & 1u) && ((bits[w + 1] & 1u) == 0)) {
-        if (--k == 0) {
-          return (w << 6) + 63;
+      if (((bits[word_index] >> 63) & 1u) &&
+          ((bits[word_index + 1] & 1u) == 0)) {
+        if (--target_pattern_rank == 0) {
+          return (word_index << 6) + 63;
         }
       }
       // full word w+1 (positions 0..62)
-      const std::uint64_t x = bits[w + 1];
-      const std::uint64_t P = (x & ~(x >> 1)) & 0x7FFFFFFFFFFFFFFFull;
-      const int c = std::popcount(P);
-      if (k <= (size_t)c) {
-        const int off = select_in_word(P, k);
-        if (off == -1) {
+      const std::uint64_t next_word = bits[word_index + 1];
+      const std::uint64_t pattern_mask =
+          (next_word & ~(next_word >> 1)) & 0x7FFFFFFFFFFFFFFFull;
+      const int count = std::popcount(pattern_mask);
+      if (target_pattern_rank <= (size_t)count) {
+        const int offset = select_in_word(pattern_mask, target_pattern_rank);
+        if (offset == -1) {
           return npos;
         }
-        return ((w + 1) << 6) + (size_t)off;
+        return ((word_index + 1) << 6) + (size_t)offset;
       }
-      k -= c;
+      target_pattern_rank -= count;
     }
 
     // boundary (w_r-1, w_r)
-    if (((bits[w_r - 1] >> 63) & 1u) && ((bits[w_r] & 1u) == 0)) {
-      if (--k == 0) {
-        return ((w_r - 1) << 6) + 63;
+    if (((bits[right_word_index - 1] >> 63) & 1u) &&
+        ((bits[right_word_index] & 1u) == 0)) {
+      if (--target_pattern_rank == 0) {
+        return ((right_word_index - 1) << 6) + 63;
       }
     }
 
     // suffix word w_r: [0..off_r]
     {
-      const int len = off_r + 1;
-      const std::uint64_t mask =
-          (len == 64) ? ~std::uint64_t(0) : ((std::uint64_t(1) << len) - 1);
-      const std::uint64_t slice = bits[w_r] & mask;
-      const int off = select_in_masked_slice(slice, len, k);
-      if (off >= 0) {
-        return (w_r << 6) + (size_t)off;
+      const int length = right_offset + 1;
+      const std::uint64_t mask = (length == 64)
+                                     ? ~std::uint64_t(0)
+                                     : ((std::uint64_t(1) << length) - 1);
+      const std::uint64_t slice = bits[right_word_index] & mask;
+      const int offset =
+          select_in_masked_slice(slice, length, target_pattern_rank);
+      if (offset >= 0) {
+        return (right_word_index << 6) + (size_t)offset;
       }
     }
     return npos;
@@ -1079,534 +1162,670 @@ class RmMTree {
     uint8_t pos_first_max;    // pos of first maximum in this byte
   };
 
+  struct LUT8Tables {
+    std::array<ByteAgg, 256> agg;
+    /**
+     * @brief For each byte and each delta ∈ [-8..+8] — the first position in
+     * the byte (0..7) where the prefix equals delta; -1 if none
+     */
+    std::array<std::array<int8_t, 17>, 256> fwd_pos;
+    /**
+     * @brief For each byte and each delta ∈ [-8..+8] — the last position in the
+     * byte (0..7) where the prefix equals delta; -1 if none
+     */
+    std::array<std::array<int8_t, 17>, 256> bwd_pos;
+  };
+
   /**
-   * @brief Returns the static lookup table of byte aggregates.
+   * @brief Unified initialization of all byte-based LUTs.
    */
-  static inline const std::array<ByteAgg, 256>& LUT8() noexcept {
-    static const std::array<ByteAgg, 256> T = [] {
-      std::array<ByteAgg, 256> t{};
-      for (int b = 0; b < 256; ++b) {
-        int cur = 0, mn = INT_MAX, mx = INT_MIN, cnt = 0, rrc = 0;
-        int pm = 0, pM = 0;
-        const auto get = [&](const int& k) {
-          return (b >> k) & 1;
+  static inline const LUT8Tables& LUT8_ALL() noexcept {
+    static const LUT8Tables tables = [] {
+      LUT8Tables lookup_tables{};
+      for (int byte_value = 0; byte_value < 256; ++byte_value) {
+        int current_excess = 0, min_prefix = INT_MAX, max_prefix = INT_MIN,
+            min_count = 0, pattern10_count = 0;
+        int first_min_position = 0, first_max_position = 0;
+        int prefixes[8];
+        const auto bit_at = [&](const int& bit_index) {
+          return (byte_value >> bit_index) & 1;
         };  // LSB-first
-        for (int k = 0; k < 8; ++k) {
-          int bit = get(k);
-          if (k + 1 < 8 && bit && get(k + 1) == 0) {
-            ++rrc;
+        for (int bit_index = 0; bit_index < 8; ++bit_index) {
+          int bit_value = bit_at(bit_index);
+          if (bit_index + 1 < 8 && bit_value && bit_at(bit_index + 1) == 0) {
+            ++pattern10_count;
+          }
+          current_excess += bit_value ? +1 : -1;
+          prefixes[bit_index] = current_excess;
+          if (current_excess < min_prefix) {
+            min_prefix = current_excess;
+            min_count = 1;
+            first_min_position = bit_index;
+          } else if (current_excess == min_prefix) {
+            ++min_count;
+          }
+          if (current_excess > max_prefix) {
+            max_prefix = current_excess;
+            first_max_position = bit_index;
           }
-          cur += bit ? +1 : -1;
-          if (cur < mn) {
-            mn = cur;
-            cnt = 1;
-            pm = k;
-          } else if (cur == mn) {
-            ++cnt;
+        }
+        ByteAgg aggregates{};
+        aggregates.excess_total = current_excess;
+        aggregates.min_prefix = (min_prefix == INT_MAX ? 0 : min_prefix);
+        aggregates.max_prefix = (max_prefix == INT_MIN ? 0 : max_prefix);
+        aggregates.min_count = min_count;
+        aggregates.pattern10_count = pattern10_count;
+        aggregates.first_bit = bit_at(0);
+        aggregates.last_bit = bit_at(7);
+        aggregates.pos_first_min = first_min_position;
+        aggregates.pos_first_max = first_max_position;
+        lookup_tables.agg[byte_value] = aggregates;
+        auto& forward_positions = lookup_tables.fwd_pos[byte_value];
+        auto& backward_positions = lookup_tables.bwd_pos[byte_value];
+        forward_positions.fill(-1);
+        backward_positions.fill(-1);
+        for (int delta = -8; delta <= 8; ++delta) {
+          for (int bit_index = 0; bit_index < 8; ++bit_index) {
+            if (prefixes[bit_index] == delta) {
+              forward_positions[delta + 8] = bit_index;
+              break;
+            }
           }
-          if (cur > mx) {
-            mx = cur;
-            pM = k;
+          for (int bit_index = 7; bit_index >= 0; --bit_index) {
+            if (prefixes[bit_index] == delta) {
+              backward_positions[delta + 8] = bit_index;
+              break;
+            }
           }
         }
-        ByteAgg a{};
-        a.excess_total = cur;
-        a.min_prefix = (mn == INT_MAX ? 0 : mn);
-        a.max_prefix = (mx == INT_MIN ? 0 : mx);
-        a.min_count = cnt;
-        a.pattern10_count = rrc;
-        a.first_bit = get(0);
-        a.last_bit = get(7);
-        a.pos_first_min = pm;
-        a.pos_first_max = pM;
-        t[b] = a;
-      }
-      return t;
+      }
+      return lookup_tables;
     }();
-    return T;
+    return tables;
   }
 
   /**
-   * @brief Extract 8 bits starting at position pos (LSB-first across words).
+   * @brief Byte-level aggregates.
    */
-  inline uint8_t get_byte(const size_t& pos) const noexcept {
-    const size_t w = pos >> 6;
-    const size_t off = pos & 63;
-    const std::uint64_t lo = bits[w] >> off;
-    if (off <= 56) {
-      return uint8_t(lo & 0xFFu);
-    }
-    const std::uint64_t hi = (w + 1 < bits.size()) ? bits[w + 1] : 0;
-    const std::uint64_t x = (lo | (hi << (64 - off))) & 0xFFu;
-    return uint8_t(x);
+  static inline const std::array<ByteAgg, 256>& LUT8() noexcept {
+    return LUT8_ALL().agg;
+  }
+
+  /**
+   * @brief Position of the first occurrence of delta within the byte.
+   */
+  static inline const std::array<std::array<int8_t, 17>, 256>&
+  LUT8_FWD_POS() noexcept {
+    return LUT8_ALL().fwd_pos;
+  }
+
+  /**
+   * @brief Position of the last occurrence of delta within the byte.
+   */
+  static inline const std::array<std::array<int8_t, 17>, 256>&
+  LUT8_BWD_POS() noexcept {
+    return LUT8_ALL().bwd_pos;
+  }
+
+  /**
+   * @brief Forward search within a single leaf over the range [@p search_start,
+   * @p search_end).
+   * @param required_delta required relative excess (relative to @p
+   * search_start).
+   * @return The position, or npos.
+   */
+  inline size_t scan_leaf_fwd(const size_t& search_start,
+                              const size_t& search_end,
+                              const int& required_delta) const noexcept {
+    if (search_start >= search_end) {
+      return npos;
+    }
+    const auto& aggregates_table = LUT8();
+    const auto& forward_lookup = LUT8_FWD_POS();
+    int current_excess = 0;
+    size_t position = search_start;
+    while (position + 8 <= search_end) {
+      const uint8_t byte_value = get_byte(position);
+      const auto& byte_aggregate = aggregates_table[byte_value];
+      const int local_need = required_delta - current_excess;
+      if (local_need >= byte_aggregate.min_prefix &&
+          local_need <= byte_aggregate.max_prefix && local_need >= -8 &&
+          local_need <= 8) {
+        const int8_t offset = forward_lookup[byte_value][local_need + 8];
+        if (offset >= 0) {
+          return position + size_t(offset);
+        }
+      }
+      current_excess += byte_aggregate.excess_total;
+      position += 8;
+    }
+
+    while (position < search_end) {
+      current_excess += bit(position) ? 1 : -1;
+      if (current_excess == required_delta) {
+        return position;
+      }
+      ++position;
+    }
+
+    return npos;
+  }
+
+  /**
+   * @brief Backward search within a single leaf over [@p block_begin, @p
+   * block_end) (does not look to the right of @p block_end).
+   * @param required_delta relative excess from @p block_begin.
+   * @param allow_right_boundary whether @p block_end may be used as a valid
+   * answer.
+   * @param global_right_border global right limit (as in descend_bwd).
+   */
+  inline size_t scan_leaf_bwd(
+      const size_t& block_begin,
+      const size_t& block_end,
+      const int& required_delta,
+      const bool& allow_right_boundary,
+      const size_t& global_right_border) const noexcept {
+    if (block_begin >= block_end) {
+      if ((block_begin < global_right_border || allow_right_boundary) &&
+          required_delta == 0) {
+        return block_begin;
+      }
+      return npos;
+    }
+
+    const int16_t* leaf_prefix_ptr =
+        &leaf_prefix[block_of(block_begin) * leaf_prefix_stride];
+    const size_t end_inclusive = block_end - block_begin;
+    if (end_inclusive > 0) {
+      const size_t position = ::find_backward_equal_i16_avx2(
+          leaf_prefix_ptr, 1, end_inclusive, required_delta, npos);
+      if (position != npos) {
+        return block_begin + position;
+      }
+    }
+
+    if ((block_begin < global_right_border || allow_right_boundary) &&
+        required_delta == 0) {
+      return block_begin;
+    }
+
+    return npos;
+  }
+
+  /**
+   * @brief Extract 8 bits starting at @p position (LSB-first across words).
+   */
+  inline uint8_t get_byte(const size_t& position) const noexcept {
+    const size_t word_index = position >> 6;
+    const size_t offset = position & 63;
+    const std::uint64_t lower_word = bits[word_index] >> offset;
+    if (offset <= 56) {
+      return uint8_t(lower_word & 0xFFu);
+    }
+    const std::uint64_t higher_word =
+        (word_index + 1 < bits.size()) ? bits[word_index + 1] : 0;
+    const std::uint64_t byte_value =
+        (lower_word | (higher_word << (64 - offset))) & 0xFFu;
+    return uint8_t(byte_value);
   }
 
   /**
    * @brief Descend to the first (leftmost) maximum with node-relative prefix
-   * equal to d.
-   * @param v Node index.
-   * @param d Target prefix within node coordinates.
-   * @param base Starting global position of node.
+   * equal to @p target_prefix.
+   * @param node_index Node index.
+   * @param target_prefix Target prefix within node coordinates.
+   * @param segment_base Starting global position of node.
    * @return Position or npos.
    */
-  size_t descend_first_max(size_t v, int d, size_t base) const noexcept {
-    const size_t leaf0 = first_leaf_index();
-    while (v < leaf0) {
-      const size_t Lc = v << 1, Rc = Lc | 1;
-      const int leftX = node_max_prefix_excess[Lc];
-      const int rightX = node_total_excess[Lc] + node_max_prefix_excess[Rc];
-      if (leftX >= rightX && leftX == d) {
-        v = Lc;
-      } else if (rightX == d) {
-        base += segment_size_bits[Lc];
-        d -= node_total_excess[Lc];
-        v = Rc;
+  size_t descend_first_max(size_t node_index,
+                           int target_prefix,
+                           size_t segment_base) const noexcept {
+    while (node_index < first_leaf_index) {
+      const size_t left_child = node_index << 1, right_child = left_child | 1;
+      const int left_max = node_max_prefix_excess[left_child];
+      const int right_max =
+          node_total_excess[left_child] + node_max_prefix_excess[right_child];
+      if (left_max >= right_max && left_max == target_prefix) {
+        node_index = left_child;
+      } else if (right_max == target_prefix) {
+        segment_base += segment_size_bits[left_child];
+        target_prefix -= node_total_excess[left_child];
+        node_index = right_child;
       } else {
         return npos;
       }
     }
 
-    const size_t Lb = base;
-    const size_t Rb = std::min(base + segment_size_bits[v], num_bits);
-    int mx;
-    size_t pos;
+    const size_t segment_begin = segment_base;
+    const size_t segment_end =
+        std::min(segment_base + segment_size_bits[node_index], num_bits);
+    int max_value;
+    size_t position;
 
-    first_max_value_pos8(Lb, Rb ? (Rb - 1) : Lb, mx, pos);
-    return (mx == d ? pos : npos);
+    first_max_value_pos8(segment_begin,
+                         segment_end ? (segment_end - 1) : segment_begin,
+                         max_value, position);
+    return (max_value == target_prefix ? position : npos);
   }
 
   /**
    * @brief Heap index of the first leaf node.
    */
-  inline size_t first_leaf_index() const noexcept {
-    return std::bit_ceil(std::max<size_t>(1, leaf_count));
-  }
+  size_t first_leaf_index = 1;
 
   /**
-   * @brief Block index containing position i.
+   * @brief Block index containing position @p position.
    */
-  size_t block_of(const size_t& i) const noexcept { return i / block_bits; }
+  size_t block_of(const size_t& position) const noexcept {
+    return position / block_bits;
+  }
 
   /**
-   * @brief Heap index of the leaf whose segment starts at block_start.
+   * @brief Heap index of the leaf whose segment starts at @p block_start.
    */
   size_t leaf_index_of(const size_t& block_start) const noexcept {
-    return first_leaf_index() + block_of(block_start);
+    return first_leaf_index + block_of(block_start);
   }
 
   /**
-   * @brief Starting global position for node v (0-indexed).
+   * @brief Starting global position for node @p node_index (0-indexed).
    * @details Walk up to compute base for internal nodes; direct for leaves.
    */
-  size_t node_base(size_t v) const noexcept {
-    const size_t leaf0 = first_leaf_index();
-    if (v >= leaf0) {
-      return (v - leaf0) * block_bits;
+  size_t node_base(size_t node_index) const noexcept {
+    if (node_index >= first_leaf_index) {
+      return (node_index - first_leaf_index) * block_bits;
     }
 
     size_t base = 0;
-    for (; v > 1; v >>= 1) {
-      if (v & 1) {
-        base += segment_size_bits[v - 1];
+    for (; node_index > 1; node_index >>= 1) {
+      if (node_index & 1) {
+        base += segment_size_bits[node_index - 1];
       }
     }
     return base;
   }
 
   /**
-   * @brief Cover a range of whole blocks [a..b] (inclusive) with O(log n)
-   * maximal nodes.
+   * @brief Cover a range of whole blocks [@p block_begin_index..@p
+   * block_end_index] (inclusive) with O(log n) maximal nodes.
    * @details Returns node indices in left-to-right order.
    */
-  std::vector<size_t> cover_blocks(const size_t& a, const size_t& b) const {
-    const size_t leaf0 = first_leaf_index();
-    size_t l = leaf0 + a;
-    size_t r = leaf0 + b;
-    std::vector<size_t> Lnodes, Rnodes;
-    while (l <= r) {
-      if ((l & 1) == 1) {
-        Lnodes.push_back(l++);
-      }
-      if ((r & 1) == 0) {
-        Rnodes.push_back(r--);
-      }
-      l >>= 1;
-      r >>= 1;
-    }
-    std::reverse(Rnodes.begin(), Rnodes.end());
-    Lnodes.insert(Lnodes.end(), Rnodes.begin(), Rnodes.end());
-    return Lnodes;
+  std::vector<size_t> cover_blocks(const size_t& block_begin_index,
+                                   const size_t& block_end_index) const {
+    size_t left_index = first_leaf_index + block_begin_index;
+    size_t right_index = first_leaf_index + block_end_index;
+    std::vector<size_t> left_nodes, right_nodes;
+    while (left_index <= right_index) {
+      if ((left_index & 1) == 1) {
+        left_nodes.push_back(left_index++);
+      }
+      if ((right_index & 1) == 0) {
+        right_nodes.push_back(right_index--);
+      }
+      left_index >>= 1;
+      right_index >>= 1;
+    }
+    std::reverse(right_nodes.begin(), right_nodes.end());
+    left_nodes.insert(left_nodes.end(), right_nodes.begin(), right_nodes.end());
+    return left_nodes;
   }
 
   /**
    * @brief Descend for fwdsearch to find first position where relative prefix
-   * equals 'need'.
+   * equals @p need.
    */
-  size_t descend_fwd(size_t v, int need, size_t base) const noexcept {
-    const size_t leaf0 = first_leaf_index();
-    while (v < leaf0) {
-      const size_t Lc = v << 1;
-      const size_t Rc = Lc | 1;
-      if (node_min_prefix_excess[Lc] <= need &&
-          need <= node_max_prefix_excess[Lc]) {
-        v = Lc;
+  size_t descend_fwd(size_t node_index,
+                     int required_delta,
+                     size_t segment_base) const noexcept {
+    while (node_index < first_leaf_index) {
+      const size_t left_child = node_index << 1;
+      const size_t right_child = left_child | 1;
+      if (node_min_prefix_excess[left_child] <= required_delta &&
+          required_delta <= node_max_prefix_excess[left_child]) {
+        node_index = left_child;
       } else {
-        need -= node_total_excess[Lc];
-        base += segment_size_bits[Lc];
-        v = Rc;
+        required_delta -= node_total_excess[left_child];
+        segment_base += segment_size_bits[left_child];
+        node_index = right_child;
       }
     }
-    const size_t Rb = std::min(base + segment_size_bits[v], num_bits);
-    int cur = 0;
-    for (size_t p = base; p < Rb; ++p) {
-      cur += bit(p) ? +1 : -1;
-      if (cur == need) {
-        return p;
-      }
-    }
-    return npos;
+    return scan_leaf_fwd(
+        segment_base,
+        std::min(segment_base + segment_size_bits[node_index], num_bits),
+        required_delta);
   }
 
   /**
    * @brief Descend for bwdsearch to return the rightmost solution.
-   * @param v Current node.
-   * @param base Left border of node v.
-   * @param need Target relative prefix inside v.
-   * @param right_border Global right limit (exclusive if !allow_rb).
-   * @param allow_rb Whether right border is allowed to match.
+   * @param node_index Current node.
+   * @param segment_base Left border of node @p node_index.
+   * @param required_delta Target relative prefix inside @p node_index.
+   * @param global_right_border Global right limit (exclusive if
+   * !@p allow_right_boundary).
+   * @param allow_right_boundary Whether right border is allowed to match.
    */
-  size_t descend_bwd(size_t v,
-                     const size_t& base,
-                     const int& need,
-                     const size_t& right_border,
-                     const bool& allow_rb) const noexcept {
-    const size_t leaf0 = first_leaf_index();
-    while (v < leaf0) {
-      const size_t Lc = v << 1;
-      const size_t Rc = Lc | 1;
-      const int need_r = need - node_total_excess[Lc];
+  size_t descend_bwd(size_t node_index,
+                     const size_t& segment_base,
+                     const int& required_delta,
+                     const size_t& global_right_border,
+                     const bool& allow_right_boundary) const noexcept {
+    while (node_index < first_leaf_index) {
+      const size_t left_child = node_index << 1;
+      const size_t right_child = left_child | 1;
+      const int required_in_right =
+          required_delta - node_total_excess[left_child];
 
       // 1) try the right child first (to capture the rightmost j)
-      if (node_min_prefix_excess[Rc] <= need_r &&
-          need_r <= node_max_prefix_excess[Rc]) {
-        const size_t ans = descend_bwd(Rc, base + segment_size_bits[Lc], need_r,
-                                       right_border, allow_rb);
-        if (ans != npos) {
-          return ans;
+      if (node_min_prefix_excess[right_child] <= required_in_right &&
+          required_in_right <= node_max_prefix_excess[right_child]) {
+        const size_t result = descend_bwd(
+            right_child, segment_base + segment_size_bits[left_child],
+            required_in_right, global_right_border, allow_right_boundary);
+        if (result != npos) {
+          return result;
         }
       }
 
       // 2) junction between children (end of the left child)
-      const size_t j_border = base + segment_size_bits[Lc];
-      if (need == node_total_excess[Lc] &&
-          (j_border < right_border || allow_rb)) {
-        return j_border;
+      const size_t junction = segment_base + segment_size_bits[left_child];
+      if (required_delta == node_total_excess[left_child] &&
+          (junction < global_right_border || allow_right_boundary)) {
+        return junction;
       }
 
       // 3) can we move left within the range?
-      if (node_min_prefix_excess[Lc] <= need &&
-          need <= node_max_prefix_excess[Lc]) {
-        v = Lc;
+      if (node_min_prefix_excess[left_child] <= required_delta &&
+          required_delta <= node_max_prefix_excess[left_child]) {
+        node_index = left_child;
         continue;
       }
 
       // None of (1)-(3) worked. The only possible point is the left border of
       // the node.
-      if (need == 0 && (base < right_border || allow_rb)) {
-        return base;
+      if (required_delta == 0 &&
+          (segment_base < global_right_border || allow_right_boundary)) {
+        return segment_base;
       }
 
       return npos;
     }
 
-    const size_t Lb = base;
-    const size_t Rb = std::min(base + segment_size_bits[v], num_bits);
-    const size_t RB = std::min(right_border, Rb);
-
-    int cur = 0;
-    for (size_t p = Lb; p < RB; ++p) {
-      cur += bit(p) ? +1 : -1;
-    }
-
-    if (allow_rb && cur == need) {
-      return RB;
-    }
-
-    for (size_t p = RB; p > Lb;) {
-      --p;
-      cur += bit(p) ? -1 : +1;
-      if (cur == need) {
-        return p;
-      }
-      if (p == 0) {
-        break;
-      }
-    }
-
-    if ((Lb < right_border || allow_rb) && cur == need) {
-      return Lb;
-    }
-
-    return npos;
+    return scan_leaf_bwd(
+        segment_base,
+        std::min(
+            global_right_border,
+            std::min(segment_base + segment_size_bits[node_index], num_bits)),
+        required_delta, allow_right_boundary, global_right_border);
   }
 
   /**
-   * @brief Descend to find first position where node-relative prefix equals d
-   * (minimum).
+   * @brief Descend to find first position where node-relative prefix equals
+   * @p target_prefix (minimum).
    */
-  size_t descend_first_min(size_t v, int d, size_t base) const noexcept {
-    const size_t leaf0 = first_leaf_index();
-    while (v < leaf0) {
-      const size_t Lc = v << 1, Rc = Lc | 1;
-      const int leftm = node_min_prefix_excess[Lc];
-      const int rightm = node_total_excess[Lc] + node_min_prefix_excess[Rc];
-      if (leftm <= rightm && leftm == d) {
-        v = Lc;
-      } else if (rightm == d) {
-        base += segment_size_bits[Lc];
-        d -= node_total_excess[Lc];
-        v = Rc;
+  size_t descend_first_min(size_t node_index,
+                           int target_prefix,
+                           size_t segment_base) const noexcept {
+    while (node_index < first_leaf_index) {
+      const size_t left_child = node_index << 1, right_child = left_child | 1;
+      const int left_min = node_min_prefix_excess[left_child];
+      const int right_min =
+          node_total_excess[left_child] + node_min_prefix_excess[right_child];
+      if (left_min <= right_min && left_min == target_prefix) {
+        node_index = left_child;
+      } else if (right_min == target_prefix) {
+        segment_base += segment_size_bits[left_child];
+        target_prefix -= node_total_excess[left_child];
+        node_index = right_child;
       } else {
         return npos;
       }
     }
 
-    const size_t Lb = base;
-    const size_t Rb = std::min(base + segment_size_bits[v], num_bits);
-    int mn;
-    size_t pos;
+    const size_t segment_begin = segment_base;
+    const size_t segment_end =
+        std::min(segment_base + segment_size_bits[node_index], num_bits);
+    int min_value;
+    size_t position;
 
-    first_min_value_pos8(Lb, Rb ? (Rb - 1) : Lb, mn, pos);
-    return (mn == d ? pos : npos);
+    first_min_value_pos8(segment_begin,
+                         segment_end ? (segment_end - 1) : segment_begin,
+                         min_value, position);
+    return (min_value == target_prefix ? position : npos);
   }
 
   /**
-   * @brief Descend to find the q-th minimum (1-based) where node-relative
-   * prefix equals d.
+   * @brief Descend to find the @p target_min_rank-th minimum (1-based) where
+   * node-relative prefix equals @p target_prefix.
    */
-  size_t descend_qth_min(size_t v,
-                         int d,
-                         size_t q,
-                         size_t base) const noexcept {
-    const size_t leaf0 = first_leaf_index();
-    while (v < leaf0) {
-      const size_t Lc = v << 1;
-      const size_t Rc = Lc | 1;
-      const int leftm = node_min_prefix_excess[Lc];
-      const int rightm = node_total_excess[Lc] + node_min_prefix_excess[Rc];
-      if (leftm == d) {
-        if (node_min_count[Lc] >= q) {
-          v = Lc;
+  size_t descend_qth_min(size_t node_index,
+                         int target_prefix,
+                         size_t target_min_rank,
+                         size_t segment_base) const noexcept {
+    while (node_index < first_leaf_index) {
+      const size_t left_child = node_index << 1;
+      const size_t right_child = left_child | 1;
+      const int left_min = node_min_prefix_excess[left_child];
+      const int right_min =
+          node_total_excess[left_child] + node_min_prefix_excess[right_child];
+      if (left_min == target_prefix) {
+        if (node_min_count[left_child] >= target_min_rank) {
+          node_index = left_child;
           continue;
         }
-        q -= node_min_count[Lc];
+        target_min_rank -= node_min_count[left_child];
       }
-      if (rightm == d) {
-        base += segment_size_bits[Lc];
-        d -= node_total_excess[Lc];
-        v = Rc;
+      if (right_min == target_prefix) {
+        segment_base += segment_size_bits[left_child];
+        target_prefix -= node_total_excess[left_child];
+        node_index = right_child;
         continue;
       }
       return npos;
     }
     return qth_min_in_block(
-        base, std::min(base + segment_size_bits[v], num_bits) - 1, q);
+        segment_base,
+        std::min(segment_base + segment_size_bits[node_index], num_bits) - 1,
+        target_min_rank);
   }
 
   /**
-   * @brief 1-based select of k-th '1' within [Lb, Rb).
+   * @brief 1-based select of @p target_one_rank-th '1' within [@p block_begin,
+   * @p block_end).
    */
-  size_t select1_in_block(const size_t& Lb,
-                          const size_t& Rb,
-                          size_t k) const noexcept {
-    size_t w_l = Lb >> 6;
-    const size_t w_r = (Rb >> 6);
-    const size_t off_l = Lb & 63;
-    const std::uint64_t mask_l =
-        (off_l ? (~std::uint64_t(0) << off_l) : ~std::uint64_t(0));
-    if (w_l == w_r) {
-      const std::uint64_t w = bits[w_l] & mask_l &
-                              ((Rb & 63) ? ((std::uint64_t(1) << (Rb & 63)) - 1)
-                                         : ~std::uint64_t(0));
-      return Lb + select_in_word(w, k);
+  size_t select1_in_block(const size_t& block_begin,
+                          const size_t& block_end,
+                          size_t target_one_rank) const noexcept {
+    size_t left_word_index = block_begin >> 6;
+    const size_t right_word_index = (block_end >> 6);
+    const size_t left_offset = block_begin & 63;
+    const std::uint64_t left_mask =
+        (left_offset ? (~std::uint64_t(0) << left_offset) : ~std::uint64_t(0));
+    if (left_word_index == right_word_index) {
+      const std::uint64_t word =
+          bits[left_word_index] & left_mask &
+          ((block_end & 63) ? ((std::uint64_t(1) << (block_end & 63)) - 1)
+                            : ~std::uint64_t(0));
+      return block_begin + select_in_word(word, target_one_rank);
     }
     // prefix
-    if (off_l) {
-      const std::uint64_t w = bits[w_l] & mask_l;
-      const int c = std::popcount(w);
-      if (k <= (size_t)c) {
-        return Lb + select_in_word(w, k);
+    if (left_offset) {
+      const std::uint64_t word = bits[left_word_index] & left_mask;
+      const int count = std::popcount(word);
+      if (target_one_rank <= (size_t)count) {
+        return block_begin + select_in_word(word, target_one_rank);
       }
-      k -= c;
-      w_l++;
+      target_one_rank -= count;
+      left_word_index++;
     }
     // full words
-    while (w_l < w_r) {
-      const std::uint64_t w = bits[w_l];
-      const int c = std::popcount(w);
-      if (k <= (size_t)c) {
-        return (w_l << 6) + select_in_word(w, k);
+    while (left_word_index < right_word_index) {
+      const std::uint64_t word = bits[left_word_index];
+      const int count = std::popcount(word);
+      if (target_one_rank <= (size_t)count) {
+        return (left_word_index << 6) + select_in_word(word, target_one_rank);
       }
-      k -= c;
-      ++w_l;
+      target_one_rank -= count;
+      ++left_word_index;
     }
     // tail
-    const size_t off_r = Rb & 63;
-    if (off_r) {
-      const std::uint64_t w = bits[w_l] & ((std::uint64_t(1) << off_r) - 1);
-      const int c = std::popcount(w);
-      if (k <= (size_t)c) {
-        return (w_l << 6) + select_in_word(w, k);
+    const size_t right_offset = block_end & 63;
+    if (right_offset) {
+      const std::uint64_t word =
+          bits[left_word_index] & ((std::uint64_t(1) << right_offset) - 1);
+      const int count = std::popcount(word);
+      if (target_one_rank <= (size_t)count) {
+        return (left_word_index << 6) + select_in_word(word, target_one_rank);
       }
     }
     return npos;
   }
 
   /**
-   * @brief 1-based select of k-th '0' within [Lb, Rb).
+   * @brief 1-based select of @p target_zero_rank-th '0' within [@p
+   * block_begin, @p block_end).
    */
-  size_t select0_in_block(const size_t& Lb,
-                          const size_t& Rb,
-                          size_t k) const noexcept {
-    if (Rb <= Lb) {
+  size_t select0_in_block(const size_t& block_begin,
+                          const size_t& block_end,
+                          size_t target_zero_rank) const noexcept {
+    if (block_end <= block_begin) {
       return npos;
     }
 
-    size_t w_l = Lb >> 6;
-    const size_t w_r = Rb >> 6;
-    const size_t off_l = Lb & 63;
+    size_t left_word_index = block_begin >> 6;
+    const size_t right_word_index = block_end >> 6;
+    const size_t left_offset = block_begin & 63;
 
-    if (w_l == w_r) {
-      const std::uint64_t mask_l =
-          (off_l ? (~std::uint64_t(0) << off_l) : ~std::uint64_t(0));
-      const std::uint64_t mask_r =
-          ((Rb & 63) ? ((std::uint64_t(1) << (Rb & 63)) - 1)
-                     : ~std::uint64_t(0));
-      const std::uint64_t w = (~bits[w_l]) & mask_l & mask_r;
-      const int off = select_in_word(w, k);
-      return (off >= 0) ? (Lb + (size_t)off) : npos;
+    if (left_word_index == right_word_index) {
+      const std::uint64_t left_mask =
+          (left_offset ? (~std::uint64_t(0) << left_offset)
+                       : ~std::uint64_t(0));
+      const std::uint64_t right_mask =
+          ((block_end & 63) ? ((std::uint64_t(1) << (block_end & 63)) - 1)
+                            : ~std::uint64_t(0));
+      const std::uint64_t word =
+          (~bits[left_word_index]) & left_mask & right_mask;
+      const int offset = select_in_word(word, target_zero_rank);
+      return (offset >= 0) ? (block_begin + (size_t)offset) : npos;
     }
 
     // prefix
-    if (off_l) {
-      const std::uint64_t w = (~bits[w_l]) & (~std::uint64_t(0) << off_l);
-      const int c = std::popcount(w);
-      if (k <= (size_t)c) {
-        const int off = select_in_word(w, k);
-        return (off >= 0) ? (Lb + (size_t)off) : npos;
+    if (left_offset) {
+      const std::uint64_t word =
+          (~bits[left_word_index]) & (~std::uint64_t(0) << left_offset);
+      const int count = std::popcount(word);
+      if (target_zero_rank <= (size_t)count) {
+        const int offset = select_in_word(word, target_zero_rank);
+        return (offset >= 0) ? (block_begin + (size_t)offset) : npos;
       }
-      k -= c;
-      ++w_l;
+      target_zero_rank -= count;
+      ++left_word_index;
     }
 
     // full words
-    while (w_l < w_r) {
-      const std::uint64_t w = ~bits[w_l];
-      const int c = std::popcount(w);
-      if (k <= (size_t)c) {
-        const int off = select_in_word(w, k);
-        return (off >= 0) ? ((w_l << 6) + (size_t)off) : npos;
+    while (left_word_index < right_word_index) {
+      const std::uint64_t word = ~bits[left_word_index];
+      const int count = std::popcount(word);
+      if (target_zero_rank <= (size_t)count) {
+        const int offset = select_in_word(word, target_zero_rank);
+        return (offset >= 0) ? ((left_word_index << 6) + (size_t)offset) : npos;
       }
-      k -= c;
-      ++w_l;
+      target_zero_rank -= count;
+      ++left_word_index;
     }
 
     // tail
-    const size_t off_r = Rb & 63;
-    if (off_r) {
-      const std::uint64_t w = (~bits[w_l]) & ((std::uint64_t(1) << off_r) - 1);
-      const int c = std::popcount(w);
-      if (k <= (size_t)c) {
-        const int off = select_in_word(w, k);
-        return (off >= 0) ? ((w_l << 6) + (size_t)off) : npos;
+    const size_t right_offset = block_end & 63;
+    if (right_offset) {
+      const std::uint64_t word =
+          (~bits[left_word_index]) & ((std::uint64_t(1) << right_offset) - 1);
+      const int count = std::popcount(word);
+      if (target_zero_rank <= (size_t)count) {
+        const int offset = select_in_word(word, target_zero_rank);
+        return (offset >= 0) ? ((left_word_index << 6) + (size_t)offset) : npos;
       }
     }
     return npos;
   }
 
   /**
-   * @brief 1-based select of k-th set bit inside a 64-bit word.
+   * @brief 1-based select of @p target_rank-th set bit inside a 64-bit word.
    * @return Bit index [0..63] or −1 if not found.
    */
-  static inline int select_in_word(std::uint64_t w, size_t k) noexcept {
-#ifdef __GNUC__
-    while (w) {
-      if (--k == 0) {
-        return __builtin_ctzll(w);
-      }
-      w &= (w - 1);
-    }
-    return -1;
-#else
-    for (int i = 0; i < 64; ++i) {
-      if ((w >> i) & 1u) {
-        if (--k == 0) {
-          return i;
-        }
+  static inline int select_in_word(std::uint64_t word,
+                                   size_t target_rank) noexcept {
+    while (word) {
+      if (--target_rank == 0) {
+        return std::countr_zero(word);
       }
+      word &= (word - 1);
     }
     return -1;
-#endif
   }
 
   /**
    * @brief Ceil division of positive integers.
    */
-  static inline size_t ceil_div(const size_t& a, const size_t& b) noexcept {
-    return (a + b - 1) / b;
+  static inline size_t ceil_div(const size_t& numerator,
+                                const size_t& denominator) noexcept {
+    return (numerator + denominator - 1) / denominator;
   }
 
   /**
-   * @brief Number of node slots needed for @p Nbits with leaf size @p Bpow2.
+   * @brief Number of node slots needed for @p bit_count with leaf size @p
+   * block_size_pow2.
    * @details Includes internal node space (heap layout) plus leaves.
    */
-  static inline size_t nodeslots_for(const size_t& Nbits,
-                                     const size_t& Bpow2) noexcept {
-    if (Nbits == 0) {
+  static inline size_t nodeslots_for(const size_t& bit_count,
+                                     const size_t& block_size_pow2) noexcept {
+    if (bit_count == 0) {
       return 0;
     }
-    size_t leaf_count = ceil_div(Nbits, Bpow2);
-    return std::bit_ceil(std::max<size_t>(1, leaf_count)) + leaf_count;
+    size_t leaf_node_count = ceil_div(bit_count, block_size_pow2);
+    return std::bit_ceil(std::max<size_t>(1, leaf_node_count)) +
+           leaf_node_count;
   }
 
   /**
    * @brief Auxiliary overhead in bytes/bitvector bytes for given parameters.
    */
-  static inline float overhead_for(const size_t& Nbits,
-                                   const size_t& Bpow2) noexcept {
+  static inline float overhead_for(const size_t& bit_count,
+                                   const size_t& block_size_pow2) noexcept {
     static constexpr size_t AUX_SLOT_BYTES =
         sizeof(uint32_t) + sizeof(int32_t) + sizeof(int32_t) + sizeof(int32_t) +
-        sizeof(uint32_t) + sizeof(uint32_t) + sizeof(uint8_t) + sizeof(uint8_t);
+        sizeof(uint32_t) + sizeof(uint32_t) + sizeof(uint8_t) +
+        sizeof(uint8_t) + sizeof(int16_t);
 
-    size_t bb = ceil_div(Nbits, 64) * 8;
-    if (bb == 0) {
+    size_t bitvector_bytes = ceil_div(bit_count, 64) * 8;
+    if (bitvector_bytes == 0) {
       return 0;
     }
-    size_t slots = nodeslots_for(Nbits, Bpow2);
-    size_t aux = slots * AUX_SLOT_BYTES;
-    return ((float)aux) / ((float)bb);
+    size_t slot_count = nodeslots_for(bit_count, block_size_pow2);
+    size_t aux_bytes = slot_count * AUX_SLOT_BYTES;
+    return ((float)aux_bytes) / ((float)bitvector_bytes);
   }
 
   /**
-   * @brief Choose minimal block size (power of two) keeping overhead ≤ cap.
-   * @details Returns 64 if cap < 0 (no constraint). Clamped to ≤16384 or Nbits.
+   * @brief Choose minimal block size (power of two) keeping overhead ≤ @p
+   * overhead_cap.
+   * @details Returns 64 if @p overhead_cap < 0 (no constraint). Clamped to
+   * ≤16384 or @p bit_count.
    */
   static inline size_t choose_block_bits_for_overhead(
-      const size_t& Nbits,
-      const float& cap) noexcept {
-    if (cap < 0.f) {
+      const size_t& bit_count,
+      const float& overhead_cap) noexcept {
+    if (overhead_cap < 0.f) {
       return 64;
     }
 
-    const size_t Bmax = std::min<size_t>(Nbits, 16384);
-    size_t block_bits = 64;
-    while (block_bits < Bmax) {
-      if (overhead_for(Nbits, block_bits) <= cap) {
+    const size_t max_block_bits = std::min<size_t>(bit_count, 16384);
+    size_t candidate_block_bits = 64;
+    while (candidate_block_bits < max_block_bits) {
+      if (overhead_for(bit_count, candidate_block_bits) <= overhead_cap) {
         break;
       }
-      block_bits <<= 1;
+      candidate_block_bits <<= 1;
     }
-    return block_bits;
+    return candidate_block_bits;
   }
 
   /**
@@ -1614,13 +1833,13 @@ class RmMTree {
    * @param leaf_block_bits Desired leaf size (0 = auto).
    * @param max_overhead Overhead cap (<0 = no cap).
    */
-  void build_from_string(const std::string& s,
+  void build_from_string(const std::string& bit_string_input,
                          const size_t& leaf_block_bits = 0,
                          const float& max_overhead = -1.0) {
-    num_bits = s.size();
+    num_bits = bit_string_input.size();
     bits.assign((num_bits + 63) / 64, 0);
     for (size_t i = 0; i < num_bits; ++i) {
-      if (s[i] == '1') {
+      if (bit_string_input[i] == '1') {
         set1(i);
       }
     }
@@ -1630,16 +1849,16 @@ class RmMTree {
   /**
    * @brief Build internal structures from 64-bit words.
    * @param words Words with LSB-first bits.
-   * @param Nbits Number of valid bits.
+   * @param bit_count Number of valid bits.
    * @param leaf_block_bits Desired leaf size (0 = auto).
    * @param max_overhead Overhead cap (<0 = no cap).
    */
   void build_from_words(const std::vector<std::uint64_t>& words,
-                        const size_t& Nbits,
+                        const size_t& bit_count,
                         const size_t& leaf_block_bits = 0,
                         const float& max_overhead = -1.0) {
     bits = words;
-    num_bits = Nbits;
+    num_bits = bit_count;
     if (bits.size() * 64 < num_bits) {
       bits.resize((num_bits + 63) / 64);
     }
@@ -1647,265 +1866,375 @@ class RmMTree {
   }
 
   /**
-   * @brief Read bit at position i (LSB-first across words).
+   * @brief Read bit at position @p position (LSB-first across words).
    */
-  inline int bit(const size_t& i) const noexcept {
-    return (bits[i >> 6] >> (i & 63)) & 1u;
+  inline int bit(const size_t& position) const noexcept {
+    return (bits[position >> 6] >> (position & 63)) & 1u;
   }
 
   /**
-   * @brief Set bit at position i to 1.
+   * @brief Set bit at position @p position to 1.
    */
-  inline void set1(const size_t& i) noexcept {
-    bits[i >> 6] |= (std::uint64_t(1) << (i & 63));
+  inline void set1(const size_t& position) noexcept {
+    bits[position >> 6] |= (std::uint64_t(1) << (position & 63));
   }
 
   /**
-   * @brief Number of ones in node v computed from size and total excess.
+   * @brief Number of ones in node @p node_index computed from size and total
+   * excess.
    */
-  inline uint32_t ones_in_node(const size_t& v) const noexcept {
-    return ((int64_t)segment_size_bits[v] + (int64_t)node_total_excess[v]) >> 1;
+  inline uint32_t ones_in_node(const size_t& node_index) const noexcept {
+    return ((int64_t)segment_size_bits[node_index] +
+            (int64_t)node_total_excess[node_index]) >>
+           1;
   }
 
   /**
-   * @brief Scan [l..r] inclusive, computing minimum value and how many times it
-   * is attained.
-   * @details Uses 8-bit LUT for speed; outputs cur at end, mn, and count.
+   * @brief Scan [@p range_begin..@p range_end] inclusive, computing minimum
+   * value and how many
+   * times it is attained.
+   * @details Uses 8-bit LUT for speed; outputs @p current_excess at end, @p
+   * min_value, and @p count.
    */
-  inline void scan_range_min_count8(size_t l,
-                                    const size_t& r,
-                                    int& cur,
-                                    int& mn,
-                                    uint32_t& cnt) const noexcept {
-    cur = 0;
-    mn = INT_MAX;
-    cnt = 0;
-    if (r < l) {
-      mn = 0;
+  inline void scan_range_min_count8(size_t range_begin,
+                                    const size_t& range_end,
+                                    int& current_excess,
+                                    int& min_value,
+                                    uint32_t& count) const noexcept {
+    current_excess = 0;
+    min_value = INT_MAX;
+    count = 0;
+    if (range_end < range_begin) {
+      min_value = 0;
       return;
     }
     // to byte alignment
-    while (l <= r && (l & 7)) {
-      cur += bit(l) ? +1 : -1;
-      if (cur < mn) {
-        mn = cur;
-        cnt = 1;
-      } else if (cur == mn) {
-        ++cnt;
+    while (range_begin <= range_end && (range_begin & 7)) {
+      current_excess += bit(range_begin) ? +1 : -1;
+      if (current_excess < min_value) {
+        min_value = current_excess;
+        count = 1;
+      } else if (current_excess == min_value) {
+        ++count;
       }
-      ++l;
+      ++range_begin;
     }
     // full bytes
-    const auto& T = LUT8();
-    while (l + 7 <= r) {
-      const auto& a = T[get_byte(l)];
-      const int cand = cur + a.min_prefix;
-      if (cand < mn) {
-        mn = cand;
-        cnt = a.min_count;
-      } else if (cand == mn) {
-        cnt += a.min_count;
-      }
-      cur += a.excess_total;
-      l += 8;
+    const auto& aggregates_table = LUT8();
+    while (range_begin + 7 <= range_end) {
+      const auto& byte_aggregate = aggregates_table[get_byte(range_begin)];
+      const int candidate = current_excess + byte_aggregate.min_prefix;
+      if (candidate < min_value) {
+        min_value = candidate;
+        count = byte_aggregate.min_count;
+      } else if (candidate == min_value) {
+        count += byte_aggregate.min_count;
+      }
+      current_excess += byte_aggregate.excess_total;
+      range_begin += 8;
     }
     // tail
-    while (l <= r) {
-      cur += bit(l) ? +1 : -1;
-      if (cur < mn) {
-        mn = cur;
-        cnt = 1;
-      } else if (cur == mn) {
-        ++cnt;
+    while (range_begin <= range_end) {
+      current_excess += bit(range_begin) ? +1 : -1;
+      if (current_excess < min_value) {
+        min_value = current_excess;
+        count = 1;
+      } else if (current_excess == min_value) {
+        ++count;
+      }
+      ++range_begin;
+    }
+    if (min_value == INT_MAX) {
+      min_value = count = 0;
+    }
+  }
+
+  /**
+   * @brief Cover the range of whole blocks [@p block_begin_index..@p
+   * block_end_index] (inclusive) with tree nodes, writing them to
+   * @p out_nodes from left to right.
+   * @return The number of nodes written.
+   */
+  inline size_t cover_blocks_collect(const size_t& block_begin_index,
+                                     const size_t& block_end_index,
+                                     size_t (&out_nodes)[64]) const noexcept {
+    if (leaf_count == 0 || block_begin_index > block_end_index) {
+      return 0;
+    }
+    size_t left_index = first_leaf_index + block_begin_index;
+    size_t right_index = first_leaf_index + block_end_index;
+    size_t left_nodes[32];
+    size_t right_nodes[32];
+    size_t left_count = 0, right_count = 0;
+    while (left_index <= right_index) {
+      if (left_index & 1) {
+        left_nodes[left_count++] = left_index++;
+      }
+      if ((right_index & 1) == 0) {
+        right_nodes[right_count++] = right_index--;
       }
-      ++l;
+      left_index >>= 1;
+      right_index >>= 1;
     }
-    if (mn == INT_MAX) {
-      mn = cnt = 0;
+    size_t out_count = 0;
+    for (size_t i = 0; i < left_count; ++i) {
+      out_nodes[out_count++] = left_nodes[i];
     }
+    while (right_count > 0) {
+      out_nodes[out_count++] = right_nodes[--right_count];
+    }
+    return out_count;
   }
 
   /**
-   * @brief Select q-th minimum (1-based) inside [l..r] inclusive in two 8-bit
-   * passes.
-   * @details First pass finds global minimum, second selects the q-th position.
+   * @brief Select @p target_min_rank-th minimum (1-based) inside [@p
+   * range_begin..@p range_end] inclusive in two 8-bit passes.
+   * @details First pass finds global minimum, second selects the requested
+   * position.
    */
-  inline size_t qth_min_in_block(const size_t& l,
-                                 const size_t& r,
-                                 size_t q) const noexcept {
-    if (r < l || q == 0) {
+  inline size_t qth_min_in_block(const size_t& range_begin,
+                                 const size_t& range_end,
+                                 size_t target_min_rank) const noexcept {
+    if (range_end < range_begin || target_min_rank == 0) {
       return npos;
     }
 
-    const auto& T = LUT8();
+    const auto& aggregates_table = LUT8();
 
-    int cur = 0, mn = INT_MAX;
-    size_t p = l;
+    int current_excess = 0, min_value = INT_MAX;
+    size_t position = range_begin;
 
-    while (p <= r && (p & 7)) {
-      cur += bit(p) ? +1 : -1;
-      if (cur < mn) {
-        mn = cur;
+    while (position <= range_end && (position & 7)) {
+      current_excess += bit(position) ? +1 : -1;
+      if (current_excess < min_value) {
+        min_value = current_excess;
       }
-      ++p;
+      ++position;
     }
-    while (p + 7 <= r) {
-      const auto& a = T[get_byte(p)];
-      mn = std::min(mn, cur + a.min_prefix);
-      cur += a.excess_total;
-      p += 8;
+    while (position + 7 <= range_end) {
+      const auto& byte_aggregate = aggregates_table[get_byte(position)];
+      min_value =
+          std::min(min_value, current_excess + byte_aggregate.min_prefix);
+      current_excess += byte_aggregate.excess_total;
+      position += 8;
     }
-    while (p <= r) {
-      cur += bit(p) ? +1 : -1;
-      if (cur < mn) {
-        mn = cur;
+    while (position <= range_end) {
+      current_excess += bit(position) ? +1 : -1;
+      if (current_excess < min_value) {
+        min_value = current_excess;
       }
-      ++p;
+      ++position;
     }
 
-    cur = 0;
-    p = l;
+    current_excess = 0;
+    position = range_begin;
 
     // to byte alignment
-    while (p <= r && (p & 7)) {
-      cur += bit(p) ? +1 : -1;
-      if (cur == mn) {
-        if (--q == 0) {
-          return p;
+    while (position <= range_end && (position & 7)) {
+      current_excess += bit(position) ? +1 : -1;
+      if (current_excess == min_value) {
+        if (--target_min_rank == 0) {
+          return position;
         }
       }
-      ++p;
+      ++position;
     }
 
     // full bytes
-    while (p + 7 <= r) {
-      const uint8_t b = get_byte(p);
-      const auto& a = T[b];
-      const int cand = cur + a.min_prefix;
-      if (cand == mn) {
-        int s = 0;
+    while (position + 7 <= range_end) {
+      const uint8_t byte_value = get_byte(position);
+      const auto& byte_aggregate = aggregates_table[byte_value];
+      const int candidate = current_excess + byte_aggregate.min_prefix;
+      if (candidate == min_value) {
+        int prefix_sum = 0;
         for (int k = 0; k < 8; ++k) {
-          s += ((b >> k) & 1u) ? +1 : -1;
-          if (s == a.min_prefix) {
-            if (--q == 0) {
-              return p + k;
+          prefix_sum += ((byte_value >> k) & 1u) ? +1 : -1;
+          if (prefix_sum == byte_aggregate.min_prefix) {
+            if (--target_min_rank == 0) {
+              return position + k;
             }
           }
         }
       }
-      cur += a.excess_total;
-      p += 8;
+      current_excess += byte_aggregate.excess_total;
+      position += 8;
     }
 
     // tail
-    while (p <= r) {
-      cur += bit(p) ? +1 : -1;
-      if (cur == mn) {
-        if (--q == 0) {
-          return p;
+    while (position <= range_end) {
+      current_excess += bit(position) ? +1 : -1;
+      if (current_excess == min_value) {
+        if (--target_min_rank == 0) {
+          return position;
         }
       }
-      ++p;
+      ++position;
     }
 
     return npos;
   }
 
   /**
-   * @brief Find first minimum value and its first position in [l..r] inclusive
-   * using 8-bit LUT.
-   * @param mn_out Output minimum value (0 if empty).
-   * @param first_pos Output position of first minimum (npos if none).
+   * @brief Performs a forward search within a single leaf.
+   * @details Starting from position @p start_position (global index) inside the
+   * leaf that begins at @p leaf_block_begin, it looks for the nearest position
+   * where the excess changes by @p delta relative to @p start_position. The
+   * function uses a SIMD-accelerated search over the leaf’s prefix-sum array.
+   * If a matching position is found, returns its global index. If not found,
+   * returns npos and sets @p leaf_delta to the net excess change from @p
+   * start_position to the end of the leaf.
    */
-  inline void first_min_value_pos8(size_t l,
-                                   const size_t& r,
-                                   int& mn_out,
-                                   size_t& first_pos) const noexcept {
-    const auto& T = LUT8();
-    int cur = 0;
-    int mn = INT_MAX;
-    first_pos = npos;
+  inline size_t leaf_fwd_bp_simd(const size_t& leaf_index,
+                                 const size_t& leaf_block_begin,
+                                 const size_t& start_position,
+                                 const int& delta,
+                                 int& leaf_delta) const noexcept {
+    const size_t leaf_length = segment_size_bits[first_leaf_index + leaf_index];
+    const int16_t* leaf_prefix_ptr =
+        &leaf_prefix[leaf_index * leaf_prefix_stride];
+    const size_t offset_in_leaf = start_position - leaf_block_begin;
+    if (offset_in_leaf > leaf_length) {
+      leaf_delta = 0;
+      return npos;
+    }
+    const int16_t start_prefix = leaf_prefix_ptr[offset_in_leaf];
+    size_t match_index = ::find_forward_equal_i16_avx2(
+        leaf_prefix_ptr, offset_in_leaf + 1, leaf_length + 1,
+        start_prefix + delta, npos);
+    if (match_index != npos) {
+      return leaf_block_begin + match_index - 1;
+    }
+    leaf_delta = leaf_prefix_ptr[leaf_length] - start_prefix;
+    return npos;
+  }
+
+  /**
+   * @brief Performs a backward search within a leaf segment for a position
+   * earlier than @p start_position such that the excess difference from @p
+   * start_position equals @p delta.
+   * @details Uses SIMD (AVX2) to scan the leaf’s prefix-sum array of excess
+   * values. Returns the absolute bit position if found; otherwise returns npos
+   * and outputs the leaf-local excess delta at @p start_position.
+   */
+  inline size_t leaf_bwd_bp_simd(const size_t& leaf_index,
+                                 const size_t& leaf_block_begin,
+                                 const size_t& start_position,
+                                 const int& delta,
+                                 int& leaf_delta) const noexcept {
+    const int16_t* leaf_prefix_ptr =
+        &leaf_prefix[leaf_index * leaf_prefix_stride];
+    const size_t offset_in_leaf = start_position - leaf_block_begin;
+    if (offset_in_leaf > segment_size_bits[first_leaf_index + leaf_index]) {
+      leaf_delta = 0;
+      return npos;
+    }
+    const int16_t start_prefix = leaf_prefix_ptr[offset_in_leaf];
+    if (offset_in_leaf > 0) {
+      const size_t match_index = ::find_backward_equal_i16_avx2(
+          leaf_prefix_ptr, 0, offset_in_leaf - 1, start_prefix + delta, npos);
+      if (match_index != npos) {
+        return leaf_block_begin + match_index;
+      }
+    }
+    leaf_delta = start_prefix;
+    return npos;
+  }
+
+  /**
+   * @brief Find first minimum value and its first position in [@p range_begin
+   * ..@p range_end] inclusive using 8-bit LUT.
+   * @param min_value_out Output minimum value (0 if empty).
+   * @param first_position Output position of first minimum (npos if none).
+   */
+  inline void first_min_value_pos8(size_t range_begin,
+                                   const size_t& range_end,
+                                   int& min_value_out,
+                                   size_t& first_position) const noexcept {
+    const auto& aggregates_table = LUT8();
+    int current_excess = 0;
+    int min_value = INT_MAX;
+    first_position = npos;
 
     // to byte allignment
-    while (l <= r && (l & 7)) {
-      cur += bit(l) ? +1 : -1;
-      if (cur < mn) {
-        mn = cur;
-        first_pos = l;
+    while (range_begin <= range_end && (range_begin & 7)) {
+      current_excess += bit(range_begin) ? +1 : -1;
+      if (current_excess < min_value) {
+        min_value = current_excess;
+        first_position = range_begin;
       }
-      ++l;
+      ++range_begin;
     }
 
     // full bytes
-    while (l + 7 <= r) {
-      const auto& a = T[get_byte(l)];
-      const int cand = cur + a.min_prefix;
-      if (cand < mn) {
-        mn = cand;
-        first_pos = l + a.pos_first_min;
+    while (range_begin + 7 <= range_end) {
+      const auto& byte_aggregate = aggregates_table[get_byte(range_begin)];
+      const int candidate = current_excess + byte_aggregate.min_prefix;
+      if (candidate < min_value) {
+        min_value = candidate;
+        first_position = range_begin + byte_aggregate.pos_first_min;
       }
-      cur += a.excess_total;
-      l += 8;
+      current_excess += byte_aggregate.excess_total;
+      range_begin += 8;
     }
 
     // tail
-    while (l <= r) {
-      cur += bit(l) ? +1 : -1;
-      if (cur < mn) {
-        mn = cur;
-        first_pos = l;
+    while (range_begin <= range_end) {
+      current_excess += bit(range_begin) ? +1 : -1;
+      if (current_excess < min_value) {
+        min_value = current_excess;
+        first_position = range_begin;
       }
-      ++l;
+      ++range_begin;
     }
 
-    mn_out = (mn == INT_MAX ? 0 : mn);
+    min_value_out = (min_value == INT_MAX ? 0 : min_value);
   }
 
   /**
-   * @brief Find first maximum value and its first position in [l..r] inclusive
-   * using 8-bit LUT.
-   * @param mx_out Output maximum value (0 if empty).
-   * @param first_pos Output position of first maximum (npos if none).
+   * @brief Find first maximum value and its first position in [@p range_begin
+   * ..@p range_end] inclusive using 8-bit LUT.
+   * @param max_value_out Output maximum value (0 if empty).
+   * @param first_position Output position of first maximum (npos if none).
    */
-  inline void first_max_value_pos8(size_t l,
-                                   const size_t& r,
-                                   int& mx_out,
-                                   size_t& first_pos) const noexcept {
-    const auto& T = LUT8();
-    int cur = 0;
-    int mx = INT_MIN;
-    first_pos = npos;
+  inline void first_max_value_pos8(size_t range_begin,
+                                   const size_t& range_end,
+                                   int& max_value_out,
+                                   size_t& first_position) const noexcept {
+    const auto& aggregates_table = LUT8();
+    int current_excess = 0;
+    int max_value = INT_MIN;
+    first_position = npos;
 
-    while (l <= r && (l & 7)) {
-      cur += bit(l) ? +1 : -1;
-      if (cur > mx) {
-        mx = cur;
-        first_pos = l;
+    while (range_begin <= range_end && (range_begin & 7)) {
+      current_excess += bit(range_begin) ? +1 : -1;
+      if (current_excess > max_value) {
+        max_value = current_excess;
+        first_position = range_begin;
       }
-      ++l;
+      ++range_begin;
     }
 
-    while (l + 7 <= r) {
-      const auto& a = T[get_byte(l)];
-      const int cand = cur + a.max_prefix;
-      if (cand > mx) {
-        mx = cand;
-        first_pos = l + a.pos_first_max;
+    while (range_begin + 7 <= range_end) {
+      const auto& byte_aggregate = aggregates_table[get_byte(range_begin)];
+      const int candidate = current_excess + byte_aggregate.max_prefix;
+      if (candidate > max_value) {
+        max_value = candidate;
+        first_position = range_begin + byte_aggregate.pos_first_max;
       }
-      cur += a.excess_total;
-      l += 8;
+      current_excess += byte_aggregate.excess_total;
+      range_begin += 8;
     }
 
-    while (l <= r) {
-      cur += bit(l) ? +1 : -1;
-      if (cur > mx) {
-        mx = cur;
-        first_pos = l;
+    while (range_begin <= range_end) {
+      current_excess += bit(range_begin) ? +1 : -1;
+      if (current_excess > max_value) {
+        max_value = current_excess;
+        first_position = range_begin;
       }
-      ++l;
+      ++range_begin;
     }
 
-    mx_out = (mx == INT_MIN ? 0 : mx);
+    max_value_out = (max_value == INT_MIN ? 0 : max_value);
   }
 
   /**
@@ -1941,8 +2270,8 @@ class RmMTree {
 #endif
 
     leaf_count = ceil_div(num_bits, block_bits);
-    const size_t leaf0 = first_leaf_index();
-    const size_t tree_size = leaf0 + leaf_count - 1;
+    first_leaf_index = std::bit_ceil(std::max<size_t>(1, leaf_count));
+    const size_t tree_size = first_leaf_index + leaf_count - 1;
     segment_size_bits.assign(tree_size + 1, 0);
     node_total_excess.assign(tree_size + 1, 0);
     node_min_prefix_excess.assign(tree_size + 1, 0);
@@ -1951,135 +2280,167 @@ class RmMTree {
     node_pattern10_count.assign(tree_size + 1, 0);
     node_first_bit.assign(tree_size + 1, 0);
     node_last_bit.assign(tree_size + 1, 0);
+    leaf_prefix_stride = block_bits + 1;
+    leaf_prefix.assign(leaf_count * leaf_prefix_stride, 0);
 
     // leaves
-    for (size_t k = 0; k < leaf_count; ++k) {
-      const size_t v = leaf0 + k;
-      const size_t Lb = k * block_bits;
-      const size_t Rb = std::min(num_bits, Lb + block_bits);
-      segment_size_bits[v] = Rb - Lb;
+    for (size_t leaf_block_index = 0; leaf_block_index < leaf_count;
+         ++leaf_block_index) {
+      const size_t leaf_node_index = first_leaf_index + leaf_block_index;
+      const size_t segment_begin = leaf_block_index * block_bits;
+      const size_t segment_end = std::min(num_bits, segment_begin + block_bits);
+      segment_size_bits[leaf_node_index] = segment_end - segment_begin;
 
-      if (Lb < Rb) {
-        node_first_bit[v] = bit(Lb);
+      if (segment_begin < segment_end) {
+        node_first_bit[leaf_node_index] = bit(segment_begin);
       }
 
-      const auto& T = LUT8();
+      const auto& aggregates_table = LUT8();
 
-      int cur = 0, mn = INT_MAX, mx = INT_MIN;
-      uint32_t mn_cnt = 0;
-      uint32_t rrc = 0;
+      int current_excess = 0, min_value = INT_MAX, max_value = INT_MIN;
+      uint32_t min_count = 0;
+      uint32_t pattern10_count = 0;
 
-      uint8_t prev_bit = 0;
+      uint8_t previous_bit = 0;
 
-      size_t p = Lb;
+      size_t position = segment_begin;
 
       // Full bytes
-      while (p + 8 <= Rb) {
-        const uint8_t b = get_byte(p);
-        const auto& a = T[b];
+      while (position + 8 <= segment_end) {
+        const uint8_t byte_value = get_byte(position);
+        const auto& byte_aggregate = aggregates_table[byte_value];
 
         // internal "10" inside the byte
-        rrc += a.pattern10_count;
+        pattern10_count += byte_aggregate.pattern10_count;
         // stitching across the boundary between the previous and current byte
         // (within the segment)
-        if (prev_bit == 1 && a.first_bit == 0) {
-          rrc++;
+        if (previous_bit == 1 && byte_aggregate.first_bit == 0) {
+          pattern10_count++;
         }
 
         // prefix min/max accounting for the current offset
-        const int cand_min = cur + a.min_prefix;
-        if (cand_min < mn) {
-          mn = cand_min;
-          mn_cnt = a.min_count;
-        } else if (cand_min == mn) {
-          mn_cnt += a.min_count;
+        const int candidate_min = current_excess + byte_aggregate.min_prefix;
+        if (candidate_min < min_value) {
+          min_value = candidate_min;
+          min_count = byte_aggregate.min_count;
+        } else if (candidate_min == min_value) {
+          min_count += byte_aggregate.min_count;
         }
 
-        mx = std::max(mx, cur + a.max_prefix);
-        cur += a.excess_total;
-        prev_bit = a.last_bit;
-        p += 8;
+        max_value =
+            std::max(max_value, current_excess + byte_aggregate.max_prefix);
+        current_excess += byte_aggregate.excess_total;
+        previous_bit = byte_aggregate.last_bit;
+        position += 8;
       }
 
       // Tail < 8 bits
-      while (p < Rb) {
-        const uint8_t b = bit(p);
-        if (prev_bit == 1 && b == 0) {
-          rrc++;
+      while (position < segment_end) {
+        const uint8_t bit_value = bit(position);
+        if (previous_bit == 1 && bit_value == 0) {
+          pattern10_count++;
         }
-        const int step = b ? +1 : -1;
-        cur += step;
-        if (cur < mn) {
-          mn = cur;
-          mn_cnt = 1;
-        } else if (cur == mn) {
-          ++mn_cnt;
+        const int step = bit_value ? +1 : -1;
+        current_excess += step;
+        if (current_excess < min_value) {
+          min_value = current_excess;
+          min_count = 1;
+        } else if (current_excess == min_value) {
+          ++min_count;
         }
-        if (cur > mx) {
-          mx = cur;
+        if (current_excess > max_value) {
+          max_value = current_excess;
         }
 
-        prev_bit = b;
-        ++p;
+        previous_bit = bit_value;
+        ++position;
       }
 
-      if (Lb < Rb) {
-        node_last_bit[v] = prev_bit;
+      if (segment_begin < segment_end) {
+        node_last_bit[leaf_node_index] = previous_bit;
       }
 
-      node_total_excess[v] = cur;
-      node_min_prefix_excess[v] = (segment_size_bits[v] == 0 ? 0 : mn);
-      node_max_prefix_excess[v] = (segment_size_bits[v] == 0 ? 0 : mx);
-      node_min_count[v] = mn_cnt;
-      node_pattern10_count[v] = (uint32_t)rrc;
+      node_total_excess[leaf_node_index] = current_excess;
+      node_min_prefix_excess[leaf_node_index] =
+          (segment_size_bits[leaf_node_index] == 0 ? 0 : min_value);
+      node_max_prefix_excess[leaf_node_index] =
+          (segment_size_bits[leaf_node_index] == 0 ? 0 : max_value);
+      node_min_count[leaf_node_index] = min_count;
+      node_pattern10_count[leaf_node_index] = (uint32_t)pattern10_count;
+      int16_t* leaf_prefix_ptr =
+          &leaf_prefix[leaf_block_index * leaf_prefix_stride];
+      int16_t prefix_accumulator = 0;
+      leaf_prefix_ptr[0] = 0;
+      for (size_t position_in_leaf = segment_begin, prefix_index = 1;
+           position_in_leaf < segment_end; ++position_in_leaf, ++prefix_index) {
+        prefix_accumulator += bit(position_in_leaf) ? 1 : -1;
+        leaf_prefix_ptr[prefix_index] = prefix_accumulator;
+      }
     }
     // internal nodes
-    for (size_t v = leaf0 - 1; v >= 1; --v) {
-      const size_t Lc = v << 1;
-      const size_t Rc = Lc | 1;
-      const bool has_l = (Lc <= tree_size) && segment_size_bits[Lc];
-      const bool has_r = (Rc <= tree_size) && segment_size_bits[Rc];
-      if (!has_l && !has_r) {
-        segment_size_bits[v] = 0;
+    for (size_t node_index = first_leaf_index - 1; node_index >= 1;
+         --node_index) {
+      const size_t left_child = node_index << 1;
+      const size_t right_child = left_child | 1;
+      const bool has_left =
+          (left_child <= tree_size) && segment_size_bits[left_child];
+      const bool has_right =
+          (right_child <= tree_size) && segment_size_bits[right_child];
+      if (!has_left && !has_right) {
+        segment_size_bits[node_index] = 0;
         continue;
       }
-      if (has_l && !has_r) {
-        segment_size_bits[v] = segment_size_bits[Lc];
-        node_total_excess[v] = node_total_excess[Lc];
-        node_min_prefix_excess[v] = node_min_prefix_excess[Lc];
-        node_max_prefix_excess[v] = node_max_prefix_excess[Lc];
-        node_min_count[v] = node_min_count[Lc];
-        node_pattern10_count[v] = node_pattern10_count[Lc];
-        node_first_bit[v] = node_first_bit[Lc];
-        node_last_bit[v] = node_last_bit[Lc];
-      } else if (!has_l && has_r) {
-        segment_size_bits[v] = segment_size_bits[Rc];
-        node_total_excess[v] = node_total_excess[Rc];
-        node_min_prefix_excess[v] = node_min_prefix_excess[Rc];
-        node_max_prefix_excess[v] = node_max_prefix_excess[Rc];
-        node_min_count[v] = node_min_count[Rc];
-        node_pattern10_count[v] = node_pattern10_count[Rc];
-        node_first_bit[v] = node_first_bit[Rc];
-        node_last_bit[v] = node_last_bit[Rc];
+      if (has_left && !has_right) {
+        segment_size_bits[node_index] = segment_size_bits[left_child];
+        node_total_excess[node_index] = node_total_excess[left_child];
+        node_min_prefix_excess[node_index] = node_min_prefix_excess[left_child];
+        node_max_prefix_excess[node_index] = node_max_prefix_excess[left_child];
+        node_min_count[node_index] = node_min_count[left_child];
+        node_pattern10_count[node_index] = node_pattern10_count[left_child];
+        node_first_bit[node_index] = node_first_bit[left_child];
+        node_last_bit[node_index] = node_last_bit[left_child];
+      } else if (!has_left && has_right) {
+        segment_size_bits[node_index] = segment_size_bits[right_child];
+        node_total_excess[node_index] = node_total_excess[right_child];
+        node_min_prefix_excess[node_index] =
+            node_min_prefix_excess[right_child];
+        node_max_prefix_excess[node_index] =
+            node_max_prefix_excess[right_child];
+        node_min_count[node_index] = node_min_count[right_child];
+        node_pattern10_count[node_index] = node_pattern10_count[right_child];
+        node_first_bit[node_index] = node_first_bit[right_child];
+        node_last_bit[node_index] = node_last_bit[right_child];
       } else {
-        segment_size_bits[v] = segment_size_bits[Lc] + segment_size_bits[Rc];
-        node_total_excess[v] = node_total_excess[Lc] + node_total_excess[Rc];
-        const int m_r = node_total_excess[Lc] + node_min_prefix_excess[Rc];
-        const int M_R = node_total_excess[Lc] + node_max_prefix_excess[Rc];
-        node_min_prefix_excess[v] = std::min(node_min_prefix_excess[Lc], m_r);
-        node_max_prefix_excess[v] = std::max(node_max_prefix_excess[Lc], M_R);
-        node_min_count[v] =
-            (node_min_prefix_excess[Lc] == node_min_prefix_excess[v]
-                 ? node_min_count[Lc]
+        segment_size_bits[node_index] =
+            segment_size_bits[left_child] + segment_size_bits[right_child];
+        node_total_excess[node_index] =
+            node_total_excess[left_child] + node_total_excess[right_child];
+        const int right_min_candidate =
+            node_total_excess[left_child] + node_min_prefix_excess[right_child];
+        const int right_max_candidate =
+            node_total_excess[left_child] + node_max_prefix_excess[right_child];
+        node_min_prefix_excess[node_index] =
+            std::min(node_min_prefix_excess[left_child], right_min_candidate);
+        node_max_prefix_excess[node_index] =
+            std::max(node_max_prefix_excess[left_child], right_max_candidate);
+        node_min_count[node_index] =
+            (node_min_prefix_excess[left_child] ==
+                     node_min_prefix_excess[node_index]
+                 ? node_min_count[left_child]
                  : 0) +
-            (m_r == node_min_prefix_excess[v] ? node_min_count[Rc] : 0);
-        node_pattern10_count[v] =
-            node_pattern10_count[Lc] + node_pattern10_count[Rc] +
-            ((node_last_bit[Lc] == 1 && node_first_bit[Rc] == 0) ? 1u : 0u);
-        node_first_bit[v] = node_first_bit[Lc];
-        node_last_bit[v] = node_last_bit[Rc];
-      }
-      if (v == 1) {
+            (right_min_candidate == node_min_prefix_excess[node_index]
+                 ? node_min_count[right_child]
+                 : 0);
+        node_pattern10_count[node_index] = node_pattern10_count[left_child] +
+                                           node_pattern10_count[right_child] +
+                                           ((node_last_bit[left_child] == 1 &&
+                                             node_first_bit[right_child] == 0)
+                                                ? 1u
+                                                : 0u);
+        node_first_bit[node_index] = node_first_bit[left_child];
+        node_last_bit[node_index] = node_last_bit[right_child];
+      }
+      if (node_index == 1) {
         break;
       }
     }
diff --git a/misc/plot_rmm.py b/misc/plot_rmm.py
index 881d50c..78e7beb 100644
--- a/misc/plot_rmm.py
+++ b/misc/plot_rmm.py
@@ -1,17 +1,18 @@
 """
 Plot RmMTree benchmark results.
 
-This script reads a CSV produced by `bench_rmm.cpp` and,
+This script reads a JSON produced by `bench_rmm` and,
 for each operation, draws a scatter plot of individual points
 and a trend line (optionally median-smoothed).
 Plots are saved as PNG files and can also be shown interactively.
 
 Examples:
-  python3 plot_rmm.py rmm_bench.csv --save-dir plots --logx --smooth 3
-  python3 plot_rmm.py rmm_bench.csv --show
+  python3 plot_rmm.py rmm_bench.json --save-dir plots --logx --smooth 3
+  python3 plot_rmm.py rmm_bench.json --show
 """
 
 import argparse
+import json
 import os
 import pandas as pd
 import matplotlib.pyplot as plt
@@ -41,20 +42,20 @@
 def main():
     ap = argparse.ArgumentParser(
         description=(
-            "Read a CSV with RmMTree benchmark results and plot time per operation "
+            "Read a JSON with RmMTree benchmark results and plot time per operation "
             "versus sequence size N for each operation."
         ),
         epilog=(
             "Examples:\n"
-            "  python3 plot_rmm.py rmm_bench.csv --save-dir plots --logx --smooth 3\n"
-            "  python3 plot_rmm.py rmm_bench.csv --show"
+            "  python3 plot_rmm.py rmm_bench.json --save-dir plots --logx --smooth 3\n"
+            "  python3 plot_rmm.py rmm_bench.json --show"
         ),
         formatter_class=argparse.RawDescriptionHelpFormatter,
     )
     ap.add_argument(
-        "csv",
-        metavar="CSV",
-        help="Path to the CSV file with results (output of bench_rmm.cpp).",
+        "json",
+        metavar="JSON",
+        help="Path to the JSON file with results (output of bench_rmm).",
     )
     ap.add_argument(
         "--save-dir",
@@ -77,7 +78,7 @@ def main():
         "--logx",
         action="store_true",
         help=(
-            "Use a logarithmic X axis (base 2). " "Handy when N grows in powers of two."
+            "Use a logarithmic X axis (base 2). Handy when N grows in powers of two."
         ),
     )
     ap.add_argument(
@@ -93,7 +94,11 @@ def main():
     args = ap.parse_args()
 
     os.makedirs(args.save_dir, exist_ok=True)
-    df = pd.read_csv(args.csv)
+    with open(args.json, "r", encoding="utf-8") as f:
+        data = json.load(f)
+    if isinstance(data, dict) and "benchmarks" in data:
+        data = data["benchmarks"]
+    df = pd.DataFrame(data)
     df = df.dropna(subset=["cpu_time", "N"])
 
     for op in OPS_ORDER:
diff --git a/src/bench_rmm.cpp b/src/bench_rmm.cpp
index 06db0db..523516e 100644
--- a/src/bench_rmm.cpp
+++ b/src/bench_rmm.cpp
@@ -3,6 +3,7 @@
 #include <algorithm>
 #include <cctype>
 #include <cmath>
+#include <cstddef>
 #include <cstdint>
 #include <iostream>
 #include <memory>
@@ -10,42 +11,42 @@
 #include <random>
 #include <set>
 #include <string>
+#include <utility>
 #include <vector>
 
 #include "rmm_tree.h"
 
-using namespace std;
 using pixie::RmMTree;
 
 struct Args {
   int min_exp = 14;
   int max_exp = 22;
-  size_t Q = 200000;
+  std::size_t Q = 200000;
   double p1 = 0.5;
-  uint64_t seed = 42;
-  size_t block_bits = 64;
+  std::uint64_t seed = 42;
+  std::size_t block_bits = 64;
   int per_octave = 10;
-  vector<size_t> explicit_sizes;
+  std::vector<std::size_t> explicit_sizes;
 };
 
 static Args args;
 
 static void parse_args_and_strip(int& argc, char**& argv) {
-  auto getv = [&](const string& key) -> optional<string> {
-    string pref = "--" + key + "=";
+  auto getv = [&](const std::string& key) -> std::optional<std::string> {
+    std::string pref = "--" + key + "=";
     for (int i = 1; i < argc; ++i) {
-      string s = argv[i];
+      std::string s = argv[i];
       if (s.rfind(pref, 0) == 0) {
         return s.substr(pref.size());
       }
     }
-    return nullopt;
+    return std::nullopt;
   };
-  auto strip = [&](const string& key) {
-    string pref = "--" + key + "=";
+  auto strip = [&](const std::string& key) {
+    std::string pref = "--" + key + "=";
     int w = 0;
     for (int i = 0; i < argc; ++i) {
-      if (i > 0 && string(argv[i]).rfind(pref, 0) == 0) {
+      if (i > 0 && std::string(argv[i]).rfind(pref, 0) == 0) {
         continue;
       }
       argv[w++] = argv[i];
@@ -54,85 +55,88 @@ static void parse_args_and_strip(int& argc, char**& argv) {
   };
 
   if (auto v = getv("min-exp")) {
-    args.min_exp = stoi(*v), strip("min-exp");
+    args.min_exp = std::stoi(*v), strip("min-exp");
   }
   if (auto v = getv("max-exp")) {
-    args.max_exp = stoi(*v), strip("max-exp");
+    args.max_exp = std::stoi(*v), strip("max-exp");
   }
   if (auto v = getv("q")) {
-    args.Q = stoull(*v), strip("q");
+    args.Q = std::stoull(*v), strip("q");
   }
   if (auto v = getv("p")) {
-    args.p1 = stod(*v), strip("p");
+    args.p1 = std::stod(*v), strip("p");
   }
   if (auto v = getv("seed")) {
-    args.seed = stoull(*v), strip("seed");
+    args.seed = std::stoull(*v), strip("seed");
   }
   if (auto v = getv("block")) {
-    args.block_bits = stoull(*v), strip("block");
+    args.block_bits = std::stoull(*v), strip("block");
   }
   if (auto v = getv("per-octave")) {
-    args.per_octave = stoi(*v), strip("per-octave");
+    args.per_octave = std::stoi(*v), strip("per-octave");
   }
 
   if (auto v = getv("sizes")) {
-    string s = *v;
+    std::string s = *v;
     strip("sizes");
-    size_t pos = 0;
+    std::size_t pos = 0;
     while (pos < s.size()) {
       while (pos < s.size() &&
-             (s[pos] == ',' || isspace((unsigned char)s[pos]))) {
+             (s[pos] == ',' || std::isspace((unsigned char)s[pos]))) {
         ++pos;
       }
-      size_t start = pos;
-      while (pos < s.size() && isdigit((unsigned char)s[pos])) {
+      std::size_t start = pos;
+      while (pos < s.size() && std::isdigit((unsigned char)s[pos])) {
         ++pos;
       }
       if (start < pos) {
-        args.explicit_sizes.push_back(stoull(s.substr(start, pos - start)));
+        args.explicit_sizes.push_back(
+            std::stoull(s.substr(start, pos - start)));
       }
     }
-    sort(args.explicit_sizes.begin(), args.explicit_sizes.end());
+    std::sort(args.explicit_sizes.begin(), args.explicit_sizes.end());
     args.explicit_sizes.erase(
-        unique(args.explicit_sizes.begin(), args.explicit_sizes.end()),
+        std::unique(args.explicit_sizes.begin(), args.explicit_sizes.end()),
         args.explicit_sizes.end());
   }
 }
 
-static string make_random_bits(size_t N, double p1, mt19937_64& rng) {
-  uniform_real_distribution<double> U(0.0, 1.0);
-  string s;
+static std::string make_random_bits(std::size_t N,
+                                    double p1,
+                                    std::mt19937_64& rng) {
+  std::uniform_real_distribution<double> U(0.0, 1.0);
+  std::string s;
   s.resize(N);
-  for (size_t i = 0; i < N; ++i) {
+  for (std::size_t i = 0; i < N; ++i) {
     s[i] = (U(rng) < p1 ? '1' : '0');
   }
   return s;
 }
 
-static vector<size_t> build_size_grid() {
+static std::vector<std::size_t> build_size_grid() {
   if (!args.explicit_sizes.empty()) {
-    vector<size_t> v = args.explicit_sizes;
-    v.erase(remove(v.begin(), v.end(), 0), v.end());
+    std::vector<std::size_t> v = args.explicit_sizes;
+    v.erase(std::remove(v.begin(), v.end(), 0), v.end());
     return v;
   }
 
   const int lo = args.min_exp, hi = args.max_exp;
   if (args.per_octave <= 0) {
-    vector<size_t> v;
+    std::vector<std::size_t> v;
     for (int e = lo; e <= hi; ++e) {
-      v.push_back(size_t(1) << e);
+      v.push_back(std::size_t(1) << e);
     }
     return v;
   }
 
   int per_octave_steps = args.per_octave;
-  set<size_t> N_cands;
-  size_t min_N = size_t(1) << lo, max_N = size_t(1) << hi;
+  std::set<std::size_t> N_cands;
+  std::size_t min_N = std::size_t(1) << lo, max_N = std::size_t(1) << hi;
 
   for (int e = lo; e <= hi; ++e) {
     for (int t = 0; t <= per_octave_steps; ++t) {
       long double x = e + (long double)t / (long double)per_octave_steps;
-      size_t N = (size_t)llround(powl(2.0L, x));
+      std::size_t N = (std::size_t)std::llround(std::pow(2.0L, x));
       if (N < min_N) {
         N = min_N;
       }
@@ -145,7 +149,7 @@ static vector<size_t> build_size_grid() {
     }
   }
 
-  vector<size_t> v(N_cands.begin(), N_cands.end());
+  std::vector<std::size_t> v(N_cands.begin(), N_cands.end());
   if (v.front() != min_N) {
     v.insert(v.begin(), min_N);
   }
@@ -156,33 +160,33 @@ static vector<size_t> build_size_grid() {
 }
 
 struct Pools {
-  vector<size_t> inds_any;
-  vector<size_t> inds;
-  vector<size_t> inds_1N;
-  vector<int> deltas;
-  vector<pair<size_t, size_t>> segs;
-
-  vector<size_t> ks1, ks0, ks10;
-  vector<size_t> minselect_q;
+  std::vector<std::size_t> inds_any;
+  std::vector<std::size_t> inds;
+  std::vector<std::size_t> inds_1N;
+  std::vector<int> deltas;
+  std::vector<std::pair<std::size_t, std::size_t>> segs;
+
+  std::vector<std::size_t> ks1, ks0, ks10;
+  std::vector<std::size_t> minselect_q;
 };
 
 struct Dataset {
-  size_t N{};
-  string bits;
+  std::size_t N{};
+  std::string bits;
   RmMTree t;
 
-  size_t cnt1{}, cnt0{}, cnt10{};
+  std::size_t cnt1{}, cnt0{}, cnt10{};
   Pools pool;
 };
 
-static vector<shared_ptr<Dataset>> keepalive;
+static std::vector<std::shared_ptr<Dataset>> keepalive;
 
-static size_t count10(const string& s) {
-  size_t c = 0;
+static std::size_t count10(const std::string& s) {
+  std::size_t c = 0;
   if (s.size() < 2) {
     return 0;
   }
-  for (size_t i = 0; i + 1 < s.size(); ++i) {
+  for (std::size_t i = 0; i + 1 < s.size(); ++i) {
     if (s[i] == '1' && s[i + 1] == '0') {
       ++c;
     }
@@ -190,8 +194,9 @@ static size_t count10(const string& s) {
   return c;
 }
 
-static Dataset build_dataset(size_t N) {
-  mt19937_64 rng(args.seed ^ (uint64_t)N * 0x9E3779B185EBCA87ull);
+static Dataset build_dataset(std::size_t N) {
+  std::mt19937_64 rng(args.seed ^
+                      static_cast<std::uint64_t>(N) * 0x9E3779B185EBCA87ull);
 
   Dataset d;
   d.N = N;
@@ -207,12 +212,12 @@ static Dataset build_dataset(size_t N) {
   d.cnt0 = N - d.cnt1;
   d.cnt10 = count10(d.bits);
 
-  const size_t L = min<size_t>(args.Q, 32768);
+  const std::size_t L = std::min<std::size_t>(args.Q, 32768);
 
-  uniform_int_distribution<size_t> ind_dist_incl(0, N ? N : 0);
-  uniform_int_distribution<size_t> ind_dist(0, N ? (N - 1) : 0);
-  uniform_int_distribution<size_t> ind_dist_1N(1, N ? N : 1);
-  uniform_int_distribution<int> d_dist(-8, +8);
+  std::uniform_int_distribution<std::size_t> ind_dist_incl(0, N ? N : 0);
+  std::uniform_int_distribution<std::size_t> ind_dist(0, N ? (N - 1) : 0);
+  std::uniform_int_distribution<std::size_t> ind_dist_1N(1, N ? N : 1);
+  std::uniform_int_distribution<int> d_dist(-8, +8);
 
   d.pool.inds_any.resize(L);
   d.pool.inds.resize(L);
@@ -220,18 +225,18 @@ static Dataset build_dataset(size_t N) {
   d.pool.deltas.resize(L);
   d.pool.segs.resize(L);
 
-  auto rand_ij = [&]() -> pair<size_t, size_t> {
+  auto rand_ij = [&]() -> std::pair<std::size_t, std::size_t> {
     if (N == 0) {
       return {0, 0};
     }
-    size_t a = ind_dist(rng), b = ind_dist(rng);
+    std::size_t a = ind_dist(rng), b = ind_dist(rng);
     if (a > b) {
-      swap(a, b);
+      std::swap(a, b);
     }
     return {a, b};
   };
 
-  for (size_t i = 0; i < L; ++i) {
+  for (std::size_t i = 0; i < L; ++i) {
     d.pool.inds_any[i] = ind_dist_incl(rng);
     d.pool.inds[i] = (N ? ind_dist(rng) : 0);
     d.pool.inds_1N[i] = (N ? ind_dist_1N(rng) : 0);
@@ -239,14 +244,14 @@ static Dataset build_dataset(size_t N) {
     d.pool.segs[i] = rand_ij();
   }
 
-  auto fill_ks = [&](size_t total, vector<size_t>& out) {
+  auto fill_ks = [&](std::size_t total, std::vector<std::size_t>& out) {
     out.resize(L);
     if (total == 0) {
-      fill(out.begin(), out.end(), 1);
+      std::fill(out.begin(), out.end(), 1);
       return;
     }
-    uniform_int_distribution<size_t> dist(1, total);
-    for (size_t i = 0; i < L; ++i) {
+    std::uniform_int_distribution<std::size_t> dist(1, total);
+    for (std::size_t i = 0; i < L; ++i) {
       out[i] = dist(rng);
     }
   };
@@ -255,13 +260,13 @@ static Dataset build_dataset(size_t N) {
   fill_ks(d.cnt10, d.pool.ks10);
 
   d.pool.minselect_q.resize(L);
-  for (size_t i = 0; i < L; ++i) {
+  for (std::size_t i = 0; i < L; ++i) {
     auto [l, r] = d.pool.segs[i];
-    size_t c = d.t.mincount(l, r);
+    std::size_t c = d.t.mincount(l, r);
     if (c == 0) {
       c = 1;
     }
-    uniform_int_distribution<size_t> uq(1, c);
+    std::uniform_int_distribution<std::size_t> uq(1, c);
     d.pool.minselect_q[i] = uq(rng);
   }
 
@@ -269,14 +274,16 @@ static Dataset build_dataset(size_t N) {
 }
 
 template <class Fn>
-static void register_op(const string& op, shared_ptr<Dataset> data, Fn&& body) {
-  auto idx_ptr = make_shared<size_t>(0);
+static void register_op(const std::string& op,
+                        std::shared_ptr<Dataset> data,
+                        Fn&& body) {
+  auto idx_ptr = std::make_shared<std::size_t>(0);
 
   auto* b = benchmark::RegisterBenchmark(
       op.c_str(), [data, idx_ptr, body](benchmark::State& state) {
         const Dataset& D = *data;
         for (auto _ : state) {
-          size_t i = (*idx_ptr)++;
+          std::size_t i = (*idx_ptr)++;
           body(state, D, i);
         }
         state.counters["N"] = static_cast<double>(D.N);
@@ -289,106 +296,106 @@ static void register_op(const string& op, shared_ptr<Dataset> data, Fn&& body) {
 
 static void register_all() {
   auto Ns = build_size_grid();
-  for (size_t N : Ns) {
-    auto data = make_shared<Dataset>(build_dataset(N));
+  for (std::size_t N : Ns) {
+    auto data = std::make_shared<Dataset>(build_dataset(N));
     keepalive.push_back(data);
 
     const auto& P = data->pool;
 
     register_op("rank1", data,
-                [&](benchmark::State&, const Dataset& D, size_t k) {
+                [&](benchmark::State&, const Dataset& D, std::size_t k) {
                   auto r = D.t.rank1(P.inds_any[k % P.inds_any.size()]);
                   benchmark::DoNotOptimize(r);
                 });
 
     register_op("rank0", data,
-                [&](benchmark::State&, const Dataset& D, size_t k) {
+                [&](benchmark::State&, const Dataset& D, std::size_t k) {
                   auto r = D.t.rank0(P.inds_any[k % P.inds_any.size()]);
                   benchmark::DoNotOptimize(r);
                 });
 
     register_op("select1", data,
-                [&](benchmark::State&, const Dataset& D, size_t k) {
+                [&](benchmark::State&, const Dataset& D, std::size_t k) {
                   auto r = D.t.select1(P.ks1[k % P.ks1.size()]);
                   benchmark::DoNotOptimize(r);
                 });
 
     register_op("select0", data,
-                [&](benchmark::State&, const Dataset& D, size_t k) {
+                [&](benchmark::State&, const Dataset& D, std::size_t k) {
                   auto r = D.t.select0(P.ks0[k % P.ks0.size()]);
                   benchmark::DoNotOptimize(r);
                 });
 
     register_op("rank10", data,
-                [&](benchmark::State&, const Dataset& D, size_t k) {
+                [&](benchmark::State&, const Dataset& D, std::size_t k) {
                   auto r = D.t.rank10(P.inds_any[k % P.inds_any.size()]);
                   benchmark::DoNotOptimize(r);
                 });
 
     register_op("select10", data,
-                [&](benchmark::State&, const Dataset& D, size_t k) {
+                [&](benchmark::State&, const Dataset& D, std::size_t k) {
                   auto r = D.t.select10(P.ks10[k % P.ks10.size()]);
                   benchmark::DoNotOptimize(r);
                 });
 
     register_op("excess", data,
-                [&](benchmark::State&, const Dataset& D, size_t k) {
+                [&](benchmark::State&, const Dataset& D, std::size_t k) {
                   auto r = D.t.excess(P.inds_any[k % P.inds_any.size()]);
                   benchmark::DoNotOptimize(r);
                 });
 
     register_op("fwdsearch", data,
-                [&](benchmark::State&, const Dataset& D, size_t k) {
+                [&](benchmark::State&, const Dataset& D, std::size_t k) {
                   auto r = D.t.fwdsearch(P.inds[k % P.inds.size()],
                                          P.deltas[k % P.deltas.size()]);
                   benchmark::DoNotOptimize(r);
                 });
 
     register_op("bwdsearch", data,
-                [&](benchmark::State&, const Dataset& D, size_t k) {
+                [&](benchmark::State&, const Dataset& D, std::size_t k) {
                   auto r = D.t.bwdsearch(P.inds_1N[k % P.inds_1N.size()],
                                          P.deltas[k % P.deltas.size()]);
                   benchmark::DoNotOptimize(r);
                 });
 
     register_op("range_min_query_pos", data,
-                [&](benchmark::State&, const Dataset& D, size_t k) {
+                [&](benchmark::State&, const Dataset& D, std::size_t k) {
                   auto [i, j] = P.segs[k % P.segs.size()];
                   auto r = D.t.range_min_query_pos(i, j);
                   benchmark::DoNotOptimize(r);
                 });
 
     register_op("range_min_query_val", data,
-                [&](benchmark::State&, const Dataset& D, size_t k) {
+                [&](benchmark::State&, const Dataset& D, std::size_t k) {
                   auto [i, j] = P.segs[k % P.segs.size()];
                   auto r = D.t.range_min_query_val(i, j);
                   benchmark::DoNotOptimize(r);
                 });
 
     register_op("range_max_query_pos", data,
-                [&](benchmark::State&, const Dataset& D, size_t k) {
+                [&](benchmark::State&, const Dataset& D, std::size_t k) {
                   auto [i, j] = P.segs[k % P.segs.size()];
                   auto r = D.t.range_max_query_pos(i, j);
                   benchmark::DoNotOptimize(r);
                 });
 
     register_op("range_max_query_val", data,
-                [&](benchmark::State&, const Dataset& D, size_t k) {
+                [&](benchmark::State&, const Dataset& D, std::size_t k) {
                   auto [i, j] = P.segs[k % P.segs.size()];
                   auto r = D.t.range_max_query_val(i, j);
                   benchmark::DoNotOptimize(r);
                 });
 
     register_op("mincount", data,
-                [&](benchmark::State&, const Dataset& D, size_t k) {
+                [&](benchmark::State&, const Dataset& D, std::size_t k) {
                   auto [i, j] = P.segs[k % P.segs.size()];
                   auto r = D.t.mincount(i, j);
                   benchmark::DoNotOptimize(r);
                 });
 
     register_op("minselect", data,
-                [&](benchmark::State&, const Dataset& D, size_t k) {
-                  const size_t idx = k % P.segs.size();
+                [&](benchmark::State&, const Dataset& D, std::size_t k) {
+                  const std::size_t idx = k % P.segs.size();
                   auto [i, j] = P.segs[idx];
                   auto q = P.minselect_q[idx];
                   auto r = D.t.minselect(i, j, q);
@@ -396,7 +403,7 @@ static void register_all() {
                 });
 
     register_op("close", data,
-                [&](benchmark::State&, const Dataset& D, size_t k) {
+                [&](benchmark::State&, const Dataset& D, std::size_t k) {
                   if (D.N == 0) {
                     benchmark::DoNotOptimize(0);
                     return;
@@ -407,13 +414,13 @@ static void register_all() {
                 });
 
     register_op("open", data,
-                [&](benchmark::State&, const Dataset& D, size_t k) {
+                [&](benchmark::State&, const Dataset& D, std::size_t k) {
                   auto r = D.t.open(P.inds_1N[k % P.inds_1N.size()]);
                   benchmark::DoNotOptimize(r);
                 });
 
     register_op("enclose", data,
-                [&](benchmark::State&, const Dataset& D, size_t k) {
+                [&](benchmark::State&, const Dataset& D, std::size_t k) {
                   auto r = D.t.enclose(P.inds_1N[k % P.inds_1N.size()]);
                   benchmark::DoNotOptimize(r);
                 });
@@ -421,12 +428,13 @@ static void register_all() {
 }
 
 int main(int argc, char** argv) {
+  benchmark::MaybeReenterWithoutASLR(argc, argv);
   parse_args_and_strip(argc, argv);
 
   auto has = [&](const char* key) {
-    string p1 = string(key) + "=";
+    std::string p1 = std::string(key) + "=";
     for (int i = 1; i < argc; ++i) {
-      string s = argv[i];
+      std::string s = argv[i];
       if (s == key || s.rfind(p1, 0) == 0) {
         return true;
       }
@@ -434,9 +442,9 @@ int main(int argc, char** argv) {
     return false;
   };
 
-  static vector<string> extra;
+  static std::vector<std::string> extra;
   if (!has("--benchmark_out_format")) {
-    extra.emplace_back("--benchmark_out_format=csv");
+    extra.emplace_back("--benchmark_out_format=json");
   }
   if (!has("--benchmark_counters_tabular")) {
     extra.emplace_back("--benchmark_counters_tabular=true");
@@ -445,7 +453,7 @@ int main(int argc, char** argv) {
     extra.emplace_back("--benchmark_time_unit=ns");
   }
 
-  static vector<char*> argv_vec;
+  static std::vector<char*> argv_vec;
   argv_vec.assign(argv, argv + argc);
   for (auto& s : extra) {
     argv_vec.push_back(s.data());