diff --git a/.gitignore b/.gitignore index b3ccf80..f0a51cf 100644 --- a/.gitignore +++ b/.gitignore @@ -2,3 +2,4 @@ build/ .vscode/ Testing/ plans/* +venv/ diff --git a/CMakeLists.txt b/CMakeLists.txt index b520326..d24fde2 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -45,4 +45,24 @@ target_link_libraries(unittests gtest gtest_main) +add_executable(bench_rmm + src/bench_rmm.cpp) +target_include_directories(bench_rmm + PUBLIC include +) +set_property(TARGET bench_rmm PROPERTY CXX_STANDARD 20) +target_link_libraries(bench_rmm + benchmark) + +add_executable(test_rmm + src/test_rmm.cpp) +target_include_directories(test_rmm + PUBLIC include + PUBLIC include/misc +) +set_property(TARGET test_rmm PROPERTY CXX_STANDARD 20) +target_link_libraries(test_rmm + gtest + gtest_main) + enable_testing() \ No newline at end of file diff --git a/README.md b/README.md index 658114b..564058f 100644 --- a/README.md +++ b/README.md @@ -40,20 +40,38 @@ This will build the library along with benchmarks and tests. After building: +### BitVector + +```bash +./unittests +``` + +### RmM Tree + ```bash -/unittests +./test_rmm ``` --- ## Running Benchmarks -Bencharks are random 50/50 0-1 bitvectors up to $2^34$ bits. +### BitVector + +Benchmarks are random 50/50 0-1 bitvectors up to $2^34$ bits. ```bash ./benchmarks ``` +### RmM Tree + +```bash +./bench_rmm +``` + +For visualization, write the CSV output to a file using `--benchmark_out=` (e.g. `./bench_rmm --benchmark_out=rmm_bench.csv`) and plot it with `misc/plot_rmm.py`. + --- ## Example Usage @@ -75,6 +93,30 @@ int main() { } ``` +```cpp +#include "rmm_tree.h" +#include +#include + +using namespace pixie; + +int main() { + // root + // ├─ A + // │ ├─ a1 + // │ └─ a2 + // ├─ B + // └─ C + // └─ c1 + std::string bits = "11101001011000"; + RmMTree t(bits); + + std::cout << "close(1): " << t.close(1) << "\n"; // expected 6 (A) + std::cout << "open(3): " << t.open(3) << "\n"; // expected 2 (a1) + std::cout << "enclose(1): " << t.enclose(1) << "\n"; // expected 0 (root) +} +``` + --- ## References diff --git a/include/misc/naive_rmm_tree.h b/include/misc/naive_rmm_tree.h new file mode 100644 index 0000000..79ecfe9 --- /dev/null +++ b/include/misc/naive_rmm_tree.h @@ -0,0 +1,279 @@ +#pragma once +#include +#include +#include +#include +#include + +class NaiveRmM { + std::vector bits; + std::size_t num_bits = 0; + + public: + static constexpr std::size_t npos = std::numeric_limits::max(); + + NaiveRmM() = default; + explicit NaiveRmM(const std::string& s) { build_from_string(s); } + NaiveRmM(const std::vector& words, std::size_t nbits) { + build_from_words(words, nbits); + } + + private: + void build_from_string(const std::string& s) { + num_bits = s.size(); + bits.resize(num_bits); + for (std::size_t i = 0; i < num_bits; i++) { + bits[i] = (s[i] == '1'); + } + } + void build_from_words(const std::vector& words, + std::size_t nbits) { + num_bits = nbits; + bits.assign(num_bits, 0); + for (std::size_t i = 0; i < num_bits; i++) { + const std::size_t w = (i >> 6); + bits[i] = static_cast( + ((w < words.size() ? words[w] : 0ull) >> (i & 63)) & 1u); + } + } + + inline int bit(std::size_t i) const { return bits[i]; } + + public: + std::size_t rank1(std::size_t i) const { + if (i > num_bits) { + i = num_bits; + } + std::size_t c = 0; + for (std::size_t p = 0; p < i; p++) { + c += (bits[p] != 0); + } + return c; + } + std::size_t rank0(std::size_t i) const { return i - rank1(i); } + + // 1-based + std::size_t select1(std::size_t k) const { + if (k == 0) { + return npos; + } + for (std::size_t p = 0; p < num_bits; p++) { + if (bits[p]) { + if (--k == 0) { + return p; + } + } + } + return npos; + } + + // 1-based + std::size_t select0(std::size_t k) const { + if (k == 0) { + return npos; + } + for (std::size_t p = 0; p < num_bits; p++) { + if (!bits[p]) { + if (--k == 0) { + return p; + } + } + } + return npos; + } + + std::size_t rank10(std::size_t i) const { + if (i <= 1) { + return 0; + } + std::size_t c = 0; + for (std::size_t p = 0; p + 1 < i; ++p) { + if (bits[p] == 1 && bits[p + 1] == 0) { + ++c; + } + } + return c; + } + + // 1-based + std::size_t select10(std::size_t k) const { + if (k == 0) { + return npos; + } + for (std::size_t p = 0; p + 1 < num_bits; ++p) { + if (bits[p] == 1 && bits[p + 1] == 0) { + if (--k == 0) { + return p; + } + } + } + return npos; + } + + int excess(std::size_t i) const { + return int(2 * static_cast(rank1(i)) - + static_cast(i)); + } + + std::size_t fwdsearch(std::size_t i, int d) const { + if (i >= num_bits) { + return npos; + } + int target = excess(i) + d; + int cur = excess(i); + for (std::size_t p = i; p < num_bits; ++p) { + cur += bits[p] ? +1 : -1; + if (cur == target) { + return p; + } + } + return npos; + } + + std::size_t bwdsearch(std::size_t i, int d) const { + if (i > num_bits) { + return npos; + } + if (i == 0) { + return npos; + } + int target = excess(i) + d; + + int cur = excess(i); + for (std::size_t p = i; p > 0;) { + --p; + cur += bits[p] ? -1 : +1; + if (cur == target) { + return p; + } + } + return npos; + } + + std::size_t range_min_query_pos(std::size_t i, std::size_t j) const { + if (i > j || j >= num_bits) { + return npos; + } + int cur = 0, mn = std::numeric_limits::max(); + std::size_t pos = npos; + for (std::size_t p = i; p <= j; ++p) { + cur += bits[p] ? +1 : -1; + if (cur < mn) { + mn = cur; + pos = p; + } + } + return pos; + } + + int range_min_query_val(std::size_t i, std::size_t j) const { + if (i > j || j >= num_bits) { + return 0; + } + int cur = 0, mn = std::numeric_limits::max(); + for (std::size_t p = i; p <= j; ++p) { + cur += bits[p] ? +1 : -1; + if (cur < mn) { + mn = cur; + } + } + return mn; + } + + std::size_t range_max_query_pos(std::size_t i, std::size_t j) const { + if (i > j || j >= num_bits) { + return npos; + } + int cur = 0, mx = std::numeric_limits::min(); + std::size_t pos = npos; + for (std::size_t p = i; p <= j; ++p) { + cur += bits[p] ? +1 : -1; + if (cur > mx) { + mx = cur; + pos = p; + } + } + return pos; + } + + int range_max_query_val(std::size_t i, std::size_t j) const { + if (i > j || j >= num_bits) { + return 0; + } + int cur = 0, mx = std::numeric_limits::min(); + for (std::size_t p = i; p <= j; ++p) { + cur += bits[p] ? +1 : -1; + if (cur > mx) { + mx = cur; + } + } + return mx; + } + + std::size_t mincount(std::size_t i, std::size_t j) const { + if (i > j || j >= num_bits) { + return 0; + } + int cur = 0, mn = std::numeric_limits::max(); + for (std::size_t p = i; p <= j; ++p) { + cur += bits[p] ? +1 : -1; + if (cur < mn) { + mn = cur; + } + } + std::size_t cnt = 0; + cur = 0; + for (std::size_t p = i; p <= j; ++p) { + cur += bits[p] ? +1 : -1; + if (cur == mn) { + ++cnt; + } + } + return cnt; + } + + // (1-based q) + std::size_t minselect(std::size_t i, std::size_t j, std::size_t q) const { + if (i > j || j >= num_bits || q == 0) { + return npos; + } + int cur = 0, mn = std::numeric_limits::max(); + for (std::size_t p = i; p <= j; ++p) { + cur += bits[p] ? +1 : -1; + if (cur < mn) { + mn = cur; + } + } + cur = 0; + for (std::size_t p = i; p <= j; ++p) { + cur += bits[p] ? +1 : -1; + if (cur == mn) { + if (--q == 0) { + return p; + } + } + } + return npos; + } + + std::size_t close(std::size_t i) const { + if (i >= num_bits) { + return npos; + } + return fwdsearch(i, -1); + } + std::size_t open(std::size_t i) const { + if (i == 0 || i > num_bits) { + return npos; + } + auto r = bwdsearch(i, 0); + return (r == npos ? npos : r + 1); + } + std::size_t enclose(std::size_t i) const { + if (i == 0 || i > num_bits) { + return npos; + } + auto r = bwdsearch(i, -2); + return (r == npos ? npos : r + 1); + } +}; diff --git a/include/rmm_tree.h b/include/rmm_tree.h new file mode 100644 index 0000000..263b668 --- /dev/null +++ b/include/rmm_tree.h @@ -0,0 +1,2089 @@ +#pragma once +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace pixie { +/** + * @brief Range min–max tree over a bitvector (LSB-first) tailored for + * balanced-parentheses (BP). + * @details + * Implements the classic rmq/rMq family and BP navigation in O(log n) + * using a perfectly balanced binary tree of blocks. Supported queries: + * - rank1 / rank0 + * - select1 / select0 + * - rank10 / select10 (starts of "10") + * - excess (prefix +1 for '1', −1 for '0') + * - fwdsearch / bwdsearch (prefix–sum search) + * - range_min_query_pos / range_min_query_val (first minimum within a range) + * - range_max_query_pos / range_max_query_val (first maximum within a range) + * - mincount / minselect (count/selection of minima) + * - close / open / enclose (BP navigation) + * + * The bitvector is LSB-first inside each 64-bit word. + */ +class RmMTree { + // ------------ bitvector ------------ + std::vector bits; // LSB-first + size_t num_bits = 0; // number of bits + + // ------------ blocking ------------ + size_t block_bits = 64; // block size (bits), leaf covers <= block_bits bits + size_t leaf_count = 0; // #leaves = ceil(num_bits/block_bits) + + // ------------ tree arrays (heap order: 1 is root) ------------ + // size of segment (in bits) covered by node + // needed for: rank1/rank0, select1/select0, select10, + // excess, fwdsearch/bwdsearch/close/open/enclose, + // range_min_query/range_max_query, minselect. + std::vector segment_size_bits; + + // node_total_excess = total excess (+1 for '1', -1 for '0') on the node + // needed for: rank1/rank0, select1/select0, excess, + // fwdsearch/bwdsearch/close/open/enclose, + // range_min_query/range_max_query, mincount/minselect. + std::vector node_total_excess; + + // node_min_prefix_excess = minimum pref-excess on the node (from 0) + // needed for: fwdsearch/bwdsearch/close/open/enclose, range_min_query, + // mincount/minselect. + std::vector node_min_prefix_excess; + + // node_max_prefix_excess = maximum pref-excess on the node (from 0) + // needed for: fwdsearch/bwdsearch/close/open/enclose, range_max_query. + std::vector node_max_prefix_excess; + + // node_min_count = number of positions where the minimum is attained + // needed for: mincount/minselect. + std::vector node_min_count; + + // node_pattern10_count = # of "10" pattern occurrences inside the node + // needed for: rank10, select10. + std::vector node_pattern10_count; + + // node_first_bit = first bit (0/1), node_last_bit = last bit (0/1) of the + // segment (to handle "10" crossing) both needed for: rank10, select10. + std::vector node_first_bit, node_last_bit; + + public: + /** + * @brief Sentinel for "not found". + */ + static constexpr size_t npos = std::numeric_limits::max(); + +#ifdef DEBUG + float built_overhead = 0.0; +#endif + + // --------- construction ---------- + + /** + * @brief Construct empty structure. + */ + RmMTree() = default; + + /** + * @brief Build from a '0'/'1' string. + * @param bp 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, + 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); + } + + /** + * @brief Build from 64-bit words (LSB-first). + * @param words Array of words holding bits LSB-first. + * @param Nbits 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). + * @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::vector& words, + size_t Nbits, + 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); + } + + // --------- 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. + */ + size_t rank1(const size_t& i) const { + if (i == 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; + } + + /** + * @brief Number of zeros in prefix [0, i). + * @details Computed as i - rank1(i). + */ + size_t rank0(const size_t& i) const { return i - rank1(i); } + + /** + * @brief 1-based select of the k-th one. + * @return Position of k-th '1' or npos if not found. + */ + size_t select1(size_t k) const { + if (k == 0 || num_bits == 0) { + return npos; + } + size_t v = 1; + if (ones_in_node(v) < k) { + 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; + } else { + k -= o_l; + base += segment_size_bits[Lc]; + v = Rc; + } + } + return select1_in_block(base, + std::min(base + segment_size_bits[v], num_bits), k); + } + + /** + * @brief 1-based select of the k-th zero. + * @return Position of k-th '0' or npos if not found. + */ + size_t select0(size_t k) const { + if (k == 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); + }; + if (zeros(v) < k) { + 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; + } else { + k -= z_l; + base += segment_size_bits[Lc]; + v = Rc; + } + } + return select0_in_block(base, + std::min(base + segment_size_bits[v], num_bits), k); + } + + /** + * @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 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; + } + prev_last = node_last_bit[v]; + } + } + const size_t Lb = blk * block_bits; + ans += rr_in_block(Lb, i); + // boundary between the last full node and the leaf tail + if (blk > 0 && i > Lb && prev_last == 1 && bit(Lb) == 0) { + ++ans; + } + return ans; + } + + /** + * @brief 1-based select of the k-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) { + return npos; + } + size_t ind = 1; + if (node_pattern10_count[ind] < k) { + 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; + continue; + } + size_t rem = k - node_pattern10_count[Lc]; + if (cross) { + if (rem == 1) { + return base + segment_size_bits[Lc] - 1; + } + --rem; + } + base += segment_size_bits[Lc]; + ind = Rc; + k = rem; + } + return select10_in_block( + base, std::min(base + segment_size_bits[ind], num_bits), k); + } + + /** + * @brief Prefix excess on [0, i): +1 for '1', −1 for '0'. + */ + inline int excess(const size_t& i) const { + return int64_t(rank1(i)) * 2 - int64_t(i); + } + + /** + * @brief Forward search: first position p ≥ i where excess(p) = excess(i) + + * d. + * @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) { + 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; + } + if (node_min_prefix_excess[v] <= need && + need <= node_max_prefix_excess[v]) { + return descend_fwd(v, need, base); + } + need -= node_total_excess[v]; + base += segment_size_bits[v]; + } + } + return npos; + } + + /** + * @brief Backward search: last position p ≤ i where excess(p) = excess(i) + + * d. + * @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) { + 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 + + // 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; + } + } + // junction between children is a separate branch (allowed only if < i) + if (need_node == node_total_excess[sib] && border < i) { + return border; + } + + // stepped over the sibling, shifted the zero point of the coordinates + need += node_total_excess[sib]; + base -= segment_size_bits[sib]; + } + v >>= 1; + } + return npos; + } + + /** + * @brief Position of the first minimum of excess on [i, j] (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) { + 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; + + int best_val = INT_MAX; + size_t best_pos = npos; + size_t chosen_node = 0; + int pref = 0, pref_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; + } + + // 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; + } + pref += node_total_excess[v]; + } + if ((r & 1) == 0) { + Rnodes[rn++] = r--; + } + 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; + } + pref += node_total_excess[v]; + } + } + + // 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 (best_pos != npos) { + return best_pos; + } + + return descend_first_min(chosen_node, best_val - pref_at_choice, + node_base(chosen_node)); + } + + /** + * @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)). + */ + int range_min_query_val(const size_t& i, const size_t& j) const { + if (i > j || j >= num_bits) { + return 0; + } + size_t p = range_min_query_pos(i, j); + if (p == npos) { + return 0; + } + return excess(p + 1) - excess(i); + } + + /** + * @brief Position of the first maximum of excess on [i, j] (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) { + 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; + + int best_val = INT_MIN; + size_t best_pos = npos; + size_t chosen_node = 0; + int pref = 0, pref_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; + } + + // 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; + } + pref += node_total_excess[v]; + } + if ((r & 1) == 0) { + Rnodes[rn++] = r--; + } + 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; + } + pref += node_total_excess[v]; + } + } + + // 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 (best_pos != npos) { + return best_pos; + } + + return descend_first_max(chosen_node, best_val - pref_at_choice, + node_base(chosen_node)); + } + + /** + * @brief Value of the maximum prefix excess on [i, j] relative to i. + */ + int range_max_query_val(const size_t& i, const size_t& j) const { + if (i > j || j >= num_bits) { + return 0; + } + size_t p = range_max_query_pos(i, j); + if (p == npos) { + return 0; + } + return excess(p + 1) - excess(i); + } + + /** + * @brief How many times the minimum prefix excess occurs on [i, j]. + */ + size_t mincount(const size_t& i, const size_t& j) const { + if (i > j || j >= 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; + + int best_val = INT_MAX; + size_t cnt = 0; + int pref = 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; + } + } + best_val = mn; + cnt = c; + pref = cur; // 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]; + } + pref += node_total_excess[v]; + } + } + + // 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; + } + } + const int cand = pref + mn; + if (cand < best_val) { + best_val = cand; + cnt = c; + } else if (cand == best_val) { + cnt += c; + } + } + return cnt; + } + + /** + * @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. + */ + size_t minselect(const size_t& i, const size_t& j, size_t q) const { + if (i > j || j >= num_bits || q == 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; + + // 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); + } else { + cur_first = 0; + mn_first = INT_MAX; + c_first = 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 + + 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; + + // 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]; + } + pref += node_total_excess[l++]; + } + if ((r & 1) == 0) { + Rnodes[rn++] = r--; + } + 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]; + } + pref += node_total_excess[v]; + } + } + + // 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) { + return npos; + } + + // prefix + if (mn_first == best_val && c_first) { + if (q <= c_first) { + return qth_min_in_block(i, end_first, q); + } + q -= c_first; + } + + // 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)); + } + q -= node_min_count[v]; + } + pref += node_total_excess[v]; + } + if (!(r & 1)) { + Rnodes[rn++] = r--; + } + 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)); + } + q -= node_min_count[v]; + } + pref += node_total_excess[v]; + } + } + + // tail + if (blk_j != blk_i && (pref + mn_last) == best_val) { + return qth_min_in_block(Lbj, j, q); + } + + return npos; + } + + // ----- parentheses navigation (BP) ----- + + /** + * @brief close(i): matching ')' for '(' at i. + * @return Position of matching ')', or npos. + */ + inline size_t close(const size_t& i) const { + if (i >= num_bits) { + return npos; + } + return fwdsearch(i, -1); + } + + /** + * @brief open(i): matching '(' for ')' at i. + * @return Position of matching '(', or npos. + */ + inline size_t open(const size_t& i) const { + // bwdsearch allows i in [1..num_bits] + if (i == 0 || i > num_bits) { + return npos; + } + const size_t r = bwdsearch(i, 0); + return (r == npos ? npos : r + 1); + } + + /** + * @brief enclose(i): opening '(' that strictly encloses position i. + * @return Position of enclosing '(', or npos. + */ + inline size_t enclose(const size_t& i) const { + if (i == 0 || i > num_bits) { + return npos; + } + const size_t r = bwdsearch(i, -2); + return (r == npos ? npos : r + 1); + } + + private: + /** + * @brief Count "10" occurrences inside a 64-bit slice of given logical + * length. + * @details Only positions fully inside the slice are counted. + */ + static inline size_t pop10_in_slice64(const std::uint64_t& slice, + const int& len) noexcept { + if (len <= 1) { + return 0; + } + std::uint64_t P = slice & ~(slice >> 1); // candidates for "10" + if (len < 64) { + P &= ((std::uint64_t(1) << (len - 1)) - 1); + } else { + P &= 0x7FFFFFFFFFFFFFFFull; + } + return (size_t)std::popcount(P); + } + + /** + * @brief Rank of ones within [Lb, Rb). + * @details Works on word boundaries; Rb may equal Lb. + */ + size_t rank1_in_block(const size_t& Lb, const size_t& Rb) const noexcept { + if (Rb <= Lb) { + 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) { + 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); + } + if (off_l) { + cnt += (size_t)std::popcount(bits[w_l] & (~std::uint64_t(0) << off_l)); + ++w_l; + } + while (w_l < w_r) { + cnt += (size_t)std::popcount(bits[w_l]); + ++w_l; + } + if (off_r) { + cnt += + (size_t)std::popcount(bits[w_r] & ((std::uint64_t(1) << off_r) - 1)); + } + return cnt; + } + + /** + * @brief Count "10" starts within [Lb, Rb). + * @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) { + 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; + + 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); + } + + // prefix word + { + const int len = 64 - off_l; + const std::uint64_t slice = bits[w_l] >> off_l; + cnt += pop10_in_slice64(slice, len); + } + // 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); + } + // 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); + } + // 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; + } + } + return cnt; + } + + /** + * @brief 1-based select of k-th "10" within [Lb, Rb). + * @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) { + 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; + + const auto select_in_masked_slice = [&](const std::uint64_t& slice, + const int& len, + const size_t& kk) noexcept -> int { + if (len <= 1) { + return -1; + } + std::uint64_t P = slice & ~(slice >> 1); + if (len < 64) { + P &= ((std::uint64_t(1) << (len - 1)) - 1); + } else { + P &= 0x7FFFFFFFFFFFFFFFull; + } + return select_in_word(P, kk); + }; + + 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; + } + + // 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; + } + k -= c; + } + + // walk interior boundaries and words + for (size_t w = w_l; w + 1 < w_r; ++w) { + // boundary between w and w+1 + if (((bits[w] >> 63) & 1u) && ((bits[w + 1] & 1u) == 0)) { + if (--k == 0) { + return (w << 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) { + return npos; + } + return ((w + 1) << 6) + (size_t)off; + } + k -= c; + } + + // 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; + } + } + + // 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; + } + } + return npos; + } + + struct ByteAgg { + int8_t excess_total; // total excess for the byte + int8_t min_prefix; // minimum prefix within the byte (from 0) + int8_t max_prefix; // maximum prefix within the byte (from 0) + uint8_t min_count; // number of positions attaining the minimum in the byte + uint8_t pattern10_count; // number of "10" patterns inside the byte + uint8_t first_bit; // first bit (LSB) + uint8_t last_bit; // last bit (MSB) + uint8_t pos_first_min; // pos of first minimum in this byte + uint8_t pos_first_max; // pos of first maximum in this byte + }; + + /** + * @brief Returns the static lookup table of byte aggregates. + */ + static inline const std::array& LUT8() noexcept { + static const std::array T = [] { + std::array 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; + }; // LSB-first + for (int k = 0; k < 8; ++k) { + int bit = get(k); + if (k + 1 < 8 && bit && get(k + 1) == 0) { + ++rrc; + } + cur += bit ? +1 : -1; + if (cur < mn) { + mn = cur; + cnt = 1; + pm = k; + } else if (cur == mn) { + ++cnt; + } + if (cur > mx) { + mx = cur; + pM = k; + } + } + 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 T; + } + + /** + * @brief Extract 8 bits starting at position pos (LSB-first across words). + */ + 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); + } + + /** + * @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. + * @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; + } 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; + + first_max_value_pos8(Lb, Rb ? (Rb - 1) : Lb, mx, pos); + return (mx == d ? pos : npos); + } + + /** + * @brief Heap index of the first leaf node. + */ + inline size_t first_leaf_index() const noexcept { + return std::bit_ceil(std::max(1, leaf_count)); + } + + /** + * @brief Block index containing position i. + */ + size_t block_of(const size_t& i) const noexcept { return i / block_bits; } + + /** + * @brief Heap index of the leaf whose segment starts at block_start. + */ + size_t leaf_index_of(const size_t& block_start) const noexcept { + return first_leaf_index() + block_of(block_start); + } + + /** + * @brief Starting global position for node v (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 base = 0; + for (; v > 1; v >>= 1) { + if (v & 1) { + base += segment_size_bits[v - 1]; + } + } + return base; + } + + /** + * @brief Cover a range of whole blocks [a..b] (inclusive) with O(log n) + * maximal nodes. + * @details Returns node indices in left-to-right order. + */ + std::vector 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 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; + } + + /** + * @brief Descend for fwdsearch to find first position where relative prefix + * equals '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; + } else { + need -= node_total_excess[Lc]; + base += segment_size_bits[Lc]; + v = Rc; + } + } + 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; + } + + /** + * @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. + */ + 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]; + + // 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; + } + } + + // 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; + } + + // 3) can we move left within the range? + if (node_min_prefix_excess[Lc] <= need && + need <= node_max_prefix_excess[Lc]) { + v = Lc; + 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; + } + + 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; + } + + /** + * @brief Descend to find first position where node-relative prefix equals d + * (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; + } 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; + + first_min_value_pos8(Lb, Rb ? (Rb - 1) : Lb, mn, pos); + return (mn == d ? pos : npos); + } + + /** + * @brief Descend to find the q-th minimum (1-based) where node-relative + * prefix equals d. + */ + 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; + continue; + } + q -= node_min_count[Lc]; + } + if (rightm == d) { + base += segment_size_bits[Lc]; + d -= node_total_excess[Lc]; + v = Rc; + continue; + } + return npos; + } + return qth_min_in_block( + base, std::min(base + segment_size_bits[v], num_bits) - 1, q); + } + + /** + * @brief 1-based select of k-th '1' within [Lb, Rb). + */ + 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); + } + // 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); + } + k -= c; + w_l++; + } + // 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); + } + k -= c; + ++w_l; + } + // 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); + } + } + return npos; + } + + /** + * @brief 1-based select of k-th '0' within [Lb, Rb). + */ + size_t select0_in_block(const size_t& Lb, + const size_t& Rb, + size_t k) const noexcept { + if (Rb <= Lb) { + return npos; + } + + size_t w_l = Lb >> 6; + const size_t w_r = Rb >> 6; + const size_t off_l = Lb & 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; + } + + // 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; + } + k -= c; + ++w_l; + } + + // 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; + } + k -= c; + ++w_l; + } + + // 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; + } + } + return npos; + } + + /** + * @brief 1-based select of k-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; + } + } + } + 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; + } + + /** + * @brief Number of node slots needed for @p Nbits with leaf size @p Bpow2. + * @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) { + return 0; + } + size_t leaf_count = ceil_div(Nbits, Bpow2); + return std::bit_ceil(std::max(1, leaf_count)) + leaf_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 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); + + size_t bb = ceil_div(Nbits, 64) * 8; + if (bb == 0) { + return 0; + } + size_t slots = nodeslots_for(Nbits, Bpow2); + size_t aux = slots * AUX_SLOT_BYTES; + return ((float)aux) / ((float)bb); + } + + /** + * @brief Choose minimal block size (power of two) keeping overhead ≤ cap. + * @details Returns 64 if cap < 0 (no constraint). Clamped to ≤16384 or Nbits. + */ + static inline size_t choose_block_bits_for_overhead( + const size_t& Nbits, + const float& cap) noexcept { + if (cap < 0.f) { + return 64; + } + + const size_t Bmax = std::min(Nbits, 16384); + size_t block_bits = 64; + while (block_bits < Bmax) { + if (overhead_for(Nbits, block_bits) <= cap) { + break; + } + block_bits <<= 1; + } + return block_bits; + } + + /** + * @brief Build internal structures from a 0/1 string. + * @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, + const size_t& leaf_block_bits = 0, + const float& max_overhead = -1.0) { + num_bits = s.size(); + bits.assign((num_bits + 63) / 64, 0); + for (size_t i = 0; i < num_bits; ++i) { + if (s[i] == '1') { + set1(i); + } + } + build(leaf_block_bits, max_overhead); + } + + /** + * @brief Build internal structures from 64-bit words. + * @param words Words with LSB-first bits. + * @param Nbits 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& words, + const size_t& Nbits, + const size_t& leaf_block_bits = 0, + const float& max_overhead = -1.0) { + bits = words; + num_bits = Nbits; + if (bits.size() * 64 < num_bits) { + bits.resize((num_bits + 63) / 64); + } + build(leaf_block_bits, max_overhead); + } + + /** + * @brief Read bit at position i (LSB-first across words). + */ + inline int bit(const size_t& i) const noexcept { + return (bits[i >> 6] >> (i & 63)) & 1u; + } + + /** + * @brief Set bit at position i to 1. + */ + inline void set1(const size_t& i) noexcept { + bits[i >> 6] |= (std::uint64_t(1) << (i & 63)); + } + + /** + * @brief Number of ones in node v 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; + } + + /** + * @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. + */ + 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; + 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; + } + ++l; + } + // 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; + } + // tail + while (l <= r) { + cur += bit(l) ? +1 : -1; + if (cur < mn) { + mn = cur; + cnt = 1; + } else if (cur == mn) { + ++cnt; + } + ++l; + } + if (mn == INT_MAX) { + mn = cnt = 0; + } + } + + /** + * @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. + */ + 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) { + return npos; + } + + const auto& T = LUT8(); + + int cur = 0, mn = INT_MAX; + size_t p = l; + + while (p <= r && (p & 7)) { + cur += bit(p) ? +1 : -1; + if (cur < mn) { + mn = cur; + } + ++p; + } + 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 (p <= r) { + cur += bit(p) ? +1 : -1; + if (cur < mn) { + mn = cur; + } + ++p; + } + + cur = 0; + p = l; + + // to byte alignment + while (p <= r && (p & 7)) { + cur += bit(p) ? +1 : -1; + if (cur == mn) { + if (--q == 0) { + return p; + } + } + ++p; + } + + // 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; + for (int k = 0; k < 8; ++k) { + s += ((b >> k) & 1u) ? +1 : -1; + if (s == a.min_prefix) { + if (--q == 0) { + return p + k; + } + } + } + } + cur += a.excess_total; + p += 8; + } + + // tail + while (p <= r) { + cur += bit(p) ? +1 : -1; + if (cur == mn) { + if (--q == 0) { + return p; + } + } + ++p; + } + + 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). + */ + 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; + + // to byte allignment + while (l <= r && (l & 7)) { + cur += bit(l) ? +1 : -1; + if (cur < mn) { + mn = cur; + first_pos = l; + } + ++l; + } + + // 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; + } + cur += a.excess_total; + l += 8; + } + + // tail + while (l <= r) { + cur += bit(l) ? +1 : -1; + if (cur < mn) { + mn = cur; + first_pos = l; + } + ++l; + } + + mn_out = (mn == INT_MAX ? 0 : mn); + } + + /** + * @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). + */ + 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; + + while (l <= r && (l & 7)) { + cur += bit(l) ? +1 : -1; + if (cur > mx) { + mx = cur; + first_pos = l; + } + ++l; + } + + 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; + } + cur += a.excess_total; + l += 8; + } + + while (l <= r) { + cur += bit(l) ? +1 : -1; + if (cur > mx) { + mx = cur; + first_pos = l; + } + ++l; + } + + mx_out = (mx == INT_MIN ? 0 : mx); + } + + /** + * @brief Build the tree arrays and per-node aggregates. + * @details Chooses block_bits honoring @p max_overhead or explicit @p + * leaf_block_bits, allocates arrays, fills leaves via LUT, and builds + * internal nodes bottom-up. + */ + void build(const size_t& leaf_block_bits, const float& max_overhead) { + // the lower clamp depends on the desired overhead fraction; otherwise use + // 64 + const size_t clamp_by_overhead = + (max_overhead >= 0.0 + ? choose_block_bits_for_overhead(num_bits, max_overhead) + : size_t(64)); + + // chosen block_bits: honor an explicit request, but not below + // clamp_by_overhead + if (leaf_block_bits == 0) { + block_bits = + std::max(clamp_by_overhead, + std::bit_ceil( + (num_bits <= 1) ? 1 : std::bit_width(num_bits - 1))); + } else { + block_bits = + std::max(clamp_by_overhead, + std::bit_ceil(std::max(1, leaf_block_bits))); + } + +#ifdef DEBUG + // finalizes the achieved overhead percentage + built_overhead = overhead_for(num_bits, block_bits); +#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; + 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); + node_max_prefix_excess.assign(tree_size + 1, 0); + node_min_count.assign(tree_size + 1, 0); + node_pattern10_count.assign(tree_size + 1, 0); + node_first_bit.assign(tree_size + 1, 0); + node_last_bit.assign(tree_size + 1, 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; + + if (Lb < Rb) { + node_first_bit[v] = bit(Lb); + } + + const auto& T = LUT8(); + + int cur = 0, mn = INT_MAX, mx = INT_MIN; + uint32_t mn_cnt = 0; + uint32_t rrc = 0; + + uint8_t prev_bit = 0; + + size_t p = Lb; + + // Full bytes + while (p + 8 <= Rb) { + const uint8_t b = get_byte(p); + const auto& a = T[b]; + + // internal "10" inside the byte + rrc += a.pattern10_count; + // stitching across the boundary between the previous and current byte + // (within the segment) + if (prev_bit == 1 && a.first_bit == 0) { + rrc++; + } + + // 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; + } + + mx = std::max(mx, cur + a.max_prefix); + cur += a.excess_total; + prev_bit = a.last_bit; + p += 8; + } + + // Tail < 8 bits + while (p < Rb) { + const uint8_t b = bit(p); + if (prev_bit == 1 && b == 0) { + rrc++; + } + const int step = b ? +1 : -1; + cur += step; + if (cur < mn) { + mn = cur; + mn_cnt = 1; + } else if (cur == mn) { + ++mn_cnt; + } + if (cur > mx) { + mx = cur; + } + + prev_bit = b; + ++p; + } + + if (Lb < Rb) { + node_last_bit[v] = prev_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; + } + // 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; + 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]; + } 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] + : 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) { + break; + } + } + } +}; + +} // namespace pixie diff --git a/misc/plot_rmm.py b/misc/plot_rmm.py new file mode 100644 index 0000000..881d50c --- /dev/null +++ b/misc/plot_rmm.py @@ -0,0 +1,135 @@ +""" +Plot RmMTree benchmark results. + +This script reads a CSV produced by `bench_rmm.cpp` 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 +""" + +import argparse +import os +import pandas as pd +import matplotlib.pyplot as plt + +OPS_ORDER = [ + "rank1", + "rank0", + "select1", + "select0", + "rank10", + "select10", + "excess", + "fwdsearch", + "bwdsearch", + "range_min_query_pos", + "range_min_query_val", + "range_max_query_pos", + "range_max_query_val", + "mincount", + "minselect", + "close", + "open", + "enclose", +] + + +def main(): + ap = argparse.ArgumentParser( + description=( + "Read a CSV 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" + ), + formatter_class=argparse.RawDescriptionHelpFormatter, + ) + ap.add_argument( + "csv", + metavar="CSV", + help="Path to the CSV file with results (output of bench_rmm.cpp).", + ) + ap.add_argument( + "--save-dir", + default="plots", + metavar="DIR", + help=( + "Directory to save PNG plots. Will be created if it doesn't exist. " + "Default: %(default)s" + ), + ) + ap.add_argument( + "--show", + action="store_true", + help=( + "Show plot windows after saving. " + "By default, plots are only written to disk." + ), + ) + ap.add_argument( + "--logx", + action="store_true", + help=( + "Use a logarithmic X axis (base 2). " "Handy when N grows in powers of two." + ), + ) + ap.add_argument( + "--smooth", + type=int, + default=0, + metavar="W", + help=( + "Median smoothing window size for the trend line. " + "0 or 1 means no smoothing. Default: %(default)s" + ), + ) + args = ap.parse_args() + + os.makedirs(args.save_dir, exist_ok=True) + df = pd.read_csv(args.csv) + df = df.dropna(subset=["cpu_time", "N"]) + + for op in OPS_ORDER: + d = df[df["name"] == op].copy() + if d.empty: + continue + d = d.groupby("N", as_index=False)["cpu_time"].median().sort_values("N") + + yplot = d["cpu_time"] + if args.smooth and args.smooth > 1: + d["ns_smooth"] = ( + d["cpu_time"] + .rolling(window=args.smooth, center=True, min_periods=1) + .median() + ) + yplot = d["ns_smooth"] + + plt.figure() + plt.scatter(d["N"], d["cpu_time"], s=8, alpha=0.3, linewidths=0) + plt.plot(d["N"], yplot, linewidth=1.5) + + if args.logx: + plt.xscale("log", base=2) + plt.xlabel("Sequence size N (bits)") + plt.ylabel("Time per operation, ns") + plt.title(f"RmMTree: {op}") + plt.grid(True, which="both", linestyle="--", alpha=0.4) + out = os.path.join(args.save_dir, f"{op}.png") + plt.savefig(out, bbox_inches="tight", dpi=160) + if not args.show: + plt.close() + print(f"[saved] {out}") + + if args.show: + plt.show() + + +if __name__ == "__main__": + main() diff --git a/src/bench_rmm.cpp b/src/bench_rmm.cpp new file mode 100644 index 0000000..06db0db --- /dev/null +++ b/src/bench_rmm.cpp @@ -0,0 +1,462 @@ +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "rmm_tree.h" + +using namespace std; +using pixie::RmMTree; + +struct Args { + int min_exp = 14; + int max_exp = 22; + size_t Q = 200000; + double p1 = 0.5; + uint64_t seed = 42; + size_t block_bits = 64; + int per_octave = 10; + vector explicit_sizes; +}; + +static Args args; + +static void parse_args_and_strip(int& argc, char**& argv) { + auto getv = [&](const string& key) -> optional { + string pref = "--" + key + "="; + for (int i = 1; i < argc; ++i) { + string s = argv[i]; + if (s.rfind(pref, 0) == 0) { + return s.substr(pref.size()); + } + } + return nullopt; + }; + auto strip = [&](const string& key) { + string pref = "--" + key + "="; + int w = 0; + for (int i = 0; i < argc; ++i) { + if (i > 0 && string(argv[i]).rfind(pref, 0) == 0) { + continue; + } + argv[w++] = argv[i]; + } + argc = w; + }; + + if (auto v = getv("min-exp")) { + args.min_exp = stoi(*v), strip("min-exp"); + } + if (auto v = getv("max-exp")) { + args.max_exp = stoi(*v), strip("max-exp"); + } + if (auto v = getv("q")) { + args.Q = stoull(*v), strip("q"); + } + if (auto v = getv("p")) { + args.p1 = stod(*v), strip("p"); + } + if (auto v = getv("seed")) { + args.seed = stoull(*v), strip("seed"); + } + if (auto v = getv("block")) { + args.block_bits = stoull(*v), strip("block"); + } + if (auto v = getv("per-octave")) { + args.per_octave = stoi(*v), strip("per-octave"); + } + + if (auto v = getv("sizes")) { + string s = *v; + strip("sizes"); + size_t pos = 0; + while (pos < s.size()) { + while (pos < s.size() && + (s[pos] == ',' || isspace((unsigned char)s[pos]))) { + ++pos; + } + size_t start = pos; + while (pos < s.size() && isdigit((unsigned char)s[pos])) { + ++pos; + } + if (start < pos) { + args.explicit_sizes.push_back(stoull(s.substr(start, pos - start))); + } + } + sort(args.explicit_sizes.begin(), args.explicit_sizes.end()); + args.explicit_sizes.erase( + 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 U(0.0, 1.0); + string s; + s.resize(N); + for (size_t i = 0; i < N; ++i) { + s[i] = (U(rng) < p1 ? '1' : '0'); + } + return s; +} + +static vector build_size_grid() { + if (!args.explicit_sizes.empty()) { + vector v = args.explicit_sizes; + v.erase(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 v; + for (int e = lo; e <= hi; ++e) { + v.push_back(size_t(1) << e); + } + return v; + } + + int per_octave_steps = args.per_octave; + set N_cands; + size_t min_N = size_t(1) << lo, max_N = 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)); + if (N < min_N) { + N = min_N; + } + if (N > max_N) { + N = max_N; + } + if (N) { + N_cands.insert(N); + } + } + } + + vector v(N_cands.begin(), N_cands.end()); + if (v.front() != min_N) { + v.insert(v.begin(), min_N); + } + if (v.back() != max_N) { + v.push_back(max_N); + } + return v; +} + +struct Pools { + vector inds_any; + vector inds; + vector inds_1N; + vector deltas; + vector> segs; + + vector ks1, ks0, ks10; + vector minselect_q; +}; + +struct Dataset { + size_t N{}; + string bits; + RmMTree t; + + size_t cnt1{}, cnt0{}, cnt10{}; + Pools pool; +}; + +static vector> keepalive; + +static size_t count10(const string& s) { + size_t c = 0; + if (s.size() < 2) { + return 0; + } + for (size_t i = 0; i + 1 < s.size(); ++i) { + if (s[i] == '1' && s[i + 1] == '0') { + ++c; + } + } + return c; +} + +static Dataset build_dataset(size_t N) { + mt19937_64 rng(args.seed ^ (uint64_t)N * 0x9E3779B185EBCA87ull); + + Dataset d; + d.N = N; + d.bits = make_random_bits(N, args.p1, rng); + + if (args.block_bits == 0) { + d.t = RmMTree(d.bits, 0, -1.0f); + } else { + d.t = RmMTree(d.bits, args.block_bits); + } + + d.cnt1 = d.t.rank1(N); + d.cnt0 = N - d.cnt1; + d.cnt10 = count10(d.bits); + + const size_t L = min(args.Q, 32768); + + uniform_int_distribution ind_dist_incl(0, N ? N : 0); + uniform_int_distribution ind_dist(0, N ? (N - 1) : 0); + uniform_int_distribution ind_dist_1N(1, N ? N : 1); + uniform_int_distribution d_dist(-8, +8); + + d.pool.inds_any.resize(L); + d.pool.inds.resize(L); + d.pool.inds_1N.resize(L); + d.pool.deltas.resize(L); + d.pool.segs.resize(L); + + auto rand_ij = [&]() -> pair { + if (N == 0) { + return {0, 0}; + } + size_t a = ind_dist(rng), b = ind_dist(rng); + if (a > b) { + swap(a, b); + } + return {a, b}; + }; + + for (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); + d.pool.deltas[i] = d_dist(rng); + d.pool.segs[i] = rand_ij(); + } + + auto fill_ks = [&](size_t total, vector& out) { + out.resize(L); + if (total == 0) { + fill(out.begin(), out.end(), 1); + return; + } + uniform_int_distribution dist(1, total); + for (size_t i = 0; i < L; ++i) { + out[i] = dist(rng); + } + }; + fill_ks(d.cnt1, d.pool.ks1); + fill_ks(d.cnt0, d.pool.ks0); + fill_ks(d.cnt10, d.pool.ks10); + + d.pool.minselect_q.resize(L); + for (size_t i = 0; i < L; ++i) { + auto [l, r] = d.pool.segs[i]; + size_t c = d.t.mincount(l, r); + if (c == 0) { + c = 1; + } + uniform_int_distribution uq(1, c); + d.pool.minselect_q[i] = uq(rng); + } + + return d; +} + +template +static void register_op(const string& op, shared_ptr data, Fn&& body) { + auto idx_ptr = make_shared(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)++; + body(state, D, i); + } + state.counters["N"] = static_cast(D.N); + state.counters["seed"] = static_cast(args.seed); + state.counters["block_bits"] = static_cast(args.block_bits); + }); + + b->Unit(benchmark::kNanosecond); +} + +static void register_all() { + auto Ns = build_size_grid(); + for (size_t N : Ns) { + auto data = make_shared(build_dataset(N)); + keepalive.push_back(data); + + const auto& P = data->pool; + + register_op("rank1", data, + [&](benchmark::State&, const Dataset& D, 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) { + 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) { + 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) { + 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) { + 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) { + 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) { + 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) { + 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) { + 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) { + 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) { + 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) { + 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) { + 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) { + 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(); + auto [i, j] = P.segs[idx]; + auto q = P.minselect_q[idx]; + auto r = D.t.minselect(i, j, q); + benchmark::DoNotOptimize(r); + }); + + register_op("close", data, + [&](benchmark::State&, const Dataset& D, size_t k) { + if (D.N == 0) { + benchmark::DoNotOptimize(0); + return; + } + auto i = P.inds[k % P.inds.size()]; + auto r = D.t.close(i); + benchmark::DoNotOptimize(r); + }); + + register_op("open", data, + [&](benchmark::State&, const Dataset& D, 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) { + auto r = D.t.enclose(P.inds_1N[k % P.inds_1N.size()]); + benchmark::DoNotOptimize(r); + }); + } +} + +int main(int argc, char** argv) { + parse_args_and_strip(argc, argv); + + auto has = [&](const char* key) { + string p1 = string(key) + "="; + for (int i = 1; i < argc; ++i) { + string s = argv[i]; + if (s == key || s.rfind(p1, 0) == 0) { + return true; + } + } + return false; + }; + + static vector extra; + if (!has("--benchmark_out_format")) { + extra.emplace_back("--benchmark_out_format=csv"); + } + if (!has("--benchmark_counters_tabular")) { + extra.emplace_back("--benchmark_counters_tabular=true"); + } + if (!has("--benchmark_time_unit")) { + extra.emplace_back("--benchmark_time_unit=ns"); + } + + static vector argv_vec; + argv_vec.assign(argv, argv + argc); + for (auto& s : extra) { + argv_vec.push_back(s.data()); + } + argv_vec.push_back(nullptr); + argc = (int)argv_vec.size() - 1; + argv = argv_vec.data(); + + benchmark::Initialize(&argc, argv); + register_all(); + benchmark::RunSpecifiedBenchmarks(); + benchmark::Shutdown(); + return 0; +} diff --git a/src/test_rmm.cpp b/src/test_rmm.cpp new file mode 100644 index 0000000..9571a91 --- /dev/null +++ b/src/test_rmm.cpp @@ -0,0 +1,643 @@ +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "naive_rmm_tree.h" +#include "rmm_tree.h" + +using std::size_t; + +static std::string bits_to_parens(const std::string& bits) { + std::string s; + s.reserve(bits.size()); + for (char c : bits) { + s.push_back(c == '1' ? '(' : ')'); + } + return s; +} + +static std::string vecbits_to_string(const std::vector& v) { + std::string s; + s.resize(v.size()); + for (size_t i = 0; i < v.size(); ++i) { + s[i] = v[i] ? '1' : '0'; + } + return s; +} + +static std::vector pack_words_lsb_first( + const std::string& bits) { + const size_t n = bits.size(); + std::vector w((n + 63) / 64, 0); + for (size_t i = 0; i < n; ++i) { + if (bits[i] == '1') { + w[i >> 6] |= (std::uint64_t(1) << (i & 63)); + } + } + return w; +} + +static std::string random_bits(std::mt19937_64& rng, size_t n) { + std::uniform_int_distribution b01(0, 1); + std::string s; + s.resize(n); + for (size_t i = 0; i < n; i++) { + s[i] = char('0' + b01(rng)); + } + return s; +} + +static std::string random_dyck_bits(std::mt19937_64& rng, size_t m) { + std::string s(2 * m, '0'); + if (m == 0) { + return s; + } + size_t opens_left = m; + size_t closes_left = m; + int h = 0; + std::bernoulli_distribution coin(0.5); + for (size_t pos = 0; pos < 2 * m; ++pos) { + if (opens_left == 0) { + s[pos] = '0'; + --closes_left; + --h; + continue; + } + if (closes_left == 0) { + s[pos] = '1'; + --opens_left; + ++h; + continue; + } + if (h == 0) { + s[pos] = '1'; + --opens_left; + ++h; + continue; + } + if ((size_t)h == closes_left) { + s[pos] = '0'; + --closes_left; + --h; + continue; + } + if (coin(rng)) { + s[pos] = '1'; + --opens_left; + ++h; + } else { + s[pos] = '0'; + --closes_left; + --h; + } + } + return s; +} + +struct Limits { + size_t CASES = 150; + size_t OPS_PER_CASE = 600; + size_t MAX_N = 20000; +}; + +static Limits load_limits_from_env() { + Limits L; + if (const char* s = std::getenv("RMM_CASES")) { + L.CASES = std::max(1, std::strtoull(s, nullptr, 10)); + } + if (const char* s = std::getenv("RMM_OPS")) { + L.OPS_PER_CASE = std::max(1, std::strtoull(s, nullptr, 10)); + } + if (const char* s = std::getenv("RMM_MAX_N")) { + L.MAX_N = std::max(1, std::strtoull(s, nullptr, 10)); + } + return L; +} + +static uint64_t choose_seed() { + if (const char* s = std::getenv("RMM_SEED")) { + return std::strtoull(s, nullptr, 10); + } + std::random_device rd; + return ((uint64_t)rd() << 32) ^ (uint64_t)rd() ^ + (uint64_t)std::chrono::high_resolution_clock::now() + .time_since_epoch() + .count(); +} + +static void run_case_and_compare(const std::string& bits, + std::mt19937_64& rng, + size_t ops_per_case, + uint64_t seed) { + const bool use_words = std::uniform_int_distribution(0, 1)(rng); + pixie::RmMTree rm; + NaiveRmM nv; + if (use_words) { + auto words = pack_words_lsb_first(bits); + rm = pixie::RmMTree(words, bits.size()); + nv = NaiveRmM(words, bits.size()); + } else { + rm = pixie::RmMTree(bits); + nv = NaiveRmM(bits); + } + + const size_t N = bits.size(); + const size_t ones = nv.rank1(N); + const size_t zeros = N - ones; + const size_t pairs10 = (N >= 2 ? nv.rank10(N) : 0); + + std::uniform_int_distribution pos_i(0, N); + std::uniform_int_distribution pos_i_nz(0, N ? N - 1 : 0); + std::uniform_int_distribution d_dist(-(int)std::min(N, 200), + (int)std::min(N, 200)); + + for (size_t q = 0; q < ops_per_case; ++q) { + int which = std::uniform_int_distribution(0, 17)(rng); + + size_t i = 0, j = 0; + if (N > 0) { + i = pos_i_nz(rng); + j = pos_i_nz(rng); + if (i > j) { + std::swap(i, j); + } + } + + auto trace_header = ::testing::Message() + << "\nseed=" << seed << " | N=" << N + << " | use_words=" << (use_words ? 1 : 0) + << "\nbits: " << bits + << "\nparens: " << bits_to_parens(bits) + << "\nwhich=" << which << " (0..17)\n"; + + switch (which) { + case 0: { // rank1 + size_t x = pos_i(rng); + auto a = nv.rank1(x), b = rm.rank1(x); + EXPECT_EQ(a, b) << trace_header << "op=rank1(" << x << ")"; + } break; + + case 1: { // rank0 + size_t x = pos_i(rng); + auto a = nv.rank0(x), b = rm.rank0(x); + EXPECT_EQ(a, b) << trace_header << "op=rank0(" << x << ")"; + } break; + + case 2: { // select1 + size_t k = std::uniform_int_distribution(0, ones + 3)(rng); + auto a = nv.select1(k), b = rm.select1(k); + EXPECT_EQ(a, b) << trace_header << "op=select1(" << k << ")"; + } break; + + case 3: { // select0 + size_t k = std::uniform_int_distribution(0, zeros + 3)(rng); + auto a = nv.select0(k), b = rm.select0(k); + EXPECT_EQ(a, b) << trace_header << "op=select0(" << k << ")"; + } break; + + case 4: { // rank10 + size_t x = + (N >= 2 ? std::uniform_int_distribution(0, N)(rng) : 0); + auto a = nv.rank10(x), b = rm.rank10(x); + EXPECT_EQ(a, b) << trace_header << "op=rank10(" << x << ")"; + } break; + + case 5: { // select10 + size_t k = std::uniform_int_distribution(0, pairs10 + 3)(rng); + auto a = nv.select10(k), b = rm.select10(k); + EXPECT_EQ(a, b) << trace_header << "op=select10(" << k << ")"; + } break; + + case 6: { // excess + size_t x = pos_i(rng); + auto a = nv.excess(x), b = rm.excess(x); + EXPECT_EQ(a, b) << trace_header << "op=excess(" << x << ")"; + } break; + + case 7: { // fwdsearch + if (N == 0) { + break; + } + size_t start = pos_i_nz(rng); + int d = d_dist(rng); + auto a = nv.fwdsearch(start, d), b = rm.fwdsearch(start, d); + EXPECT_EQ(a, b) << trace_header << "op=fwdsearch(" << start << "," << d + << ")"; + } break; + + case 8: { // bwdsearch + if (N == 0) { + break; + } + size_t start = pos_i_nz(rng); + int d = d_dist(rng); + auto a = nv.bwdsearch(start, d), b = rm.bwdsearch(start, d); + EXPECT_EQ(a, b) << trace_header << "op=bwdsearch(" << start << "," << d + << ")"; + } break; + + case 9: { // range_min_query_pos + if (N == 0) { + break; + } + auto a = nv.range_min_query_pos(i, j), b = rm.range_min_query_pos(i, j); + EXPECT_EQ(a, b) << trace_header << "op=range_min_query_pos(" << i << "," + << j << ")"; + } break; + + case 10: { // range_min_query_val + if (N == 0) { + break; + } + auto a = nv.range_min_query_val(i, j), b = rm.range_min_query_val(i, j); + EXPECT_EQ(a, b) << trace_header << "op=range_min_query_val(" << i << "," + << j << ")"; + } break; + + case 11: { // mincount + if (N == 0) { + break; + } + auto a = nv.mincount(i, j), b = rm.mincount(i, j); + EXPECT_EQ(a, b) << trace_header << "op=mincount(" << i << "," << j + << ")"; + } break; + + case 12: { // minselect + if (N == 0) { + break; + } + size_t cnt = nv.mincount(i, j); + size_t k = cnt == 0 + ? 1 + : std::uniform_int_distribution(1, cnt + 1)(rng); + auto a = nv.minselect(i, j, k), b = rm.minselect(i, j, k); + EXPECT_EQ(a, b) << trace_header << "op=minselect(" << i << "," << j + << "," << k << ")"; + } break; + + case 13: { // range_max_query_pos + if (N == 0) { + break; + } + auto a = nv.range_max_query_pos(i, j), b = rm.range_max_query_pos(i, j); + EXPECT_EQ(a, b) << trace_header << "op=range_max_query_pos(" << i << "," + << j << ")"; + } break; + + case 14: { // range_max_query_val + if (N == 0) { + break; + } + auto a = nv.range_max_query_val(i, j), b = rm.range_max_query_val(i, j); + EXPECT_EQ(a, b) << trace_header << "op=range_max_query_val(" << i << "," + << j << ")"; + } break; + + case 15: { // close + if (N == 0) { + break; + } + size_t x = pos_i_nz(rng); + auto a = nv.close(x), b = rm.close(x); + EXPECT_EQ(a, b) << trace_header << "op=close(" << x << ")"; + } break; + + case 16: { // open + if (N == 0) { + break; + } + size_t x = pos_i(rng); + auto a = nv.open(x), b = rm.open(x); + EXPECT_EQ(a, b) << trace_header << "op=open(" << x << ")"; + } break; + + case 17: { // enclose + if (N == 0) { + break; + } + size_t x = pos_i(rng); + auto a = nv.enclose(x), b = rm.enclose(x); + EXPECT_EQ(a, b) << trace_header << "op=enclose(" << x << ")"; + } break; + } + } +} + +class RmMRandomTest : public ::testing::Test { + protected: + uint64_t seed{}; + std::mt19937_64 rng; + Limits L; + void SetUp() override { + L = load_limits_from_env(); + seed = choose_seed(); + rng.seed(seed); + std::cerr << "[ RmMRandomTest ] seed=" << seed << " CASES=" << L.CASES + << " OPS=" << L.OPS_PER_CASE << " MAX_N=" << L.MAX_N << "\n"; + } +}; + +TEST_F(RmMRandomTest, RandomMixedBits) { + std::uniform_int_distribution len_u(1, (int)L.MAX_N); + for (size_t t = 0; t < L.CASES; ++t) { + const size_t n = (size_t)len_u(rng); + const std::string bits = random_bits(rng, n); + run_case_and_compare(bits, rng, L.OPS_PER_CASE, seed); + } +} + +TEST_F(RmMRandomTest, RandomDyckBits) { + std::uniform_int_distribution len_even(0, (int)(L.MAX_N / 2)); + for (size_t t = 0; t < L.CASES; ++t) { + size_t m = 1 + (size_t)len_even(rng); + std::string bits = random_dyck_bits(rng, m); + if (bits.empty()) { + bits = "10"; + } + run_case_and_compare(bits, rng, L.OPS_PER_CASE, seed); + } +} + +TEST_F(RmMRandomTest, ShortInputs) { + for (size_t n = 1; n <= 8; ++n) { + const size_t total = (n ? (1ull << n) : 1ull); + for (size_t mask = 0; mask < total; ++mask) { + std::string bits; + bits.resize(n); + for (size_t i = 0; i < n; ++i) { + bits[i] = ((mask >> i) & 1ull) ? '1' : '0'; + } + { + SCOPED_TRACE(::testing::Message() + << "short-inputs:string bits=" << bits); + run_case_and_compare(bits, rng, /*ops_per_case=*/200, + /*seed=*/0xC0FFEEull); + } + { + SCOPED_TRACE(::testing::Message() + << "short-inputs:words bits=" << bits); + run_case_and_compare(bits, rng, /*ops_per_case=*/200, + /*seed=*/0xC0FFEEull); + } + } + } +} + +static void run_case_string_vs_words(const std::string& bits, + std::mt19937_64& rng) { + const size_t N = bits.size(); + pixie::RmMTree rm_s(bits); + auto words = pack_words_lsb_first(bits); + pixie::RmMTree rm_w(words, N); + + std::uniform_int_distribution pos_i(0, N); + std::uniform_int_distribution pos_i_nz(0, N ? N - 1 : 0); + std::uniform_int_distribution d_dist(-(int)std::min(N, 200), + (int)std::min(N, 200)); + + for (int t = 0; t < 128; ++t) { + int which = std::uniform_int_distribution(0, 10)(rng); + size_t i = 0, j = 0; + if (N > 0) { + i = pos_i_nz(rng); + j = pos_i_nz(rng); + if (i > j) { + std::swap(i, j); + } + } + switch (which) { + case 0: { + size_t x = pos_i(rng); + EXPECT_EQ(rm_s.rank1(x), rm_w.rank1(x)); + break; + } + case 1: { + size_t x = pos_i(rng); + EXPECT_EQ(rm_s.rank0(x), rm_w.rank0(x)); + break; + } + case 2: { + size_t k = std::uniform_int_distribution(0, N + 3)(rng); + EXPECT_EQ(rm_s.select1(k), rm_w.select1(k)); + break; + } + case 3: { + size_t k = std::uniform_int_distribution(0, N + 3)(rng); + EXPECT_EQ(rm_s.select0(k), rm_w.select0(k)); + break; + } + case 4: { + size_t x = + (N >= 2 ? std::uniform_int_distribution(0, N)(rng) : 0); + EXPECT_EQ(rm_s.rank10(x), rm_w.rank10(x)); + break; + } + case 5: { + size_t k = std::uniform_int_distribution(0, N)(rng); + EXPECT_EQ(rm_s.select10(k), rm_w.select10(k)); + break; + } + case 6: { + size_t x = pos_i(rng); + EXPECT_EQ(rm_s.excess(x), rm_w.excess(x)); + break; + } + case 7: { + if (N == 0) { + break; + } + size_t start = pos_i_nz(rng); + int d = d_dist(rng); + EXPECT_EQ(rm_s.fwdsearch(start, d), rm_w.fwdsearch(start, d)); + break; + } + case 8: { + if (N == 0) { + break; + } + size_t start = pos_i_nz(rng); + int d = d_dist(rng); + EXPECT_EQ(rm_s.bwdsearch(start, d), rm_w.bwdsearch(start, d)); + break; + } + case 9: { + if (N == 0) { + break; + } + EXPECT_EQ(rm_s.range_min_query_pos(i, j), + rm_w.range_min_query_pos(i, j)); + break; + } + case 10: { + if (N == 0) { + break; + } + EXPECT_EQ(rm_s.range_max_query_pos(i, j), + rm_w.range_max_query_pos(i, j)); + break; + } + } + } +} + +TEST_F(RmMRandomTest, WordsVsString) { + std::uniform_int_distribution len_u(0, (int)L.MAX_N); + for (size_t t = 0; t < L.CASES; ++t) { + const size_t n = (size_t)len_u(rng); + const std::string bits = random_bits(rng, n); + run_case_string_vs_words(bits, rng); + } +} + +TEST(RmMEdgeCases, EmptyInput) { + pixie::RmMTree rm(std::string{}); + NaiveRmM nv(std::string{}); + EXPECT_EQ(rm.rank1(0), nv.rank1(0)); + EXPECT_EQ(rm.rank0(0), nv.rank0(0)); + EXPECT_EQ(rm.rank10(0), nv.rank10(0)); + EXPECT_EQ(rm.select1(1), nv.select1(1)); + EXPECT_EQ(rm.select0(1), nv.select0(1)); + EXPECT_EQ(rm.fwdsearch(0, 0), nv.fwdsearch(0, 0)); + EXPECT_EQ(rm.bwdsearch(0, 0), nv.bwdsearch(0, 0)); + EXPECT_EQ(rm.range_min_query_pos(0, 0), nv.range_min_query_pos(0, 0)); + EXPECT_EQ(rm.range_max_query_pos(0, 0), nv.range_max_query_pos(0, 0)); +} + +TEST(RmMTreeStress, LongRandom) { + Limits L; + L.OPS_PER_CASE = 2000; + L.MAX_N = 65536; + + if (const char* s = std::getenv("RMM_CASES")) { + L.CASES = std::max(1, std::strtoull(s, nullptr, 10)); + } + if (const char* s = std::getenv("RMM_OPS")) { + L.OPS_PER_CASE = std::max(1, std::strtoull(s, nullptr, 10)); + } + if (const char* s = std::getenv("RMM_MAX_N")) { + L.MAX_N = std::max(1, std::strtoull(s, nullptr, 10)); + } + + const uint64_t seed = choose_seed(); + std::mt19937_64 rng(seed); + size_t LOG_EVERY = 10; + if (const char* s = std::getenv("RMM_LOG_EVERY")) { + size_t v = std::strtoull(s, nullptr, 10); + if (v) { + LOG_EVERY = v; + } + } + + std::cerr << "[ LongRandom ] seed=" << seed << " CASES=" << L.CASES + << " OPS=" << L.OPS_PER_CASE << " MAX_N=" << L.MAX_N + << " LOG_EVERY=" << LOG_EVERY << "\n"; + + std::uniform_int_distribution coin(0, 1); + std::uniform_int_distribution len_u(1, (int)L.MAX_N); + std::uniform_int_distribution len_even(0, (int)(L.MAX_N / 2)); + + size_t total_ops = 0; + + for (size_t iter = 1; iter <= L.CASES; ++iter) { + std::string bits; + if (coin(rng) == 0) { + bits = random_bits(rng, (size_t)len_u(rng)); + } else { + size_t m = 1 + (size_t)len_even(rng); + bits = random_dyck_bits(rng, m); + if (bits.empty()) { + bits = "10"; + } + } + + run_case_and_compare(bits, rng, L.OPS_PER_CASE, seed); + total_ops += L.OPS_PER_CASE; + + if (iter % LOG_EVERY == 0) { + std::cerr << "[ LongRandom ] iter=" << iter << " total_ops=" << total_ops + << " last_N=" << bits.size() << " ok\n"; + } + } +} + +TEST(RmMTest, RankBasic) { + std::vector bits = {0b10110}; + pixie::RmMTree rm(bits, 5); + + EXPECT_EQ(rm.rank1(0), 0); // No bits + EXPECT_EQ(rm.rank1(1), 0); // 0 + EXPECT_EQ(rm.rank1(2), 1); // 10 + EXPECT_EQ(rm.rank1(3), 2); // 110 + EXPECT_EQ(rm.rank1(4), 2); // 0110 + EXPECT_EQ(rm.rank1(5), 3); // 10110 +} + +TEST(RmMTest, RankWithZeros) { + std::vector bits = {0}; + pixie::RmMTree rm(bits, 5); + + for (size_t i = 0; i <= 5; i++) { + EXPECT_EQ(rm.rank1(i), 0); + } +} + +TEST(RmMTest, SelectBasic) { + std::vector bits = {0b1100010110010110}; + pixie::RmMTree rm(bits, 16); + + EXPECT_EQ(rm.select1(1), 1); + EXPECT_EQ(rm.select1(2), 2); + EXPECT_EQ(rm.select1(3), 4); + EXPECT_EQ(rm.select1(4), 7); + EXPECT_EQ(rm.select1(5), 8); + EXPECT_EQ(rm.select1(6), 10); + EXPECT_EQ(rm.select1(7), 14); + EXPECT_EQ(rm.select1(8), 15); +} + +TEST(RmMTest, MainRankTest) { + std::mt19937_64 rng(42); + std::vector bits(65536 * 32); + for (size_t i = 0; i < 65536 * 32; i++) { + bits[i] = rng(); + } + + size_t rm_size = 65536 * 32 * 64; + pixie::RmMTree rm(bits, rm_size); + size_t rank = 0; + for (size_t i = 0; i < rm_size; ++i) { + ASSERT_EQ(rank, rm.rank1(i)); + rank += (bits[i >> 6] >> (i & 63)) & 1u; + } +} + +TEST(RmMTest, MainSelectTest) { + std::mt19937_64 rng(42); + std::vector bits(65536 * 32); + for (size_t i = 0; i < 65536 * 32; i++) { + bits[i] = rng(); + } + + size_t rm_size = 65536 * 32 * 64; + pixie::RmMTree rm(bits, rm_size); + size_t rank = 0; + + for (size_t i = 0; i < rm_size; ++i) { + if ((bits[i >> 6] >> (i & 63)) & 1u) { + ASSERT_EQ(rm.select1(++rank), i); + ASSERT_EQ(rm.rank1(i), rank - 1); + ASSERT_EQ(rm.rank1(i + 1), rank); + } + } +}