Wavelet matrix (tools/wavelet_matrix.hpp)

View this file on GitHub
View document part on GitHub
Last update: 2025-02-10 20:05:30+09:00
Include: #include "tools/wavelet_matrix.hpp"

It is a data structure which can process various queries at high speed for $n$ weighted points $(x_0, y_0, w_0), (x_1, y_1, w_1), \ldots, (x_{n - 1}, y_{n - 1}, w_{n - 1})$ on a two-dimensional plane.

License

Author

anqooqie

Constructor

wavelet_matrix<T> wm();

It creates an empty two-dimensional plane.

Constraints

None

Time Complexity

$O(1)$

size

std::size_t wm.size();

It returns the number of points on the plane.

Constraints

None

Time Complexity

$O(1)$

add_point

std::size_t wm.add_point(T x, T y);

It adds a point $(x, y)$ to the plane. It returns an integer $i$ such that this is the $i$-th ($0$ indexed) point that is added.

Constraints

wm.build() has not been called ever.

Time Complexity

$O(1)$ amortized

get_point

std::pair<T, T> wm.get_point(std::size_t i)

It returns the $i$-th ($0$ indexed) point.

Constraints

$i < n$ where $n$ is wm.size()

Time Complexity

$O(1)$

points

const std::vector<std::pair<T, T>>& wm.points();

It returns all the points on the plane. The points are ordered in the same order as added by add_point.

Constraints

None

Time Complexity

$O(1)$

build

std::vector<std::vector<std::size_t>> wm.build();

It internally creates the data structure called as wavelet matrix.

It also returns the matrix $A_{h, j}$ with $O(\log n)$ rows and $n$ columns. Each row of $A_h$ is the permutation of $(0, 1, \ldots, n - 1)$. You can initialize an auxiliary data structure by $A_{h, j}$, for querying weighted points.

Example

// Given n weighted points, answer the sum of w_i such that l <= x_i < r and d <= y_i < u.

tools::wavelet_matrix<int> wm;
for (int i = 0; i < n; ++i) {
  wm.add_point(x[i], y[i]);
}

std::vector<tools::online_cumsum<int>> cumsums;
for (const auto& A_h : wm.build()) {
  cumsums.emplace_back(n);
  for (int j = 0; j < n; ++j) {
    cumsums.back()[j] = w[A_h[j]];
  }
}

ll answer = 0;
for (const auto& [h, jl, jr] : wm.range_prod(l, r, u)) {
  answer += cumsums[h].sum(jl, jr);
}
for (const auto& [h, jl, jr] : wm.range_prod(l, r, d)) {
  answer -= cumsums[h].sum(jl, jr);
}
std::cout << answer << '\n';

Constraints

wm.build() has not been called ever.

Time Complexity

$O(n \log n)$

apply

std::vector<std::pair<std::size_t, std::size_t>> wm.apply(std::size_t i);

It has no side effects, but just returns pairs $(0, j_0), (1, j_1), \ldots, (H - 1, j_{H - 1})$ where $H$ is the number of rows of the returned matrix from wm.build(). You can update the auxiliary data structure by the returned pairs, for querying weighted points.

Example

// Answer Q queries.
// (type 1 query) 1 x y w: Add a point (x, y) at weight w.
// (type 2 query) 2 l r d u: Print the sum of w_i such that l <= x_i < r and d <= y_i < u.

std::queue<int> query_types;
tools::wavelet_matrix<int> wm;
std::vector<int> weights;
std::queue<std::tuple<int, int, int, int>> query2;

int Q;
std::cin >> Q;
for (int q = 0; q < Q; ++q) {
  int t;
  std::cin >> t;
  query_types.push(t);
  if (t == 1) {
    int x, y, w;
    wm.add_point(x, y);
    weights.push_back(w);
  } else {
    int l, r, d, u;
    std::cin >> l >> r >> d >> u;
    query2.emplace(l, r, d, u);
  }
}

auto fws = std::vector(wm.build().size(), atcoder::fenwick_tree<ll>(wm.size()));

int i = 0;
while (!query_types.empty()) {
  const auto t = query_types.front();
  query_types.pop();
  if (t == 1) {
    for (const auto& [h, j] : wm.apply(i)) {
      fws[h].add(j, weights[i]);
    }
    ++i;
  } else {
    const auto [l, r, d, u] = query2.front();
    query2.pop();

    ll answer = 0;
    for (const auto& [h, jl, jr] : wm.range_prod(l, r, u)) {
      answer += fws[h].sum(jl, jr);
    }
    for (const auto& [h, jl, jr] : wm.range_prod(l, r, d)) {
      answer -= fws[h].sum(jl, jr);
    }
    std::cout << answer << '\n';
  }
}

Constraints

wm.build() has been called ever.
$i < n$

Time Complexity

$O(\log n)$

kth_smallest

std::size_t wm.kth_smallest(T l, T r, std::size_t k);

It returns $i$ such that $l \leq x_i < r$ and $|\{j \mid l \leq x_j < r \land (y_j, j) < (y_i, i)) \}| = k$.

Constraints

wm.build() has been called ever.
$l \leq r$
$k < |\{i \mid l \leq x_i < r\}|$

Time Complexity

$O(\log n)$

kth_largest

std::size_t wm.kth_largest(T l, T r, std::size_t k);

It returns $i$ such that $l \leq x_i < r$ and $|\{j \mid l \leq x_j < r \land (y_i, i) < (y_j, j)) \}| = k$.

Constraints

wm.build() has been called ever.
$l \leq r$
$k < |\{i \mid l \leq x_i < r\}|$

Time Complexity

$O(\log n)$

range_prod

std::vector<std::tuple<std::size_t, std::size_t, std::size_t>> wm.range_prod(T l, T r, T u);

It has no side effects, but just returns tuples $(0, l_0, r_0), (1, l_1, r_1), \ldots, (H - 1, l_{H - 1}, r_{H - 1})$ where $H$ is the number of rows of the returned matrix from wm.build(). You can answer the query about weighted points using the auxiliary data structure and the returned tuples.

Example

See the examples of build and apply.

Constraints

wm.build() has been called ever.
$l \leq r$

Time Complexity

$O(\log n)$

range_freq

(1) std::size_t wm.range_freq(T l, T r);
(2) std::size_t wm.range_freq(T l, T r, T u);
(3) std::size_t wm.range_freq(T l, T r, T d, T u);

(1)
- It returns $|\{i \mid l \leq x_i < r \}|$.
(2)
- It returns $|\{i \mid l \leq x_i < r \land y_i < u \}|$.
(3)
- It returns $|\{i \mid l \leq x_i < r \land d \leq y_i < u \}|$.

Constraints

(1), (2)
- wm.build() has been called ever.
- $l \leq r$
(3)
- wm.build() has been called ever.
- $l \leq r$
- $d \leq u$

Time Complexity

$O(\log n)$

prev_value

std::optional<T> wm.prev_value(T l, T r, T u);

It returns the maximum $y_i$ such that $l \leq x_i < r$ and $y_i < u$. If such the $y_i$ does not exist, it returns std::nullopt.

Constraints

wm.build() has been called ever.
$l \leq r$

Time Complexity

$O(\log n)$

next_value

std::optional<T> wm.next_value(T l, T r, T d);

It returns the minimum $y_i$ such that $l \leq x_i < r$ and $d \leq y_i$. If such the $y_i$ does not exist, it returns std::nullopt.

Constraints

wm.build() has been called ever.
$l \leq r$

Time Complexity

$O(\log n)$

prev_points

std::pair<typename std::vector<std::size_t>::const_iterator, typename std::vector<std::size_t>::const_iterator> prev_points(T l, T r, T u)

Let us denote the returned iterator pair by begin and end. std::vector<std::size_t>(begin, end) would be $\{i \mid l \leq x_i < r \land y_i = y\}$ sorted in ascending order, where $y$ is wm.prev_value(l, r, u).

Constraints

wm.build() has been called ever.
$l \leq r$

Time Complexity

$O(\log n)$

next_points

std::pair<typename std::vector<std::size_t>::const_iterator, typename std::vector<std::size_t>::const_iterator> next_points(T l, T r, T d)

Constraints

wm.build() has been called ever.
$l \leq r$

Time Complexity

$O(\log n)$

Depends on

Required by

Range count distinct (tools/range_count_distinct.hpp)

Verified with

Code

#ifndef TOOLS_WAVELET_MATRIX_HPP
#define TOOLS_WAVELET_MATRIX_HPP

#include <vector>
#include <utility>
#include <cstddef>
#include <cassert>
#include <algorithm>
#include <iterator>
#include <array>
#include <tuple>
#include <optional>
#include "tools/bit_vector.hpp"
#include "tools/lower_bound.hpp"
#include "tools/less_by_first.hpp"
#include "tools/floor_log2.hpp"

namespace tools {
  template <typename T>
  class wavelet_matrix {
  private:
    ::std::vector<::std::pair<T, T>> m_ps;
    ::std::vector<::std::pair<T, ::std::size_t>> m_xs;
    ::std::vector<T> m_ys;
    ::std::vector<::std::size_t> m_is;
    ::std::vector<::tools::bit_vector> m_bvs;

    ::std::size_t iid(const ::std::size_t i) const {
      return ::tools::lower_bound(this->m_xs.begin(), this->m_xs.end(), ::std::make_pair(this->m_ps[i].first, i));
    }
    ::std::size_t xid(const T& x) const {
      return ::tools::lower_bound(this->m_xs.begin(), this->m_xs.end(), ::std::make_pair(x, ::std::size_t{}), ::tools::less_by_first());
    }
    ::std::size_t yid(const T& y) const {
      return ::tools::lower_bound(this->m_ys.begin(), this->m_ys.end(), y);
    }
    ::std::size_t lg() const {
      return this->m_bvs.size();
    }
    bool built() const {
      return this->lg() > 0;
    }

  public:
    wavelet_matrix() = default;
    wavelet_matrix(const ::tools::wavelet_matrix<T>&) = default;
    wavelet_matrix(::tools::wavelet_matrix<T>&&) = default;
    ~wavelet_matrix() = default;
    ::tools::wavelet_matrix<T>& operator=(const ::tools::wavelet_matrix<T>&) = default;
    ::tools::wavelet_matrix<T>& operator=(::tools::wavelet_matrix<T>&&) = default;

    ::std::size_t size() const {
      return this->m_ps.size();
    }

    ::std::size_t add_point(const T& x, const T& y) {
      assert(!this->built());
      this->m_ps.emplace_back(x, y);
      return this->size() - 1;
    }

    ::std::pair<T, T> get_point(const ::std::size_t i) const {
      assert(i < this->size());
      return this->m_ps[i];
    }

    const ::std::vector<::std::pair<T, T>>& points() const {
      return this->m_ps;
    }

    ::std::vector<::std::vector<::std::size_t>> build() {
      assert(!this->built());

      this->m_xs.reserve(this->size());
      for (::std::size_t i = 0; i < this->size(); ++i) {
        this->m_xs.emplace_back(this->m_ps[i].first, i);
      }
      ::std::sort(this->m_xs.begin(), this->m_xs.end());

      this->m_ys.reserve(this->size());
      ::std::transform(this->m_ps.begin(), this->m_ps.end(), ::std::back_inserter(this->m_ys), [](const auto& p) { return p.second; });
      ::std::sort(this->m_ys.begin(), this->m_ys.end());
      this->m_ys.erase(::std::unique(this->m_ys.begin(), this->m_ys.end()), this->m_ys.end());

      const auto n = ::std::max<::std::size_t>(this->size(), 1);
      this->m_bvs.assign(::tools::floor_log2(::std::max<::std::size_t>(this->m_ys.size(), 1)) + 1, ::tools::bit_vector(n));

      ::std::vector<::std::size_t> cur;
      cur.reserve(n);
      ::std::transform(this->m_xs.begin(), this->m_xs.end(), ::std::back_inserter(cur), [&](const auto& p) { return this->yid(this->m_ps[p.second].second); });
      cur.resize(n);
      ::std::vector<::std::size_t> nxt(n);

      auto res = ::std::vector(this->lg() + 1, ::std::vector<::std::size_t>(n));
      ::std::transform(this->m_xs.begin(), this->m_xs.end(), res.back().begin(), [&](const auto& p) { return p.second; });

      for (::std::size_t h = this->lg(); h --> 0;) {
        for (::std::size_t i = 0; i < n; ++i) {
          if ((cur[i] >> h) & 1) {
            this->m_bvs[h].set(i);
          }
        }
        this->m_bvs[h].build();
        ::std::array<::std::size_t, 2> offsets = {0, this->m_bvs[h].zeros()};
        for (::std::size_t i = 0; i < n; ++i) {
          auto& offset = offsets[this->m_bvs[h].get(i)];
          nxt[offset] = cur[i];
          res[h][offset] = res[h + 1][i];
          ++offset;
        }
        ::std::swap(cur, nxt);
      }

      this->m_is = res.front();
      res.pop_back();
      return res;
    }

    ::std::vector<::std::pair<::std::size_t, ::std::size_t>> apply(::std::size_t i) const {
      assert(this->built());
      assert(i < this->size());

      ::std::vector<::std::pair<::std::size_t, ::std::size_t>> res(this->lg());
      i = this->iid(i);
      for (::std::size_t h = this->lg(); h --> 0;) {
        i = this->m_bvs[h].get(i) ? this->m_bvs[h].zeros() + this->m_bvs[h].rank1(i) : this->m_bvs[h].rank0(i);
        res[h] = ::std::make_pair(h, i);
      }
      return res;
    }

    ::std::size_t kth_smallest(const T& l, const T& r, ::std::size_t k) const {
      assert(this->built());
      assert(l <= r);

      auto lid = this->xid(l);
      auto rid = this->xid(r);

      assert(k < rid - lid);

      for (::std::size_t h = this->lg(); h --> 0;) {
        const auto l0 = this->m_bvs[h].rank0(lid);
        const auto r0 = this->m_bvs[h].rank0(rid);
        if (k < r0 - l0) {
          lid = l0;
          rid = r0;
        } else {
          k -= r0 - l0;
          lid += this->m_bvs[h].zeros() - l0;
          rid += this->m_bvs[h].zeros() - r0;
        }
      }

      return this->m_is[lid + k];
    }

    ::std::size_t kth_largest(const T& l, const T& r, const ::std::size_t k) const {
      assert(this->built());
      assert(l <= r);

      const auto lid = this->xid(l);
      const auto rid = this->xid(r);

      assert(k < rid - lid);

      return this->kth_smallest(l, r, rid - lid - k - 1);
    }

    ::std::vector<::std::tuple<::std::size_t, ::std::size_t, ::std::size_t>> range_prod(const T& l, const T& r, const T& u) const {
      assert(this->built());
      assert(l <= r);

      auto lid = this->xid(l);
      auto rid = this->xid(r);
      const auto uid = this->yid(u);

      ::std::vector<::std::tuple<::std::size_t, ::std::size_t, ::std::size_t>> res(this->lg());
      for (::std::size_t h = this->lg(); h --> 0;) {
        const auto l0 = this->m_bvs[h].rank0(lid);
        const auto r0 = this->m_bvs[h].rank0(rid);
        if ((uid >> h) & 1) {
          res[h] = ::std::make_tuple(h, l0, r0);
          lid += this->m_bvs[h].zeros() - l0;
          rid += this->m_bvs[h].zeros() - r0;
        } else {
          res[h] = ::std::make_tuple(h, 0, 0);
          lid = l0;
          rid = r0;
        }
      }
      return res;
    }

    ::std::size_t range_freq(const T& l, const T& r) const {
      assert(this->built());
      assert(l <= r);

      return this->xid(r) - this->xid(l);
    }

    ::std::size_t range_freq(const T& l, const T& r, const T& u) const {
      assert(this->built());
      assert(l <= r);

      auto lid = this->xid(l);
      auto rid = this->xid(r);
      const auto uid = this->yid(u);

      ::std::size_t res = 0;
      for (::std::size_t h = this->lg(); h --> 0;) {
        const auto l0 = this->m_bvs[h].rank0(lid);
        const auto r0 = this->m_bvs[h].rank0(rid);
        if ((uid >> h) & 1) {
          res += r0 - l0;
          lid += this->m_bvs[h].zeros() - l0;
          rid += this->m_bvs[h].zeros() - r0;
        } else {
          lid = l0;
          rid = r0;
        }
      }
      return res;
    }

    ::std::size_t range_freq(const T& l, const T& r, const T& d, const T& u) const {
      assert(this->built());
      assert(l <= r);
      assert(d <= u);

      return this->range_freq(l, r, u) - this->range_freq(l, r, d);
    }

    ::std::optional<T> prev_value(const T& l, const T& r, const T& u) const {
      assert(this->built());
      assert(l <= r);

      const auto k = this->range_freq(l, r, u);
      return k == 0 ? ::std::nullopt : ::std::make_optional(this->m_ps[this->kth_smallest(l, r, k - 1)].second);
    }

    ::std::optional<T> next_value(const T& l, const T& r, const T& d) const {
      assert(this->built());
      assert(l <= r);

      const auto k = this->range_freq(l, r, d);
      return k == this->range_freq(l, r) ? ::std::nullopt : ::std::make_optional(this->m_ps[this->kth_smallest(l, r, k)].second);
    }

    ::std::pair<typename ::std::vector<::std::size_t>::const_iterator, typename ::std::vector<::std::size_t>::const_iterator> prev_points(const T& l, const T& r, const T& u) const {
      assert(this->built());
      assert(l <= r);

      auto lid = this->xid(l);
      auto rid = this->xid(r);
      auto k = this->range_freq(l, r, u);
      if (k == 0) return ::std::make_pair(this->m_is.cend(), this->m_is.cend());
      --k;

      for (::std::size_t h = this->lg(); h --> 0;) {
        const auto l0 = this->m_bvs[h].rank0(lid);
        const auto r0 = this->m_bvs[h].rank0(rid);
        if (k < r0 - l0) {
          lid = l0;
          rid = r0;
        } else {
          k -= r0 - l0;
          lid += this->m_bvs[h].zeros() - l0;
          rid += this->m_bvs[h].zeros() - r0;
        }
      }

      return ::std::make_pair(this->m_is.cbegin() + lid, this->m_is.cbegin() + rid);
    }

    ::std::pair<typename ::std::vector<::std::size_t>::const_iterator, typename ::std::vector<::std::size_t>::const_iterator> next_points(const T& l, const T& r, const T& d) const {
      assert(this->built());
      assert(l <= r);

      auto lid = this->xid(l);
      auto rid = this->xid(r);
      auto k = this->range_freq(l, r, d);
      if (k == rid - lid) return ::std::make_pair(this->m_is.cend(), this->m_is.cend());

      for (::std::size_t h = this->lg(); h --> 0;) {
        const auto l0 = this->m_bvs[h].rank0(lid);
        const auto r0 = this->m_bvs[h].rank0(rid);
        if (k < r0 - l0) {
          lid = l0;
          rid = r0;
        } else {
          k -= r0 - l0;
          lid += this->m_bvs[h].zeros() - l0;
          rid += this->m_bvs[h].zeros() - r0;
        }
      }

      return ::std::make_pair(this->m_is.cbegin() + lid, this->m_is.cbegin() + rid);
    }
  };
}

#endif

#line 1 "tools/wavelet_matrix.hpp"



#include <vector>
#include <utility>
#include <cstddef>
#include <cassert>
#include <algorithm>
#include <iterator>
#include <array>
#include <tuple>
#include <optional>
#line 1 "tools/bit_vector.hpp"



#include <cstdint>
#line 6 "tools/bit_vector.hpp"
#include <immintrin.h>

// Source: https://nyaannyaan.github.io/library/data-structure-2d/wavelet-matrix.hpp.html
// License: CC0 1.0 Universal
// Author: Nyaan

namespace tools {
  class bit_vector {
  private:
    using u32 = ::std::uint32_t;
    using i64 = ::std::int64_t;
    using u64 = ::std::uint64_t;

    static constexpr u32 w = 64;
    ::std::vector<u64> m_block;
    ::std::vector<u32> m_count;
    u32 m_size, m_zeros;

  public:
    bit_vector() {}
    explicit bit_vector(int _n) { init(_n); }
    __attribute__((optimize("O3", "unroll-loops"))) void init(int _n) {
      m_size = m_zeros = _n;
      m_block.resize(m_size / w + 1, 0);
      m_count.resize(m_block.size(), 0);
    }

    u32 size() const { return m_size; }
    inline u32 get(u32 i) const { return u32(m_block[i / w] >> (i % w)) & 1u; }
    inline void set(u32 i) { m_block[i / w] |= 1LL << (i % w); }

    __attribute__((target("popcnt"))) void build() {
      for (u32 i = 1; i < m_block.size(); ++i)
        m_count[i] = m_count[i - 1] + _mm_popcnt_u64(m_block[i - 1]);
      m_zeros = rank0(m_size);
    }

    u32 zeros() const { return m_zeros; }
    inline u32 rank0(u32 i) const { return i - rank1(i); }
    __attribute__((target("bmi2,popcnt"))) inline u32 rank1(u32 i) const {
      return m_count[i / w] + _mm_popcnt_u64(_bzhi_u64(m_block[i / w], i % w));
    }
  };
}


#line 1 "tools/lower_bound.hpp"



#line 6 "tools/lower_bound.hpp"

namespace tools {

  template <class ForwardIterator, class T>
  auto lower_bound(ForwardIterator first, ForwardIterator last, const T& value) {
    return ::std::distance(first, ::std::lower_bound(first, last, value));
  }

  template <class ForwardIterator, class T, class Compare>
  auto lower_bound(ForwardIterator first, ForwardIterator last, const T& value, Compare comp) {
    return ::std::distance(first, ::std::lower_bound(first, last, value, comp));
  }
}


#line 1 "tools/less_by_first.hpp"



#line 5 "tools/less_by_first.hpp"

namespace tools {

  class less_by_first {
  public:
    template <class T1, class T2>
    bool operator()(const ::std::pair<T1, T2>& x, const ::std::pair<T1, T2>& y) const {
      return x.first < y.first;
    }
  };
}


#line 1 "tools/floor_log2.hpp"



#line 1 "tools/bit_width.hpp"



#include <bit>
#line 6 "tools/bit_width.hpp"
#include <type_traits>
#line 1 "tools/is_integral.hpp"



#line 5 "tools/is_integral.hpp"

namespace tools {
  template <typename T>
  struct is_integral : ::std::is_integral<T> {};

  template <typename T>
  inline constexpr bool is_integral_v = ::tools::is_integral<T>::value;
}


#line 1 "tools/is_signed.hpp"



#line 5 "tools/is_signed.hpp"

namespace tools {
  template <typename T>
  struct is_signed : ::std::is_signed<T> {};

  template <typename T>
  inline constexpr bool is_signed_v = ::tools::is_signed<T>::value;
}


#line 1 "tools/make_unsigned.hpp"



#line 5 "tools/make_unsigned.hpp"

namespace tools {
  template <typename T>
  struct make_unsigned : ::std::make_unsigned<T> {};

  template <typename T>
  using make_unsigned_t = typename ::tools::make_unsigned<T>::type;
}


#line 10 "tools/bit_width.hpp"

namespace tools {
  template <typename T>
  constexpr int bit_width(T) noexcept;

  template <typename T>
  constexpr int bit_width(const T x) noexcept {
    static_assert(::tools::is_integral_v<T> && !::std::is_same_v<::std::remove_cv_t<T>, bool>);
    if constexpr (::tools::is_signed_v<T>) {
      assert(x >= 0);
      return ::tools::bit_width<::tools::make_unsigned_t<T>>(x);
    } else {
      return ::std::bit_width(x);
    }
  }
}


#line 6 "tools/floor_log2.hpp"

namespace tools {
  template <typename T>
  constexpr T floor_log2(T x) noexcept {
    assert(x > 0);
    return ::tools::bit_width(x) - 1;
  }
}


#line 17 "tools/wavelet_matrix.hpp"

namespace tools {
  template <typename T>
  class wavelet_matrix {
  private:
    ::std::vector<::std::pair<T, T>> m_ps;
    ::std::vector<::std::pair<T, ::std::size_t>> m_xs;
    ::std::vector<T> m_ys;
    ::std::vector<::std::size_t> m_is;
    ::std::vector<::tools::bit_vector> m_bvs;

    ::std::size_t iid(const ::std::size_t i) const {
      return ::tools::lower_bound(this->m_xs.begin(), this->m_xs.end(), ::std::make_pair(this->m_ps[i].first, i));
    }
    ::std::size_t xid(const T& x) const {
      return ::tools::lower_bound(this->m_xs.begin(), this->m_xs.end(), ::std::make_pair(x, ::std::size_t{}), ::tools::less_by_first());
    }
    ::std::size_t yid(const T& y) const {
      return ::tools::lower_bound(this->m_ys.begin(), this->m_ys.end(), y);
    }
    ::std::size_t lg() const {
      return this->m_bvs.size();
    }
    bool built() const {
      return this->lg() > 0;
    }

  public:
    wavelet_matrix() = default;
    wavelet_matrix(const ::tools::wavelet_matrix<T>&) = default;
    wavelet_matrix(::tools::wavelet_matrix<T>&&) = default;
    ~wavelet_matrix() = default;
    ::tools::wavelet_matrix<T>& operator=(const ::tools::wavelet_matrix<T>&) = default;
    ::tools::wavelet_matrix<T>& operator=(::tools::wavelet_matrix<T>&&) = default;

    ::std::size_t size() const {
      return this->m_ps.size();
    }

    ::std::size_t add_point(const T& x, const T& y) {
      assert(!this->built());
      this->m_ps.emplace_back(x, y);
      return this->size() - 1;
    }

    ::std::pair<T, T> get_point(const ::std::size_t i) const {
      assert(i < this->size());
      return this->m_ps[i];
    }

    const ::std::vector<::std::pair<T, T>>& points() const {
      return this->m_ps;
    }

    ::std::vector<::std::vector<::std::size_t>> build() {
      assert(!this->built());

      this->m_xs.reserve(this->size());
      for (::std::size_t i = 0; i < this->size(); ++i) {
        this->m_xs.emplace_back(this->m_ps[i].first, i);
      }
      ::std::sort(this->m_xs.begin(), this->m_xs.end());

      this->m_ys.reserve(this->size());
      ::std::transform(this->m_ps.begin(), this->m_ps.end(), ::std::back_inserter(this->m_ys), [](const auto& p) { return p.second; });
      ::std::sort(this->m_ys.begin(), this->m_ys.end());
      this->m_ys.erase(::std::unique(this->m_ys.begin(), this->m_ys.end()), this->m_ys.end());

      const auto n = ::std::max<::std::size_t>(this->size(), 1);
      this->m_bvs.assign(::tools::floor_log2(::std::max<::std::size_t>(this->m_ys.size(), 1)) + 1, ::tools::bit_vector(n));

      ::std::vector<::std::size_t> cur;
      cur.reserve(n);
      ::std::transform(this->m_xs.begin(), this->m_xs.end(), ::std::back_inserter(cur), [&](const auto& p) { return this->yid(this->m_ps[p.second].second); });
      cur.resize(n);
      ::std::vector<::std::size_t> nxt(n);

      auto res = ::std::vector(this->lg() + 1, ::std::vector<::std::size_t>(n));
      ::std::transform(this->m_xs.begin(), this->m_xs.end(), res.back().begin(), [&](const auto& p) { return p.second; });

      for (::std::size_t h = this->lg(); h --> 0;) {
        for (::std::size_t i = 0; i < n; ++i) {
          if ((cur[i] >> h) & 1) {
            this->m_bvs[h].set(i);
          }
        }
        this->m_bvs[h].build();
        ::std::array<::std::size_t, 2> offsets = {0, this->m_bvs[h].zeros()};
        for (::std::size_t i = 0; i < n; ++i) {
          auto& offset = offsets[this->m_bvs[h].get(i)];
          nxt[offset] = cur[i];
          res[h][offset] = res[h + 1][i];
          ++offset;
        }
        ::std::swap(cur, nxt);
      }

      this->m_is = res.front();
      res.pop_back();
      return res;
    }

    ::std::vector<::std::pair<::std::size_t, ::std::size_t>> apply(::std::size_t i) const {
      assert(this->built());
      assert(i < this->size());

      ::std::vector<::std::pair<::std::size_t, ::std::size_t>> res(this->lg());
      i = this->iid(i);
      for (::std::size_t h = this->lg(); h --> 0;) {
        i = this->m_bvs[h].get(i) ? this->m_bvs[h].zeros() + this->m_bvs[h].rank1(i) : this->m_bvs[h].rank0(i);
        res[h] = ::std::make_pair(h, i);
      }
      return res;
    }

    ::std::size_t kth_smallest(const T& l, const T& r, ::std::size_t k) const {
      assert(this->built());
      assert(l <= r);

      auto lid = this->xid(l);
      auto rid = this->xid(r);

      assert(k < rid - lid);

      for (::std::size_t h = this->lg(); h --> 0;) {
        const auto l0 = this->m_bvs[h].rank0(lid);
        const auto r0 = this->m_bvs[h].rank0(rid);
        if (k < r0 - l0) {
          lid = l0;
          rid = r0;
        } else {
          k -= r0 - l0;
          lid += this->m_bvs[h].zeros() - l0;
          rid += this->m_bvs[h].zeros() - r0;
        }
      }

      return this->m_is[lid + k];
    }

    ::std::size_t kth_largest(const T& l, const T& r, const ::std::size_t k) const {
      assert(this->built());
      assert(l <= r);

      const auto lid = this->xid(l);
      const auto rid = this->xid(r);

      assert(k < rid - lid);

      return this->kth_smallest(l, r, rid - lid - k - 1);
    }

    ::std::vector<::std::tuple<::std::size_t, ::std::size_t, ::std::size_t>> range_prod(const T& l, const T& r, const T& u) const {
      assert(this->built());
      assert(l <= r);

      auto lid = this->xid(l);
      auto rid = this->xid(r);
      const auto uid = this->yid(u);

      ::std::vector<::std::tuple<::std::size_t, ::std::size_t, ::std::size_t>> res(this->lg());
      for (::std::size_t h = this->lg(); h --> 0;) {
        const auto l0 = this->m_bvs[h].rank0(lid);
        const auto r0 = this->m_bvs[h].rank0(rid);
        if ((uid >> h) & 1) {
          res[h] = ::std::make_tuple(h, l0, r0);
          lid += this->m_bvs[h].zeros() - l0;
          rid += this->m_bvs[h].zeros() - r0;
        } else {
          res[h] = ::std::make_tuple(h, 0, 0);
          lid = l0;
          rid = r0;
        }
      }
      return res;
    }

    ::std::size_t range_freq(const T& l, const T& r) const {
      assert(this->built());
      assert(l <= r);

      return this->xid(r) - this->xid(l);
    }

    ::std::size_t range_freq(const T& l, const T& r, const T& u) const {
      assert(this->built());
      assert(l <= r);

      auto lid = this->xid(l);
      auto rid = this->xid(r);
      const auto uid = this->yid(u);

      ::std::size_t res = 0;
      for (::std::size_t h = this->lg(); h --> 0;) {
        const auto l0 = this->m_bvs[h].rank0(lid);
        const auto r0 = this->m_bvs[h].rank0(rid);
        if ((uid >> h) & 1) {
          res += r0 - l0;
          lid += this->m_bvs[h].zeros() - l0;
          rid += this->m_bvs[h].zeros() - r0;
        } else {
          lid = l0;
          rid = r0;
        }
      }
      return res;
    }

    ::std::size_t range_freq(const T& l, const T& r, const T& d, const T& u) const {
      assert(this->built());
      assert(l <= r);
      assert(d <= u);

      return this->range_freq(l, r, u) - this->range_freq(l, r, d);
    }

    ::std::optional<T> prev_value(const T& l, const T& r, const T& u) const {
      assert(this->built());
      assert(l <= r);

      const auto k = this->range_freq(l, r, u);
      return k == 0 ? ::std::nullopt : ::std::make_optional(this->m_ps[this->kth_smallest(l, r, k - 1)].second);
    }

    ::std::optional<T> next_value(const T& l, const T& r, const T& d) const {
      assert(this->built());
      assert(l <= r);

      const auto k = this->range_freq(l, r, d);
      return k == this->range_freq(l, r) ? ::std::nullopt : ::std::make_optional(this->m_ps[this->kth_smallest(l, r, k)].second);
    }

    ::std::pair<typename ::std::vector<::std::size_t>::const_iterator, typename ::std::vector<::std::size_t>::const_iterator> prev_points(const T& l, const T& r, const T& u) const {
      assert(this->built());
      assert(l <= r);

      auto lid = this->xid(l);
      auto rid = this->xid(r);
      auto k = this->range_freq(l, r, u);
      if (k == 0) return ::std::make_pair(this->m_is.cend(), this->m_is.cend());
      --k;

      for (::std::size_t h = this->lg(); h --> 0;) {
        const auto l0 = this->m_bvs[h].rank0(lid);
        const auto r0 = this->m_bvs[h].rank0(rid);
        if (k < r0 - l0) {
          lid = l0;
          rid = r0;
        } else {
          k -= r0 - l0;
          lid += this->m_bvs[h].zeros() - l0;
          rid += this->m_bvs[h].zeros() - r0;
        }
      }

      return ::std::make_pair(this->m_is.cbegin() + lid, this->m_is.cbegin() + rid);
    }

    ::std::pair<typename ::std::vector<::std::size_t>::const_iterator, typename ::std::vector<::std::size_t>::const_iterator> next_points(const T& l, const T& r, const T& d) const {
      assert(this->built());
      assert(l <= r);

      auto lid = this->xid(l);
      auto rid = this->xid(r);
      auto k = this->range_freq(l, r, d);
      if (k == rid - lid) return ::std::make_pair(this->m_is.cend(), this->m_is.cend());

      for (::std::size_t h = this->lg(); h --> 0;) {
        const auto l0 = this->m_bvs[h].rank0(lid);
        const auto r0 = this->m_bvs[h].rank0(rid);
        if (k < r0 - l0) {
          lid = l0;
          rid = r0;
        } else {
          k -= r0 - l0;
          lid += this->m_bvs[h].zeros() - l0;
          rid += this->m_bvs[h].zeros() - r0;
        }
      }

      return ::std::make_pair(this->m_is.cbegin() + lid, this->m_is.cbegin() + rid);
    }
  };
}