M.Hiroi's Home Page

Algorithms with Python

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

[ PrevPage | Python | NextPage ]

はじめに

今回は ratio-2 による三分割法の改良方法について説明します。

●ratio-2 による三分割法の改良

三分割法に ratio-2 を適用する場合、2 つの方法が考えられます。ひとつは、最初の三分割を ratio-2 で行う方法です。もうひとつは最初の三分割を ratio-1 で行い、そのあとで ratio-2 を適用する方法です。ここでは三分割法のプログラムを改造することで簡単に対応できる、後者の方法を採用することにします。

また、この方法でも ratio-2 に三分割法 (Si > Si+2, Si = Si+2, Si < Si+2) を適用できると考えられますが、Si = Si+2 の部分に他の Type と依存関係があるため、サフィックスの順序を決定できない場合があります。この部分を場合分けしてソートすることもできますが、プログラムがかなり複雑になると思われます。そこで、今回の ratio-2 は二分割にとどめます。ratio-2 の三分割は今後の課題にしたいと思います。

ratio-1 で三分割したあと ratio-2 を適用する場合、ソートするタイプで場合分けが異なります。次の図を見てください。

Type C をソートする場合、Type C を ratio-2 で二分割します。この場合、文字 a から順番にソートしていくので、Si > Si+2 (Type D) と Si <= Si+2 (Type E) に分けて、Type E をソートします。Type D は今までと同様にソート結果を使ってサフィックスの順序を決定できます。

Type A をソートする場合、Type A を ratio-2 で二分割します。この場合、文字 e からソートしていくので、Si >= Si+2 (Type D) と Si < Si+2 (Type E) に分けて、Type D をソートします。Type E はコピーによりサフィックスの順序を決定できます。

この処理は下図のように表すことができます。文字 b, c, d の部分だけを示します。

Type C に ratio-2 を適用すると、区間の前半が Type D で後半が Type E に分けらます。ソートするのは Type E だけで、Type D のサフィックスの順序はソート結果から決定できます。

たとえば、区間 bc, bd の Type D は bca... や bda... というサフィックスになるので、文字 a がソートされていればサフィックスの位置を決定できます。文字 a の Type C は Type E しかないので、この区間はすべてソートされます。同様に、区間 cd の Type D は cda... や cdb... というサフィックスになるので、文字 a と b がソートされていればサフィックスの順序を決定できます。次の図を見てください。

区間 cd の Type D で、区間 cda のインデックスは区間 a からセットされ、区間 cdb のインデックスは区間 b からセットされます。前から順番にセットしていくので、Type D の区間も前から順番にセットされていくことに注意してください。

Type A に ratio-2 を適用する場合も同じです。今度は区間の前半である Type D をソートして、区間の後半である Type E はソート結果によりサフィックスの順序が決定されます。

たとえば、区間 db, dc の Type E は dbe... や dce... というサフィックスになるので、記号 e がソートされていればサフィックスの順序を決定できます。記号 e の Type A は Type D しかないので、この区間はすべてソートされます。同様に、区間 cb の Type E は cbd... や cbe... というサフィックスになるので、記号 e と d がソートされていればサフィックスの順序を決定できます。次の図を見てください。

区間 cb の Type E で、区間 cbe のインデックスは区間 e からセットされ、区間 cbd の部分は区間 d からセットされます。後ろの文字からソートしていくので、Type E の区間も後ろからセットされていくことに注意してください。

●プログラムの作成

それでは、ratio-2 を適用した三分割法のプログラムを作りましょう。関数 suffix_sort() は三分割法 ratio-1 と同じです。Type C をソートする関数 sort_type_c() は次のようになります。

リスト : 三分割法(ratio-2) Type C をソート

def sort_type_c():
    global start2
    start2 = [0] * (SIZE + 1)
    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]
            start2[(x << 8) + y] = low
            m = low
            for n in range(low, high):
                i = idx[n]
                if buff[i] > buff[i + 2]:
                    idx[n] = idx[m]
                    idx[m] = i
                    m += 1
            if high - m > 1:
                mqsort(m, high - 1, 2)
    # Type A, B をセット
    set_type_ab()

Type C をソートする場合、Type C の区間を Type D と Type E に分割してから Type E の区間をソートします。Type D は前から後ろへインデックスをセットしていくので、Type D の先頭の位置を配列 start2 に格納しておきます。あとは文字 buff[i] と buff[i + 2] を比較して、データを Type D と Type E に分割します。そして、Type E のデータをマルチキークイックソート (mqsort) でソートします。

Type A をソートする関数 sort_type_a() では、記号の比較を buff[i] <= buff[i + 2] とし、Type E に Si = Si+2 のデータを含めます。そして、Type D をソートして Type E をコピーするので、start2 にはType E の最後尾の位置をセットして、インデックスを後ろから前へセットしていきます。

セット処理も難しくありません。次のリストを見てください。

リスト : 三分割法(ratio-2) Type A, B をセット

def set_type_ab():
    start = [0] * 256
    end = [0] * 256
    # 最後尾のデータをセット
    z1 = -1
    i = data_size
    if buff[i - 1] >= buff[i]:
        c = (buff[i - 1] << 8) + buff[i]
        idx[count_sum[c]] = data_size - 1
        z1 = c
    z2 = -1
    if buff[i - 2] > buff[i] and buff[i - 2] < buff[i - 1]:
        c = (buff[i - 2] << 8) + buff[i - 1]
        idx[count_sum[c]] = i - 2
        z2 = c
    #
    for i in range(0, 256):
        for j in range(i, 256):
            start[j] = count_sum[(j << 8) + i]
            if (j << 8) + i == z1: start[j] += 1
            if (j << 8) + i == z2: start[j] += 1
            end[j] = count_sum[(j << 8) + i + 1] - 1
        # 前から後ろへ
        j = count_sum[i << 8]
        while j < start[i]:
            # ratio-1
            x = idx[j]
            if x > 0 and buff[x - 1] >= buff[x]:
                c = buff[x - 1]
                idx[start[c]] = x - 1
                start[c] += 1
            # ratio-2
            if x > 1 and buff[x - 2] > buff[x] and buff[x - 2] < buff[x - 1]:
                c = (buff[x - 2] << 8) + buff[x - 1]
                idx[start2[c]] = x - 2
                start2[c] += 1
            j += 1
        # 後ろから前へ
        j = count_sum[(i + 1) << 8] - 1
        work = []
        while j > end[i]:
            # ratio-1
            x = idx[j]
            if x > 0 and buff[x - 1] >= buff[x]:
                c = buff[x - 1]
                idx[end[c]] = x - 1
                end[c] -= 1
            # ratio-2
            if x > 1 and buff[x - 2] > buff[x] and buff[x - 2] < buff[x - 1]:
                work.append(x - 2)
            j -= 1
        # work からコピーする
        while len(work) > 0:
            x = work.pop()
            c = (buff[x] << 8) + buff[x + 1]
            idx[start2[c]] = x
            start2[c] += 1

リストは少々長いですが、今までのセット処理に ratio-2 用の処理が追加されているだけです。最初に最後尾のデータをセットします。変数 z1 には ratio-1 でインデックスをセットした区間を、変数 z2 では ratio-2 でインデックスをセットした区間をセットします。そして、start にデータをセットするとき、z1, z2 と同じ区間であれば start の値を +1 するようにします。

Type E をソートして Type A, B, D のインデックスをセットする場合、前から後ろへセットするときは簡単です。buff[x - 2] > buff[x] and buff[x - 2] < buff[x - 1] を満たせば ratio-2 のデータです。start2 の位置へインデックスをセットします。後ろから前へセットする場合は、ratio-2 のデータを配列 work にセットします。work はスタックとして使用します。そのあとで、work からデータを取り出して、start2 の位置へインデックスをセットします。

Type D をソートして Type B, C, E へコピーする関数 set_type_bc() の場合は、インデックスを前から後ろへセットするときに work を使います。そして、ratio-2 のデータの判定が buff[x - 2] < buff[x] and buff[x - 2] > buff[x - 1] になることに注意してください。

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

●評価結果

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

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

  ファイル名      サイズ    LS        TS1       TS2      ABC'     ABC2
  ------------------------------------------------------------------------
  alice29.txt    152,089   0.965     0.669     0.605     0.608     0.556
                                   (74,618)  (53.507)  (67,459)  (48,384)
  asyoulik.txt   125,179   0.760     0.497     0.438     0.494     0.425
                                   (62,477)  (42,567)  (58,834)  (40,898)
  cp.html         24,603   0.141     0.149     0.123     0.146     0.152
                                   (12,901)   (8.981)  (11,181)   (7,814)
  fields.c        11,150   0.081     0.083     0.080     0.082     0.088
                                    (6,251)   (4,230)   (4,699)   (3,156)
  grammar.lsp      3,721   0.039     0.054     0.049     0.051     0.070
                                    (2,231)   (1,679)   (1,357)   (1,078)
  kennedy.xls  1,029,744   6.964     5.806     4.697     3.332     2.624
                                  (623,975) (467,929) (404,858) (201,689)
  lcet10.txt     426,754   3.142     2.218     1.977     1.679     1.567
                                  (219,172) (157,115) (194,732) (137,506)
  plrabn12.txt   481,861   3.430     2.158     1.835     1.972     1.798
                                  (234,061) (162,222) (224,466) (158,938)
  ptt5           513,216   6.015     -----     -----     1.119     1.067
                                                       (34,601)  (29,788)
  sum             38,240   0.277     0.538     0.398     0.337     0.278
                                   (23,496)  (18,559)  (13,175)  (10,187)
  xargs.1          4,227   0.039     0.042     0.046     0.053     0.070
                                    (2,053)   (1,351)   (1,981)   (1,333)
  ------------------------------------------------------------------------
  合計         2,810,784  21.853     -----     -----    11.054     8.695


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

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

三分割法 ratio-2 (ABC2) は、二段階ソート法 (TS1, TS2) や三分割法 ratio-1 (ABC') よりもソートした個数が少なくなり、その分だけ実行速度が速くなるファイルが多くなりました。とくに、kennedy.xls は ratio-2 の効果がとても大きく、ここまで速くなるとは M.Hiroi も予想していませんでした。三分割法に ratio-2 を適用した効果は十分に出ていると思います。ただし、ratio-2 のプログラムは少し複雑になるので、小さなファイルでは実行速度が遅くなる場合もあります。

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

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

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

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

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

三分割法 ratio-2 にランクソート (ABC2R) を適用すると、小さなテキストファイルは ABC2, ABC'R よりも遅くなりますが、それ以外のファイルでは速くなりました。repeat*.txt のように繰り返しの多いデータは、ランクソートの効果により実行速度は速くなりましたが、ratio-2 のプログラムが複雑になる分だけ、ratio-1 よりも少し遅くなるようです。

ランクソート法は、三分割法 ratio-2 でも最悪時の性能を改善する効果があるのは確かなようですが、小さなファイルではその効果は少ないようです。興味のある方は、もっと大きなファイルで試してみてください。

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

●参考文献, URL

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

●プログラムリスト1

#
# mksa43.py : Suffix Array の構築 (two-stage sort, 三分割法 ratio-2)
#
#             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
    # 最後尾のデータをセット
    z1 = -1
    i = data_size
    if buff[i - 1] >= buff[i]:
        c = (buff[i - 1] << 8) + buff[i]
        idx[count_sum[c]] = data_size - 1
        z1 = c
    z2 = -1
    if buff[i - 2] > buff[i] and buff[i - 2] < buff[i - 1]:
        c = (buff[i - 2] << 8) + buff[i - 1]
        idx[count_sum[c]] = i - 2
        z2 = c
    #
    for i in range(0, 256):
        for j in range(i, 256):
            start[j] = count_sum[(j << 8) + i]
            if (j << 8) + i == z1: start[j] += 1
            if (j << 8) + i == z2: start[j] += 1
            end[j] = count_sum[(j << 8) + i + 1] - 1
        # 前から後ろへ
        j = count_sum[i << 8]
        while j < start[i]:
            # ratio-1
            x = idx[j]
            if x > 0 and buff[x - 1] >= buff[x]:
                c = buff[x - 1]
                idx[start[c]] = x - 1
                start[c] += 1
            # ratio-2
            if x > 1 and buff[x - 2] > buff[x] and buff[x - 2] < buff[x - 1]:
                c = (buff[x - 2] << 8) + buff[x - 1]
                idx[start2[c]] = x - 2
                start2[c] += 1
            j += 1
        # 後ろから前へ
        j = count_sum[(i + 1) << 8] - 1
        work = []
        while j > end[i]:
            # ratio-1
            x = idx[j]
            if x > 0 and buff[x - 1] >= buff[x]:
                c = buff[x - 1]
                idx[end[c]] = x - 1
                end[c] -= 1
            # ratio-2
            if x > 1 and buff[x - 2] > buff[x] and buff[x - 2] < buff[x - 1]:
                work.append(x - 2)
            j -= 1
        # work からコピーする
        while len(work) > 0:
            x = work.pop()
            c = (buff[x] << 8) + buff[x + 1]
            idx[start2[c]] = x
            start2[c] += 1

# Type C をソート
def sort_type_c():
    global start2
    start2 = [0] * (SIZE + 1)
    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]
            start2[(x << 8) + y] = low
            m = low
            for n in range(low, high):
                i = idx[n]
                if buff[i] > buff[i + 2]:
                    idx[n] = idx[m]
                    idx[m] = i
                    m += 1
            if high - m > 1:
                mqsort(m, high - 1, 2)
    # Type A, B をセット
    set_type_ab()

# Type B, C をセット
def set_type_bc():
    start = [0] * 256
    end = [0] * 256
    # 最後尾のデータをセット
    i = data_size
    z1 = -1
    if buff[i - 1] <= buff[i]:
        c = (buff[i - 1] << 8) + buff[i]
        idx[count_sum[c]] = i - 1
        z1 = c
    z2 = -1
    if buff[i - 2] < buff[i] and buff[i - 2] > buff[i - 1]:
        c = (buff[i - 2] << 8) + buff[i - 1]
        idx[count_sum[c]] = i - 1
        z2 = 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 == z1: start[j] += 1
            if (j << 8) + i == z2: start[j] += 1
            end[j] = count_sum[(j << 8) + i + 1] - 1
        # 前から後ろへ
        j = count_sum[i << 8]
        work = []
        while j < start[i]:
            # ratio-1
            x = idx[j]
            if x > 0 and buff[x - 1] <= buff[x]:
                c = buff[x - 1]
                idx[start[c]] = x - 1
                start[c] += 1
            # ratio-2
            if x > 1 and buff[x - 2] < buff[x] and buff[x - 2] > buff[x - 1]:
                work.append(x - 2)
            j += 1
        # 後ろから前へ
        j = count_sum[(i + 1) << 8] - 1
        while j > end[i]:
            # ratio-1
            x = idx[j]
            if x > 0 and buff[x - 1] <= buff[x]:
                c = buff[x - 1]
                idx[end[c]] = x - 1
                end[c] -= 1
            # ratio-2
            if x > 1 and buff[x - 2] < buff[x] and buff[x - 2] > buff[x - 1]:
                c = (buff[x - 2] << 8) + buff[x - 1]
                idx[start2[c]] = x - 2
                start2[c] -= 1
            j -= 1
        # work からコピー
        while len(work):
            x = work.pop()
            c = (buff[x] << 8) + buff[x + 1]
            idx[start2[c]] = x
            start2[c] -= 1


# Type A をソート
def sort_type_a():
    global start2
    start2 = [0] * (SIZE + 1)
    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]
            start2[(x << 8) + y] = high - 1
            m = low
            for n in range(low, high):
                i = idx[n]
                if buff[i] >= buff[i + 2]:
                    idx[n] = idx[m]
                    idx[m] = i
                    m += 1
            if m - low > 1:
                mqsort(low, m - 1, 2)
    # Type B, C をセット
    set_type_bc()

# 三分割法 (ratio-2)
def suffix_sort():
    global count, 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)
        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

#
# mksa44.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
    # 最後尾のデータをセット
    z1 = -1
    i = data_size
    if buff[i - 1] >= buff[i]:
        c = (buff[i - 1] << 8) + buff[i]
        idx[count_sum[c]] = data_size - 1
        z1 = c
    z2 = -1
    if buff[i - 2] > buff[i] and buff[i - 2] < buff[i - 1]:
        c = (buff[i - 2] << 8) + buff[i - 1]
        idx[count_sum[c]] = i - 2
        z2 = c
    #
    for i in range(0, 256):
        for j in range(i, 256):
            start[j] = count_sum[(j << 8) + i]
            if (j << 8) + i == z1: start[j] += 1
            if (j << 8) + i == z2: start[j] += 1
            end[j] = count_sum[(j << 8) + i + 1] - 1
        # 前から後ろへ
        j = count_sum[i << 8]
        while j < start[i]:
            # ratio-1
            x = idx[j]
            if x > 0 and buff[x - 1] >= buff[x]:
                c = buff[x - 1]
                idx[start[c]] = x - 1
                start[c] += 1
            # ratio-2
            if x > 1 and buff[x - 2] > buff[x] and buff[x - 2] < buff[x - 1]:
                c = (buff[x - 2] << 8) + buff[x - 1]
                idx[start2[c]] = x - 2
                start2[c] += 1
            j += 1
        # 後ろから前へ
        j = count_sum[(i + 1) << 8] - 1
        work = []
        while j > end[i]:
            # ratio-1
            x = idx[j]
            if x > 0 and buff[x - 1] >= buff[x]:
                c = buff[x - 1]
                idx[end[c]] = x - 1
                end[c] -= 1
            # ratio-2
            if x > 1 and buff[x - 2] > buff[x] and buff[x - 2] < buff[x - 1]:
                work.append(x - 2)
            j -= 1
        # work からコピーする
        while len(work) > 0:
            x = work.pop()
            c = (buff[x] << 8) + buff[x + 1]
            idx[start2[c]] = x
            start2[c] += 1

# Type C をソート
def sort_type_c():
    global start2
    start2 = [0] * (SIZE + 1)
    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]
            start2[(x << 8) + y] = low
            m = low
            for n in range(low, high):
                i = idx[n]
                if buff[i] > buff[i + 2]:
                    idx[n] = idx[m]
                    idx[m] = i
                    m += 1
            if high - m > 1:
                for z in range(low, m):
                    rank[idx[z]] = m - 1
                mqsort(m, high - 1, 2)
    # Type A, B をセット
    set_type_ab()

# Type B, C をセット
def set_type_bc():
    start = [0] * 256
    end = [0] * 256
    # 最後尾のデータをセット
    i = data_size
    z1 = -1
    if buff[i - 1] <= buff[i]:
        c = (buff[i - 1] << 8) + buff[i]
        idx[count_sum[c]] = i - 1
        z1 = c
    z2 = -1
    if buff[i - 2] < buff[i] and buff[i - 2] > buff[i - 1]:
        c = (buff[i - 2] << 8) + buff[i - 1]
        idx[count_sum[c]] = i - 1
        z2 = 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 == z1: start[j] += 1
            if (j << 8) + i == z2: start[j] += 1
            end[j] = count_sum[(j << 8) + i + 1] - 1
        # 前から後ろへ
        j = count_sum[i << 8]
        work = []
        while j < start[i]:
            # ratio-1
            x = idx[j]
            if x > 0 and buff[x - 1] <= buff[x]:
                c = buff[x - 1]
                idx[start[c]] = x - 1
                start[c] += 1
            # ratio-2
            if x > 1 and buff[x - 2] < buff[x] and buff[x - 2] > buff[x - 1]:
                work.append(x - 2)
            j += 1
        # 後ろから前へ
        j = count_sum[(i + 1) << 8] - 1
        while j > end[i]:
            # ratio-1
            x = idx[j]
            if x > 0 and buff[x - 1] <= buff[x]:
                c = buff[x - 1]
                idx[end[c]] = x - 1
                end[c] -= 1
            # ratio-2
            if x > 1 and buff[x - 2] < buff[x] and buff[x - 2] > buff[x - 1]:
                c = (buff[x - 2] << 8) + buff[x - 1]
                idx[start2[c]] = x - 2
                start2[c] -= 1
            j -= 1
        # work からコピー
        while len(work):
            x = work.pop()
            c = (buff[x] << 8) + buff[x + 1]
            idx[start2[c]] = x
            start2[c] -= 1

# Type A をソート
def sort_type_a():
    global start2
    start2 = [0] * (SIZE + 1)
    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]
            start2[(x << 8) + y] = high - 1
            m = low
            for n in range(low, high):
                i = idx[n]
                if buff[i] >= buff[i + 2]:
                    idx[n] = idx[m]
                    idx[m] = i
                    m += 1
            if m - low > 1:
                mqsort(low, m - 1, 2)
    # Type B, C をセット
    set_type_bc()

# 三分割法 (ratio-2)
def suffix_sort():
    global count, 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)
        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 月 22 日
改訂 2022 年 9 月 25 日

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

[ PrevPage | Python | NextPage ]