セグメント木をあきらめた人のための平方分割

この記事はCompetitive Programming Advent Calendar 2016(その2)の12月15日の記事です。

www.adventar.org

はじめに

私はRMQのような典型的なセグメント木までは比較的容易に実装できるのですが,それよりも難しいクエリに対応したセグメント木を書くのは難しく感じています。とはいえ「区間に対するクエリを処理するための道具」としてセグメント木を書けることが前提であるような場合もあるため放っておくわけにはいきません。そのため,私はセグメント木を書くことは諦めて平方分割を使っています。

平方分割を使うことには以下のようなメリットがあると考えています。

  • 再帰的に考える必要がない

セグメント木はlog N層の構造であり再帰的に考える必要があるのに対し,平方分割は2層しかないので考えやすいです。

セグメント木のような二分木をデバッグするのはなかなか大変です。平方分割では配列を2つprintするだけで状態が分かるので,デバッグも比較的容易です。

  • セグメント木と比べて意外と遅くない

セグメント木と平方分割を比べると,計算量的にはO(log N)からO(√N)になり悪化しています。しかしながら,平方分割ではメモリへのアクセスが連続的であるためキャッシュが効きやすく,上手に実装したセグメント木と比べても3〜5倍遅くなる程度で収まることが多いです。この差が致命的になる場合もあるのかもしれませんが,セグメント木は実装により速度に幅が出やすいので,時間制限に余裕を持たせている問題が多いように感じています。

この記事では,簡単なクエリに対応する平方分割の説明からはじめ,少し応用的なクエリが要求されたときに私がどのようなことに気をつけているかについて説明したいと思います。

基本事項

「平方分割って何?」という方は以下のスライドのp.19〜をご参照ください。

www.slideshare.net

1点に対する変更クエリ・区間に対する質問クエリ

Range Sum Query

http://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=DSL_2_B

この問題では2種類のクエリに対応する必要があります。

  • add(i,x): a[i]にxを加算する。
  • getSum(s,t): a[s],a[s+1],…,a[t]の合計値を出力する。

ここから平方分割を書くにあたって,「どのような補助データが必要か」「クエリに対する処理はどのようにすればよいか」を考える必要があります。この問題のように「1点に対する変更クエリ」と「区間に対する質問クエリ」に対応した平方分割を書くときには,考えた処理が以下の3つの条件を満たすかどうか確認するようにしています。クエリ処理はこれらの条件を満たさなければならないという前提を持って考えると,必要な補助データ構造やクエリの処理が思いつきやすいです。

  • 変更クエリはO(√N)か?
  • 変更クエリのあとで,1つのバケット全体に対する質問クエリにO(1)で答えられるか?
  • 変更クエリのあとで,バケットの中に対する質問クエリにO(√N)で答えられるか?

これらの条件を満たす処理が実現できれば,バケットの個数は√N個なので,どのような質問クエリに対してもO(1)×√N+O(√N)×2で全体としてO(√N)です。

この条件を満たさなければならないことを注意すると,私は次のように考えます。

バケットの中に対する質問クエリにO(√N)で答えるという制約は単純にループを回せば達成できるけど,1つのバケット全体に対する質問クエリにO(1)で答えないといけないという制約が厳しいな。補助データ構造を作って,1つのバケット全体の合計値をどこかのタイミングで保存しておくことができればO(1)で答えられそうだ。そういえば,変更クエリは単純な方法だとO(1)でできるけどO(√N)時間かけても大丈夫みたいだな。O(√N)時間かけたら何ができるだろう? このタイミングでバケット全体の合計値を再計算するとO(√N)時間だからうまくいきそうだ。

ここでは,次のように処理することでクエリに対応することを考えましょう。

区間を√N個のバケットに分割し,バケット中の合計値を bucketSum という配列に保存する。値の変更時には,対応するバケットの合計値を再計算して bucketSum を更新する。質問クエリでは,バケット全体に対しては bucketSum の値を返し,バケットの中に対しては単純にループを回して合計値を計算する。

変更クエリ add(7,2) の処理の例を見てみましょう。クエリ処理前の状態が以下のようなものだとします。

f:id:kujira16:20161214221706p:plain

値を変更します。

f:id:kujira16:20161214221712p:plain

対応するバケットの合計値を計算します。

f:id:kujira16:20161214221724p:plain

bucketSum を更新します。

f:id:kujira16:20161214221733p:plain

getSum(2,9) という質問クエリでは,[4,7]の区間は「バケット全体に対する質問クエリ」,[2,3], [8,9]の区間は「バケットの中に対する質問クエリ」ですので,sum(2, 6, 16, 8, 4) を計算して合計値は36です。

f:id:kujira16:20161215140334p:plain

先ほど示した条件を満たしているか確認しましょう。

  • 変更クエリはO(√N)か?

1個のバケットに含まれる要素の数はO(√N)ですので,対応するバケットの合計値を再計算して bucketSum を更新する処理はO(√N)です。

  • 変更クエリのあとで,1つのバケット全体に対する質問クエリにO(1)で答えられるか?

バケット全体に対しては bucketSum の値を返せば良いだけなのでO(1)です。

  • 変更クエリのあとで,バケットの中に対する質問クエリにO(√N)で答えられるか?

単純にループを回して合計値を計算することでO(√N)です。

ソースコード例は以下のようになります。

http://judge.u-aizu.ac.jp/onlinejudge/review.jsp?rid=2123752

#define int long long
const int sqrtN = 512;
struct SqrtDecomposition {
  int N, K;
  vector<int> data;
  vector<int> bucketSum;
  SqrtDecomposition(int n) : N(n) {
    K = (N + sqrtN - 1) / sqrtN;
    data.assign(K * sqrtN, 0);
    bucketSum.assign(K, 0);
  }
  void add(int x, int y) {
    int k = x / sqrtN;
    data[x] += y;
    int sum = 0;
    for (int i = k * sqrtN; i < (k + 1) * sqrtN; ++i) {
      sum += data[i];
    }
    bucketSum[k] = sum;
  }
  // [x, y)
  int getSum(int x, int y) {
    int sum = 0;
    for (int k = 0; k < K; ++k) {
      int l = k * sqrtN, r = (k + 1) * sqrtN;
      if (r <= x || y <= l)
        continue;
      if (x <= l && r <= y) {
        sum += bucketSum[k];
      } else {
        for (int i = max(x, l); i < min(y, r); ++i) {
          sum += data[i];
        }
      }
    }
    return sum;
  }
};

Range Minimum Query

http://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=DSL_2_A

この問題では2種類のクエリに対応する必要があります。

  • update(i,x): a[i]をxに変更する。
  • find(s,t): a[s],a[s+1],…,a[t]の最小値を出力する。

ここでは,次のように処理することでクエリに対応することを考えましょう。

区間を√N個のバケットに分割し,バケット中の最小値を bucketMin という配列に保存する。値の変更時には,対応するバケットの最小値を再計算して bucketMin を更新する。質問クエリでは,バケット全体に対しては bucketMin の値を返し,バケットの中に対しては単純にループを回して最小値を計算する。

変更クエリ update(9,1) の処理の例を見てみましょう。クエリ処理前の状態が以下のようなものだとします。

f:id:kujira16:20161214183800p:plain

値を変更します。

f:id:kujira16:20161214183809p:plain

対応するバケットの最小値を計算します。

f:id:kujira16:20161214183815p:plain

bucketMin を更新します。

f:id:kujira16:20161214183822p:plain

find(6,13) という質問クエリでは,[8,11]の区間は「バケット全体に対する質問クエリ」,[6,7], [12,13]の区間は「バケットの中に対する質問クエリ」ですので,min(3, 0, 1, 9, 2) を計算して最小値は0です。

f:id:kujira16:20161215140614p:plain

先ほど示した条件を満たしているか確認しましょう。

  • 変更クエリはO(√N)か?

1個のバケットに含まれる要素の数はO(√N)ですので,対応するバケットの最小値を再計算して bucketMin を更新する処理はO(√N)です。

  • 変更クエリのあとで,1つのバケット全体に対する質問クエリにO(1)で答えられるか?

バケット全体に対しては bucketMin の値を返せば良いだけなのでO(1)です。

  • 変更クエリのあとで,バケットの中に対する質問クエリにO(√N)で答えられるか?

単純にループを回して最小値を計算することでO(√N)です。

ソースコード例は以下のようになります。

http://judge.u-aizu.ac.jp/onlinejudge/review.jsp?rid=2123670

#define int long long
const int INF = (1LL << 31) - 1;
const int sqrtN = 512;
struct SqrtDecomposition {
  int N, K;
  vector<int> data;
  vector<int> bucketMin;
  SqrtDecomposition(int n) : N(n) {
    K = (N + sqrtN - 1) / sqrtN;
    data.assign(K * sqrtN, INF);
    bucketMin.assign(K, INF);
  }
  void update(int x, int y) {
    int k = x / sqrtN;
    int minVal = data[x] = y;
    for (int i = k * sqrtN; i < (k + 1) * sqrtN; ++i) {
      minVal = min(minVal, data[i]);
    }
    bucketMin[k] = minVal;
  }
  // [x, y)
  int find(int x, int y) {
    int minVal = INF;
    for (int k = 0; k < K; ++k) {
      int l = k * sqrtN, r = (k + 1) * sqrtN;
      if (r <= x || y <= l)
        continue;
      if (x <= l && r <= y) {
        minVal = min(minVal, bucketMin[k]);
      } else {
        for (int i = max(x, l); i < min(y, r); ++i) {
          minVal = min(minVal, data[i]);
        }
      }
    }
    return minVal;
  }
};

区間に対する変更クエリ・1点に対する質問クエリ

Range Add Query

http://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=DSL_2_E

  • add(s,t,x): a[s],a[s+1],…,a[t] にxを加算する。
  • get(i): a[i]の値を出力する。

この問題のように「区間に対する変更クエリ」と「1点に対する質問クエリ」に対応した平方分割を書くときには,考えた処理が以下の3つの条件を満たすかどうか確認するようにしています。

  • 1つのバケット全体に対する変更クエリはO(1)か?
  • バケットの中に対する変更クエリはO(√N)か?
  • 変更クエリのあとで,質問クエリにO(√N)で答えられるか?

この条件を満たさなければならないことを注意すると,次のように考えます。

バケットの中に対する変更クエリをO(√N)で行うという制約は単純にループを回す方法で問題ないけど,1つのバケット全体に対する変更クエリをO(1)で行うという制約が厳しいから,バケット全体に対する変更は補助データ構造だけに書き込むようにしないといけないな。ということは,バケット全体に足された値を補助データ構造に持っておけば良さそうだな。get(i)を処理する時に,iを含むバケット全体に足された値を追加で加算するようにすれば答えも合いそうだ!

ここでは,次のように処理することでクエリに対応することを考えましょう。

区間を√N個のバケットに分割し,バケット中のすべての要素に加算される値は bucketAdd という配列に保存する。質問クエリでは,その要素に加算された値と,対応するバケットの bucketAdd の値の和を答えにする

変更クエリ add(9,10,2) の処理の例を見てみましょう。クエリ処理前の状態が以下のようなものだとします。

f:id:kujira16:20161215142450p:plain

この区間への変更は「「バケットの中に対する変更クエリ」なので単純にループを回します。

f:id:kujira16:20161215142456p:plain

get(i) クエリでは data[i] の値と,対応するバケットのbucketAddの値を足した値を答えにします。

f:id:kujira16:20161215142501p:plain

続いて,変更クエリ add(3,8,1) の処理の例を見てみましょう。[4,7]の区間は「バケット全体に対する変更クエリ」なので,bucketAdd のみに書き込みます。[3,3], [8,8] の区間は「バケットの中に対する変更クエリ」なので単純にループを回します。

f:id:kujira16:20161215142504p:plain

get(i) クエリでは data[i] の値と,対応するバケットのbucketAddの値を足した値を答えにします。

f:id:kujira16:20161215142507p:plain

先ほど示した条件を満たしているか確認しましょう。

  • 1つのバケット全体に対する変更クエリはO(1)か?

bucketAdd に加算するだけなのでO(1)です。

  • バケットの中に対する変更クエリはO(√N)か?

バケット内の要素の数はO(√N)なので,単純にループを回してもO(√N)です。

  • 変更クエリのあとで,質問クエリにO(√N)で答えられるか?

databucketAddの値を足すだけなのでO(1)です。

ソースコード例は以下のようになります。

http://judge.u-aizu.ac.jp/onlinejudge/review.jsp?rid=2123808

#define int long long
const int sqrtN = 512;
struct SqrtDecomposition {
  int N, K;
  vector<int> data;
  vector<int> bucketAdd;
  SqrtDecomposition(int n) : N(n) {
    K = (N + sqrtN - 1) / sqrtN;
    data.assign(K * sqrtN, 0);
    bucketAdd.assign(K, 0);
  }
  // [s, t)
  void add(int s, int t, int x) {
    for (int k = 0; k < K; ++k) {
      int l = k * sqrtN, r = (k + 1) * sqrtN;
      if (r <= s || t <= l)
        continue;
      if (s <= l && r <= t) {
        bucketAdd[k] += x;
      } else {
        for (int i = max(s, l); i < min(t, r); ++i) {
          data[i] += x;
        }
      }
    }
  }
  int get(int i) {
    int k = i / sqrtN;
    return data[i] + bucketAdd[k];
  }
};

Range Update Query

http://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=DSL_2_D

この問題では2種類のクエリに対応する必要があります。

  • update(s,t,x): a[s],a[s+1],…,a[t] をxに変更する。
  • find(i): a[i] の値を出力する。

今回の「区間の値をまとめて変更する」という操作のように,操作の種類によっては「遅延伝搬」や「遅延評価」「遅延更新」「遅延セグメント木」などと呼ばれるテクニックが必要になることがあります*1。もう少し詳しく言うと,「操作の順番を入れ替えると結果が変わってしまう操作」では遅延伝搬テクが必要です。

可算は操作の順番を入れ替えても結果が変わりません。+3してから+5しても,+5してから+3しても同じです。そのため「区間に対する可算」では遅延伝搬テクを使わなくても解けます*2

一方で,「値をまとめて変更する」といった操作では操作の順番を入れ替えると結果が変わってしまいます。値を3に変更してから5に変更するのと,5に変更してから3に変更するのとでは結果が変わってしまいます。

今回の場合,updateクエリは「操作の順番を入れ替えると結果が変わってしまう操作」に該当するので,遅延伝搬テクを使います。

この場合も「区間に対する変更クエリ」と「1点に対する質問クエリ」に対応しなければならないという点では同じなので,やはり以下の3つの条件は満たす必要があります。

  • 1つのバケット全体に対する変更クエリはO(1)か?
  • バケットの中に対する変更クエリはO(√N)か?
  • 変更クエリのあとで,質問クエリにO(√N)で答えられるか?

そのため,以下のように考えます。

1つのバケット全体に対する変更クエリをO(1)で行うという制約が厳しいからバケット全体に対する変更は補助データ構造だけに書き込むようにしないといけないな。でも data の内容と矛盾が生じてしまうから,どこかのタイミングで data の内容と同期を取らないといけないな。特に,質問クエリと「バケットの中に対する変更クエリ」の前には同期を取っておかないと大変なことになりそうだ。1つのバケットについて同期を取るのはO(√N)時間だから,質問クエリと「バケットの中に対する変更クエリ」の直前に,対応するバケットについてのみ同期を取れば条件を満たせそうだな。

ここでは,次のように処理することでクエリに対応することを考えましょう。

区間を√N個のバケットに分割し,バケット中のすべての要素の値を変更するときは lazyUpdate という配列に保存する。質問クエリやバケットの中に対する変更クエリを処理するときには,対応するバケットについてのみ lazyUpdate を伝搬させてから処理を行う。

例を見てみましょう。初期状態が以下のようなものだとします。

f:id:kujira16:20161215154045p:plain

update(9,10,2) はバケットの中に対する変更クエリですが,対応するバケットの lazyUpdate には値が入っていないので,data をそのまま変更します。

f:id:kujira16:20161215154051p:plain

get(10) を処理したいのですが,lazyUpdate には値が入っていないので data の値を返すだけです。

f:id:kujira16:20161215154056p:plain

続いて,update(2,14,3) について考えてみましょう。[4,7], [8, 11] の区間はバケット全体に対する変更クエリなので,変更後の値を lazyUpdate に書き込みます。[2, 3], [12, 14] の区間はバケットの中に対する変更クエリなので,data をそのまま変更します。

f:id:kujira16:20161215154100p:plain

get(5) を処理したいとしましょう。今回は lazyUpdate に値が入っているので,値を伝搬させる必要があります。

f:id:kujira16:20161215154103p:plain

対応しているバケットについて lazyUpdate の値を data に伝搬させました。

f:id:kujira16:20161215154106p:plain

あとは data の値を返すだけです。

f:id:kujira16:20161215154109p:plain

「バケットの中に対する変更クエリ」の場合にも伝搬は必要です。update(6,10,4) を処理したい場合,[8,10] の区間はバケットの中に対する変更クエリですが,lazyUpdate に値が入っているので値を伝搬させる必要があります。

f:id:kujira16:20161215154112p:plain

対応しているバケットについて lazyUpdate の値を data に伝搬させました。

f:id:kujira16:20161215154115p:plain

あとは data をそのまま変更するだけです。

f:id:kujira16:20161215154119p:plain

3つの条件を満たしているか確認しましょう。

  • 1つのバケット全体に対する変更クエリはO(1)か?

lazyUpdate に書き込むだけなので O(1) です。

  • バケットの中に対する変更クエリはO(√N)か?

lazyUpdate に値が入っていたら値を伝搬させる必要がありますが,これはO(√N)時間かかります。data の変更は単純にループを回してO(√N)です。合計するとO(√N)です。

  • 変更クエリのあとで,質問クエリにO(√N)で答えられるか?

lazyUpdate に値が入っていたら値を伝搬させる必要がありますが,これはO(√N)時間かかります。値の取得はO(1)です。合計するとO(√N)です。

ソースコード例は以下のようになります。

http://judge.u-aizu.ac.jp/onlinejudge/review.jsp?rid=2124136

#define int long long
const int INF = (1LL << 31) - 1;
const int sqrtN = 512;
struct SqrtDecomposition {
  int N, K;
  vector<int> data;
  vector<bool> lazyFlag;
  vector<int> lazyUpdate;
  SqrtDecomposition(int n) : N(n) {
    K = (N + sqrtN - 1) / sqrtN;
    data.assign(K * sqrtN, INF);
    lazyFlag.assign(K, false);
    lazyUpdate.assign(K, 0);
  }
  void eval(int k) {
    if (lazyFlag[k]) {
      lazyFlag[k] = false;
      for (int i = k * sqrtN; i < (k + 1) * sqrtN; ++i) {
        data[i] = lazyUpdate[k];
      }
    }
  }
  // [s, t)
  void update(int s, int t, int x) {
    for (int k = 0; k < K; ++k) {
      int l = k * sqrtN, r = (k + 1) * sqrtN;
      if (r <= s || t <= l)
        continue;
      if (s <= l && r <= t) {
        lazyFlag[k] = true;
        lazyUpdate[k] = x;
      } else {
        eval(k);
        for (int i = max(s, l); i < min(t, r); ++i) {
          data[i] = x;
        }
      }
    }
  }
  int find(int i) {
    int k = i / sqrtN;
    eval(k);
    return data[i];
  }
};

区間に対する変更クエリ・区間に対する質問クエリ

RSQ and RAQ

http://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=DSL_2_G

この問題では2種類のクエリに対応する必要があります。

  • add(s,t,x): a[s],a[s+1],…,a[t] にxを加算する。
  • getSum(s,t): a[s],a[s+1],…,a[t] の合計値を出力する。

遅延伝搬のところで説明しましたが,加算は操作の順番を入れ替えても結果が変わらないので,遅延伝搬を使わなくても解けます。今回は遅延伝搬を使わない方法で解いてみましょう。

この問題のように「区間に対する変更クエリ」と「区間に対する質問クエリ」に対応した平方分割を書くときには,考えた処理が以下の4つの条件を満たすかどうか確認するようにしています。

  • バケットの中に対する変更クエリはO(√N)か?
  • バケットの中に対する質問クエリにO(√N)で答えられるか?
  • 1つのバケット全体に対する変更クエリはO(1)か?
  • 1つのバケット全体に対する質問クエリにO(1)で答えられるか?

そのため,以下のように考えます。

1つのバケット全体に対する変更クエリをO(1)で行うという制約が厳しいからバケット全体に加算された値は補助データ構造だけに書き込むようにしないといけないな。あと,1つのバケット全体に対する質問クエリにO(1)で答えるという制約もあるから,1つのバケット全体の合計値も別の補助データ構造に保存しておく必要があるな。

というわけで,補助データ構造が2つ必要です。バケット全体に加算された値は bucketAdd,バケット全体の合計値は bucketSum という配列に格納することにしましょう。

例を見てみましょう。初期状態が以下のようなものだとします。

f:id:kujira16:20161215173641p:plain

add(3,8,2) を処理します。[4,7] の区間はバケット全体に対する変更,[3,3], [8,8] の区間はバケットの中に対する変更です。バケット全体に対する変更はO(1)である必要があるので bucketAdd だけに書き込みます。バケットの中に対する変更クエリはO(√N)かけてよいので単純にループを回して処理しますが,あとでバケット全体に対する質問クエリにO(1)で答えるために,バケット全体の合計値の再計算も行う必要があります。

f:id:kujira16:20161215173644p:plain

getSum(6,7) を処理しましょう。これはバケットの中に対する質問なのでO(√N)かけてよく,単純にループを回せば良いです。ただし,どの要素にも bucketAdd の値が加算されているので,bucketAddの値×要素数も含める必要があります。

f:id:kujira16:20161215173648p:plain

最後に getSum(0,8) を処理しましょう。[0,3], [4,7] の区間はバケット全体に対する質問,[8,8]の区間はバケットの中に対する質問です。バケット全体に対する質問はO(1)で答える必要があるので,bucketSum の値を使います。また,バケット内の要素には bucketAdd の値が加算されているので,bucketAddの値×要素数も含める必要があります。

f:id:kujira16:20161215173652p:plain

条件を満たしているか確認しましょう。

  • バケットの中に対する変更クエリはO(√N)か?

data を変更したあと,対応するバケットの bucketSum を再計算する必要があるのでO(√N)です。

  • バケットの中に対する質問クエリにO(√N)で答えられるか?

data の合計値の計算がO(√N), bucketAdd×要素数の計算がO(1)なのでO(√N)です。

  • 1つのバケット全体に対する変更クエリはO(1)か?

bucketAddに値を足すだけなのでO(1)です。

  • 1つのバケット全体に対する質問クエリにO(1)で答えられるか?

bucketSumbucketAdd×要素数の計算をするだけなのでO(1)です。

ソースコード例は以下のようになります。

http://judge.u-aizu.ac.jp/onlinejudge/review.jsp?rid=2124559

#define int long long
const int sqrtN = 512;
struct SqrtDecomposition {
  int N, K;
  vector<int> data;
  vector<int> bucketSum;
  vector<int> bucketAdd;
  SqrtDecomposition(int n) : N(n) {
    K = (N + sqrtN - 1) / sqrtN;
    data.assign(K * sqrtN, 0);
    bucketSum.assign(K, 0);
    bucketAdd.assign(K, 0);
  }
  // [s, t)
  void add(int s, int t, int x) {
    for (int k = 0; k < K; ++k) {
      int l = k * sqrtN, r = (k + 1) * sqrtN;
      if (r <= s || t <= l)
        continue;
      if (s <= l && r <= t) {
        bucketAdd[k] += x;
      } else {
        for (int i = max(s, l); i < min(t, r); ++i) {
          data[i] += x;
          bucketSum[k] += x;
        }
      }
    }
  }
  // [s, t)
  int getSum(int s, int t) {
    int sum = 0;
    for (int k = 0; k < K; ++k) {
      int l = k * sqrtN, r = (k + 1) * sqrtN;
      if (r <= s || t <= l)
        continue;
      if (s <= l && r <= t) {
        sum += bucketSum[k] + bucketAdd[k] * sqrtN;
      } else {
        for (int i = max(s, l); i < min(t, r); ++i) {
          sum += data[i] + bucketAdd[k];
        }
      }
    }
    return sum;
  }
};

RMQ and RUQ

http://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=DSL_2_F

この問題では2種類のクエリに答える必要があります。

  • update(s,t,x): a[s],a[s+1],…,a[t] をxに変更する。
  • find(s,t): a[s],a[s+1],…,a[t]の最小値を出力する。

「値の変更」というクエリがありますが,これは操作の順番を入れ替えると結果が変わってしまうので遅延伝搬を使うと解きやすいです。

必要な条件を確認しましょう。

  • バケットの中に対する変更クエリはO(√N)か?
  • バケットの中に対する質問クエリにO(√N)で答えられるか?
  • 1つのバケット全体に対する変更クエリはO(1)か?
  • 1つのバケット全体に対する質問クエリにO(1)で答えられるか?

以下のように考えます。

1つのバケット全体に対する変更クエリがO(1)になる必要があるから,変更後の値は補助データ構造だけに書き込むようにしないといけないな。1つのバケット全体に対する質問クエリにもO(1)で答えないといけないから,どこかのタイミングでバケットの最小値を別の補助データ構造に書き込む必要があるな。

この場合もやはり補助データ構造が2つ必要です。バケット全体の変更後の値は lazyUpdate,バケット全体の最小値は bucketMin という配列に格納することにしましょう。

例を見てみましょう。初期状態が以下のようなものだとします。

f:id:kujira16:20161215200534p:plain

update(2,9,5) というクエリを処理しましょう。[4,7]の区間はバケット全体に対する変更,[2,3], [8,9]の区間はバケットの中に対する変更です。

バケット全体に対する変更はO(1)である必要があるので lazyUpdate だけに書き込みます。また,あとでバケット全体に対する質問クエリにO(1)で答えるためにbucketMinも書き換えておきましょう。

バケットの中に対する変更クエリはO(√N)かけてよいので単純にループを回して処理しますが,あとでバケット全体に対する質問クエリにO(1)で答えるために,bucketMin の再計算も行う必要があります。

f:id:kujira16:20161215200541p:plain

find(0,8) というクエリを処理しましょう。[0,3], [4,7]の区間はバケット全体に対する質問なので bucketMin を返しておけばよく,[8,8]の区間はバケットの中に対する質問なのでO(√N)でループを回して最小値を計算すればよいです。

f:id:kujira16:20161215200547p:plain

次に update(6,9,2) というクエリを処理しましょう。[6,7], [8,9]の区間はバケットの中に対する質問ですが,[6,7]の区間に対応する updateLazy に値が入っているので値を伝搬させる必要があります。[8,9]の区間updateLazy に値が入っていないので値の伝搬は必要ありません。

f:id:kujira16:20161215200550p:plain

updateLazy を伝搬させました。

f:id:kujira16:20161215200554p:plain

単純にループを回して値の変更を行います。bucketMinの再計算を忘れないようにしてください。

f:id:kujira16:20161215200556p:plain

条件を満たしているか確認しましょう。

  • バケットの中に対する変更クエリはO(√N)か?

lazyUpdateに値が入っていれば,値を伝搬させる必要があります。これはO(√N)です。その後,単純にループを回して値の変更を行ってから,bucketMinの再計算を行います。これはO(√N)です。合計するとO(√N)です。

  • バケットの中に対する質問クエリにO(√N)で答えられるか?

lazyUpdateに値が入っていれば,値を伝搬させる必要があります。これはO(√N)です。その後,単純にループを回して最小値を計算するのにO(√N)かかります。合計するとO(√N)です。

  • 1つのバケット全体に対する変更クエリはO(1)か?

lazyUpdateに値を書き込むだけなのでO(1)です。

  • 1つのバケット全体に対する質問クエリにO(1)で答えられるか?

bucketMinを参照するだけなのでO(1)です。

ソースコード例は以下のようになります。

http://judge.u-aizu.ac.jp/onlinejudge/review.jsp?rid=2124563

#define int long long
const int INF = (1LL << 31) - 1;
const int sqrtN = 512;
struct SqrtDecomposition {
  int N, K;
  vector<int> data;
  vector<int> bucketMin;
  vector<bool> lazyFlag;
  vector<int> lazyUpdate;
  SqrtDecomposition(int n) : N(n) {
    K = (N + sqrtN - 1) / sqrtN;
    data.assign(K * sqrtN, INF);
    bucketMin.assign(K, INF);
    lazyFlag.assign(K, false);
    lazyUpdate.assign(K, 0);
  }
  void eval(int k) {
    if (lazyFlag[k]) {
      lazyFlag[k] = false;
      for (int i = k * sqrtN; i < (k + 1) * sqrtN; ++i) {
        data[i] = lazyUpdate[k];
      }
    }
  }
  // [s, t)
  void update(int s, int t, int x) {
    for (int k = 0; k < K; ++k) {
      int l = k * sqrtN, r = (k + 1) * sqrtN;
      if (r <= s || t <= l)
        continue;
      if (s <= l && r <= t) {
        bucketMin[k] = x;
        lazyFlag[k] = true;
        lazyUpdate[k] = x;
      } else {
        eval(k);
        for (int i = max(s, l); i < min(t, r); ++i) {
          data[i] = x;
        }
        int &minVal = bucketMin[k] = INF;
        for (int i = l; i < r; ++i) {
          minVal = min(minVal, data[i]);
        }
      }
    }
  }
  // [s, t)
  int find(int s, int t) {
    int minVal = INF;
    for (int k = 0; k < K; ++k) {
      int l = k * sqrtN, r = (k + 1) * sqrtN;
      if (r <= s || t <= l)
        continue;
      if (s <= l && r <= t) {
        minVal = min(minVal, bucketMin[k]);
      } else {
        eval(k);
        for (int i = max(s, l); i < min(t, r); ++i) {
          minVal = min(minVal, data[i]);
        }
      }
    }
    return minVal;
  }
};

応用例

Flipping Parentheses

http://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=1351

平方分割に関連するところだけ解説します。この問題では次の2種類のクエリに答える必要があります。

  • add(s,t,x): a[s],a[s+1],…,a[t] にxを加算する。
  • find(i): a[j],a[j+1],…,a[i]の最小値が2未満となる最大のjを求める。

区間に対する変更クエリ」と「区間に対する質問クエリ」が必要だということがわかります。ブロック全体に対する加算をO(1)で行うための補助データ構造 bucketAdd と,ブロック全体の最小値をO(1)で求めるための補助データ構造 bucketMin が必要です。

add(s,t,x)クエリの処理は特に解説するところがないので省略します。問題はfind(i)クエリの処理で,[j,i]の最小値を求める関数を作ってjを二分探索するような方法を使うと平方分割では O(√N・log N) かかってしまいます。N=3*105を代入してみると104くらいになってしまい,間に合いそうにありません。

find(i)を効率的に処理するためには,バケットを右から左に走査していき,明らかに条件を満たさないバケットはスキップするという方法を使います。この方法によりO(√N)時間で処理を実現できます。いま,以下のような状態で find(12) を求めたいとしましょう。

f:id:kujira16:20161215211337p:plain

bucketMinを見ることにより,[8,11] の区間の最小値は2だということが分かります。求めたいのは「a[j…i]の最小値が2未満となる最小のj」ですが,[8,11]の左端からiまでの最小値はmin([8,12])=2ですので,[8,11]の区間内に「a[j…i]の最小値が2未満」となるものは存在しないことが分かります。よって,[8,11]のバケットはスキップできます。

次に[4,7] の区間を見ると,bucketMinを見ることによりmin([4,7])=1であることが分かります。よって,この区間内に「a[j…i]の最小値が2未満」となるものが存在することがわかります。よって,このバケットを右から左へ走査して,最小値が2未満となるjを見つければよいです。

f:id:kujira16:20161215211343p:plain

バケットの左端を調べる処理はO(√N)回,「a[j…i]の最小値が2未満」となるものが存在するバケットの中を詳しく調べる処理にO(√N)時間かかるので,findクエリは合計でO(√N)時間で行えます。

セグメントツリーにチャレンジしたくなったら

@pekempeyさんが分かりやすい可視化をしてくれています。

https://pekempey.github.io/lazy_propagation_segment_tree/lazy_propagation_segment_tree.html

アドベントカレンダーで同じ日の @kobae964 さんが,セグメントツリーに適用可能な演算についてのサーベイを書いてくれています。

koba-e964.hatenablog.com

謝辞

記事を確認していただいた @yazaten さん,@yurahuna さんには多くの有意義な助言をいただきましたことを感謝いたします。

*1:実は遅延更新を使わなくても実現可能だそうです https://pekempey.github.io/lazy_propagation_segment_tree/lazy_propagation_segment_tree.html

*2:使っても解けます

Jupyter Notebookの小ネタ (display, tqdm, SSHポート転送)

この記事は jupyter notebook Advent Calendar 2016 の12月13日の記事です。

qiita.com

Jupyter Notebookを使う上で,この機能を知らない人を見つけたら教えてあげたくなるような小ネタをいくつか紹介します。

IPython.display.display

Jupyter Notebook上でpandasのDataFrameを表示すると,HTMLのtable要素として見やすく表示してくれます。活用している方も多いかと思います。

df[df['f'] == 'C']df[df['f'] == 'T'] を表示したいのですが,Jupyter Notebookで表示されるのはそれぞれのセルで最後に評価されたものだけです。そのため,表示したいDataFrameごとにセルを分ける必要があります。

gist.github.com

しかしながら,このように複数のセルに分けてしまうと,DataFrameに前処理を行った場合などにいくつものセルを再実行するのは面倒です。このような場合には IPython.display.display を使うのが便利です。

gist.github.com

以下のページで存在を知りました。

stackoverflow.com

tqdm

Iterableを渡すだけでプログレスバーを表示してくれるtqdmというPythonのライブラリがあります。

https://tqdm.github.io/

コマンドライン上で使うとこんな感じです。

import time
from tqdm import tqdm

s = 0
for i in tqdm(range(20)):
    s += i
    time.sleep(0.1)
print(s)

www.youtube.com

このライブラリはコマンドライン上でしか使えないと思い込んでいたのですが,ドキュメントをよく読むとIPython/Jupyter Integrationなるものが追加されていました。これは以下のように使えます。

import time
from tqdm import tqdm_notebook as tqdm

s = 0
for i in tqdm(range(20)):
    s += i
    time.sleep(0.1)
print(s)

www.youtube.com

SSHポート転送

知っている人にとっては「それ当たり前じゃん!」となってしまうかもしれませんが,小ネタが2つだけだと物足りなかったので紹介させてください。

手元のコンピュータの性能が貧弱な場合には,クラウド上のサーバーや研究室の計算機などでJupyter Notebookを使いたいことがあります。

しかしながら,何も設定していない状態だと,Jupyter Notebookを起動したサーバーのIPアドレスの8888ポートをブラウザに打ち込んでも接続できないと思います。

これはJupyter Notebookのデフォルト設定で127.0.0.1以外に8888ポートを公開していないことが原因なのですが,これを公開設定にするだけではセキュリティ的に問題があります。そのため,接続元のIPアドレスに制限をかけたり,パスワードをかけたりするための設定が必要です。扱っているデータによってはSSL/TLSの設定も必要かもしれません。これらの設定については以下のページで解説されています。

Running a notebook server — Jupyter Notebook 5.0.0.dev documentation

しかしながら,これらの設定を行うのはちょっと面倒です。もっと簡単な方法として,SSHのポート転送機能を使う方法があります。私はSSHでサーバに接続するとき,以下のようにしています。

$ ssh -L8888:localhost:8888 username@hostname

-L8888:localhost:8888 の部分は,手元のコンピュータの8888ポートへの接続を,SSHでの接続先から見てlocalhostの8888ポートに転送するという設定です。SSHでの接続先でJupyter Notebookを起動してから,ブラウザで手元のコンピュータのlocalhost:8888に接続すれば繋がります。Jupyter Notebookから見ればlocalhostからのアクセスとなるためアクセス制限にも引っかかりません。

サーバー上でJupyter Notebookを使う際には,ほかにも以下のようなことに気をつけています。

  • 「ブラウザが起動できませんでした」という警告が出るのが目障りなので,jupyter notebook --no-browser としてブラウザが起動しないようにしています。
  • SSHの接続が切断されてしまうとJupyter Notebookが停止してしまい,再接続しても復活しないので,何かしらの工夫が必要です。私の場合はターミナルマルチプレクサとしてtmuxを使うことでこれを防いでいます。tmuxをインストールすることが難しい環境では nohup & を使うのが良いでしょう(参考: http://www.glamenv-septzen.net/view/854)。

ACM-ICPC 2016アジア地区つくば大会 参加記

結果

5完26位でした。

Results and Standings | ACM-ICPC 2016 Asia Tsukuba Regional

メンバー

会場に行くまで

奈良に新幹線なんてあるわけないじゃん!

移動時間の見積もりで,1時間に2本のバスに乗り遅れたときのバッファとみどりの窓口が混雑していたときのバッファを乗せていたら,13:00時到着の予定で6時30分に家を出ることになりました。結局バッファを使うことがなかったので,12時すぎに到着したのですが,つくば ↔︎ けいはんな に5時間半かかるのはちょっと憂鬱ですね…

駅前では筑波大の研究室のロボットが動いていました。

練習セッション

事前に System Trial DVD が提供されていたので,滞りなく終わりました。

練習問題を解いたあとは印刷を試してみたりClarを試してみたり "I want to go to restroom." の練習をしてみたりしました。

そのあとは実行環境のスタックサイズを確認したり,ループがどれくらいの回数回せるかを確認したかったのですが,キーボードに慣れるためにタイピング練習をしていたら練習セッションの時間が終わってしまいました。去年までのJava Challengeの時間はキーボード練習の時間を兼ねていた考えると,練習セッションの時間はあと15分くらいは欲しかったかもしれません。

去年のキーボードは配列がイマイチで選手からも不評だったのですが,今年のキーボードはアンケートの内容を反映してもらえたのか,去年よりも拡大に快適でした。

Welcome Party

会場はアヒルボートの博士号がある池に近い場所だったらしかったので見に行けたら良かったのですが,良いタイミングが見つからず見に行けませんでした。

宿泊

東横イン最高!!去年はあまり眠れず競技時間中に頭が回らなくて困ったことになったのですが,今年は快適でした。

コンテスト

A

解いた

B

時間がかかったけど解いた

C

その位置に到達できる製造ラインの番号のminとmaxを更新していくDPみたいなものをめんふぃむが提案してくれて10分くらいで解けた。

D

基数をBとして

(aの出現回数*B0 + bの出現回数*B1 + cの出現回数*B2 + … + zの出現回数*B25) % M

をハッシュにした。ハッシュ値をstd::set<int64_t> に格納して O(N2 log N) にするとTLEするし,ハッシュをstd::unordered_setに格納して O(N2) にするとハッシュがコンフリクトしてWAになるしペナルティがたくさん加算されて辛い思いをした。O(N2 log N) の定数倍を改善すると通った。想定解は O(N2) らしい。ふーん…

↓ AOJでは通るけど本番環境では通るか自信なし

https://gist.github.com/arosh/61b3753f1d747b6f251e90f34bac2299#file-d-cpp

G

ビットの配列の上で多倍長の足し算を実装できればいいね,ということで,以下のクエリを実現するデータ構造を平方分割で実装するという愚かな方針を選んでしまった。

  • set1(p) … pビット目を1にする
  • set0range(l,r) … [l,r) の範囲のビットを全て0にする
  • search(p) … pよりも左側にあって,最も右にある0の位置を返す
  • check0(p) … pよりも右側のビットが全て0であるか確認する

常識的に考えて,他のチームはかなり早く通しているのだから,もっと簡単な方針がないか検討するべきだった。

↓実際に通したコード

https://gist.github.com/arosh/61b3753f1d747b6f251e90f34bac2299#file-g-cpp

E

40分しか残っていない状況で実装を始めた。

終了数分前に -1+1 を -(1+1) と解釈してしまうバグが見つかり,どこを直せば良いのかすぐに見つけられず無念のまま競技時間が終わってしまい,E(p)F(p) に直せば良いだけじゃん,ということがわかって絶望していたのだが,そのほかにも 01 みたいな数値リテラルは無効にするという処理を入れていなかったので,1バイト直せばACするというのは嘘だった。

↓ あと数分あれば通せたコード

https://gist.github.com/arosh/61b3753f1d747b6f251e90f34bac2299#file-e-cpp

@@ -113,7 +113,7 @@ Result F(iter &p) {
   }
   else if(*p == '-') {
     ++p;
-    Result e = E(p);
+    Result e = F(p);
     if(!e.valid) return fail;
     e.value *= -1;
     return e;
@@ -125,10 +125,13 @@ Result N(iter &p) {
     return fail;
   }
   int n = 0;
+  int i = 0;
   while(isdigit(*p)) {
     n *= 2;
+    if(i > 0 && n == 0) return fail;
     n += *p - '0';
     ++p;
+    ++i;
   }
   // DEBUG(n);
   return Result(true, n);

表彰式

13*2位という理由で(?)IBMから副賞を貰いました。順位はイマイチだったけど,今までICPCで一度も賞を貰ったことがなかったので,最後にもらえて嬉しいです。

最後に

5年前にICPCの出場歴が無い大学から参加を始めて,最初は国内予選60位という順位でしたが,最終的にアジア地区大会に4回も出ることができて良かったです。

続き