ds/offline_query/rectangle_add_rectangle_sum.hpp

• View this file on GitHub
• Last update: 2025-12-07 20:35:27+09:00
• Include: #include "ds/offline_query/rectangle_add_rectangle_sum.hpp"
• Link: View error logs on GitHub Actions

Depends on

• alg/monoid/add.hpp
• ds/fenwicktree/fenwicktree.hpp
• ds/index_compression.hpp
• ds/offline_query/coeffient_query_2d.hpp

Verified with

• test/2_library_checker/data_structure/static_rectangle_add_rectangle_sum.test.cpp
• test/3_yukicoder/1490.test.cpp

Code

#include "ds/offline_query/coeffient_query_2d.hpp"

template <typename T, typename XY>
struct Rectangle_Add_Rectangle_Sum {
  Coefficient_Query_2D<2, 2, T, XY> CQ;
  void add_query(XY x1, XY x2, XY y1, XY y2, T w) {
    CQ.add_query(x1, y1, w), CQ.add_query(x1, y2, -w);
    CQ.add_query(x2, y1, -w), CQ.add_query(x2, y2, w);
  }
  void sum_query(XY x1, XY x2, XY y1, XY y2) {
    --x1, --y1, --x2, --y2;
    CQ.sum_query(x1, y1), CQ.sum_query(x1, y2);
    CQ.sum_query(x2, y1), CQ.sum_query(x2, y2);
  }
  vc<T> calc() {
    vc<T> tmp = CQ.calc();
    int Q = len(tmp) / 4;
    vc<T> res(Q);
    FOR(q, Q) {
      res[q] = tmp[4 * q] - tmp[4 * q + 1] - tmp[4 * q + 2] + tmp[4 * q + 3];
    }
    return res;
  }
};

#line 1 "ds/index_compression.hpp"
template <typename T>
struct Index_Compression_DISTINCT_SMALL {
  int mi, ma;
  vc<T> dat;
  vc<T> build(vc<int> X) {
    mi = 0, ma = -1;
    if (!X.empty()) mi = MIN(X), ma = MAX(X);
    dat.assign(ma - mi + 2, 0);
    for (auto& x : X) dat[x - mi + 1]++;
    FOR(i, len(dat) - 1) dat[i + 1] += dat[i];
    for (auto& x : X) {
      x = dat[x - mi]++;
    }
    FOR_R(i, 1, len(dat)) dat[i] = dat[i - 1];
    dat[0] = 0;
    return X;
  }
  int size() { return len(dat); }
  int operator()(ll x) { return dat[clamp<ll>(x - mi, 0, ma - mi + 1)]; }
};

template <typename T>
struct Index_Compression_SAME_SMALL {
  int mi, ma;
  vc<T> dat;
  vc<T> build(vc<T> X) {
    mi = 0, ma = -1;
    if (!X.empty()) mi = MIN(X), ma = MAX(X);
    dat.assign(ma - mi + 2, 0);
    for (auto& x : X) dat[x - mi + 1] = 1;
    FOR(i, len(dat) - 1) dat[i + 1] += dat[i];
    for (auto& x : X) {
      x = dat[x - mi];
    }
    return X;
  }
  int size() { return len(dat); }
  int operator()(ll x) { return dat[clamp<ll>(x - mi, 0, ma - mi + 1)]; }
};

template <typename T>
struct Index_Compression_SAME_LARGE {
  vc<T> dat;
  vc<int> build(vc<T> X) {
    vc<int> I = argsort(X);
    vc<int> res(len(X));
    for (auto& i : I) {
      if (!dat.empty() && dat.back() == X[i]) {
        res[i] = len(dat) - 1;
      } else {
        res[i] = len(dat);
        dat.eb(X[i]);
      }
    }
    dat.shrink_to_fit();
    return res;
  }
  int size() { return len(dat); }
  int operator()(T x) { return LB(dat, x); }
};

template <typename T>
struct Index_Compression_DISTINCT_LARGE {
  vc<T> dat;
  vc<int> build(vc<T> X) {
    vc<int> I = argsort(X);
    vc<int> res(len(X));
    for (auto& i : I) {
      res[i] = len(dat), dat.eb(X[i]);
    }
    dat.shrink_to_fit();
    return res;
  }
  int size() { return len(dat); }
  int operator()(T x) { return LB(dat, x); }
};

template <typename T, bool SMALL>
using Index_Compression_DISTINCT =
    typename std::conditional<SMALL, Index_Compression_DISTINCT_SMALL<T>,
                              Index_Compression_DISTINCT_LARGE<T>>::type;
template <typename T, bool SMALL>
using Index_Compression_SAME =
    typename std::conditional<SMALL, Index_Compression_SAME_SMALL<T>,
                              Index_Compression_SAME_LARGE<T>>::type;

// SAME: [2,3,2] -> [0,1,0]
// DISTINCT: [2,2,3] -> [0,2,1]
// build で列を圧縮してくれる. そのあと
// (x): lower_bound(X,x) をかえす
template <typename T, bool SAME, bool SMALL>
using Index_Compression =
    typename std::conditional<SAME, Index_Compression_SAME<T, SMALL>,
                              Index_Compression_DISTINCT<T, SMALL>>::type;
#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 MX = 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);
  }
  void set(int k, E x) { add(k, G::op(G::inverse(prod(k, 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 3 "ds/offline_query/coeffient_query_2d.hpp"

// A, B：定数
// Sparse Laurent Polynomial f(x,y) を与える
// [x^py^q] f(x,y)/(1-x)^A(1-y)^B をたくさん求める
// O(AB N logN) 時間
template <int A, int B, typename T, bool STATIC>
struct Coefficient_Query_2D {
  struct Mono {
    using value_type = array<T, A * B>;
    using X = value_type;
    static X op(X x, X y) {
      FOR(i, A * B) x[i] += y[i];
      return x;
    }
    static constexpr X unit() { return X{}; }
    static constexpr bool commute = 1;
  };
  vc<tuple<ll, ll, T, int>> query;

  int nsum = 0;
  Coefficient_Query_2D() {}
  void add_query(ll x, ll y, T c) {
    if (c != T(0)) query.eb(x, y, c, -1);
  }
  void sum_query(ll p, ll q) { query.eb(p, q, 0, nsum++); }

  // オーバーフローなどの状況次第で書き換える
  template <int n>
  void comb_array(ll x, array<T, n>& S) {
    static_assert(n < 4);
    if constexpr (n == 1) S = {T(1)};
    if constexpr (n == 2) S = {T(1), T(x)};
    if constexpr (n == 3) S = {T(1), T(x), T(x * (x - 1) / 2)};
  }
  template <int n>
  void coef_array(ll b, array<T, n>& S) {
    static_assert(n < 4);
    // [t^x]t^b(1-t)^{-n} を binom(x,k) の線形結合で表す係数
    if constexpr (n == 1) S = {T(1)};
    if constexpr (n == 2) S = {T(1 - b), T(1)};
    if constexpr (n == 3) S = {T((b - 1) * (b - 2) / 2), T(2 - b), T(1)};
  }

  vc<T> ANS;
  bool done = false;
  void calc_static(const vc<int>& ADD_I, vc<int>& GET_I) {
    if (ADD_I.empty() || GET_I.empty()) return;
    Index_Compression<ll, true, false> IY;
    {
      vc<ll> tmp;
      for (int q : ADD_I) {
        auto [a, b, w, qid] = query[q];
        if (qid == -1) tmp.eb(b);
      }
      IY.build(tmp);
    }

    FenwickTree<Mono> bit(len(IY));

    array<T, A> CX;
    array<T, B> CY;
    array<T, A * B> tmp;

    int ptr = 0;
    for (int q : GET_I) {
      auto [a, b, w, qid] = query[q];
      while (ptr < len(ADD_I) && (get<0>(query[ADD_I[ptr]])) <= a) {
        int q = ADD_I[ptr++];
        auto [a, b, w, qid] = query[q];
        coef_array<A>(a, CX);
        coef_array<B>(b, CY);
        FOR(i, A) FOR(j, B) tmp[B * i + j] = CX[i] * CY[j] * w;
        bit.add(IY(b), tmp);
      }
      comb_array<A>(a, CX);
      comb_array<B>(b, CY);
      // calc query
      tmp = bit.prod(IY(b + 1));
      T ans = 0;
      FOR(i, A) FOR(j, B) ans += tmp[B * i + j] * CX[i] * CY[j];
      ANS[qid] += ans;
    }
  }

  vc<T> calc() {
    assert(!done);
    done = 1;
    ANS.assign(nsum, 0);
    int Q = len(query);
    auto comp = [&](int i, int j) -> bool {
      return (get<0>(query[i])) < (get<0>(query[j]));
    };
    if (STATIC) {
      vc<int> ADD, GET;
      FOR(i, Q) { (get<3>(query[i]) == -1 ? ADD : GET).eb(i); }
      sort(all(ADD), comp);
      sort(all(GET), comp);
      calc_static(ADD, GET);
      return ANS;
    }
    auto dfs = [&](auto& dfs, int L, int R) -> pair<vc<int>, vc<int>> {
      vc<int> ADD, GET;
      if (R == L + 1) {
        (get<3>(query[L]) == -1 ? ADD : GET).eb(L);
        return {ADD, GET};
      }
      int M = (L + R) / 2;
      auto [ADD1, GET1] = dfs(dfs, L, M);
      auto [ADD2, GET2] = dfs(dfs, M, R);
      calc_static(ADD1, GET2);
      ADD.resize(len(ADD1) + len(ADD2));
      GET.resize(len(GET1) + len(GET2));
      merge(all(ADD1), all(ADD2), ADD.begin(), comp);
      merge(all(GET1), all(GET2), GET.begin(), comp);
      return {ADD, GET};
    };
    dfs(dfs, 0, Q);
    return ANS;
  }
};
#line 2 "ds/offline_query/rectangle_add_rectangle_sum.hpp"

template <typename T, typename XY>
struct Rectangle_Add_Rectangle_Sum {
  Coefficient_Query_2D<2, 2, T, XY> CQ;
  void add_query(XY x1, XY x2, XY y1, XY y2, T w) {
    CQ.add_query(x1, y1, w), CQ.add_query(x1, y2, -w);
    CQ.add_query(x2, y1, -w), CQ.add_query(x2, y2, w);
  }
  void sum_query(XY x1, XY x2, XY y1, XY y2) {
    --x1, --y1, --x2, --y2;
    CQ.sum_query(x1, y1), CQ.sum_query(x1, y2);
    CQ.sum_query(x2, y1), CQ.sum_query(x2, y2);
  }
  vc<T> calc() {
    vc<T> tmp = CQ.calc();
    int Q = len(tmp) / 4;
    vc<T> res(Q);
    FOR(q, Q) {
      res[q] = tmp[4 * q] - tmp[4 * q + 1] - tmp[4 * q + 2] + tmp[4 * q + 3];
    }
    return res;
  }
};

Back to top page