「整列 (sorting)」はある規則に従ってデータを順番に並べ換える操作です。たとえば、データが整数であれば大きい順に並べる、もしくは小さい順に並べます。一般に、このような操作を「ソート」と呼んでいます。ソートは昔から研究されている分野で、優秀なアルゴリズムが確立しています。Yet Another Clang Problems (2) では、バブルソート、選択ソート、単純挿入ソートを取り上げました。要素数を N とすると、これらは N2 に比例する遅いソートです。今回はもっと速いソートを紹介しましょう。
ソートは大きく分けると、内部ソート (internal sort) と外部ソート (external sort) にわかれます。内部ソートはすべてのデータをメモリに読み込んでソートします。外部ソートはメモリにすべて読み込むことができない巨大なデータをソートするときに使われ、外部記憶装置に途中経過を記憶させながらソートします。
今回作成するプログラムは内部ソートで、データは浮動小数点数 (double) とします。データは配列に格納します。データは乱数のほかにも、山型データ (中央のデータがいちばん大きく端にいくほど小さいデータになる)、ソート済みのデータ (昇順)、逆順にソートされたデータ (降順) の 4 種類を用意して、実行時間を比較してみましょう。
次に、ソートの安定性について説明します。安定 (stable) なソートとは、ソートキーが等しい場合、入力された順番が崩れないソートのことです。逆に、不安定なソートとは、入力された順番とは異なる結果になるソートのことです。次の図を見てください。
元のデータ 安定なソート 不安定なソート
------------------------------------------------
(123, 'abc') (123, 'abc') (789, 'abc')
(456, 'def') (789, 'abc') (123, 'abc')
(789, 'abc') (456, 'def') (456, 'def')
図 : ソートの安定性
2 番目の要素をキーにソートした場合、1 行目と 3 行目はソートキーが等しいですね。安定なソートであればソート結果が 1, 3, 2 の順番になり、入力時の位置関係が保たれています。不安定なソートはソート結果が、たとえば 3, 1, 2 の順番となり、入力時の位置関係が崩れます。
実用的なアプリケーションを作成する場合、ソートの安定性が重要になる場合があります。一般に、単純なソート（遅いソート）は安定で、複雑なソート（高速なソート）は安定ではない、という傾向があります。
次はテストデータを作成する関数を作りましょう。
リスト : テストデータの作成
// 昇順
void make_sequence(double buff[], int size)
{
for (int i = 0; i < size; i++) buff[i] = i;
}
// 降順
void make_reverse(double buff[], int size)
{
for (int i = 0, j = size - 1; i < size; i++, j--)
buff[i] = j;
}
// 山型
void make_yama(double buff[], int size)
{
int j = 0;
for (int i = 0; i < size / 2; i++, j += 2)
buff[i] = j;
j--;
for (int i = size / 2; i < size; i++, j -= 2)
buff[i] = j;
}
どの関数も簡単なので説明は不要でしょう。次は乱数データの作成ですが、その前に乱数について簡単に説明します。
私たちが適当な数字を決める場合、サイコロを使えば 1 から 6 までの数字を簡単に決めることができます。たとえば、サイコロを振って出た目を記録したら次のようになったとしましょう。
5, 2, 1, 2, 6, 3, 4, 3, 1, 5, .....
サイコロの目はイカサマをしないかぎり、出る確率が 1 / 6 で規則性はまったくありません。したがって、数字の出る順番には規則性はなく、まったくでたらめになります。いま 2 が出たから次は 1 が出るとか、3, 4 と続いたから次は 5 が出るなどのように、前に出た数字から次に出る数字を予測することはできないのです。このように、でたらめに並んだ数列を「乱数列 (random sequence)」といい、乱数列の中の一つ一つの数字を「乱数 (random numbers)」といいます。
コンピュータは決められた手順（プログラム）を高速に実行することは得意なのですが、まったくでたらめの数を作れといわれると、とたんに困ってしまいます。そこで、何かしらの数式をプログラムして、それを実行することで乱数を発生させます。厳密にいえば乱数ではありませんが、それを乱数としてみなして使うことにするのです。このような乱数を「疑似乱数 (pseudo-random numbers)」といいます。
Ｃ言語の標準ライブラリには乱数を生成する関数 rand があります。
<stdlib.h> int rand(void);
この関数は 0 から RAND_MAX (32767 以上の値) までの乱数を生成します。M.Hiroi が使用している Clang の場合、RAND_MAX の値は 2147483647 で、int の最大値 (INT_MAX) と同じです。
それでは実際に試してみましょう。
リスト : 乱数の生成 (rand.c)
#include <stdio.h>
#include <stdlib.h>
int main(void)
{
for (int i = 0; i < 8; i++)
printf("%d\n", rand());
return 0;
}
$ clang rand.c $ ./a.out 1804289383 846930886 1681692777 1714636915 1957747793 424238335 719885386 1649760492
ところで、この乱数を使うためにはちょっとしたテクニックが必要になります。実は、このプログラムを実行すると、必ずこの乱数列が生成されます。一般に、擬似乱数を生成する場合、種 (seed : シード) となる数値が必要になります。Ｃ言語の場合、シードは大域変数に定義されているのが普通なので、初期値は 0 になります。その直後に rand を実行すれば、必ず同じ乱数列が生成されるというわけです。
異なる乱数列を生成するには、シードに違う値を設定すればいいのです。シードに値を設定するには関数 srand を使います。
<stdlib.h> void srand(unsigned int seed);
この設定に rand を使うことはできませんね。いちばん簡単な方法は、シードの設定に現在時刻を使うことです。Ｃ言語の場合、関数 time で現在時刻を求めることができます。
<time.h> time_t time(time_t *timer);
引数 timer が NULL でなければ、現在時刻を *timer にセットします。time_t は整数値なので、返り値を srand に渡してシードを設定すると、プログラムを実行するたびに違う乱数列を発生させることができます。今回は srand を使わなくてもいいのですが、簡単なゲームを作る場合はシードを設定してください。
次は乱数データを作成するため、配列をシャッフルする関数 shuffle を作ります。
リスト : 配列のシャッフル
void shuffle(double *buff, int size)
{
for (int i = size - 1; i > 0; i--) {
int j = rand() % (i + 1);
double tmp = buff[i];
buff[i] = buff[j];
buff[j] = tmp;
}
}
プログラムは簡単です。0 から i 番目の要素を乱数で選んで、それと i 番目の要素を交換します。rand() % (i + 1) で 0 から i までの乱数を得ることができます。次に、0 から i - 1 番目の要素を乱数で選んで、それと i - 1 番目の要素を交換します。つまり、配列の前半部分 (0 番目から i 番目) が未選択の要素で、後半部分 (i + 1 番目から size - 1 番目) が乱数で選んだ要素になります。未選択の要素がなくなったら処理を終了します。この方法を「Fisher-Yates shuffle アルゴリズム」といいます。
make_sequence で昇順のデータを作成し、そのあとで shuffle を実行すれば、乱数データを生成することができます。
次はソートの実行時間を計測するプログラムを作りましょう。時間計測には、ライブラリ関数 clock を使うと簡単です。
<time.h> clock_t clock(void);
関数 clock は、プログラム実行開始からの経過時間（プロセッサ時間）を取得します。返り値の clock_t は整数型です。関数を呼び出す前と呼び出し後の差分をとれば、関数の実行時間を求めることができます。秒数に変換する場合は、マクロ CLOCKS_PER_SEC で割り算します。
プログラムは次のようになります。
リスト : ソートのテスト
#define N 80000
double buff[N];
bool check(int size)
{
for (int i = 1; i < size; i++) {
if (buff[i -1] > buff[i]) {
printf("sort error\n");
return false;
}
}
return true;
}
void test(void (*func)(double *, int))
{
printf(" 個数 : 乱数 : 昇順 : 逆順 : 山型 \n");
printf("-------+-------+-------+-------+-------\n");
for (int i = 10000; i <= N; i *= 2) {
clock_t s;
printf("%6d : ", i);
make_sequence(buff, i);
shuffle(buff, i);
s = clock();
(*func)(buff, i);
printf("%.3f : ", (double)(clock() - s) / CLOCKS_PER_SEC);
check(i);
make_sequence(buff, i);
s = clock();
(*func)(buff, i);
printf("%.3f : ", (double)(clock() - s) / CLOCKS_PER_SEC);
check(i);
make_reverse(buff, i);
s = clock();
(*func)(buff, i);
printf("%.3f : ", (double)(clock() - s) / CLOCKS_PER_SEC);
check(i);
make_yama(buff, i);
s = clock();
(*func)(buff, i);
printf("%.3f\n", (double)(clock() - s) / CLOCKS_PER_SEC);
check(i);
}
}
関数 test の引数 func にはソート関数を渡します。実行時間は整数のまま割り算すると 0 になる場合があるので、double に「型変換 (キャスト)」してから CLOCKS_PER_SEC で割り算します。キャストは値の前に変換したいデータ型をカッコで囲んで付けるだけです。あとは特に難しいところないでしょう。
まず最初に、バブルソート、選択ソート、単純挿入ソートの実行時間を計測しましょう。プログラムと実行結果を下記に示します。
リスト : バブルソート、選択ソート、単純挿入ソート
void buble_sort(double *buff, int size)
{
for (int i = 0; i < size; i++) {
for (int j = size - 1; j > i; j--) {
if (buff[j] < buff[j - 1]) {
double temp = buff[j];
buff[j] = buff[j - 1];
buff[j - 1] = temp;
}
}
}
}
void select_sort(double *buff, int size)
{
for (int i = 0; i < size - 1; i++) {
double min = buff[i];
int n = i;
for (int j = i + 1; j < size; j++) {
if (buff[j] < min) {
min = buff[j];
n = j;
}
}
buff[n] = buff[i];
buff[i] = min;
}
}
void insert_sort(double *buff, int size)
{
for (int i = 1; i < size; i++) {
int j;
double temp = buff[i];
for (j = i - 1; j >= 0 && temp < buff[j]; j--) {
buff[j + 1] = buff[j];
}
buff[j + 1] = temp;
}
}
表 : バブルソート (単位 : 秒) 個数 : 乱数 : 昇順 : 逆順 : 山型 -------+-------+-------+-------+------- 10000 : 0.148 : 0.036 : 0.091 : 0.064 20000 : 0.679 : 0.148 : 0.367 : 0.257 40000 : 2.985 : 0.589 : 1.467 : 1.026 80000 :12.626 : 2.376 : 5.869 : 4.120
表 : 選択ソート (単位 : 秒) 個数 : 乱数 : 昇順 : 逆順 : 山型 -------+-------+-------+-------+------- 10000 : 0.075 : 0.072 : 0.073 : 0.073 20000 : 0.295 : 0.291 : 0.292 : 0.293 40000 : 1.173 : 1.166 : 1.169 : 1.166 80000 : 4.674 : 4.667 : 4.666 : 4.667
表 : 単純挿入ソート (単位 : 秒) 個数 : 乱数 : 昇順 : 逆順 : 山型 -------+-------+-------+-------+------- 10000 : 0.021 : 0.000 : 0.043 : 0.023 20000 : 0.084 : 0.000 : 0.164 : 0.082 40000 : 0.342 : 0.000 : 0.663 : 0.333 80000 : 1.334 : 0.000 : 2.677 : 1.342
乱数データではバブルソートが一番遅い結果になりました。選択ソートはデータの種類によらず一定の性能を発揮しています。それから、挿入ソートは昇順データを高速にソートできることがわかります。与えられたデータが大まかにでもソートされていれば、2 番目のループで繰り返す回数が少なくなり、挿入ソートでも高速にソートすることができます。
選択ソートと挿入ソートを比較すると、データの比較回数は挿入ソートのほうが少なくなり (平均すると約半分になる)、データの移動回数は選択ソートの方が少なくなります。今回は挿入ソートの方が速くなりましたが、使用するプログラミング言語やデータの種類によっては、結果が異なる場合もあるでしょう。興味のある方はいろいろ試してみてください。
シェルソート (shell sort) は挿入ソートの改良版ともいえる方法です。最初は遠く離れた要素間でソートを開始し、徐々に間隔を狭めていきます。最後は隣り合った要素間でソートします。つまり、単純挿入ソートと同じになります。
間隔が大きいときは要素の個数が少ないので、単純なアルゴリズムでもソートにかかる時間は少なくてすみます。間隔が小さくなると要素の個数は多くなりますが、大まかにソートされているので挿入ソートでも高速にソートすることが可能です。
9 5 3 7 6 4 2 8 最初の状態
9 6 間隔を 4 で分割する
5 4
3 8
7 2
6 9 各群をソートする
4 5
3 8
2 7
6 3 9 8 間隔を 2 で分割する
4 2 5 7
3 6 8 9 各群をソートする
2 4 5 7
3 2 6 4 8 5 9 7 間隔を 1 で分割する（単純挿入ソートと同じ）
2 3 4 5 6 7 8 9 ソート完了
図 : シェルソート
プログラムは次のようになります。
リスト : シェルソート
void shell_sort(double *buff, int size)
{
for (int gap = size / 2; gap > 0; gap /= 2) {
for (int i = gap; i < size; i++) {
double temp = buff[i];
int j = i - gap;
for (; j >=0 && temp < buff[j]; j -= gap)
buff[j + gap] = buff[j];
buff[j + gap] = temp;
}
}
}
最初のループで間隔を徐々に狭めていきます。ここでは単純に 2 で割っていくことにしました。次のループで比較する要素を取り出します。最後のループでこの要素を挿入する位置を探索します。このときの探索は隣り合った要素ではなく gap 離れた要素を比較します。
2 番目のループでは、各群を並列にソートしていることに注意してください。群のひとつの要素を取り出して位置を決めたら、次の群の要素を取り出して位置を決めています。最後に gap は 1 になるので、挿入ソートと同じになりソートが完了します。
それでは実行結果を示します。
表 : shell sort の結果 (単位 : 秒) 個数 : 乱数 : 昇順 : 逆順 : 山型 --------+-------+-------+-------+------- 1000000 : 0.214 : 0.024 : 0.041 : 0.031 2000000 : 0.460 : 0.052 : 0.086 : 0.068 4000000 : 1.038 : 0.110 : 0.192 : 0.149 8000000 : 2.338 : 0.241 : 0.375 : 0.307
データ数を 100 倍していることに注意してください。結果を見ればおわかりのように、シェルソートは速いですね。gap を常に奇数になるようにすると、シェルソートの実行速度はデータの個数 n の 1.5 乗に比例します。また、Knuth 先生によると、gap の値に次の数列を用いると、シェルソートは n の 1.25 乗に比例するそうです。
gap = ..., 121, 40, 13, 4, 1
この数列は 3 倍して 1 を加えることで得られる数列を逆にしたものです。これをプログラムすると、次のようになります。
リスト : シェルソートの改良版
void shell_sort_k(double *buff, int size)
{
int gap = 1;
while (gap < size / 9) gap = gap * 3 + 1;
for (; gap > 0; gap /= 3) {
for (int i = gap; i < size; i++) {
double temp = buff[i];
int j = i - gap;
for (; j >= 0 && temp < buff[j]; j -= gap)
buff[j + gap] = buff[j];
buff[j + gap] = temp;
}
}
}
それでは実行結果を示します。
表 : shell sort (改良版) の結果 (単位 : 秒) 個数 : 乱数 : 昇順 : 逆順 : 山型 --------+-------+-------+-------+------- 1000000 : 0.182 : 0.016 : 0.022 : 0.022 2000000 : 0.395 : 0.032 : 0.052 : 0.048 4000000 : 0.874 : 0.074 : 0.108 : 0.102 8000000 : 1.930 : 0.152 : 0.232 : 0.229
少しですが速くなっていることがわかります。シェルソートは実装が簡単で、極端に要素数が大きくなければ十分実用になるソートです。なお、シェルソートは不安定なソートです。
コムソート (comb sort) はバブルソートの改良版といえる方法です。コームソートと呼ばれることもあります。今野氏の日記 (http://www.bsdclub.org/~motoyuki/d/d200005c.html#26-1) より引用します。
ソートの様子は間隔があいた櫛 (comb) ですいていくのに似ている。という感じである。 1.3 というのは数多くのデータを sort して実験的に得た数値とのこと。
- 最初に間隔が (要素数 / 1.3) の櫛で一度すく。
- 櫛の間隔を 1 / 1.3 しながら繰り返す。
- 櫛の間隔が 1 になったら bubble sort
出典は 『日経バイト 1991 年 11 月号』 とのことです。コムソートの考え方はシェルソートとよく似ています。シェルソートは挿入ソートの改良版といえる方法で、最初は遠く離れた要素間でソートを開始し、徐々に間隔を狭めていきます。そして、最後には隣り合った要素間でソートします。
シェルソートは実装が簡単で、要素数が極端に多くなければ十分実用になるソートですが、平均すればあとで説明するクイックソートよりも遅くなります。コムソートはシェルソートと同様の手法で改良を行っているようですが、それでどの程度の速さになるのかとても興味があります。
さっそくプログラムを作ってみましょう。次のリストを見てください。
リスト : コムソート
void comb_sort(double *buff, int size)
{
int gap = size;
bool done = false;
while (gap > 1 || !done) {
gap = (gap * 10) / 13;
if (!gap)
gap = 1;
else if (gap == 9 || gap == 10)
gap = 11;
done = true;
for (int i = 0; i < size - gap; i++) {
if (buff[i] > buff[i + gap]) {
double temp = buff[i];
buff[i] = buff[i + gap];
buff[i + gap] = temp;
done = false;
}
}
}
}
gap の間隔でデータを比較していきます。データの順序が逆ならば交換するところはバブルソートと同じです。gap を 1 / 1.3 ずつ狭めていき、gap が 1 になったらバブルソートと同じになります。ここで、変数 done を終了判定のフラグとして使っていることに注意してください。ここがコムソートを理解するポイントです。
バブルソートの場合、データの交換が行われないときはソートが完了しています。つまり、done が true のままであれば、ソートを終了していいわけです。櫛をすくことでデータはおおまかにソートされるため、バブルソートの途中でソートが完了する可能性は極めて高くなるでしょう。これがコムソートの原理です。
それから、間隔 gap が 9, 10 の場合には強制的に 11 にすることで、速度が向上するように改良したコムソートを comb sort 11 というそうです。今回のプログラムは comb sort 11 を使っています。
それでは実行結果を示します。
表 : comb sort の結果 (単位 : 秒) 個数 : 乱数 : 昇順 : 逆順 : 山型 --------+-------+-------+-------+------- 1000000 : 0.154 : 0.032 : 0.041 : 0.047 2000000 : 0.329 : 0.076 : 0.088 : 0.105 4000000 : 0.693 : 0.155 : 0.193 : 0.219 8000000 : 1.496 : 0.327 : 0.410 : 0.460
バブルソートや挿入ソートに比べると、コムソートの方が断然速いですね。乱数データではシェルソート (改良版) よりも速くなりました。他のデータでは、シェルソートよりも少し遅くなるようです。
クイックソートはある値を基準にして、要素をそれより大きいものと小さいものの 2 つに分割していくことでソートを行います。2 つに分けた各々の区間を同様に分割して 2 つの区間に分けます。最後は区間の要素がひとつになってソートが完了します。
9 5 3 7 6 4 2 8 最初の状態
9 5 3 7 6 4 2 8 7 を枢軸にして左側から 7 以上の値を探し、
L R 右側から 7 以下の値を探す。
2 5 3 7 6 4 9 8 交換する
L R
2 5 3 7 6 4 9 8 検索する
L R
2 5 3 4 6 7 9 8 交換する
L R
2 5 3 4 6 7 9 8 検索する。R と L が交差したら分割終了。
R L
[2 5 3 4 6] [7 9 8] この 2 つの区間について再び同様な分割を行う
図 : クイックソート
基準になる値のことを「枢軸 (pivot)」といいます。枢軸は要素の中から適当な値を選びます。今回は中央にある要素を選ぶことにしましょう。上図を見てください。左側から枢軸 7 以上の要素を探し、左側から 7 以下の要素を探します。探索のときは枢軸が番兵の役割を果たすので、ソート範囲外の要素を探索することはありません。見つけたらお互いの要素を交換します。探索位置が交差したら分割は終了です。
あとは同じ手順を分割した 2 つの区間に適用します。これは再帰定義を使えば簡単に実現できます。分割した区間の要素数が 1 になったときが再帰の停止条件になります。プログラムは次のようになります。
リスト : クイックソート (1)
void qsort_sub(double *buff, int low, int high)
{
double pivot = buff[low + (high - low) / 2];
int i = low;
int j = high;
while (true) {
while (pivot > buff[i]) i++;
while (pivot < buff[j]) j--;
if (i >= j) break;
double temp = buff[i];
buff[i] = buff[j];
buff[j] = temp;
i++;
j--;
}
if (low < i - 1) qsort_sub(buff, low, i - 1);
if (high > j + 1) qsort_sub(buff, j + 1, high);
}
void quick_sort(double *buff, int size)
{
qsort_sub(buff, 0, size - 1);
}
実際の処理は関数 qsort_sub で行います。引数 buff がソートする配列、low が区間の下限値、high が区間の上限値です。qsort_sub は buff の low から high までの区間をソートします。最初に、区間の中央にあるデータを枢軸として選びます。そして、最初の while ループで pivot を基準にして区間を 2 つに分けます。
次の while ループで、左側から枢軸以上の要素を探しています。ここでは枢軸以上という条件を、枢軸より小さい間は探索位置を進める、というように置き換えています。同様に次の while ループで右側から枢軸以下の要素を探します。お互いの探索位置 i, j が交差したら分割は終了です。break 文で while ループから脱出します。そうでなければお互いの要素を交換します。交換したあとは i と j の値を更新しておくことを忘れないでください。
そして、分割した区間に対して quick_sort を再帰呼び出しします。このとき要素数をチェックして、2 個以上ある場合に再帰呼び出しを行います。この停止条件を忘れると正常に動作しません。ご注意ください。
クイックソートは、枢軸の選び方で効率が大きく左右されます。区間の中間値を枢軸に選ぶと、区間をほぼ半分に分割することができます。この場合がいちばん効率が良く、データ数を N とすると N * log N に比例する時間でソートすることができます。
逆に、区間での最大値または最小値を枢軸に選ぶと、その要素と残りの要素の 2 つに分割にされることになります。これが最悪の場合で、分割のたびに最大値もしくは最小値を選ぶと、実行時間は要素数の 2 乗に比例することになります。つまり、挿入ソートと同じくらい遅いソートになるのです。それだけでなく、要素数が多くなるとスタックがオーバーフローする危険性もあります。
今回は区間の中央に位置する要素を枢軸としたので、中央付近に大きい要素があるデータが最悪の場合にあてはまります。つまり、山型データがこのプログラムでは最悪の結果になるのです。
それでは実行結果を示します。
表 : quick sort (1) の結果 (単位 : 秒) 個数 : 乱数 : 昇順 : 逆順 : 山型 --------+-------+-------+-------+------- 1000000 : 0.116 : 0.016 : 0.016 : ----- 2000000 : 0.235 : 0.033 : 0.035 : ----- 4000000 : 0.492 : 0.071 : 0.072 : ----- 8000000 : 1.035 : 0.145 : 0.151 : -----
山型データは 1000000 個しか試していませんが、時間がかかりすぎて途中で実行を中止しました。データ数を増やすとスタックオーバーフローするかもしれません。その他のデータでは、とても高速にソートすることができました。やっぱり、クイックソートは速いですね。なお、クイックソートは不安定なソートです。
それでは、クイックソートのプログラム (1) を改良してみましょう。まずは枢軸の選び方を工夫します。区間の中からいくつかの要素を選び、その中で真ん中の値を持つ要素を枢軸とします。たくさんの要素を選ぶとそれだけ最悪の枢軸を選ぶ危険性は減少しますが、真ん中の値を選ぶのに時間がかかってしまいます。そこで、実際には数個の要素を選んで、その中で真ん中の値を持つ要素を枢軸とする場合が多いようです。
次に再帰呼び出しをループに展開します。2 つに分割した区間の長い方を自分で用意したスタックに積み、短い区間からソートしていきます。そうすると、スタックの深さは要素数を N とすると log2 N 程度におさまります。たとえば、40 億個の要素をソートする場合でも、スタックの大きさは 32 程度あれば十分です。ひとつの区間の処理が終わったら、スタックから次の区間を取り出して処理を行います。これを図に示すと次のようになります。
スタック
[ A ][ B ] ( ) : 分割する。スタックは最初は空
(A) : A をスタックに積む
[ C ][D] (A) : B を分割する
[D] (C A) : C をスタックに積む
[ ] (C A) : D をソートする
[ C ][ ] (A) : C をスタックから取り出す
[ E ][F][ ] (A) : 分割する
[F][ ] (E A) : E をスタックに積む
[ ][ ] (E A) : F をソートする
[ E ] (A) : E をスタックより取り出す
・・・・・あとはスタックが空になるまで繰り返す・・・・・
図 : クイックソートでのスタックの動作
最後に、要素数が少なくなったらクイックソートを打ち切ります。そして、最後に挿入ソートで仕上げます。データ数が少ない場合は、クイックソートよりも単純なソートアルゴリズムの方が高速です。また、全体をソートする場合でも、すでに大まかにソートされているので、挿入ソートでも高速にソートすることができます。
まずは枢軸を選択する関数 median3 を作ります。
リスト : 枢軸の選択
int median3(int a, int b, int c)
{
if (a > b) {
if (b > c) {
return b;
} else if (a < c) {
return a;
} else {
return c;
}
} else {
if (b < c) {
return b;
} else if (a < c) {
return c;
} else {
return a;
}
}
}
引数 a, b, c の中から中央値を選んで返します。プログラムは特に難しいところはないでしょう。
ところで、枢軸の選択を工夫しても、その方法に対して最悪のデータが存在します。しかしながら、最悪の場合でも再帰呼び出しの削除とスタックの積み方の工夫により、エラーが発生することはありません。ただし、実行時間が遅くなるのは仕方がないところです。
それでは、クイックソート本体を作ります。次のリストを見てください。
リスト : クイックソート (2)
#define LIMIT 32
#define STACK_SIZE 32
void quick_sort(double *buff, int size)
{
int low_stack[STACK_SIZE];
int high_stack[STACK_SIZE];
int sp = 0, low = 0, high = size - 1;
while (true) {
if(high - low <= LIMIT){
if(!sp) break;
low = low_stack[--sp];
high = high_stack[sp];
} else {
double pivot = median3(buff[low],
buff[low + (high - low) / 2],
buff[high]);
int i = low, j = high;
while (true) {
while (pivot > buff[i]) i++;
while (pivot < buff[j]) j--;
if (i >= j) break;
double temp = buff[i];
buff[i] = buff[j];
buff[j] = temp;
i++;
j--;
}
if (i - low > high - j) {
low_stack[sp] = low;
high_stack[sp++] = i - 1;
low = j + 1;
} else {
low_stack[sp] = j + 1;
high_stack[sp++] = high;
high = i - 1;
}
}
}
insert_sort(buff, size);
}
最初に、high - low の値をチェックします。LIMIT (32) 以下になったらクイックソートを打ち切ります。スタック stack にデータがなければ break で while ループを脱出して、最後に挿入ソート (insert_sort) で仕上げます。スタックには区間の上限値と下限値がセットされています。スタックにデータがある場合は、それを取り出して low と high にセットし、処理を繰り返します。
区間が LIMIT より大きい場合はクイックソートを行います。区間の分割は今までのプログラムと同じです。枢軸は先頭の位置 low、中央の位置 (low + high) / 2、最後尾の位置 high にある要素を選ぶことにしました。median3 を呼び出して中央値を選びます。
そして、長い方の区間をスタックにセットして、短い方の区間からクイックソートします。前半部分 (i - low) が長い場合は、(low, i - 1) をスタックに積んで (j + 1, high) をソートします。low の値を j + 1 に書き換えて処理を繰り返します。後半部分が長い場合は、(j + 1, high) をスタックに積んで (low, i - 1) をソートします。この場合は high を i - 1 に書き換えて処理を繰り返します。
それでは実行結果を示します。
表 : quick sort (2) の結果 (単位 : 秒) 個数 : 乱数 : 昇順 : 逆順 : 山型 --------+-------+-------+-------+------- 1000000 : 0.101 : 0.011 : 0.012 : 0.067 2000000 : 0.208 : 0.024 : 0.025 : 0.151 4000000 : 0.431 : 0.050 : 0.053 : 0.329 8000000 : 0.916 : 0.106 : 0.117 : 0.720
非再帰版のほうが速くなりました。山型データもソートすることができますが、シェルソートやコムソートよりも遅くなりました。ちなみに、9 個の中から中央値を選ぶ median-of-9 という方法を使うと、他のデータは少し遅くなりますが、山型データはもっと速くなります。
median-of-9 は 9 つの要素から枢軸を選ぶ方法です。実際に 9 つの要素の中央値を求めているわけではありませんが、3 つの要素の中央値を求める方法 (median-of-3) よりも実行速度が速くなる場合があるようです。
枢軸を選択するプログラムは次のようになります。
リスト : 枢軸の選択
double median9(double *buff, int low, int high)
{
int m2 = (high - low) / 2;
int m4 = m2 / 2;
int m8 = m4 / 2;
int a = buff[low];
int b = buff[low + m8];
int c = buff[low + m4];
int d = buff[low + m2 - m8];
int e = buff[low + m2];
int f = buff[low + m2 + m8];
int g = buff[high - m4];
int h = buff[high - m8];
int i = buff[high];
return median3(median3(a,b,c), median3(d,e,f), median3(g,h,i));
}
関数 median9 は区間 (low, high) から 9 つの要素を選びます。区間を (0, 1) とすると、0, 1/8, 1/4, 3/8, 1/2, 5/8, 3/4, 7/8, 1 の位置にある要素を選びます。次に、9 つの要素を 3 つのグループ (0, 1/8, 1/4), (3/18, 1/2, 5/8), (3/4, 7/8, 1) に分けて、おのおののグループの中央値を median3 で求めます。さらに、その 3 つから中央値を median3 で求め、その値が枢軸となります。今回の方法は 9 つの要素の中央値を選択しているわけではありませんが、これでも十分に効果を発揮するようです。
あとは関数 quick_sort で median3 のかわりに median9 を呼び出して枢軸を選ぶだけです。それでは実行結果を示します。
表 : median-of-9 の結果 (単位 : 秒) 個数 : 乱数 : 昇順 : 逆順 : 山型 -------+-------+-------+-------+------- 1000000 : 0.105 : 0.009 : 0.010 : 0.029 2000000 : 0.215 : 0.021 : 0.022 : 0.061 4000000 : 0.458 : 0.044 : 0.046 : 0.125 8000000 : 0.942 : 0.092 : 0.101 : 0.264
乱数データでは median-of-9 のほうが少し遅くなりましたが、山型データでは median-of-9 のほうが 2 倍以上速くなりました。median-of-9 は少ないコストで最悪のケースを回避する優れた方法だと思います。もちろん、median-of-9 でも最悪のケースが存在するはずですが、最悪のケースに遭遇する確率は median-of-3 よりも低くなると思います。興味のある方はいろいろ試してみてください。
ご参考までに、ライブラリ関数 qsort の実行結果を示します。
リスト : qsort を使う場合
// 比較関数
int compare(const double *x, const double *y)
{
if (*x < *y) return -1;
if (*x > *y) return 1;
return 0;
}
void quick_sort(double *buff, int size)
{
qsort(buff, size, sizeof(double), (int(*)(const void*, const void*))compare);
}
表 : qsort の結果 (単位 : 秒) 個数 : 乱数 : 昇順 : 逆順 : 山型 --------+-------+-------+-------+------- 1000000 : 0.175 : 0.040 : 0.056 : 0.049 2000000 : 0.357 : 0.079 : 0.100 : 0.091 4000000 : 0.729 : 0.165 : 0.215 : 0.194 8000000 : 1.534 : 0.354 : 0.452 : 0.422
qsort は汎用的に作られているので、double 専用の quick sort より遅くなるのは仕方がないでしょう。
なお、これらの結果は M.Hiroi のコーディング、実行したマシン、プログラミング言語などの環境に大きく依存しています。また、これらの環境だけではなく、データの種類によっても実行時間はかなり左右されます。興味のある方は、いろいろなデータをご自分の環境で試してみてください。