Preset wavelet matrix (tools/preset_wavelet_matrix.hpp)

• View this file on GitHub
• View document part on GitHub
• Last update: 2025-07-27 03:09:14+09:00
• Include: #include "tools/preset_wavelet_matrix.hpp"

It is a pair of a wavelet matrix and auxiliary data structures, preconfigured for commonly used operations. Given $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, it can update $w_i$ for any $i$ and answer the sum of $w_i$ such that $l \leq x_i < r$ and $d \leq y_i < u$ hold.

License

• CC0

Author

• anqooqie

Constructor

(1) preset_wavelet_matrix<T, U, Updatable> wm();
(2) preset_wavelet_matrix<T, G, Updatable> wm();

• (1)
  • It is identical to preset_wavelet_matrix<T, tools::group::plus<U>, Updatable> wm();
• (2)
  • It creates a two-dimensional plane with $0$ weighted points.
  • The weight of each point is on a given commutative group $G$.
  • If Updatable is true, the weights can be updated even after wm.build() has been called in exchange for worse time complexity.

Constraints

• (1)
  • tools::is_group_v<U> does not hold.
• (2)
  • tools::is_group_v<G> holds.
  • For all $x$ in typename G::T and $y$ in typename G::T, G::op(x, y) $=$ G::op(y, x).
  • For all $x$ in typename G::T, $y$ in typename G::T and $z$ in typename G::T, G::op(G::op(x, y), z) $=$ G::op(x, G::op(y, z)).
  • For all $x$ in typename G::T, G::op(x, G::e()) $= x$.
  • For all $x$ in typename G::T, G::op(x, G::inv(x)) $=$ G::e().

Time Complexity

• $O(1)$

size

int wm.size();

It returns the number of weighted points on the plane.

Constraints

• None

Time Complexity

• $O(1)$

add_point

int wm.add_point(T x, T y, typename G::T w);

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

Constraints

• wm.build() has not been called ever.

Time Complexity

• $O(1)$ amortized

get_point

std::tuple<T, T, typename G::T> wm.get_point(int i);

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

Constraints

• $0 \leq i < n$ where $n$ is wm.size()

Time Complexity

• $O(1)$

points

std::vector<std::tuple<T, T, typename G::T>> wm.points();

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

Constraints

• None

Time Complexity

• $O(n)$

build

void wm.build();

It internally creates the data structure called as wavelet matrix.

Constraints

• wm.build() has not been called ever.

Time Complexity

• $O(n \log n)$

set_weight

void wm.set_weight(int i, typename G::T w);

It updates the weight of the $i$-th ($0$ indexed) weighted point to $w$.

Constraints

• <Updatable> is true.
• $0 \leq i < n$

Time Complexity

• (Before wm.build() has been called): $O(1)$
• (After wm.build() has been called): $O\left((\log n)^2\right)$

kth_smallest

int wm.kth_smallest(T l, T r, int 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$
• $0 \leq k < |\{i \mid l \leq x_i < r\}|$

Time Complexity

• $O(\log n)$

kth_largest

int wm.kth_largest(T l, T r, int 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$
• $0 \leq k < |\{i \mid l \leq x_i < r\}|$

Time Complexity

• $O(\log n)$

range_sum

(1) typename G::T wm.range_sum(T l, T r, T u);
(2) typename G::T wm.range_sum(T l, T r, T d, T u);

• (1)
  • It returns the sum of $w_i$ such that $l \leq x_i < r$ and $y_i < u$.
• (2)
  • It returns the sum of $w_i$ such that $l \leq x_i < r$ and $d \leq y_i < u$.

Note that the addition is defined by the commutative group $G$.

Constraints

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

Time Complexity

• (<Updatable> is true): $O\left((\log n)^2\right)$
• (<Updatable> is false): $O(\log n)$

range_freq

(1) int wm.range_freq(T l, T r);
(2) int wm.range_freq(T l, T r, T u);
(3) int 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::ranges::subrange<std::vector<int>::const_iterator> prev_points(T l, T r, T u)

It returns $\{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::ranges::subrange<std::vector<int>::const_iterator> next_points(T l, T r, T d)

It returns $\{i \mid l \leq x_i < r \land y_i = y\}$ sorted in ascending order, where $y$ is wm.next_value(l, r, d).

Constraints

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

Time Complexity

• $O(\log n)$

Depends on

• Bit vector (tools/bit_vector.hpp)
• Bit width required to represent a given integer (tools/bit_width.hpp)
• $\left\lfloor \log_2(x) \right\rfloor$ (tools/floor_log2.hpp)
• Typical groups (tools/group.hpp)
• std::is_arithmetic extended for tools::int128_t and tools::uint128_t (tools/is_arithmetic.hpp)
• Check whether T is a group (tools/is_group.hpp)
• std::is_integral extended for tools::int128_t and tools::uint128_t (tools/is_integral.hpp)
• std::is_signed extended for tools::int128_t and tools::uint128_t (tools/is_signed.hpp)
• std::less by first (tools/less_by_first.hpp)
• std::ranges::lower_bound, but returns index (tools/lower_bound.hpp)
• std::make_unsigned extended for tools::int128_t and tools::uint128_t (tools/make_unsigned.hpp)
• Wavelet matrix (tools/wavelet_matrix.hpp)

Verified with

• tests/preset_wavelet_matrix/static.test.cpp
• tests/preset_wavelet_matrix/updatable.test.cpp

Code

#ifndef TOOLS_PRESET_WAVELET_MATRIX_HPP
#define TOOLS_PRESET_WAVELET_MATRIX_HPP

#include <cassert>
#include <optional>
#include <ranges>
#include <tuple>
#include <type_traits>
#include <vector>
#include "atcoder/segtree.hpp"
#include "tools/group.hpp"
#include "tools/is_group.hpp"
#include "tools/wavelet_matrix.hpp"

namespace tools {
  template <typename T, typename X, bool Updatable>
  class preset_wavelet_matrix {
    using G = ::std::conditional_t<::tools::is_group_v<X>, X, ::tools::group::plus<X>>;
    using U = typename G::T;

    ::tools::wavelet_matrix<T> m_wm;
    ::std::vector<U> m_ws;
    ::std::vector<::std::conditional_t<Updatable, ::atcoder::segtree<U, G::op, G::e>, ::std::vector<U>>> m_aux;

    bool built() const {
      return !this->m_aux.empty();
    }

  public:
    preset_wavelet_matrix() = default;

    int size() const {
      return this->m_wm.size();
    }

    int add_point(const T& x, const T& y, const U& w) {
      assert(!this->built());

      this->m_ws.push_back(w);
      return this->m_wm.add_point(x, y);
    }

    ::std::tuple<T, T, U> get_point(const int i) const {
      assert(0 <= i && i < this->size());

      const auto [x, y] = this->m_wm.get_point(i);
      return {x, y, this->m_ws[i]};
    }

    ::std::vector<::std::tuple<T, T, U>> points() const {
      ::std::vector<::std::tuple<T, T, U>> res;
      res.reserve(this->size());
      for (int i = 0; i < this->size(); ++i) {
        res.push_back(this->get_point(i));
      }
      return res;
    }

    void build() {
      assert(!this->built());

      for (const auto& A_h : this->m_wm.build()) {
        if constexpr (Updatable) {
          ::std::vector<U> v(this->size());
          for (int j = 0; j < this->size(); ++j) {
            v[j] = this->m_ws[A_h[j]];
          }
          this->m_aux.emplace_back(v);
        } else {
          this->m_aux.emplace_back();
          this->m_aux.back().reserve(this->size() + 1);
          this->m_aux.back().push_back(G::e());
          for (int j = 0; j < this->size(); ++j) {
            this->m_aux.back().push_back(G::op(this->m_aux.back().back(), this->m_ws[A_h[j]]));
          }
        }
      }
    }

    template <bool SFINAE = Updatable> requires (SFINAE)
    void set_weight(const int i, const U& w) {
      assert(0 <= i && i < this->size());

      this->m_ws[i] = w;
      if (this->built()) {
        for (const auto& [h, j] : this->m_wm.apply(i)) {
          this->m_aux[h].set(j, w);
        }
      }
    }

    int kth_smallest(const T& l, const T& r, const int k) const {
      assert(this->built());
      assert(l <= r);
      assert(k >= 0);

      return this->m_wm.kth_smallest(l, r, k);
    }

    int kth_largest(const T& l, const T& r, const int k) const {
      assert(this->built());
      assert(l <= r);
      assert(k >= 0);

      return this->m_wm.kth_largest(l, r, k);
    }

    U range_sum(const T& l, const T& r, const T& u) const {
      assert(this->built());
      assert(l <= r);

      U res = G::e();
      for (const auto& [h, jl, jr] : this->m_wm.range_prod(l, r, u)) {
        if constexpr (Updatable) {
          res = G::op(res, this->m_aux[h].prod(jl, jr));
        } else {
          res = G::op(res, G::op(G::inv(this->m_aux[h][jl]), this->m_aux[h][jr]));
        }
      }

      return res;
    }

    U range_sum(const T& l, const T& r, const T& d, const T& u) const {
      assert(this->built());
      assert(l <= r);
      assert(d <= u);

      return G::op(this->range_sum(l, r, u), G::inv(this->range_sum(l, r, d)));
    }

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

      return this->m_wm.range_freq(l, r);
    }

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

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

    int 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->m_wm.range_freq(l, r, d, u);
    }

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

      return this->m_wm.prev_value(l, r, u);
    }

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

      return this->m_wm.next_value(l, r, d);
    }

    ::std::ranges::subrange<::std::vector<int>::const_iterator> prev_points(const T& l, const T& r, const T& u) const {
      assert(this->built());
      assert(l <= r);

      return this->m_wm.prev_points(l, r, u);
    }

    ::std::ranges::subrange<::std::vector<int>::const_iterator> next_points(const T& l, const T& r, const T& d) const {
      assert(this->built());
      assert(l <= r);

      return this->m_wm.next_points(l, r, d);
    }
  };
}

#endif

#line 1 "tools/preset_wavelet_matrix.hpp"



#include <cassert>
#include <optional>
#include <ranges>
#include <tuple>
#include <type_traits>
#include <vector>
#line 1 "lib/ac-library/atcoder/segtree.hpp"



#include <algorithm>
#line 6 "lib/ac-library/atcoder/segtree.hpp"
#include <functional>
#line 8 "lib/ac-library/atcoder/segtree.hpp"

#line 1 "lib/ac-library/atcoder/internal_bit.hpp"



#ifdef _MSC_VER
#include <intrin.h>
#endif

#if __cplusplus >= 202002L
#include <bit>
#endif

namespace atcoder {

namespace internal {

#if __cplusplus >= 202002L

using std::bit_ceil;

#else

// @return same with std::bit::bit_ceil
unsigned int bit_ceil(unsigned int n) {
    unsigned int x = 1;
    while (x < (unsigned int)(n)) x *= 2;
    return x;
}

#endif

// @param n `1 <= n`
// @return same with std::bit::countr_zero
int countr_zero(unsigned int n) {
#ifdef _MSC_VER
    unsigned long index;
    _BitScanForward(&index, n);
    return index;
#else
    return __builtin_ctz(n);
#endif
}

// @param n `1 <= n`
// @return same with std::bit::countr_zero
constexpr int countr_zero_constexpr(unsigned int n) {
    int x = 0;
    while (!(n & (1 << x))) x++;
    return x;
}

}  // namespace internal

}  // namespace atcoder


#line 10 "lib/ac-library/atcoder/segtree.hpp"

namespace atcoder {

#if __cplusplus >= 201703L

template <class S, auto op, auto e> struct segtree {
    static_assert(std::is_convertible_v<decltype(op), std::function<S(S, S)>>,
                  "op must work as S(S, S)");
    static_assert(std::is_convertible_v<decltype(e), std::function<S()>>,
                  "e must work as S()");

#else

template <class S, S (*op)(S, S), S (*e)()> struct segtree {

#endif

  public:
    segtree() : segtree(0) {}
    explicit segtree(int n) : segtree(std::vector<S>(n, e())) {}
    explicit segtree(const std::vector<S>& v) : _n(int(v.size())) {
        size = (int)internal::bit_ceil((unsigned int)(_n));
        log = internal::countr_zero((unsigned int)size);
        d = std::vector<S>(2 * size, e());
        for (int i = 0; i < _n; i++) d[size + i] = v[i];
        for (int i = size - 1; i >= 1; i--) {
            update(i);
        }
    }

    void set(int p, S x) {
        assert(0 <= p && p < _n);
        p += size;
        d[p] = x;
        for (int i = 1; i <= log; i++) update(p >> i);
    }

    S get(int p) const {
        assert(0 <= p && p < _n);
        return d[p + size];
    }

    S prod(int l, int r) const {
        assert(0 <= l && l <= r && r <= _n);
        S sml = e(), smr = e();
        l += size;
        r += size;

        while (l < r) {
            if (l & 1) sml = op(sml, d[l++]);
            if (r & 1) smr = op(d[--r], smr);
            l >>= 1;
            r >>= 1;
        }
        return op(sml, smr);
    }

    S all_prod() const { return d[1]; }

    template <bool (*f)(S)> int max_right(int l) const {
        return max_right(l, [](S x) { return f(x); });
    }
    template <class F> int max_right(int l, F f) const {
        assert(0 <= l && l <= _n);
        assert(f(e()));
        if (l == _n) return _n;
        l += size;
        S sm = e();
        do {
            while (l % 2 == 0) l >>= 1;
            if (!f(op(sm, d[l]))) {
                while (l < size) {
                    l = (2 * l);
                    if (f(op(sm, d[l]))) {
                        sm = op(sm, d[l]);
                        l++;
                    }
                }
                return l - size;
            }
            sm = op(sm, d[l]);
            l++;
        } while ((l & -l) != l);
        return _n;
    }

    template <bool (*f)(S)> int min_left(int r) const {
        return min_left(r, [](S x) { return f(x); });
    }
    template <class F> int min_left(int r, F f) const {
        assert(0 <= r && r <= _n);
        assert(f(e()));
        if (r == 0) return 0;
        r += size;
        S sm = e();
        do {
            r--;
            while (r > 1 && (r % 2)) r >>= 1;
            if (!f(op(d[r], sm))) {
                while (r < size) {
                    r = (2 * r + 1);
                    if (f(op(d[r], sm))) {
                        sm = op(d[r], sm);
                        r--;
                    }
                }
                return r + 1 - size;
            }
            sm = op(d[r], sm);
        } while ((r & -r) != r);
        return 0;
    }

  private:
    int _n, size, log;
    std::vector<S> d;

    void update(int k) { d[k] = op(d[2 * k], d[2 * k + 1]); }
};

}  // namespace atcoder


#line 1 "tools/group.hpp"



#include <cstddef>
#line 1 "tools/is_arithmetic.hpp"



#line 5 "tools/is_arithmetic.hpp"

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

  template <typename T>
  inline constexpr bool is_arithmetic_v = ::tools::is_arithmetic<T>::value;
}


#line 7 "tools/group.hpp"

namespace tools {
  namespace group {
    template <typename G>
    class bit_xor {
      using VR = ::std::conditional_t<::tools::is_arithmetic_v<G> && sizeof(G) <= sizeof(::std::size_t), const G, const G&>;

    public:
      using T = G;
      static T op(VR x, VR y) {
        return x ^ y;
      }
      static T e() {
        return T(0);
      }
      static T inv(VR x) {
        return x;
      }
    };

    template <typename G>
    class multiplies {
      using VR = ::std::conditional_t<::tools::is_arithmetic_v<G> && sizeof(G) <= sizeof(::std::size_t), const G, const G&>;

    public:
      using T = G;
      static T op(VR x, VR y) {
        return x * y;
      }
      static T e() {
        return T(1);
      }
      static T inv(VR x) {
        return e() / x;
      }
    };

    template <typename G>
    class plus {
      using VR = ::std::conditional_t<::tools::is_arithmetic_v<G> && sizeof(G) <= sizeof(::std::size_t), const G, const G&>;

    public:
      using T = G;
      static T op(VR x, VR y) {
        return x + y;
      }
      static T e() {
        return T(0);
      }
      static T inv(VR x) {
        return -x;
      }
    };
  }
}


#line 1 "tools/is_group.hpp"



#line 5 "tools/is_group.hpp"
#include <utility>

namespace tools {

  template <typename G, typename = void>
  struct is_group : ::std::false_type {};

  template <typename G>
  struct is_group<G, ::std::enable_if_t<
    ::std::is_same_v<typename G::T, decltype(G::op(::std::declval<typename G::T>(), ::std::declval<typename G::T>()))> &&
    ::std::is_same_v<typename G::T, decltype(G::e())> &&
    ::std::is_same_v<typename G::T, decltype(G::inv(::std::declval<typename G::T>()))>
  , void>> : ::std::true_type {};

  template <typename G>
  inline constexpr bool is_group_v = ::tools::is_group<G>::value;
}


#line 1 "tools/wavelet_matrix.hpp"



#line 5 "tools/wavelet_matrix.hpp"
#include <array>
#line 7 "tools/wavelet_matrix.hpp"
#include <iterator>
#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/floor_log2.hpp"



#line 1 "tools/bit_width.hpp"



#include <bit>
#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 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/lower_bound.hpp"



#line 8 "tools/lower_bound.hpp"

namespace tools {

  template <
    ::std::random_access_iterator I,
    ::std::sentinel_for<I> S,
    typename Proj = ::std::identity,
    typename T = ::std::remove_cvref_t<::std::invoke_result_t<Proj&, ::std::iter_value_t<I>&>>,
    ::std::indirect_strict_weak_order<const T*, ::std::projected<I, Proj>> Comp = ::std::ranges::less
  >
  constexpr ::std::iter_difference_t<I> lower_bound(const I first, const S last, const T& value, const Comp comp = {}, const Proj proj = {}) {
    return ::std::ranges::distance(first, ::std::ranges::lower_bound(first, last, value, comp, proj));
  }

  template <
    ::std::ranges::random_access_range R,
    typename Proj = ::std::identity,
    typename T = ::std::remove_cvref_t<::std::invoke_result_t<Proj&, ::std::ranges::range_value_t<R>&>>,
    ::std::indirect_strict_weak_order<const T*, ::std::projected<::std::ranges::iterator_t<R>, Proj>> Comp = ::std::ranges::less
  >
  constexpr ::std::ranges::range_difference_t<R> lower_bound(R&& r, const T& value, const Comp comp = {}, const Proj proj = {}) {
    return ::std::ranges::distance(::std::ranges::begin(r), ::std::ranges::lower_bound(r, value, comp, proj));
  }
}


#line 17 "tools/wavelet_matrix.hpp"

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

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

  public:
    wavelet_matrix() = default;

    int size() const {
      return this->m_ps.size();
    }

    int 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 int i) const {
      assert(0 <= i && 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<int>> build() {
      assert(!this->built());

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

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

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

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

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

      for (int h = this->lg(); h --> 0;) {
        for (int i = 0; i < n; ++i) {
          if ((cur[i] >> h) & 1) {
            this->m_bvs[h].set(i);
          }
        }
        this->m_bvs[h].build();
        ::std::array<int, 2> offsets = {0, static_cast<int>(this->m_bvs[h].zeros())};
        for (int 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<int, int>> apply(int i) const {
      assert(this->built());
      assert(0 <= i && i < this->size());

      ::std::vector<::std::pair<int, int>> res(this->lg());
      i = this->iid(i);
      for (int 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;
    }

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

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

      assert(0 <= k && k < rid - lid);

      for (int 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];
    }

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

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

      assert(0 <= k && k < rid - lid);

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

    ::std::vector<::std::tuple<int, int, int>> 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<int, int, int>> res(this->lg());
      for (int 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;
    }

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

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

    int 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);

      int res = 0;
      for (int 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;
    }

    int 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::ranges::subrange<::std::vector<int>::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::ranges::subrange(this->m_is.cend(), this->m_is.cend());
      --k;

      for (int 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::ranges::subrange(this->m_is.cbegin() + lid, this->m_is.cbegin() + rid);
    }

    ::std::ranges::subrange<::std::vector<int>::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::ranges::subrange(this->m_is.cend(), this->m_is.cend());

      for (int 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::ranges::subrange(this->m_is.cbegin() + lid, this->m_is.cbegin() + rid);
    }
  };
}


#line 14 "tools/preset_wavelet_matrix.hpp"

namespace tools {
  template <typename T, typename X, bool Updatable>
  class preset_wavelet_matrix {
    using G = ::std::conditional_t<::tools::is_group_v<X>, X, ::tools::group::plus<X>>;
    using U = typename G::T;

    ::tools::wavelet_matrix<T> m_wm;
    ::std::vector<U> m_ws;
    ::std::vector<::std::conditional_t<Updatable, ::atcoder::segtree<U, G::op, G::e>, ::std::vector<U>>> m_aux;

    bool built() const {
      return !this->m_aux.empty();
    }

  public:
    preset_wavelet_matrix() = default;

    int size() const {
      return this->m_wm.size();
    }

    int add_point(const T& x, const T& y, const U& w) {
      assert(!this->built());

      this->m_ws.push_back(w);
      return this->m_wm.add_point(x, y);
    }

    ::std::tuple<T, T, U> get_point(const int i) const {
      assert(0 <= i && i < this->size());

      const auto [x, y] = this->m_wm.get_point(i);
      return {x, y, this->m_ws[i]};
    }

    ::std::vector<::std::tuple<T, T, U>> points() const {
      ::std::vector<::std::tuple<T, T, U>> res;
      res.reserve(this->size());
      for (int i = 0; i < this->size(); ++i) {
        res.push_back(this->get_point(i));
      }
      return res;
    }

    void build() {
      assert(!this->built());

      for (const auto& A_h : this->m_wm.build()) {
        if constexpr (Updatable) {
          ::std::vector<U> v(this->size());
          for (int j = 0; j < this->size(); ++j) {
            v[j] = this->m_ws[A_h[j]];
          }
          this->m_aux.emplace_back(v);
        } else {
          this->m_aux.emplace_back();
          this->m_aux.back().reserve(this->size() + 1);
          this->m_aux.back().push_back(G::e());
          for (int j = 0; j < this->size(); ++j) {
            this->m_aux.back().push_back(G::op(this->m_aux.back().back(), this->m_ws[A_h[j]]));
          }
        }
      }
    }

    template <bool SFINAE = Updatable> requires (SFINAE)
    void set_weight(const int i, const U& w) {
      assert(0 <= i && i < this->size());

      this->m_ws[i] = w;
      if (this->built()) {
        for (const auto& [h, j] : this->m_wm.apply(i)) {
          this->m_aux[h].set(j, w);
        }
      }
    }

    int kth_smallest(const T& l, const T& r, const int k) const {
      assert(this->built());
      assert(l <= r);
      assert(k >= 0);

      return this->m_wm.kth_smallest(l, r, k);
    }

    int kth_largest(const T& l, const T& r, const int k) const {
      assert(this->built());
      assert(l <= r);
      assert(k >= 0);

      return this->m_wm.kth_largest(l, r, k);
    }

    U range_sum(const T& l, const T& r, const T& u) const {
      assert(this->built());
      assert(l <= r);

      U res = G::e();
      for (const auto& [h, jl, jr] : this->m_wm.range_prod(l, r, u)) {
        if constexpr (Updatable) {
          res = G::op(res, this->m_aux[h].prod(jl, jr));
        } else {
          res = G::op(res, G::op(G::inv(this->m_aux[h][jl]), this->m_aux[h][jr]));
        }
      }

      return res;
    }

    U range_sum(const T& l, const T& r, const T& d, const T& u) const {
      assert(this->built());
      assert(l <= r);
      assert(d <= u);

      return G::op(this->range_sum(l, r, u), G::inv(this->range_sum(l, r, d)));
    }

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

      return this->m_wm.range_freq(l, r);
    }

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

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

    int 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->m_wm.range_freq(l, r, d, u);
    }

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

      return this->m_wm.prev_value(l, r, u);
    }

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

      return this->m_wm.next_value(l, r, d);
    }

    ::std::ranges::subrange<::std::vector<int>::const_iterator> prev_points(const T& l, const T& r, const T& u) const {
      assert(this->built());
      assert(l <= r);

      return this->m_wm.prev_points(l, r, u);
    }

    ::std::ranges::subrange<::std::vector<int>::const_iterator> next_points(const T& l, const T& r, const T& d) const {
      assert(this->built());
      assert(l <= r);

      return this->m_wm.next_points(l, r, d);
    }
  };
}

Back to top page