From 65434b3d744f4b1272449148cd2f9fa157263809 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=D0=A0=D0=BE=D0=B1=D0=B5=D1=80=D1=82=20=D0=A1=D0=BC=D0=B0?= =?UTF-8?q?=D0=B9=D1=82?= Date: Sat, 22 Nov 2025 19:00:06 +0300 Subject: [PATCH 1/4] faster searches --- README.md | 5 +- include/rmm_tree.h | 560 ++++++++++++++++++++++++++++++++++----------- src/bench_rmm.cpp | 4 +- 3 files changed, 430 insertions(+), 139 deletions(-) diff --git a/README.md b/README.md index 564058f..f1404b0 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,8 @@ 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=` (e.g. `./bench_rmm --benchmark_out=rmm_bench.csv`) and plot it with `misc/plot_rmm.py`. +Results print to stdout as CSV by default. Redirect to a file (e.g. `./bench_rmm > rmm_bench.csv`) and visualize it with `misc/plot_rmm.py`. +For a human-readable console output, use `./bench_rmm --benchmark_format=console`. --- diff --git a/include/rmm_tree.h b/include/rmm_tree.h index 263b668..b2cf622 100644 --- a/include/rmm_tree.h +++ b/include/rmm_tree.h @@ -8,6 +8,9 @@ #include #include #include +#if defined(__AVX2__) +#include +#endif namespace pixie { /** @@ -71,6 +74,12 @@ class RmMTree { // segment (to handle "10" crossing) both needed for: rank10, select10. std::vector node_first_bit, node_last_bit; + // leaf_prefix[k * leaf_prefix_stride + j] = excess(Lb + j) - excess(Lb), j in + // [0..leaf_size], Lb = k * block_bits. needed for: + // fwdsearch/bwdsearch/close/open/enclose + std::vector leaf_prefix; + size_t leaf_prefix_stride = 0; + public: /** * @brief Sentinel for "not found". @@ -133,9 +142,10 @@ class RmMTree { 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); + size_t nodes_buf[64]; + const size_t cnt = cover_blocks_collect(0, blk - 1, nodes_buf); + for (size_t j = 0; j < cnt; ++j) { + ans += ones_in_node(nodes_buf[j]); } } const size_t Lb = blk * block_bits; @@ -163,8 +173,7 @@ class RmMTree { return npos; } size_t base = 0; - const size_t leaf0 = first_leaf_index(); - while (v < leaf0) { + while (v < first_leaf_index) { const size_t Lc = v << 1, Rc = Lc | 1; const uint32_t o_l = ones_in_node(Lc); if (o_l >= k) { @@ -195,8 +204,7 @@ class RmMTree { return npos; } size_t base = 0; - const size_t leaf0 = first_leaf_index(); - while (v < leaf0) { + while (v < first_leaf_index) { const size_t Lc = v << 1, Rc = Lc | 1; const size_t z_l = zeros(Lc); if (z_l >= k) { @@ -255,8 +263,7 @@ class RmMTree { return npos; } size_t base = 0; - const size_t leaf0 = first_leaf_index(); - while (ind < leaf0) { + while (ind < first_leaf_index) { 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; @@ -297,29 +304,24 @@ class RmMTree { 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; - } + int delta_leaf = 0; + const size_t pos_leaf = leaf_fwd_bp_simd(blk, Lb, i, d, delta_leaf); + if (pos_leaf != npos) { + return pos_leaf; } - // 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); + int need = d - delta_leaf; + if (blk + 1 < leaf_count) { + size_t nodes_buf[64]; + const size_t cnt = + cover_blocks_collect(blk + 1, leaf_count - 1, nodes_buf); size_t base = (blk + 1) * block_bits; - for (const size_t& v : nodes) { + for (size_t j = 0; j < cnt; ++j) { + const size_t v = nodes_buf[j]; if (need == 0) { return base; } @@ -344,30 +346,18 @@ class RmMTree { 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; + int delta_leaf = 0; // excess(i) - excess(Lb) + const size_t res = leaf_bwd_bp_simd(blk, Lb, i, d, delta_leaf); + if (res != npos) { + return res; } - // 2) climb up - int need = target - excess_at_Lb; + // need = target - excess(Lb) = excess(i) + d - excess(Lb) = delta_leaf + d + int need = delta_leaf + d; size_t v = leaf_index_of(Lb); size_t base = Lb; while (v > 1) { @@ -437,10 +427,9 @@ class RmMTree { } // 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 l = first_leaf_index + (blk_i + 1); + size_t r = first_leaf_index + (blk_j - 1); size_t Rnodes[64]; int rn = 0; @@ -545,10 +534,9 @@ class RmMTree { } // 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 l = first_leaf_index + (blk_i + 1); + size_t r = first_leaf_index + (blk_j - 1); size_t Rnodes[64]; int rn = 0; @@ -724,9 +712,8 @@ class RmMTree { 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 l = first_leaf_index + blk_i + 1; + size_t r = first_leaf_index + blk_j - 1; size_t Rnodes[64]; int rn = 0; @@ -791,8 +778,8 @@ class RmMTree { // middle pref = cur_first; if (blk_i + 1 <= blk_j - 1) { - l = leaf0 + (blk_i + 1); - r = leaf0 + (blk_j - 1); + l = first_leaf_index + (blk_i + 1); + r = first_leaf_index + (blk_j - 1); rn = 0; while (l <= r) { if (l & 1) { @@ -1079,15 +1066,30 @@ class RmMTree { uint8_t pos_first_max; // pos of first maximum in this byte }; + struct LUT8Tables { + std::array 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, 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, 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& LUT8() noexcept { - static const std::array T = [] { - std::array t{}; + static inline const LUT8Tables& LUT8_ALL() noexcept { + static const LUT8Tables T = [] { + LUT8Tables tab{}; 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; + int pref[8]; const auto get = [&](const int& k) { return (b >> k) & 1; }; // LSB-first @@ -1097,6 +1099,7 @@ class RmMTree { ++rrc; } cur += bit ? +1 : -1; + pref[k] = cur; if (cur < mn) { mn = cur; cnt = 1; @@ -1119,13 +1122,138 @@ class RmMTree { a.last_bit = get(7); a.pos_first_min = pm; a.pos_first_max = pM; - t[b] = a; + tab.agg[b] = a; + auto& fp = tab.fwd_pos[b]; + auto& bp = tab.bwd_pos[b]; + fp.fill(-1); + bp.fill(-1); + for (int d = -8; d <= 8; ++d) { + for (int k = 0; k < 8; ++k) { + if (pref[k] == d) { + fp[d + 8] = k; + break; + } + } + for (int k = 7; k >= 0; --k) { + if (pref[k] == d) { + bp[d + 8] = k; + break; + } + } + } } - return t; + return tab; }(); return T; } + /** + * @brief Byte-level aggregates. + */ + static inline const std::array& LUT8() noexcept { + return LUT8_ALL().agg; + } + + /** + * @brief Position of the first occurrence of delta within the byte. + */ + static inline const std::array, 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, 256>& + LUT8_BWD_POS() noexcept { + return LUT8_ALL().bwd_pos; + } + + /** + * @brief Forward search within a single leaf over the range [from, to). + * @param need required relative excess (relative to `from`). + * @return The position, or npos. + */ + inline size_t scan_leaf_fwd(const size_t& from, + const size_t& to, + const int& need) const noexcept { + if (from >= to) { + return npos; + } + const auto& Agg = LUT8(); + const auto& Pos = LUT8_FWD_POS(); + int cur = 0; + size_t pos = from; + while (pos + 8 <= to) { + const uint8_t b = get_byte(pos); + const auto& a = Agg[b]; + const int local = need - cur; + if (local >= a.min_prefix && local <= a.max_prefix && local >= -8 && + local <= 8) { + const int8_t offset = Pos[b][local + 8]; + if (offset >= 0) { + return pos + size_t(offset); + } + } + cur += a.excess_total; + pos += 8; + } + + while (pos < to) { + cur += bit(pos) ? 1 : -1; + if (cur == need) { + return pos; + } + ++pos; + } + + return npos; + } + + /** + * @brief Backward search within a single leaf over [Lb, RB) (does not look to + * the right of RB). + * @param need relative excess from Lb. + * @param allow_rb whether RB may be used as a valid answer. + * @param right_border global right limit (as in descend_bwd). + */ + inline size_t scan_leaf_bwd(const size_t& Lb, + const size_t& RB, + const int& need, + const bool& allow_rb, + const size_t& right_border) const noexcept { + if (Lb >= RB) { + if ((Lb < right_border || allow_rb) && need == 0) { + return Lb; + } + return npos; + } + + const int16_t* lp = &leaf_prefix[block_of(Lb) * leaf_prefix_stride]; + const size_t end_incl = RB - Lb; + if (end_incl > 0) { +#if defined(__AVX2__) + const size_t pos = find_backward_equal_i16_avx2(lp, 1, end_incl, need); + if (pos != npos) { + return Lb + pos; + } +#else + for (size_t i = end_incl; i > 0; --i) { + if (lp[i] == need) { + return Lb + i; + } + } +#endif + } + + if ((Lb < right_border || allow_rb) && need == 0) { + return Lb; + } + + return npos; + } + /** * @brief Extract 8 bits starting at position pos (LSB-first across words). */ @@ -1150,8 +1278,7 @@ class RmMTree { * @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) { + while (v < first_leaf_index) { 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]; @@ -1178,9 +1305,7 @@ class RmMTree { /** * @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)); - } + size_t first_leaf_index = 1; /** * @brief Block index containing position i. @@ -1191,7 +1316,7 @@ class RmMTree { * @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); + return first_leaf_index + block_of(block_start); } /** @@ -1199,9 +1324,8 @@ class RmMTree { * @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; + if (v >= first_leaf_index) { + return (v - first_leaf_index) * block_bits; } size_t base = 0; @@ -1219,9 +1343,8 @@ class RmMTree { * @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; + size_t l = first_leaf_index + a; + size_t r = first_leaf_index + b; std::vector Lnodes, Rnodes; while (l <= r) { if ((l & 1) == 1) { @@ -1243,8 +1366,7 @@ class RmMTree { * 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) { + while (v < first_leaf_index) { const size_t Lc = v << 1; const size_t Rc = Lc | 1; if (node_min_prefix_excess[Lc] <= need && @@ -1256,15 +1378,8 @@ class RmMTree { 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; + return scan_leaf_fwd(base, std::min(base + segment_size_bits[v], num_bits), + need); } /** @@ -1280,8 +1395,7 @@ class RmMTree { const int& need, const size_t& right_border, const bool& allow_rb) const noexcept { - const size_t leaf0 = first_leaf_index(); - while (v < leaf0) { + while (v < first_leaf_index) { const size_t Lc = v << 1; const size_t Rc = Lc | 1; const int need_r = need - node_total_excess[Lc]; @@ -1319,35 +1433,10 @@ class RmMTree { 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( + base, + std::min(right_border, std::min(base + segment_size_bits[v], num_bits)), + need, allow_rb, right_border); } /** @@ -1355,8 +1444,7 @@ class RmMTree { * (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) { + while (v < first_leaf_index) { 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]; @@ -1388,8 +1476,7 @@ class RmMTree { int d, size_t q, size_t base) const noexcept { - const size_t leaf0 = first_leaf_index(); - while (v < leaf0) { + while (v < first_leaf_index) { const size_t Lc = v << 1; const size_t Rc = Lc | 1; const int leftm = node_min_prefix_excess[Lc]; @@ -1529,24 +1616,13 @@ class RmMTree { * @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); + return std::countr_zero(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 } /** @@ -1576,7 +1652,8 @@ class RmMTree { 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); + 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) { @@ -1725,6 +1802,42 @@ class RmMTree { } } + /** + * @brief Cover the range of whole blocks [a..b] (inclusive) with tree nodes, + * writing them to `out` from left to right. + * @return The number of nodes written. + */ + inline size_t cover_blocks_collect(const size_t& a, + const size_t& b, + size_t (&out)[64]) const noexcept { + if (leaf_count == 0 || a > b) { + return 0; + } + size_t l = first_leaf_index + a; + size_t r = first_leaf_index + b; + size_t Lnodes[32]; + size_t Rnodes[32]; + size_t Ln = 0, Rn = 0; + while (l <= r) { + if (l & 1) { + Lnodes[Ln++] = l++; + } + if ((r & 1) == 0) { + Rnodes[Rn++] = r--; + } + l >>= 1; + r >>= 1; + } + size_t ind = 0; + for (size_t i = 0; i < Ln; ++i) { + out[ind++] = Lnodes[i]; + } + while (Rn > 0) { + out[ind++] = Rnodes[--Rn]; + } + return ind; + } + /** * @brief Select q-th minimum (1-based) inside [l..r] inclusive in two 8-bit * passes. @@ -1811,6 +1924,174 @@ class RmMTree { return npos; } + /** + * @brief Efficiently searches for the first occurrence of a 16-bit value in + * the range [begin, end_excl) using AVX2 when available. + * @details Loads 16 consecutive int16_t elements (256 bits) per iteration. + * Compares them against the 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 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) noexcept { +#if defined(__AVX2__) + 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(arr + i)), + vtarget)); + if (mask) { + return i + (std::countr_zero(mask) >> 1); + } + } + for (; i < n; ++i) { + if (arr[i] == target) { + return i; + } + } + return npos; +#else + for (size_t i = begin; i < end_excl; ++i) { + if (arr[i] == target) { + return i; + } + } + return npos; +#endif + } + + /** + * @brief Performs a backward search for a 16-bit value in a given range. + * @details Scans the array segment [begin .. 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 `target`, or `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) noexcept { + if (begin > end_incl) { + return npos; + } +#if defined(__AVX2__) + 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(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; + } + } + return npos; + } else { + for (size_t i = end_incl + 1; i > begin;) { + --i; + if (arr[i] == target) { + return i; + } + } + return npos; + } +#else + for (size_t i = end_incl + 1; i > begin;) { + --i; + if (arr[i] == target) { + return i; + } + } + return npos; +#endif + } + + /** + * @brief Performs a forward search within a single leaf. + * @details Starting from position `i` (global index) inside the leaf [Lb .. + * Lb+len), it looks for the nearest position where the excess changes by `d` + * relative to `i`. The function uses a SIMD-accelerated search over the + * leaf’s prefix-sum array `lp`. If a matching position is found, returns its + * global index. If not found, returns `npos` and sets `delta_leaf` to the net + * excess change from `i` to the end of the leaf. + */ + inline size_t leaf_fwd_bp_simd(const size_t& leaf_ind, + const size_t& Lb, + const size_t& i, + const int& d, + int& delta_leaf) const noexcept { + const size_t len = segment_size_bits[first_leaf_index + leaf_ind]; + const int16_t* lp = &leaf_prefix[leaf_ind * leaf_prefix_stride]; + const size_t pos = i - Lb; + if (pos > len) { + delta_leaf = 0; + return npos; + } + const int16_t start_pref = lp[pos]; + size_t match_ind = + find_forward_equal_i16_avx2(lp, pos + 1, len + 1, start_pref + d); + if (match_ind != npos) { + return Lb + match_ind - 1; + } + delta_leaf = lp[len] - start_pref; + return npos; + } + + /** + * @brief Performs a backward search within a leaf segment for a position `p < + * i` such that the excess difference from `i` equals `d`. + * @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 position `i`. + */ + inline size_t leaf_bwd_bp_simd(const size_t& leaf_ind, + const size_t& Lb, + const size_t& i, + const int& d, + int& delta_leaf) const noexcept { + const int16_t* lp = &leaf_prefix[leaf_ind * leaf_prefix_stride]; + const size_t pos = i - Lb; + if (pos > segment_size_bits[first_leaf_index + leaf_ind]) { + delta_leaf = 0; + return npos; + } + const int16_t start_pref = lp[pos]; + if (pos > 0) { + const size_t match_ind = + find_backward_equal_i16_avx2(lp, 0, pos - 1, start_pref + d); + if (match_ind != npos) { + return Lb + match_ind; + } + } + delta_leaf = start_pref; + return npos; + } + /** * @brief Find first minimum value and its first position in [l..r] inclusive * using 8-bit LUT. @@ -1941,8 +2222,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(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,10 +2232,12 @@ 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 v = first_leaf_index + 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; @@ -2032,9 +2315,16 @@ class RmMTree { 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; + int16_t* lp = &leaf_prefix[k * leaf_prefix_stride]; + int16_t acc = 0; + lp[0] = 0; + for (size_t pos = Lb, i = 1; pos < Rb; ++pos, ++i) { + acc += bit(pos) ? 1 : -1; + lp[i] = acc; + } } // internal nodes - for (size_t v = leaf0 - 1; v >= 1; --v) { + for (size_t v = first_leaf_index - 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]; diff --git a/src/bench_rmm.cpp b/src/bench_rmm.cpp index 06db0db..568ad68 100644 --- a/src/bench_rmm.cpp +++ b/src/bench_rmm.cpp @@ -435,8 +435,8 @@ int main(int argc, char** argv) { }; static vector extra; - if (!has("--benchmark_out_format")) { - extra.emplace_back("--benchmark_out_format=csv"); + if (!has("--benchmark_format")) { + extra.emplace_back("--benchmark_format=csv"); } if (!has("--benchmark_counters_tabular")) { extra.emplace_back("--benchmark_counters_tabular=true"); From 1d90d18122cce6b1065f5257a5c333138786434c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=D0=A0=D0=BE=D0=B1=D0=B5=D1=80=D1=82=20=D0=A1=D0=BC=D0=B0?= =?UTF-8?q?=D0=B9=D1=82?= Date: Sat, 29 Nov 2025 17:00:47 +0300 Subject: [PATCH 2/4] refactoring --- README.md | 3 +- include/bits.h | 105 ++++++++++++++++ include/rmm_tree.h | 307 +++++++++++++++------------------------------ misc/plot_rmm.py | 27 ++-- src/bench_rmm.cpp | 214 ++++++++++++++++--------------- 5 files changed, 335 insertions(+), 321 deletions(-) diff --git a/README.md b/README.md index f1404b0..516fa64 100644 --- a/README.md +++ b/README.md @@ -70,8 +70,7 @@ Benchmarks are random 50/50 0-1 bitvectors up to $2^{34}$ bits. ./bench_rmm ``` -Results print to stdout as CSV by default. Redirect to a file (e.g. `./bench_rmm > rmm_bench.csv`) and visualize it with `misc/plot_rmm.py`. -For a human-readable console output, use `./bench_rmm --benchmark_format=console`. +For visualization, write the JSON output to a file using `--benchmark_out=` (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..a85b8fe 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 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(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 `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(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 b2cf622..3d5ae0c 100644 --- a/include/rmm_tree.h +++ b/include/rmm_tree.h @@ -1,4 +1,6 @@ #pragma once +#include + #include #include #include @@ -8,9 +10,8 @@ #include #include #include -#if defined(__AVX2__) -#include -#endif + +#include "bits.h" namespace pixie { /** @@ -132,8 +133,8 @@ class RmMTree { // --------- 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 i). + * @details Returns 0 for @p i == 0. */ size_t rank1(const size_t& i) const { if (i == 0) { @@ -155,14 +156,14 @@ class RmMTree { } /** - * @brief Number of zeros in prefix [0, i). - * @details Computed as i - rank1(i). + * @brief Number of zeros in prefix [0, @p i). + * @details Computed as @p i - rank1(@p 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. + * @brief 1-based select of the @p k-th one. + * @return Position of @p k-th '1' or npos if not found. */ size_t select1(size_t k) const { if (k == 0 || num_bits == 0) { @@ -189,8 +190,8 @@ class RmMTree { } /** - * @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 k-th zero. + * @return Position of @p k-th '0' or npos if not found. */ size_t select0(size_t k) const { if (k == 0 || num_bits == 0) { @@ -220,8 +221,8 @@ class RmMTree { } /** - * @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 j || j >= num_bits) { @@ -501,7 +505,8 @@ class RmMTree { } /** - * @brief Position of the first maximum of excess on [i, j] (inclusive). + * @brief Position of the first maximum of excess on [@p i, @p 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 { @@ -593,7 +598,8 @@ class RmMTree { } /** - * @brief Value of the maximum prefix excess on [i, j] relative to i. + * @brief Value of the maximum prefix excess on [@p i, @p j] relative to + * @p i. */ int range_max_query_val(const size_t& i, const size_t& j) const { if (i > j || j >= num_bits) { @@ -607,7 +613,7 @@ class RmMTree { } /** - * @brief How many times the minimum prefix excess occurs on [i, j]. + * @brief How many times the minimum prefix excess occurs on [@p i, @p j]. */ size_t mincount(const size_t& i, const size_t& j) const { if (i > j || j >= num_bits) { @@ -681,8 +687,9 @@ class RmMTree { } /** - * @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 q-th (1-based) occurrence of the minimum on + * [@p i, @p j]. + * @return Position or npos if @p 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) { @@ -823,7 +830,7 @@ class RmMTree { // ----- parentheses navigation (BP) ----- /** - * @brief close(i): matching ')' for '(' at i. + * @brief close(@p i): matching ')' for '(' at @p i. * @return Position of matching ')', or npos. */ inline size_t close(const size_t& i) const { @@ -834,7 +841,7 @@ class RmMTree { } /** - * @brief open(i): matching '(' for ')' at i. + * @brief open(@p i): matching '(' for ')' at @p i. * @return Position of matching '(', or npos. */ inline size_t open(const size_t& i) const { @@ -847,7 +854,7 @@ class RmMTree { } /** - * @brief enclose(i): opening '(' that strictly encloses position i. + * @brief enclose(@p i): opening '(' that strictly encloses position @p i. * @return Position of enclosing '(', or npos. */ inline size_t enclose(const size_t& i) const { @@ -861,8 +868,8 @@ class RmMTree { 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 len. + * @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 { @@ -879,8 +886,8 @@ class RmMTree { } /** - * @brief Rank of ones within [Lb, Rb). - * @details Works on word boundaries; Rb may equal Lb. + * @brief Rank of ones within [@p Lb, @p Rb). + * @details Works on word boundaries; @p Rb may equal @p Lb. */ size_t rank1_in_block(const size_t& Lb, const size_t& Rb) const noexcept { if (Rb <= Lb) { @@ -913,7 +920,7 @@ class RmMTree { } /** - * @brief Count "10" starts within [Lb, Rb). + * @brief Count "10" starts within [@p Lb, @p Rb). * @details Accounts for cross-word boundaries. */ size_t rr_in_block(const size_t& Lb, const size_t& Rb) const noexcept { @@ -961,7 +968,7 @@ class RmMTree { } /** - * @brief 1-based select of k-th "10" within [Lb, Rb). + * @brief 1-based select of @p k-th "10" within [@p Lb, @p Rb). * @return Position or npos. */ size_t select10_in_block(const size_t& Lb, @@ -1171,8 +1178,8 @@ class RmMTree { } /** - * @brief Forward search within a single leaf over the range [from, to). - * @param need required relative excess (relative to `from`). + * @brief Forward search within a single leaf over the range [@p from, @p to). + * @param need required relative excess (relative to @p from). * @return The position, or npos. */ inline size_t scan_leaf_fwd(const size_t& from, @@ -1212,10 +1219,10 @@ class RmMTree { } /** - * @brief Backward search within a single leaf over [Lb, RB) (does not look to - * the right of RB). - * @param need relative excess from Lb. - * @param allow_rb whether RB may be used as a valid answer. + * @brief Backward search within a single leaf over [@p Lb, @p RB) (does not + * look to the right of @p RB). + * @param need relative excess from @p Lb. + * @param allow_rb whether @p RB may be used as a valid answer. * @param right_border global right limit (as in descend_bwd). */ inline size_t scan_leaf_bwd(const size_t& Lb, @@ -1233,18 +1240,11 @@ class RmMTree { const int16_t* lp = &leaf_prefix[block_of(Lb) * leaf_prefix_stride]; const size_t end_incl = RB - Lb; if (end_incl > 0) { -#if defined(__AVX2__) - const size_t pos = find_backward_equal_i16_avx2(lp, 1, end_incl, need); + const size_t pos = + ::find_backward_equal_i16_avx2(lp, 1, end_incl, need, npos); if (pos != npos) { return Lb + pos; } -#else - for (size_t i = end_incl; i > 0; --i) { - if (lp[i] == need) { - return Lb + i; - } - } -#endif } if ((Lb < right_border || allow_rb) && need == 0) { @@ -1255,7 +1255,8 @@ class RmMTree { } /** - * @brief Extract 8 bits starting at position pos (LSB-first across words). + * @brief Extract 8 bits starting at position @p pos (LSB-first across + * words). */ inline uint8_t get_byte(const size_t& pos) const noexcept { const size_t w = pos >> 6; @@ -1271,7 +1272,7 @@ class RmMTree { /** * @brief Descend to the first (leftmost) maximum with node-relative prefix - * equal to d. + * equal to @p d. * @param v Node index. * @param d Target prefix within node coordinates. * @param base Starting global position of node. @@ -1308,19 +1309,19 @@ class RmMTree { size_t first_leaf_index = 1; /** - * @brief Block index containing position i. + * @brief Block index containing position @p 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. + * @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); } /** - * @brief Starting global position for node v (0-indexed). + * @brief Starting global position for node @p v (0-indexed). * @details Walk up to compute base for internal nodes; direct for leaves. */ size_t node_base(size_t v) const noexcept { @@ -1338,8 +1339,8 @@ class RmMTree { } /** - * @brief Cover a range of whole blocks [a..b] (inclusive) with O(log n) - * maximal nodes. + * @brief Cover a range of whole blocks [@p a..@p 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 { @@ -1363,7 +1364,7 @@ class RmMTree { /** * @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 { while (v < first_leaf_index) { @@ -1385,9 +1386,9 @@ class RmMTree { /** * @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 base Left border of node @p v. + * @param need Target relative prefix inside @p v. + * @param right_border Global right limit (exclusive if !@p allow_rb). * @param allow_rb Whether right border is allowed to match. */ size_t descend_bwd(size_t v, @@ -1440,8 +1441,8 @@ class RmMTree { } /** - * @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 d (minimum). */ size_t descend_first_min(size_t v, int d, size_t base) const noexcept { while (v < first_leaf_index) { @@ -1469,8 +1470,8 @@ class RmMTree { } /** - * @brief Descend to find the q-th minimum (1-based) where node-relative - * prefix equals d. + * @brief Descend to find the @p q-th minimum (1-based) where node-relative + * prefix equals @p d. */ size_t descend_qth_min(size_t v, int d, @@ -1501,7 +1502,7 @@ class RmMTree { } /** - * @brief 1-based select of k-th '1' within [Lb, Rb). + * @brief 1-based select of @p k-th '1' within [@p Lb, @p Rb). */ size_t select1_in_block(const size_t& Lb, const size_t& Rb, @@ -1550,7 +1551,7 @@ class RmMTree { } /** - * @brief 1-based select of k-th '0' within [Lb, Rb). + * @brief 1-based select of @p k-th '0' within [@p Lb, @p Rb). */ size_t select0_in_block(const size_t& Lb, const size_t& Rb, @@ -1612,7 +1613,7 @@ class RmMTree { } /** - * @brief 1-based select of k-th set bit inside a 64-bit word. + * @brief 1-based select of @p 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 { @@ -1665,8 +1666,9 @@ class RmMTree { } /** - * @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 cap. + * @details Returns 64 if @p cap < 0 (no constraint). Clamped to ≤16384 or + * @p Nbits. */ static inline size_t choose_block_bits_for_overhead( const size_t& Nbits, @@ -1724,30 +1726,31 @@ class RmMTree { } /** - * @brief Read bit at position i (LSB-first across words). + * @brief Read bit at position @p 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. + * @brief Set bit at position @p 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. + * @brief Number of ones in node @p 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. + * @brief Scan [@p l..@p r] inclusive, computing minimum value and how many + * times it is attained. + * @details Uses 8-bit LUT for speed; outputs @p cur at end, @p mn, and + * @p cnt. */ inline void scan_range_min_count8(size_t l, const size_t& r, @@ -1803,8 +1806,8 @@ class RmMTree { } /** - * @brief Cover the range of whole blocks [a..b] (inclusive) with tree nodes, - * writing them to `out` from left to right. + * @brief Cover the range of whole blocks [@p a..@p b] (inclusive) with tree + * nodes, writing them to @p out from left to right. * @return The number of nodes written. */ inline size_t cover_blocks_collect(const size_t& a, @@ -1839,9 +1842,10 @@ class RmMTree { } /** - * @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 q-th minimum (1-based) inside [@p l..@p r] inclusive in + * two 8-bit passes. + * @details First pass finds global minimum, second selects the @p q-th + * position. */ inline size_t qth_min_in_block(const size_t& l, const size_t& r, @@ -1924,121 +1928,14 @@ class RmMTree { return npos; } - /** - * @brief Efficiently searches for the first occurrence of a 16-bit value in - * the range [begin, end_excl) using AVX2 when available. - * @details Loads 16 consecutive int16_t elements (256 bits) per iteration. - * Compares them against the 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 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) noexcept { -#if defined(__AVX2__) - 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(arr + i)), - vtarget)); - if (mask) { - return i + (std::countr_zero(mask) >> 1); - } - } - for (; i < n; ++i) { - if (arr[i] == target) { - return i; - } - } - return npos; -#else - for (size_t i = begin; i < end_excl; ++i) { - if (arr[i] == target) { - return i; - } - } - return npos; -#endif - } - - /** - * @brief Performs a backward search for a 16-bit value in a given range. - * @details Scans the array segment [begin .. 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 `target`, or `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) noexcept { - if (begin > end_incl) { - return npos; - } -#if defined(__AVX2__) - 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(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; - } - } - return npos; - } else { - for (size_t i = end_incl + 1; i > begin;) { - --i; - if (arr[i] == target) { - return i; - } - } - return npos; - } -#else - for (size_t i = end_incl + 1; i > begin;) { - --i; - if (arr[i] == target) { - return i; - } - } - return npos; -#endif - } - /** * @brief Performs a forward search within a single leaf. - * @details Starting from position `i` (global index) inside the leaf [Lb .. - * Lb+len), it looks for the nearest position where the excess changes by `d` - * relative to `i`. The function uses a SIMD-accelerated search over the - * leaf’s prefix-sum array `lp`. If a matching position is found, returns its - * global index. If not found, returns `npos` and sets `delta_leaf` to the net - * excess change from `i` to the end of the leaf. + * @details Starting from position @p i (global index) inside the leaf [@p Lb + * .. Lb+len), it looks for the nearest position where the excess changes by + * @p d relative to @p i. 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 delta_leaf to the net + * excess change from @p i to the end of the leaf. */ inline size_t leaf_fwd_bp_simd(const size_t& leaf_ind, const size_t& Lb, @@ -2053,8 +1950,8 @@ class RmMTree { return npos; } const int16_t start_pref = lp[pos]; - size_t match_ind = - find_forward_equal_i16_avx2(lp, pos + 1, len + 1, start_pref + d); + size_t match_ind = ::find_forward_equal_i16_avx2(lp, pos + 1, len + 1, + start_pref + d, npos); if (match_ind != npos) { return Lb + match_ind - 1; } @@ -2063,11 +1960,11 @@ class RmMTree { } /** - * @brief Performs a backward search within a leaf segment for a position `p < - * i` such that the excess difference from `i` equals `d`. + * @brief Performs a backward search within a leaf segment for a position + * earlier than @p i such that the excess difference from @p i equals @p d. * @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 position `i`. + * and outputs the leaf-local excess delta at position @p i. */ inline size_t leaf_bwd_bp_simd(const size_t& leaf_ind, const size_t& Lb, @@ -2083,7 +1980,7 @@ class RmMTree { const int16_t start_pref = lp[pos]; if (pos > 0) { const size_t match_ind = - find_backward_equal_i16_avx2(lp, 0, pos - 1, start_pref + d); + ::find_backward_equal_i16_avx2(lp, 0, pos - 1, start_pref + d, npos); if (match_ind != npos) { return Lb + match_ind; } @@ -2093,8 +1990,8 @@ class RmMTree { } /** - * @brief Find first minimum value and its first position in [l..r] inclusive - * using 8-bit LUT. + * @brief Find first minimum value and its first position in [@p l..@p 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). */ @@ -2143,8 +2040,8 @@ class RmMTree { } /** - * @brief Find first maximum value and its first position in [l..r] inclusive - * using 8-bit LUT. + * @brief Find first maximum value and its first position in [@p l..@p 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). */ 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 568ad68..523516e 100644 --- a/src/bench_rmm.cpp +++ b/src/bench_rmm.cpp @@ -3,6 +3,7 @@ #include #include #include +#include #include #include #include @@ -10,42 +11,42 @@ #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; + 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 explicit_sizes; + std::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 + "="; + auto getv = [&](const std::string& key) -> std::optional { + 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 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 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 build_size_grid() { +static std::vector build_size_grid() { if (!args.explicit_sizes.empty()) { - vector v = args.explicit_sizes; - v.erase(remove(v.begin(), v.end(), 0), v.end()); + std::vector 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 v; + std::vector 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 N_cands; - size_t min_N = size_t(1) << lo, max_N = size_t(1) << hi; + std::set 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 build_size_grid() { } } - vector v(N_cands.begin(), N_cands.end()); + std::vector v(N_cands.begin(), N_cands.end()); if (v.front() != min_N) { v.insert(v.begin(), min_N); } @@ -156,33 +160,33 @@ static vector build_size_grid() { } struct Pools { - vector inds_any; - vector inds; - vector inds_1N; - vector deltas; - vector> segs; - - vector ks1, ks0, ks10; - vector minselect_q; + std::vector inds_any; + std::vector inds; + std::vector inds_1N; + std::vector deltas; + std::vector> segs; + + std::vector ks1, ks0, ks10; + std::vector 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> keepalive; +static std::vector> 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(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(args.Q, 32768); + const std::size_t L = std::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); + std::uniform_int_distribution ind_dist_incl(0, N ? N : 0); + std::uniform_int_distribution ind_dist(0, N ? (N - 1) : 0); + std::uniform_int_distribution ind_dist_1N(1, N ? N : 1); + std::uniform_int_distribution 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 { + auto rand_ij = [&]() -> std::pair { 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& out) { + auto fill_ks = [&](std::size_t total, std::vector& out) { out.resize(L); if (total == 0) { - fill(out.begin(), out.end(), 1); + std::fill(out.begin(), out.end(), 1); return; } - uniform_int_distribution dist(1, total); - for (size_t i = 0; i < L; ++i) { + std::uniform_int_distribution 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 uq(1, c); + std::uniform_int_distribution uq(1, c); d.pool.minselect_q[i] = uq(rng); } @@ -269,14 +274,16 @@ static Dataset build_dataset(size_t N) { } template -static void register_op(const string& op, shared_ptr data, Fn&& body) { - auto idx_ptr = make_shared(0); +static void register_op(const std::string& op, + std::shared_ptr data, + Fn&& body) { + auto idx_ptr = std::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)++; + std::size_t i = (*idx_ptr)++; body(state, D, i); } state.counters["N"] = static_cast(D.N); @@ -289,106 +296,106 @@ static void register_op(const string& op, shared_ptr data, Fn&& body) { static void register_all() { auto Ns = build_size_grid(); - for (size_t N : Ns) { - auto data = make_shared(build_dataset(N)); + for (std::size_t N : Ns) { + auto data = std::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) { + [&](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 extra; - if (!has("--benchmark_format")) { - extra.emplace_back("--benchmark_format=csv"); + static std::vector extra; + if (!has("--benchmark_out_format")) { + 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 argv_vec; + static std::vector argv_vec; argv_vec.assign(argv, argv + argc); for (auto& s : extra) { argv_vec.push_back(s.data()); From 694b03a9ec1f89dbafd689df3e123eb45154a97d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=D0=A0=D0=BE=D0=B1=D0=B5=D1=80=D1=82=20=D0=A1=D0=BC=D0=B0?= =?UTF-8?q?=D0=B9=D1=82?= Date: Sat, 29 Nov 2025 17:35:29 +0300 Subject: [PATCH 3/4] annihilation of short variable names --- include/rmm_tree.h | 2718 +++++++++++++++++++++++--------------------- 1 file changed, 1446 insertions(+), 1272 deletions(-) diff --git a/include/rmm_tree.h b/include/rmm_tree.h index 3d5ae0c..30b2ad9 100644 --- a/include/rmm_tree.h +++ b/include/rmm_tree.h @@ -75,8 +75,9 @@ class RmMTree { // segment (to handle "10" crossing) both needed for: rank10, select10. std::vector node_first_bit, node_last_bit; - // leaf_prefix[k * leaf_prefix_stride + j] = excess(Lb + j) - excess(Lb), j in - // [0..leaf_size], Lb = k * block_bits. needed for: + // 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 leaf_prefix; size_t leaf_prefix_stride = 0; @@ -100,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). @@ -124,704 +125,775 @@ class RmMTree { * leaf_block_bits, (3) set to ceil_pow2(log2(num_bits)). */ explicit RmMTree(const std::vector& 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, @p i). - * @details Returns 0 for @p i == 0. + * @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) { - size_t nodes_buf[64]; - const size_t cnt = cover_blocks_collect(0, blk - 1, nodes_buf); - for (size_t j = 0; j < cnt; ++j) { - ans += ones_in_node(nodes_buf[j]); - } - } - 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, @p i). - * @details Computed as @p i - rank1(@p 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 @p k-th one. - * @return Position of @p 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; - while (v < first_leaf_index) { - 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 @p k-th zero. - * @return Position of @p 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; - while (v < first_leaf_index) { - 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, @p i). - * @details Counts p where bit[p]==1 and bit[p+1]==0 with p+1<@p 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 @p 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; - while (ind < first_leaf_index) { - 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, @p 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 ≥ @p i where excess(p) = - * excess(@p i) + @p 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; } // 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 delta_leaf = 0; - const size_t pos_leaf = leaf_fwd_bp_simd(blk, Lb, i, d, delta_leaf); - if (pos_leaf != npos) { - return pos_leaf; - } - - int need = d - delta_leaf; - if (blk + 1 < leaf_count) { - size_t nodes_buf[64]; - const size_t cnt = - cover_blocks_collect(blk + 1, leaf_count - 1, nodes_buf); - size_t base = (blk + 1) * block_bits; - for (size_t j = 0; j < cnt; ++j) { - const size_t v = nodes_buf[j]; - 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 ≤ @p i where excess(p) = - * excess(@p i) + @p 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; } // 1) scan inside the block - const size_t blk = block_of(i - 1); - const size_t Lb = blk * block_bits; - int delta_leaf = 0; // excess(i) - excess(Lb) - const size_t res = leaf_bwd_bp_simd(blk, Lb, i, d, delta_leaf); - if (res != npos) { - return res; - } - - // need = target - excess(Lb) = excess(i) + d - excess(Lb) = delta_leaf + d - int need = delta_leaf + d; - 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 [@p i, @p j] + * @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 - if (blk_i + 1 <= blk_j - 1) { - size_t l = first_leaf_index + (blk_i + 1); - size_t r = first_leaf_index + (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 [@p i, @p j] relative to - * @p i. - * @details Equivalent to min_{t in [@p i..@p j]} (excess(t+1) - excess(@p - * 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 [@p i, @p j] + * @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 - if (blk_i + 1 <= blk_j - 1) { - size_t l = first_leaf_index + (blk_i + 1); - size_t r = first_leaf_index + (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 [@p i, @p j] relative to - * @p 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 [@p i, @p 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 @p q-th (1-based) occurrence of the minimum on - * [@p i, @p j]. - * @return Position or npos if @p 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 - size_t l = first_leaf_index + blk_i + 1; - size_t r = first_leaf_index + 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 = first_leaf_index + (blk_i + 1); - r = first_leaf_index + (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; @@ -830,232 +902,249 @@ class RmMTree { // ----- parentheses navigation (BP) ----- /** - * @brief close(@p i): matching ')' for '(' at @p 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(@p i): matching '(' for ')' at @p 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(@p i): opening '(' that strictly encloses position @p 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 @p len. + * 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 [@p Lb, @p Rb). - * @details Works on word boundaries; @p Rb may equal @p 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 [@p Lb, @p 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 @p k-th "10" within [@p Lb, @p 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; @@ -1091,67 +1180,68 @@ class RmMTree { * @brief Unified initialization of all byte-based LUTs. */ static inline const LUT8Tables& LUT8_ALL() noexcept { - static const LUT8Tables T = [] { - LUT8Tables tab{}; - 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; - int pref[8]; - const auto get = [&](const int& k) { - return (b >> k) & 1; + 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; } - cur += bit ? +1 : -1; - pref[k] = cur; - if (cur < mn) { - mn = cur; - cnt = 1; - pm = k; - } else if (cur == mn) { - ++cnt; + 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 (cur > mx) { - mx = cur; - pM = k; + if (current_excess > max_prefix) { + max_prefix = current_excess; + first_max_position = bit_index; } } - 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; - tab.agg[b] = a; - auto& fp = tab.fwd_pos[b]; - auto& bp = tab.bwd_pos[b]; - fp.fill(-1); - bp.fill(-1); - for (int d = -8; d <= 8; ++d) { - for (int k = 0; k < 8; ++k) { - if (pref[k] == d) { - fp[d + 8] = k; + 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; } } - for (int k = 7; k >= 0; --k) { - if (pref[k] == d) { - bp[d + 8] = k; + for (int bit_index = 7; bit_index >= 0; --bit_index) { + if (prefixes[bit_index] == delta) { + backward_positions[delta + 8] = bit_index; break; } } } } - return tab; + return lookup_tables; }(); - return T; + return tables; } /** @@ -1178,129 +1268,144 @@ class RmMTree { } /** - * @brief Forward search within a single leaf over the range [@p from, @p to). - * @param need required relative excess (relative to @p from). + * @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& from, - const size_t& to, - const int& need) const noexcept { - if (from >= to) { + 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& Agg = LUT8(); - const auto& Pos = LUT8_FWD_POS(); - int cur = 0; - size_t pos = from; - while (pos + 8 <= to) { - const uint8_t b = get_byte(pos); - const auto& a = Agg[b]; - const int local = need - cur; - if (local >= a.min_prefix && local <= a.max_prefix && local >= -8 && - local <= 8) { - const int8_t offset = Pos[b][local + 8]; + 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 pos + size_t(offset); + return position + size_t(offset); } } - cur += a.excess_total; - pos += 8; + current_excess += byte_aggregate.excess_total; + position += 8; } - while (pos < to) { - cur += bit(pos) ? 1 : -1; - if (cur == need) { - return pos; + while (position < search_end) { + current_excess += bit(position) ? 1 : -1; + if (current_excess == required_delta) { + return position; } - ++pos; + ++position; } return npos; } /** - * @brief Backward search within a single leaf over [@p Lb, @p RB) (does not - * look to the right of @p RB). - * @param need relative excess from @p Lb. - * @param allow_rb whether @p RB may be used as a valid answer. - * @param right_border global right limit (as in descend_bwd). + * @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& Lb, - const size_t& RB, - const int& need, - const bool& allow_rb, - const size_t& right_border) const noexcept { - if (Lb >= RB) { - if ((Lb < right_border || allow_rb) && need == 0) { - return Lb; + 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* lp = &leaf_prefix[block_of(Lb) * leaf_prefix_stride]; - const size_t end_incl = RB - Lb; - if (end_incl > 0) { - const size_t pos = - ::find_backward_equal_i16_avx2(lp, 1, end_incl, need, npos); - if (pos != npos) { - return Lb + pos; + 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 ((Lb < right_border || allow_rb) && need == 0) { - return Lb; + if ((block_begin < global_right_border || allow_right_boundary) && + required_delta == 0) { + return block_begin; } return npos; } /** - * @brief Extract 8 bits starting at position @p pos (LSB-first across - * words). + * @brief Extract 8 bits starting at @p position (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); + 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 @p 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 { - while (v < first_leaf_index) { - 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); } /** @@ -1309,9 +1414,11 @@ class RmMTree { size_t first_leaf_index = 1; /** - * @brief Block index containing position @p 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 @p block_start. @@ -1321,307 +1428,336 @@ class RmMTree { } /** - * @brief Starting global position for node @p 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 { - if (v >= first_leaf_index) { - return (v - first_leaf_index) * 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 [@p a..@p 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 cover_blocks(const size_t& a, const size_t& b) const { - size_t l = first_leaf_index + a; - size_t r = first_leaf_index + 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; + std::vector 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 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 @p need. */ - size_t descend_fwd(size_t v, int need, size_t base) const noexcept { - while (v < first_leaf_index) { - 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; } } - return scan_leaf_fwd(base, std::min(base + segment_size_bits[v], num_bits), - need); + 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 @p v. - * @param need Target relative prefix inside @p v. - * @param right_border Global right limit (exclusive if !@p 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 { - while (v < first_leaf_index) { - 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; } return scan_leaf_bwd( - base, - std::min(right_border, std::min(base + segment_size_bits[v], num_bits)), - need, allow_rb, right_border); + 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 - * @p d (minimum). + * @p target_prefix (minimum). */ - size_t descend_first_min(size_t v, int d, size_t base) const noexcept { - while (v < first_leaf_index) { - 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 @p q-th minimum (1-based) where node-relative - * prefix equals @p 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 { - while (v < first_leaf_index) { - 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 @p k-th '1' within [@p Lb, @p 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 @p k-th '0' within [@p Lb, @p 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 @p 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 { - while (w) { - if (--k == 0) { - return std::countr_zero(w); + 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); } - w &= (w - 1); + word &= (word - 1); } return -1; } @@ -1629,63 +1765,67 @@ class RmMTree { /** * @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(1, leaf_count)) + leaf_count; + size_t leaf_node_count = ceil_div(bit_count, block_size_pow2); + return std::bit_ceil(std::max(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(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 ≤ @p cap. - * @details Returns 64 if @p cap < 0 (no constraint). Clamped to ≤16384 or - * @p 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(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(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; } /** @@ -1693,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); } } @@ -1709,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& 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); } @@ -1726,203 +1866,209 @@ class RmMTree { } /** - * @brief Read bit at position @p 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 @p 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 @p 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 [@p l..@p r] inclusive, computing minimum value and how many + * @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 cur at end, @p mn, and - * @p cnt. + * @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; } - ++l; + ++range_begin; } - if (mn == INT_MAX) { - mn = cnt = 0; + if (min_value == INT_MAX) { + min_value = count = 0; } } /** - * @brief Cover the range of whole blocks [@p a..@p b] (inclusive) with tree - * nodes, writing them to @p out from left to right. + * @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& a, - const size_t& b, - size_t (&out)[64]) const noexcept { - if (leaf_count == 0 || a > b) { + 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 l = first_leaf_index + a; - size_t r = first_leaf_index + b; - size_t Lnodes[32]; - size_t Rnodes[32]; - size_t Ln = 0, Rn = 0; - while (l <= r) { - if (l & 1) { - Lnodes[Ln++] = l++; + 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 ((r & 1) == 0) { - Rnodes[Rn++] = r--; + if ((right_index & 1) == 0) { + right_nodes[right_count++] = right_index--; } - l >>= 1; - r >>= 1; + left_index >>= 1; + right_index >>= 1; } - size_t ind = 0; - for (size_t i = 0; i < Ln; ++i) { - out[ind++] = Lnodes[i]; + size_t out_count = 0; + for (size_t i = 0; i < left_count; ++i) { + out_nodes[out_count++] = left_nodes[i]; } - while (Rn > 0) { - out[ind++] = Rnodes[--Rn]; + while (right_count > 0) { + out_nodes[out_count++] = right_nodes[--right_count]; } - return ind; + return out_count; } /** - * @brief Select @p q-th minimum (1-based) inside [@p l..@p r] inclusive in - * two 8-bit passes. - * @details First pass finds global minimum, second selects the @p q-th + * @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; @@ -1930,160 +2076,165 @@ class RmMTree { /** * @brief Performs a forward search within a single leaf. - * @details Starting from position @p i (global index) inside the leaf [@p Lb - * .. Lb+len), it looks for the nearest position where the excess changes by - * @p d relative to @p i. 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 delta_leaf to the net - * excess change from @p i to the end of the 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 size_t leaf_fwd_bp_simd(const size_t& leaf_ind, - const size_t& Lb, - const size_t& i, - const int& d, - int& delta_leaf) const noexcept { - const size_t len = segment_size_bits[first_leaf_index + leaf_ind]; - const int16_t* lp = &leaf_prefix[leaf_ind * leaf_prefix_stride]; - const size_t pos = i - Lb; - if (pos > len) { - delta_leaf = 0; + 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_pref = lp[pos]; - size_t match_ind = ::find_forward_equal_i16_avx2(lp, pos + 1, len + 1, - start_pref + d, npos); - if (match_ind != npos) { - return Lb + match_ind - 1; + 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; } - delta_leaf = lp[len] - start_pref; + 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 i such that the excess difference from @p i equals @p d. + * 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 position @p i. + * and outputs the leaf-local excess delta at @p start_position. */ - inline size_t leaf_bwd_bp_simd(const size_t& leaf_ind, - const size_t& Lb, - const size_t& i, - const int& d, - int& delta_leaf) const noexcept { - const int16_t* lp = &leaf_prefix[leaf_ind * leaf_prefix_stride]; - const size_t pos = i - Lb; - if (pos > segment_size_bits[first_leaf_index + leaf_ind]) { - delta_leaf = 0; + 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_pref = lp[pos]; - if (pos > 0) { - const size_t match_ind = - ::find_backward_equal_i16_avx2(lp, 0, pos - 1, start_pref + d, npos); - if (match_ind != npos) { - return Lb + match_ind; + 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; } } - delta_leaf = start_pref; + leaf_delta = start_prefix; return npos; } /** - * @brief Find first minimum value and its first position in [@p l..@p 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 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 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 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 [@p l..@p 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); } /** @@ -2133,140 +2284,163 @@ class RmMTree { leaf_prefix.assign(leaf_count * leaf_prefix_stride, 0); // leaves - for (size_t k = 0; k < leaf_count; ++k) { - const size_t v = first_leaf_index + 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; - int16_t* lp = &leaf_prefix[k * leaf_prefix_stride]; - int16_t acc = 0; - lp[0] = 0; - for (size_t pos = Lb, i = 1; pos < Rb; ++pos, ++i) { - acc += bit(pos) ? 1 : -1; - lp[i] = acc; + 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 = first_leaf_index - 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; } } From decc9f3ced3ffd40d215eb37e7989a0a7a122917 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=D0=A0=D0=BE=D0=B1=D0=B5=D1=80=D1=82=20=D0=A1=D0=BC=D0=B0?= =?UTF-8?q?=D0=B9=D1=82?= Date: Sat, 29 Nov 2025 17:49:11 +0300 Subject: [PATCH 4/4] small comment refactoring --- include/bits.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/include/bits.h b/include/bits.h index a85b8fe..15f0e6c 100644 --- a/include/bits.h +++ b/include/bits.h @@ -414,7 +414,7 @@ void rank_32x8(const uint8_t* x, uint8_t* result) { * 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 npos if the value is not found. + * @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, @@ -456,7 +456,7 @@ static inline size_t find_forward_equal_i16_avx2(const int16_t* arr, * 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 `npos` if no match is found. + * 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,