接尾辞配列とその応用

接尾辞配列

単一文字列の部分文字列に関する辞書順問題に一般的に利用できます。

アルゴリズムの流れ

まず、文字列のすべての接尾辞をソートします。

定義:\(sa[i]\) はすべての接尾辞をソートした後、第 \(i\) 番目に小さい接尾辞のインデックスを表します。\(rk[i]\) は接尾辞 \(i\) のランクを表します。

倍増法と基数ソート(\(O(n\log n)\))を採用します(簡略化/エラー回避のためクイックソート(\(O(n\log^2 n )\))を使用することも可能です)。

長さ \(w\) の部分文字列のランク \(rk_w[1..n]\) が既知であると仮定します(つまり、\(rk_w[i]\) は \(s[i..min(i+w−1,n)]\) が $ {s[x..min(x+w−1,n)] x \in [1,n] } $ の中で何番目に小さいかを表します)、次に、\(rk_w[i]\) を第一キー、\(rk_w[i+w]\) を第二キー(\(i+w>n\) の場合は \(rk_w[i+w]\) を無限小に設定)としてソートを行うことで、\(rk_{2w}[1..n]\) を求めることができます。

基数ソートについては、各キーを優先度の高い順にバケツソートを行い、時間計算量は値域に依存します。これにより、\(O(n\log n)\) のクイックソートを \(O(n)\) に最適化します。

アルゴリズムの可視化を参考に理解を深めることができます。

元のコード

scanf("%s", s + 1);
  n = strlen(s + 1);
  m = max(n, 300);
  //mは値域のサイズを表します
  for (i = 1; i <= n; ++i) ++cnt[rk[i] = s[i]];
  for (i = 1; i <= m; ++i) cnt[i] += cnt[i - 1];
  for (i = n; i >= 1; --i) sa[cnt[rk[i]]--] = i;
  for (w = 1; w < n; w <<= 1) {
    //rk[i+w]で初めてバケツソートし、結果をidに保存
    memset(cnt, 0, sizeof(cnt));
    for (i = 1; i <= n; ++i) id[i] = sa[i];
    for (i = 1; i <= n; ++i) ++cnt[rk[id[i] + w]];
    for (i = 1; i <= m; ++i) cnt[i] += cnt[i - 1];
    for (i = n; i >= 1; --i) sa[cnt[rk[id[i] + w]]--] = id[i];
    //rk[i]で二回目のバケツソートし、結果をsaに保存
    memset(cnt, 0, sizeof(cnt));
    for (i = 1; i <= n; ++i) id[i] = sa[i];
    for (i = 1; i <= n; ++i) ++cnt[rk[id[i]]];
    for (i = 1; i <= m; ++i) cnt[i] += cnt[i - 1];
    for (i = n; i >= 1; --i) sa[cnt[rk[id[i]]]--] = id[i];
    /*
    上記の二回のバケツソートは以下の簡潔なクイックソートに変換できます
    bool cmp(int x,int y){
		return rk[x] == rk[y] ? rk[x + w] < rk[y + w] : rk[x] < rk[y];
	}
    sort(sa + 1, sa + n + 1, cmp)	
    */
    memcpy(oldrk, rk, sizeof(rk));
    //saの値に基づいてrkの値を更新(以前のランクが同じものを更新)
    for (p = 0, i = 1; i <= n; ++i) {
      if (oldrk[sa[i]] == oldrk[sa[i - 1]] &&
          oldrk[sa[i] + w] == oldrk[sa[i - 1] + w]) {
        rk[sa[i]] = p;
      } else {
        rk[sa[i]] = ++p;
      }
    }
  }

最適化:

\(rk[i+w]\) のソートは不要です。なぜなら、すでに \(sa[1\sim n]\) がわかっているからです。\(i+w\) が \(n\) を超える場合は無限小として処理し、残りの第 \(i\) 位は \(sa[i]-p\) とすればよいです。

値域の範囲を更新できます。最後の部分で現在の値域の範囲 \(p\) を求めているので、それを使って \(m\) を更新します。

値域がすでに \(n\) である、つまり \(rk\) がすべて異なる場合は、終了できます。

その他、いくつかの特殊な最適化テクニック(連続したメモリ領域にアクセスすると高速化されるというもの)があります。

とにかく、最終的なコードは以下のようになります。

  int i, m = 300, p, w;
  scanf("%s", s + 1);
  n = strlen(s + 1);
  for (i = 1; i <= n; ++i) ++cnt[rk[i] = s[i]];
  for (i = 1; i <= m; ++i) cnt[i] += cnt[i - 1];
  for (i = n; i >= 1; --i) sa[cnt[rk[i]]--] = i;
  for (w = 1;; w <<= 1, m = p) {  // 最適化2
    //最適化1
    for (p = 0, i = n; i > n - w; --i) id[++p] = i;
    for (i = 1; i <= n; ++i)
      if (sa[i] > w) id[++p] = sa[i] - w;
    memset(cnt, 0, sizeof(cnt));
    //最適化4
    for (i = 1; i <= n; ++i) ++cnt[px[i] = rk[id[i]]];
    for (i = 1; i <= m; ++i) cnt[i] += cnt[i - 1];
    for (i = n; i >= 1; --i) sa[cnt[px[i]]--] = id[i];
    memcpy(oldrk, rk, sizeof(rk));
    for (p = 0, i = 1; i <= n; ++i)
      //最適化4
      rk[sa[i]] = cmp(sa[i], sa[i - 1], w) ? p : ++p;
    //最適化3
    if (p == n) {
      for (int i = 1; i <= n; ++i) sa[rk[i]] = i;
      break;
    }
  }


応用

アルゴリズムの説明はここまでですが、一つ問題に気づきました。接尾辞配列の用途が少し少ないのではないでしょうか...接尾辞の辞書順を処理するだけで、他には使えないようにも見えます...

実は、接尾辞配列で本当に役立つ重要な配列がここで登場します。

\(height\)配列

定義:\(height[i]=lcp(sa[i],sa[i−1])\)、つまり第$ i $位の接尾辞とその前の位の接尾辞の最長共通接頭辞です。

この定義の利点や、何ができるか考えてみましょう。

これまでに苦労してすべての接尾辞をソートしましたが、下図のようになります。

例として、ランク6の接尾辞を考えてみましょう(すべての\(height\)が求まっていると仮定します)。この接尾辞が他のすべての接尾辞と接頭辞がいくつ同じか知りたい場合、どうすればよいでしょうか?

まず、ランク5の接尾辞とは\(height[6]=2\)個の部分が同じです。ランク4の接尾辞とは1個の部分が同じ、ランク4のもので1個、ランク3のもので0個、ランク2のもので0個同じです。ランク\(i\)の接頭辞と\(j\)の接頭辞の共通部分は\(min_{j\le k\le i}\{height[k]\}\)であることがわかります(\(j> i\)の場合も同様です)。なぜでしょうか?

辞書順の定義を考えてみてください。\(s_1\)と\(s_2\)の辞書順が小さいということは、\(s_1[1\sim i]\)と\(s_2[1\sim i]\)が同じで、\(s_1[i+1]\)と\(s_2[i+1]\)が異なることを意味します。もし\(height\)が一定の値以上であれば、その長さの接頭辞はずっと同じままであることがわかります。もし特定の値\(k\)より小さくなった場合、\(s[k]\)以降は再び一致することはありません。

次に、これら\(n\)個の接尾辞が持つ性質を考えてみましょう。

各\(height\)を求めた後、実際には\(n\)個の接尾辞の等価情報を圧縮しています。そして、\(n\)個の接尾辞の各接頭辞は、すべての部分文字列を完全に表現できるため、部分文字列の等価問題を解く上で大きな助けとなります!

これが\(height\)配列が接尾辞配列で応用される実質的な意味です。次に、\(height\)配列をどのように求めるかを考えてみましょう(ついに求め方の説明を始めます...)。

現在\(height[rk[i]]\)を取得したい場合(\(height[i]\)は\(lcp(sa[i],sa[i−1])\)を表します)、実際には\(height[rk[i-1]]\)を利用できます。ここで補題\(height[rk[i]]≥height[rk[i−1]]−1\)を示します。

式が理解できなくても心配は不要です。図を参考に想像してみてください(下図の青い線と赤い線は接尾辞\(i-1\)と接尾辞\(sa[rk[i-1]-1]\)が等しい部分を表します)。

\(height[rk[i]]\)にはすでに「最低保証」があることがわかります。つまり、黒い下線と緑の下線の部分は必ず等しいですが、\(sa[rk[i]-1]\)が\(x+1\)であるとは限りません。

実際には問題ありません。辞書順を考慮すると、文字列\(s_i\)と\(s_x\)(\(s_x<s_i\))があり、すでに\(s_i\)と\(s_x\)の最初の\(y\)文字が等しい場合、\(s_x< s_j<s_i\)であり、かつ\(s_j\)と\(s_i\)の等しい接頭辞の長さが\(y\)より小さいような\(s_j\)は存在しないことがわかります。

したがって、\(height[rk[i]]\)を求める際には、すでに\(height[rk[i-1]]-1\)の「最低保証」があるため、接尾辞\(sa[rk[i]-1]\)と接尾辞\(i\)の接頭辞の等しい長さはそれより長くなる可能性がありますが、最低保証の長さより小さくなることはありません。

\(height\)を求めるコードは、\(sa\)を求めた後に続けます。

for (i = 1, k = 0; i <= n; ++i) {
  if (k) --k;
  while (s[i + k] == s[sa[rk[i] - 1] + k]) ++k;
  ht[rk[i]] = k;  // heightは長いためhtと省略
}

これが接尾辞配列の基本的な流れです(それほど理解しやすいでしょう...)。

次に、いくつかの古典的な問題を紹介します。

。。。。。。。

タグ: 接尾辞配列 文字列アルゴリズム 基数ソート データ構造

5月24日 06:50 投稿