GNU C Libraryのqsortはクイックソートではないの？

# はじめに

筆者の下記の記事で Microsoft の uCRT の `qsort` のソースを読んでみた。

https://qiita.com/TETSURO1999/items/ff12ba697a1ae79cbf9f

しかしながら，一般的には Free Software Foundation の GNU C library（glibc）のほうがメジャーな存在かもしれない。

# glibc の qsort をオンライン環境で実行してみる

わざわざ GCC をインストールしなくてもオンライン環境で検証できるのは助かる[^1]。今回は下記のサイトを使ってみた。

[wandbox.org](https://wandbox.org/)

テストプログラムを下記に示す。[以前の記事](https://qiita.com/TETSURO1999/items/68422ff100dbecf2ff1a)で使ったものと同じである。
※日本語のコメントは削除した。

```c++:test.c（テストプログラム）
#include <stdio.h>
#include <stdlib.h>
#define countof(a)  (sizeof(a)/sizeof(a[0]))
typedef	int (*SORTFUNC)(const void *, const void *);
static  int table[] = {9, 15, 16, 14, 11, 12,  7, 13, 10, 0, 8, 5, 6, 1, 4, 3, 2};
static	int	count = 0;
static	int	sortfunc(int *p, int *q) {
	count++;
	fprintf(stdout, "%d\t%d\t%d\n", *p, *q, *p - *q);
	return *p - *q;
}
int main() {
    qsort(table, countof(table), sizeof(table[0]), (SORTFUNC)sortfunc);
    fprintf(stdout, "%d\n", count);
    return 0;
}
```

実行結果を以下に示す。驚くべきことに比較関数の呼び出し回数が46回と uCRT の78回よりも少ない。最初の三回の比較要素を比べても，uCRT では先頭と末尾，中央の三点でピボットとなるべき基準値をサンプリングしていたが，glibc では明らかに様相が異なる。

```plantuml
@startuml
<style>
root {
  margin 0
  padding 0
}
title {
  padding 0 16 16 16
  roundcorner 16
  horizontalalignment left
  backgroundcolor #202020
  fontcolor #e3e3e3
  fontstyle normal
  fontname monospaced
  fontsize 14
  backgroundcolor #202020
}
</style>
sprite TITLE <svg width="640" height="21">
<path d="M0 0 v24 a4 4 0 0 0 4 4 h260 a 4 4 0 0 0 4 -4 v-24" fill="#606060" />
<text x="8" y="20" font-size="14" fill="#e3e3e3">test.c（テストプログラム）の実行結果</text>
</svg>
title
<$TITLE>

9	15	-6	<color:#3bbfbf>#  0 番目と  1 番目の比較</color>
16	14	2	<color:#3bbfbf>#  2 番目と  3 番目の比較</color>
9	14	-5
15	14	1
15	16	-1
11	12	-1	<color:#3bbfbf>#  4 番目と  5 番目の比較</color>
7	13	-6	<color:#3bbfbf>#  6 番目と  7 番目の比較</color>
11	7	4
11	13	-2
12	13	-1
9	7	2
9	11	-2
14	11	3
14	12	2
14	13	1
10	0	10	<color:#3bbfbf>#  8 番目と  9 番目の比較</color>
8	5	3	<color:#3bbfbf># 10 番目と 11 番目の比較</color>
0	5	-5
10	5	5
10	8	2
6	1	5	<color:#3bbfbf># 12 番目と 13 番目の比較</color>
3	2	1	<color:#3bbfbf># 15 番目と 16 番目の比較</color>
4	2	2
4	3	1
1	2	-1
6	2	4
6	3	3
6	4	2
0	1	-1
5	1	4
5	2	3
5	3	2
5	4	1
5	6	-1
8	6	2
7	0	7
7	1	6
7	2	5
7	3	4
7	4	3
7	5	2
7	6	1
7	8	-1
9	8	1
9	10	-1
11	10	1
46
end title
@enduml
```

比較関数の呼び出し回数を比べると，.NET Framework の106回の半分以下である。

<table>
<caption>表１　比較関数の呼び出し回数</caption>
<thead>
<tr><th style="text-align:center">ライブラリ</th>
    <th style="text-align:center">回数</th></tr>
</thead>
<tbody>
<tr></tr>
<tr><th>Array.Sort (.NET Framework)</th>
    <td style="text-align:right">106</td></tr>
<tr></tr>
<tr><th>qsort (Microsoft uCRT)</th>
    <td style="text-align:right">78</td></tr>
<tr></tr>
<tr><th>qsort (GNU C Library)</th>
    <td style="text-align:right">46</td></tr>
</tbody>
</table>

結論から言うと glibc の `qsort` はマージソートのようだ[^2]。

# glibc の qsort の意訳版を作る

glibc の qsort.c のソースコード[^3]を読むとはっきりする。とはいえ，全文を掲載しても理解を妨げるので，アルゴリズムを自分なりに咀嚼し，整数配列のソートとして新たに作り直したものを掲載する。大して長くは無いのでまずは全文を掲載しよう。

```c++:int 型配列のソートに作り直したマージソート
/* 001 */ #include <stdlib.h>
/* 002 */ typedef int     (*COMPARE)(int, int);
/* 003 */ #define TEMP_SIZE   256
/* 004 */ static  void    sort(int *buf, int *a, int len, COMPARE compare) {
/* 005 */     if(len < 2) return;
/* 006 */     int n1 = len / 2;
/* 007 */     int n2 = len - n1;
/* 008 */     int *p1 = a;
/* 009 */     int *p2 = &a[n1];
/* 010 */     sort(buf, p1, n1, compare);
/* 011 */     sort(buf, p2, n2, compare);
/* 012 */     int *q = buf;
/* 013 */     for(;;) {
/* 014 */         if(*p1 <= *p2) {
/* 015 */             *q++ = *p1++;
/* 016 */             if(--n1 == 0) break;
/* 017 */         } else {
/* 018 */             *q++ = *p2++;
/* 019 */             if(--n2 == 0) break;
/* 020 */         }
/* 021 */     }
/* 022 */     while(n1 > 0) {
/* 023 */         *q++ = *p1++;
/* 024 */         n1--;
/* 025 */     }
/* 026 */     for(int i = 0; i < len - n2; i++)
/* 027 */         a[i] = buf[i];
/* 028 */ }
/* 029 */ void    qsort_glibc(int *a, int len, COMPARE compare) {
/* 030 */     if(len < 2) return;
/* 031 */     if(len <= TEMP_SIZE) {
/* 032 */         int tmp[TEMP_SIZE];
/* 033 */         sort(tmp, a, len, compare);
/* 034 */     } else {
/* 035 */         int *buf = (int*)malloc(sizeof(int) * len);
/* 036 */         if(buf == (int*)NULL) return;
/* 037 */         sort(buf, a, len, compare);
/* 038 */         free(buf);
/* 039 */     }
/* 040 */ }
```

# アルゴリズム解説

それでは意訳コードを頭から解説していきたい。

## 1. ソート関数本体の冒頭部分

ソート関数本体 `qsort_glibc()` の冒頭部分である。引数は `int` 型の配列 `a[]` と配列の要素数 `len`，比較関数 `compare` の三つのみである。比較関数の型は標準ライブラリのものと少し変えており，二つの `int` 型の整数を直接比較する。

マージソートはソート対象の配列の他に作業用の記憶領域を必要とする外部ソートである。配列の要素数 `len` が所定の大きさ `TEMP_SIZE` 以下であれば作業用領域をスタック上に確保するが，`TEMP_SIZE` より大きければヒープ上に確保する。


なお，本自作コードではメモリの確保に失敗すると何もしないでリターンするが，本物はヒープソート（ヒープ領域のヒープとは無関係で，二分ヒープ木を用いて並べ替えを行うソート）を用いてリトライするよう作られている。本実装では簡単のためヒープソートは省略した。

```c++:ソート関数本体の冒頭部分
/* 001 */ #include <stdlib.h>
/* 002 */ typedef int     (*COMPARE)(int, int);
/* 003 */ #define TEMP_SIZE   256

/* 029 */ void    qsort_glibc(int *a, int len, COMPARE compare) {
/* 030 */     if(len < 2) return;
/* 031 */     if(len <= TEMP_SIZE) {
/* 032 */         int tmp[TEMP_SIZE];
/* 033 */         sort(tmp, a, len, compare);
/* 034 */     } else {
/* 035 */         int *buf = (int*)malloc(sizeof(int) * len);
/* 036 */         if(buf == (int*)NULL) return;
/* 037 */         sort(buf, a, len, compare);
/* 038 */         free(buf);
/* 039 */     }
/* 040 */ }
```

## 2. マージソートエンジン・冒頭部

次はマージソートエンジンの冒頭部である。引数は作業用領域 `buf[]` とソート対象の配列 `a[]` と配列の要素数 `len`，比較関数 `compare` の四つである。作業用領域は十分確保されているのでオーバーフローを気にしなくてよい。

配列の要素数 `len` が 2 個未満のときは何もしないでリターンする。2 個以上の場合は，配列を二つに分けてそれぞれにマージソートを行う。すなわちマージソートエンジンは再帰関数として作られている。

```c++:マージソートエンジン・冒頭部
/* 004 */ static  void    sort(int *buf, int *a, int len, COMPARE compare) {
/* 005 */     if(len < 2) return;
/* 006 */     int n1 = len / 2;
/* 007 */     int n2 = len - n1;
/* 008 */     int *p1 = a;
/* 009 */     int *p2 = &a[n1];
/* 010 */     sort(buf, p1, n1, compare);
/* 011 */     sort(buf, p2, n2, compare);
```

## 3. マージソートエンジン・マージ部

二つの配列を作業用領域 `buf[]` 上でマージした後，元の配列 `a[]` に書き戻す部分である。ここの22～27行でちょっとした工夫が行われている。

```c++:マージソートエンジン・マージ部
/* 012 */     int *q = buf;
/* 013 */     for(;;) {
/* 014 */         if(*p1 <= *p2) {
/* 015 */             *q++ = *p1++;
/* 016 */             if(--n1 == 0) break;
/* 017 */         } else {
/* 018 */             *q++ = *p2++;
/* 019 */             if(--n2 == 0) break;
/* 020 */         }
/* 021 */     }
/* 022 */     while(n1 > 0) {
/* 023 */         *q++ = *p1++;
/* 024 */         n1--;
/* 025 */     }
/* 026 */     for(int i = 0; i < len - n2; i++)
/* 027 */         a[i] = buf[i];
/* 028 */ }
```

上記22～27行の部分を教科書通りに作れば下記のようになるはずだ。つまり，上記の22～27行は転送領域が不足しているように見えるのだ。

```c++
while(n1 > 0) {
    *q++ = *p1++; n1--;
}
while(n2 > 0) {
    *q++ = *p2++; n2--;
}
for(int i = 0; i < len; i++)
    a[i] = buf[i];
```

これは簡単な図を描いてみると分かる。ソート対象の配列 `a[]` を前半部と後半部に分け，マージしながら作業用領域 `buf[]` に転送していく。21行目終了時点で転送済みの領域を<font color=lightgray>■</font>で示す。

21行目終了時点で `n1 > 0` のときの場合を図１に示す。このとき残った領域もいったん作業用領域 `buf[]` の残った領域に転送してから配列 `a[]` に書き戻す必要がある。

```plantuml
@startuml
<style>
root {
  margin 0
  padding 0
}
title {
  fontStyle normal
  fontSize 16
}
</style>
sprite FIGURE <svg width="640" height="160">
<path d="M32  20 h272 v1 h-272" fill="black" />
<path d="M32  20 h10 v1 h-10" fill="black" transform="rotate(20,32,20)"/>
<path d="M32  20 h10 v1 h-10" fill="black" transform="rotate(-20,32,20)"/>
<path d="M608 20 h-272 v1 h272" fill="black" />
<path d="M608 20 h-10 v1 h10" fill="black" transform="rotate(20,608,20)"/>
<path d="M608 20 h-10 v1 h10" fill="black" transform="rotate(-20,608,20)"/>
<text x="310" y="24" font-size="14" fill="black">len</text>
<path d="M256 76 h20 v1 h-20" fill="black" />
<path d="M256 76 h10 v1 h-10" fill="black" transform="rotate(20,256,76)"/>
<path d="M256 76 h10 v1 h-10" fill="black" transform="rotate(-20,256,76)"/>
<path d="M320 76 h-20 v1 h20" fill="black" />
<path d="M320 76 h-10 v1 h10" fill="black" transform="rotate(20,320,76)"/>
<path d="M320 76 h-10 v1 h10" fill="black" transform="rotate(-20,320,76)"/>
<text x="278" y="82" font-size="14" fill="black">n1</text>
<text x="12" y="53" font-size="14" fill="black">a</text>
<path d="M32  32 v32 h224 v-32" fill="lightgray" />
<path d="M320 32 v32 h288 v-32" fill="lightgray" />
<path d="M32  32 h576 v1 h-576" fill="black" />
<path d="M32  64 h576 v1 h-576" fill="black" />
<path d="M32  32 v32 h1 v-32" fill="black" />
<path d="M320 32 v32 h1 v-32" fill="black" />
<path d="M608 32 v32 h1 v-32" fill="black" />
<text x="4" y="117" font-size="14" fill="black">buf</text>
<path d="M32  96 v32 h512 v-32" fill="lightgray" />
<path d="M32  96 h576 v1 h-576" fill="black" />
<path d="M32 128 h576 v1 h-576" fill="black" />
<path d="M32  96 v32 h1 v-32" fill="black" />
<path d="M608 96 v32 h1 v-32" fill="black" />
<path d="M544 140 h20 v1 h-20" fill="black" />
<path d="M544 140 h10 v1 h-10" fill="black" transform="rotate(20,544,140)"/>
<path d="M544 140 h10 v1 h-10" fill="black" transform="rotate(-20,544,140)"/>
<path d="M608 140 h-20 v1 h20" fill="black" />
<path d="M608 140 h-10 v1 h10" fill="black" transform="rotate(20,608,140)"/>
<path d="M608 140 h-10 v1 h10" fill="black" transform="rotate(-20,608,140)"/>
<text x="566" y="146" font-size="14" fill="black">n1</text>
</svg>
title
<$FIGURE>
図１　21行目終了時点で n1 > 0 の場合
end title
@enduml
```

一方，21行目終了時点で `n2 > 0` のときの場合を図２に示す。このとき残った領域はそのまま残してもよいことが分かる。

```plantuml
@startuml
<style>
root {
  margin 0
  padding 0
}
title {
  fontStyle normal
  fontSize 16
}
</style>
sprite FIGURE <svg width="640" height="160">
<path d="M32  20 h272 v1 h-272" fill="black" />
<path d="M32  20 h10 v1 h-10" fill="black" transform="rotate(20,32,20)"/>
<path d="M32  20 h10 v1 h-10" fill="black" transform="rotate(-20,32,20)"/>
<path d="M608 20 h-272 v1 h272" fill="black" />
<path d="M608 20 h-10 v1 h10" fill="black" transform="rotate(20,608,20)"/>
<path d="M608 20 h-10 v1 h10" fill="black" transform="rotate(-20,608,20)"/>
<text x="310" y="24" font-size="14" fill="black">len</text>
<path d="M544 76 h20 v1 h-20" fill="black" />
<path d="M544 76 h10 v1 h-10" fill="black" transform="rotate(20,544,76)"/>
<path d="M544 76 h10 v1 h-10" fill="black" transform="rotate(-20,544,76)"/>
<path d="M608 76 h-20 v1 h20" fill="black" />
<path d="M608 76 h-10 v1 h10" fill="black" transform="rotate(20,608,76)"/>
<path d="M608 76 h-10 v1 h10" fill="black" transform="rotate(-20,608,76)"/>
<text x="566" y="82" font-size="14" fill="black">n2</text>
<text x="12" y="53" font-size="14" fill="black">a</text>
<path d="M32  32 v32 h512 v-32" fill="lightgray" />
<path d="M32  32 h576 v1 h-576" fill="black" />
<path d="M32  64 h576 v1 h-576" fill="black" />
<path d="M32  32 v32 h1 v-32" fill="black" />
<path d="M320 32 v32 h1 v-32" fill="black" />
<path d="M608 32 v32 h1 v-32" fill="black" />
<text x="4" y="117" font-size="14" fill="black">buf</text>
<path d="M32  96 v32 h512 v-32" fill="lightgray" />
<path d="M32  96 h576 v1 h-576" fill="black" />
<path d="M32 128 h576 v1 h-576" fill="black" />
<path d="M32  96 v32 h1 v-32" fill="black" />
<path d="M608 96 v32 h1 v-32" fill="black" />
<path d="M544 140 h20 v1 h-20" fill="black" />
<path d="M544 140 h10 v1 h-10" fill="black" transform="rotate(20,544,140)"/>
<path d="M544 140 h10 v1 h-10" fill="black" transform="rotate(-20,544,140)"/>
<path d="M608 140 h-20 v1 h20" fill="black" />
<path d="M608 140 h-10 v1 h10" fill="black" transform="rotate(20,608,140)"/>
<path d="M608 140 h-10 v1 h10" fill="black" transform="rotate(-20,608,140)"/>
<text x="566" y="146" font-size="14" fill="black">n2</text>
</svg>
title
<$FIGURE>
図２　21行目終了時点で n2 > 0 の場合
end title
@enduml
```

ということで，転送量を少しでも減らせるよう細かい工夫が行われているようだ。

また，配列 `a[]` を前半部と後半部に分割する際，配列の要素数が奇数だと等分割できない。非均等に分割するのであれば，後半部の要素数を多くした方が得だと分かるだろう。これは後半部の要素数のほうが多いと，後半部の要素が残ってしまう確率が高くなるが，その場合，残った要素を転送しなくても済むからだ。実際，配列の要素数 `len` が奇数の場合，前半部の要素数 `n1` よりも後半部の要素数 `n2` がのほうが一つだけ多くなるよう分割される。

```c++:配列の前後分割
/* 006 */     int n1 = len / 2;
/* 007 */     int n2 = len - n1;
```

# パフォーマンス比較

さて，前回の記事で紹介した uCRT を参考にして作成した `qsort_ucrt()` と本記事で紹介した glibc を参考にして作成した `qsort_glibc()` ではどっちが速いだろうか？

試験条件を以下に示す。

- 200,000,000 個の 32bit 整数を昇順ソートする
- データの順序は以下の４種類
  - ゼロ
  - 昇順：0 → 199,999,999
  - 降順：199,999,999 → 0
  - 乱数：0～199,999,999 の範囲で一様乱数（重複あり）
- 測定環境：Windows10 Pro 22H2，Core i7-13700，5.1GHz，32GB
- 開発環境：Visual Studio 2022 Community Edition，32bit版コンパイラ
- 比較関数は以下に示すように整数の昇順ソートである。
  ```c++
  static	int		compare(int p, int q) {
	return p - q;
  }
  ```

ソート方式の比較結果を以下に示す。あれれ？期待に反してマージソートの `qsort_glibc()` のほうが遅いではないか！！

<table>
<caption>表２　ソート方式の比較</caption>
<thead>
<tr><th style="text-align:center">順序</th>
    <th style="text-align:center">項目</th>
    <th style="text-align:center"><code>qsort_ucrt()</code></th>
    <th style="text-align:center"><code>qsort_glibc()</code></th></tr>
</thead>
<tbody>
<tr></tr>
<tr><th rowspan=5>ゼロ</th>
    <th>時間</th>
    <td style="text-align:right">0.090[秒]</td>
    <td style="text-align:right">2.660[秒]</td></tr>
<tr></tr>
<tr><th>比較回数</th>
    <td style="text-align:right">399,999,999[回]</td>
    <td style="text-align:right">2,728,894,208[回]</td></tr>
<tr></tr>
<tr><th>転送回数</th>
    <td style="text-align:right">0[回]</td>
    <td style="text-align:right">5,457,788,416[回]</td></tr>
<tr></tr>
<tr><th rowspan=5>昇順</th>
    <th>時間</th>
    <td style="text-align:right">1.433[秒]</td>
    <td style="text-align:right">2.656[秒]</td></tr>
<tr></tr>
<tr><th>比較回数</th>
    <td style="text-align:right">5,363,792,412[回]</td>
    <td style="text-align:right">2,728,894,208[回]</td></tr>
<tr></tr>
<tr><th>転送回数</th>
    <td style="text-align:right">265,782,274[回]</td>
    <td style="text-align:right">5,457,788,416[回]</td></tr>
<tr></tr>
<tr><th rowspan=5>降順</th>
    <th>時間</th>
    <td style="text-align:right">1.475[秒]</td>
    <td style="text-align:right">3.914[秒]</td></tr>
<tr></tr>
<tr><th>比較回数</th>
    <td style="text-align:right">5,363,792,412[回]</td>
    <td style="text-align:right">2,802,670,336[回]</td></tr>
<tr></tr>
<tr><th>転送回数</th>
    <td style="text-align:right">465,782,278[回]</td>
    <td style="text-align:right">11,063,129,088[回]</td></tr>
<tr></tr>
<tr><th rowspan=5>乱数</th>
    <th>時間</th>
    <td style="text-align:right">15.505[秒]</td>
    <td style="text-align:right">17.996[秒]</td></tr>
<tr></tr>
<tr><th>比較回数</th>
    <td style="text-align:right">6,400,427,154[回]</td>
    <td style="text-align:right">5,265,836,886[回]</td></tr>
<tr></tr>
<tr><th>転送回数</th>
    <td style="text-align:right">2,736,238,624[回]</td>
    <td style="text-align:right">10,752,119,342[回]</td></tr>
</tbody>
</table>


<font color=#5B9BD5>■</font>qsort_ucrt，<font color=#ED7D31>■</font>qsort_glibc</font>

```plantuml
@startuml
<style>
root {
  margin 0
  padding 0
}
title {
  fontStyle normal
  fontSize 17
}
</style>
sprite FIGURE <svg width="560" height="540">
<g fill="lightgray">
<path d="M60 20 h480 v1 h-480">
<path d="M60 140 h480 v1 h-480">
<path d="M60 260 h480 v1 h-480">
<path d="M60 380 h480 v1 h-480">
<path d="M60 500 h480 v1 h-480">
</g>
<text x="45" y="506" font-size="16" fill="black">0</text>
<text x="45" y="386" font-size="16" fill="black">5</text>
<text x="35" y="266" font-size="16" fill="black">10</text>
<text x="35" y="146" font-size="16" fill="black">15</text>
<text x="35" y="26"  font-size="16" fill="black">20</text>
<text x="105" y="520" font-size="16" fill="black">ゼロ</text>
<text x="225" y="520" font-size="16" fill="black">昇順</text>
<text x="345" y="520" font-size="16" fill="black">降順</text>
<text x="465" y="520" font-size="16" fill="black">乱数</text>
<g fill="#5B9BD5">
<g fill="#5B9BD5">
<path d="M76 500 h40 v-2 h-40">
<path d="M196 500 h40 v-34 h-40">
<path d="M316 500 h40 v-35 h-40">
<path d="M436 500 h40 v-372 h-40">
</g>
<g fill="#ED7D31">
<path d="M124 500 h40 v-64 h-40">
<path d="M244 500 h40 v-64 h-40">
<path d="M364 500 h40 v-94 h-40">
<path d="M484 500 h40 v-432 h-40">
</g>
</svg>
title
<$FIGURE>
図３　ソート時間の比較（単位は秒）
end title
@enduml
```

ということで比較回数・転送回数も比べてみた。なお `qsort_ucrt()` では１回のスワップあたり，データ転送２回分とカウントしている。比較回数はマージソートである `qsort_glibc()` のほうが安定して少ないのは流石であるが，転送回数は圧倒的に多い。これが基本的にマージソートの遅い理由ではある。しかし，同じ `qsort_glibc()` の降順と乱数で転送回数は大差ないにも関わらず，計算時間は乱数のほうが４倍近く要していることに注目されたい。降順のほうが圧倒的に速いのは，おそらくメモリアクセスがシーケンシャルアクセスに近いせいでキャッシュメモリに乗り易いからだろうと考える。一方，乱数のほうはメモリアクセスがランダムアクセスに近いせいでキャッシュミスヒットが連続しているのだろう。

```plantuml
@startuml
<style>
root {
  margin 0
  padding 0
}
title {
  fontStyle normal
  fontSize 17
}
</style>
sprite FIGURE <svg width="340" height="530">
<path d="M40 490 h280 v1 h-280" fill="lightgray" />
<path d="M40 410 h280 v1 h-280" fill="lightgray" />
<path d="M40 330 h280 v1 h-280" fill="lightgray" />
<path d="M40 250 h280 v1 h-280" fill="lightgray" />
<path d="M40 170 h280 v1 h-280" fill="lightgray" />
<path d="M40  90 h280 v1 h-280" fill="lightgray" />
<path d="M40  10 h280 v1 h-280" fill="lightgray" />

<path d="M55 490 h20 v-16 h-20" fill="#5B9BD5">
<path d="M125 490 h20 v-215 h-20" fill="#5B9BD5">
<path d="M195 490 h20 v-215 h-20" fill="#5B9BD5">
<path d="M265 490 h20 v-256 h-20" fill="#5B9BD5">
<path d="M79 490 h20 v-109 h-20" fill="#ED7D31">
<path d="M149 490 h20 v-109 h-20" fill="#ED7D31">
<path d="M219 490 h20 v-112 h-20" fill="#ED7D31">
<path d="M289 490 h20 v-211 h-20" fill="#ED7D31">

<text x="25" y="496" font-size="16" fill="black">0</text>
<text x="15" y="416" font-size="16" fill="black">20</text>
<text x="15" y="336" font-size="16" fill="black">40</text>
<text x="15" y="256" font-size="16" fill="black">60</text>
<text x="15" y="176" font-size="16" fill="black">80</text>
<text x="5" y="96"  font-size="16" fill="black">100</text>
<text x="5" y="16"  font-size="16" fill="black">120</text>
<text x="60"  y="510" font-size="16" fill="black">ゼロ</text>
<text x="130" y="510" font-size="16" fill="black">昇順</text>
<text x="200" y="510" font-size="16" fill="black">降順</text>
<text x="270" y="510" font-size="16" fill="black">乱数</text>
</svg>
title
<$FIGURE>
図４　比較回数の比較（単位は億回）
end title
@enduml
```

```plantuml
@startuml
<style>
root {
  margin 0
  padding 0
}
title {
  fontStyle normal
  fontSize 17
}
</style>
sprite FIGURE <svg width="340" height="530">
<path d="M40 490 h280 v1 h-280" fill="lightgray" />
<path d="M40 410 h280 v1 h-280" fill="lightgray" />
<path d="M40 330 h280 v1 h-280" fill="lightgray" />
<path d="M40 250 h280 v1 h-280" fill="lightgray" />
<path d="M40 170 h280 v1 h-280" fill="lightgray" />
<path d="M40  90 h280 v1 h-280" fill="lightgray" />
<path d="M40  10 h280 v1 h-280" fill="lightgray" />

<path d="M55 490 h20 v-0 h-20" fill="#5B9BD5">
<path d="M125 490 h20 v-11 h-20" fill="#5B9BD5">
<path d="M195 490 h20 v-19 h-20" fill="#5B9BD5">
<path d="M265 490 h20 v-109 h-20" fill="#5B9BD5">
<path d="M79 490 h20 v-218 h-20" fill="#ED7D31">
<path d="M149 490 h20 v-218 h-20" fill="#ED7D31">
<path d="M219 490 h20 v-443 h-20" fill="#ED7D31">
<path d="M289 490 h20 v-430 h-20" fill="#ED7D31">

<text x="25" y="496" font-size="16" fill="black">0</text>
<text x="15" y="416" font-size="16" fill="black">20</text>
<text x="15" y="336" font-size="16" fill="black">40</text>
<text x="15" y="256" font-size="16" fill="black">60</text>
<text x="15" y="176" font-size="16" fill="black">80</text>
<text x="5" y="96"  font-size="16" fill="black">100</text>
<text x="5" y="16"  font-size="16" fill="black">120</text>
<text x="60"  y="510" font-size="16" fill="black">ゼロ</text>
<text x="130" y="510" font-size="16" fill="black">昇順</text>
<text x="200" y="510" font-size="16" fill="black">降順</text>
<text x="270" y="510" font-size="16" fill="black">乱数</text>
</svg>
title
<$FIGURE>
図５　転送回数の比較（単位は億回）
end title
@enduml
```

# 考察

`qsort_ucrt()` に比べて `qsort_glibc()` が振るわなかった理由を考察する。

1. `qsort_glibc()` はマージソートであり，そもそもデータ転送量が多いため，キャッシュメモリに収まり切らない大規模な問題では性能が出ないのかもしれない。

2. 今回のプログラムはいずれも比較関数を簡素な整数比較にしたことに加え，ソート関数と比較関数を一緒にコンパイルしたことにより，最適化によって比較関数がソート関数内にインライン展開され，比較関数を呼び出すオーバーヘッドが軽減された。このせいで比較関数の呼び出し回数の多いクイックソートである `qsort_ucrt()` に有利な条件になっている。

3. `qsort_glibc()` は普通に再帰関数として作成しているが，`qsort_ucrt()` は再帰関数ではなくループ処理で済ませており，再帰関数の呼び出しオーバーヘッドを減らすなど細かい工夫が施されている。

4. `qsort_glibc()` は配列の要素数が 2 個になってもマージソートを行っているが，`qsort_ucrt()` は分割された配列の要素数が 8 個以下になると選択ソートに切り替えるなど細かい工夫が施されている。

# オリジナル作者の意見

なぜ，glibc の `qsort()` ではマージソートを採用したのか？という疑問に対して，オリジナル作者の mikehaertel 氏が Reddit で次のように述べている[^4]。

> I am the (original) author of the merge sort code. Someone else wrote the quick sort code. The reason we made qsort() default to using merge sort was, indeed, that it always makes fewer compares than quick sort. We thought that in the common case, people would sort small objects (numbers, pointers) rather than large objects, and that therefore the performance would be dominated by comparisons rather than data movement.

> 私はマージソートのコードの（オリジナルの）作者です。クイックソートのコードは他の人が書きました。qsort() のデフォルトをマージソートにしたのはクイックソートよりも比較回数を少なくするためです。一般的なケースでは人々は大きなオブジェクトではなく，数値やポインタのような小さなオブジェクトをソートするでしょうから，データ転送よりも比較のほうが性能上重要だと考えたのです。

> In 1988, the dominant cost was usually function call overhead for calling the compare function. Nowdays, due to out-of-order processors with multiported caches and call/return branch predictors, the function call overhead is less noticeable, and has probably been replaced by branch misprediction costs. Assuming random data, you'd expect the compare to go either way about half the time. Many compare functions have branches, and of course so does the sort code itself. So, for different reasons than in 1988, performance is still usually dominated by the number of compares.

> 1988年当時，比較関数を呼び出すための関数コールのオーバーヘッドが主要なコストでした。現在ではマルチポートキャッシュとコール／リターン分岐予測を備えたアウトオブオーダーのプロセッサにより関数コールのオーバーヘッドは目立たなくなっており，おそらくは分岐予測ミスのコストが代わりに支配的になっています。ランダムなデータを仮定すると，比較は約半分の確率でどちらに転ぶか分かりません。多くの比較関数には分岐があり，もちろんソートのコード自体にも分岐があります。そのため1988年当時とは異なった理由で今でも比較回数が性能上重要なのです。

1988年といえば，当時先端の EWS ですら CPU クロックはせいぜい 25～33MHz 程度。メインメモリアクセスに対するウェイトは小さく，データ転送の負荷は少なかったものと推察される。なので1988年当時はマージソートを採用するのが正しかったかもしれない。

とはいえ，Reddit で上記の回答が行われたのは今から18年前なので，上記の回答における「現在」や「今」とは2007年のことである。2007年当時の CPU は既にクロック3GHzに達しており，1988年に比べてデータ転送の負荷は飛躍的に増えているはずなのだ。おそらく2007年の時点では，もうマージソートの優位性は失われていたと思う。

# 付録～ソースコードおよびメイクファイル

<details><summary><b>ベンチ用メイン</b></summary>

```c++:bench.c
#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <string.h>
#include <time.h>
#define	MACROOF(arg)	#arg
#define	STRINGOF(macro)	MACROOF(macro)
typedef	int		(*COMPARE)(int, int);
extern	void	SORTFUNC(int *a, int len, COMPARE compare);
uint32_t	xorshift( void ) {
static	uint64_t	x = 1;
	x ^= x << 13;
	x ^= x >>  7;
	x ^= x << 17;
	return (uint32_t)x;
}
#ifdef	COUNT
uint64_t	move_count = 0;
uint64_t	comp_count = 0;
#endif
static	int		compare(int p, int q) {
#ifdef	COUNT
	comp_count++;
#endif
	return p - q;
}
int	main(int argc, char *argv[]) {
	if(argc < 3) {
		fprintf(stderr, "Usage: bench(.exe) [data] [size] (debug)\n");
		fprintf(stderr, " data: zero, ascend, descend, random\n");
		fprintf(stderr, " func: " STRINGOF(SORTFUNC) "\n");
		return -1;
	}
	int	size = atoi(argv[2]);
	int	*table = malloc(sizeof(int) * size);
	if(table == (int*)NULL) {
		fprintf(stderr, "cannot allocate size: %d\n", size);
		return -1;
	}
	switch(argv[1][0]) {
	  case 'z': for(int i = 0; i < size; i++) table[i] = 0;					break;
	  case 'a':	for(int i = 0; i < size; i++) table[i] = i;					break;
	  case 'd':	for(int i = 0; i < size; i++) table[i] = size - i - 1;		break;
	  case 'r':	for(int i = 0; i < size; i++) table[i] = xorshift() % size;	break;
	  default:
		fprintf(stderr, "illegal data: %s\n", argv[1]);
		return -1;
	}
#ifdef	COUNT
	SORTFUNC(table, size, compare);
	fprintf(stdout, "comp_count: %lld\n", comp_count);
	fprintf(stdout, "move_count: %lld\n", move_count);
#else
	clock_t	t1 = clock();
	SORTFUNC(table, size, compare);
	clock_t	t2 = clock();
	fprintf(stdout, "time: %.3f\n", (double)(t2 - t1) / CLOCKS_PER_SEC);
#endif
	if(argc >= 4 && argv[3][0] == 'd')
		for(int i = 0; i < size; i++) fprintf(stdout, "%d\n", table[i]);
	return 0;
}
```

</details>

<details><summary><b>uCRT 版 qsort の意訳版</b></summary>

```c++:qsort_ucrt.c
#ifdef	COUNT
#include <stdint.h>
extern	uint64_t	move_count;
#endif
typedef	int		(*COMPARE)(int, int);
#define	STACK_SIZE	30
static	void	swap(int *p, int *q) {
#ifdef	COUNT
	move_count += 2;
#endif
	int	temp = *p;
	*p = *q;
	*q = temp;
}
void	qsort_ucrt(int *a, int len, COMPARE compare) {
	if(len < 2) return;
	int	top_stack[STACK_SIZE];
	int	end_stack[STACK_SIZE];
	int	stack = 0;
	int	top = 0;
	int	end = len - 1;
LOOP:
	int	n = end - top + 1;
	if(n <= 8) {
		for(int e = end; e > top; e--) {
			int	max = top;
			for(int i = top + 1; i <= e; i++)
				if(compare(a[i], a[max]) > 0) max = i;
			swap(&a[max], &a[e]);
		}
	} else {
		int	mid = top + n / 2;
		int	t = top;
		int	e = end;
		if(compare(a[top], a[mid]) > 0) swap(&a[top], &a[mid]);
		if(compare(a[top], a[end]) > 0) swap(&a[top], &a[end]);
		if(compare(a[mid], a[end]) > 0) swap(&a[mid], &a[end]);
		for(;;) {
			if(t <  mid) while(++t <  mid && compare(a[t], a[mid]) <= 0);
			if(t >= mid) while(++t <= end && compare(a[t], a[mid]) <= 0);
			while(--e > mid && compare(a[e], a[mid]) > 0);
			if(t > e) break;
			swap(&a[t], &a[e]);
			if(e == mid) mid = t;
		}
		e++;
		if(e >  mid) while(--e > mid && compare(a[e], a[mid]) == 0);
		if(e <= mid) while(--e > top && compare(a[e], a[mid]) == 0);
		if(e - top >= end - t) {
			if(top < e) {
				top_stack[stack] = top;
				end_stack[stack] = e;
				stack++;
			}
			if(t < end) {
				top = t;
				goto LOOP;
			}
		} else {
			if(t < end) {
				top_stack[stack] = t;
				end_stack[stack] = end;
				stack++;
			}
			if(top < e) {
				end = e;
				goto LOOP;
			}
		}
	}
	if(--stack >= 0) {
		top = top_stack[stack];
		end = end_stack[stack];
		goto LOOP;
	}
}
```

</details>

<details><summary><b>glibc 版 qsort の意訳版</b></summary>

```c++:qsort_glib.c
#ifdef	COUNT
#include <stdint.h>
extern	uint64_t	move_count;
#endif
#include <stdlib.h>
typedef	int		(*COMPARE)(int, int);
#define	TEMP_SIZE	256
static	void	sort(int *buf, int *a, int len, COMPARE compare) {
	if(len < 2) return;
	int	n1 = len / 2;
	int	n2 = len - n1;
	int	*p1 = a;
	int	*p2 = &a[n1];
	sort(buf, p1, n1, compare);
	sort(buf, p2, n2, compare);
	int	*q = buf;
	for(;;) {
#ifdef	COUNT
		move_count++;
#endif
		if(compare(*p1, *p2) <= 0) {
			*q++ = *p1++;
			if(--n1 == 0) break;
		} else {
			*q++ = *p2++;
			if(--n2 == 0) break;
		}
	}
#ifdef	COUNT
	move_count += n1 + len - n2;
#endif
	while(n1 > 0) {
		*q++ = *p1++;
		n1--;
	}
	for(int i = 0; i < len - n2; i++)
		a[i] = buf[i];
}
void	qsort_glibc(int *a, int len, COMPARE compare) {
	if(len < 2) return;
	if(len <= TEMP_SIZE) {
		int	tmp[TEMP_SIZE];
		sort(tmp, a, len, compare);
	} else {
		int	*buf = (int*)malloc(sizeof(int) * len);
		if(buf == (int*)NULL) return;
		sort(buf, a, len, compare);
		free(buf);
	}
}
```

</details>

<details><summary><b>メイクファイル</b></summary>

```:Makefile
CC=cl
CFLAGS=/c /O2 /W4 /arch:AVX2 /GL /GS-
LFLAGS=/GL
LIBS=

OUTDIRS=count count\ucrt count\glibc \
        clock clock\ucrt clock\glibc

target: $(OUTDIRS) \
        count\ucrt\bench.exe count\glibc\bench.exe \
        clock\ucrt\bench.exe clock\glibc\bench.exe

$(OUTDIRS):
	@if not exist $@ mkdir $@

count\ucrt\bench.exe:  count\ucrt\bench.obj  count\ucrt\qsort_ucrt.obj
count\glibc\bench.exe: count\glibc\bench.obj count\glibc\qsort_glibc.obj
clock\ucrt\bench.exe:  clock\ucrt\bench.obj  clock\ucrt\qsort_ucrt.obj
clock\glibc\bench.exe: clock\glibc\bench.obj clock\glibc\qsort_glibc.obj

count\ucrt\bench.exe \
count\glibc\bench.exe \
clock\ucrt\bench.exe \
clock\glibc\bench.exe:
	$(CC) $** /Fe$* $(LFLAGS) $(LIBS)

.c{count\ucrt}.obj::
	$(CC) $(CFLAGS) $< /Focount\ucrt\ /DSORTFUNC=qsort_ucrt /DCOUNT

.c{count\glibc}.obj::
	$(CC) $(CFLAGS) $< /Focount\glibc\ /DSORTFUNC=qsort_glibc /DCOUNT

.c{clock\ucrt}.obj::
	$(CC) $(CFLAGS) $< /Foclock\ucrt\ /DSORTFUNC=qsort_ucrt

.c{clock\glibc}.obj::
	$(CC) $(CFLAGS) $< /Foclock\glibc\ /DSORTFUNC=qsort_glibc

clean:
	if exist count rd /s /q count
	if exist clock rd /s /q clock
```

</details>

上記４つのファイルを同じフォルダに保存して `nmake` を実行すると，下記のプログラムが作られる。

<table>
<caption>表３　作成プログラム一覧</caption>
<thead>
<tr><th></th>
    <th style="text-align:center">uCRT 版</th>
    <th style="text-align:center">glibc 版</th></tr>
</thead>
<tbody>
<tr></tr>
<tr><th>比較・転送回数計測用</th>
    <td><code>clock\ucrt\bench.exe</code></td>
    <td><code>clock\glibc\bench.exe</code></td></tr>
<tr></tr>
<tr><th>時間計測用</th>
    <td><code>count\ucrt\bench.exe</code></td>
    <td><code>count\glibc\bench.exe</code></td></tr>
</tbody>
</table>

# 次回

次回に続く。

- [最強コンパイラを使えばマージソートは遅くない？ - Qiita](https://qiita.com/TETSURO1999/items/0cc093f427531407d2d4)
- [マージソートを非再帰型に書き直す～誰でも思いつくアイディアだが～ - Qiita](https://qiita.com/TETSURO1999/items/23bdd493975cbd8a5085)
- [マージソートを非再帰型に書き直す２～さらなる改善を試みる～ - Qiita](https://qiita.com/TETSURO1999/items/783aa5732b9f19a615ca)
- [マージソートを非再帰型に書き直す３～もう限界かもしれない～ - Qiita](https://qiita.com/TETSURO1999/items/5aa007604198b2ff7b06)
- [マージソートを非再帰型に書き直す４～限界突破のアイディアを思いつく～ - Qiita](https://qiita.com/TETSURO1999/items/17d751555fcf9b275518)

<!-- 参考文献 -->
[^1]: [オンラインでプログラミングして実行できるサイト - Qiita](https://qiita.com/tttamaki/items/2b009aa957cfa4895d50)

[^2]: [qsortの実装はなぜマージソート？ 高速な安定ソートとは？ (qs14) - Qiita](https://qiita.com/t-kawa/items/d62ee8c20af652ac87fe)

[^3]: [GNU C Library 2.41 stdlib/qsort.c - glibc.git](https://sourceware.org/git/?p=glibc.git;a=blob_plain;f=stdlib/qsort.c;hb=e5893e6349541d871e8a25120bca014551d13ff5)

[^4]: [Ask Reddit: Why does the GNU libc use mergesort to implement qsort()?](https://www.reddit.com/r/programming/comments/1jxq4/ask_reddit_why_does_the_gnu_libc_use_mergesort_to/?rdt=44001)