接尾辞配列 (suffix array) の続きです。今回は Yuta Mori さん が考案された improved two-stage suffix sorting algorithm という方法を紹介します。このページでは略して ITS と呼ぶことにしましょう。ITS は二段階ソート法を改良したアルゴリズムです。二段階ソート法の改良には、前回と前々回で説明した三分割法 (ratio-1, ratio-2) がありますが、ソートするデータ数を比較すると、ITS の方が三分割法 ratio-2 よりも少なくなる場合が多いようです。そして、その分だけ高速に suffix array を構築することができます。
ITS のポイントはサフィックスのタイプ分けにあります。基本的な考え方は簡単で、三分割法にも適用することができます。そこで、まず最初に ITS のアイデアを拝借して、三分割法を改良してみましょう。そのあとで、Yuta Mori さんが考案された ITS を説明します。
三分割法は i 番目の記号 Si と i + 1 番目の記号 Si+1 の大小関係を使って、サフィックスを次に示すように Type A, B, C の 3 種類に分割する方法です。
Type A : Si > Si+1 Type B : Si = Si+1 Type C : Si < Si+1
Type C をソートする場合、三分割法 ratio-2 は Type C の区間を Type D (Si > Si+2) と Type E (Si <:= Si+2) に分けて、Type E だけをソートしました。ITS の基本的な考え方は、タイプ分けを行うところは三分割法 ratio-2 と同じですが、Si+1 と Si+2 を比較してタイプ分けを行うところが異なります。次の図を見てください。
┌─┬─┬─┬─┬─┬─┐
│ │ａ│ｂ│ｃ│ｄ│ｅ│Si+1
├─┼─┼─┼─┼─┼─┤ 1. Type C をソートする場合
│ａ│Ｂ│Ｃ│Ｃ│Ｃ│Ｃ│ Type C を二分割する
├─┼─┼─┼─┼─┼─┤
│ｂ│Ａ│Ｂ│Ｃ│Ｃ│Ｃ│ Type D : Si+1 >= Si+2 (ソートする)
├─┼─┼─┼─┼─┼─┤ Type E : Si+1 < Si+2
│ｃ│Ａ│Ａ│Ｂ│Ｃ│Ｃ│
├─┼─┼─┼─┼─┼─┤
│ｄ│Ａ│Ａ│Ａ│Ｂ│Ｃ│
├─┼─┼─┼─┼─┼─┤
│ｅ│Ａ│Ａ│Ａ│Ａ│Ｂ│
└─┴─┴─┴─┴─┴─┘
Si
図 : 三分割法の改良 (1)
記号は {a, b, c, d, e} の 5 種類とします。Type C をソートする場合、Type C を Si+1 と Si+2 の大小関係で二分割します。このページでは、Si+1 >= Si+2 を Type D とし、Si+1 < Si+2 を Type E とします。Type D をソートすると、その結果を使って Type E のサフィックスの順序を決定することができます。
サフィックスのタイプ分けは下図のように表すことができます。
┌─┬─────┬─────┬─────┬─────┬─────┐
│ │ ａ │ ｂ │ ｃ │ ｄ │ ｅ │Si+1
├─┼─────┼─┬───┼──┬──┼───┬─┼─────┤
│ａ│ Ｂ │Ｄ│ Ｅ │ Ｄ │ Ｅ │ Ｄ │Ｅ│ Ｄ │
├─┼─────┼─┴───┼──┼──┼───┼─┼─────┤
│ｂ│ Ａ │ Ｂ │ Ｄ │ Ｅ │ Ｄ │Ｅ│ Ｄ │
├─┼─────┼─────┼──┴──┼───┼─┼─────┤
│ｃ│ Ａ │ Ａ │ Ｂ │ Ｄ │Ｅ│ Ｄ │
├─┼─────┼─────┼─────┼───┴─┼─────┤
│ｄ│ Ａ │ Ａ │ Ａ │ Ｂ │ Ｄ │
├─┼─────┼─────┼─────┼─────┼─────┤
│ｅ│ Ａ │ Ａ │ Ａ │ Ａ │ Ｂ │
└─┴─────┴─────┴─────┴─────┴─────┘
Si
Type D をソート
Type E はセット
図 : 三分割法の改良 (2)
たとえば、記号 Si+1 が e の場合、e よりも大きな記号がないので、ee を除いたすべての区間 (ae, be, ce, de) は Type D になります。Si+1 が d の場合は、ade, bde, cde の区間が Type E で、それ以外の区間が Type D になります。同様に、Si+1 が b, c の場合も Type D, E に分けることができます。
次に、Type D の区間をソートして、Type D の区間を後ろから順番に調べていきます。まず最初に区間 de を調べます。この区間は Type D で全てソートされています。この区間のサフィックスの順序が決まると、サフィックス ade..., bde..., cde... の順序を決定することができます。つまり、区間 ad, bd, cd の Type E の区間にインデックスをセットすることができます。
これで、記号 c の Type C の区間 (cd, ce) がソートされます。次は、区間 ce, cd の結果を使って、区間 ac, bc の Type E にインデックスをセットします。これにより、記号 b の Type C の区間 (bc, bd, be) がソートされます。最後に、区間 bc, bd, be の結果を使って、区間 ab の Type E にインデックスをセットします。
このように、Si+1 と Si+2 の大小関係を使ってタイプ分けを行っても、Type D をソートするだけで Type E の区間のインデックスをセットすることができます。あとは、三分割法と同様に Type C の結果から Type A, B の区間のインデックスをセットすることができます。
それではプログラムを作りましょう。最初に先頭 2 記号で分布数えソートを行うのは今までと同じです。次に、Type C の区間を Type D, E に分離して、Type D の区間をソートします。プログラムは次のようになります。
リスト : Type C をソート
def sort_type_c():
# Type C を Type D, E を分離して Type D をソート
for x in xrange(0, 256):
for y in xrange(x + 1, 256):
c = (x << 8) + y
low = count_sum[c]
high = count_sum[c + 1]
m = low
for n in xrange(low, high):
i = idx[n]
if buff[i + 1] >= buff[i + 2]:
idx[n] = idx[m]
idx[m] = i
m += 1
if m - low > 1:
mqsort(low, m - 1, 2)
#
# Type E のセット
#
# 累積度数表をコピー
end = count_sum[:]
# 記号 0, 255 の区間はスキャンする必要なし
for x in xrange(254, 0, -1):
for y in xrange(255, x, -1):
c = (x << 8) + y
for n in xrange(count_sum[c + 1] - 1, count_sum[c] - 1, -1):
i = idx[n]
if i > 0 and buff[i - 1] < buff[i] < buff[i + 1]:
# Type E をセット
c1 = (buff[i - 1] << 8) + buff[i] + 1
end[c1] -= 1
idx[end[c1]] = i - 1
Type C をソートする場合、Type C の区間を Type D と Type E に分割してから Type D の区間をソートします。三分割法 ratio-2 とは違って、記号 buff[i + 1] と buff[i + 2] を比較して、データを Type D と Type E に分けます。このとき、Type D のサフィックスを区間の前半に集めることに注意してください。そして、Type D のデータをマルチキークイックソート (mqsort) でソートします。
次に、Type E の区間にインデックスをセットします。最初に累積度数表 count_sum をコピーして変数 end にセットします。count_sum は Type A, B のインデックスをセットする処理 (関数 set_type_ab) でも使うので、ここで値を書き換えてはいけません。インデックスは後ろからセットしていくので、大きな記号の区間から順番に調べていきます。記号 255 は Type A, B の区間しかないので調べる必要はありません。また、記号 0 の区間はインデックスがセットされるだけなので、この区間も調べる必要はありません。
あとは、count_sum から区間の上限値と下限値を求めて、後ろのサフィックスから順番に調べていきます。buff[i - 1] < buff[i] < buff[i + 1] であれば Type E のサフィックスなので、その区間の後ろからインデックス i - 1 をセットします。このとき、後ろからインデックスをセットしていくので、end[c1] を -1 することをお忘れなく。
Type A, B のインデックスをセットする関数 set_type_ab は三分割法と同じです。あとは特に難しいところはないでしょう。説明は割愛いたしますので、詳細は プログラムリスト１ をお読みください。
それでは、実際に Canterbury Corpus で配布されているテストデータ The Canterbury Corpus で suffix array を作成してみましょう。結果は suffix_sort の時間で、ファイルの入出力処理は含んでいません。結果は次にようになりました。
表 : suffix array の作成 (suffix_sort の時間)
ファイル名 サイズ ABC1 ABC2 ABC3
--------------------------------------------------------
alice29.txt 152,089 2.67 2.52 2.42
(67,459) (48,384) (47,620)
asyoulik.txt 125,179 2.19 2.02 1.93
(58,834) (40,898) (39,645)
cp.html 24,603 0.59 0.58 0.58
(11,181) (7,814) (7,213)
fields.c 11,150 0.30 0.32 0.38
(4,699) (3,156) (3,363)
grammar.lsp 3,721 0.16 0.22 0.25
(1,357) (1,078) (1,032)
kennedy.xls 1,029,744 16.54 13.06 18.27
(404,858) (201,689) (402,276)
lcet10.txt 426,754 8.41 7.46 6.86
(194,732) (137,506) (131,505)
plrabn12.txt 481,861 8.86 7.92 7.09
(224,466) (158,938) (143,021)
ptt5 513,216 4.66 5.17 4.43
(34,601) (29,788) (27,431)
sum 38,240 1.83 1.46 1.57
(13,175) (10,187) (9,245)
xargs.1 4,227 0.17 0.22 0.25
(1,981) (1,333) (1,285)
--------------------------------------------------------
合計 2,810,784 46.38 40.95 44.03
# 時間の単位 : 秒, () : ソートしたデータ数
実行環境 : Windows XP, celeron 1.40 GHz, Python 2.4.2
ABC1 は三分割法 ratio-1、ABC2 は三分割法 ratio-2、ABC3 が今回の改良版のプログラムです。ABC3 は ABC1 よりもソートするデータ数が少なくなり、その分だけ実行時間も速くなるファイルが多くなりました。テキストファイルの場合、ソートしたデータ数は ABC2 よりも少なくなり、ABC3 の方が速くなりました。バイナリファイルや小さなファイルでは、それほど大きな効果はないようです。それでも、記号 Si+1 と Si+2 でタイプを分ける方法 (ABC3) は、三分割法 ratio-2 と同様に大きな効果がある方法だと思います。
なお、実行時間の結果は M.Hiroi のコーディング、実行したマシン、プログラミング言語などの環境に大きく依存しています。また、これらの環境だけではなく、データの種類によっても実行時間はかなり左右されます。興味のある方は、いろいろなデータをご自分の環境で試してみてください。
次は improved two-stage suffix sorting algorithm (ITS) について説明します。ITS はサフィックスを次のように Type A, B, B* の 3 つに分けて、Type B* の区間だけをソートします。
Type A (1) Si > Si+1
(2) Si = Si+1 and Si+1 = Type A
Type B (1) Si < Si+1
(2) Si = Si+1 and Si+1 = Type B
Type B* (1) Si < Si+1 and Si+1 = Type A
ITS は Si のタイプを Si と Si+1 の大小関係と Si+1 のタイプで決定します。Si+1 のタイプを使って Si のタイプを分けるので、Si と Si+1 を比較するだけではタイプを決定することはできません。この場合、記号を後ろから比較していくと簡単にタイプを決めることができます。
三分割法と比較すると、タイプ分けが複雑になったように思われますが、基本的な考え方に大きな違いはありません。Type A, B の条件 (1) は三分割法の Type A, C と同じです。このあと (2) の条件で、Si+1 のタイプをチェックすることで、Si = Si+1 の区間を Type A, B に分割します。
そして、Type B* の条件は、Si+1 と Si+2 を比較してタイプを分ける処理に相当します。この処理が ITS のポイントです。Si+1 が Type A の場合、Si+1 >= Si+2 が成立します。ここで、Si+1 = Si+2 の場合を考えてみましょう。三分割法改良版は単純に Si+1 と Si+2 を比較するだけなので、Si+1 = Si+2 のデータは全てソートしますが、ITS では Type B* と B に分けることができます。次の図を見てください。
a b c c d e a b c c a b c $
-----------------------------------------
B B B B B A B B A A B B A /
* * *
図 : ITS のタイプ分け
最後尾のデータ c は Type A になります。ここで "bcc" の部分に注目してください。最初の "bcc" では、記号 b のタイプは B になりますが、次の "bcc" では B* になります。
cc が Type B の場合、その後ろには c が続くか、それよりも大きな記号になります。逆にcc が Type A の場合は、その後ろには c が続くか、それよりも小さな記号になります。したがって、同じ bcc でも cc が Type A であれば、cc が Type B のサフィックスよりも小さいことがわかります。つまり Si+1 = Si+2 の区間で、小さい方を Type B* に、大きい方を Type B に分けることができるのです。これを図に示すと次のようになります。
┌─┬─────┬─────┬─────┬─────┬─────┐
│ │ ａ │ ｂ │ ｃ │ ｄ │ ｅ │Si+1
├─┼─────┼─┬───┼──┬──┼───┬─┼─────┤
│ａ│ Ｂ │B*│ Ｂ │ B* │ Ｂ │ B* │Ｂ│ B* │
├─┼─────┼─┼───┼──┼──┼───┼─┼─────┤
│ｂ│ Ａ │Ａ│ Ｂ │ B* │ Ｂ │ B* │Ｂ│ B* │
├─┼─────┼─┴───┼──┼──┼───┼─┼─────┤
│ｃ│ Ａ │ Ａ │ Ａ │ Ｂ │ B* │Ｂ│ B* │
├─┼─────┼─────┼──┴──┼───┼─┼─────┤
│ｄ│ Ａ │ Ａ │ Ａ │ Ａ │Ｂ│ B* │
├─┼─────┼─────┼─────┼───┴─┼─────┤
│ｅ│ Ａ │ Ａ │ Ａ │ Ａ │ Ａ │
└─┴─────┴─────┴─────┴─────┴─────┘
Si
Type B* をソート
Type A, B はセット
図 : ITS のタイプ分け
たとえば、記号 Si+1 が e の場合、e よりも大きな記号がないので、ee を除いたすべての区間 (ae, be, ce, de) は Type B* になります。Si+1 が d の場合は、ade, bde, cde の区間と add, bdd, cdd の大きい部分が Type B で、それ以外の区間が Type B* になります。また、区間 dd は dde と ddd の大きい部分が Type B で残りの部分が Type A になります。同様に、Si+1 が b, c の場合も Type B*, B に分けることができます。区間 aa は Type B だけになりますが、最後尾の記号が a の場合は、区間の先頭が Type A になります。
次に、Type B* の区間をソートして、Type B* の区間を後ろから調べていきます。最初に、区間 de を調べます。この区間は Type B* で全てソートされています。この区間のサフィックスの順序が決まると、サフィックス ade..., bde..., cde..., dde... の順序を決定することができます。つまり、区間 ad, bd, cd の Type B の区間にインデックスをセットすることができます。
ただし、これだけでは *dd の大きな部分のインデックスがセットされません。この場合、区間 dd の Type B にセットされたインデックスを使って、*dd の大きい部分のインデックスをセットします。区間 de だけではなく、区間 dd の Type B の部分も調べることに注意してください。これで、記号 c の 区間 cd, ce がソート完了になります。
次は、区間 ce, cd の結果を使って、区間 ac, bc, cc の Type B にインデックスをセットします。これにより、記号 b の区間 bc, bd, be がソートされます。最後に、区間 bc, bd, be の結果を使って、区間 ab の Type B にインデックスをセットします。そして、その結果を使って残りの区間 aa にインデックスをセットします。
これで Type B にインデックスをセットすることができました。あとは Type A のインデックスをセットすればよいわけです。ITS は Si+1 = Si+2 の部分を Type B* と Type B に分けることができるので、三分割法改良版よりもソートするデータ数を少なくすることができます。
それではプログラムを作りましょう。ITS の場合、最初に分布数えソートを行う必要はありません。累積度数表 count_sum があれば、Type B* を分離してソートすることができます。次のリストを見てください。
リスト : Type B* のソート
def sort_type_b_star():
b_star_end = count_sum[:]
# 最後尾のデータは Type A
prev = 1
# Type B* の分離
for x in xrange(data_size - 2, -1, -1):
r = buff[x] - buff[x + 1]
if r == 0: continue
if r < 0 and prev > 0:
# Type B* をセット
c = (buff[x] << 8) + buff[x + 1]
idx[b_star_end[c]] = x
b_star_end[c] += 1
prev = r
# Type B* をソート
for x in xrange(0, 256):
for y in xrange(x + 1, 256):
c = (x << 8) + y
low = count_sum[c]
high = b_star_end[c]
if high - low > 1:
mqsort(low, high - 1, 2)
関数 sort_type_b_star は Type B* を分離してソートします。最初に count_sum をコピーして b_star_end にセットします。Type B* のインデックスを idx にセットするときは b_star_end を使います。そして、count_sum と b_star_end で Type B* の区間を表します。
次に Type B* を分離します。prev は一つ後ろの記号のタイプを表していて、正の整数が Type A で、負の整数が Type B になります。そして、buff を後ろから順番に調べていきます。buff[x] - buff[x + 1] を r にセットします。r が 0 の場合は、prev の値はそのままです。r が負の値でかつ prev が正の値であれば、そのインデックスは Type B* です。区間を変数 c にセットして、b_star_end[c] の位置にインデックスをセットします。そして、prev の値を更新します。
それから Type B* の区間をソートします。Type B* があるのは x < y の区間です。count_sum から下限値 low を、b_star_end から上限値 high を求めて、マルチキークイックソート (mqsort) でソートします。
次は Type A, B にインデックスをセットする関数 set_type_ab を作ります。次のリストを見てください。
リスト : Type A, B をセット
def set_type_ab():
end = [0] * 256
type_a_end = [0] * 256
# Type B をセット
for i in xrange(254, -1, -1):
for j in xrange(0, i + 1):
end[j] = count_sum[(j << 8) + i + 1] - 1
# 後ろから前へ
j = count_sum[(i + 1) << 8] - 1
while j > end[i]:
x = idx[j]
if x > 0 and buff[x - 1] <= buff[x]:
c = buff[x - 1]
idx[end[c]] = x - 1
end[c] -= 1
j -= 1
type_a_end[i] = end[i]
type_a_end[255] = count_sum[256 << 8] - 1
#
# Type A をセット
#
# 最後尾のデータをセット
if buff[data_size - 1] >= buff[data_size]:
c = (buff[data_size - 1] << 8) + buff[data_size]
idx[count_sum[c]] = data_size - 1
count_sum[c] += 1
# 前から Type A をセットする
for x in xrange(data_size):
i = idx[x]
if i > 0 and buff[i - 1] >= buff[i]:
c = (buff[i - 1] << 8) + buff[i]
if count_sum[c] <= type_a_end[buff[i - 1]]:
idx[count_sum[c]] = i - 1
count_sum[c] += 1
Type B は後ろからインデックスをセットするので、区間の最後尾の位置を配列 end にセットします。end の使い方は三分割法のプログラムと同じです。記号 i の区間を後ろから前へ探索しますが、探索する範囲は記号 i の区間の最後尾から Type B の区間の end[i] までになります。、
ここで、Type B の区間にインデックスがセットされると end[i] の値が減ることに注意してください。これで、探索の範囲が広がるわけです。そして、Type B にセットされたインデックスを利用して、他のサフィックスのインデックスをセットすることができます。Type B にインデックスをセットしたら、Type B の先頭位置 end[i] を type_a_end[i] にセットします。これは Type A のインデックスをセットするときに使います。
次に、Type A のインデックスをセットします。最初に、最後尾の記号のインデックスをセットします。次に、Type A のインデックスを探してセットします。この処理は二段階ソート法とほとんど同じですが、buff[i - 1] == buff[i] の区間は Type A と Type B に分かれているところが違います。これを区別するために type_a_end を使います。インデックスをセットする位置 count_sum[c] が type_a_end[buff[i - 1]] 以下であれば、その位置にインデックスをセットします。そうでなければ、その区間の Type A のインデックスはセット完了したので何もしません。
あとはとくに難しいところはないと思います。説明は割愛いたしますので、詳細は プログラムリスト２ をお読みください。
それでは、実際に Canterbury Corpus で配布されているテストデータ The Canterbury Corpus で suffix array を作成してみましょう。結果は suffix_sort の時間で、ファイルの入出力処理は含んでいません。結果は次にようになりました。
表 : suffix array の作成 (suffix_sort の時間)
ファイル名 サイズ ABC1 ABC2 ABC3 ITS
-------------------------------------------------------------------
alice29.txt 152,089 2.67 2.52 2.42 2.14
(67,459) (48,384) (47,620) (46,334)
asyoulik.txt 125,179 2.19 2.02 1.93 1.72
(58,834) (40,898) (39,645) (39,420)
cp.html 24,603 0.59 0.58 0.58 0.42
(11,181) (7,814) (7,213) (7,175)
fields.c 11,150 0.30 0.32 0.38 0.22
(4,699) (3,156) (3,363) (3,114)
grammar.lsp 3,721 0.16 0.22 0.25 0.11
(1,357) (1,078) (1,032) (891)
kennedy.xls 1,029,744 16.54 13.06 18.27 17.70
(404,858) (201,689) (402,276) (393,559)
lcet10.txt 426,754 8.41 7.46 6.86 6.42
(194,732) (137,506) (131,505) (129,139)
plrabn12.txt 481,861 8.86 7.92 7.09 6.68
(224,466) (158,938) (143,021) (142,114)
ptt5 513,216 4.66 5.17 4.43 4.86
(34,601) (29,788) (27,431) (27,010)
sum 38,240 1.83 1.46 1.57 1.40
(13,175) (10,187) (9,245) (9,134)
xargs.1 4,227 0.17 0.22 0.25 0.11
(1,981) (1,333) (1,285) (1,281)
-------------------------------------------------------------------
合計 2,810,784 46.38 40.95 44.03 41.78
# 時間の単位 : 秒, () : ソートしたデータ数
実行環境 : Windows XP, celeron 1.40 GHz, Python 2.4.2
ABC1 は三分割法 ratio-1、ABC2 は三分割法 ratio-2、ABC3 が今回の三分割法改良版、ITS が improved two-stage suffix sorting algorithm のプログラムです。ソートしたデータ数は kennedy.xls を除いた全てのファイルで ITS が一番少なくなりました。その結果、kennedy.xls と ptt5 以外のファイルで、ITS が一番速くなりました。ITS の効果はとても高いことがわかります。ITS は優秀なアルゴリズムだと思います。
ところで、ITS にランクソート法を適用することもできます。説明は割愛いたしますので、詳細は プログラムリスト３ をお読みください。結果は次のようになりました。
表 : suffix array の作成 (suffix_sort の時間)
ファイル名 サイズ LS ABC1R ABC2R ITSR
---------------------------------------------------
alice29.txt 152,089 4.10 2.73 2.56 2.28
asyoulik.txt 125,179 3.25 2.28 2.10 1.86
cp.html 24,603 0.56 0.48 0.51 0.38
fields.c 11,150 0.28 0.27 0.31 0.20
grammar.lsp 3,721 0.12 0.16 0.21 0.11
kennedy.xls 1,029,744 31.14 18.19 13.64 18.94
lcet10.txt 426,754 12.62 8.34 7.45 6.65
plrabn12.txt 481,861 14.52 9.43 8.46 7.35
ptt5 513,216 24.33 4.19 4.68 4.73
sum 38,240 0.89 0.60 0.67 0.53
xargs.1 4,227 0.13 0.17 0.22 0.12
---------------------------------------------------
合計 2,810,784 91.94 46.84 40.81 43.15
ファイル名 サイズ LS ABC1R ABC2R ITSR
------------------------------------------------
repeat1.txt 2,600 0.14 0.17 0.22 0.12
repeat2.txt 5,200 0.25 0.25 0.30 0.20
repeat4.txt 10,400 0.47 0.44 0.52 0.41
repeat8.txt 20,800 0.95 0.99 1.07 0.96
# 時間の単位 : 秒
実行環境 : Windows XP, celeron 1.40 GHz, Python 2.4.2
ITSR が ITS にランクソートを適用した結果です。テキストファイルの場合、ITS にランクソートを適用した効果はほとんどないようです。バイナリファイルの場合、データによってはランクソートが効果的なようで、sum の実行速度は速くなりました。特に repeat*.txt のように繰り返しの多いデータは、ランクソートの効果により実行速度はとても速くなります。ランクソート法は、ITS でも最悪時の性能を改善する効果があるようです。
なお、これらの結果は M.Hiroi のコーディング、実行したマシン、コンパイラなどの環境に大きく依存しています。また、これらの環境だけではなく、データの種類によっても実行時間はかなり左右されます。興味のある方は、いろいろなデータをご自分の環境で試してみてください。
今回のプログラムを作成するに当たり、Yuta Mori さんの Web サイト white page で公開されているプログラムとドキュメントを参考にさせていただきました。素晴らしいプログラムとドキュメントを公開されている Yuta Mori さんに感謝いたします。
# coding: utf-8
#
# mksa51.py : Suffix Array の構築 (三分割法の改良)
#
# Copyright (C) 2007 Makoto Hiroi
#
import time, sys, os.path
from array import *
# 定数
SIZE = 256 * 256
LIMIT = 10
# 記号の取得
def getcode(x):
if x < data_size: return buff[x]
return -1
# ソート用比較関数
def compare(x, y, n):
for i in xrange(n, data_size - max(x, y)):
r = buff[x + i] - buff[y + i]
if r != 0: return r
return y - x
# 挿入ソート
def insert_sort(low, high, n):
for i in xrange(low + 1, high + 1):
temp = idx[i]
j = i - 1
while j >= low and compare(temp, idx[j], n) < 0:
idx[j + 1] = idx[j]
j -= 1
idx[j + 1] = temp
# 枢軸の選択
def select_pivot(low, high, n):
m = (high - low) / 4
a = getcode(idx[low + m] + n)
b = getcode(idx[low + m * 2] + n)
c = getcode(idx[low + m * 3] + n)
if a > b:
tmp = a
a = b
b = tmp
if b > c:
b = c
if a > b: b = a
return b
# マルチキークイックソート
def mqsort(low, high, n = 0):
while True:
if high - low <= LIMIT:
insert_sort(low, high, n)
return
# 枢軸
p = select_pivot(low, high, n)
# 4 分割
i = m1 = low
j = m2 = high
while True:
while i <= j:
k = getcode(idx[i] + n) - p
if k > 0: break
if k == 0:
idx[i], idx[m1] = idx[m1], idx[i]
m1 += 1
i += 1
while i <= j:
k = getcode(idx[j] + n) - p
if k < 0: break
if k == 0:
idx[j], idx[m2] = idx[m2], idx[j]
m2 -= 1
j -= 1
if i > j: break
idx[i], idx[j] = idx[j], idx[i]
i += 1
j -= 1
# 等しいデータを中央に集める
for k in xrange(min(m1 - low, i - m1)):
idx[low + k], idx[j - k] = idx[j - k], idx[low + k]
m1 = low + (i - m1)
for k in xrange(min(high - m2, m2 - j)):
idx[i + k], idx[high - k] = idx[high - k], idx[i + k]
m2 = high - (m2 - j) + 1
mqsort(low, m1 - 1, n)
mqsort(m2, high, n)
if m1 >= m2: break
low = m1
high = m2 - 1
n += 1
# Type C をソート
def sort_type_c():
# Type C を Type D, E を分離して Type D をソート
for x in xrange(0, 256):
for y in xrange(x + 1, 256):
c = (x << 8) + y
low = count_sum[c]
high = count_sum[c + 1]
m = low
for n in xrange(low, high):
i = idx[n]
if buff[i + 1] >= buff[i + 2]:
idx[n] = idx[m]
idx[m] = i
m += 1
if m - low > 1:
mqsort(low, m - 1, 2)
#
# Type E のセット
#
# 累積度数表をコピー
end = count_sum[:]
# 記号 0, 255 の区間はスキャンする必要なし
for x in xrange(254, 0, -1):
for y in xrange(255, x, -1):
c = (x << 8) + y
for n in xrange(count_sum[c + 1] - 1, count_sum[c] - 1, -1):
i = idx[n]
if i > 0 and buff[i - 1] < buff[i] < buff[i + 1]:
# Type E をセット
c1 = (buff[i - 1] << 8) + buff[i] + 1
end[c1] -= 1
idx[end[c1]] = i - 1
# Type A, B をセット
def set_type_ab():
start = [0] * 256
end = [0] * 256
# 最後尾のデータをセット
z = -1
if buff[data_size - 1] >= buff[data_size]:
c = (buff[data_size - 1] << 8) + buff[data_size]
idx[count_sum[c]] = data_size - 1
z = c
#
for i in xrange(0, 256):
for j in xrange(i, 256):
start[j] = count_sum[(j << 8) + i]
if (j << 8) + i == z: start[j] += 1
end[j] = count_sum[(j << 8) + i + 1] - 1
# 前から後ろへ
j = count_sum[i << 8]
while j < start[i]:
x = idx[j]
if x > 0 and buff[x - 1] >= buff[x]:
c = buff[x - 1]
idx[start[c]] = x - 1
start[c] += 1
j += 1
# 後ろから前へ
j = count_sum[(i + 1) << 8] - 1
while j > end[i]:
x = idx[j]
if x > 0 and buff[x - 1] >= buff[x]:
c = buff[x - 1]
idx[end[c]] = x - 1
end[c] -= 1
j -= 1
# 三分割法の改良
def suffix_sort():
global count_sum
# 分布数えソート
count = [0] * SIZE
count_sum = [0] * (SIZE + 1)
for x in xrange(data_size):
count[(buff[x] << 8) + buff[x + 1]] += 1
for x in xrange(1, SIZE + 1):
count_sum[x] = count[x - 1] + count_sum[x - 1]
for x in xrange(1, SIZE):
count[x] += count[x - 1]
for x in xrange(data_size - 1, -1, -1):
c = (buff[x] << 8) + buff[x + 1]
count[c] -= 1
idx[count[c]] = x
# Type C をソート
sort_type_c()
# Type A, B をセット
set_type_ab()
# suffix array の構築
def make_suffix_array(name):
global buff, idx, data_size
# 入力
data_size = os.path.getsize(name)
fin = open(name, 'rb')
buff = array('B')
buff.fromfile(fin, data_size)
buff.append(0)
buff.append(0)
idx = array('L')
for _ in xrange(data_size): idx.append(0)
# ソート
s = time.clock()
suffix_sort()
e = time.clock()
print "%.3f" % (e - s)
# 出力
name1 = os.path.basename(name) + '.idx'
fout = open(name1, 'wb')
idx.tofile(fout)
fin.close()
fout.close()
#
make_suffix_array(sys.argv[1])
# coding: utf-8
#
# mksa52.py : Suffix Array の構築 (improved two-stage suffix sorting)
#
# Copyright (C) 2007 Makoto Hiroi
#
import time, sys, os.path
from array import *
# 定数
SIZE = 256 * 256
LIMIT = 10
# 記号の取得
def getcode(x):
if x < data_size: return buff[x]
return -1
# ソート用比較関数
def compare(x, y, n):
for i in xrange(n, data_size - max(x, y)):
r = buff[x + i] - buff[y + i]
if r != 0: return r
return y - x
# 挿入ソート
def insert_sort(low, high, n):
for i in xrange(low + 1, high + 1):
temp = idx[i]
j = i - 1
while j >= low and compare(temp, idx[j], n) < 0:
idx[j + 1] = idx[j]
j -= 1
idx[j + 1] = temp
# 枢軸の選択
def select_pivot(low, high, n):
m = (high - low) / 4
a = getcode(idx[low + m] + n)
b = getcode(idx[low + m * 2] + n)
c = getcode(idx[low + m * 3] + n)
if a > b:
tmp = a
a = b
b = tmp
if b > c:
b = c
if a > b: b = a
return b
# マルチキークイックソート
def mqsort(low, high, n = 0):
while True:
if high - low <= LIMIT:
insert_sort(low, high, n)
return
# 枢軸
p = select_pivot(low, high, n)
# 4 分割
i = m1 = low
j = m2 = high
while True:
while i <= j:
k = getcode(idx[i] + n) - p
if k > 0: break
if k == 0:
idx[i], idx[m1] = idx[m1], idx[i]
m1 += 1
i += 1
while i <= j:
k = getcode(idx[j] + n) - p
if k < 0: break
if k == 0:
idx[j], idx[m2] = idx[m2], idx[j]
m2 -= 1
j -= 1
if i > j: break
idx[i], idx[j] = idx[j], idx[i]
i += 1
j -= 1
# 等しいデータを中央に集める
for k in xrange(min(m1 - low, i - m1)):
idx[low + k], idx[j - k] = idx[j - k], idx[low + k]
m1 = low + (i - m1)
for k in xrange(min(high - m2, m2 - j)):
idx[i + k], idx[high - k] = idx[high - k], idx[i + k]
m2 = high - (m2 - j) + 1
mqsort(low, m1 - 1, n)
mqsort(m2, high, n)
if m1 >= m2: break
low = m1
high = m2 - 1
n += 1
# Type B* のソート
def sort_type_b_star():
b_star_end = count_sum[:]
# 最後尾のデータは Type A
prev = 1
# Type B* の分離
for x in xrange(data_size - 2, -1, -1):
r = buff[x] - buff[x + 1]
if r == 0: continue
if r < 0 and prev > 0:
# Type B* をセット
c = (buff[x] << 8) + buff[x + 1]
idx[b_star_end[c]] = x
b_star_end[c] += 1
prev = r
# Type B* をソート
for x in xrange(0, 256):
for y in xrange(x + 1, 256):
c = (x << 8) + y
low = count_sum[c]
high = b_star_end[c]
if high - low > 1:
mqsort(low, high - 1, 2)
# Type A, B をセット
def set_type_ab():
end = [0] * 256
type_a_end = [0] * 256
# Type B をセット
for i in xrange(254, -1, -1):
for j in xrange(0, i + 1):
end[j] = count_sum[(j << 8) + i + 1] - 1
# 後ろから前へ
j = count_sum[(i + 1) << 8] - 1
while j > end[i]:
x = idx[j]
if x > 0 and buff[x - 1] <= buff[x]:
c = buff[x - 1]
idx[end[c]] = x - 1
end[c] -= 1
j -= 1
type_a_end[i] = end[i]
type_a_end[255] = count_sum[256 << 8] - 1
# Type A をセット
# 最後尾のデータをセット
if buff[data_size - 1] >= buff[data_size]:
c = (buff[data_size - 1] << 8) + buff[data_size]
idx[count_sum[c]] = data_size - 1
count_sum[c] += 1
# 前から Type A をセットする
for x in xrange(data_size):
i = idx[x]
if i > 0 and buff[i - 1] >= buff[i]:
c = (buff[i - 1] << 8) + buff[i]
if count_sum[c] <= type_a_end[buff[i - 1]]:
idx[count_sum[c]] = i - 1
count_sum[c] += 1
# improved two-stage suffix sorting
def suffix_sort():
global count_sum
# 累積度数表を求める
count = [0] * SIZE
count_sum = [0] * (SIZE + 1)
for x in xrange(data_size):
count[(buff[x] << 8) + buff[x + 1]] += 1
for x in xrange(1, SIZE + 1):
count_sum[x] = count[x - 1] + count_sum[x - 1]
# improved two-stage sort
# Type B* のソート
sort_type_b_star()
# Type A, B のセット
set_type_ab()
# suffix array の構築
def make_suffix_array(name):
global buff, idx, data_size
# 入力
data_size = os.path.getsize(name)
fin = open(name, 'rb')
buff = array('B')
buff.fromfile(fin, data_size)
buff.append(0)
idx = array('L')
for _ in xrange(data_size): idx.append(0)
# ソート
s = time.clock()
suffix_sort()
e = time.clock()
print "%.3f" % (e - s)
# 出力
name1 = os.path.basename(name) + '.idx'
fout = open(name1, 'wb')
idx.tofile(fout)
fin.close()
fout.close()
#
make_suffix_array(sys.argv[1])
# coding: utf-8
#
# mksa53.py : Suffix Array の構築
# (Improved two-stage suffix sorting + ランクソート)
#
# Copyright (C) 2007 Makoto Hiroi
#
import time, sys, os.path
from array import *
# 定数
SIZE = 256 * 256
LIMIT = 10
# ランクの取得
def getrank(x):
if x < data_size: return rank[x]
return -1
# ソート用比較関数
def compare(x, y, n):
for i in xrange(n, data_size - max(x, y), 2):
r = rank[x + i] - rank[y + i]
if r != 0: return r
return y - x
# 挿入ソート
def insert_sort(low, high, n):
for i in xrange(low + 1, high + 1):
temp = idx[i]
j = i - 1
while j >= low and compare(temp, idx[j], n) < 0:
idx[j + 1] = idx[j]
j -= 1
idx[j + 1] = temp
# ランクの設定
for i in xrange(low, high + 1):
rank[idx[i]] = i
# 枢軸の選択
def select_pivot(low, high, n):
m = (high - low) / 4
a = getrank(idx[low + m] + n)
b = getrank(idx[low + m * 2] + n)
c = getrank(idx[low + m * 3] + n)
if a > b:
tmp = a
a = b
b = tmp
if b > c:
b = c
if a > b: b = a
return b
# マルチキークイックソート
def mqsort(low, high, n = 0):
stack = []
while True:
if high - low <= LIMIT:
insert_sort(low, high, n)
if len(stack) == 0: break
low, high, n = stack.pop()
continue
#
p = select_pivot(low, high, n)
# 4 分割
i = m1 = low
j = m2 = high
while True:
while i <= j:
k = getrank(idx[i] + n) - p
if k > 0: break
if k == 0:
idx[i], idx[m1] = idx[m1], idx[i]
m1 += 1
i += 1
while i <= j:
k = getrank(idx[j] + n) - p
if k < 0: break
if k == 0:
idx[j], idx[m2] = idx[m2], idx[j]
m2 -= 1
j -= 1
if i > j: break
idx[i], idx[j] = idx[j], idx[i]
i += 1
j -= 1
# 等しいデータを中央に集める
for k in xrange(min(m1 - low, i - m1)):
idx[low + k], idx[j - k] = idx[j - k], idx[low + k]
m1 = low + (i - m1)
for k in xrange(min(high - m2, m2 - j)):
idx[i + k], idx[high - k] = idx[high - k], idx[i + k]
m2 = high - (m2 - j) + 1
if m2 <= high: stack.append((m2, high, n))
if m1 <= m2 - 1: stack.append((m1, m2 - 1, n + 2))
# < の部分から処理する
high = m1 - 1
# Type B* のソート
def sort_type_b_star():
b_star_end = count_sum[:]
# 最後尾のデータは Type A
prev = 1
# Type B* の分離
for x in xrange(data_size - 2, -1, -1):
r = buff[x] - buff[x + 1]
if r == 0: continue
if r < 0 and prev > 0:
# Type B* をセット
c = (buff[x] << 8) + buff[x + 1]
idx[b_star_end[c]] = x
b_star_end[c] += 1
prev = r
# Type B* をソート
for x in xrange(0, 256):
for y in xrange(x + 1, 256):
c = (x << 8) + y
low = count_sum[c]
high = b_star_end[c]
if high - low > 1:
mqsort(low, high - 1, 2)
# Type A, B をセット
def set_type_ab():
end = [0] * 256
type_a_end = [0] * 256
# Type B をセット
for i in xrange(254, -1, -1):
for j in xrange(0, i + 1):
end[j] = count_sum[(j << 8) + i + 1] - 1
# 後ろから前へ
j = count_sum[(i + 1) << 8] - 1
while j > end[i]:
x = idx[j]
if x > 0 and buff[x - 1] <= buff[x]:
c = buff[x - 1]
idx[end[c]] = x - 1
end[c] -= 1
j -= 1
type_a_end[i] = end[i]
type_a_end[255] = count_sum[256 << 8] - 1
# Type A をセット
# 最後尾のデータをセット
if buff[data_size - 1] >= buff[data_size]:
c = (buff[data_size - 1] << 8) + buff[data_size]
idx[count_sum[c]] = data_size - 1
count_sum[c] += 1
# 前から Type A をセットする
for x in xrange(data_size):
i = idx[x]
if i > 0 and buff[i - 1] >= buff[i]:
c = (buff[i - 1] << 8) + buff[i]
if count_sum[c] <= type_a_end[buff[i - 1]]:
idx[count_sum[c]] = i - 1
count_sum[c] += 1
# improved two-stage suffix sorting
def suffix_sort():
global count_sum
# 累積度数表とランクの生成
count = [0] * SIZE
count_sum = [0] * (SIZE + 1)
for x in xrange(data_size):
count[(buff[x] << 8) + buff[x + 1]] += 1
for x in xrange(1, SIZE + 1):
count_sum[x] = count[x - 1] + count_sum[x - 1]
for x in xrange(data_size):
c = (buff[x] << 8) + buff[x + 1]
rank[x] = count_sum[c + 1] - 1
# improved two-stage sort
# Type B* のソート
sort_type_b_star()
# Type A, B のセット
set_type_ab()
#
def make_suffix_array(name):
global buff, idx, data_size, rank
# 入力
data_size = os.path.getsize(name)
fin = open(name, 'rb')
buff = array('B')
buff.fromfile(fin, data_size)
buff.append(0)
idx = array('L')
rank = array('L')
for _ in xrange(data_size):
idx.append(0)
rank.append(0)
# ソート
s = time.clock()
suffix_sort()
e = time.clock()
print "%.3f" % (e - s)
# 出力
name1 = os.path.basename(name) + '.idx'
fout = open(name1, 'wb')
idx.tofile(fout)
fin.close()
fout.close()
#
make_suffix_array(sys.argv[1])