M.Hiroi's Home Page

Algorithms with Python

接尾辞配列 (suffix array) [3]

[ PrevPage | Python | NextPage ]

はじめに

今回は「二段階ソート法 (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 つに分けることができます。このため、三分割法を「四分割法」と呼ぶこともあります。

それでは詳しく説明しましょう。次の図を見てください。

記号は {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 のサフィックスの順序を決める処理が少し複雑になります。次の図を見てください。

区間 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 range(256):
        for y in range(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 range(data_size):
        count[(buff[x] << 8) + buff[x + 1]] += 1
    for x in range(1, SIZE + 1):
        count_sum[x] = count[x - 1] + count_sum[x - 1]
    for x in range(1, SIZE):
        count[x] += count[x - 1]
    for x in range(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 range(0, 256):
        for j in range(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 の区間にインデックスを正しくセットすることができます。

あとはとくに難しいところはないと思います。説明は割愛いたしますので、詳細は プログラムリスト1 をお読みください。

●評価結果

それでは、実際に Canterbury Corpus で配布されているテストデータ The Canterbury Corpus で suffix array を作成してみましょう。結果は suffix_sort の時間で、ファイルの入出力処理は含んでいません。結果は次にようになりました。

            表 : suffix array の作成 (suffix_sort の時間)

  ファイル名      サイズ    LS      TS1             TS2              ABC
  -----------------------------------------------------------------------------------
  alice29.txt    152,089   0.965   0.669  (74,618)  0.605  (53.507)  0.608  (67,459)
  asyoulik.txt   125,179   0.760   0.497  (62,477)  0.438  (42,567)  0.494  (58,834)
  cp.html         24,603   0.141   0.149  (12,901)  0.123   (8.981)  0.169  (11,541)
  fields.c        11,150   0.081   0.083   (6,251)  0.080   (4,230)  0.084   (5,023)
  grammar.lsp      3,721   0.039   0.054   (2,231)  0.049   (1,679)  0.059   (1,824)
  kennedy.xls  1,029,744   6.964   5.806 (623,975)  4.697 (467,929)  4.391 (460,331)
  lcet10.txt     426,754   3.142   2.218 (219,172)  1.977 (157,115)  1.679 (194,732)
  plrabn12.txt   481,861   3.430   2.158 (234,061)  1.835 (162,222)  1.972 (224,466)
  ptt5           513,216   6.015   -----            -----            1.144  (40,306)
  sum             38,240   0.277   0.538  (23,496)  0.398  (18,559)  0.398  (14,560)
  xargs.1          4,227   0.039   0.042   (2,053)  0.046   (1,351)  0.056   (1,978)
  -----------------------------------------------------------------------------------
  合計         2,810,784  21.853   -----            -----           11.054

時間の単位 : 秒, () : ソートしたデータ数
実行環境 : Python 3.8.10, Ubuntu 20.04 LTS (WSL1), Intel Core i5-6200U 2.30GHz

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 のコーディング、実行したマシン、コンパイラなどの環境に大きく依存しています。また、これらの環境だけではなく、データの種類によっても実行時間はかなり左右されます。興味のある方は、いろいろなデータをご自分の環境で試してみてください。

●三分割法の改良

二段階ソート法は小さな記号から順番にソートしていきます。それでは、逆に大きな記号からソートしていくと、どうなるのでしょうか。次の図を見てください。

この場合、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 range(data_size):
        count[(buff[x] << 8) + buff[x + 1]] += 1
    for x in range(1, SIZE + 1):
        count_sum[x] = count[x - 1] + count_sum[x - 1]
    for x in range(1, SIZE):
        count[x] += count[x - 1]
    for x in range(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 range(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 について調べてみてください。

あとはとくに難しいところはないと思います。説明は割愛いたしますので、詳細は プログラムリスト2 をお読みください。

●評価結果

それでは、実際に Canterbury Corpus で配布されているテストデータ The Canterbury Corpus で suffix array を作成してみましょう。結果は suffix_sort の時間で、ファイルの入出力処理は含んでいません。結果は次にようになりました。

            表 : suffix array の作成 (suffix_sort の時間)

  ファイル名      サイズ    LS        TS1       TS2       ABC      ABC'
  ------------------------------------------------------------------------
  alice29.txt    152,089   0.965     0.669     0.605     0.608     0.608
                                   (74,618)  (53.507)  (67,459)  (67,459)
  asyoulik.txt   125,179   0.760     0.497     0.438     0.494     0.494
                                   (62,477)  (42,567)  (58,834)  (58,834)
  cp.html         24,603   0.141     0.149     0.123     0.169     0.146
                                   (12,901)   (8.981)  (11,541)  (11,181)
  fields.c        11,150   0.081     0.083     0.080     0.084     0.082
                                    (6,251)   (4,230)   (5,023)   (4,699)
  grammar.lsp      3,721   0.039     0.054     0.049     0.059     0.051
                                    (2,231)   (1,679)   (1,824)   (1,357)
  kennedy.xls  1,029,744   6.964     5.806     4.697     4.391     3.332
                                  (623,975) (467,929) (460,331) (404,858)
  lcet10.txt     426,754   3.142     2.218     1.977     1.679     1.679
                                  (219,172) (157,115) (194,732) (194,732)
  plrabn12.txt   481,861   3.430     2.158     1.835     1.972     1.972
                                  (234,061) (162,222) (224,466) (224,466)
  ptt5           513,216   6.015     -----     -----     1.144     1.119
                                                       (40,306)  (34,601)
  sum             38,240   0.277     0.538     0.398     0.398     0.337
                                   (23,496)  (18,559)  (14,560)  (13,175)
  xargs.1          4,227   0.039     0.042     0.046     0.056     0.053
                                 (2,053)    (1,351)    (1,978)    (1,981)
  ------------------------------------------------------------------------
  合計         2,810,784  21.853     -----     -----    11.054     9.873


  ファイル名   サイズ    LS    ABC'
  ----------------------------------------
  repeat1.txt   2,600  0.043  0.152  (99)
  repeat2.txt   5,200  0.075  0.274 (199)
  repeat4.txt  10,400  0.124  0.980 (399)
  repeat8.txt  20,800  0.241  3.849 (799)

時間の単位 : 秒, () : ソートしたデータ数
実行環境 : Python 3.8.10, Ubuntu 20.04 LTS (WSL1), Intel Core i5-6200U 2.30GHz

LS が Larsson Sadakane 法、TS1 が二段階ソート法 ratio-1、TS2 が二段階ソート法 ratio-2、ABC が三分割法、ABC' が三分割法改良版の結果です。個数が少ない Type をソートする三分割法の改良は kennedy.xls で大きな効果がありました。ソートするデータ数が少なくなった分だけ、実行速度は高速になりました。

個数が少ないタイプをソートする方法は、二段階ソート法を改善する簡単で効果的な方法だと思います。ただし、繰り返しの多いデータには大きな効果はないようです。

ところで、三分割法にランクソート法を適用することもできます。説明は割愛いたしますので、詳細は プログラムリスト3 をお読みください。結果は次のようになりました。

  ファイル名      サイズ    LS     TSR1    TSR2    ABC'    ABC'R
  ---------------------------------------------------------------
  alice29.txt    152,089   0.965   0.575   0.553   0.608   0.610
  asyoulik.txt   125,179   0.760   0.463   0.405   0.494   0.508
  cp.html         24,603   0.141   0.110   0.121   0.146   0.118
  fields.c        11,150   0.081   0.077   0.068   0.082   0.078
  grammar.lsp      3,721   0.039   0.042   0.049   0.051   0.057
  kennedy.xls  1,029,744   6.964   4.829   4.356   3.332   3.753
  lcet10.txt     426,754   3.142   1.875   1.637   1.697   1.730
  plrabn12.txt   481,861   3.430   2.039   1.748   1.972   1.895
  ptt5           513,216   6.015   -----   -----   1.119   0.961
  sum             38,240   0.277   0.202   0.184   0.337   0.148
  xargs.1          4,227   0.039   0.047   0.051   0.053   0.055
  ---------------------------------------------------------------
  合計         2,810,784  21.853   -----   -----   9.873   9.913

  ファイル名   サイズ    LS     TSR1    TSR2    ABC'    ABC'R
  ------------------------------------------------------------
  repeat1.txt   2,600   0.043   0.060   0.063   0.152   0.052
  repeat2.txt   5,200   0.075   0.095   0.107   0.274   0.074
  repeat4.txt  10,400   0.124   0.198   0.196   0.980   0.119
  repeat8.txt  20,800   0.241   0.434   0.433   3.849   0.198

時間の単位 : 秒
実行環境 : Python 3.8.10, Ubuntu 20.04 LTS (WSL1), Intel Core i5-6200U 2.30GHz

ABC'R が三分割法改良版にランクソートを適用した結果です。テキストファイルの場合、三分割法にランクソートを適用した効果はほとんどありません。バイナリファイルの場合、データによってはランクソートが効果的なようで、sum の実行速度はとても速くなりました。

特に repeat*.txt のように繰り返しの多いデータは、ランクソートの効果により実行速度はとても速くなりました。ランクソート法は、三分割法でも最悪時の性能を改善する効果があるようです。興味のある方は、もっと大きなファイルで試してみてください。

なお、これらの結果は M.Hiroi のコーディング、実行したマシン、コンパイラなどの環境に大きく依存しています。また、これらの環境だけではなく、データの種類によっても実行時間はかなり左右されます。興味のある方は、いろいろなデータをご自分の環境で試してみてください。

●参考文献, URL

  1. 儘田真吾, 『suffix array の構築 --- 二段階ソート法とその改良』, (リンク切れ)
  2. 広井誠, 『高性能圧縮ツール bsrc の改良 bsrc2 (前編)』, Interface 2004 年 10 月号, CQ出版社

●プログラムリスト1

#
# mksa4.py : Suffix Array の構築 (two-stage sort, 三分割法)
#
#            Copyright (C) 2007-2022 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 range(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 range(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 range(min(m1 - low, i - m1)):
            idx[low + k], idx[j - k] = idx[j - k], idx[low + k]
        m1 = low + (i - m1)
        for k in range(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 range(0, 256):
        for j in range(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 range(256):
        for y in range(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 range(data_size):
        count[(buff[x] << 8) + buff[x + 1]] += 1
    for x in range(1, SIZE + 1):
        count_sum[x] = count[x - 1] + count_sum[x - 1]
    for x in range(1, SIZE):
        count[x] += count[x - 1]
    for x in range(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)
    buff = array('B')
    with open(name, 'rb') as fin:
        buff.fromfile(fin, data_size)
        buff.append(0)
    idx = array('I')
    for _ in range(data_size): idx.append(0)
    # ソート
    s = time.time()
    suffix_sort()
    e = time.time()
    print('{:.3f}'.format(e - s))
    # 出力
    name1 = os.path.basename(name) + '.idx'
    with open(name1, 'wb') as fout:
        idx.tofile(fout)

# 実行
make_suffix_array(sys.argv[1])

●プログラムリスト2

#
# mksa41.py : Suffix Array の構築 (two-stage sort, 三分割法)
#
#            Copyright (C) 2007-2022 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 range(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 range(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 range(min(m1 - low, i - m1)):
            idx[low + k], idx[j - k] = idx[j - k], idx[low + k]
        m1 = low + (i - m1)
        for k in range(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 range(0, 256):
        for j in range(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 range(256):
        for y in range(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 range(255, -1, -1):
        for j in range(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 range(255, -1, -1):
        for y in range(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 range(data_size):
        count[(buff[x] << 8) + buff[x + 1]] += 1
    for x in range(1, SIZE + 1):
        count_sum[x] = count[x - 1] + count_sum[x - 1]
    for x in range(1, SIZE):
        count[x] += count[x - 1]
    for x in range(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 range(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)
    buff = array('B')
    with open(name, 'rb') as fin:
        buff.fromfile(fin, data_size)
        buff.append(0)
    idx = array('I')
    for _ in range(data_size): idx.append(0)
    # ソート
    s = time.time()
    suffix_sort()
    e = time.time()
    print('{:.3f}'.format(e - s))
    # 出力
    name1 = os.path.basename(name) + '.idx'
    with open(name1, 'wb') as fout:
        idx.tofile(fout)

# 実行
make_suffix_array(sys.argv[1])

●プログラムリスト3

#
# mksa42.py : Suffix Array の構築 (two-stage sort, 三分割法, ランクソート)
#
#            Copyright (C) 2007-2022 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 range(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 range(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 range(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 range(min(m1 - low, i - m1)):
            idx[low + k], idx[j - k] = idx[j - k], idx[low + k]
        m1 = low + (i - m1)
        for k in range(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 range(0, 256):
        for j in range(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 range(256):
        for y in range(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 range(255, -1, -1):
        for j in range(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 range(255, -1, -1):
        for y in range(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 range(data_size):
        count[(buff[x] << 8) + buff[x + 1]] += 1
    for x in range(1, SIZE + 1):
        count_sum[x] = count[x - 1] + count_sum[x - 1]
    for x in range(1, SIZE):
        count[x] += count[x - 1]
    for x in range(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 range(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)
    buff = array('B')
    with open(name, 'rb') as fin:
        buff.fromfile(fin, data_size)
        buff.append(0)
    idx = array('I')
    rank = array('I')
    for _ in range(data_size):
        idx.append(0)
        rank.append(0)
    # ソート
    s = time.time()
    suffix_sort()
    e = time.time()
    print('{:.3f}'.format(e - s))
    # 出力
    name1 = os.path.basename(name) + '.idx'
    with open(name1, 'wb') as fout:
        idx.tofile(fout)

# 実行
make_suffix_array(sys.argv[1])

初版 2007 年 7 月 21 日
改訂 2022 年 9 月 25 日

Copyright (C) 2007-2022 Makoto Hiroi
All rights reserved.

[ PrevPage | Python | NextPage ]