典型的な条件のもとでは、ソートに時間計算量O(n log n)かかるのはよく知られた話です。そして実装には、例えばC/C++でいう qsort()
/ std::sort()
のような標準ライブラリを使うのがよくある方法だと思います。
ただし、ソート対象や条件に何らかの追加の仮定を置いてよいなら、(理論的あるいは実用的な観点で)より高速な実装が得られることがあります。
例えば値の分布によっては基数ソートやバケットソートが有用なことがあるかもしれません。
この記事で話題にするのは、「ソート対象の要素数が少なかったら?」です。
ソーティングネットワーク
要素数nが事前に分かっているとき、直感的には「何番目の要素と何番目の要素を比較していけばソートが完成するのか、最適な手順はnによって完全に決まるのではないか?」というアイディアが思い浮かびます。
この議論はソーティングネットワークという分野で研究されているようです。
この記事では主に実装を扱うため、ソーティングネットワークの詳細には触れませんが、こんなようなものです。
- 大小を比較して交換するコンポーネント(コンパレータ)を入力に対して並べることで、入力をソートした結果を出力するもの(ソーティングネットワーク)を作れる
- nが小さい時の最適なソーティングネットワーク(コンパレータの個数が最小のもの)は既知となっている
- 例としてn=8のとき、19個7段のコンパレータがあればソーティングネットワークが構成できる
つまり、この最適なソーティングネットワークを使ってソートすれば、とても効率が良いのでは? ということが考えられます。
※ ここの議論は少し粗雑で、比較結果に応じて次の比較対象を変えてよいのであれば、比較回数をより少なく(n=8のとき16回に)できることが知られているようです。一方で、ソーティングネットワークは比較結果に関わらず、次の比較対象は固定されています。このことを無記憶性という言い方をするようです。比較回数は少ないほうが本来嬉しいはずですが、ソーティングネットワークは高効率な実装が可能でした(この性質がなければ難しかったのではと思います)。
実装して評価してみる
お題として、8要素の32ビット整数のソート問題を考えます。言語はC/C++1とします。
まずコンパレータを実装します。
inline void network_sort_swap(int32_t *v, int a, int b)
{
if (v[a] > v[b]) {
// swap
auto tmp = v[a];
v[a] = v[b];
v[b] = tmp;
}
}
具体的なソーティングネットワークの構成はソルバがあるので、ソルバの出力を使うことにします。フォーマットの違いを揃えるためにマクロを噛ませています。
#define SWAP(a, b) network_sort_swap(v, a, b)
void network_sort_8n(int32_t *v)
{
SWAP(0, 1); SWAP(2, 3); SWAP(4, 5); SWAP(6, 7);
SWAP(0, 2); SWAP(1, 3); SWAP(4, 6); SWAP(5, 7);
SWAP(1, 2); SWAP(5, 6); SWAP(0, 4); SWAP(3, 7);
SWAP(1, 5); SWAP(2, 6);
SWAP(1, 4); SWAP(3, 6);
SWAP(2, 4); SWAP(3, 5);
SWAP(3, 4);
}
さて、大量にnetwork_sort_8n()
を呼び出すベンチマークを書きます。
// https://gcc.gnu.org/onlinedocs/gcc-8.2.0/gcc/Machine-Constraints.html#Machine-Constraints
inline uint64_t rdtsc()
{
uint32_t tickl, tickh;
__asm__ __volatile__("rdtsc":"=a"(tickl),"=d"(tickh));
return (static_cast<uint64_t>(tickh) << 32) | tickl;
}
// xorshift: https://ja.wikipedia.org/wiki/Xorshift
uint32_t rnd(void)
{
static uint32_t y = 2463534242;
y = y ^ (y << 13); y = y ^ (y >> 17);
return y = y ^ (y << 15);
}
void benchmark()
{
constexpr const int N = 80 * 1000 * 1000;
// 領域vを確保し乱数で埋める
auto v = static_cast<int32_t*>(calloc(sizeof(int32_t), N));
for (int i = 0; i < N; i++) {
v[i] = rnd();
}
auto t0 = rdtsc();
for (int n = 0; n < N; n += 8) {
network_sort_8n(v + n);
//std::sort(v + n, v + n + 8);
}
auto t1 = rdtsc();
printf("%lld [Mcycles]\n", (t1 - t0)/1000/1000);
}
- ソート対象の要素数は80M(8000万個)の乱数です
-
rdtsc()
から取得できるのはクロックカウンタです。性能測定区間をrdtsc()
で挟み差を取ることで、処理にかかったクロックサイクル数が出せます。 -
network_sort_8n()
の比較対象として今回使うのはstd::sort()
です
benchmark()
を実行することでソーティングネットワークの評価をしていきます。
評価結果
手元の少し古いMacBook Pro(Late 2013)で実行しました。処理系はClang 7.0です。
項目 | 結果(Mcycles2) |
---|---|
network_sort_8n() + Clang(-O1) | 2928 |
network_sort_8n() + Clang(-O3) | 1582 |
std::sort() + Clang(-O1) | 1986 |
std::sort() + Clang(-O3) | 1961 |
1961/1582 ≒ 1.24ということは、20%ほど高速化しました。
実装の改善方法を考える
n=8に特化した割に2割改善じゃあんまりおいしくない!! ということで改善方法を考えます。
ソーティングネットワークは処理時間のほぼ全てがnetwork_sort_swap()
にかかっているはずです。network_sort_swap()
を再掲します。
inline void network_sort_swap(int32_t *v, int a, int b)
{
if (v[a] > v[b]) {
// swap
auto tmp = v[a];
v[a] = v[b];
v[b] = tmp;
}
}
遅そうなのがif文です。評価ではランダムな入力値の比較をしているため、if文で実行がどちらに分岐するかは予想困難なはずで、したがって分岐予測ミスのコストが非常に大きいと推測できます。
ということで、なんとかしてif文を消すことを考えます。
この箇所でやりたいのは、v[a]
とv[b]
がそれぞれ2値の小さい方・大きい方になっていてほしいということです。言い換えると2値のmin, maxが代入されればよいということなので、次のような書き換えができます。
inline void network_sort_swap(int32_t *v, int a, int b)
{
auto va = v[a];
auto vb = v[b];
v[a] = std::min(va, vb);
v[b] = std::max(va, vb);
}
これも結局は内部でif使ってるのではという思いになりつつも、再評価します。
項目 | 結果(Mcycles2) |
---|---|
network_sort_8n() + Clang(-O1) | 1273 |
network_sort_8n() + Clang(-O3) | 390 |
かなり良くなりました。1961/390 ≒ 5.0 ということは、std::sort()
と比べて5倍ぐらい速くなりました。network_sort_swap
は 390/(80/8*19) ≒ 2 なので、平均2cycle/call程度しかかかっていないと考えられます。とても効率が良いと言えそうです。
アセンブリ出力を確認すると、Clangはstd::min
, std::max
をcmov系の命令で実現することで、分岐命令を使わないコード生成をしてくれたようです3。
理想的な出力になっているか?
Clang7.0では、network_sort_swap()
は次のアセンブリコードになるようです(最適化コンパイルしたもの)。
.globl __Z17network_sort_swapPiii
.align 4, 0x90
__Z17network_sort_swapPiii: ## @_Z17network_sort_swapPiii
.cfi_startproc
## BB#0:
push rbp
Ltmp0:
.cfi_def_cfa_offset 16
Ltmp1:
.cfi_offset rbp, -16
mov rbp, rsp
Ltmp2:
.cfi_def_cfa_register rbp
movsxd r8, esi
mov eax, dword ptr [rdi + 4*r8]
movsxd rdx, edx
mov esi, dword ptr [rdi + 4*rdx]
cmp esi, eax
mov ecx, eax
cmovle ecx, esi
mov dword ptr [rdi + 4*r8], ecx
cmp eax, esi
cmovl eax, esi
mov dword ptr [rdi + 4*rdx], eax
pop rbp
ret
.cfi_endproc
cmp
が何故か2回あることに気づきます。cmovはmovの条件を自由に指定できるので、1つでいいはずです。これはおそらくstd::min
/ std::max
の実装で大小比較していることから来ています。
次に考えたのは、「min, maxの実装で条件分岐が共通部分式にできることにコンパイラは気づいていないのでは?」ということでした。そこで、条件式を明示的に揃えてみます。
inline void network_sort_swap(int32_t *v, int a, int b)
{
auto va = v[a];
auto vb = v[b];
v[a] = va < vb ? va : vb;
v[b] = va < vb ? vb : va;
}
対応するアセンブリ出力です。
.globl __Z17network_sort_swapPiii
.align 4, 0x90
__Z17network_sort_swapPiii: ## @_Z17network_sort_swapPiii
.cfi_startproc
## BB#0:
push rbp
Ltmp0:
.cfi_def_cfa_offset 16
Ltmp1:
.cfi_offset rbp, -16
mov rbp, rsp
Ltmp2:
.cfi_def_cfa_register rbp
movsxd rax, esi
mov ecx, dword ptr [rdi + 4*rax]
movsxd r8, edx
mov edx, dword ptr [rdi + 4*r8]
cmp ecx, edx
mov esi, edx
cmovle esi, ecx
mov dword ptr [rdi + 4*rax], esi
cmovl ecx, edx
mov dword ptr [rdi + 4*r8], ecx
pop rbp
ret
.cfi_endproc
1つのcmp
の結果を使ってcmovを2回することでmin, maxを得る、期待した最適化になっています。
項目 | 結果(Mcycles2) |
---|---|
network_sort_8n() + Clang(-O3) | 340 |
性能的にも、さらに15%ほどの改善になっていることが分かります。
紆余曲折を経ながら3種類のnetwork_sort_swap()
を実装しました。非常に微妙な実装の違いながら、実行効率には明らかな差が出るという、興味深い結果となりました。
この結果はコンパイラにも強く依存します。Clangはcmovを使ったコード生成をしてくれましたが、GCCは8.2ではしてくれず、trunkではやってくれるようです。
まとめ
- ソーティングネットワークの実装を通して、特化したソートアルゴリズムが典型的なソートアルゴリズムより高速になる例を見てみました
- ちょっとした実装のチューニングで大きな効果を得られるケースにも遭遇しました
Appendix
余談
筆者がソーティングネットワークを知ったのはLet Over Lambdaを読んだときです。難解な記述が多かった(理解できない箇所がたくさんあった)のですが、面白い本でした。
ソースコード
参考リンク
- http://pages.ripco.net/~jgamble/nw.html
- https://stackoverflow.com/questions/2786899/fastest-sort-of-fixed-length-6-int-array
- https://stackoverflow.com/q/8673860