ei1333の日記

ぺこい

Mo's algorithm

うしのアルゴリズムです. もーもー.

(なんかアルゴリズムの概要というよりは, Moを使う問題のソースコードをたくさん貼る場所になってます…)

間違っている可能性があるので, ゆるして><

Mo’s algorithm (Query Square Root Decomposition)

概要

以下の  {3} つの条件を満たしていれば Mo を適用できる可能性がある.

  1. 配列の要素が不変.
  2. クエリを先読みできる(オフライン).
  3. 区間  {[l, r)} の結果から区間  {[l + 1, r)},  {[l - 1, r)},  {[l, r - 1)},  {[l, r + 1)} の結果を容易に計算できる.

アルゴリズムの流れとしては以下の通り. 非常にシンプル.

  1. 区間 {\sqrt N} 個のブロックに分割する.
  2. 全てのクエリを, それぞれのクエリの左端が属するブロックで昇順ソート. 左端が同じもの同士では, 右端の値で昇順ソートする.
  3. クエリごとに左端と右端を尺取りのように移動させて, 区間を伸縮させる.

事前にクエリを適切に並び替えておき, あとは区間の移動を尺取りして伸縮させることで,  {1} 回あたりの伸縮の計算量を  {O(\alpha)} とすると計算量が  {O(\alpha (N + Q) \sqrt N)} になるアルゴリズムといえる.

計算量

計算量の解析は比較的容易. 左端と右端, それぞれポインタの移動回数を見ればよい(微妙な定数倍は無視しています).

  • 左端:  {1} つのクエリあたり最大  {\sqrt N} 回移動する.  {Q} 個のクエリがあるので全体で  {O(Q \sqrt N)} 回. f:id:ei1333:20170908145752p:plain
  • 右端: それぞれのブロックごとに最大  {N} 回移動する.  {\sqrt N} 個のブロックがあるので全体で  {O(N \sqrt N)} 回. f:id:ei1333:20170908144703p:plain

したがって, 伸縮  {1} 回あたりの計算量を  {O(\alpha)} とすると, 全体の計算量は  {O(\alpha (N + Q) \sqrt N)} となることがわかる.  {O(\alpha)} {O(1)} だったり  {O(\log N)} だったりする.

実装例

struct Mo
{
  vector< int > left, right, order;
  vector< bool > v;
  int width;
  int nl, nr, ptr;

  Mo(int n) : width((int) sqrt(n)), nl(0), nr(0), ptr(0), v(n) {}

  void insert(int l, int r) /* [l, r) */
  {
    left.push_back(l);
    right.push_back(r);
  }

  /* ソート */
  void build()
  {
    order.resize(left.size());
    iota(begin(order), end(order), 0);
    sort(begin(order), end(order), [&](int a, int b)
    {
      if(left[a] / width != left[b] / width) return left[a] < left[b];
      return right[a] < right[b];
    });
  }

  /* クエリを 1 つぶんすすめて, クエリのidを返す */
  int process()
  {
    if(ptr == order.size()) return (-1);
    const auto id = order[ptr];
    while(nl > left[id]) distribute(--nl);
    while(nr < right[id]) distribute(nr++);
    while(nl < left[id]) distribute(nl++);
    while(nr > right[id]) distribute(--nr);
    return (order[ptr++]);
  }

  inline void distribute(int idx)
  {
    v[idx].flip();
    if(v[idx]) add(idx);
    else del(idx);
  }

  void add(int idx);

  void del(int idx);
};

実装例についての補足

それぞれのクエリの区間  {[l, r)} を順に指定して insert(l, r) を呼び出す. そのあと build() を呼び出すと, 追加されたクエリを  {\sqrt{N}}区間に分けてソートする.

process() を呼び出すことで, いい感じに区間が進められて  {1} つぶんのクエリを処理できる(戻り値としてもともとのクエリの要素番号を返す). 問題に応じて実装すべきなのは add() と del() で, それぞれ  {\mathrm{idx}} 番目の要素が加わる時,  {\mathrm{idx}} 番目の要素が消える時に呼び出される.


Mo を適用可能な例題と解法をあげる.

DQUERY - D-query

www.spoj.com

長さ  {n(1 \le n \le 3 \times 10^4)} の数列  {a_i(1 \le a_i \le 10^6)} が与えられます. 区間  {[l_i, r_i] (1 \le l_i \le r_i \le n)} の値の種類数を求める  {q(1 \le q \le 2 \times 10^5)} 個のクエリに答えてください.

区間の伸縮の操作がとても容易.

現在の状態として以下の  {2} つを持っておく.

  1.  {\mathrm{cnt}[i] := } 区間内の値  {i} の個数
  2.  {\mathrm{sum} := } 区間内の種類数

要素  {\mathrm{idx}} を加える時  {\mathrm{cnt}[a_i]} {1} を足す. このときもともとの  {a_i} の出現数が  {0} であれば  {\mathrm{sum}} {1} 加える. 要素  {\mathrm{idx}} を除く時  {\mathrm{cnt}[a_i]} から  {1} を引く. 同様に  {a_i} の出現数が  {0} になるとき  {\mathrm{sum}} {1} 減らす.

区間の伸縮の操作が  {O(1)} でできたため, 全体の計算量は  {O((N + Q) \sqrt N)} となる.

(以下 Mo の部分は同じなので省略します)

int N, A[300000], Q;
int ans[200000];
int cnt[1000001], sum;

void Mo::add(int idx)
{
  if(cnt[A[idx]]++ == 0) ++sum;
}

void Mo::del(int idx)
{
  if(--cnt[A[idx]] == 0) --sum;
}

int main()
{
  scanf("%d", &N);
  for(int i = 0; i < N; i++) {
    scanf("%d", &A[i]);
  }
  scanf("%d", &Q);
  Mo mo(N);
  for(int i = 0; i < Q; i++) {
    int a, b;
    scanf("%d %d", &a, &b);
    mo.insert(--a, b);
  }
  mo.build();
  for(int i = 0; i < Q; i++) {
    ans[mo.process()] = sum;
  }
  for(int i = 0; i < Q; i++) {
    printf("%d\n", ans[i]);
  }
}

Mo mo っておもしろくないですか? おもしろくないですね.

FREQUENT - Frequent values

www.spoj.com

長さ  {n(1 \le n \le 3 \times 10^5)} の数列  {a_i(-10^5 \le a_i \le 10^5)} が与えられます. 区間  {[l_i, r_i] (1 \le l_i \le r_i \le n)} の最頻値の出現回数を求める  {q(1 \le q \le 2 \times 10^5)} 個のクエリに答えてください.

これも  {O(1)} でできる. 現在の状態として以下の  {3} つを持っておく.

  1.  {\mathrm{sum}[i] := } 区間内の値  {i} の個数
  2.  {\mathrm{reach}[i] := } 区間内に  {i} 個ある値の個数
  3.  {\mathrm{now} := } 最頻値の出現回数

要素  {\mathrm{idx}} を加える時  {\mathrm{sum}[a_i]} {1} を足す. このあと  {\mathrm{reach}[\mathrm{sum}[a_i]]} {1} かつ  {\mathrm{sum}[a_i]} {\mathrm{now}} が大きければ, 最頻値の出現回数が更新されたことを意味している. 削除も同じようにいい感じに処理をする.

int N, M, A[100000], ans[100000];
int reach[100001], sum[200001], now;

void Mo::add(int idx)
{
  reach[sum[A[idx]]]--;
  ++sum[A[idx]];
  reach[sum[A[idx]]]++;
  if(sum[A[idx]] > now) ++now;
}

void Mo::del(int idx)
{
  reach[sum[A[idx]]]--;
  if(sum[A[idx]] == now && reach[sum[A[idx]]] == 0) --now;
  --sum[A[idx]];
  reach[sum[A[idx]]]++;
}

int main()
{
  while(scanf("%d", &N), N) {
    memset(sum, 0, sizeof(sum));
    memset(reach, 0, sizeof(reach));
    now = 0;

    scanf("%d", &M);
    for(int i = 0; i < N; i++) {
      scanf("%d", &A[i]);
      A[i] += 100000;
    }
    Mo mo(N);
    for(int i = 0; i < M; i++) {
      int x, y;
      scanf("%d %d", &x, &y);
      mo.insert(--x, y);
    }
    mo.build();
    for(int i = 0; i < M; i++) {
      ans[mo.process()] = now;
    }
    for(int i = 0; i < M; i++) {
      printf("%d\n", ans[i]);
    }
  }
}

木上の Mo

DFSして木の頂点を適当な順序で並べておくことで, 木に対する様々なオフラインクエリに答えることができる.

部分木に対するクエリ

部分木は頂点をDFSの行きがけ順(pre-order) に並べておくと, 以下のように配列の区間に落とすことが可能

f:id:ei1333:20170911201752p:plain

具体的な実装としては, 現在何個の頂点を訪れたかという変数を更新しながら, それぞれの頂点に訪れたとき, また戻るタイミングでその変数の情報を保存していけば良い(上の図は閉区間でしたが, 以下では半開区間の実装になっています).

vector< int > tree[100000];
int order[100000], in[100000], out[100000], ptr; // [in, out)

void dfs(int idx, int par = -1)
{
  order[ptr] = idx;
  in[idx] = ptr++;
  for(auto &to : tree[idx]) if(to != par) dfs(to, idx);
  out[idx] = ptr;
}

D. Tree and Queries

codeforces.com

 {n} 頂点の木があって 頂点  {i} の色は  {c_i} です. 頂点  {v_j} の部分木の頂点をみたときに,  {k_j} 個以上ある色の個数を求める  {m} 個のクエリに答えてください.

現在の状態として以下の  {2} つを持っておく.

  •  {\mathrm{color}[i] := } 区間内で色  {i} が現れた個数
  •  {\mathrm{sum}[i] := } 区間内で  {i} 個以上ある色の個数

 {\mathrm{color}[i]} が 1 つ増減した時, 変化するのは  {\mathrm{sum}[\mathrm{color}[i]]} だけなので更新が  {O(1)} で行える.

vector< int > tree[100000];
int rev[100000], in[100000], out[100000], ptr; // [in, out)
int N, M, C[100000], V[100000], K[100000];
int color[100000], sum[100001];
int ans[100000];

void Mo::add(int idx)
{
  ++sum[++color[C[rev[idx]]]];
}

void Mo::del(int idx)
{
  --sum[color[C[rev[idx]]]--];
}

void dfs(int idx, int par = -1)
{
  rev[ptr] = idx;
  in[idx] = ptr++;
  for(auto &to : tree[idx]) if(to != par) dfs(to, idx);
  out[idx] = ptr;
}

int main()
{
  scanf("%d %d", &N, &M);
  for(int i = 0; i < N; i++) {
    scanf("%d", &C[i]);
    --C[i];
  }
  for(int i = 1; i < N; i++) {
    int a, b;
    scanf("%d %d", &a, &b);
    --a, --b;
    tree[a].emplace_back(b);
    tree[b].emplace_back(a);
  }
  dfs(0);

  Mo mo(N);
  for(int i = 0; i < M; i++) {
    scanf("%d %d", &V[i], &K[i]);
    --V[i];
    mo.insert(in[V[i]], out[V[i]]);
  }
  mo.build();

  for(int i = 0; i < M; i++) {
    int idx = mo.process();
    ans[idx] = sum[K[idx]];
  }

  for(int i = 0; i < M; i++) {
    printf("%d\n", ans[i]);
  }
}

パスに対するクエリ

パスに対するクエリに対しても適用可能である.

まずオイラーツアーをして, 訪問順に頂点を並べる. たとえば以下の図の通りになる. 行った時と戻った時両方について見ることに注意する. またそれぞれの辺の端点に対して  {1} つずつ消費するので  {2n-1} 個の大きさの配列が必要.

f:id:ei1333:20170907200655p:plain

ある  {2} 頂点  {u, v} を取ってきて, 区間  {[in_u, in_v]} に着目してみる.

f:id:ei1333:20170907202031p:plain

この区間に着目すると, パス上にある辺は奇数回, パス上にない辺は偶数回出現するようになっている. 奇数回目の訪問時にデータを追加し, 偶数回目の訪問時にデータを削除する操作を行うことで, パス上にない頂点を相殺できることを意味している. この操作をそのまま Mo に適用すればそのままできる.

辺に対するクエリははいだったが, 頂点に対するクエリに対してもこの応用で実現可能. 任意の頂点間の LCA(最小共通祖先) を  {O(\log n)} くらいで求められると便利なので準備しておく(ここではダブリングで求めることとする).

それぞれの辺のコストを深さが深い頂点のコストと対応付けるようにすると, 以下のような感じになって上手くいくことがわかる. LCA の頂点だけ例外処理みたいなことをする.

f:id:ei1333:20170907202411p:plain

f:id:ei1333:20170907202424p:plain

(Mo ではないけど, これと同様のテクを使って Tree | Aizu Online Judge とか highway: 高速道路 (Highway) - 2010年 日本情報オリンピック春合宿OJ | AtCoder とか解ける.(話が脱線してる気がしてきた, なんか木に対する有用テク流出書みたいになってませんか?なってませんね))

以下のコードでは, 普通の Mo を拡張してパス(頂点)に対するクエリにマッチするように書き直している.

(実装が険しいね, えーん, LCAが間違っていて1時間溶かした)

struct Mo_Tree_Vertex
{
  vector< vector< int > > g, parent;
  vector< int > vs, in;

  vector< int > left, right, order, lca, dep;
  vector< bool > v;
  int width;
  int nl, nr, ptr;

  Mo_Tree_Vertex(int n) : width((int) sqrt(2 * n - 1)), nl(0), nr(0), ptr(0), g(n), in(n), dep(n), v(n)
  {
    const auto lg = (int) (log2(n) + 1);
    parent.resize(lg, vector< int >(n));
    vs.reserve(2 * n - 1);
  }

  void add_edge(int x, int y)
  {
    g[x].push_back(y);
    g[y].push_back(x);
  }

  void dfs(int idx, int depth, int par)
  {
    in[idx] = (int) vs.size();
    dep[idx] = depth;
    parent[0][idx] = par;
    vs.push_back(idx);
    for(auto &to : g[idx]) {
      if(to == par) continue;
      dfs(to, depth + 1, idx);
      vs.push_back(to); // 本来のEuler-tourはidxだが, 今回は深い方の頂点をとりたいのでto
    }
  }

  int getlca(int u, int v)
  {
    if(dep[u] > dep[v]) swap(u, v);
    for(int k = 0; k < parent.size(); k++) {
      if(((dep[v] - dep[u]) >> k) & 1) {
        v = parent[k][v];
      }
    }
    if(u == v) return (u);
    for(int k = (int) parent.size() - 1; k >= 0; k--) {
      if(parent[k][u] != parent[k][v]) {
        u = parent[k][u];
        v = parent[k][v];
      }
    }
    return (parent[0][u]);
  }

  void insert(int x, int y)
  {
    if(in[x] > in[y]) swap(x, y);
    left.push_back(in[x] + 1);
    right.push_back(in[y] + 1);
    lca.push_back(getlca(x, y));
  }

  void prepare()
  {
    dfs(0, 0, -1);
    for(int k = 0; k + 1 < parent.size(); k++) {
      for(int i = 0; i < parent[k].size(); i++) {
        if(parent[k][i] == -1) parent[k + 1][i] = -1;
        else parent[k + 1][i] = parent[k][parent[k][i]];
      }
    }
  }

  void build()
  {
    order.resize(left.size());
    iota(begin(order), end(order), 0);
    sort(begin(order), end(order), [&](int a, int b)
    {
      if(left[a] / width != left[b] / width) return left[a] < left[b];
      return right[a] < right[b];
    });
  }

  int process()
  {
    if(ptr == order.size()) return (-1);
    if(ptr > 0) distribute(lca[order[ptr - 1]]); // 前のクエリで追加したLCAを削除
    const auto id = order[ptr];
    while(nl > left[id]) distribute(vs[--nl]);
    while(nr < right[id]) distribute(vs[nr++]);
    while(nl < left[id]) distribute(vs[nl++]);
    while(nr > right[id]) distribute(vs[--nr]);
    distribute(lca[id]); // LCA の頂点を例外として追加する
    return (order[ptr++]);
  }

  /* 偶数回目の訪問時は削除, 奇数回目の訪問時は追加を呼び出す */
  inline void distribute(int vertex)
  {
    v[vertex].flip();
    if(v[vertex]) add(vertex);
    else del(vertex);
  }

  void add(int idx);

  void del(int idx);
};

COT2 - Count on a tree II

www.spoj.com

 {n} 頂点の木があって 頂点  {i} には値  {v_i} が書かれています. 頂点  {a_i, b_i} のパス上の頂点に書かれた値の種類数を求める  {q} 個のクエリに答えてください.

配列のときの種類数の求め方と同様.

int N, M, X[40000];
int ans[100000];
int cnt[40000], sum;

void Mo_Tree_Vertex::add(int idx)
{
  if(cnt[X[idx]]++ == 0) ++sum;
}

void Mo_Tree_Vertex::del(int idx)
{
  if(--cnt[X[idx]] == 0) --sum;
}

int main()
{
  scanf("%d %d", &N, &M);

  vector< int > vs(N);
  for(int i = 0; i < N; i++) {
    scanf("%d", &X[i]);
    vs[i] = X[i];
  }
  sort(begin(vs), end(vs));
  vs.erase(unique(begin(vs), end(vs)), end(vs));
  for(int i = 0; i < N; i++) {
    X[i] = lower_bound(begin(vs), end(vs), X[i]) - begin(vs);
  }

  Mo_Tree_Vertex mo(N);
  for(int i = 1; i < N; i++) {
    int a, b;
    scanf("%d %d", &a, &b);
    mo.add_edge(--a, --b);
  }
  mo.prepare();

  for(int i = 0; i < M; i++) {
    int a, b;
    scanf("%d %d", &a, &b);
    mo.insert(--a, --b);
  }
  mo.build();
  for(int i = 0; i < M; i++) {
    ans[mo.process()] = sum;
  }
  for(int i = 0; i < M; i++) {
    printf("%d\n", ans[i]);
  }
}

Bubble Cup X - Finals I. Dating

codeforces.com

 {n} 頂点の木があります. 頂点  {i} には男性か女性が住んでいて, その人は値  {v_i} が好きです. 頂点  {a_i, b_i} のパス上にある頂点のみを取り出して, 同じ値が好きな男女でペアを組めるとき, 考えられるペアの個数を求める  {q} 個のクエリに答えてください.

値の種類ごとに独立に数えることができることに着目すると, 種類数の問題とほとんど同じ問題に帰着できる.

ある頂点が加わるときにはもう一方の同じ値を持つ頂点数を答えに加えればよくて, 逆も然り.

typedef long long int64;

int N, Q, A[100000], B[100000];
int cnt1[100000], cnt2[100000];
int64 sum, ans[100000];

void Mo_Tree_Vertex::add(int idx)
{
  if(A[idx] == 0) {
    cnt1[B[idx]]++;
    sum += cnt2[B[idx]];
  } else {
    cnt2[B[idx]]++;
    sum += cnt1[B[idx]];
  }
}

void Mo_Tree_Vertex::del(int idx)
{
  if(A[idx] == 0) {
    cnt1[B[idx]]--;
    sum -= cnt2[B[idx]];
  } else {
    cnt2[B[idx]]--;
    sum -= cnt1[B[idx]];
  }
}

int main()
{
  scanf("%d", &N);

  for(int i = 0; i < N; i++) scanf("%d", &A[i]);
  vector< int > vs(N);
  for(int i = 0; i < N; i++) scanf("%d", &B[i]), vs[i] = B[i];

  sort(begin(vs), end(vs));
  vs.erase(unique(begin(vs), end(vs)), end(vs));
  for(int i = 0; i < N; i++) {
    B[i] = lower_bound(begin(vs), end(vs), B[i]) - begin(vs);
  }

  Mo_Tree_Vertex mo(N);
  for(int i = 1; i < N; i++) {
    int a, b;
    scanf("%d %d", &a, &b);
    mo.add_edge(--a, --b);
  }
  mo.prepare();

  scanf("%d", &Q);
  for(int i = 0; i < Q; i++) {
    int a, b;
    scanf("%d %d", &a, &b);
    mo.insert(--a, --b);
  }
  mo.build();
  for(int i = 0; i < Q; i++) {
    ans[mo.process()] = sum;
  }
  for(int i = 0; i < Q; i++) {
    printf("%lld\n", ans[i]);
  }
}

The L-th Number

The L-th Number | Aizu Online Judge

 {n} 頂点の木があって 頂点  {i} には値  {v_i} が書かれています. 頂点  {a_i, b_i} のパス上の頂点に書かれた値をソートしたときの  {k_i} 番目の値を求める  {q} 個のクエリに答えてください.

座標圧縮したあとBITを使う. 伸縮は  {O(\log n)} で,  {k} 番目の値を求めるのはBIT上の二分探索を使えば  {O(\log n)}. 全体で  {O(q \sqrt n \log n)}. 計算量がヤバな気がするけど BIT は定数だと信じればlogは消えるので神頼みをすれば解ける. (僕にはタプリスさん(@ei13333) という神ならぬ天使がついています(は? 2.57secらしい

補足すると, BIT で  {i} 番目の要素に  {x} 足す操作は,  {[i, \infty)} {x} 足す操作とみなせる. BIT の要素は単調増加になっており,  {[1, x]} の和が  {k} を超えるような最小の要素位置  {x} を二分探索できる.

typedef long long int64;

template< class T >
struct BinaryIndexedTree
{
  vector< T > data;

  BinaryIndexedTree(int sz)
  {
    data.assign(++sz, 0);
  }

  T sum(int k)
  {
    T ret = 0;
    for(++k; k > 0; k -= k & -k) ret += data[k];
    return (ret);
  }

  void add(int k, T x)
  {
    for(++k; k < data.size(); k += k & -k) data[k] += x;
  }
};

template< class T >
struct AdditionalBIT : BinaryIndexedTree< T >
{
  int curr;

  AdditionalBIT(int sz) : BinaryIndexedTree< T >(sz)
  {
    curr = 1;
    while(curr <= sz) curr <<= 1;
  }

  int lower_bound(T w)
  {
    if(w <= 0) return (0);
    int i = 0;
    for(int k = curr; k > 0; k >>= 1) {
      if(i + k < BinaryIndexedTree< T >::data.size() && BinaryIndexedTree< T >::data[i + k] < w) {
        w -= BinaryIndexedTree< T >::data[i + k];
        i += k;
      }
    }
    return (i);
  }
};

int N, Q, X[100000];
int A[100000], B[100000], K[100000];
int64 ans[100000];
AdditionalBIT< int > bit(100000);

void Mo_Tree_Vertex::add(int idx)
{
  bit.add(X[idx], 1);
}

void Mo_Tree_Vertex::del(int idx)
{
  bit.add(X[idx], -1);
}

int main()
{
  scanf("%d %d", &N, &Q);
  vector< int > vs(N);
  for(int i = 0; i < N; i++) scanf("%d", &X[i]), vs[i] = X[i];
  sort(begin(vs), end(vs));
  vs.erase(unique(begin(vs), end(vs)), end(vs));
  for(int i = 0; i < N; i++) {
    X[i] = lower_bound(begin(vs), end(vs), X[i]) - begin(vs);
  }
  Mo_Tree_Vertex mo(N);
  for(int i = 1; i < N; i++) {
    int a, b;
    scanf("%d %d", &a, &b);
    mo.add_edge(--a, --b);
  }
  mo.prepare();
  for(int i = 0; i < Q; i++) {
    scanf("%d %d %d", &A[i], &B[i], &K[i]);
    mo.insert(--A[i], --B[i]);
  }
  mo.build();
  for(int i = 0; i < Q; i++) {
    int idx = mo.process();
    ans[idx] = vs[bit.lower_bound(K[idx])];
  }
  for(int i = 0; i < Q; i++) {
    printf("%lld\n", ans[i]);
  }
}

D - 閉路

abc014.contest.atcoder.jp

明らかに無駄だが, Mo を貼っても出来る. 練習用にはいいかも.

int sum;
void Mo_Tree_Vertex::add(int idx) { sum++; }
void Mo_Tree_Vertex::del(int idx) { sum--; }

int main()
{
  int N, Q, ans[100000];
  scanf("%d", &N);
  Mo_Tree_Vertex mo(N);
  for(int i = 1; i < N; i++) {
    int x, y;
    scanf("%d %d", &x, &y);
    mo.add_edge(--x, --y);
  }
  mo.prepare();
  scanf("%d", &Q);
  for(int i = 0; i < Q; i++) {
    int a, b;
    scanf("%d %d", &a, &b);
    mo.insert(--a, --b);
  }
  mo.build();
  for(int i = 0; i < Q; i++) {
    ans[mo.process()] = sum;
  }
  for(int i = 0; i < Q; i++) {
    printf("%d\n", ans[i]);
  }
}

おまけ

その 1: 定数倍改善

もともとのソートでは, 区間の左端が属するブロックが同じならば右端の昇順になるようにソートしていた. ここで偶数番目のブロックを反転, つまり奇数番目のブロックは昇順, 偶数番目のブロックは降順となるようにソートすると, ブロックの切り替わり時に削除のために一番左に戻ることがなくなって, 効率が上がることが期待される.

(2.57 sec -> 1.52 sec, すごくはやいね(実行速度の比較は L番目の数字 を用いています(比較しやすかったので)))

sort(begin(order), end(order), [&](int a, int b)
{
 if(left[a] / width != left[b] / width) return left[a] < left[b];
 return bool((right[a] < right[b]) ^ (left[a] / width % 2));
});

その 2: 時空間 Mo

Twitterでふと目にしたのでちょっとだけ(ちょっとだけでわない). (参考: Update query on Mo's Algorithm - Codeforces など)

Mo’s を応用して, 配列の要素の更新もやってしまおうという考え方. アルゴリズムの概略と計算量は以下に示す. 計算量がアなので配列の長さ  {n} とクエリの数  {q = n} で同一視している.

  1. 配列を  {n^{\frac 1 3}} 個のブロックに分割する. 当然  {1} ブロックあたりの長さは  {n^{\frac 2 3}} となる.
  2. 普通の Mo と同じように, 左端のブロックで昇順ソート, 同じならば右端のブロックで降順ソートする(右端のブロックというところが異なる). 右端のブロックも同じならばクエリが与えられた時刻でソートする.
  3. 現在の状態を表す変数として,  {nl, nr, time} を持つ.  {nl} は左端,  {nr} は右端,  {time} は時刻を表す.
  4.  {nl} が動く回数を考える.  {1} ブロックの長さが  {n^{\frac 2 3}} だったので最大で  {n n^{\frac 2 3} = n^{\frac 5 3} } 回移動する. f:id:ei1333:20170910215955p:plain
  5.  {nr} が動く回数を考える.  {1} ブロックの長さが  {n^{\frac 2 3}} だったので最大で  {n n^{\frac 2 3} = n^{\frac 5 3} } 回移動する. またこれとは別に  {n^{\frac 1 3}} 通りの  {nl} (青色の塗りつぶし)に対して, 長さ  {n} の移動を行うので  {n^{\frac 4 3} } 回移動する. f:id:ei1333:20170910220012p:plain
  6.  {time} が動く回数を考える.  {nl, nr} の組み合わせは  {n^{\frac 1 3} n^{\frac 1 3} = n^{\frac 2 3}} 通りある. それぞれに対して  {n} 回移動するので  {n^{\frac 5 3} }. (なんか誤解を招きそうな図になってますが, 上の図の1行分の操作が下のアです) f:id:ei1333:20170910220022p:plain
  7. すべてのポインタの動く回数が  { n^{\frac 5 3} } で抑えられたため,  {1} 回あたりの伸縮の計算量を  {\alpha} とすると全体で  {O(\alpha n^{\frac 5 3})} となる.

更新ができるといっても, ヤバい更新とかはヤバくなるのでア(シンプルな  {1} 点更新みたいなのだと出来る可能性が高い, 定数時間で更新するための前計算とかが大事な気がする, 前に戻る時の操作に多分前計算が必要).

pekempey さんによる区間更新の例 (なるほど…) pekempey.hatenablog.com

struct Mo_3D
{
  vector< int > left, right, index, order;
  vector< bool > v;
  int width;
  int nl, nr, time, ptr;

  Mo_3D(int n) : width((int) pow(n, 2.0 / 3.0)), nl(0), nr(0), ptr(0), time(0), v(n) {}

  void insert(int idx, int l, int r) /* [l, r) */
  {
    left.push_back(l);
    right.push_back(r);
    index.push_back(idx);
  }

  void build()
  {
    order.resize(left.size());
    iota(begin(order), end(order), 0);
    sort(begin(order), end(order), [&](int a, int b)
    {
      if(left[a] / width != left[b] / width) return left[a] < left[b];
      if(right[a] / width != right[b] / width) return bool((right[a] < right[b]) ^ (left[a] / width % 2));
      return bool((index[a] < index[b]) ^ (right[a] / width % 2));
    });
  }

  int process()
  {
    if(ptr == order.size()) return (-1);
    const auto id = order[ptr];
    while(time < index[id]) addquery(time++);
    while(time > index[id]) delquery(--time);
    while(nl > left[id]) distribute(--nl);
    while(nr < right[id]) distribute(nr++);
    while(nl < left[id]) distribute(nl++);
    while(nr > right[id]) distribute(--nr);
    return (index[order[ptr++]]);
  }

  inline void distribute(int idx)
  {
    v[idx].flip();
    if(v[idx]) add(idx);
    else del(idx);
  }

  void add(int idx);

  void del(int idx);

  void addquery(int idx);

  void delquery(int idx);
};

クエリを  {1} 個進めるための addquery() と  {1} 個戻すための delquery() の実装が追加で必要になる.

addquery() と delquery() の処理で, 変更する要素が現在のクエリの範囲  {[nl, nr)} に入っているかどうかで虚無な場合分けが発生しそうだが, 楽にする手段として  {v[idx]} が true のときに 1 回 distribute(idx) を呼び出して範囲から切り離した後に要素を更新して, もう一度 distribute(idx) を呼び出すような実装が考えられる(要素の一部を切り離せないアのときは仕方ないね).

これから問題解くの非常に虚無なんですが, それが求められている気がするので頑張ります(?).

Range Sum Query

Range Sum Query (RSQ) | Aizu Online Judge

区間の伸縮は要素の値をえいするだけで, 変更要素もまあ適当にはいなのでともに  {O(1)}. 全体で  {O(n^{\frac 5 3})}.

1sec なので非常に厳しい気持ちになるが, 信じればぎりぎり間に合う.

int N, Q, A[100000], B[100000], C[100000];
int D[100000];
int ans[100000], sum;

void Mo_3D::add(int idx)
{
  sum += D[idx];
}

void Mo_3D::del(int idx)
{
  sum -= D[idx];
}

void Mo_3D::addquery(int idx)
{
  if(A[idx] == 1) return; // タイプ1のクエリは無視(それはそう
  if(v[B[idx]]) { // 区間に入っている場合は取り消してから処理
    distribute(B[idx]);
    D[B[idx]] += C[idx];
    distribute(B[idx]);
  } else {
    D[B[idx]] += C[idx];
  }
}

void Mo_3D::delquery(int idx)
{
  if(A[idx] == 1) return;
  if(v[B[idx]]) {
    distribute(B[idx]);
    D[B[idx]] -= C[idx];
    distribute(B[idx]);
  } else {
    D[B[idx]] -= C[idx];
  }
}

int main()
{
  scanf("%d %d", &N, &Q);

  Mo_3D mo(N);
  for(int i = 0; i < Q; i++) {
    scanf("%d %d %d", &A[i], &B[i], &C[i]);
    --B[i];
    if(A[i] == 1) mo.insert(i, B[i], C[i]);
  }

  mo.build();

  int idx;
  while((idx = mo.process()) != -1) {
    ans[idx] = sum;
  }

  for(int i = 0; i < Q; i++) {
    if(A[i] == 1) printf("%d\n", ans[i]);
  }
}

C. Goodbye Souvenir

Problem - C - Codeforces

長さ  {n(1 \le n \le 10^5)} の数列  {a_i(1 \le a_i \le n)} が与えられます. 以下の  {m(1 \le m \le 10^5)} 個のクエリを順に処理してください.
① 要素  {p} {x} に更新する
区間  {[l_i, r_i] (1 \le l_i \le r_i \le n)} を取り出したとき, 区間 {1} 回以上出現するすべての値  {x} について (区間内の最後の出現位置 - 区間内の最初の出現位置) を求めそれを足し合わせたものを求める

現在の状態として以下の  {3} つを持っておく.

  •  {\mathrm{first}[i] := } 区間内で値  {i} が最初に現れた位置 (出現してなければ -1)
  •  {\mathrm{last}[i] := } 区間内で値  {i} が最後に現れた位置 (出現してなければ -1)
  •  {\mathrm{ans} :=} 答え(つまり配列 last の和 - first の和)

要素を加えるときは必要ならば first[i], last[i], ans を更新する. 要素をけずるとき, first[i], last[i] が削る idx を指していれば, それを適切なもの(1個前の出現か1個後の出現か出現しなくなるかのいずれか)に更新する. 出現の情報は前計算しておくことで  {O(1)} で求めることが出来る.

また要素を更新するときは, 直前直後の同じ値の要素と自分(の1個前の出現位置, 1個後の出現位置の情報)なので高々  {5} 箇所, これらを前計算しておくと  {O(1)} で求めることができる.

ココらへんの問題になると Mo の恩恵を受けられているような気がする.

typedef long long int64;

int N, M, A[100005], buff[100005];

set< int > line[100005];
int init_pv[100005], init_nt[100005];
int X[100005], Y[100005], Z[100005];

vector< pair< int, int > > update_pv[100000], update_nt[100000];
int first[100001], last[100001];

int64 ret[100001], ans;

void Mo_3D::add(int idx)
{
  if(first[A[idx]] == -1 || (first[A[idx]] != -1 && idx < first[A[idx]])) {
    if(first[A[idx]] != -1) ans += first[A[idx]];
    ans -= idx;
    first[A[idx]] = idx;
  }
  if(last[A[idx]] == -1 || (last[A[idx]] != -1 && idx > last[A[idx]])) {
    if(last[A[idx]] != -1) ans -= last[A[idx]];
    ans += idx;
    last[A[idx]] = idx;
  }
}

void Mo_3D::del(int idx)
{
  if(first[A[idx]] == idx) {
    ans += idx;
    if(init_nt[idx] != -1 && init_nt[idx] < nr) {
      first[A[idx]] = init_nt[idx];
      ans -= first[A[idx]];
    } else {
      first[A[idx]] = -1;
    }
  }
  if(last[A[idx]] == idx) {
    ans -= idx;
    if(init_pv[idx] != -1 && nl <= init_pv[idx]) {
      last[A[idx]] = init_pv[idx];
      ans += last[A[idx]];
    } else {
      last[A[idx]] = -1;
    }
  }
}

void Mo_3D::addquery(int idx)
{
  if(X[idx] == 2) return;

  if(A[Y[idx]] == Z[idx]) return;

  if(v[Y[idx]]) {
    distribute(Y[idx]);
    for(auto &q : update_pv[idx]) swap(init_pv[q.first], q.second);
    for(auto &q : update_nt[idx]) swap(init_nt[q.first], q.second);
    swap(A[Y[idx]], Z[idx]);
    distribute(Y[idx]);
  } else {
    for(auto &q : update_pv[idx]) swap(init_pv[q.first], q.second);
    for(auto &q : update_nt[idx]) swap(init_nt[q.first], q.second);
    swap(A[Y[idx]], Z[idx]);
  }
}

void Mo_3D::delquery(int idx)
{
  addquery(idx);
}

int main()
{
  scanf("%d %d", &N, &M);
  for(int i = 0; i < N; i++) {
    scanf("%d", &A[i]);
    line[A[i]].emplace(i);
    buff[i] = A[i];
  }
  for(int i = 0; i < M; i++) {
    scanf("%d %d %d", &X[i], &Y[i], &Z[i]);
    --Y[i];
  }

  memset(init_pv, -1, sizeof(init_pv));
  memset(init_nt, -1, sizeof(init_nt));
  for(int i = 0; i < N; i++) {
    auto p = line[A[i]].lower_bound(i);
    if(p != begin(line[A[i]])) init_pv[i] = *prev(p);
    if(next(p) != end(line[A[i]])) init_nt[i] = *next(p);
  }

  Mo_3D mo(N);
  for(int i = 0; i < M; i++) {

    if(X[i] == 1) {
      {
        auto p = line[A[Y[i]]].lower_bound(Y[i]);
        const int pv = (p != begin(line[A[Y[i]]]) ? *prev(p) : -1);
        const int nt = (next(p) != end(line[A[Y[i]]]) ? *next(p) : -1);
        if(~pv) update_nt[i].emplace_back(pv, nt);
        if(~nt) update_pv[i].emplace_back(nt, pv);
        line[A[Y[i]]].erase(Y[i]);
      }
      {
        A[Y[i]] = Z[i];
        line[A[Y[i]]].emplace(Y[i]);
        auto p = line[A[Y[i]]].lower_bound(Y[i]);
        const int pv = (p != begin(line[A[Y[i]]]) ? *prev(p) : -1);
        const int nt = (next(p) != end(line[A[Y[i]]]) ? *next(p) : -1);
        if(~pv) update_nt[i].emplace_back(pv, Y[i]);
        if(~nt) update_pv[i].emplace_back(nt, Y[i]);
        update_pv[i].emplace_back(Y[i], pv);
        update_nt[i].emplace_back(Y[i], nt);
      }
    } else {
      mo.insert(i, Y[i], Z[i]);
    }
  }

  mo.build();

  int id;
  memset(first, -1, sizeof(first));
  memset(last, -1, sizeof(last));
  for(int i = 0; i < N; i++) A[i] = buff[i];
  while((id = mo.process()) != -1) ret[id] = ans;

  for(int i = 0; i < M; i++) {
    if(X[i] == 2) printf("%lld\n", ret[i]);
  }
}

問題探すのが面倒なのでこれで終わり(これは見せかけで, 定数倍も厳しい上に実装も厳しいので体力のほうがね)

参考資料

www.hackerearth.com

blog.anudeep2011.com

最後に

もう。