ds/offline_query/point_add_rectangle_sum.hpp

• View this file on GitHub
• Last update: 2024-04-09 15:17:41+09:00
• Include: #include "ds/offline_query/point_add_rectangle_sum.hpp"

Depends on

• alg/monoid/add.hpp
• ds/fenwicktree/fenwicktree.hpp

Verified with

• test/library_checker/datastructure/rectangle_sum_sweep.test.cpp
• test/yukicoder/2161.test.cpp

Code

#include "ds/fenwicktree/fenwicktree.hpp"


template <typename AbelGroup, typename XY, bool SMALL_X = false>
struct Point_Add_Rectangle_Sum {
  using G = typename AbelGroup::value_type;
  using Point = tuple<XY, XY, G>;
  vector<Point> point;
  vector<tuple<XY, XY, XY, XY>> rect;

  Point_Add_Rectangle_Sum() {}

  void add_query(XY x, XY y, G w) { point.eb(x, y, w); }
  void sum_query(XY xl, XY xr, XY yl, XY yr) { rect.eb(xl, xr, yl, yr); }

  vector<G> calc() {
    int N = point.size(), Q = rect.size();
    if (N == 0 || Q == 0) return vector<G>(Q, AbelGroup::unit());
    // X 方向の座圧

    int NX = 0;
    if (!SMALL_X) {
      sort(all(point),
           [&](auto &x, auto &y) -> bool { return get<0>(x) < get<0>(y); });
      vc<XY> keyX;
      keyX.reserve(N);
      for (auto &&[a, b, c]: point) {
        if (len(keyX) == 0 || keyX.back() != a) { keyX.eb(a); }
        a = len(keyX) - 1;
      }
      for (auto &&[xl, xr, yl, yr]: rect) {
        xl = LB(keyX, xl);
        xr = LB(keyX, xr);
      }
      NX = len(keyX);
    }
    if (SMALL_X) {
      XY mx = infty<XY>;
      for (auto &&[x, y, g]: point) chmin(mx, x);
      for (auto &&[x, y, g]: point) x -= mx, chmax(NX, x + 1);
      for (auto &&[xl, xr, yl, yr]: rect) {
        xl -= mx, xr -= mx;
        xl = max(0, min<int>(xl, NX));
        xr = max(0, min<int>(xr, NX));
      }
    }

    vc<tuple<XY, int, int, int>> event(Q + Q);
    FOR(q, Q) {
      auto &[xl, xr, yl, yr] = rect[q];
      event[2 * q] = {yl, xl, xr, 2 * q};
      event[2 * q + 1] = {yr, xl, xr, 2 * q + 1};
    }
    sort(all(point),
         [&](auto &x, auto &y) -> bool { return get<1>(x) < get<1>(y); });
    sort(all(event),
         [&](auto &x, auto &y) -> bool { return get<0>(x) < get<0>(y); });

    FenwickTree<AbelGroup> bit(NX);
    vc<G> res(Q, AbelGroup::unit());
    int j = 0;
    for (auto &&[y, xl, xr, qq]: event) {
      while (j < N && get<1>(point[j]) < y) {
        bit.add(get<0>(point[j]), get<2>(point[j]));
        ++j;
      }
      G g = bit.sum(xl, xr);
      int q = qq / 2;
      if (qq % 2 == 0) g = AbelGroup::inverse(g);
      res[q] = AbelGroup::op(res[q], g);
    }
    return res;
  }
};

#line 2 "alg/monoid/add.hpp"

template <typename E>
struct Monoid_Add {
  using X = E;
  using value_type = X;
  static constexpr X op(const X &x, const X &y) noexcept { return x + y; }
  static constexpr X inverse(const X &x) noexcept { return -x; }
  static constexpr X power(const X &x, ll n) noexcept { return X(n) * x; }
  static constexpr X unit() { return X(0); }
  static constexpr bool commute = true;
};
#line 3 "ds/fenwicktree/fenwicktree.hpp"

template <typename Monoid>
struct FenwickTree {
  using G = Monoid;
  using E = typename G::value_type;
  int n;
  vector<E> dat;
  E total;

  FenwickTree() {}
  FenwickTree(int n) { build(n); }
  template <typename F>
  FenwickTree(int n, F f) {
    build(n, f);
  }
  FenwickTree(const vc<E>& v) { build(v); }

  void build(int m) {
    n = m;
    dat.assign(m, G::unit());
    total = G::unit();
  }
  void build(const vc<E>& v) {
    build(len(v), [&](int i) -> E { return v[i]; });
  }
  template <typename F>
  void build(int m, F f) {
    n = m;
    dat.clear();
    dat.reserve(n);
    total = G::unit();
    FOR(i, n) { dat.eb(f(i)); }
    for (int i = 1; i <= n; ++i) {
      int j = i + (i & -i);
      if (j <= n) dat[j - 1] = G::op(dat[i - 1], dat[j - 1]);
    }
    total = prefix_sum(m);
  }

  E prod_all() { return total; }
  E sum_all() { return total; }
  E sum(int k) { return prefix_sum(k); }
  E prod(int k) { return prefix_prod(k); }
  E prefix_sum(int k) { return prefix_prod(k); }
  E prefix_prod(int k) {
    chmin(k, n);
    E ret = G::unit();
    for (; k > 0; k -= k & -k) ret = G::op(ret, dat[k - 1]);
    return ret;
  }
  E sum(int L, int R) { return prod(L, R); }
  E prod(int L, int R) {
    chmax(L, 0), chmin(R, n);
    if (L == 0) return prefix_prod(R);
    assert(0 <= L && L <= R && R <= n);
    E pos = G::unit(), neg = G::unit();
    while (L < R) { pos = G::op(pos, dat[R - 1]), R -= R & -R; }
    while (R < L) { neg = G::op(neg, dat[L - 1]), L -= L & -L; }
    return G::op(pos, G::inverse(neg));
  }

  vc<E> get_all() {
    vc<E> res(n);
    FOR(i, n) res[i] = prod(i, i + 1);
    return res;
  }

  void add(int k, E x) { multiply(k, x); }
  void multiply(int k, E x) {
    static_assert(G::commute);
    total = G::op(total, x);
    for (++k; k <= n; k += k & -k) dat[k - 1] = G::op(dat[k - 1], x);
  }

  template <class F>
  int max_right(const F check, int L = 0) {
    assert(check(G::unit()));
    E s = G::unit();
    int i = L;
    // 2^k 進むとダメ
    int k = [&]() {
      while (1) {
        if (i % 2 == 1) { s = G::op(s, G::inverse(dat[i - 1])), i -= 1; }
        if (i == 0) { return topbit(n) + 1; }
        int k = lowbit(i) - 1;
        if (i + (1 << k) > n) return k;
        E t = G::op(s, dat[i + (1 << k) - 1]);
        if (!check(t)) { return k; }
        s = G::op(s, G::inverse(dat[i - 1])), i -= i & -i;
      }
    }();
    while (k) {
      --k;
      if (i + (1 << k) - 1 < len(dat)) {
        E t = G::op(s, dat[i + (1 << k) - 1]);
        if (check(t)) { i += (1 << k), s = t; }
      }
    }
    return i;
  }

  // check(i, x)
  template <class F>
  int max_right_with_index(const F check, int L = 0) {
    assert(check(L, G::unit()));
    E s = G::unit();
    int i = L;
    // 2^k 進むとダメ
    int k = [&]() {
      while (1) {
        if (i % 2 == 1) { s = G::op(s, G::inverse(dat[i - 1])), i -= 1; }
        if (i == 0) { return topbit(n) + 1; }
        int k = lowbit(i) - 1;
        if (i + (1 << k) > n) return k;
        E t = G::op(s, dat[i + (1 << k) - 1]);
        if (!check(i + (1 << k), t)) { return k; }
        s = G::op(s, G::inverse(dat[i - 1])), i -= i & -i;
      }
    }();
    while (k) {
      --k;
      if (i + (1 << k) - 1 < len(dat)) {
        E t = G::op(s, dat[i + (1 << k) - 1]);
        if (check(i + (1 << k), t)) { i += (1 << k), s = t; }
      }
    }
    return i;
  }

  template <class F>
  int min_left(const F check, int R) {
    assert(check(G::unit()));
    E s = G::unit();
    int i = R;
    // false になるところまで戻る
    int k = 0;
    while (i > 0 && check(s)) {
      s = G::op(s, dat[i - 1]);
      k = lowbit(i);
      i -= i & -i;
    }
    if (check(s)) {
      assert(i == 0);
      return 0;
    }
    // 2^k 進むと ok になる
    // false を維持して進む
    while (k) {
      --k;
      E t = G::op(s, G::inverse(dat[i + (1 << k) - 1]));
      if (!check(t)) { i += (1 << k), s = t; }
    }
    return i + 1;
  }

  int kth(E k, int L = 0) {
    return max_right([&k](E x) -> bool { return x <= k; }, L);
  }
};
#line 2 "ds/offline_query/point_add_rectangle_sum.hpp"

template <typename AbelGroup, typename XY, bool SMALL_X = false>
struct Point_Add_Rectangle_Sum {
  using G = typename AbelGroup::value_type;
  using Point = tuple<XY, XY, G>;
  vector<Point> point;
  vector<tuple<XY, XY, XY, XY>> rect;

  Point_Add_Rectangle_Sum() {}

  void add_query(XY x, XY y, G w) { point.eb(x, y, w); }
  void sum_query(XY xl, XY xr, XY yl, XY yr) { rect.eb(xl, xr, yl, yr); }

  vector<G> calc() {
    int N = point.size(), Q = rect.size();
    if (N == 0 || Q == 0) return vector<G>(Q, AbelGroup::unit());
    // X 方向の座圧

    int NX = 0;
    if (!SMALL_X) {
      sort(all(point),
           [&](auto &x, auto &y) -> bool { return get<0>(x) < get<0>(y); });
      vc<XY> keyX;
      keyX.reserve(N);
      for (auto &&[a, b, c]: point) {
        if (len(keyX) == 0 || keyX.back() != a) { keyX.eb(a); }
        a = len(keyX) - 1;
      }
      for (auto &&[xl, xr, yl, yr]: rect) {
        xl = LB(keyX, xl);
        xr = LB(keyX, xr);
      }
      NX = len(keyX);
    }
    if (SMALL_X) {
      XY mx = infty<XY>;
      for (auto &&[x, y, g]: point) chmin(mx, x);
      for (auto &&[x, y, g]: point) x -= mx, chmax(NX, x + 1);
      for (auto &&[xl, xr, yl, yr]: rect) {
        xl -= mx, xr -= mx;
        xl = max(0, min<int>(xl, NX));
        xr = max(0, min<int>(xr, NX));
      }
    }

    vc<tuple<XY, int, int, int>> event(Q + Q);
    FOR(q, Q) {
      auto &[xl, xr, yl, yr] = rect[q];
      event[2 * q] = {yl, xl, xr, 2 * q};
      event[2 * q + 1] = {yr, xl, xr, 2 * q + 1};
    }
    sort(all(point),
         [&](auto &x, auto &y) -> bool { return get<1>(x) < get<1>(y); });
    sort(all(event),
         [&](auto &x, auto &y) -> bool { return get<0>(x) < get<0>(y); });

    FenwickTree<AbelGroup> bit(NX);
    vc<G> res(Q, AbelGroup::unit());
    int j = 0;
    for (auto &&[y, xl, xr, qq]: event) {
      while (j < N && get<1>(point[j]) < y) {
        bit.add(get<0>(point[j]), get<2>(point[j]));
        ++j;
      }
      G g = bit.sum(xl, xr);
      int q = qq / 2;
      if (qq % 2 == 0) g = AbelGroup::inverse(g);
      res[q] = AbelGroup::op(res[q], g);
    }
    return res;
  }
};

Back to top page