接尾辞配列 (suffix array) の続きです。今回は「二段階ソート法 (two-stage sort) 」の改良版である「三分割法」という方法を取り上げます。
前回説明したように、二段階ソート法はソートするデータ数を減らす画期的なアルゴリズムです。しかしながら、二段階ソート法にも弱点があります。一般に、バイナリデータはあまり得意ではありません。特に、同じ記号が多数続くデータ（連長の多いデータ）や繰り返しの多いデータでは極端に遅くなります。
連長の多いデータは、二段階ソート法を改良することで処理速度を高速にすることができます。儘田真吾氏 (参考文献 1) によると、この改良方法を「三分割法」と呼ぶそうです。三分割法は 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
三分割法は Si と Si+1 が等しい場合を Type B とし、Si < Si+1 を Type C とします。そして、Type C のサフィックスをソートするだけで、Type A と Type B のサフィックスの順序を決定することができます。Type A のサフィックスの順序は Type C により決定されますが、Type B のサフィックスの順序は Type A から決定される部分と Type C から決定される部分の 2 つに分けることができます。このため、三分割法を「四分割法」と呼ぶこともあります。
それでは詳しく説明しましょう。次の図を見てください。
（Ａ）Ｔｙｐｅ （Ｂ）ソートとセットの順番
┌─┬─┬─┬─┬─┬─┐ ┌─┬─┬─┬─┬─┬─┐
│ │ａ│ｂ│ｃ│ｄ│ｅ│Si+1 │ │ａ│ｂ│ｃ│ｄ│ｅ│Si+1
├─┼─┼─┼─┼─┼─┤ ├─┼─┼─┼─┼─┼─┤
│ａ│Ｂ│Ｃ│Ｃ│Ｃ│Ｃ│ │ａ│││─┼─┼─┼→│１
├─┼─┼─┼─┼─┼─┤ ├─┼┼┼─┼─┼─┼─┤
│ｂ│Ａ│Ｂ│Ｃ│Ｃ│Ｃ│ │ｂ│││││─┼─┼→│２
├─┼─┼─┼─┼─┼─┤ ├─┼┼┼┼┼─┼─┼─┤
│ｃ│Ａ│Ａ│Ｂ│Ｃ│Ｃ│ │ｃ│││││││─┼→│３
├─┼─┼─┼─┼─┼─┤ ├─┼┼┼┼┼┼┼─┼─┤
│ｄ│Ａ│Ａ│Ａ│Ｂ│Ｃ│ │ｄ│││││││││→│４
├─┼─┼─┼─┼─┼─┤ ├─┼┼┼┼┼┼┼┼┼─┤
│ｅ│Ａ│Ａ│Ａ│Ａ│Ｂ│ │ｅ│↓│↓│↓│↓│↓│
└─┴─┴─┴─┴─┴─┘ └─┴─┴─┴─┴─┴─┘
Si Si ５ ６ ７ ８ ９
→ : ソート区間
↓ : セット区間
図 : 三分割法による区間分け
記号は {a, b, c, d, e} の 5 種類とします。三分割法の場合、区間の Type は上図 (A) のように分けることができます。三分割法は Type C の区間をソートし、その結果を使って Type A と Type B の区間のサフィックスの順序を決定します。
図 (B) を見てください。まず最初に、二段階ソート法と同様に Type C の区間をソートします(図 5 (B) 1, 2, 3, 4)。次に Type A と Type B のインデックスをセットします。このとき、Type B のサフィックスの順序を決める処理が少し複雑になります。次の図を見てください。
│──── Set ───→│←─── Set ────│
┌───┬───┬─┬─┬─┬─┬───┬───┐
│ Ａ │ Ａ │B1│B2│B3│B4│ Ｃ │ Ｃ │
└───┴───┴─┴─┴─┴─┴───┴───┘
│- ca -│- cb -│----- cc -----│- cd -│- ce -│
B1 は ca, cb からセットされ、B2 は B1 からセットされる
B4 は ce, cd からセットされ、B3 は B4 からセットされる
図 : Type B の位置決定
区間 c を考えてみましょう。この場合、区間 ca, cb が Type A, 区間 cc が Type B, 区間 cd, ce が Type C になります。Type A の区間 ca, cb は、二段階ソート法と同様に区間 a と b のセット処理により、既にソートされています。あとは Type B の区間 cc のサフィックスの位置を決定します。
区間 cc の前半は区間 ca, cb のセット処理によりインデックスがセットされます。これは二段階ソート法と同じですが、これだけではなく、区間 cc の後半は区間 cd, ce からインデックスをセットする必要があります。このとき、上図に示すように区間 ce, cd は後方からセットしていきます。つまり、区間の先頭だけではなく、後方からも Type A, B のサフィックスを探して、インデックスをセットするのです。
ここが三分割法を理解するポイントです。区間 cd の先頭から Type B のサフィックスを探してインデックスをセットしようとすると、区間 cc でインデックスをセットする位置が決定できないのです。後方からセットしていくことで、区間 cc には後ろからインデックスをセットすることができます。つまり、区間 ce, cd から区間 cc の B4 にインデックスがセットされ、区間 cc の B4 から B3 にインデックスがセットされるわけです。これで、区間 cc のサフィックスの順序を決定することができます。
上図の場合、区間 a には Type A の区間はありません。この場合、Type C の区間をソートしたあと、その区間の後ろから Type A, B のサフィックスを探します。したがって、Type B の区間 aa は後ろからインデックスがセットされます。
このように、三分割法は記号 S の区間 SS をソートする必要がありません。したがって、連長の多いデータに三分割法を適用すれば、極めて高速に suffix array を構築することができます。
それでは、三分割法のプログラムを作りましょう。次のリストを見てください。
リスト : 三分割法
# Type C をソート
def sort_type_c():
for x in xrange(256):
for y in xrange(x + 1, 256):
low = count_sum[(x << 8) + y]
high = count_sum[(x << 8) + y + 1] - 1
if high - low >= 1:
mqsort(low, high, 2)
# Type A, B をセット
set_type_ab()
# 三分割法
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()
最初に、先頭 2 記号で分布数えソートを行います。この処理は今までと同じです。関数 sort_type_c は Type C の区間 (x > y) をソートします。そのあとで、関数 set_type_ab を呼び出して Type A, B のインデックスをセットします。
次は関数 set_type_ab を作ります。
リスト : 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):
# start, end をセット
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
関数 set_type_ab は、インデックスのセットを行う区間の先頭を配列 start に、最後尾を配列 end にセットします。今までのように、累積度数表 count_sum を直接書き換えるとプログラムは正常に動作しません。start と end の使い方がポイントになります。
まず最初に、最後尾のデータのインデックスをセットします。このとき、count_sum の値を +1 することはできないので、インデックスをセットした区間を変数 z に格納しておきます。そして、start にデータをセットするとき、z と同じ区間であれば start の値を +1 するようにします。
次に、Type A, B のインデックスをセットするため、記号 i の区間を前から後ろへ探索します。まず、start と end に Type A, B の区間の上限値と下限値をセットします。そして、前から後ろへインデックスをセットします。探索する範囲は記号 i の区間の先頭から Type B の区間の start[i] までになります。
ここで、Type B の区間にインデックスがセットされると start[i] の値も増えることに注意してください。これで、探索の範囲が広がるわけです。そして、Type B にセットされたインデックスを利用して、他のサフィックスのインデックスをセットすることができます。同様に、後ろから前へ探索するときは end[i] の値が減少するので、探索の範囲が広がります。これで、Type B の区間にインデックスを正しくセットすることができます。
あとはとくに難しいところはないと思います。説明は割愛いたしますので、詳細は プログラムリスト１ をお読みください。
それでは、実際に Canterbury Corpus で配布されているテストデータ The Canterbury Corpus で suffix array を作成してみましょう。結果は suffix_sort の時間で、ファイルの入出力処理は含んでいません。結果は次にようになりました。
表 : suffix array の作成 (suffix_sort の時間) ファイル名 サイズ LS TS1 TS2 ABC ------------------------------------------------------------------------------ alice29.txt 152,089 4.10 2.95 (74,618) 2.69 (53.507) 2.67 (67,459) asyoulik.txt 125,179 3.25 2.29 (62,477) 2.04 (42,567) 2.19 (58,834) cp.html 24,603 0.56 0.61 (12,901) 0.55 (8.981) 0.59 (11,541) fields.c 11,150 0.28 0.34 (6,251) 0.33 (4,230) 0.32 (5,023) grammar.lsp 3,721 0.12 0.16 (2,231) 0.18 (1,679) 0.18 (1,824) kennedy.xls 1,029,744 31.14 26.82 (623,975) 22.85 (467,929) 19.59 (460,331) lcet10.txt 426,754 12.62 10.59 (219,172) 9.17 (157,115) 8.41 (194,732) plrabn12.txt 481,861 14.52 9.28 (234,061) 8.05 (162,222) 8.86 (224,466) ptt5 513,216 24.33 ---- ---- 4.87 (40,306) sum 38,240 0.89 2.84 (23,496) 2.17 (18,559) 2.07 (14,560) xargs.1 4,227 0.13 0.14 (2,053) 0.17 (1,351) 0.17 (1,978) ------------------------------------------------------------------------------ 合計 2,810,784 91.94 ---- ---- 49.92 # 時間の単位 : 秒, () : ソートしたデータ数 実行環境 : Windows XP, celeron 1.40 GHz, Python 2.4.2
LS が Larsson Sadakane 法、TS1 が二段階ソート法 ratio-1、TS2 が二段階ソート法 ratio-2、ABC が三分割法の結果です。テキストファイルの場合、三分割法は ratio-1 よりソートしたデータ数が少なくなり、その分だけ高速になります。ratio-2 と比較すると、ソートしたデータ数は多くなりますが、ファイルによっては ratio-2 よりも速くなる場合があります。
バイナリファイルは ratio-2 よりも少し速くなるようです。特に、ptt5 のような同じ記号が連続しているデータの場合、ratio-1 と ratio-2 は計測不能でした。このようなデータでも、三分割法は高速に suffix array を作成することができます。Larsson Sadakane 法より数倍も速いのですから、三分割法の効果は十分に出ていると思います。ただし、繰り返しの多いデータでは三分割法でも極端に遅くなります。
なお、これらの結果は M.Hiroi のコーディング、実行したマシン、コンパイラなどの環境に大きく依存しています。また、これらの環境だけではなく、データの種類によっても実行時間はかなり左右されます。興味のある方は、いろいろなデータをご自分の環境で試してみてください。
二段階ソート法は小さな記号から順番にソートしていきます。それでは、逆に大きな記号からソートしていくと、どうなるのでしょうか。次の図を見てください。
Ｔｙｐｅ ソートとセットの順番
┌─┬─┬─┬─┬─┬─┐ ┌─┬─┬─┬─┬─┬─┐
│ │ａ│ｂ│ｃ│ｄ│ｅ│Si+1 │ │ａ│ｂ│ｃ│ｄ│ｅ│Si+1
├─┼─┼─┼─┼─┼─┤ ├─┼─┼─┼─┼─┼─┤
│ａ│Ｂ│Ｃ│Ｃ│Ｃ│Ｃ│ │ａ│↓│││││││││
├─┼─┼─┼─┼─┼─┤ ├─┼─┼┼┼┼┼┼┼┼┤
│ｂ│Ａ│Ｂ│Ｃ│Ｃ│Ｃ│ │ｂ│→│↓│││││││４
├─┼─┼─┼─┼─┼─┤ ├─┼─┼─┼┼┼┼┼┼┤
│ｃ│Ａ│Ａ│Ｂ│Ｃ│Ｃ│ │ｃ│─┼→│↓│││││３
├─┼─┼─┼─┼─┼─┤ ├─┼─┼─┼─┼┼┼┼┤
│ｄ│Ａ│Ａ│Ａ│Ｂ│Ｃ│ │ｄ│─┼─┼→│↓│││２
├─┼─┼─┼─┼─┼─┤ ├─┼─┼─┼─┼─┼┼┤
│ｅ│Ａ│Ａ│Ａ│Ａ│Ｂ│ │ｅ│─┼─┼─┼→│↓│１
└─┴─┴─┴─┴─┴─┘ └─┴─┴─┴─┴─┴─┘
Si Si ９ ８ ７ ６ ５
→ : ソート区間
↓ : セット区間
図 : 三分割法で Type A をソートする場合
この場合、Type A の区間をソートすることになります。そして、その結果を使って Type B, C のサフィックスの順序を決定します。インデックスのセットは大きな記号から順番に行います。まず、記号 e の区間を探索して区間 ee, de, ce, be, ae にインデックスをセットします(図 (B) 5)。次に、記号 d の区間を探索して、区間 dd, cd, bd, ad にインデックスをセットします(図 (B) 6)。あとはこれを繰り返すだけです(図 6 (B) 7, 8, 9)。
このように、二段階ソート法 (三分割法) は Type A の区間をソートしても実現することができます。そうであれば、Type A と Type C の個数を求めて、個数が少ないタイプをソートした方がよいでしょう。Type A の個数が少ない場合は、今までよりも処理時間を短縮できる可能性があります。
プログラムは次のようになります。
リスト : 三分割法の改良
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 A, C の個数をカウントする
type_a = 0
type_c = 0
for i in xrange(256):
type_a += count_sum[(i << 8) + i] - count_sum[i << 8]
type_c += count_sum[(i + 1) << 8] - count_sum[(i << 8) + i + 1]
# 少ない方をソート
if type_a > type_c:
sort_type_c()
else:
sort_type_a()
Type A の個数を変数 type_a に、Type C の個数を変数 type_c に求めます。記号を i とすると、Type A は区間 (i, 0) から (i, i - 1) までになります。同様に、Type C は区間 (i, i + 1) から (i, 0xff) までなので、累積度数表を使って簡単に計算することができます。
あとは、変数 type_a と type_c を比較して、type_c が少なければ関数 sort_type_c を呼び出して Type C の区間をソートします。type_a が少なければ関数 sort_type_a を呼び出して、Type A の区間をソートします。sort_type_a は記号 0xff から Type A の区間をソートし、関数 set_type_bc を呼び出して Type B, C のインデックスのセット処理を行います。
このように、三分割法は小さい記号からソートしても、大きな記号からソートしても成立します。それでは、個数が少ない記号からソートする方法も考えられるでしょう。実は、bzip2 の作者 Julian Seward が提案した Copy Mehtod がそうなのです。Copy Method は二段階ソート法と同様にサフィックス間の関係を利用したアルゴリズムで、個数の少ない記号から順番にソートし、その結果を使って他の区間のサフィックスの順序を決定します。興味のある方は Copy Method について調べてみてください。
あとはとくに難しいところはないと思います。説明は割愛いたしますので、詳細は プログラムリスト２ をお読みください。
それでは、実際に Canterbury Corpus で配布されているテストデータ The Canterbury Corpus で suffix array を作成してみましょう。結果は suffix_sort の時間で、ファイルの入出力処理は含んでいません。結果は次にようになりました。
表 : suffix array の作成 (suffix_sort の時間)
ファイル名 サイズ LS TS1 TS2 ABC ABC'
--------------------------------------------------------------------------
alice29.txt 152,089 4.10 2.95 2.69 2.67 2.67
(74,618) (53.507) (67,459) (67,459)
asyoulik.txt 125,179 3.25 2.29 2.04 2.19 2.19
(62,477) (42,567) (58,834) (58,834)
cp.html 24,603 0.56 0.61 0.55 0.59 0.59
(12,901) (8.981) (11,541) (11,181)
fields.c 11,150 0.28 0.34 0.33 0.32 0.30
(6,251) (4,230) (5,023) (4,699)
grammar.lsp 3,721 0.12 0.16 0.18 0.18 0.16
(2,231) (1,679) (1,824) (1,357)
kennedy.xls 1,029,744 31.14 26.82 22.85 19.59 16.54
(623,975) (467,929) (460,331) (404,858)
lcet10.txt 426,754 12.62 10.59 9.17 8.41 8.41
(219,172) (157,115) (194,732) (194,732)
plrabn12.txt 481,861 14.52 9.28 8.05 8.86 8.86
(234,061) (162,222) (224,466) (224,466)
ptt5 513,216 24.33 ---- ---- 4.87 4.66
(40,306) (34,601)
sum 38,240 0.89 2.84 2.17 2.07 1.83
(23,496) (18,559) (14,560) (13,175)
xargs.1 4,227 0.13 0.14 0.17 0.17 0.17
(2,053) (1,351) (1,978) (1,981)
--------------------------------------------------------------------------
合計 2,810,784 91.94 ---- ---- 49.92 46.38
ファイル名 サイズ LS ABC'
----------------------------------------
repeat1.txt 2,600 0.14 0.44 (99)
repeat2.txt 5,200 0.25 1.32 (199)
repeat4.txt 10,400 0.47 4.74 (399)
repeat8.txt 20,800 0.95 18.25 (799)
# 時間の単位 : 秒, () : ソートしたデータ数
実行環境 : Windows XP, celeron 1.40 GHz, Python 2.4.2
LS が Larsson Sadakane 法、TS1 が二段階ソート法 ratio-1、TS2 が二段階ソート法 ratio-2、ABC が三分割法、ABC' が三分割法改良版の結果です。個数が少ない Type をソートする三分割法の改良は kennedy.xls で大きな効果がありました。ソートするデータ数が少なくなった分だけ、実行速度は高速になりました。
個数が少ないタイプをソートする方法は、二段階ソート法を改善する簡単で効果的な方法だと思います。ただし、繰り返しの多いデータには大きな効果はないようです。
ところで、三分割法にランクソート法を適用することもできます。説明は割愛いたしますので、詳細は プログラムリスト３ をお読みください。結果は次のようになりました。
表 : suffix array の作成 (suffix_sort の時間)
ファイル名 サイズ LS TSR1 TSR2 ABC' ABC'R
----------------------------------------------------------
alice29.txt 152,089 4.10 2.96 2.48 2.67 2.73
asyoulik.txt 125,179 3.25 2.37 1.94 2.19 2.28
cp.html 24,603 0.56 0.49 0.45 0.59 0.48
fields.c 11,150 0.28 0.28 0.28 0.30 0.27
grammar.lsp 3,721 0.12 0.15 0.17 0.16 0.16
kennedy.xls 1,029,744 31.14 26.50 20.75 16.54 18.19
lcet10.txt 426,754 12.62 9.33 7.54 8.41 8.34
plrabn12.txt 481,861 14.52 9.69 7.75 8.86 9.43
ptt5 513,216 24.33 ---- ---- 4.66 4.19
sum 38,240 0.89 0.93 0.85 1.83 0.60
xargs.1 4,227 0.13 0.14 0.17 0.17 0.17
----------------------------------------------------------
合計 2,810,784 91.94 ---- ---- 46.38 46.84
ファイル名 サイズ LS ABC' ABC'R
------------------------------------------
repeat1.txt 2,600 0.14 0.44 0.17
repeat2.txt 5,200 0.25 1.32 0.25
repeat4.txt 10,400 0.47 4.74 0.44
repeat8.txt 20,800 0.95 18.25 0.99
# 時間の単位 : 秒
実行環境 : Windows XP, celeron 1.40 GHz, Python 2.4.2
ABC'R が三分割法改良版にランクソートを適用した結果です。テキストファイルの場合、三分割法にランクソートを適用した効果はほとんどありません。バイナリファイルの場合、データによってはランクソートが効果的なようで、sum の実行速度はとても速くなりました。特に repeat*.txt のように繰り返しの多いデータは、ランクソートの効果により実行速度はとても速くなりました。ランクソート法は、三分割法でも最悪時の性能を改善する効果があるようです。それでも、Larsson Sadakane 法の方が少し速いですね。興味のある方は、もっと大きなファイルで試してみてください。
なお、これらの結果は M.Hiroi のコーディング、実行したマシン、コンパイラなどの環境に大きく依存しています。また、これらの環境だけではなく、データの種類によっても実行時間はかなり左右されます。興味のある方は、いろいろなデータをご自分の環境で試してみてください。
# coding: utf-8
#
# mksa4.py : Suffix Array の構築 (two-stage sort, 三分割法)
#
# 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 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
# Type C をソート
def sort_type_c():
for x in xrange(256):
for y in xrange(x + 1, 256):
low = count_sum[(x << 8) + y]
high = count_sum[(x << 8) + y + 1] - 1
if high - low >= 1:
mqsort(low, high, 2)
#
set_type_ab()
# 三分割法
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()
# 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
#
# mksa41.py : Suffix Array の構築 (two-stage sort, 三分割法)
#
# 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 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
# Type C をソート
def sort_type_c():
for x in xrange(256):
for y in xrange(x + 1, 256):
low = count_sum[(x << 8) + y]
high = count_sum[(x << 8) + y + 1] - 1
if high - low >= 1:
mqsort(low, high, 2)
#
set_type_ab()
# Type B, C をセット
def set_type_bc():
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(255, -1, -1):
for j in xrange(0, i + 1):
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
# Type A をソート
def sort_type_a():
for x in xrange(255, -1, -1):
for y in xrange(0, x):
low = count_sum[(x << 8) + y]
high = count_sum[(x << 8) + y + 1] - 1
if high - low >= 1:
mqsort(low, high, 2)
#
set_type_bc()
# 三分割法
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 A, C の個数をカウントする
type_a = 0
type_c = 0
for i in xrange(256):
type_a += count_sum[(i << 8) + i] - count_sum[i << 8]
type_c += count_sum[(i + 1) << 8] - count_sum[(i << 8) + i + 1]
# 少ない方をソート
if type_a > type_c:
sort_type_c()
else:
sort_type_a()
# 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
#
# mksa42.py : Suffix Array の構築 (two-stage sort, 三分割法, ランクソート)
#
# 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 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
# Type C をソート
def sort_type_c():
for x in xrange(256):
for y in xrange(x + 1, 256):
low = count_sum[(x << 8) + y]
high = count_sum[(x << 8) + y + 1] - 1
if high - low >= 1:
mqsort(low, high, 2)
#
set_type_ab()
# Type B, C をセット
def set_type_bc():
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(255, -1, -1):
for j in xrange(0, i + 1):
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
# Type A をソート
def sort_type_a():
for x in xrange(255, -1, -1):
for y in xrange(0, x):
low = count_sum[(x << 8) + y]
high = count_sum[(x << 8) + y + 1] - 1
if high - low >= 1:
mqsort(low, high, 2)
#
set_type_bc()
# 三分割法
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
rank[x] = count_sum[c + 1] - 1
# Type A, C の個数をカウントする
type_a = 0
type_c = 0
for i in xrange(256):
type_a += count_sum[(i << 8) + i] - count_sum[i << 8]
type_c += count_sum[(i + 1) << 8] - count_sum[(i << 8) + i + 1]
# 少ない方をソート
if type_a > type_c:
sort_type_c()
else:
sort_type_a()
# suffix array の構築
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])