SA-IS 法のメモ

suffix array を線形時間で構築する SA-IS 法についてのメモです。

SA-IS 法の解説はわりと世の中にいっぱいありますが、実際のプログラムにする方法がよくわからなったり、どうしてそれでうまく行くのか書かれてなくて気持ち悪かったりするものが多くて、自分の望む感じのものがありませんでした。アルゴリズムは当然プログラムで見たいし、厳密な証明は要らないけど直観的な理由は知りたい。

ということで、自分なりの理解を書いてみました。

suffix array とは

文字列の suffix とは、文字列の途中から最後までの部分文字列のことです。題材の文字列を "mmiissiissiippii" とすると、次の 17 個の部分文字列です。

 0 mmiissiissiippii$
 1 miissiissiippii$
 2 iissiissiippii$
 3 issiissiippii$
 4 ssiissiippii$
 5 siissiippii$
 6 iissiippii$
 7 issiippii$
 8 ssiippii$
 9 siippii$
10 iippii$
11 ippii$
12 ppii$
13 pii$
14 ii$
15 i$
16 $

最後の $ は番兵です。

これを普通に辞書順にソートします。文字には普通に全順序があり、番兵は他のどの文字よりも小さいとします。

16 $
15 i$
14 ii$
10 iippii$
 6 iissiippii$
 2 iissiissiippii$
11 ippii$
 7 issiippii$
 3 issiissiippii$
 1 miissiissiippii$
 0 mmiissiissiippii$
13 pii$
12 ppii$
 9 siippii$
 5 siissiippii$
 8 ssiippii$
 4 ssiissiippii$

このときのインデックスの列、[16, 15, 14, 10, 6, 2, 11, 7, 3, 1, 0, 13, 12, 9, 5, 8, 4] を suffix array といいます。今回説明しませんが、suffix array は検索や圧縮などでいろいろと役立つ数列です。

suffix array を求めるプログラムを Ruby でナイーブに書くと、これだけです。

s = "mmiissiissiippii".bytes + [0]

p (0...s.size).sort_by {|i| s[i..-1] }
#=> [16, 15, 14, 10, 6, 2, 11, 7, 3, 1, 0, 13, 12, 9, 5, 8, 4]

このプログラムは、ソートに O(n log n) 、しかも各比較に O(n) もかかるので、全体で O(n^2 log n) もかかります。SA-IS はこれを O(n) でやってしまいます。すごい。

SA-IS 法の超概要

流れとしては、次のような感じのアルゴリズムです。

1. 適当な「種」を元に、induced sort というのをやる → 間違った suffix array が得られる
2. 間違った suffix array を使って、正しい「種」を得る
3. 正しい「種」を元に、induced sort をもう一度やる → 正しい suffix array が得られる

「種」が何かはあとで説明するとして、最初の目標は induced sort というものを理解することです。しかし induced sort を理解するには、「L 型と S 型」「LMS」「ビン」という概念をまず理解する必要があります。順に説明します。

L 型と S 型

文字列 s のインデックス i から始まる suffix を s[i..] と書くことにします。s[0..] = "mmiissiissiippii$" です。

各インデックスを、次のルールで L 型と S 型に分けます。

• 最後（番兵）のインデックスは S åž‹
• s[i..] > s[i+1..] だったらインデックス i は L 型（i は i+1 より Larger）
• s[i..] < s[i+1..] だったらインデックス i は S 型（i は i+1 より Smaller）

「インデックス i が L 型」と言ったり、「s[i..] が L 型」と言ったりします。なお、すべての suffix は長さが異なるので、s[i..] == s[i+1..] になることはありません。

例題の文字列を分類すると次の通り。

mmiissiissiippii$
LLSSLLSSLLSSLLLLS

これを計算するには、文字列を逆順に走査して決めていけばいいです。基本的に s[i] と s[i+1] の文字を比べて決めていきます。

# 題材の文字列
s = "mmiissiissiippii".bytes + [0]
k = 256 # 文字種の数

# L 型か S 型か
t = [nil] * s.size

# 最後は S
t[-1] = :S

(s.size - 2).downto(0) do |i|
  # s[i] < s[i+1] なら明らかに s[i..] < s[i+1..] => i は S 型
  # s[i] > s[i+1] なら明らかに s[i..] > s[i+1..] => i は L 型
  # s[i] == s[i+1] の場合、s[i..] <=> s[i+1..] の比較結果は
  # s[i+1..] <=> s[i+2..] に準ずる (つまり t[i + 1] と同じ)
  case
  when s[i] < s[i + 1] then t[i] = :S
  when s[i] > s[i + 1] then t[i] = :L
  else                      t[i] = t[i + 1]
  end
end

LMS と LMS 部分文字列

インデックス i - 1 が L 型、i が S 型になっているとき、i を LMS（Left-most S, 最も左の S）と言います。

例題で言えば、2 番目、6 番目、10 番目、16 番目が LMS 。印をつけた位置。

          1111111
01234567890123456
mmiissiissiippii$
LLSSLLSSLLSSLLLLS
  ^   ^   ^     ^ <= LMS のインデックス

ついでに、LMS から次の LMS までの部分文字列を LMS 部分文字列と言います。これはだいぶ後で出てきます。

mmiissiissiippii$
LLSSLLSSLLSSLLLLS
  ^   ^   ^     ^
  iissi           <= LMS 部分文字列
      iissi       <= LMS 部分文字列
          iippii$ <= LMS 部分文字列
                $ <= LMS 部分文字列

ビン

まず、文字列 s と同じ長さの配列 sa を用意します。ここに suffix array の数列を入れるのが目標です。

# 作業領域
sa = [nil] * s.size

冒頭の suffix array を見ると、i で始まる suffix は 1 番目から 8 番目、m で始まる suffix は 9 番目から 10 番目、というように、suffix の先頭の文字ごとに範囲が決まっています。この範囲を、文字の「ビン」といいます。たとえば、「文字 i のビン」は 1 番目から 8 番目です。

このビンは、文字列 s の中の各文字の出現回数を数えればわかります。

# s に出現する文字種ごとのカウント
bin = [0] * (k + 1)
s.each {|ch| bin[ch + 1] += 1 }

# 文字種を累積することで bin のインデックスを決める
sum = 0
bin = bin.map {|v| sum += v }

これによって、文字 ch で始まる suffix は、sa の bin[ch] 番目から bin[ch+1] - 1 番目のどこかに入れればよい、ということになります。

それから、L 型の suffix と S 型の suffix が同じビンに入る場合（つまり先頭の文字が同じ場合）、必ず L 型の方を前に入れることになります。理由は→*1

ということで、配列 sa を次のように表現することにします。

 0 $ S : --
 1 i L : --
 2 i L : --
 3 i S : --
 4 i S : --
 5 i S : --
 6 i S : --
 7 i S : --
 8 i S : --
 9 m L : --
10 m L : --
11 p L : --
12 p L : --
13 s L : --
14 s L : --
15 s L : --
16 s L : --

左端の数字は sa 内のインデックス、その右の文字はビン、さらにその右の文字は L 型か S 型かを表してます。たとえば、i で始まる S 型の suffix は 3 番目から 8 番目のどこかに入ります。induced sort では、このルールを満たすように suffix のインデックスを入れていきます。

induced sort

いよいよ最初の目標であった induced sort の説明に入ります。induced sort は 3 つの段階からなります。

1. LMS のインデックスをビンの終わりから詰めていく
2. L 型のインデックスをビンの頭から詰めていく
3. S 型のインデックスを（LMS を含めて再度）ビンの終わりから詰めていく

とりあえずアルゴリズムを説明して、それからそれがだいたい何をやっているか説明します。

induced sort: ステップ 1

ステップ 1 は、LMS のインデックスを詰めていきます。LMS の 2 、6 、10 、16 を、それぞれ i 、i 、i 、$ のビンに終わりから入れます。

 0 $ S : 16 $
 1 i L : --
 2 i L : --
 3 i S : --
 4 i S : --
 5 i S : --
 6 i S :  2 iissiissiippii$
 7 i S :  6 iissiippii$
 8 i S : 10 iippii$
 9 m L : --
10 m L : --
11 p L : --
12 p L : --
13 p L : --
14 p L : --
15 s L : --
16 s L : --

右端に、インデックスから始まる suffix を参考に書いています。

ステップ 1 を Ruby で書くとこんな感じ。

# インデックス i が LMS かどうか
def lms?(t, i)
  i > 0 && t[i - 1] == :L && t[i] == :S
end

# LMS のインデックスだけ取り出したリスト（「種」）
lmss = (0...s.size).select {|i| lms?(t, i) }

# ステップ 1: LMS のインデックスをビンの終わりの方から入れる

count = [0] * k # ビンごとにすでに挿入した数をカウントする

# 挿入する順序は適当
lmss.reverse_each do |i|
  ch = s[i]
  # ch を入れるビンの終わり (bin[ch + 1] - 1) から詰めていれる
  sa[bin[ch + 1] - 1 - count[ch]] = i
  count[ch] += 1
end

実はこのとき、LMS をどういう順序で挿入するかが、超概要で言っていた「種」です。正しい順で LMS を挿入すれば、induced sort によって正しい suffix array が計算できてしまいます。しかし、この時点で正しい順序はわからないので、適当に 2 、6 、10 と並べました。正しくは上から 10 、6 、2 の順にならないといけないので、不正解です。最終的な結果も間違ったものになります。しかし不思議なことに、その間違った結果から、正しい「種」を抽出できるのです。詳しくは後で。

induced sort: ステップ 2

ステップ 2 は、sa を正順に走査して L 型の suffix を埋めていきます。

まず、sa に入っている最初のインデックスは 16 です。その 1 つ前のインデックス 15 は L 型なので、ステップ 2 で扱う対象です。s[15] は i なので i のビンに入れます。このとき、ビンの頭に詰めます。

 0 $ S : 16 $    <===== 今見ているところ
 1 i L : 15 i$   <===== 挿入位置
 2 i L : --
 3 i S : --
 4 i S : --
 5 i S : --
 6 i S :  2 iissiissiippii$
 7 i S :  6 iissiippii$
 8 i S : 10 iippii$
 9 m L : --
10 m L : --
11 p L : --
12 p L : --
13 p L : --
14 p L : --
15 s L : --
16 s L : --

このとき、「挿入位置」は必ず「今見ているところ」より後に来ます。理由→*2

次は、入れたばかりの 15 です。14 も L 型なので入れます。これを繰り返していくと、最終的に次のようになります。

 0 $ S : 16 $
 1 i L : 15 i$
 2 i L : 14 ii$
 3 i S : --
 4 i S : --
 5 i S : --
 6 i S :  2 iissiissiippii$
 7 i S :  6 iissiippii$
 8 i S : 10 iippii$
 9 m L :  1 miissiissiippii$
10 m L :  0 mmiissiissiippii$
11 p L : 13 pii$
12 p L : 12 ppii$
13 s L :  5 siissiippii$
14 s L :  9 siippii$
15 s L :  4 ssiissiippii$
16 s L :  8 ssiippii$

ステップ 2 を Ruby で書くとこんな感じ。

# ステップ 2: sa を昇順に走査していく

count = [0] * k # ビンごとにすでに挿入した数をカウントする

sa.each do |i|
  next if i == nil
  next if i == 0
  next if t[i - 1] == :S

  # sa に入っているインデックス i について、i - 1 が L 型であるとき、
  # 文字 s[i - 1] に対応するビンに i - 1 を頭から詰めていれる
  ch = s[i - 1]
  sa[bin[ch] + count[ch]] = i - 1
  count[ch] += 1
end

induced sort: ステップ 3

ステップ 3 は、sa を逆順に捜査して S 型の suffix を埋めていきます。LMS も S 型の一種なので埋めなおします。

まず、sa の最後に入っているインデックスは 8 です。その 1 つ前のインデックス 7 は S 型なので埋めます。s[7] は i なので i のビンに入れます。このとき、ビンの終わりから詰めていれます。

 0 $ S : 16 $
 1 i L : 15 i$
 2 i L : 14 ii$
 3 i S : --
 4 i S : --
 5 i S : --
 6 i S :  2 iissiissiippii$
 7 i S :  6 iissiippii$
 8 i S :  7 issiippii$   <===== 挿入位置
 9 m L :  1 miissiissiippii$
10 m L :  0 mmiissiissiippii$
11 p L : 13 pii$
12 p L : 12 ppii$
13 s L :  5 siissiippii$
14 s L :  9 siippii$
15 s L :  4 ssiissiippii$
16 s L :  8 ssiippii$    <===== 今見ているところ

ステップ 2 と同じように、「挿入位置」は必ず「今見ているところ」より前に来ます。理由も同様です。

なお、もともと入っていた 10 を書きつぶしていることに注意。最初から入っていた LMS は（"$" を除いて）いったん全部書きつぶされ、その後で再挿入されます。理由は次の節で。

これを繰り返すと、最終的に次のようになります。

 0 $ S : 16 $
 1 i L : 15 i$
 2 i L : 14 ii$
 3 i S : 10 iippii$
 4 i S :  2 iissiissiippii$
 5 i S :  6 iissiippii$
 6 i S : 11 ippii$
 7 i S :  3 issiissiippii$
 8 i S :  7 issiippii$
 9 m L :  1 miissiissiippii$
10 m L :  0 mmiissiissiippii$
11 p L : 13 pii$
12 p L : 12 ppii$
13 s L :  5 siissiippii$
14 s L :  9 siippii$
15 s L :  4 ssiissiippii$
16 s L :  8 ssiippii$

確かに LMS の 2 と 6 と 10 が（元と異なる位置に）挿入されています。

こうして得られた sa は、ほぼソートされているように見えますが、一部間違っています。例えば、s[4..] = "ssiissiippii$" が 15 番目、s[8..] = "ssiippii$" が 16 番目に来ていますが、この順序は逆ですね。これは後で直します。

Ruby はこちら。

# ステップ 3: sa を逆順に走査していく

count = [0] * k # ビンごとにすでに挿入した数をカウントする

sa.reverse_each do |i|
  next if i == nil
  next if i == 0
  next if t[i - 1] == :L

  # sa に入っているインデックス i について、i - 1 が S 型であるとき、
  # 文字 s[i - 1] に対応するビンに i - 1 を終わりから詰めていれる
  ch = s[i - 1]
  sa[bin[ch + 1] - 1 - count[ch]] = i - 1 # 上書きすることもある
  count[ch] += 1
end

induced sort の結果の考察

induced sort をしたとき、その結果はどんな性質を満たしているでしょうか。

ステップ 1 はビンソートなので、少なくとも各 suffix の最初の文字についてはちゃんとソートできています。しかし、各 suffix の 2 文字目以降はソートされていないかもしれません。次の図はステップ 1 の結果で、ソートされていないかもしれないところを ... で省略して表示しています。

 0 $ S : 16 $
 1 i L : --
 2 i L : --
 3 i S : --
 4 i S : --
 5 i S : --
 6 i S :  2 i...
 7 i S :  6 i...
 8 i S : 10 i...
 9 m L : --
10 m L : --
11 p L : --
12 p L : --
13 s L : --
14 s L : --
15 s L : --
16 s L : --

ステップ 2 では、すでに入っている suffix より 1 つ長い suffix を入れていきます。このとき、先頭 n 文字についてソートされていた suffix を 1 つ伸ばした suffix は、先頭 n + 1 文字がソートされていることが保証されます。理由→*3

これを繰り返すと、ステップ 2 完了時点で次のようになります。

 0 $ S : 16 $
 1 i L : 15 i$
 2 i L : 14 ii$
 3 i S : --
 4 i S : --
 5 i S : --
 6 i S :  2 i...
 7 i S :  6 i...
 8 i S : 10 i...
 9 m L :  1 mi...
10 m L :  0 mmi...
11 p L : 13 pii$
12 p L : 12 ppii$
13 s L :  5 si...
14 s L :  9 si...
15 s L :  4 ssi...
16 s L :  8 ssi...

ステップ 3 はステップ 2 と全く同様です。

 0 $ S : 16 $
 1 i L : 15 i$
 2 i L : 14 ii$
 3 i S : 10 iippii$
 4 i S :  2 iissi...
 5 i S :  6 iissi...
 6 i S : 11 ippii$
 7 i S :  3 issi...
 8 i S :  7 issi...
 9 m L :  1 mi...
10 m L :  0 mmi...
11 p L : 13 pii$
12 p L : 12 ppii$
13 s L :  5 si...
14 s L :  9 si...
15 s L :  4 ssi...
16 s L :  8 ssi...

ということで、... になっていない部分についてはソートされていることが保証されます。これが induced sort によって言えること。

さて、もしもステップ 1 の段階で、LMS が正しい順序（「種」）で挿入されていたとします。ステップ 2 とステップ 3 の「先頭 n 文字についてソートされている suffix を元に、1 つ長い suffix を先頭 n + 1 文字についてソートされた状態で挿入する」という性質から、正しい「種」から始めて最終的に得られる sa は、すべての suffix の全体についてソートされていることになります。それはすなわち、正しい suffix array です。

ということで、どうにか正しい「種」を得るというのが残る課題です。実はこれは、上記の「間違った induced sort の結果」を使って取り出すことができるのですが、その前に正しい「種」とは何かを考察します。

正しい「種」

正しい「種」を取り出すナイーブな実装としては、LMS の suffix を列挙して、辞書順にソートすればいいです。

LMS の suffix を列挙したもの。

 2: iissiissiippii$
 6: iissiippii$
10: iippii$
16: $

↓ソート

16: $
10: iippii$
 6: iissiippii$
 2: iissiissiippii$

この [16, 10, 6, 2] が求める「種」です。induced sort のステップ 1 で、この「種」を逆順に走査して、各ビンの最後から詰めていくと次のようになります。

 0 $ S : 16 $
 1 i L : --
 2 i L : --
 3 i S : --
 4 i S : --
 5 i S : --
 6 i S : 10 iippii$
 7 i S :  6 iissiippii$
 8 i S :  2 iissiissiippii$
 9 m L : --
10 m L : --
11 p L : --
12 p L : --
13 p L : --
14 p L : --
15 s L : --
16 s L : --

この状態でステップ 2 と 3 を実行すれば、めでたく正しい suffix array が得られます。

ただ、ナイーブな実装では O(n) が達成できません。そこで、このソートに SA-IS 法を再帰的に使うというのが、SA-IS 法の肝です。

SA-IS 法の再帰

LMS の suffix の列を SA-IS 法でソートするために、suffix の文字単位を「粗く」するのがポイントです。

"mmiissiissiippii$" を LMS 部分文字列（LMS を導入したときについでに説明してました）で分割すると、["iissi", "iissi", "iippii$", "$"] になります *4 。ここで、"iissi" や "iippii$" を 1 つの「文字」とみなし、この列を「文字列」と考えます。そして、この「文字列」のすべての suffix を並べてみます。

0: ["iissi", "iissi", "iippii$", "$"] ( 2: iissiissiippii$)
1: ["iissi", "iippii$", "$"]          ( 6: iissiippii$)
2: ["iippii$", "$"]                   (10: iippii$)
3: ["$"]                              (16: $)

最後の "(2: iissiissiippii$)" や "(16: $)" は、各 suffix が元の文字列のどの suffix に対応するかを表しています。

各「文字」の順序関係（"$" < "iippii$" < "iissi"）に注意しつつ、この「文字列」の suffix の列をソートします。

3: ["$"]                              (16: $)
2: ["iippii$", "$"]                   (10: iippii$)
1: ["iissi", "iippii$", "$"]          ( 6: iissiippii$)
0: ["iissi", "iissi", "iippii$", "$"] ( 2: iissiissiippii$)

注目すべきは、元文字列のインデックスの列。めでたく、[16, 10, 6, 2] という正しい「種」が得られています。

ところで、今おこなったソートは、"iissi" などを「文字」とみなした「文字列」の suffix array を作ったのと同じです。よって SA-IS 法が再帰的に適用できます。再帰すると全体で O(n) にならなくなりそうですが、「文字列」の長さは元の文字列の長さに比べて半分未満になってるので、O(n + n/2 + n/4 + ...) = O(2n) = O(n) ということで、線形時間が保たれます。すごい。

ちなみに、このアルゴリズム構成技法を再帰減速（recursive slowdown）といいます。余談ですが、これを遅延評価でやった暗黙的再帰減速（implicit recursive slowdown）という技法もあります。こういうのが面白いなと思った人は、次の本がとても面白いかもしれません。

純粋関数型データ構造

posted with amazlet at 18.01.30

Chris Okasaki
KADOKAWA (2017-04-28)
売り上げランキング: 82,013

Amazon.co.jpで詳細を見る

「文字」に番号を振る

さて、"iissi" や "iippii$" を 1 つの「文字」として扱うと言いましたが、実際にこういう部分文字列を切り出すと、比較に O(n) かかってしまうのでダメです。そこで、次のように番号を振ることを考えます。

"$"       => 0
"iippii$" => 1
"iissi"   => 2

この番号付けは、"$" < "iippii$" < "iissi" という順序関係を維持しています。そして、番号同士なら比較が O(1) でできます。よって、番号で置き換えてから SA-IS 法を再帰すれば OK です。

ただ、この番号を振る方法は自明ではありません。下手をしたら番号を振るために O(n^2) とかかかってしまいます。

この番号付けのために、1 回目の induced sort の（間違った）結果が使えます。再掲。

 0 $ S : 16 $
 1 i L : 15 i$
 2 i L : 14 ii$
 3 i S : 10 iippii$
 4 i S :  2 iissi...
 5 i S :  6 iissi...
 6 i S : 11 ippii$
 7 i S :  3 issi...
 8 i S :  7 issi...
 9 m L :  1 mi...
10 m L :  0 mmi...
11 p L : 13 pii$
12 p L : 12 ppii$
13 s L :  5 si...
14 s L :  9 si...
15 s L :  4 ssi...
16 s L :  8 ssi...

ここから、LMS の suffix だけを取り出してみます。

 0 $ S : 16 $
 3 i S : 10 iippii$
 4 i S :  2 iissi...
 5 i S :  6 iissi...

... を除くと、LMS 部分文字列ばかり出てきました。これは偶然ではありません。理由→*5

よって、この情報を活用することで、さっきのような番号付けを実現できます。隣り合う LMS 部分文字列を比べて、異なる場合は別の番号を、同じ場合は同じ番号を順に振っていけば OK です。具体的には次のような感じ。

• "$" に番号 0 を振っておく
• "$" と "iippii$" は異なるので、"iippii$" に番号 1 を振る
• "iippii$" と "iissi" は異なるので、"iissi" に番号 2 を振る
• "iissi" と "iissi" は同じなので、新しい番号は振らない

これで所望の番号付けができました。ちなみに、このときは LMS 部分文字列を直接比較しますが、文字の比較回数の合計が O(n) なので、問題ありません。

実装

実際にほしいのは、単なる番号付けではなく、["iissi", "iissi", "iippii$", "$"] の各 LMS 部分文字列を番号で置き換えたもの、つまり [2, 2, 1, 0] です。これを一気に求めます。

# induced sort の結果から LMS の suffix だけ取り出す
sa = sa.select {|i| lms?(t, i) }

# LMS のインデックス i に対して番号 nums[i] を振る
nums = []

# sa[0] の LMS は $ と決まっているので、番号 0 を振っておく
nums[sa[0]] = num = 0

# 隣り合う LMS を比べる
sa.each_cons(2) do |i, j|
  diff = false

  # 隣り合う LMS 部分文字列の比較
  s.size.times do |d|
    if s[i + d] != s[j + d] || lms?(t, i + d) != lms?(t, j + d)
      # LMS 部分文字列の範囲で異なる文字があった
      diff = true
      break
    elsif d > 0 && (lms?(t, i + d) || lms?(t, j + d))
      # LMS 部分文字列の終端に至った
      break
    end
  end

  # 隣り合う LMS 部分文字列が異なる場合は、番号を増やす
  num += 1 if diff
  
  # LMS のインデックス j に番号 num を振る
  nums[j] = num
end

# nums の中に出現する番号のみを並べると、LMS 部分文字列を番号に置き換えた列が得られる
nums = nums.compact

sa.each_cons(2) の中で sa.size.times を使っているので、一瞬 O(n^2) の二重ループのようにも見えますが、よく考えるとループの実行回数は最大でも元の文字列全体を 1 回走査するのと同じなので、O(n) に収まります。

さて、これで得られた nums に SA-IS 法を再帰適用します。ただ、もし隣り合う LMS 部分文字列が全部異なるものだった場合、つまり番号の重複がない場合は、いちいち再帰するまでもなく、ビンソート（各ビンのサイズは 1）の考え方で suffix array を直接求めることができます。

if num + 1 < nums.size
  # 番号の重複があるので再帰
  sa = sa_is(nums, num + 1)
else
  # 番号の重複がない場合、suffix array を容易に求められる
  sa = []
  nums.each_with_index {|ch, i| sa[ch] = i }
end

そして、これで得られる sa は上記で言う [3, 2, 1, 0] に相当する数列なので、これを LMS インデックスの列 [16, 10, 6, 2] に変換します。

lmss = (0...s.size).select {|i| lms?(t, i) }
seed = sa.map {|i| lmss[i] }

そして、この「種」を元に induced sort を再度行えば、ついに完了です。長かった。

まとめ

suffix array を O(n) で作るアルゴリズム SA-IS 法を説明しました。

世にある解説が自分にはいまいちわかりにくいものしか見つからなかったので書いてみましたが、やっぱりこれも他の人にとってはわかりにくいものになっているのかもしれません。コメントをくれたら改良を考えます。

おまけ：プログラム全体

# インデックス i が LMS かどうか
def lms?(t, i)
  i > 0 && t[i - 1] == :L && t[i] == :S
end

def induced_sort(s, k, t, lmss)
  # 作業領域
  sa = [nil] * s.size

  # s に出現する文字種ごとのカウント
  bin = [0] * (k + 1)
  s.each {|ch| bin[ch + 1] += 1 }

  # 文字種を累積することで bin のインデックスを決める
  sum = 0
  bin = bin.map {|v| sum += v }
  
  # ステップ 1: LMS のインデックスをビンの終わりの方から入れる

  count = [0] * k # ビンごとにすでに挿入した数をカウントする

  lmss.reverse_each do |i|
    ch = s[i]
    # ch を入れるビンの終わり (bin[ch + 1] - 1) から詰めていれる
    sa[bin[ch + 1] - 1 - count[ch]] = i
    count[ch] += 1
  end

  # ステップ 2: sa を昇順に走査していく

  count = [0] * k # ビンごとにすでに挿入した数をカウントする

  sa.each do |i|
    next if i == nil
    next if i == 0
    next if t[i - 1] == :S

    # sa に入っているインデックス i について、i - 1 が L 型であるとき、
    # 文字 s[i - 1] に対応するビンに i - 1 を頭から詰めていれる
    ch = s[i - 1]
    sa[bin[ch] + count[ch]] = i - 1
    count[ch] += 1
  end

  # ステップ 3: sa を逆順に走査していく

  count = [0] * k # ビンごとにすでに挿入した数をカウントする

  sa.reverse_each do |i|
    next if i == nil
    next if i == 0
    next if t[i - 1] == :L

    # sa に入っているインデックス i について、i - 1 が S 型であるとき、
    # 文字 s[i - 1] に対応するビンに i - 1 を終わりから詰めていれる
    ch = s[i - 1]
    sa[bin[ch + 1] - 1 - count[ch]] = i - 1 # 上書きすることもある
    count[ch] += 1
  end

  sa
end

def sa_is(s, k)
  # L 型か S 型かのテーブル
  t = [nil] * s.size

  # 最後は S
  t[-1] = :S

  (s.size - 2).downto(0) do |i|
    # s[i] < s[i+1] なら明らかに s[i..] < s[i+1..] => i は S 型
    # s[i] > s[i+1] なら明らかに s[i..] > s[i+1..] => i は L 型
    # s[i] == s[i+1] の場合、s[i..] <=> s[i+1..] の比較結果は
    # s[i+1..] <=> s[i+2..] に準ずる (つまり t[i + 1] と同じ)
    case
    when s[i] < s[i + 1] then t[i] = :S
    when s[i] > s[i + 1] then t[i] = :L
    else                      t[i] = t[i + 1]
    end
  end

  # LMS のインデックスだけを集めた配列
  lmss = (0...s.size).select {|i| lms?(t, i) }

  # 適当な「種」：seed = lmss.shuffle でもよい
  seed = lmss

  # 1 回目の induced sort
  sa = induced_sort(s, k, t, seed)

  # induced sort の結果から LMS の suffix だけ取り出す
  sa = sa.select {|i| lms?(t, i) }

  # LMS のインデックス i に対して番号 nums[i] を振る
  nums = []

  # sa[0] の LMS は $ と決まっているので、番号 0 を振っておく
  nums[sa[0]] = num = 0

  # 隣り合う LMS を比べる
  sa.each_cons(2) do |i, j|
    diff = false

    # 隣り合う LMS 部分文字列の比較
    s.size.times do |d|
      if s[i + d] != s[j + d] || lms?(t, i + d) != lms?(t, j + d)
        # LMS 部分文字列の範囲で異なる文字があった
        diff = true
        break
      elsif d > 0 && (lms?(t, i + d) || lms?(t, j + d))
        # LMS 部分文字列の終端に至った
        break
      end
    end

    # 隣り合う LMS 部分文字列が異なる場合は、番号を増やす
    num += 1 if diff
  
    # LMS のインデックス j に番号 num を振る
    nums[j] = num
  end

  # nums の中に出現する番号のみを並べると、LMS 部分文字列を番号に置き換えた列が得られる
  nums = nums.compact

  if num + 1 < nums.size
    # 番号の重複があるので再帰
    sa = sa_is(nums, num + 1)
  else
    # 番号の重複がない場合、suffix array を容易に求められる
    sa = []
    nums.each_with_index {|ch, i| sa[ch] = i }
  end

  # 正しい「種」
  seed = sa.map {|i| lmss[i] }

  # 2 回目の induced sort
  sa = induced_sort(s, k, t, seed)

  sa
end

# 題材の文字列
s = "mmiissiissiippii".bytes + [0]
k = 256 # 文字種の数
p sa_is(s, k)

なお、このプログラムは Ruby らしく大変富豪的に書かれています（配列作りまくり）。元論文の C プログラムは、ビンの位置以外の配列をすべて SA という 1 つの配列の再利用で解決していてかっこいいです。

出典

• Ge Nong, Sen Zhang, and Wai Hong Chan. Two Efficient Algorithms for Linear Time Suffix Array Construction

追記（2018-02-10）：コメントでの指摘を受けてバグ修正。「隣り合う LMS 部分文字列の比較」のループ終了条件が間違ってました。

*1:L 型は、1 文字目が 2 文字目より Larger です（例："ba...")。S 型は、1 文字目が 2 文字目より Smaller です（例："bc..."）。明らかに "ba..." < "bc..." なので、確かに L 型は S 型より前に来てます。1 文字目と 2 文字目が同じになっている文字列の場合もまあ同じように証明できます。

*2:インデックス i に対して i-1 が L 型のときだけ挿入する。L 型の定義から s[i-1..] > s[i..] 。つまり s[i-1..] の挿入位置は s[i..] より後になる。

*3:インデックス 10 の suffix "i..." は先頭 1 文字についてソートされています。今回は、インデックス 10 を元に、インデックス 9 の suffix "si..." が sa[15] の位置に挿入されました。これに注目してみます。もしも "si..." を s[15] に入れることが間違いだとしたら、sa[15] は文字 s のビン内なので、"si..." より大きいか、または小さい文字列が入るべきだったことになります。もし仮に、s[15] の位置には "si..." より大きい文字列（たとえば "sz..."）が入るべきだったとすると、s[13] から s[14] はすでに埋まっていることから、"si..." はいったいどこに行けばいいのかわからない、ということになるので矛盾。s[15] の位置に "si..." より小さい文字列（たとえば "sa..."）が入るべきだったとすると、a の文字のビンの中に "a..." のような文字列があったはずです。そしてステップ 2 では上から順番に処理しているので、その文字列が先に処理され、sa[15] にはすでに "sa..." が入っていたはずです。しかし実際には sa[15] にそういう文字列が入っていなかったので、矛盾。ということで、インデックス 5 の suffix "si..." が入るのが正しいということになります。

*4:先頭の "mm" は LMS 部分文字列ではないので捨てます。また、LMS 部分文字列の切れ目で 1 文字重複してます。

*5:特定の suffix に注目して induced sort の動きを見てみます。ステップ 1 で sa[8] (sa の 8 番目) に入れられたインデックス 10 に注目してみましょう。ステップ 2 では、インデックス 10 を元に、L 型インデックス 9 が s[14] に入れられました。そのインデックス 9 を元に、L 型インデックス 8 が s[16] に入れられました。次のインデックス 7 は S 型なので、ステップ 2 はここで終わりです。ステップ 3 では、インデックス 8 を元に、S 型インデックス 7 が s[8] に入れられました。インデックス 7 を元に、S 型インデックス 6 が s[5] に入れられました。次のインデックス 5 は L 型なので、ステップ 3 はここで終わり。要するに、LMS から始めてその前の L 型を順次追加し、さらにその前の S 型を順次追加する、スタートの LMS からみて 1 つ前の LMS にたどり着いた時点で終了、という動きになっています。LMS から 1 つ前の LMS までの範囲（すなわち LMS 部分文字列）についてソートされる、ということになります。