matsutaku-library
This documentation is automatically generated by competitive-verifier/competitive-verifier
test/yosupo/rectangle_sum.test.cpp
- View this file on GitHub
- Last update: 2024-07-25 03:39:00+09:00
- Problem: https://judge.yosupo.jp/problem/rectangle_sum
Code
#define PROBLEM "https://judge.yosupo.jp/problem/rectangle_sum"
// #define IGNORE
#include "include/mtl/ordinal_range_search.hpp"
#include "include/mtl/fenwick_tree.hpp"
#include <iostream>
#include <bitset>
#include <bits/stdc++.h>
using namespace std;
int main() {
int n,q; cin>>n>>q;
vector<array<int,3>> N(n);
int max_idx = 0;
auto chmax = [&](int& a, int b) { a = max(a, b); };
for (int i = 0; i < n; i++) {
int x,y,w; cin>>x>>y>>w;
N[i] = {x,y,w};
chmax(max_idx, max(x, y));
}
vector<array<int,4>> Q(q);
for (int i = 0; i < q; i++) {
int l,d,r,u; cin>>l>>d>>r>>u;
Q[i] = {l,d,r,u};
chmax(max_idx, max(r, u));
}
ORS<int, long long> ors;
for (auto [x,y,w] : N) ors.add({x, y, w});
FenwickTree<long long> rsq;
ors.build(
[&](size_t size) {
rsq = FenwickTree<long long>(size);
},
[&](size_t i, int w) {
rsq.add(i, w);
});
for (auto [l,d,r,u] : Q) {
cout << ors.sum(l, r, d, u, [&](size_t l, size_t r) { return rsq.sum(l, r); }) << endl;
}
}#line 1 "test/yosupo/rectangle_sum.test.cpp"
#define PROBLEM "https://judge.yosupo.jp/problem/rectangle_sum"
// #define IGNORE
#line 2 "include/mtl/bit_manip.hpp"
#include <cstdint>
#include <cassert>
#if __cplusplus >= 202002L
#ifndef MTL_CPP20
#define MTL_CPP20
#endif
#include <bit>
#endif
namespace bm {
/// Count 1s for each 8 bits
inline constexpr uint64_t popcnt_e8(uint64_t x) {
x = (x & 0x5555555555555555) + ((x>>1) & 0x5555555555555555);
x = (x & 0x3333333333333333) + ((x>>2) & 0x3333333333333333);
x = (x & 0x0F0F0F0F0F0F0F0F) + ((x>>4) & 0x0F0F0F0F0F0F0F0F);
return x;
}
/// Count 1s
inline constexpr unsigned popcnt(uint64_t x) {
#ifdef MTL_CPP20
return std::popcount(x);
#else
return (popcnt_e8(x) * 0x0101010101010101) >> 56;
#endif
}
/// Alias to mtl::popcnt(x)
constexpr unsigned popcount(uint64_t x) {
return popcnt(x);
}
/// Count trailing 0s. s.t. *11011000 -> 3
inline constexpr unsigned ctz(uint64_t x) {
#ifdef MTL_CPP20
return std::countr_zero(x);
#else
return popcnt((x & (-x)) - 1);
#endif
}
/// Alias to mtl::ctz(x)
constexpr unsigned countr_zero(uint64_t x) {
return ctz(x);
}
/// Count trailing 1s. s.t. *11011011 -> 2
inline constexpr unsigned cto(uint64_t x) {
#ifdef MTL_CPP20
return std::countr_one(x);
#else
return ctz(~x);
#endif
}
/// Alias to mtl::cto(x)
constexpr unsigned countr_one(uint64_t x) {
return cto(x);
}
inline constexpr unsigned ctz8(uint8_t x) {
return x == 0 ? 8 : popcnt_e8((x & (-x)) - 1);
}
/// [00..0](8bit) -> 0, [**..*](not only 0) -> 1
inline constexpr uint8_t summary(uint64_t x) {
constexpr uint64_t hmask = 0x8080808080808080ull;
constexpr uint64_t lmask = 0x7F7F7F7F7F7F7F7Full;
auto a = x & hmask;
auto b = x & lmask;
b = hmask - b;
b = ~b;
auto c = (a | b) & hmask;
c *= 0x0002040810204081ull;
return uint8_t(c >> 56);
}
/// Extract target area of bits
inline constexpr uint64_t bextr(uint64_t x, unsigned start, unsigned len) {
uint64_t mask = len < 64 ? (1ull<<len)-1 : 0xFFFFFFFFFFFFFFFFull;
return (x >> start) & mask;
}
/// 00101101 -> 00111111 -count_1s-> 6
inline constexpr unsigned log2p1(uint8_t x) {
if (x & 0x80)
return 8;
uint64_t p = uint64_t(x) * 0x0101010101010101ull;
p -= 0x8040201008040201ull;
p = ~p & 0x8080808080808080ull;
p = (p >> 7) * 0x0101010101010101ull;
p >>= 56;
return p;
}
/// 00101100 -mask_mssb-> 00100000 -to_index-> 5
inline constexpr unsigned mssb8(uint8_t x) {
assert(x != 0);
return log2p1(x) - 1;
}
/// 00101100 -mask_lssb-> 00000100 -to_index-> 2
inline constexpr unsigned lssb8(uint8_t x) {
assert(x != 0);
return popcnt_e8((x & -x) - 1);
}
/// Count leading 0s. 00001011... -> 4
inline constexpr unsigned clz(uint64_t x) {
#ifdef MTL_CPP20
return std::countl_zero(x);
#else
if (x == 0)
return 64;
auto i = mssb8(summary(x));
auto j = mssb8(bextr(x, 8 * i, 8));
return 63 - (8 * i + j);
#endif
}
/// Alias to mtl::clz(x)
constexpr unsigned countl_zero(uint64_t x) {
return clz(x);
}
/// Count leading 1s. 11110100... -> 4
inline constexpr unsigned clo(uint64_t x) {
#ifdef MTL_CPP20
return std::countl_one(x);
#else
return clz(~x);
#endif
}
/// Alias to mtl::clo(x)
constexpr unsigned countl_one(uint64_t x) {
return clo(x);
}
inline constexpr unsigned clz8(uint8_t x) {
return x == 0 ? 8 : 7 - mssb8(x);
}
inline constexpr uint64_t bit_reverse(uint64_t x) {
x = ((x & 0x00000000FFFFFFFF) << 32) | ((x & 0xFFFFFFFF00000000) >> 32);
x = ((x & 0x0000FFFF0000FFFF) << 16) | ((x & 0xFFFF0000FFFF0000) >> 16);
x = ((x & 0x00FF00FF00FF00FF) << 8) | ((x & 0xFF00FF00FF00FF00) >> 8);
x = ((x & 0x0F0F0F0F0F0F0F0F) << 4) | ((x & 0xF0F0F0F0F0F0F0F0) >> 4);
x = ((x & 0x3333333333333333) << 2) | ((x & 0xCCCCCCCCCCCCCCCC) >> 2);
x = ((x & 0x5555555555555555) << 1) | ((x & 0xAAAAAAAAAAAAAAAA) >> 1);
return x;
}
/// Check if x is power of 2. 00100000 -> true, 00100001 -> false
constexpr bool has_single_bit(uint64_t x) noexcept {
#ifdef MTL_CPP20
return std::has_single_bit(x);
#else
return x != 0 && (x & (x - 1)) == 0;
#endif
}
/// Bit width needs to represent x. 00110110 -> 6
constexpr int bit_width(uint64_t x) noexcept {
#ifdef MTL_CPP20
return std::bit_width(x);
#else
return 64 - clz(x);
#endif
}
/// Ceil power of 2. 00110110 -> 01000000
constexpr uint64_t bit_ceil(uint64_t x) {
#ifdef MTL_CPP20
return std::bit_ceil(x);
#else
if (x == 0) return 1;
return 1ull << bit_width(x - 1);
#endif
}
/// Floor power of 2. 00110110 -> 00100000
constexpr uint64_t bit_floor(uint64_t x) {
#ifdef MTL_CPP20
return std::bit_floor(x);
#else
if (x == 0) return 0;
return 1ull << (bit_width(x) - 1);
#endif
}
} // namespace bm
#line 3 "include/mtl/fenwick_tree.hpp"
#include <cstddef>
#include <vector>
template <class T>
class FenwickTree {
private:
std::vector<T> tree_;
public:
FenwickTree() = default;
explicit FenwickTree(size_t size) : tree_(size+1) {}
size_t size() const { return tree_.size()-1; }
template <class Iter>
explicit FenwickTree(Iter begin, Iter end) : FenwickTree(std::distance(begin, end)) {
size_t i = 1;
for (auto it = begin; it != end; ++it) {
tree_[i] = tree_[i] + *it;
auto j = i + (i&(-i));
if (j < tree_.size())
tree_[j] = tree_[j] + tree_[i];
++i;
}
}
template<class V>
void add(size_t index, const V& x) {
for (size_t i = index+1; i < tree_.size(); i += i&(-i))
tree_[i] = tree_[i] + x;
}
T sum(size_t index) const {
T sum = 0;
for (size_t i = index+1; i > 0; i -= i&(-i))
sum = sum + tree_[i];
return sum;
}
T range_sum(size_t l, size_t r) const {
auto sl = l > 0 ? sum(l-1) : 0;
auto sr = r > 0 ? sum(r-1) : 0;
return sr - sl;
}
/// @brief Alias of range_sum(l, r)
T sum(size_t l, size_t r) const {
return range_sum(l, r);
}
template<class V>
size_t lower_bound(const V& _sum) const {
size_t ret = 0;
T s = 0;
for (int k = 63-bm::clz(size()); k >= 0; k--) {
size_t j = ret | (1ull<<k);
if (j < tree_.size() and s + tree_[j] < _sum) {
s = s + tree_[j];
ret = j;
}
}
return ret;
}
};
#line 3 "include/mtl/succinct/select.hpp"
#include <array>
struct select64 {
using size_type = uint8_t;
static struct make_select_table {
using bitmap_type = uint8_t;
std::array<std::array<size_type, 9>, 1<<8> tb;
make_select_table() {
for (int i = 0; i < 1<<8; i++) {
int c = 0;
int x = i;
tb[i].fill(8);
for (int j = 0; j < 8; j++) {
if (x & 1)
tb[i][++c] = j;
x >>= 1;
}
}
}
size_type operator()(bitmap_type bitmap, size_type ith) const {
return tb[bitmap][ith];
}
} sel_tb;
template<bool B>
static constexpr size_type select(size_type ith, uint64_t bitmap) { // 0-indexed
assert(ith < 64);
ith++; // to 1-index
// Find byte
uint64_t w = bitmap;
if constexpr (!B) w = ~w;
auto _bs = (uint64_t) bm::popcnt_e8(w) * 0x0101010101010100ull;
auto bs = (const uint8_t*) &_bs;
size_type b = 0;
auto o = ith;
/* Following bit-manipulates code is same as ... */
// {
// auto d = 8;
// while (d > 1) {
// auto c = b + d/2;
// if (bs[c] < o)
// b = c;
// d /= 2;
// }
// }
{
uint64_t x = (uint64_t) o * 0x0101010101010101ull;
uint64_t bmx = (_bs | 0x8080808080808080ull) - x;
uint64_t y = ((bmx & 0x8080808080808080ull) * 0x02040810204081ull) >> (64-8);
size_type nb = bm::ctz8(y) - 1;
// assert(b == nb);
b = nb;
}
// Calc select
auto r = o - bs[b];
uint8_t byte = ((const uint8_t*)(&w))[b];
assert(r and r <= (size_type)bm::popcnt(byte));
return b * 8 + sel_tb(byte, r);
}
static constexpr size_type select1(size_type ith, uint64_t bitmap) {
return select<1>(ith, bitmap);
}
static constexpr size_type select0(size_type ith, uint64_t bitmap) {
return select<0>(ith, bitmap);
}
};
typename select64::make_select_table select64::sel_tb;
#line 6 "include/mtl/succinct/bv.hpp"
#include <bitset>
#include <iostream>
#if __cpp_concepts >= 202002L
#include <concepts>
template<class T>
concept ConstructableBV = requires(T t, size_t s) {
{ t.size() } -> std::same_as<size_t>;
{ t.word_size() } -> std::same_as<size_t>;
{ t.get_word(s) } -> std::convertible_to<uint64_t>;
};
#endif
template<
#if __cpp_concepts >= 202002L
ConstructableBV
#else
class
#endif
BitVec, unsigned WordSize>
struct BV {
static_assert(WordSize <= 64, "WordSize must be <= 64");
static constexpr unsigned S = WordSize;
static constexpr unsigned S_PER_L = 8;
static constexpr unsigned L = S*S_PER_L;
using bitvec_type = BitVec;
const bitvec_type* bm_;
std::vector<uint64_t> _r, _s, _zs;
BV() = default;
explicit BV(const bitvec_type* bm) {
build(bm);
}
void set_ptr(const bitvec_type* bm) {
bm_ = bm;
}
void build(const bitvec_type* bm) {
set_ptr(bm);
const auto num_l = (bm->size() + L-1) / L;
_r.assign((num_l + 1) * 2, 0);
_s.clear();
_s.push_back(0);
_zs.clear();
_zs.push_back(0);
uint64_t sum = 0;
for (size_t l = 0; l < num_l; ++l) {
auto psum = sum;
uint64_t sum_l = 0;
for (size_t s = 0; s < S_PER_L; ++s) {
if (l * S_PER_L + s < bm->word_size())
sum_l += bm::popcnt(bm->get_word(l * S_PER_L + s));
if (s < S_PER_L-1)
_r[l * 2 + 1] |= sum_l << ((7-(s+1)) * 9);
}
sum += sum_l;
_r[(l + 1) * 2] = sum;
if (psum / L != sum / L) {
_s.push_back(l);
}
if ((L*l - psum) / L != (L*(l+1) - sum) / L) {
_zs.push_back(l);
}
}
_s.push_back(num_l);
_zs.push_back(num_l);
}
template<bool B>
size_t get_l(size_t l) const {
auto b = _r[l*2];
return B ? b : l * L - b;
}
static constexpr size_t s_off(size_t s) {
return (7-s) * 9;
}
template<bool B>
size_t get_s(size_t l, size_t s) const {
auto b = (_r[l*2+1] >> s_off(s)) % (1ull<<9);
return B ? b : s * S - b;
}
uint64_t mask(size_t width) const {
return width == 64 ? ~0ull : (1ull << width) - 1;
}
size_t rank1(size_t i) const {
auto l = i / L;
auto s = i % L / S;
auto q = i / S;
auto r = i % S;
assert(bm_ != nullptr);
auto w = bm_->get_word(q) & mask(r);
return get_l<1>(l) +
get_s<1>(l, s) +
bm::popcnt(w);
}
size_t rank0(size_t i) const {
return i - rank1(i);
}
template<bool B>
size_t rank(size_t i) const {
if constexpr (B)
return rank1(i);
else
return rank0(i);
}
static struct l_block_cap_mask {
uint64_t mask;
constexpr l_block_cap_mask() : mask(0) {
for (unsigned i = 0; i < S_PER_L; i++) {
uint64_t cap = i * S;
mask |= cap << s_off(i);
}
}
} l_block_cap;
template<bool B>
size_t select(size_t ith) const {
auto n = ith+1; // to 1-indexed
if (n > rank<B>(bm_->size()))
return bm_->size();
// Find L block
auto& idx = B ? _s : _zs;
size_t l = idx[n / L];
{
auto r = idx[n / L + 1] + 1;
while (l+1 < r) {
auto c = l + (r-l)/2;
if (get_l<B>(c) < n)
l = c;
else
r = c;
}
}
// Find S block
size_t s = 0;
auto m = n - get_l<B>(l);
/* Following bit-manipulates code is same as ... */
// {
// auto d = 8;
// while (d > 1) {
// auto c = s + d/2;
// if (get_s<B>(l, c) < m)
// s = c;
// d /= 2;
// }
// }
{
uint64_t x = (uint64_t) (m-1) * 0x0040201008040201ull;
uint64_t a = _r[l*2+1];
if constexpr (!B)
a = l_block_cap.mask - a; // to 0s sum
uint64_t xda = x - a;
uint64_t sm = 0x4020100804020100ull;
uint64_t ok = (~x | a) & sm;
uint64_t ng = (~x & a) & sm;
uint64_t y = ((x ^ a ^ xda) & ok) | ng;
y = y * 0x0001010101010101ull >> (64-1-7);
auto id = bm::clz8(y)-1;
auto clo = bm::clz((~xda << 1 | 1) << (9*id));
auto ns = id + (clo ? (clo - 1) / 9 : 0);
s = ns;
}
// Calc select
auto o = m - get_s<B>(l, s);
uint64_t w = bm_->get_word(l * S_PER_L + s);
return l * L +
s * S +
select64::select<B>(o-1, w);
}
size_t select1(size_t ith) const {
return select<1>(ith);
}
size_t select0(size_t ith) const {
return select<0>(ith);
}
};
template<class BitVec, unsigned WordSize>
typename BV<BitVec, WordSize>::l_block_cap_mask BV<BitVec, WordSize>::l_block_cap;
#line 8 "include/mtl/succinct/bitmap.hpp"
#include <iterator>
auto t = std::iterator_traits<std::vector<bool>::iterator>::iterator_category();
/// Bitmap is likes std::vector<bool> with advanced operations
struct Bitmap {
using value_type = bool;
using W = uint64_t;
std::vector<W> arr;
size_t sz;
static constexpr size_t required_word_size(size_t n) {
return (n+63) / 64;
}
static constexpr W word_filled_by(bool bit) {
return bit ? 0xFFFFFFFFFFFFFFFF : 0;
}
explicit Bitmap(size_t n = 0, bool bit = false)
: arr(required_word_size(n), word_filled_by(bit)), sz(n) {}
template<typename It>
Bitmap(It begin, It end) : sz(0) {
using trait = std::iterator_traits<It>;
using iterator_category = typename trait::iterator_category;
static_assert(std::is_base_of<std::input_iterator_tag, iterator_category>::value, "");
static_assert(std::is_convertible<typename trait::value_type, bool>::value, "");
if constexpr (std::is_base_of<std::random_access_iterator_tag, iterator_category>::value) {
arr.reserve(required_word_size(std::distance(begin, end)));
}
for (auto it = begin; it != end; ++it)
push_back((bool)*it);
}
size_t size() const { return sz; }
bool empty() const { return size() == 0; }
void push_back(bool bit) {
auto r = sz % 64;
if (r == 0) {
arr.push_back((W)bit);
} else {
if (bit)
arr.back() |= 1ull << r;
else
arr.back() &= ~(1ull << r);
}
++sz;
}
void pop_back() {
--sz;
}
void resize(size_t new_size, bool bit=false) { // TODO: fix when bit = true
auto old_size = size();
sz = new_size;
if (new_size < old_size) {
return;
}
arr.resize(required_word_size(new_size), word_filled_by(bit));
auto r = old_size % 64;
if (r) {
W mask = (1ull << r) - 1;
if (bit)
arr[old_size / 64] |= ~mask;
else
arr[old_size / 64] &= mask;
}
}
void assign(size_t new_size, bool bit) {
sz = new_size;
arr.assign(required_word_size(new_size), word_filled_by(bit));
}
void reserve(size_t reserved_size) {
arr.reserve(required_word_size(reserved_size));
}
struct const_reference;
struct reference;
template<bool>
struct _iterator;
using const_iterator = _iterator<true>;
using iterator = _iterator<false>;
const_reference operator[](size_t i) const {
return const_reference(arr.data() + i / 64, 1ull << (i % 64));
}
reference operator[](size_t i) {
return {arr.data() + i / 64, 1ull << (i % 64)};
}
const_reference get(size_t i) const {
return operator[](i);
}
/**
* Usable without pre-set required size
*/
void set(size_t i, bool b) {
if (i >= size())
resize(i + 1);
operator[](i) = b;
}
/**
* No build process is needed
*/
void build() const {}
void move_or_build(Bitmap&& src) {
*this = std::move(src);
}
const_iterator begin() const { return const_iterator(arr.data(), 0); }
iterator begin() { return iterator(arr.data(), 0); }
const_iterator cbegin() const { return begin(); }
const_iterator end() const { return const_iterator(arr.data() + sz / 64, sz % 64); }
iterator end() { return iterator(arr.data() + sz / 64, sz % 64); }
const_iterator cend() const { return end(); }
template<bool Const>
struct reference_base {
using _pointer = typename std::conditional<Const, const W*, W*>::type;
using _iterator = typename std::conditional<Const, const_iterator, iterator>::type;
_pointer ptr;
W mask;
reference_base(_pointer ptr, W mask) : ptr(ptr), mask(mask) {}
reference_base(const reference_base&) = delete;
reference_base& operator=(const reference_base&) = delete;
reference_base(reference_base&&) noexcept = default;
reference_base& operator=(reference_base&&) noexcept = default;
inline operator bool() const {
return (*ptr & mask) != 0;
}
inline bool operator==(bool r) const {
return (bool) *this == r;
}
inline friend bool operator==(bool l, const reference_base& r) {
return r == l;
}
inline bool operator!=(bool r) const {
return (bool) *this != r;
}
inline friend bool operator!=(bool l, const reference_base& r) {
return r != l;
}
_iterator operator&() const {
return {ptr, bm::ctz(mask)};
}
std::ostream& operator<<(std::ostream& os) const {
return os << (bool) *this;
}
};
struct const_reference : public reference_base<true> {
using _base = reference_base<true>;
const_reference(_base::_pointer ptr, W mask) : _base(ptr, mask) {}
const_reference(const reference& rhs) : _base(rhs.ptr, rhs.mask) {}
};
struct reference : public reference_base<false> {
using _base = reference_base<false>;
reference(_base::_pointer ptr, W mask) : _base(ptr, mask) {}
inline reference& operator=(bool bit) {
if (bit)
*ptr |= mask;
else
*ptr &= ~mask;
return *this;
}
inline reference& operator|=(bool bit) {
if (bit)
*ptr |= mask;
return *this;
}
inline reference& operator&=(bool bit) {
if (!bit)
*ptr &= ~mask;
return *this;
}
template<bool C>
inline reference& operator=(const reference_base<C>& rhs) {
return *this = (bool) rhs;
}
reference(const reference& rhs) = delete;
reference(reference&& rhs) noexcept = default;
reference& operator=(reference&& rhs) noexcept = default;
std::istream& operator>>(std::istream& is) {
bool b;
is >> b;
operator=(b);
return is;
}
};
template<bool Const>
struct _iterator {
using iterator_category = std::random_access_iterator_tag;
using value_type = bool;
using difference_type = std::ptrdiff_t;
using pointer = iterator;
using reference = typename std::conditional<Const, const_reference, Bitmap::reference>::type;
using _pointer = typename std::conditional<Const, const W*, W*>::type;
_pointer ptr;
unsigned ctz;
_iterator(_pointer ptr, unsigned ctz) : ptr(ptr), ctz(ctz) {}
inline reference operator*() const {
return reference(ptr, 1ull << ctz);
}
template<bool C>
inline bool operator==(const _iterator<C>& r) const {
return ptr == r.ptr and ctz == r.ctz;
}
template<bool C>
inline bool operator!=(const _iterator<C>& r) const {
return ptr != r.ptr or ctz != r.ctz;
}
inline _iterator& operator++() {
if (ctz % 64 < 63) {
++ctz;
} else {
++ptr;
ctz = 0;
}
return *this;
}
inline _iterator operator++(int) {
_iterator ret = *this;
operator++();
return ret;
}
inline _iterator& operator--() {
if (ctz % 64 > 0) {
--ctz;
} else {
--ptr;
ctz = 63;
}
return *this;
}
inline _iterator operator--(int) {
_iterator ret = *this;
operator--();
return ret;
}
inline _iterator& operator+=(std::ptrdiff_t d) {
if (d < 0)
return operator-=(-d);
auto r = d % 64;
if (ctz + r < 64) {
ptr += d / 64;
ctz += r;
} else {
ptr += d / 64 + 1;
ctz = (ctz + r) - 64;
}
return *this;
}
inline _iterator operator+(std::ptrdiff_t d) const {
return _iterator(*this) += d;
}
inline _iterator& operator-=(std::ptrdiff_t d) {
if (d < 0)
return operator+=(-d);
auto r = d % 64;
if (r <= ctz) {
ptr -= d / 64;
ctz -= r;
} else {
ptr -= d / 64 + 1;
ctz = (ctz + 64 - r) - 64;
}
return *this;
}
inline _iterator operator-(std::ptrdiff_t d) const {
return _iterator(*this) -= d;
}
inline reference operator[](size_t i) const {
return *(*this + i);
}
};
void range_set(size_t b, size_t e, uint64_t x) {
if (b >= e) return;
auto r = b % 64;
auto w = e-b;
auto mask = w < 64 ? (1ull << w) - 1 : ~0ull;
assert(x <= mask);
arr[b/64] = (arr[b/64] & ~(mask << r)) | x << r;
if (mask + r > 64) {
arr[b/64+1] = (arr[b/64+1] & ~(mask >> (64-r))) | x >> (64-r);
}
}
uint64_t range_get(size_t b, size_t e) const {
if (b >= e) return 0;
assert(e-b <= 64);
auto r = b % 64;
auto w = e-b;
auto mask = w < 64 ? (1ull << w) - 1 : ~0ull;
auto x = arr[b/64] >> r;
if (w + r > 64)
x |= arr[b/64+1] << (64-r);
return x & mask;
}
const uint64_t& get_word(size_t wi) const {
return arr[wi];
}
size_t word_size() const {
return arr.size();
}
// using rank_select_type = BV<Bitmap, 64>;
};
#line 4 "include/mtl/succinct/ty.hpp"
#include <limits>
#line 6 "include/mtl/succinct/ty.hpp"
#include <algorithm>
#line 8 "include/mtl/succinct/ty.hpp"
/**
* @brief TY: Store increasing sequence of integers.
* Memory needs for store nth integers O(n log d) bits
* which d is max diff of consecutive elements.
*/
template<class T, class DiffType = int16_t>
class TY {
using value_type = T;
static constexpr auto block_size = sizeof(value_type) * 8;
using diff_value_type = DiffType;
static constexpr unsigned max_diff = std::numeric_limits<diff_value_type>::max();
private:
std::vector<value_type> head;
std::vector<diff_value_type> diff;
public:
TY() = default;
template<class It>
TY(It first, It last) {
assert(std::is_sorted(first, last));
reserve(std::distance(first, last));
for (auto it = first; it != last; it++) {
push_back(*it);
}
}
size_t size() const {
return head.size() + diff.size();
}
bool empty() const { return size() == 0; }
void reserve(size_t n) {
head.reserve((n + block_size - 1) / block_size);
diff.reserve(n / block_size * (block_size - 1) + n % block_size);
}
template<class... Args>
void emplace_back(Args&&... args) {
if (size() % block_size == 0) {
head.emplace_back(std::forward<Args>(args)...);
} else {
value_type v(std::forward<Args>(args)...);
assert(v >= head.back());
assert(v - head.back() <= (value_type)max_diff);
diff.push_back((diff_value_type)(v - head.back()));
}
}
void push_back(const value_type& v) {
if (size() % block_size == 0) {
head.push_back(v);
} else {
assert(v >= head.back());
assert(v - head.back() <= (value_type)max_diff);
diff.push_back(v - head.back());
}
}
void push_back(value_type&& v) {
emplace_back(std::move(v));
}
value_type get(size_t i) const {
if (i % block_size == 0)
return head[i / block_size];
else
return head[i / block_size] +
(value_type)diff[i / block_size * (block_size-1) + i % block_size - 1];
}
value_type operator[](size_t i) const { return get(i); }
value_type front() const { return get(0); }
value_type back() const { return get(size()-1); }
};
#line 7 "include/mtl/succinct/rrr.hpp"
#include <map>
#line 13 "include/mtl/succinct/rrr.hpp"
constexpr unsigned need_bits(uint64_t n) {
return bm::bit_width(n);
}
template<unsigned N>
struct BinomialTable {
static_assert(N < 64,
"Too large N for BinomialTable. N must be less than 64");
using number_type = uint64_t;
using binom_table_type = std::array<std::array<number_type, N+1>, N+1>;
static constexpr binom_table_type make_binomial_table() {
binom_table_type tb{};
tb[0][0] = 1;
for (size_t i = 1; i <= N; i++) {
tb[i][0] = tb[i-1][0];
for (size_t j = 1; j <= i; j++)
tb[i][j] = tb[i-1][j-1] + tb[i-1][j];
}
return tb;
}
static constexpr binom_table_type _binom_tb = make_binomial_table();
static constexpr number_type binom(size_t n, size_t k) {
assert(n <= N and k <= N);
return _binom_tb[n][k];
}
};
template<unsigned N>
constexpr typename BinomialTable<N>::binom_table_type BinomialTable<N>::_binom_tb;
template<class Def>
struct RRRTable {
using def = Def;
static constexpr unsigned s_size = def::s_size;
using s_type = typename def::s_type;
static constexpr unsigned n_bits = def::n_bits;
using binomial_table_type = BinomialTable<s_size>;
using number_type = typename binomial_table_type::number_type;
static constexpr s_type get_int(unsigned n, number_type k, unsigned bits = s_size) {
s_type res = 0;
const auto offset = bits;
s_type mask = ((s_type(1)<<bits)-1);
auto nn = n;
unsigned i = 0;
/*
Binary search time B = ceil(\log_2 w)
Expected length of consecutive zeros Ez = \sum j binom(w-j, nn) / binom(w, nn), j=1..w-nn
Expected length of consecutive ones Eo = \sum j binom(w-j, nn-j) / binom(w, nn), j=1..nn
Approximate simple function from Ez > B to be nn <= min(20, w-1)
Approximate simple function from Eo > B to be nn > min(40, w)
*/
// TODO: When nn > 40, use binary search to find length of consecutive ones
for (; i < offset and nn > 20; i++) {
auto w = s_size - i;
if (nn == w) {
res |= ((s_type(1)<<nn)-1) << i;
return res & mask;
}
if (nn == w-1) {
res |= (((s_type(1)<<w)-1) ^ (s_type(1) << k)) << i;
return res & mask;
}
// Linear search
auto binom = binomial_table_type::binom(w-1, nn);
if (k >= binom) {
res |= s_type(1) << i;
k -= binom;
nn--;
}
}
for (; i < offset and nn > 1; i++) {
auto w = s_size - i;
if (nn == w) {
res |= ((s_type(1)<<nn)-1) << i;
return res & mask;
}
if (nn == w-1) {
res |= (((s_type(1)<<w)-1) ^ (s_type(1) << k)) << i;
return res & mask;
}
// Binary search to find length of consecutive zeros
auto l = i, r = offset+1;
while (l+1<r) {
auto c = l+(r-l)/2;
if (k < binomial_table_type::binom(s_size-c, nn))
l = c;
else
r = c;
}
if (l < offset) {
res |= s_type(1) << l;
k -= binomial_table_type::binom(s_size-l-1, nn);
nn--;
}
i = l;
}
if (nn == 1) {
res |= s_type(1) << (s_size-1-k);
return res & mask;
}
if (nn == 0)
return res;
if (k >= binomial_table_type::binom(s_size-offset-1, nn))
res |= s_type(1) << offset;
return res & ((s_type(1)<<bits)-1);
}
static constexpr bool get_bit(unsigned n, number_type k, unsigned offset) {
auto nn = n;
unsigned i = 0;
/*
Binary search time B = ceil(\log_2 w)
Expected length of consecutive zeros Ez = \sum j binom(w-j, nn) / binom(w, nn), j=1..w-nn
Expected length of consecutive ones Eo = \sum j binom(w-j, nn-j) / binom(w, nn), j=1..nn
Approximate simple function from Ez > B to be nn <= min(20, w-1)
Approximate simple function from Eo > B to be nn > min(40, w)
*/
// TODO: When nn > 40, use binary search to find length of consecutive ones
for (; i < offset and nn > 20; i++) {
auto w = s_size - i;
if (nn == w) {
return 1;
}
if (nn == w-1) {
return offset != i+k;
}
// linear search
auto binom = binomial_table_type::binom(w-1, nn);
if (k >= binom) {
k -= binom;
nn--;
}
}
for (; i < offset and nn > 1; i++) {
auto w = s_size - i;
if (nn == w) {
return 1;
}
if (nn == w-1) {
return offset != i+k;
}
// binary search
auto l = i, r = offset+1;
while (l+1<r) {
auto c = l+(r-l)/2;
if (k < binomial_table_type::binom(s_size-c, nn))
l = c;
else
r = c;
}
if (l < offset) {
k -= binomial_table_type::binom(s_size-l-1, nn);
nn--;
}
i = l;
}
if (nn == 1)
return offset == s_size-1-k;
if (nn == 0)
return 0;
return k >= binomial_table_type::binom(s_size-offset-1, nn);
}
static constexpr number_type get_number_for_popcnt(s_type bitmap, unsigned pc) {
number_type number = 0;
auto m = bitmap;
auto n = pc;
while (m) {
auto i = bm::ctz(m);
number += binomial_table_type::binom(s_size-i-1, n);
n--;
m ^= (s_type(1)<<i);
}
return number;
}
static constexpr number_type number_size(unsigned n) {
return binomial_table_type::binom(s_size, n);
}
static constexpr std::pair<unsigned, number_type> get_pc_and_number(s_type bitmap) {
unsigned pc = bm::popcnt(bitmap);
auto number = pc <= s_size-pc ? get_number_for_popcnt(bitmap, pc)
: (number_size(pc)-1-get_number_for_popcnt(
~bitmap & ((s_type(1)<<s_size)-1), s_size-pc));
return std::make_pair(pc, number);
}
using number_bits_table_type = std::array<unsigned, s_size+1>;
static constexpr number_bits_table_type make_number_bits_table() {
number_bits_table_type tb{};
for (unsigned i = 0; i <= s_size; i++) {
tb[i] = need_bits(number_size(i)-1);
}
return tb;
}
static constexpr number_bits_table_type n_len = make_number_bits_table();
static constexpr unsigned number_bits(unsigned n) {
assert(n <= s_size);
return n_len[n];
}
};
template<class Def>
constexpr typename RRRTable<Def>::number_bits_table_type RRRTable<Def>::n_len;
template<unsigned SSize, class SType>
struct RRRDefinition {
static constexpr unsigned s_size = SSize;
using s_type = SType;
static constexpr unsigned n_bits = need_bits(s_size);
};
/**
* @brief Succinct bit vector in memory of B(n, u) + O(u log log n / log n) bits
* where u is number of bits and n is number of 1s
*/
template<
unsigned SSize = 63,
class SType = uint64_t,
class MapType = std::map<size_t, SType>
>
struct RRR {
using def = RRRDefinition<SSize, SType>;
using s_type = typename def::s_type;
using rrr_table_type = RRRTable<def>;
using map_type = MapType;
using ty_type = TY<size_t>;
map_type s_map;
ty_type heads;
Bitmap bm;
RRR() = default;
void set(size_t i, bool b) {
if (b)
s_map[i/def::s_size] |= (s_type)1<<(i%def::s_size);
else
s_map[i/def::s_size] &= ~((s_type)1<<(i%def::s_size));
}
void build() {
size_t h = 0;
size_t pq = 0;
auto block_count = s_map.empty() ? 0 : std::prev(s_map.end())->first+1;
heads.reserve(block_count);
for (auto qm : s_map) {
auto qidx = qm.first;
auto mask = qm.second;
while (pq < qidx) {
heads.push_back(h);
auto w = def::n_bits;
bm.resize(h+w);
bm.range_set(h, h+w, 0);
h += w;
pq++;
}
heads.push_back(h);
auto np = rrr_table_type::get_pc_and_number(mask);
auto n = np.first;
auto p = np.second;
assert(rrr_table_type::get_int(n, p) == mask);
auto w = def::n_bits + rrr_table_type::number_bits(n);
bm.resize(h+w);
bm.range_set(h, h+def::n_bits, n);
bm.range_set(h+def::n_bits, h+w, p);
assert(bm.range_get(h, h+def::n_bits) == n);
assert(bm.range_get(h+def::n_bits, h+w) == p);
h += w;
pq++;
}
s_map.clear();
}
void move_or_build(RRR&& src) {
*this = std::move(src);
}
void move_or_build(const Bitmap& bm) {
for (size_t i = 0; i < bm.size(); i += def::s_size) {
auto w = bm.range_get(i, std::min(i+def::s_size, bm.size()));
if (w or i+def::s_size >= bm.size()) s_map.emplace(i/def::s_size, w);
}
build();
}
bool get_bit(size_t si, unsigned off) const {
if (si >= heads.size())
return false;
auto a = heads.get(si);
auto b = a+def::n_bits;
auto n = bm.range_get(a, b);
auto p = bm.range_get(b, b+rrr_table_type::number_bits(n));
return rrr_table_type::get_bit(n, p, off);
}
s_type get_mask(size_t si) const {
if (si >= heads.size())
return 0;
auto a = heads.get(si);
auto b = a+def::n_bits;
auto n = bm.range_get(a, b);
auto p = bm.range_get(b, b+rrr_table_type::number_bits(n));
return rrr_table_type::get_int(n, p);
}
uint64_t get_word(size_t si) const {
return get_mask(si);
}
size_t word_size() const {
return heads.size();
}
size_t size() const {
return heads.size() * def::s_size;
}
bool empty() const {
return size() == 0;
}
bool get(size_t i) const {
return get_bit(i/def::s_size, i%def::s_size);
}
};
#line 5 "include/mtl/succinct/traits.hpp"
template<class T>
struct RankSelectTraits : std::false_type {};
template<>
struct RankSelectTraits<Bitmap> {
using rank_select_type = BV<Bitmap, 64>;
};
template<unsigned SSize, class SType, class MapType>
struct RankSelectTraits<RRR<SSize, SType, MapType>> {
using rank_select_type = BV<RRR<SSize, SType, MapType>, SSize>;
};
#line 10 "include/mtl/succinct/bit_vector.hpp"
struct BitVector {
using bitmap_type = Bitmap;
bitmap_type bm;
using rs_type = typename RankSelectTraits<bitmap_type>::rank_select_type;
rs_type rs_support;
// std::vector<uint64_t> _r, _s, _zs;
BitVector() = default;
explicit BitVector(size_t size) : bm(size) {}
BitVector(size_t size, bool bit) : bm(size, bit) {}
BitVector(const BitVector& rhs) : bm(rhs.bm), rs_support(rhs.rs_support) {
rs_support.set_ptr(&bm);
}
BitVector& operator=(const BitVector& rhs) {
bm = rhs.bm;
rs_support = rhs.rs_support;
rs_support.set_ptr(&bm);
return *this;
}
BitVector(BitVector&& rhs) noexcept :
bm(std::move(rhs.bm)),
rs_support(std::move(rhs.rs_support)) {
rs_support.set_ptr(&bm);
}
BitVector& operator=(BitVector&& rhs) noexcept {
bm = std::move(rhs.bm);
rs_support = std::move(rhs.rs_support);
rs_support.set_ptr(&bm);
return *this;
}
template<typename It>
BitVector(It begin, It end) : bm(begin, end) {
build();
}
size_t size() const { return bm.size(); }
bool empty() const { return bm.empty(); }
void push_back(bool bit) { bm.push_back(bit); }
void resize(size_t new_size, bool bit = false) { bm.resize(new_size, bit); }
void assign(size_t new_size, bool bit) { bm.assign(new_size, bit); }
void reserve(size_t reserved_size) { bm.reserve(reserved_size); }
bitmap_type& bitmap() { return bm; }
const bitmap_type& bitmap() const { return bm; }
bitmap_type::const_reference operator[](size_t i) const { return bm[i]; }
bitmap_type::reference operator[](size_t i) { return bm[i]; }
bitmap_type::const_iterator begin() const { return bm.begin(); }
bitmap_type::const_iterator end() const { return bm.end(); }
bitmap_type::iterator begin() { return bm.begin(); }
bitmap_type::iterator end() { return bm.end(); }
bitmap_type::const_iterator cbegin() const { return bm.cbegin(); }
bitmap_type::const_iterator cend() const { return bm.cend(); }
void build() {
rs_support.build(&bm);
}
inline size_t rank(size_t i) const {
return rs_support.rank1(i);
}
size_t rank1(size_t i) const {
return rs_support.rank1(i);
}
size_t rank0(size_t i) const {
return rs_support.rank0(i);
}
template<bool B>
size_t select(size_t n) const {
return rs_support.select<B>(n);
}
size_t select1(size_t n) const {
return rs_support.select<1>(n);
}
size_t select0(size_t n) const {
return rs_support.select<0>(n);
}
};
#line 9 "include/mtl/succinct/wavelet_matrix.hpp"
#include <queue>
#include <tuple>
template<class T, class BitmapType = Bitmap, unsigned Height = 0>
struct WaveletMatrix {
static constexpr unsigned H = 64 - bm::clz(std::numeric_limits<T>::max());
size_t n,h;
using bitmap_type = BitmapType;
bitmap_type B;
using rs_type = typename RankSelectTraits<bitmap_type>::rank_select_type;
rs_type rs_b;
std::vector<size_t> RO,Z;
WaveletMatrix() = default;
template<typename It>
WaveletMatrix(It begin, It end)
: n(std::distance(begin, end)),
h(Height == 0 ? std::max(1u, 64 - bm::clz(n ? *max_element(begin, end) : 0)) : H),
B(),
rs_b(),
RO(h+1),
Z(h)
{
using trait = std::iterator_traits<It>;
static_assert(std::is_base_of<std::input_iterator_tag, typename trait::iterator_category>::value, "");
static_assert(std::is_convertible<typename trait::value_type, T>::value, "");
if (n == 0) return;
assert(*min_element(begin, end) >= 0);
std::vector<T> S(begin, end), z, o;
z.reserve(n);
o.reserve(n);
Bitmap _B(n*h, false);
auto bit = _B.begin();
for (int k = h-1; k >= 0; k--) {
for (size_t i = 0; i < n; i++) {
bool b = S[i] >> k & 1;
*bit++ = b;
if (!b)
z.push_back(S[i]);
else
o.push_back(S[i]);
}
Z[k] = n*(h-1-k+1) + z.size();
auto j = n;
while (!o.empty()) {
S[--j] = o.back();
o.pop_back();
}
while (!z.empty()) {
S[--j] = z.back();
z.pop_back();
}
assert(j == 0);
}
B.move_or_build(std::move(_B));
rs_b.build(&B);
for (size_t i = 0; i <= h; i++)
RO[i] = rs_b.rank1(n * i);
}
WaveletMatrix(const WaveletMatrix& rhs) noexcept :
n(rhs.n),
h(rhs.h),
B(rhs.B),
rs_b(rhs.rs_b),
RO(rhs.RO),
Z(rhs.Z)
{
rs_b.set_ptr(&B);
}
WaveletMatrix& operator=(const WaveletMatrix& rhs) noexcept {
n = rhs.n;
h = rhs.h;
B = rhs.B;
rs_b = rhs.rs_b;
rs_b.set_ptr(&B);
RO = rhs.RO;
Z = rhs.Z;
return *this;
}
WaveletMatrix(WaveletMatrix&& rhs) noexcept :
n(std::move(rhs.n)),
h(std::move(rhs.h)),
B(std::move(rhs.B)),
rs_b(std::move(rhs.rs_b)),
RO(std::move(rhs.RO)),
Z(std::move(rhs.Z))
{
rs_b.set_ptr(&B);
}
WaveletMatrix& operator=(WaveletMatrix&& rhs) noexcept {
n = std::move(rhs.n);
h = std::move(rhs.h);
B = std::move(rhs.B);
rs_b = std::move(rhs.rs_b);
rs_b.set_ptr(&B);
RO = std::move(rhs.RO);
Z = std::move(rhs.Z);
return *this;
}
inline size_t _child0(size_t level, size_t i) const {
return i + n + RO[level] - rs_b.rank1(i);
}
inline size_t _child1(size_t level, size_t i) const {
return n*(level+2) + rs_b.rank1(i) - RO[level+1];
}
inline size_t child(size_t level, size_t i, bool bit) const {
return !bit ? _child0(level, i) : _child1(level, i);
}
std::pair<size_t, size_t> _child_tie0(size_t level, size_t l, size_t r) const {
return std::make_pair(_child0(level, l), _child0(level, r));
}
std::pair<size_t, size_t> _child_tie1(size_t level, size_t l, size_t r) const {
return std::make_pair(_child1(level, l), _child1(level, r));
}
std::pair<size_t, size_t> child_tie(size_t level, size_t l, size_t r, bool bit) const {
return !bit ? _child_tie0(level, l, r) : _child_tie1(level, l, r);
}
inline size_t _parent0(size_t level, size_t i) const {
return rs_b.select0(i - n - RO[level]);
}
inline size_t _parent1(size_t level, size_t i) const {
return rs_b.select1(i - Z[h-1-level] + RO[level]);
}
inline size_t parent(size_t level, size_t i, bool bit) const {
return !bit ? _parent0(level, i) : _parent1(level, i);
}
T get(size_t i) const {
T c = 0;
size_t j = i;
for (int k = h-1; k > 0; k--) {
bool b = B[j];
j = child(h-1-k, j, b);
if (b)
c |= 1ull<<k;
}
if (B[j])
c |= 1u;
return c;
}
size_t range_rank(T c, size_t l, size_t r) const {
for (int k = h-1; k >= 0 and l < r; k--) {
std::tie(l,r) = child_tie(h-1-k, l, r, (c >> k) & 1u);
}
return r - l;
}
size_t rank(T c, size_t i) const {
return range_rank(c, 0, i);
}
std::tuple<size_t, size_t, size_t> rank_3way(T c, size_t l, size_t r) const {
size_t lt = 0, gt = 0;
for (int k = h-1; k >= 0; k--) {
size_t pr = r - l;
if (pr == 0)
break;
if (((c >> k) & 1u) == 0) {
std::tie(l,r) = _child_tie0(h-1-k, l, r);
gt += pr - (r - l);
} else {
std::tie(l,r) = _child_tie1(h-1-k, l, r);
lt += pr - (r - l);
}
}
return std::make_tuple(lt, r - l, gt);
}
/// Get frequency of values which (x <= value < y) in S[l,r).
size_t range_freq(size_t l, size_t r, T x, T y) const {
if (l == r) return 0;
size_t freq = 0;
std::queue<std::tuple<size_t,size_t, T>> qs;
qs.emplace(l, r, T(0));
while (!qs.empty()) {
size_t _l,_r;
T c;
std::tie(_l,_r,c) = qs.front();
qs.pop();
assert(_l < _r);
size_t level = _l/n;
int shift = h-1-level;
T clo = c, chi = c | ((1ull<<(shift+1))-1);
if (chi < x or y <= clo)
continue;
if (x <= clo and chi < y) {
freq += _r - _l;
continue;
}
assert(level < h);
size_t nl,nr;
std::tie(nl,nr) = child_tie(level, _l, _r, 0);
if (nl < nr)
qs.emplace(nl, nr, c);
std::tie(nl,nr) = child_tie(level, _l, _r, 1);
if (nl < nr)
qs.emplace(nl, nr, c | (1ull << shift));
}
return freq;
}
size_t range_select(T c, size_t l, size_t r, size_t i) const {
if (r - l <= i)
return n;
for (int k = h-1; k >= 0; k--) {
std::tie(l,r) = child_tie(h-1-k, l, r, (c >> k) & 1u);
if (r - l <= i)
return n;
}
size_t j = l+i;
for (size_t k = 0; k < h; k++) {
j = parent(h-1-k, j, (c >> k) & 1u);
assert((bool)((c>>k)&1u) == B[j]);
}
return j;
}
size_t select(T c, size_t i) const {
return range_select(c, 0, n, i);
}
/// Get kth (0-indexed) smallest value in S[l,r).
T quantile(size_t l, size_t r, size_t k) const {
assert(r - l > k);
T c = 0;
for (int d = h-1; d > 0; d--) {
auto os = rs_b.rank1(r) - rs_b.rank1(l);
auto zs = r - l - os;
if (k < zs) {
std::tie(l,r) = _child_tie0(h-1-d, l, r);
} else {
c |= 1ull << d;
k -= zs;
std::tie(l,r) = _child_tie1(h-1-d, l, r);
}
assert(l < r);
}
auto os = rs_b.rank1(r) - rs_b.rank1(l);
auto zs = r - l - os;
if (k >= zs) {
c |= 1ull;
}
return c;
}
/// Get tuples (value, frequency) of the most k frequently occurring values in S[l,r).
std::vector<std::pair<T, size_t>> top_k(size_t l, size_t r, size_t k) const {
std::vector<std::pair<T, size_t>> ret;
std::priority_queue<std::tuple<size_t, size_t, T>> qs;
qs.emplace(r-l, l, 0);
while (!qs.empty()) {
size_t range, s;
T c;
std::tie(range, s, c) = qs.top();
qs.pop();
auto level = s/n;
if (level == h) {
ret.emplace_back(c, range);
if (ret.size() == k)
break;
} else {
size_t _l, _r;
for (int b = 0; b < 2; b++) {
std::tie(_l,_r) = child_tie(level, s, s+range, b);
if (_l != _r)
qs.emplace(_r-_l, _l, c | (uint64_t(b) << (h-1-level)));
}
}
}
return ret;
}
/// Get sum of S[l,r) in O(min(r-l, \sigma) log \sigma) times.
template<typename U=T>
U sum(size_t l, size_t r) const {
U ret = 0;
for (auto p : top_k(l, r, r-l))
ret += (U) p.first * p.second;
return ret;
}
/// Get k tuples of (value, frequency) that value satisfies condition (x <= value < y) in S[l,r) from the smallest (or largest).
/// The complexity is O(k log \sigma).
template<bool ASCENDING, bool VALUE_RANGE = true>
std::vector<std::pair<T, size_t>>
range_list_k(size_t l, size_t r, size_t k, T x, T y) const {
std::vector<std::pair<T, size_t>> ret;
std::queue<std::tuple<size_t, size_t, T>> qs;
qs.emplace(r-l, l, T(0));
size_t range,s;
T c;
while (!qs.empty()) {
std::tie(range,s,c) = qs.top();
qs.pop();
auto level = s/n;
if (level == h) {
assert(!VALUE_RANGE or (x <= c and c < y));
ret.emplace_back(c, range);
if (ret.size() == k)
break;
} else {
auto shift = (h-1-level);
for (int b = ASCENDING ? 0 : 1;
ASCENDING ? b < 2 : b >= 0;
ASCENDING ? b++ : b--) {
T nc = (c << 1) | b;
if (VALUE_RANGE and (nc < (x >> shift) or ((y-1) >> shift) < nc))
continue;
size_t nl,nr;
std::tie(nl,nr) = child_tie(level, s, s+range, b);
if (nl != nr)
qs.emplace(nr-nl, nl, nc);
}
}
}
return ret;
}
/// Get tuples of (value, frequency) that value satisfies condition (x <= value < y) in S[l,r).
/// The complexity is O(k log \sigma).
std::vector<std::pair<T, size_t>> range_list(size_t l, size_t r, T x, T y) const {
return range_list_k<true>(l, r, r - l, x, y);
}
/// Get k tuples of (value, frequency) that value satisfies condition (x <= value < y) in S[l,r) from the largest.
/// The complexity is O(k log \sigma).
std::vector<std::pair<T, size_t>> range_max_k(size_t l, size_t r, size_t k) const {
return range_list_k<false, false>(l, r, k, 0, 0);
}
/// Get k tuples of (value, frequency) that value satisfies condition (x <= value < y) in S[l,r) from the smallest.
// The complexity is O(k log \sigma).
std::vector<std::pair<T, size_t>> range_min_k(size_t l, size_t r, size_t k) const {
return range_list_k<true, false>(l, r, k, 0, 0);
}
/// Get tuples (value, frequency of T1, frequency of T2) that commonly occur between T1=S[l1,r1) and T2=S[l2,r2).
std::vector<std::tuple<T, size_t, size_t>> intersect(size_t l1, size_t r1, size_t l2, size_t r2) const {
std::vector<std::tuple<T, size_t, size_t>> ret;
std::queue<std::pair<std::array<size_t,4>, T>> qs;
qs.emplace({l1,r1,l2,r2}, 0);
std::array<size_t,4> nrs{};
while (!qs.empty()) {
const auto& rs = qs.front().first;
T c = qs.front().second;
auto level = rs[0]/n;
if (level == h) {
ret.emplace_back(c, rs[1]-rs[0], rs[3]-rs[2]);
} else {
for (int b = 0; b < 2; b++) {
for (int i = 0; i < 4; i++)
nrs[i] = child(level, rs[i], b);
if (nrs[0] != nrs[1] and nrs[2] != nrs[3])
qs.emplace(nrs, c | (uint64_t(b) << (h-1-level)));
}
}
qs.pop();
}
return ret;
}
};
#line 8 "include/mtl/succinct/binary_set.hpp"
template<class BitmapType = RRR<>>
class BinarySet {
public:
using value_type = uint32_t;
private:
using bitmap_type = BitmapType;
using rs_type = typename RankSelectTraits<bitmap_type>::rank_select_type;
bitmap_type bm_;
rs_type rs_;
size_t size_;
public:
BinarySet() = default;
BinarySet(const BinarySet& rhs) : bm_(rhs.bm_), rs_(rhs.rs_), size_(rhs.size_) {
rs_.set_ptr(&bm_);
}
BinarySet& operator=(const BinarySet& rhs) {
bm_ = rhs.bm_;
rs_ = rhs.rs_;
size_ = rhs.size_;
rs_.set_ptr(&bm_);
return *this;
}
BinarySet(BinarySet&& rhs) : bm_(std::move(rhs.bm_)), rs_(std::move(rhs.rs_)), size_(std::move(rhs.size_)) {
rs_.set_ptr(&bm_);
}
BinarySet& operator=(BinarySet&& rhs) {
bm_ = std::move(rhs.bm_);
rs_ = std::move(rhs.rs_);
size_ = std::move(rhs.size_);
rs_.set_ptr(&bm_);
return *this;
}
~BinarySet() = default;
void insert(value_type x) {
bm_.set(x, 1);
}
void build() {
bm_.build();
rs_.build(&bm_);
size_ = rs_.rank1(bm_.size());
}
size_t size() const { return size_; }
bool empty() const { return size() == 0;}
bool contains(value_type x) const {
return bm_.get(x);
}
value_type lower_bound(value_type x) const {
if (contains(x)) return x;
return get(rank(x));
}
value_type upper_bound(value_type x) const {
return lower_bound(x+1);
}
size_t rank(value_type x) const {
return x < bm_.size() ? rs_.rank1(x) : size_;
}
value_type get(size_t i) const {
return i < size_ ? rs_.select1(i) : bm_.size();
}
};
template<class BitmapType = RRR<>>
class BinaryMultiset {
public:
using value_type = uint32_t;
static constexpr value_type ValueMax = std::numeric_limits<value_type>::max();
private:
using bitmap_type = BitmapType;
using rs_type = typename RankSelectTraits<bitmap_type>::rank_select_type;
bitmap_type bm_;
rs_type rs_;
size_t size_;
std::vector<value_type> values_;
public:
BinaryMultiset() = default;
BinaryMultiset(const BinaryMultiset& rhs) : bm_(rhs.bm_), rs_(rhs.rs_), size_(rhs.size_), values_(rhs.values_) {
rs_.set_ptr(&bm_);
}
BinaryMultiset& operator=(const BinaryMultiset& rhs) {
bm_ = rhs.bm_;
rs_ = rhs.rs_;
rs_.set_ptr(&bm_);
size_= rhs.size_;
values_ = rhs.values_;
return *this;
}
BinaryMultiset(BinaryMultiset&& rhs) noexcept
: bm_(std::move(rhs.bm_)), rs_(std::move(rhs.rs_)), size_(std::move(rhs.size_)), values_(std::move(rhs.values_)) {
rs_.set_ptr(&bm_);
}
BinaryMultiset& operator=(BinaryMultiset&& rhs) noexcept {
bm_ = std::move(rhs.bm_);
rs_ = std::move(rhs.rs_);
size_ = std::move(rhs.size_);
values_ = std::move(rhs.values_);
rs_.set_ptr(&bm_);
return *this;
}
~BinaryMultiset() = default;
void insert(value_type x) {
values_.push_back(x);
}
void build() {
if (!std::is_sorted(values_.begin(), values_.end()))
std::sort(values_.begin(), values_.end());
for (size_t i = 0; i < values_.size(); i++)
bm_.set((unsigned long long)values_[i] + i, 1);
size_ = values_.size();
values_.clear();
values_.shrink_to_fit();
bm_.build();
rs_.build(&bm_);
}
size_t size() const { return size_; }
bool empty() const { return size() == 0;}
size_t count(value_type x) const {
auto s = rs_.select0(x);
return s - (x == 0 ? 0 : (rs_.select0(x-1)+1));
}
bool contains(value_type x) const {
auto s = rs_.select0(x);
return s and bm_.get(s-1);
}
value_type lower_bound(value_type x) const {
return get(rank(x));
}
value_type upper_bound(value_type x) const {
return lower_bound(x+1);
}
size_t rank(value_type x) const {
return x <= (bm_.size()-size_) ? x == 0 ? 0 : rs_.select0(x-1) - (x-1) : size_;
}
value_type get(size_t i) const {
return rs_.select1(i) - i;
}
};
#line 10 "include/mtl/ordinal_range_search.hpp"
/**
* @brief Ordinal Range Search
* Maintain 2-dimential points (x,y) and their weights.
*/
template<class T, class W = void,
class BitmapType = RRR<>
>
class ORS {
public:
using index_type = T;
using weight_type = W;
// static constexpr index_type IndexLimit = std::numeric_limits<index_type>::max();
static constexpr bool kUseWeight = !std::is_same_v<weight_type, void>;
using point_data = std::conditional_t<
!kUseWeight,
std::tuple<index_type, index_type>,
std::tuple<index_type, index_type, weight_type>>;
using x_index_multiset = BinaryMultiset<BitmapType>;
using y_index_set = BinarySet<BitmapType>;
using wm_type = WaveletMatrix<int>;
private:
std::vector<point_data> ps;
x_index_multiset X;
y_index_set Y;
wm_type wm;
// using range_sum_type = FenwickTree<weight_type>;
// range_sum_type rsq;
index_type value_of_ith_x(size_t i) const {
return X.get(i);
}
size_t index_of_lower_bound_x(const T& x) const {
return X.rank(x);
}
std::pair<size_t, size_t> range_on_wm_lower_bound_x(const T& x) const {
auto l = index_of_lower_bound_x(x);
auto lb = value_of_ith_x(l);
auto c = X.count(lb);
return std::make_pair(l, l+c);
}
index_type value_of_ith_y(size_t i) const {
return Y.get(i);
}
size_t index_of_lower_bound_y(const T& y) const {
return Y.rank(y);
}
public:
ORS() = default;
void add(const point_data& point) {
ps.push_back(point);
}
void add(point_data&& point) {
ps.emplace_back(std::move(point));
}
template<class Resize, class Add>
void build(Resize resize = [](size_t){}, Add add = [](size_t, W) {}) {
std::sort(ps.begin(), ps.end(),
[](auto l, auto r) {
return std::get<0>(l) != std::get<0>(r) ? std::get<0>(l) < std::get<0>(r)
: std::get<1>(l) < std::get<1>(r);
});
for (size_t i = 0; i < ps.size(); i++) {
auto x = std::get<0>(ps[i]);
auto y = std::get<1>(ps[i]);
X.insert(x);
Y.insert(y);
}
X.build();
Y.build();
if constexpr (!kUseWeight) {
std::vector<int> S(ps.size());
for (size_t i = 0; i < ps.size(); i++) {
auto y = std::get<1>(ps[i]);
S[i] = index_of_lower_bound_y(y);
}
ps.clear();
ps.shrink_to_fit();
wm = wm_type(S.begin(), S.end());
} else {
std::vector<int> S(ps.size());
std::vector<std::pair<int, W>> SW(ps.size()), swz, swo;
for (size_t i = 0; i < ps.size(); i++) {
auto y = std::get<1>(ps[i]);
auto w = std::get<2>(ps[i]);
S[i] = index_of_lower_bound_y(y);
SW[i] = {S[i], w};
}
ps.clear();
ps.shrink_to_fit();
wm = wm_type(S.begin(), S.end());
// rsq = range_sum_type(wm.n * (wm.h+1));
resize(wm.n * (wm.h+1));
for (int k = wm.h-1; k >= 0; k--) {
for (size_t i = 0; i < wm.n; i++) {
if (((SW[i].first >> k) & 1u) == 0)
swz.push_back(SW[i]);
else
swo.push_back(SW[i]);
add(wm.n * (wm.h-1-k) + i, SW[i].second);
}
size_t j = wm.n-1;
while (!swo.empty()) {
SW[j--] = swo.back();
swo.pop_back();
}
while (!swz.empty()) {
SW[j--] = swz.back();
swz.pop_back();
}
}
for (size_t i = 0; i < wm.n; i++) {
add(wm.n * wm.h + i, SW[i].second);
}
}
}
private:
std::tuple<int, int, int, int>
compress_idx(const index_type& xl,
const index_type& xh,
const index_type& yl,
const index_type& yh) const {
return std::make_tuple(
index_of_lower_bound_x(xl),
index_of_lower_bound_x(xh),
index_of_lower_bound_y(yl),
index_of_lower_bound_y(yh)
);
}
template<typename F, typename G>
void _traverse(const index_type& xl,
const index_type& xh,
const index_type& yl,
const index_type& yh,
F internal_fn,
G leaf_fn) const {
size_t i,j;
int a,b;
std::tie(i,j,a,b) = compress_idx(xl, xh, yl, yh);
std::queue<std::tuple<size_t,size_t,T>> qs;
qs.emplace(i,j,T(0));
while (!qs.empty()) {
size_t l,r;
T c;
std::tie(l,r,c) = qs.front(); qs.pop();
size_t level = l/wm.n;
int shift = wm.h-1-level;
T clo = c, chi = c | ((1ull<<(shift+1))-1);
if (chi < a or b <= clo)
continue;
if (a <= clo and chi < b) {
internal_fn(l, r, c);
continue;
}
if (level == wm.h) {
leaf_fn(l, r, c);
continue;
}
size_t nl,nr;
std::tie(nl,nr) = wm._child_tie0(level, l, r);
qs.emplace(nl, nr, c);
std::tie(nl,nr) = wm._child_tie1(level, l, r);
qs.emplace(nl, nr, c | (1ull<<shift));
}
}
public:
size_t freq(const index_type& xl,
const index_type& xh,
const index_type& yl,
const index_type& yh) const {
size_t i,j;
int a,b;
std::tie(i,j,a,b) = compress_idx(xl, xh, yl, yh);
return wm.range_freq(i, j, a, b);
}
struct identity_add {
auto operator()(size_t) const {}
};
template<class Add = identity_add>
void weight_add(const index_type& x, const index_type& y, Add add = Add()) {
size_t l, r;
std::tie(l,r) = range_on_wm_lower_bound_x(x);
auto c = index_of_lower_bound_y(y);
for (int k = wm.h-1; k >= 0; k--) {
std::tie(l,r) = wm.child_tie(wm.h-1-k, l, r, (c >> k) & 1u);
}
assert(l != r);
size_t j = l;
for (size_t k = 0; k < wm.h; k++) {
add(j);
j = wm.parent(wm.h-1-k, j, (c >> k) & 1u);
}
add(j);
}
struct plus_op {
template<class U, class V>
auto operator()(const U& a, const V& b) const {
return a+b;
}
};
struct plus_identity {
auto operator()() const {
return 0;
}
};
template<class RangeSum, class Sum = plus_op, class Id = plus_identity>
auto sum(
const index_type& xl,
const index_type& xh,
const index_type& yl,
const index_type& yh,
RangeSum range_sum = [](size_t l, size_t r) { return r-l; },
Sum sum_op = Sum(),
Id id = Id()
) const -> decltype(range_sum(0,0)) {
using sum_type = decltype(range_sum(0,0));
sum_type ret = id();
_traverse(xl, xh, yl, yh,
[&](size_t l, size_t r, auto _) { ret = sum_op(ret, range_sum(l, r)); },
[](auto,auto,auto){});
return ret;
}
std::vector<std::tuple<index_type, index_type>>
enumerate(const index_type& xl,
const index_type& xh,
const index_type& yl,
const index_type& yh) const {
std::vector<std::tuple<index_type, index_type>> ret;
_traverse(xl, xh, yl, yh,
[](auto,auto,auto){},
[&](size_t l, size_t r, T c) {
for (size_t idx = l; idx < r; idx++) {
auto jdx = idx;
for (size_t k = 0; k < wm.h; k++)
jdx = wm.parent(wm.h-1-k, jdx, (c >> k) & 1u);
ret.emplace_back(value_of_ith_x(jdx), value_of_ith_y(c));
}
});
return ret;
}
};
#line 7 "test/yosupo/rectangle_sum.test.cpp"
#include <bits/stdc++.h>
using namespace std;
int main() {
int n,q; cin>>n>>q;
vector<array<int,3>> N(n);
int max_idx = 0;
auto chmax = [&](int& a, int b) { a = max(a, b); };
for (int i = 0; i < n; i++) {
int x,y,w; cin>>x>>y>>w;
N[i] = {x,y,w};
chmax(max_idx, max(x, y));
}
vector<array<int,4>> Q(q);
for (int i = 0; i < q; i++) {
int l,d,r,u; cin>>l>>d>>r>>u;
Q[i] = {l,d,r,u};
chmax(max_idx, max(r, u));
}
ORS<int, long long> ors;
for (auto [x,y,w] : N) ors.add({x, y, w});
FenwickTree<long long> rsq;
ors.build(
[&](size_t size) {
rsq = FenwickTree<long long>(size);
},
[&](size_t i, int w) {
rsq.add(i, w);
});
for (auto [l,d,r,u] : Q) {
cout << ors.sum(l, r, d, u, [&](size_t l, size_t r) { return rsq.sum(l, r); }) << endl;
}
}
Test cases
| Env | Name | Status | Elapsed | Memory |
|---|---|---|---|---|
| g++ | example_00 |
AC |
6 ms | 3 MB |
| g++ | max_random_00 |
AC |
2051 ms | 256 MB |
| g++ | max_random_01 |
AC |
2209 ms | 256 MB |
| g++ | max_random_02 |
AC |
2191 ms | 253 MB |
| g++ | n_131072_00 |
AC |
931 ms | 33 MB |
| g++ | random_00 |
AC |
1567 ms | 235 MB |
| g++ | random_01 |
AC |
1655 ms | 243 MB |
| g++ | random_02 |
AC |
1142 ms | 210 MB |
| g++ | small_00 |
AC |
8 ms | 4 MB |
| g++ | small_01 |
AC |
6 ms | 4 MB |
| g++ | small_02 |
AC |
6 ms | 3 MB |
| g++ | xy_concentrate_00 |
AC |
1804 ms | 245 MB |
| g++ | xy_concentrate_01 |
AC |
1825 ms | 244 MB |
| g++ | xy_concentrate_02 |
AC |
2289 ms | 244 MB |
| clang++ | example_00 |
AC |
6 ms | 3 MB |
| clang++ | max_random_00 |
AC |
2211 ms | 254 MB |
| clang++ | max_random_01 |
AC |
2168 ms | 256 MB |
| clang++ | max_random_02 |
AC |
2500 ms | 253 MB |
| clang++ | n_131072_00 |
AC |
1167 ms | 32 MB |
| clang++ | random_00 |
AC |
1611 ms | 235 MB |
| clang++ | random_01 |
AC |
1686 ms | 244 MB |
| clang++ | random_02 |
AC |
1306 ms | 211 MB |
| clang++ | small_00 |
AC |
8 ms | 4 MB |
| clang++ | small_01 |
AC |
6 ms | 3 MB |
| clang++ | small_02 |
AC |
7 ms | 3 MB |
| clang++ | xy_concentrate_00 |
AC |
1939 ms | 244 MB |
| clang++ | xy_concentrate_01 |
AC |
1928 ms | 247 MB |
| clang++ | xy_concentrate_02 |
AC |
1978 ms | 244 MB |