SA-IS: SuffixArray線形構築

『Linear Suffix Array Construction by Almost Pure Induced-Sorting』*1という論文を参考にして、Induced-Sortingを用いたSuffixArrayの線形構築アルゴリズム(SA-IS)Common Lispで実装した。
以下、そのソースコードとメモ書き。

線形構築

汎用ソートアルゴリズムと今回実装したSA-ISアルゴリズムでのSuffixArray構築時間の簡単な比較。
百文字〜一千万文字の文字数のテキストを入力として与えて、その構築時間を計測。

百文字千文字一万文字十万文字百万文字一千万文字
汎用ソート*20.000 sec0.002 sec0.023 sec0.358 sec5.604 sec94.196 sec
SA-IS0.031 sec0.120 sec0.135 sec0.335 sec2.580 sec29.616 sec
SA-ISの方は今回の実装が雑なため小さい入力に対しては効率が悪くなっているが、文字数が十万文字以上の場合の処理時間は、入力サイズに対してほぼ線形となっている。
対して、組み込みの汎用ソート関数(おそらくクイックソート)を用いた場合は、入力サイズが大きくなるにつれて処理時間が、だいたい(クリックソートの処理オーダーである)NlogNに比例する形で伸びていっている。

注意書き(2010/12/21)

以降の内容には(明確な)間違いが含まれていたので、以下の訂正記事も参照のこと。※ 基底ケースでのSuffixArrayの求め方が間違っていた
SA-IS: SuffixArray線形構築: 訂正 - sileの日記

アルゴリズム概略(メモ)

アルゴリズムの概略。
自分の理解の整理メモ。
※ 間違っている可能性があるので、正しく知りたい人はオリジナルの論文を参照のこと
※ 仮に間違っていないとしても、かなり自分流の解釈 & 不十分な記述となっているので、正しく知りたい人はオリジナルの論文を参照のこと


SA-ISアルゴリズムの大まかな流れ。

  • 1) 入力文字列をLMS部分文字列単位に分割する ※ LMS部分文字列が何かは後述
  • 2) LMS部分文字列単位でソート(induced-sort)を行う
  • 3) 全てのLMS部分文字列が互いにユニークかどうかをチェックする
  • 4-a) ユニークなら
    • 5-a) そこからSuffixArrayを求める(induced-sort)
  • 4-b) ユニークでない(重複する部分文字列がある)なら、
    • 5-b) 以下の方法によって新たな入力文字列を生成する。
      • ソート済みの各LMS部分文字列に対して、その出現順に昇順に番号を割り振る
      • その際に、互いに等しいLMS部分文字列には同じ番号を振る
      • 入力文字列の各LMS部分文字列を、上で割り振った番号で置換する
      • => もともとの入力文字列がLMS部分文字列のソート順番号によって簡約化される
    • 6-b) その文字列に対して再帰的にSA-ISを適用する
    • 7-b) その返り値からSuffixArrayを求める(induced-sort)

もっと単純化してしまえば、「入力文字列からLMS部分文字列の列を求める」と「それに対してinduced-sortを適用する」の二つの処理の繰り返し、となる。
最初に求めたLMS部分文字列列の全ての要素が互いにユニークなら、一度のinduced-sortでSuffixArrayが計算可能。
重複がある場合は、二つ(もしくはそれ以上)のLMS部分文字列間に優劣をつける必要があるので、SA-ISを再帰的に適用する。
※ 返ってきた SuffixArray == 重複する要素が正しくソートされたLMS部分文字列列 なので、それに対してもう一度induced-sortを行えば、もともとの入力に対する適切なSuffixArrayが得られる。


LMS部分文字列列を求めるのとそれに対するinduced-sortは、どちらも入力文字列に対して線形時間で行える。
かつ、SA-ISの再帰的な呼び出しの際に、新たな入力文字列は、もともとの入力文字列の半分以下のサイズになることが保証されているので、全体の処理時間も最初の入力テキストに対して線形となる。
※ 新たな入力文字列のサイズ == もともの入力文字列のLMS部分文字列列のサイズ


LMS部分文字列:
LMS部分文字列が何か。
論文中に出てくる"mmiissiissiippii$"という文字列を例とする。
※ 末尾の"$"はヌル文字を表す。SA-ISを自然に記述するには、入力文字列の末尾に番兵値(入力文字列中で最も小さい文字)が必要。


入力文字列から、LMS部分文字列を求めるには、まず以下のルールに従って、各文字をSLかのどちらかのタイプに分類する必要がある。

  • 1) 文字列[i]の文字は、文字列[i+1]の文字より小さいならSとなる。 ※ smallのS ?
  • 2) 文字列[i]の文字は、文字列[i+1]の文字より大きいならLとなる。 ※ largeのL ?
  • 3) 文字列[i]の文字が、文字列[i+1]の文字と等しい場合は、文字列[i+1]のタイプが、文字列[i]のタイプとなる。
  • 4) 最後の文字(= ヌル文字)は、常にSとなる。

以下はこのルールに従って分類した結果。

"mmiissiissiippii$"
 LLSSLLSSLLSSLLLLS

LMSは"LeftMost S-character"の略(確か…)で、「一番左側」つまりLタイプの文字の直後にあるSタイプの文字、はLMS文字となる。
例で云えば、以下で'@'が付いている文字がLMS文字。

"mmiissiissiippii$"
 LLSSLLSSLLSSLLLLS
   @   @   @     @

そして、二つのLMS文字に挟まれた区間(ただし両端のLMS文字自身も含まれる)がLMS部分文字列となる。
※ 例外として、一番最後のヌル文字は、それ単独でLMS部分文字列。
なので、入力の"mmiissiissiippii$"は、以下の四つのLMS部分文字列に分割される。

SSLLS   "iissi"
SSLLS   "iissi"    ; "iissi"は二度目 (= ユニークではない)
SSLLLLS "ippii$"
S       "$"

LMS部分文字列の特徴としては、以下のような点が挙げられる。

  • LMS部分文字列のサイズは必ず2以上となる。※ 末尾のヌル文字を除く
    • このため、SA-ISの再帰的適用の際の入力は、もともとの入力のサイズの半分以下となる。
  • LMS部分文字列は全て「1個以上のS」+「1個以上のL」+「一個のS」という形式になっている。※ 末尾のヌル文字を除く
    • つまり、末尾の文字があって、その前方により大きい文字が続き、さらにその前方により小さい文字が続く、という形なっている。
    • この性質がinduced-sortの際には重要。多分。
  • 連接するLMS部分文字列の先頭と末尾は共有される。
    • 現在のLMS部分文字列の末尾文字 == 次のLMS部分文字列の先頭文字
    • つまり、LMS部分文字列の先頭文字の前方には、「1個以上のS」+「1個以上のL」、という形式の文字列が続くことになる。
      • ※ ただし、一番初めのLMS部分文字列だけは例外で、「0個以上のS」+「1個以上のL」、となる。

次に記述するinduced-sort及びSA-ISの再帰的な適用は、このLMS部分文字列列を対象に行われる。


induced-sort:
induced-sort。
バケツの用意と、三回の走査。


バケツの用意:
induced-sortに際して、まずは入力文字列の各文字用の箱(バケツ)を用意する。

入力: "mmiissiissiippii$"

;; 入力の各文字用のバケツを用意する。
;; - バケツは辞書順(文字の値に昇順)に並べる。
;; - バケツに入れる際に、
;;  -- Lタイプの文字(で始まるsuffix)は、先頭から挿入する  ※ まず先頭に挿入、次は二番目に挿入、...
;;  -- Sタイプの文字(で始まるsuffix)は、末尾から挿入する ※ まず末尾に挿入、次は末尾から二番目に挿入、...
;; - '_'は要素未挿入を表す初期値のつもり
;; - バケツには入力文字列の各suffixの開始位置を格納する
"$": L[_]S
"i": L[_,_,_,_,_,_,_,_]S
"m": L[_,_]S
"p": L[_,_]S
"s": L[_,_,_,_]S

コメントにある通り、各文字用のバケツは辞書順に並んでいる。
つまり、各バケツに入力文字列中の文字を挿入していくだけで、各suffixの最初の文字レベルでは、ソートが行えるということになる。
※ もちろん、これだけでは先頭文字が等しいsuffix同士の順番が正しくソートされるとは限らないので、駄目だけど。


三回の走査とバケツへの挿入:
induced-sortの残り。
三ステップからなる。

  • 1) 各LMS部分文字列の先頭文字(の位置)をバケツに挿入する。
    • 先頭文字は全てSタイプなので、各文字に対応するバケツの末尾から挿入していく。
    • 先頭文字が同じとなるLMS部分文字列同士の挿入順序は任意で構わない。
      • もし、この挿入順序が適切なら、結果として得られるバケツの中身は、入力文字列に対するSuffixArrayを表すようになる
      • '適切'とは、先頭文字から始まるsuffix同士がソート順になるようにバケツに格納されている、ということ
      • SA-IS再帰的適用後は、LMS部分文字列が'適切'に並んだ列が得られるので、その順序を保持したまま、この1を行えば、最終的に得られるバケツの中身 == 入力文字列のSuffixArray、となる
  • 2) 各バケツを順に、先頭要素から走査し、もし走査中の要素が示す文字の一つ前の文字のタイプがLなら、その文字(の位置)をバケツに挿入する。
    • つまり、入力文字列[バケツ[n][i]-1].type == L、なら、バケツ[n][i]-1、の文字をバケツに挿入する。
    • この段階で挿入される文字のタイプは当然全てL
    • かつ、タイプがLの全ての文字は、この段階でバケツに挿入されることになる。 ※ 後述
  • 3) 各バケツを逆順に、末尾要素から走査し、もし走査中の要素が示す文字の一つ前の文字のタイプがSなら、その文字(の位置)をバケツに挿入する。
    • つまり、入力文字列[バケツ[n][i]-1].type == S、なら、バケツ[n][i]-1、の文字をバケツに挿入する。
    • ただし、この走査前に1の段階で挿入したSタイプの文字(に対応する位置)をバケツから除外しておく。
      • ※ ヌル文字以外
    • この段階で挿入される文字のタイプは当然全てS
    • かつ、タイプがSの全ての文字は、この段階でバケツに挿入されることになる。 ※ 後述

つまり、まず各LMS部分文字列の先頭文字(の位置)をバケツに挿入し、次にその前方に続くLタイプの文字(の位置)を全てバケツに挿入し、最後にさらにその前方に続くSタイプの文字(の位置)をバケツに挿入する、ということを行っている。
※ ただし、最後のSタイプ文字挿入ステップの前に、その時点でバケツに格納されているSタイプの文字(= 最初に挿入したLMS部分文字列の先頭文字)は、除去しておく必要がある


各LMS部分文字列が互いにユニークだとした場合、最終的に得られるバケツの中身が、入力文字列に対するSuffixArrayとなる。


何故この三ステップでSuffixArrayが求まるか、というと、だいたい以降のような感じになると思う。
まず、一ステップの段階で、各LMSの先頭文字だけを見れば、辞書順に並んでいることになる。※ 各バケツが辞書順に並んでいるため

;; 一ステップを終えた段階でのバケツの中身
;; - バケツ内の数字は、各文字の入力文字列内での位置  ※ 各suffixの開始位置でもある
;; - 各数字の下に並んでいる文字列は、その位置から次のLMS文字までの部分文字列
;;
;; - 先頭文字だけを見れば、辞書順に並んでいる。
"$"[16] "i"[__,__,__,__,__,02,06,10] "m"[__,__] "p"[__,__] "s"[__,__,__,__]
     $                      i  i  i
                            i  i  i
                            s  s  p
                            s  s  p
                            i  i  i
                                  i
                                  $

;; 入力文字列とLMS
00 01 02 03 04 05 06 07 08 09 10 11 12 13 14 15 16
 m  m  i  i  s  s  i  i  s  s  i  i  p  p  i  i  $
 L  L  S  S  L  L  S  S  L  L  S  S  L  L  L  L  S
       @           @           @                 @


次に第二ステップ。
バケツ中の要素を先頭から走査し、入力文字列[バケツ[n][i]-1].type == L、ならバケツ[n][i]-1にある文字(の位置)をバケツに挿入する。
第一ステップでバケツに格納される文字('@'が付いてる文字)は、全てLMS部分文字列の先頭文字で、その前方にはLタイプの文字が続いている。
また、Lタイプの文字は「一つ後ろの文字と等しいか、それよりも大きい文字」であるため、バケツを先頭から走査しつつ、要素の挿入を行っていけば、自然と全てのLタイプの文字(の位置)がバケツに格納されることになる。

で、その要素挿入では、以下の二つの性質が常に満たされるので、結果としてLタイプの文字(の位置)はソート順でバケツに挿入されることになる。

  • a) 値が小さい文字は、より前方のバケツに挿入される
    • これはバケツの並び順自体の性質
    • また、バケツの並び順およびLタイプの文字の定義上、現在走査中の位置より前方に、新たに要素が挿入されることはない
      • つまり、新たに挿入された要素は、必ず一度走査される (ので、連鎖的に全てのLタイプの文字が挿入・走査される)
  • b) 後ろに続く文字(文字列/suffix)がより小さい文字は、同じバケツ内のより前方に挿入される
    • これはバケツの走査順により保証される
      • 「バケツを先頭から走査する」ということは、「より小さいsuffixから順に走査する」ということに等しい
      • そのため、「より小さいsuffix」を接尾に持つLタイプの文字は、より先に走査され、より先にバケツ内に挿入されることになる
      • 新たに挿入される文字は、常にLタイプなので、それを追加したことにより、バケツ内の順序が狂うことはない
        • 言い換えれば、新たに挿入される文字は常に「現在走査中の文字より、等しいかより大きい文字」であり、それが表すsuffixは常に「現在走査中のsuffixより少しだけ大きいsuffix」となる(ので、より後方に挿入されるのは適切)

以下は、第二ステップの挿入過程を表した図。

;; 入力文字列とLMS
00 01 02 03 04 05 06 07 08 09 10 11 12 13 14 15 16
 m  m  i  i  s  s  i  i  s  s  i  i  p  p  i  i  $
 L  L  S  S  L  L  S  S  L  L  S  S  L  L  L  L  S
       @           @           @                 @

;;; 第二ステップでの挿入過程
;;; ※ '#'は現在走査中の要素を、'!'は新たに挿入された要素を表す

;; 1) input[16-1].type == L なので 15 が挿入される
     #       !
"$"[16] "i"[15,__,__,__,__,02,06,10] "m"[__,__] "p"[__,__] "s"[__,__,__,__]
     $       i              i  i  i
             $              i  i  i
                            s  s  p
                            s  s  p
                            i  i  i
                                  i
                                  $

;; 2) input[15-1].type == L なので 14 が挿入される
             #  !
"$"[16] "i"[15,14,__,__,__,02,06,10] "m"[__,__] "p"[__,__] "s"[__,__,__,__]
     $       i  i           i  i  i
             $  i           i  i  i
                $           s  s  p
                            s  s  p
                            i  i  i
                                  i
                                  $

;; 3)
                #                                    !
"$"[16] "i"[15,14,__,__,__,02,06,10] "m"[__,__] "p"[13,__] "s"[__,__,__,__]
     $       i  i           i  i  i                  p
             $  i           i  i  i                  i
                $           s  s  p                  i
                            s  s  p                  $
                            i  i  i
                                  i
                                  $

;; 4)
                            #             !
"$"[16] "i"[15,14,__,__,__,02,06,10] "m"[01,__] "p"[13,__] "s"[__,__,__,__]
     $       i  i           i  i  i       m          p
             $  i           i  i  i       i          i
                $           s  s  p                  i
                            s  s  p                  $
                            i  i  i
                                  i
                                  $

;; 5)
                               #                                !
"$"[16] "i"[15,14,__,__,__,02,06,10] "m"[01,__] "p"[13,__] "s"[05,__,__,__]
     $       i  i           i  i  i       m          p          s
             $  i           i  i  i       i          i          i
                $           s  s  p                  i
                            s  s  p                  $
                            i  i  i
                                  i
                                  $
;; 9)
                                  #                                !
"$"[16] "i"[15,14,__,__,__,02,06,10] "m"[01,__] "p"[13,__] "s"[05,09,__,__]
     $       i  i           i  i  i       m          p          s  s
             $  i           i  i  i       i          i          i  i
                $           s  s  p                  i
                            s  s  p                  $
                            i  i  i
                                  i
                                  $


;; 10)
                                          #  !
"$"[16] "i"[15,14,__,__,__,02,06,10] "m"[01,00] "p"[13,__] "s"[05,09,__,__]
     $       i  i           i  i  i       m  m       p          s  s
             $  i           i  i  i       i  m       i          i  i
                $           s  s  p          i       i
                            s  s  p                  $
                            i  i  i
                                  i
                                  $


;; 11) input[0-1]の文字はないので、挿入なし
                                             #
"$"[16] "i"[15,14,__,__,__,02,06,10] "m"[01,00] "p"[13,__] "s"[05,09,__,__]
     $       i  i           i  i  i       m  m       p          s  s
             $  i           i  i  i       i  m       i          i  i
                $           s  s  p          i       i
                            s  s  p                  $
                            i  i  i
                                  i
                                  $

;; 12)
                                                     #  !
"$"[16] "i"[15,14,__,__,__,02,06,10] "m"[01,00] "p"[13,12] "s"[05,09,__,__]
     $       i  i           i  i  i       m  m       p  p       s  s
             $  i           i  i  i       i  m       i  p       i  i
                $           s  s  p          i       i  i
                            s  s  p                  $  i
                            i  i  i                     $
                                  i
                                  $


;; 13) input[12-1].type == S なので、挿入なし
                                                        #
"$"[16] "i"[15,14,__,__,__,02,06,10] "m"[01,00] "p"[13,12] "s"[05,09,__,__]
     $       i  i           i  i  i       m  m       p  p       s  s
             $  i           i  i  i       i  m       i  p       i  i
                $           s  s  p          i       i  i
                            s  s  p                  $  i
                            i  i  i                     $
                                  i
                                  $

;; 14)
                                                                #     !
"$"[16] "i"[15,14,__,__,__,02,06,10] "m"[01,00] "p"[13,12] "s"[05,09,04,__]
     $       i  i           i  i  i       m  m       p  p       s  s  s
             $  i           i  i  i       i  m       i  p       i  i  s
                $           s  s  p          i       i  i             i
                            s  s  p                  $  i
                            i  i  i                     $
                                  i
                                  $

;; 15)
                                                                   #     !
"$"[16] "i"[15,14,__,__,__,02,06,10] "m"[01,00] "p"[13,12] "s"[05,09,04,08]
     $       i  i           i  i  i       m  m       p  p       s  s  s  s
             $  i           i  i  i       i  m       i  p       i  i  s  s
                $           s  s  p          i       i  i             i  i
                            s  s  p                  $  i
                            i  i  i                     $
                                  i
                                  $


最後は第三ステップ。
これは、バケツの走査順その他を第二ステップとは逆にしただけなので、説明は省略。


ここまでは、入力文字列から求めた各LMS部分文字列が互いにユニークだとした場合の話で、それらの間に重複があると、少し話が複雑になる。
※ 例となっている"mmiissiissiippii$"には、LMS部分文字列"iissi"が二つ含まれており、重複がある
上に書いた三ステップでは、見ての通り(?)、バケツへの挿入時に、入力文字列を各LMS部分文字列単位でしか扱っていない。
そのため、仮に二つのLMS部分文字列が全く等しいとすると、それら(が表すsuffix)の間の大小関係を決定する方法は存在せず、結果として得られる順序は、第一ステップでのバケツへの挿入順*3に依存することになる。
その問題を解消するために用いるのが、次に書くSA-ISの再帰的な適用。


SA-IS再帰的適用、およびその際に用いる入力文字列:
入力文字列から求めた各LMS部分文字列間に重複がある場合の解消法。
SA-ISの再帰的適用。


既に書いた通り、問題となっているのは、互いに重複するLMS部分文字列同士の位置関係。
そのどちらが前に来て、どちらが後ろに来るべきか、を知る必要がある。
これは、その重複する部分文字列同士だけを見ていては判断不可能なので、後続する部分文字列同士の大小関係を考慮しなければならない。
入力文字列の各文字をSタイプとLタイプに割り振った時と同様に、部分文字列AとBが等しい場合は、その互いの後続の部分文字列CとDを比較して、CがDよりも小さいなら、AをBよりも小さいと判断し、その逆ならBがAより小さいと判断することになる。


各LMS部分文字列の大小関係自体は、induced-sortを行った段階で得られているので、それが利用できる。
また、ここで必要なのはその大小関係だけであり、各部分文字列の具体的な中身は不要なので、もともとの入力文字列は、その各LMS部分文字列をinduced-sort後の出現位置で置換したものに、簡約化することが可能。
※ 最初のLMS部分文字列はソート後はA番目に位置して、次のLMS部分文字列はB番目に位置して、...、ということ情報があれば良い
※ そのため、簡約化後の文字列 == LMS部分文字列列のサイズ、となる


そして、上で書いたような「両者が等しい場合は、それぞれの後続の要素を見て、その大小関係を判断する」といったようなことは、つまりはSuffixArrayを求める(ために行っている)ことと等しいので、簡約化された入力文字列に対してSA-ISアルゴリズム再帰的に適用することで、互いに重複するLMS部分文字列の位置関係を決定することが可能、となる。


以下、その具体的な手順等。

;;; a) 入力文字列とLMS
00 01 02 03 04 05 06 07 08 09 10 11 12 13 14 15 16
 m  m  i  i  s  s  i  i  s  s  i  i  p  p  i  i  $
 L  L  S  S  L  L  S  S  L  L  S  S  L  L  L  L  S
       @           @           @                 @

;;; b) 第三ステップ終了後のバケツの中身
;;; '@'は、各LMS部分文字列の先頭位置を表す
"$"[16] "i"[15,14,10,02,06,11,03,07] "m"[01,00] "p"[13,12] "s"[05,09,04,08]
     @             @  @  @

;;; 新しい入力文字列の作り方とSA-ISの再帰的適用
;;;
;; 1) もともとの入力文字列からLMS文字の位置のみを抜き出す
=> [02, 06, 10, 16]

;; 2) 1で得た各要素を、バケツ内での出現順序で置換する
;;    ※ ただし、互いに等しいLMS部分文字列の先頭位置には、同じ番号を割り振る
;;
;;    - ソート後のバケツを見ると、[16,10,02,06]、の順番でLMS部分文字列の先頭が並んでいることが分かる
;;    - また、02と06から始まるLMS部分文字列は互いに等しいので、同じ番号を割り振る必要がある
;;     -- そのため16には0番を、10には1番を、02と06には2番を割り振る
;;
;; 結果として、以下の列が得られる。
=> [2, 2, 1, 0]   ; = LMS部分文字列列["iissi":2番目, "iissi":2番目, "ippii$":1番目, "$":0番目]を、そのソート位置で簡約化した列

;; 3) 上の列([2, 2, 1, 0])を、新たな入力文字列として、SA-ISを適用する
;;    ※ 実際には"文字列"というよりは数値列だが、論文に倣い"文字列"と呼称する
;;    
=> SA-IS([2,2,1,0]) -> [3, 2, 1, 0] 

;; 4) 結果として得られたSuffixArrayは、もともとのLMS部分文字列同士の正しい位置関係(ソート順)を示しているので、
;;    その順番を反映する形で、第一ステップで、バケツにLMS部分文字列の先頭文字を格納すれば、
;;    最終に得られるバケツの中身 == もともとの入力文字列に対するSuffixArray となる。
;;
;;    ※ 実際には、今回重複があったのは02と06で始まるLMS部分文字列だけだったので、
;;       その二つの順番のみが意味を持つ。
;;       上で返って来た列を見ると、02は四番目に、06は三番目に、位置するようになっているので、
;;       第一ステップの挿入時に、06が02よりも前方になるようにすれば良い。
;;
;;    以下、SA-IS再帰的適用の返り値を反映した、induced-sortの第一ステップ終了後のバケツの中身
"$"[16] "i"[__,__,__,__,__,10,06,02] "m"[__,__] "p"[__,__] "s"[__,__,__,__]
     $                      i  i  i
                            i  i  i 
                            p  s  s 
                            p  s  s 
                            i  i  i 
                            i       
                            $      


以上。
細かいところはいろいろ省略(or 簡略化)しているけど、それは冒頭の論文(or ソースコード)を参照、ということで。

実装

残りは実装ソースコード
最初の方にも書いた通り、雑なコードで最適化も全然されていない*4が、一応正常に動作はする。
※ 汎用ソート関数を用いたSuffixArray構築関数も末尾に掲載

;; suffix(末尾文字列)用構造体
(defstruct suffix
  (str "" :type string)   ; ソース文字列
  (pos 0  :type fixnum))  ; suffixの先頭位置

;; suffix構造体用の表示関数: suffixの十文字目までを表示する
(defmethod print-object ((o suffix) stream)
  (with-slots (str pos) o
    (format stream "#S~S" (subseq str pos (min (+ pos 10) (length str))))))


;;;;;;;;;;;;
;; SA-IS関数
(defun sais (str)
  (let ((str (make-sais-string str)))          ; 入力文字列をSA-IS用に変換する
    (conv-suffix-array str (sais-impl str))))  

;; 入力文字列をSA-ISで扱いやすいように変換する
;; 以下の二つの処理を施す
;; - ヌル終端にする
;; - character型から、整数型に変換する
(defun make-sais-string (src &aux (len (length src)))
  (let ((new (make-array (1+ len) :element-type 'fixnum :initial-element 0)))
    (dotimes (i len new)
      ;; XXX: src文字列が#\Nullを含んでいることもありえるので、全てに1を加算すべき
      (setf (aref new i) (char-code (char src i))))))

;; sais-impl関数の戻り値(= SuffixArrayの各suffixの先頭位置の配列)を、
;; suffix構造体の配列に変換する
(defun conv-suffix-array (sais-str suffix-starts)
  (let ((str (map 'string #'code-char sais-str)))
    (map 'vector (lambda (p) (make-suffix :str str :pos p)) suffix-starts)))


;;;;;;;;;;;;;;;;;;;;;;;;
;; 実際にSA-ISを行う関数
(defun sais-impl (str)
  (let* ((types (classify str))                ; 入力文字列をSとLに分類する
         (lms-blocks (calc-lms-blocks types))) ; LMS部分文字列列を求める
    (multiple-value-bind (s1 sa1 uniq?) (induced-sort str types lms-blocks)  ; induced-sort
      (if uniq?
          ;; 全てのLMS部分文字列がユニークなら、induced-sort関数の返り値をそのまま返す
          sa1
        ;; 重複する部分文字列がある場合は、簡約化された入力(s1)に対してsais-implを再帰的に適用し、
        ;; その結果に対して、もう一度induced-sortを行う。 ※ ただし、やることが若干ことなるので、上とは関数は別
        (induce-sa str types (sais-impl s1) lms-blocks)))))

;; 入力文字列の各文字をSタイプとLタイプに分類する
;; タイプはビット配列で表現し0ならSタイプ、1ならLタイプ、となる
(defun classify (str &aux (len (length str)))
  (let ((types (make-array len :element-type 'bit)))
    (setf (bit types (1- len)) 0)         ; 末尾のヌル文字は常にSタイプ
    (loop FOR i FROM (- len 2) DOWNTO 0
          FOR cur  = (aref str i)
          FOR next = (aref str (1+ i))
          FOR next-suffix = (bit types (1+ i))
      DO
      (setf (bit types i) (cond ((< cur next) 0)   ; 後続の文字より小さいならS
                                ((> cur next) 1)   ; 後続の文字より大きいならL
                                (t next-suffix)))) ; 等しい場合は、後続の文字のタイプが、自身のタイプとなる
    types))

;; posにある文字が、LMS文字かどうかを判定する関数
(declaim (inline lms?))
(defun lms? (types pos)
  (and (plusp pos)
       (> (bit types (1- pos)) (bit types pos))))  ; 自身がSタイプで、一つ前の文字がLタイプなら、LMS文字となる

;; LMS部分文字列を求める
;; 返り値は、各LMS部分文字列の開始位置を格納した配列
(defun calc-lms-blocks (types)
  (loop FOR i FROM 1 BELOW (length types)
        WHEN (lms? types i)
    COLLECT i INTO list
    FINALLY
    (return (coerce list 'vector))))


;;;;;;;;;;;;;;;;;;;;;
;; 各文字用のバケツを表す構造体
(defstruct (bucket-pos
            (:constructor make-bucket-pos (&key start last &aux (l-cur start) (s-cur last))))
  (start 0 :type fixnum)   ; バケツの開始位置
  (l-cur 0 :type fixnum)   ; Lタイプの要素を次に挿入する位置。startから始まる(増加)
  (s-cur 0 :type fixnum)   ; Sタイプの要素を次に挿入する位置。lastから始まる(減少)
  (last  0 :type fixnum))  ; バケツの最終位置

;; バケツ全体を表す構造体
(defstruct buckets
  (index #() :type simple-vector)  ; 文字からバケツ(bucket-pos)へのインデックス
  (elems #() :type simple-vector)) ; 全てのバケツの要素を実際に格納する配列

;; バケツの初期化
(defun init-buckets (s)
  (let ((limit (1+ (loop FOR x ACROSS s MAXIMIZE x))))  ; 入力文字列中の最大値の取得
    ;; 各文字の出現頻度をカウントする
    (loop WITH freq = (make-array limit :element-type 'fixnum :initial-element 0)
          FOR x ACROSS s
      DO
      (incf (aref freq x))

      FINALLY
      ;; 全てのバケツ用の格納領域と、文字からバケツを参照するためのインデックスを作成する
      (let* ((bkt (make-buckets :elems (make-array (length s) :initial-element -1)
                                :index (make-array limit)))
             (idx (buckets-index bkt)))
        (loop WITH offset OF-TYPE fixnum = 0
              FOR i FROM 0 BELOW limit
              FOR cnt = (aref freq i)
          DO
          ;; 文字i用のバケツを作成する
          (setf (aref idx i)
                (make-bucket-pos :start offset :last (1- (+ offset cnt))))
          (incf offset cnt))
        (return bkt)))))

;; バケツに要素を挿入する
(declaim (inline bucket-put))
(defun bucket-put (type char value buckets)
  (with-slots (index elems) buckets
    (with-slots (l-cur s-cur) (aref index char)  ; 文字charに対応するバケツを取得する
      (ecase type
        (:s (setf (aref elems s-cur) value)  ; Sタイプの文字なら、バケツの末尾側から順に文字を挿入
            (decf s-cur))
        (:l (setf (aref elems l-cur) value)  ; Lタイプの文字なら、バケツの先頭側から順に文字を挿入
            (incf l-cur))))))


;;;;;;;;;;;;;;;;
;; induce関数
(defun induce (first-char-pos-generator s ty bkt)
  ;; 1] ステップ一# LMS部分文字列の先頭文字をバケツに挿入する
  ;;    ※ 先頭文字の挿入順序によって、結果が変わってくるので、
  ;;       その具体的な順序は、first-char-pos-generator関数で制御する
  (loop FOR i = (funcall first-char-pos-generator)
        WHILE i
    DO
    (bucket-put :s (aref s i) i bkt))

  ;; 2] ステップ二# バケツを先頭から走査し、入力文字列中のLタイプの文字をバケツに挿入する
  (loop FOR i ACROSS (buckets-elems bkt)
        WHEN (> i 0)
    DO
    (when (= (aref ty (1- i)) 1)
      (bucket-put :l (aref s (1- i)) (1- i) bkt)))

  ;; 3-a] ステップ三の準備# ステップ一で挿入したSタイプの文字を上書き可能なように、s-curを初期化する
  (loop FOR pos ACROSS (buckets-index bkt)
    DO
    (with-slots (s-cur last) pos
      (setf s-cur last)))
 
  ;; 3-b] ステップ三# バケツを末尾から走査し、入力文字列中のSタイプの文字をバケツに挿入する
  (loop FOR x FROM (1- (length (buckets-elems bkt))) DOWNTO 0
        FOR i = (aref (buckets-elems bkt) x)
        WHEN (> i 0)
    DO
    (when (= (aref ty (1- i)) 0)
      (bucket-put :s (aref s (1- i)) (1- i) bkt)))
 
  (buckets-elems bkt))

;; SA-IS再帰呼び出し結果に対するinduced-sort
;;
;; s1はsais-impl再帰呼び出し結果のSuffixArray。
;; induced-sortのステップ一で、各LMS部分文字列先頭文字を
;; 挿入する順序を格納している。
(defun induce-sa (s ty s1 lms &aux (bkt (init-buckets s)) (i (length s1)))
  (induce (lambda ()
            (when (plusp i)
              (aref lms (aref s1 (decf i)))))
          s ty bkt))

;; 初めに適用されるinduced-sort関数
(defun induced-sort (s ty lms &aux (bkt (init-buckets s)))
  ;; まずはバケツに入力文字列の全ての文字(の位置)を挿入し、induced-sortを行う
  (let ((sa (let ((i (length lms)))
              (induce (lambda ()
                        ;; ステップ一での挿入順は任意: 以下では後方のLMS部分文字列(の先頭文字)から順に挿入している
                        (when (plusp i)
                          (aref lms (decf i))))
                      s ty bkt))))

    ;; 入力文字列の各LMS部分文字列を、その先頭文字のソート後の出現順序で簡約化した文字列(s1)を作成する
    ;; - そのために、まずはもともとの入力文字列と同じサイズの文字列を用意し、
    ;;    LMS部分文字列の先頭文字だけを、上記出現順序で置換する、ということを行う。
    ;;    ※ 先頭文字以外は-1にセットしておく
    (let ((s1 (make-array (length s) :element-type 'fixnum :initial-element -1))
          (order 0))
      (setf (aref s1 (aref sa 0)) order)  ; 末尾のヌル文字のみのLMS部分文字列の出現順序は常に0となる
      (loop WITH prev-pos = (aref sa 0)
            FOR i FROM 1 BELOW (length sa)
            FOR pos = (aref sa i)
            WHEN (lms? ty pos)
        DO
        (when (not (lmss= s ty prev-pos pos))
          ;; 現在のLMS部分文字列が既出でない場合は、出現順序をインクリメントする
          (incf order))
        (setf (aref s1 pos) order)  ; 出現順序をセット
        (setf prev-pos pos))

      ;; LMS部分文字列の先頭文字に対応する要素以外をs1から除外する
      (setf s1 (delete -1 s1))
     
      ;; (簡約化された入力文字列=s1 induced-sortの結果 s1の各要素がユニークかどうか) を返す
      (values s1 sa (= order (1- (length lms)))))))

;; 二つのLMS部分文字列が等しいかどうかを判定する
(declaim (inline lmss=))
(defun lmss= (s ty start1 start2)
  (flet ((== (i1 i2)
           (and (= (aref s i1) (aref s i2))     ; 二つの文字が等しい
                (= (bit ty i1) (bit ty i2)))))  ; かつ、二つの文字のタイプが等しい
    (declare (inline ==))

    (and (== start1 start2)
         (do ((i1 (1+ start1) (1+ i1))
              (i2 (1+ start2) (1+ i2)))
             ((and (lms? ty i1)
                   (lms? ty i2))
              t)
           (unless (== i1 i2)
             (return nil))))))


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; 汎用ソート関数によるSuffixArray構築
(defun sa (s &aux (s (format nil "~A~C" s #\Null)))
  (let ((ary (make-array (length s))))
    (dotimes (i (length s))
      (setf (aref ary i) (make-suffix :str s :pos i)))
    
    (sort ary (lambda (a b)
                (string< s s :start1 (suffix-pos a) :start2 (suffix-pos b))))))

*1:また『Two Efficient Algorithms for Linear Suffix Array Construction』という論文には、ほとんど同じ内容でアルゴリズムが説明されており、かつC++ソースコードも掲載されている。

*2:common-lisp:sort関数。おそらくクイックソート

*3:これは、素直な実装では、入力文字列中のLMS部分文字列の出現順によって決定される

*4:加えて未整理