こんにちは、junkawaです。

本記事では、ブロックソートアルゴリズムの「巡回行列のソート」についてソースコードを交えて紹介します。

前の記事

bzip2を読む はじめに
bzip2を読む ブロックソート1

前回のおさらい

ブロックソートアルゴリズムは
　1. 入力データを1字ずつずらして巡回行列を作る
　2. 巡回行列の行をソートする
　3. ソート後の行列の最後の列を出力とする

2. 巡回行列のソートでは、下記の1. 2. を行います。

　1. 先頭の2バイトでソート
　2. 先頭バイトのシンボルの出現頻度が少ない順に下記を行う
　　(3バイト目以降のソート)
　　a. 先頭2バイトが異なる値、の行をソート
　　b. 先頭2バイトが同じ値、の行をソート

本記事では、「1. 先頭の2バイトでソート」について紹介します。

先頭の2バイトでソート

計数ソート (Counting sort) を使用して、先頭2バイトでソートします。

abracadabraの例だと、下記の右側がソート結果です。

ここで、注目してほしいのが、3バイト目以降はソートされていない点です。

右側の2行目 abracadabra、3行目 abraabracad は、3バイト目以降も含めてソートした場合には3行目のabraabracadの方が上にきます。

3バイト目以降のソートは、続く別のアルゴリズムで行います。これについては別の記事で紹介します。

データ構造

入力： char block[nlbock] (nblockは900,000)

各行は入力データを1バイトずつずらしたものなので、入力データの先頭からのインデックス(以後、行の先頭ブロック番号と呼びます)で各行を表現できます。

行の先頭ブロック番号を保持する int ptr[nblock]を使って、巡回行列を表現できます。

ptrを使うことで、block[nblock][nblock]の領域を確保する必要なく巡回行列を表現できます。

出力：int ptr[nblock]

ソート後の行列を、ptr で表現します。

ソートした順に、行の先頭ブロック番号をptrにセットします。

一時データ：頻度テーブル ftab[65536]

処理が進むにつれ、異なる用途で使用されます。

1) 各行の先頭2バイト(0〜65535)の頻度を保存する
　ftab[10] = 20
　は、巡回行列の中で先頭2バイトが 10 (1バイト目が0で2バイト目が10)となる行は20個ある、を表す。

2) 計数ソートで、ソート結果を保存する ptr のインデックス番号
　ftab[10] = 1020
　は、先頭2バイトが10となる行の先頭ブロック番号を保存する場所は、ptr[1019] 〜 ptr[1000] 、を表す。

処理の概要

1. 各行の先頭2バイト(0〜65535)の頻度を数え、ftabに保存
2. 頻度を積算して、ftabに上書き保存
3. 積算した値を ptr のインデックスに利用して、行の先頭ブロック番号を保存する

ソースコード紹介

https://github.com/junkawa/bzip2/tree/master/bzip2-1.0.6
mainSort()@blocksort.c L769 〜 L830

初期化

   769:    /*-- set up the 2-byte frequency table --*/
   770:    for (i = 65536; i >= 0; i--) ftab[i] = 0;
   771: 

先頭2バイト(0〜65535)の頻度テーブルを0で初期化する。

各行の先頭2バイト(0〜65535)の頻度を数え、ftabに保存

   772:    j = block[0] << 8;

block[0]<<8は、776行目で>>8されて、block[0]に戻して使用されます。

先頭2バイトの値を、(block[i]<<8) | block[i+1] で求めますが、blcok最後のblock[nblock-1] の次の値を巡回してblock[0]とするための処理です。776行目参照。

   773:    i = nblock-1;

先頭2バイトを、(block[nblock-1]<<8) | block[0] 、(block[nblock-2]<<8) | block[nblock-1]、とblock配列の後ろから数えていくので、for ループカウンタ i を配列最後の nblock-1 としています。

   774:    for (; i >= 3; i -= 4) {

4バイト単位で処理をしています。おそらく、キャッシュを利用した高速化のためだと思われます。

   775:       quadrant[i] = 0;

quadrantを0で初期化します。
quadrantは、 3バイト目以降のソートで、後のソートを高速化するためにソート結果を保存するために利用します。

   776:       j = (j >> 8) | ( ((UInt16)block[i]) << 8);

j は先頭2バイトの値となります。
block配列の先頭側(低位アドレス)の値が2バイトの上位バイトになります。

for ループ初回 (i=nblock-1)の時は、j = block[nblock-1]<<8 | block[0] となります。
for ループ2回目(i=nblock-5)の時は、j = block[nblock-5]<<8 | block[nblock-4] となります。

   777:       ftab[j]++;

先頭2バイト j をインデックスとして、頻度テーブル ftab をインクリメントして更新します。
こうすることで、ftab[j] は、先頭2バイトが j となる行の数(頻度)を表します。

   778:       quadrant[i-1] = 0;
   779:       j = (j >> 8) | ( ((UInt16)block[i-1]) << 8);
   780:       ftab[j]++;

上記と同様の内容です。
for ループ初回 (i=nblock-1)の時は、j = block[nblock-2]<<8 | block[nblock-1] となります。

   781:       quadrant[i-2] = 0;
   782:       j = (j >> 8) | ( ((UInt16)block[i-2]) << 8);
   783:       ftab[j]++;

上記と同様の内容です。
for ループ初回 (i=nblock-1)の時は、j = block[nblock-3]<<8 | block[nblock-2] となります。

  784:       quadrant[i-3] = 0;
  785:       j = (j >> 8) | ( ((UInt16)block[i-3]) << 8);
  786:       ftab[j]++;

上記と同様の内容です。
for ループ初回 (i=nblock-1)の時は、j = block[nblock-4]<<8 | block[nblock-3] となります。

   787:    }

   788:    for (; i >= 0; i--) {
   789:       quadrant[i] = 0;
   790:       j = (j >> 8) | ( ((UInt16)block[i]) << 8);
   791:       ftab[j]++;
   792:    }

nblockが4で割り切れない値の場合、このループが使用されます。

i = 2の時、j = block[2]<<8 | block[3] 。
i = 1の時、j = block[1]<<8 | block[2] 。
i = 0の時、j = block[0]<<8 | block[1] 。

今回の記事に関係ない処理

   794:    /*-- (emphasises close relationship of block & quadrant) --*/
   795:    for (i = 0; i < BZ_N_OVERSHOOT; i++) {
   796:       block   [nblock+i] = block[i];
   797:       quadrant[nblock+i] = 0;
   798:    }
   799: 

block[nblock+i] の部分は、後のソートを高速化するためにソート結果を保存するために利用します。

@bzlib_private.h

   186: #define BZ_N_RADIX 2
   187: #define BZ_N_QSORT 12
   188: #define BZ_N_SHELL 18
   189: #define BZ_N_OVERSHOOT (BZ_N_RADIX + BZ_N_QSORT + BZ_N_SHELL + 2)

BZ_N_OVERSHOOT は 34 となります。
BZ_N_RADIX、BZ_N_QSORT、BZ_N_SHELL、及びBZ_N_OVERSHOOTの詳細な意図については調査中です。

block、BZ_N_OVERSHOOT、quadrantの関係を図示します。

   802:    /*-- Complete the initial radix sort --*/
   803:    for (i = 1; i <= 65536; i++) ftab[i] += ftab[i-1];
   804: 

ftab[0]はそのまま、ftab[1]は、ftab[0]とftab[1]の合計。ftab[2]は、ftab[0]とftab[1]とftab[2]の合計。同様に、積算してきます。
ftab[65535]は、先頭2バイトが0〜65535となる全ての場合の合計数なので行数に等しくなり、nblockとなります。

積算した値を ptr のインデックスに利用して、行の先頭ブロック番号を保存する

積算後の ftab は、ソート結果を格納するptr内の位置を指します。
先頭2バイトが i となる行のブロック番号は、ptr[ ftab[i]-1 ] にセットします。

上記図では、ブロック番号15の行(入力データの先頭から15文字ずらした所から始まる行)は、先頭(block[15])が 0、2バイト目(block[16])が2なので、先頭2バイトの値は、2 (0<<8 | 2)となります。積算後 ftab [2]は 17となっており、ptr の ftab[2]-1 (=16)の位置へ、ブロック番号である15をセットします(ptr[16] = 15)。

では、この処理をコードで紹介します。

   805:    s = block[0] << 8;
   806:    i = nblock-1;

頻度を計測した時と同じく、もう一度、先頭2バイトを取得します。
前と同じく、i=nblock-1とし、後ろ(block[nblock-1])から順に先頭2バイトの値を取得します。

   807:    for (; i >= 3; i -= 4) {
   808:       s = (s >> 8) | (block[i] << 8);

先頭2バイトの値は s となります。

for ループ初回 (i=nblock-1)の時は、s = block[nblock-1]<<8 | block[0] となります。for ループ2回目(i=nblock-5)の時は、s = block[nblock-5]<<8 | block[nblock-4] となります。

   809:       j = ftab[s] -1;

先頭2バイトsに対応する、積算後のftab[s] から 1引いた値を j とします。

   810:       ftab[s] = j;

ftab[s] を j (ftab[s]-1)で更新します。これは、次に先頭2バイトが今回と同じ s となった時に、ptr[j] でブロック番号をセットする位置を前方にずらすためです。

   811:       ptr[j] = i;

ここで、ptrにブロック番号iをセットします。for ループ初回(i=nblock-1)の時は、ptr[j] = nblock-1 となります。

ftab[s]の値は(先頭2バイト) sが小さいほど小さくなります(803行目でftabの値が積算されているため)。ブロック番号iはptr[ ftab[s]-1 ]にセットされるので、先頭2バイトsが小さい値となる行のブロック番号iは、ptrの先頭の方(ftab[s]-1が小さい場所)にセットされます。結果として、ブロック番号は先頭2バイトで昇順にソートされ、ptrに格納されます。

備考
この計数ソートの実装は安定ソートだと思います。
ブロック番号は i = nblock-1 と後ろから探してきており、ptr にセットする位置も後ろから(ftab[s] はデクリメントされるため)セットしているため、先頭2バイトが同じ値となる行は、行番号(ブロック番号)の小さいものが prt の低い位置にセットします。
上記図の例では、先頭2バイトが0となるブロック番号100と211が、ptr[8]に100、ptr[9]に211がセットされています。
このように、ソート前の巡回行列のブロック番号の順序がソート後も保たれているので、安定ソートと言えます。

   812:       s = (s >> 8) | (block[i-1] << 8);
   813:       j = ftab[s] -1;
   814:       ftab[s] = j;
   815:       ptr[j] = i-1;

上記と同様の内容です。
for ループ初回 (i=nblock-1)の時は、s = block[nblock-2]<<8 | block[nblock-1] となります。
ここで、ptrにブロック番号i-1をセットします。for ループ初回(i=nblock-1)の時は、ptr[j] = nblock-2 となります。

   816:       s = (s >> 8) | (block[i-2] << 8);
   817:       j = ftab[s] -1;
   818:       ftab[s] = j;
   819:       ptr[j] = i-2;

上記と同様の内容です。
for ループ初回 (i=nblock-1)の時は、s = block[nblock-3]<<8 | block[nblock-2] となります。
ここで、ptrにブロック番号i-2をセットします。for ループ初回(i=nblock-1)の時は、ptr[j] = nblock-3 となります。

   820:       s = (s >> 8) | (block[i-3] << 8);
   821:       j = ftab[s] -1;
   822:       ftab[s] = j;
   823:       ptr[j] = i-3;
   824:    }

上記と同様の内容です。
for ループ初回 (i=nblock-1)の時は、s = block[nblock-4]<<8 | block[nblock-3] となります。
ここで、ptrにブロック番号i-3をセットします。for ループ初回(i=nblock-1)の時は、ptr[j] = nblock-4 となります。

   825:    for (; i >= 0; i--) {
   826:       s = (s >> 8) | (block[i] << 8);
   827:       j = ftab[s] -1;
   828:       ftab[s] = j;
   829:       ptr[j] = i;
   830:    }

nblockが4で割り切れない値の場合、このループが使用されます。

i = 2の時、s = block[2]<<8 | block[3] 。ptr[j] = 2。
i = 1の時、s = block[1]<<8 | block[2] 。ptr[j] = 1。
i = 0の時、s = block[0]<<8 | block[1] 。ptr[j] = 0。

以上で、ブロック番号nblock-1からブロック番号0までの全ての行について、計数ソートの結果がptrに保存されました。

まとめ

本記事では、ブロックソートアルゴリズムの「巡回行列のソート」について、先頭2バイトでソートする部分を紹介しました。