ei1333の日記

ぺこい

Point/Rectangle Add Rectangle Sum

はい

実は ei1333 の日記はデータ構造の記事が少なかったりします

ライブラリにしておこうね！

Static Point Add Rectangle Sum

judge.yosupo.jp

問題設定

{W \times H} のマス目があります。左から {x (0 \leq x \lt W)} 列目、下から {y (0 \leq y \lt H)} 行目のマスを {(x, y)} で表します。

最初に、次の {N} 個の Add クエリが与えられます。

  • {Add(x, y, w)}: マス {(x, y)} の重みに {w} を加える

このとき、次の {Q} 個の Sum クエリを処理したいです。

  • {Sum(l, d, r, u)}: {l \leq x \lt r} かつ {d \leq y \lt u} を満たす全てのマス {(x, y)} についての重みの総和 を求める

解法

平面走査 + BIT により {O(N \log N + Q \log Q)} で解くことができます。

1. Sum クエリの分解

{Sum(r, u)} を {0 \leq x \lt r} かつ {0 \leq y \lt u} を満たす全てのマス {(x, y)} についての重みの総和 とします。

このとき、 {Sum(l, d, r, u) = Sum(r, u) + Sum(l, d) - Sum(l, u) - Sum(r, d)} が成立します。

{Sum(r, u)} が解ければ良いので 以降では Sum クエリをこの形式で扱います。

2. 寄与をかんがえる

{Sum(r,u)} に対する {Add(x,y,w)} の寄与を考えます。

{x \lt r} かつ {y \lt u} を満たす {Add(x,y,w)} について {w} です。

3. 平面走査

{Add(x, y, w)} と {Sum(r, u)} を {x, r} の昇順に処理します。( {x = r} の場合は Sum クエリを優先します）

{dp(y)} を {0 \leq x \lt r} を満たす全ての {Add(x, y, w)} についての総和とします。

{Add(x, y, w)} では {dp(y)} に {w} を足します。

{Sum(r, u)} では {\displaystyle \sum_{i=0}^{u-1} dp(i)} を求めます。

とけました。

4. BIT による高速化

一点への加算クエリとprefixの総和を求めるクエリを効率的に処理できればよいです。

BITをつかいます。

{O(Q \log H)} です。

ところで、BITを使う関係でこっそり縦方向は座標圧縮しています。

実装例

実装例はこちら！！！！！！　BITはACLを使ってます

template< typename T, typename C >
vector< C > point_add_rectangle_sum(vector< tuple< T, T, C > > add,
                                    vector< tuple< T, T, T, T > > query) {
  static_assert(is_integral< T >::value, "template parameter T must be integral type");

  int N = (int) add.size();
  int Q = (int) query.size();
  using add_t = tuple< T, T, C >;
  using query_t = tuple< T, int >;

  sort(add.begin(), add.end(), [](const add_t &a, const add_t &b) {
    return get< 1 >(a) < get< 1 >(b);
  });
  vector< T > ys;
  ys.reserve(N);
  for(auto&[x, y, w]: add) {
    if(ys.empty() or ys.back() != y) ys.emplace_back(y);
    y = (int) ys.size() - 1;
  }
  ys.shrink_to_fit();
  vector< query_t > qs;
  qs.reserve(2 * Q);
  for(int i = 0; i < Q; i++) {
    auto&[l, d, r, u] = query[i];
    d = lower_bound(ys.begin(), ys.end(), d) - ys.begin();
    u = lower_bound(ys.begin(), ys.end(), u) - ys.begin();
    qs.emplace_back(l, i);
    qs.emplace_back(r, i + Q);
  }
  sort(qs.begin(), qs.end(), [](const query_t &a, const query_t &b) {
    return get< 0 >(a) < get< 0 >(b);
  });
  sort(add.begin(), add.end(), [](const add_t &a, const add_t &b) {
    return get< 0 >(a) < get< 0 >(b);
  });
  vector< C > ans(Q);
  int j = 0;
  atcoder::fenwick_tree< C > bit(ys.size());
  for(auto&[x, i]: qs) {
    while(j < N and get< 0 >(add[j]) < x) {
      bit.add(get< 1 >(add[j]), get< 2 >(add[j]));
      ++j;
    }
    if(i < Q) {
      ans[i] -= bit.sum(get< 1 >(query[i]), get< 3 >(query[i]));
    } else {
      i -= Q;
      ans[i] += bit.sum(get< 1 >(query[i]), get< 3 >(query[i]));
    }
  }
  return ans;
}

Static Rectangle Add Rectangle Sum

Point Add ではなく Rectangle Add の場合を考えます。

問題設定

{W \times H} のマス目があります。左から {x (0 \leq x \lt W)} 列目、下から {y (0 \leq y \lt H)} 行目のマスを {(x, y)} で表します。

最初に、次の {N} 個の Rectangle Add クエリが与えられます。

  • {Add(l, d, r, u, w)}: {l \leq x \lt r} かつ {d \leq y \lt u} を満たす全てのマス {(x, y)} について重み {w} を加える

このとき、次の {Q} 個の Rectangle Sum クエリに答えたいです。

  • {Sum(l, d, r, u)}: {l \leq x \lt r} かつ {d \leq y \lt u} を満たす全てのマス {(x, y)} についての重みの総和 を求める

解法

Point Add Rectangle Sum と同様の考え方で解くことができます。

1. Sum クエリと Add クエリの分解

{Sum(r, u)} を {0 \leq x \lt r} かつ {0 \leq y \lt u} を満たす全てのマス {(x, y)} についての重みの総和 とします。

このとき、 {Sum(l, d, r, u) = Sum(r, u) + Sum(l, d) - Sum(l, u) - Sum(r, d)} が成立します。（Point Add と全く同じ）

{Add(l, d, w)} を {l \leq x \lt W} かつ {r \leq y \lt H} を満たす全てのマス {(x, y)} について {w} を加算するクエリ とします。

このとき、 {Add(l, d, r, u, w)} は {Add(l, d, w)}, {Add(r, u, w)}, {Add(l, u, -w)}, {Add(r, d, -w)} をそれぞれ処理することと同じです。

{Add(l, d, w)} と {Sum(r, u)} を処理できれば良いので、以降ではクエリをこの形式で扱います。

2. 寄与をかんがえる

{Sum(r, u)} に対する {Add(l, d, w)} の寄与を考えます。

{l \lt r} かつ {d \lt u} を満たす {Add(l, d, w)} について {w(r - l)(u - d)} です。

3. 平面走査

{Add(l, d, w)} と {Sum(r, u)} を {l, r} の昇順に処理します。( {l = r} の場合は Sum クエリを優先します）

{dp(d)} を {0 \leq l \lt r} を満たす全ての {Add(l,d,w)} についての {w(r - l)(u - d)} の総和とします。

{Add(l, d, w)} では {dp(d)} に {w(r - l)(u - d)} を足します。

{Sum(r, u)} では {\displaystyle \sum_{i=0}^{u-1} dp(i)} を求めます。

解けました。

4. 変数分離

よく見ると {dp(d)} に {w(r - l)(u - d)} を足すのは、このあとに来る Sum クエリの {r, u} に依存しているため不可能です。

{w(r - l)(u - d) = w(ru + ld - rd - lu)} です。

{Add(l, d, w)} に依存する係数 {l, d, w} と {Sum(r, u)} に依存する係数 {r, u} を分けて考えます。

以下の総和が Sum クエリの答えです。

  • Add クエリで {wld} を足して、Sum クエリで {1} を掛けた値
  • Add クエリで {w} を足して、Sum クエリで {ru} を掛けた値
  • Add クエリで {-wd} を足して、Sum クエリで {r} を掛けた値
  • Add クエリで {-wl} を足して、Sum クエリで {u} を掛けた値

4個それぞれで BIT を持てば良いです。

{O(Q \log H)} です。（Hは {y} の種類数）

実装例

実装例はこちら！！！！！！！！

コードを見ていただけるとわかるんですが、定数倍が香ばしいことになってます

template< typename T, typename C >
vector< C > rectangle_add_rectangle_sum(vector< tuple< T, T, T, T, C > > add,
                                        vector< tuple< T, T, T, T > > query) {
  static_assert(is_integral< T >::value, "template parameter T must be integral type");

  int N = (int) add.size();
  int Q = (int) query.size();
  using query_t = tuple< T, int >;

  vector< T > ys;
  ys.reserve(2 * N);
  for(auto&[l, d, r, u, w]: add) {
    ys.emplace_back(d);
    ys.emplace_back(u);
  }
  sort(ys.begin(), ys.end());
  ys.erase(unique(ys.begin(), ys.end()), ys.end());

  vector< query_t > ps;
  ps.reserve(2 * N);
  for(int i = 0; i < N; i++) {
    auto&[l, d, r, u, w] = add[i];
    d = lower_bound(ys.begin(), ys.end(), d) - ys.begin();
    u = lower_bound(ys.begin(), ys.end(), u) - ys.begin();
    ps.emplace_back(l, i);
    ps.emplace_back(r, i + N);
  }
  vector< query_t > qs;
  qs.reserve(2 * Q);
  vector< int > did(Q), uid(Q);
  for(int i = 0; i < Q; i++) {
    auto&[l, d, r, u] = query[i];
    did[i] = lower_bound(ys.begin(), ys.end(), d) - ys.begin();
    uid[i] = lower_bound(ys.begin(), ys.end(), u) - ys.begin();
    qs.emplace_back(l, i);
    qs.emplace_back(r, i + Q);
  }
  sort(ps.begin(), ps.end(), [](const query_t &a, const query_t &b) {
    return get< 0 >(a) < get< 0 >(b);
  });
  sort(qs.begin(), qs.end(), [](const query_t &a, const query_t &b) {
    return get< 0 >(a) < get< 0 >(b);
  });
  vector< C > ans(Q);
  int j = 0;
  atcoder::fenwick_tree< C > bit1(ys.size());
  atcoder::fenwick_tree< C > bit2(ys.size());
  atcoder::fenwick_tree< C > bit3(ys.size());
  atcoder::fenwick_tree< C > bit4(ys.size());
  for(auto&[x, i]: qs) {
    while(j < 2 * N and get< 0 >(ps[j]) < x) {
      int k = get< 1 >(ps[j]);
      T l, d1, d2;
      C w;
      if(k < N) {
        tie(l, d1, ignore, d2, w) = add[k];
      } else {
        tie(ignore, d2, l, d1, w) = add[k - N];
      }
      bit1.add(d1, w * l * ys[d1]);
      bit2.add(d1, w);
      bit3.add(d1, -w * ys[d1]);
      bit4.add(d1, -w * l);
      bit1.add(d2, -w * l * ys[d2]);
      bit2.add(d2, -w);
      bit3.add(d2, w * ys[d2]);
      bit4.add(d2, w * l);
      ++j;
    }
    T r, d, u;
    if(i < Q) {
      tie(r, d, ignore, u) = query[i];
      ans[i] -= bit1.sum(0, uid[i]);
      ans[i] -= bit2.sum(0, uid[i]) * r * u;
      ans[i] -= bit3.sum(0, uid[i]) * r;
      ans[i] -= bit4.sum(0, uid[i]) * u;
      ans[i] += bit1.sum(0, did[i]);
      ans[i] += bit2.sum(0, did[i]) * r * d;
      ans[i] += bit3.sum(0, did[i]) * r;
      ans[i] += bit4.sum(0, did[i]) * d;
    } else {
      i -= Q;
      tie(ignore, d, r, u) = query[i];
      ans[i] += bit1.sum(0, uid[i]);
      ans[i] += bit2.sum(0, uid[i]) * r * u;
      ans[i] += bit3.sum(0, uid[i]) * r;
      ans[i] += bit4.sum(0, uid[i]) * u;
      ans[i] -= bit1.sum(0, did[i]);
      ans[i] -= bit2.sum(0, did[i]) * r * d;
      ans[i] -= bit3.sum(0, did[i]) * r;
      ans[i] -= bit4.sum(0, did[i]) * d;
    }
  }
  return ans;
}

BITを4つもつよりもBITに長さ4の配列をのせたほうがはやいかもしれん（だれか計測してくれ～）

本質

Static Point Add Rectangle Sum と見比べるとわかるように、Add クエリはマス {(l,d)} にそれぞれ {wld, w, wd, wl} を足すクエリを処理することと等価です。

本質2

平面上にいくつかの点があります。ある点 {(x,y)} が与えられたときに、左下にある全ての点 {(a, b)} {(a \lt x, b \lt y)} に対する寄与の総和を求めたいです。各点の寄与が {x,y,a,b} を用いた式で書けるときは変数分離して平面走査 + BIT で解けます。

Dynamic Add Rectangle Sum

加算クエリが動的に与えられる問題を考えます。

ここでは特に Dynamic Point Add Rectangle Sum を考えます。

Rectangle Add も Point Add に帰着させれば等価です。

judge.yosupo.jp

問題設定

{W \times H} のマス目があります。左から {x (0 \leq x \lt W)} 列目、下から {y (0 \leq y \lt H)} 行目のマスを {(x, y)} で表します。

次の {Q} 個のクエリを処理したいです。

  • {Add(x, y, w)}: マス {(x, y)} の重みに {w} を加える
  • {Sum(l, d, r, u)}: {l \leq x \lt r} かつ {d \leq y \lt u} を満たす全てのマス {(x, y)} についての重みの総和 を求める

解法

分割統治すると Static Add になるため {O(Q \log Q \log H)} で解けます（完）

分割統治

{solve(l, r)} を時刻 {[l, r)} に与えられた Sum クエリを処理する関数としましょう。ただしこのときに考慮する Add クエリは時刻 {[l, r)} で与えられたクエリのみとします。

{\displaystyle m = \frac {l + r} {2}} として {solve(l, m), solve(m, r)} を呼び出します。

このとき 時刻 {[l, m), [m, r)} 内で完結するクエリは全て処理済みです。

残りの分は、時刻 {[l, m)} に与えられた Add クエリの寄与を 時刻 {[m,r)} に与えられた Sum クエリに足すことです。

独立しているため Static Add Rectangle Sum を解けばよいです。

実装例

template< typename T, typename C >
vector< C > dynamic_point_add_rectangle_sum(
    vector< variant< tuple< T, T, C >, tuple< T, T, T, T > > > query) {
  using add_t = tuple< T, T, C >;
  using query_t = tuple< T, T, T, T >;
  int Q = (int) query.size();
  queue< pair< int, int > > que;
  que.emplace(0, Q);
  vector< int > rev(Q);
  int ptr = 0;
  for(int i = 0; i < Q; i++) {
    if(holds_alternative< query_t >(query[i])) {
      rev[i] = ptr++;
    }
  }
  vector< C > ans(ptr);
  while(not que.empty()) {
    auto[l, r] = que.front();
    que.pop();
    if(l + 1 >= r) continue;
    int m = (l + r) / 2;
    que.emplace(l, m);
    que.emplace(m, r);
    vector< add_t > add;
    for(int k = l; k < m; k++) {
      if(holds_alternative< add_t >(query[k])) {
        add.emplace_back(get< add_t >(query[k]));
      }
    }
    vector< query_t > sum;
    for(int k = m; k < r; k++) {
      if(holds_alternative< query_t >(query[k])) {
        sum.emplace_back(get< query_t >(query[k]));
      }
    }
    auto sub = point_add_rectangle_sum(add, sum);
    for(int k = m, t = 0; k < r; k++) {
      if(holds_alternative< query_t >(query[k])) {
        ans[rev[k]] += sub[t++];
      }
    }
  }
  return ans;
}

マージソートみたいにすると定数倍改善になったりする？　わかりません

おまけ

オンラインでクエリを処理する2D SegmentTreeなどのデータ構造もありますが、こわいデータ構造わかりません（オンラインで求める必要がある場合があるときはびっくりしながら使いますが、オフラインの場合は分割統治のほうが定数倍が良いと思っています）

なんかおもしろいことを

かきたいとおもったがおもいつかなかった