【作ってみた】二次元 Disjoint Sparse Table【なにこれ】

導入

訳あってSparse Tableについて調べていたのだが、こんな記事を見つけてしまった

この記事では「二次元のRmQクエリを高速に処理をする方法」として書かれているが、演算部分を抽象化することにより一般的な冪等半群のクエリを処理することができる二次元 Sparse Tableとなる。

二次元 Sparse Tableができるなら、二次元 Disjoint Sparse Tableもできるんじゃないか?

ひとまず検索

こういうのはだいたい先人が開発しているものである。別に車輪の再発名をすることに対して抵抗はないが、よりよい実装があるなら見てみたい。とりあえずGoogle先生に聞いてみることとしよう。

「二次元 Disjoint Sparse Table」や「2D Disjoint Sparse Table」など、それっぽい単語を入れてみる


出ない


なぜだ?

というかそもそも、2D Sparse Tableの時点でこどふぉのブログ一つ以外にそれっぽい記事が見当たらない。まぁ確かにこれが必要になるような問題ならみんな二次元セグ木を生やすのかもしれない...

それとも、そもそも二次元に拡張することが困難なのか?そんなことはないように思えるのだが...

書いてみる

仕方ないから1から自分で書くことにする。

上記のブログを参考にまずは2D Sparse Tableを書いてみて感覚をつかむ。これはお手本があるのですんなりと行った。tableの値をマージしていくのが面白かった。

ここからが本題、Disjointの方だ。tableの大きさはさっきと同じで  \log n \times n \times \log m \times m としてみる。

すると、  \rm{table _ {0,i}} i 行目の要素についての 1D Disjoint Sparse Table を生やせばよさそうだ

 i にのループを回して生やしてみる ここまで  O(nm\log m)

そしたら次は生やしたDisjoint Sparse Tableを上手くマージして二次元にする

うーん、これはどうすればいいのだろう? Sparse Tableに比べてマージの仕方が少しだけ複雑なので、脳内だけで考えると結構混乱してしまった。この辺は数学力が足りていないのだろう。

1D Disjoint Sparse Tableで区間長1の部分には引数に渡した元の配列の要素を入れるが、ここにさっき生やしまくったDisjoint Sparse Tableを使えばいいことに気づいた。これは行ける気がする。

とりあえずそれっぽいものが書けた。最後にクエリの処理をする方法を考えよう。

やはり縦横をそれぞれ二つに分けて合計四つの領域に分けるのがいいのだろう。 しかし、Disjoint Sparse Tableは長さ1の区間は別処理をしなければいけない。

いい書き方が思いつかなかったので、丁寧に場合分けを書く。(といってもせいぜい四通りだ)

tableのどの部分を使うかについては、1D Disjoint Sparse Tableと全く同じ方法を使う。xorを取って、最上位bitの位置を求めるための前計算をする。

その前計算を保存しておく配列は、次元は増えても一次元ずつ処理していくので結局は変わらないわけだ。

そんなわけで完成してしまった。

template<typename T>
struct DisjointSparseTable2D{
  using F = function<T(T, T)>;
  const F f;
  vector<vector<vector<vector<T>>>> table;
  vector<int> lookup;

  DisjointSparseTable2D(const vector<vector<T>> &v, const F &f): f(f){
    const int n = v.size(), m = v[0].size();
    const int nb = 32 - __builtin_clz(n), mb = 32 - __builtin_clz(m);
    table.assign(nb, vector(n, vector(mb, vector<T>(m))));
    for(int rj = 0; rj < n; rj++){

      for(int cj = 0; cj < m; cj++) table[0][rj][0][cj] = v[rj][cj];

      for(int ci = 1; ci < mb; ci++){
        int shift = 1 << ci;
        for(int cj = 0; cj < m; cj += shift << 1){
          int t = min(cj+shift, m);
          table[0][rj][ci][t-1] = v[rj][t-1];
          for(int k = t-2; k >= cj; k--){
            table[0][rj][ci][k] = f(v[rj][k], table[0][rj][ci][k+1]);
          }
          if(m <= t) break;
          table[0][rj][ci][t] = v[rj][t];
          int r = min(t+shift, m);
          for(int k = t+1; k < r; k++){
            table[0][rj][ci][k] = f(table[0][rj][ci][k-1], v[rj][k]);
          }
        }
      }
    }

    for(int ri = 1; ri < nb; ri++){
      int rshift = 1 << ri;
      for(int rj = 0; rj < n; rj += rshift << 1){
        int rt = min(rj+rshift, n);
        for(int ci = 0; ci < mb; ci++){
          for(int cj = 0; cj < m; cj++){
            table[ri][rt-1][ci][cj] = table[0][rt-1][ci][cj];
            for(int k = rt-2; k >= rj; k--){
              table[ri][k][ci][cj] = f(table[0][k][ci][cj], table[ri][k+1][ci][cj]);
            }
          }
        }
        if(n <= rt) break;
        for(int ci = 0; ci < mb; ci++){
          for(int cj = 0; cj < m; cj++){
            table[ri][rt][ci][cj] = table[0][rt][ci][cj];
            int r = min(rt+rshift, n);
            for(int k = rt+1; k < r; k++){
              table[ri][k][ci][cj] = f(table[ri][k-1][ci][cj], table[0][k][ci][cj]);
            }
          }
        }
      }
    }
    const int b = max(nb, mb);
    lookup.assign(1 << b, 0);
    for(int i = 2; i < (1 << b); i++){
      lookup.at(i) = lookup.at(i>>1) + 1;
    }
  }

  T query(int x1, int x2, int y1, int y2) const {
    x2--; y2--;
    if(x1 >= x2){
      if(y1 >= y2) return table[0][x1][0][y1];
      int p = lookup[y1 ^ y2];
      return f(table[0][x1][p][y1], table[0][x1][p][y2]);
    }
    int np = lookup[x1 ^ x2];
    if(y1 >= y2) return f(table[np][x1][0][y1], table[np][x2][0][y1]);
    int mp = lookup[y1 ^ y2];
    T R1 = f(table[np][x1][mp][y1], table[np][x1][mp][y2]);
    T R2 = f(table[np][x2][mp][y1], table[np][x2][mp][y2]);
    return f(R1, R2);
  }
};

元の要素になる二次元配列と、マージする関数を引数に入れる設計です。 ei1333さんのライブラリの影響を多分に受けています(というかパクり

使ってみる

さすがにバグが怖い

元となっている資料があるわけでもないし、何かよくわからない四次元配列を使った計算をしているんだから当然だ。

verifyができる問題を探すが、静的な二次元の区間minを要求される問題にそもそも出会ったことがない。Library-Checkerさんにもなさそう、困った。

見つけた~

さすがAOJ、何でもある

submit...

f:id:suu_0313:20210319011136p:plain
https://judge.u-aizu.ac.jp/onlinejudge/review.jsp?rid=5299736#2

MLE!?

悲しい...

けど Case#1 は通っているので多分バグってはないと思います

追記:
テストケースをダウンロードして手元環境で動かしてみました。
出力結果はあっていたので大丈夫だと思います。

僕の実装を置いておきます

おわりに

ここまでありがとうございました。

この実装のバグや、もっといい実装、いい感じにverifyができる問題など見つけた方がいらっしゃいましたらご一報いただけると嬉しいです。

追記:
思考過程書きまくってたらろくに解説を書いていないことに気づきました。
そのうち加筆します。