std::sortより特定条件で高速な8入力ソーティングネットワークを作ってみる

典型的な条件のもとでは、ソートに時間計算量O(n log n)かかるのはよく知られた話です。そして実装には、例えばC/C++でいう `qsort()` / `std::sort()` のような標準ライブラリを使うのがよくある方法だと思います。

ただし、ソート対象や条件に何らかの追加の仮定を置いてよいなら、（理論的あるいは実用的な観点で）より高速な実装が得られることがあります。
例えば値の分布によっては[基数ソート](https://ja.wikipedia.org/wiki/%E5%9F%BA%E6%95%B0%E3%82%BD%E3%83%BC%E3%83%88)や[バケットソート](https://ja.wikipedia.org/wiki/%E3%83%90%E3%82%B1%E3%83%83%E3%83%88%E3%82%BD%E3%83%BC%E3%83%88)が有用なことがあるかもしれません。

この記事で話題にするのは、「ソート対象の要素数が少なかったら？」です。

## ソーティングネットワーク

要素数nが事前に分かっているとき、直感的には「何番目の要素と何番目の要素を比較していけばソートが完成するのか、最適な手順はnによって完全に決まるのではないか？」というアイディアが思い浮かびます。
この議論はソーティングネットワークという分野で研究されているようです。

- [ソーティングネットワーク - Wikipedia](https://ja.wikipedia.org/wiki/%E3%82%BD%E3%83%BC%E3%83%86%E3%82%A3%E3%83%B3%E3%82%B0%E3%83%8D%E3%83%83%E3%83%88%E3%83%AF%E3%83%BC%E3%82%AF)

この記事では主に実装を扱うため、ソーティングネットワークの詳細には触れませんが、こんなようなものです。

- 大小を比較して交換するコンポーネント（コンパレータ）を入力に対して並べることで、入力をソートした結果を出力するもの（ソーティングネットワーク）を作れる
- nが小さい時の最適なソーティングネットワーク（コンパレータの個数が最小のもの）は既知となっている
  - 例としてn=8のとき、19個7段のコンパレータがあればソーティングネットワークが構成できる

つまり、この最適なソーティングネットワークを使ってソートすれば、とても効率が良いのでは？ ということが考えられます。

_※ ここの議論は少し粗雑で、比較結果に応じて次の比較対象を変えてよいのであれば、比較回数をより少なく（n=8のとき16回に）できることが知られているようです。一方で、ソーティングネットワークは比較結果に関わらず、次の比較対象は固定されています。このことを無記憶性という言い方をするようです。比較回数は少ないほうが本来嬉しいはずですが、ソーティングネットワークは高効率な実装が可能でした（この性質がなければ難しかったのではと思います）。_

## 実装して評価してみる

お題として、8要素の32ビット整数のソート問題を考えます。言語はC/C++[^c-cpp]とします。

[^c-cpp]: CともC++とも言えないような、渾然としたコードになっています。C++11以降の処理系ならコンパイル通るはず。

まずコンパレータを実装します。

```cpp
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;
  }
}
```

具体的なソーティングネットワークの構成は[ソルバがある](http://pages.ripco.net/~jgamble/nw.html)ので、ソルバの出力を使うことにします。フォーマットの違いを揃えるためにマクロを噛ませています。

```cpp
#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()`を呼び出すベンチマークを書きます。

```cpp
// 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です。

| 項目                               | 結果(Mcycles[^Mcycles]) |
|-----------------------------------|------------------------|
| network_sort_8n() + Clang(-O1)    | 2928                   |
| network_sort_8n() + Clang(-O3)    | 1582                   |
| std::sort()       + Clang(-O1)    | 1986                   |
| std::sort()       + Clang(-O3)    | 1961                   |

[^Mcycles]: mega cyclesつまり100万クロックサイクルです。

1961/1582 ≒ 1.24ということは、20%ほど高速化しました。

## 実装の改善方法を考える

n=8に特化した割に2割改善じゃあんまりおいしくない！！ ということで改善方法を考えます。

ソーティングネットワークは処理時間のほぼ全てが`network_sort_swap()`にかかっているはずです。`network_sort_swap()`を再掲します。

```cpp
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が代入されればよいということなので、次のような書き換えができます。

```cpp
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使ってるのではという思いになりつつも、再評価します。

| 項目                               | 結果(Mcycles[^Mcycles]) |
|-----------------------------------|------------------------|
| 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系の命令](https://www.jaist.ac.jp/iscenter-new/mpc/altix/altixdata/opt/intel/vtune/doc/users_guide/mergedProjects/analyzer_ec/mergedProjects/reference_olh/mergedProjects/instructions/instruct32_hh/vc35.htm)で実現することで、分岐命令を使わないコード生成をしてくれたようです[^use-cmov]。

[^use-cmov]: [Intelのアーキテクチャ最適化マニュアル](https://software.intel.com/sites/default/files/managed/9e/bc/64-ia-32-architectures-optimization-manual.pdf)によると、cmovを使った分岐命令の削除は、分岐が予測不能である場合に検討すべきとされているようです。まさに今回のような状況には当てはまるように思います。

## 理想的な出力になっているか？

Clang7.0では、`network_sort_swap()`は次のアセンブリコードになるようです（最適化コンパイルしたもの）。

```nasm
        .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の実装で条件分岐が共通部分式にできることにコンパイラは気づいていないのでは？」ということでした。そこで、条件式を明示的に揃えてみます。

```cpp
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;
}
```

対応するアセンブリ出力です。

```nasm
  .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を得る、期待した最適化になっています。

| 項目                               | 結果(Mcycles[^Mcycles]) |
|-----------------------------------|------------------------|
| network_sort_8n() + Clang(-O3)    | 340                    |

性能的にも、さらに15%ほどの改善になっていることが分かります。

紆余曲折を経ながら3種類の`network_sort_swap()`を実装しました。非常に微妙な実装の違いながら、実行効率には明らかな差が出るという、興味深い結果となりました。
この結果はコンパイラにも強く依存します。Clangはcmovを使ったコード生成をしてくれましたが、GCCは8.2ではしてくれず、trunkではやってくれるようです。

## まとめ

- ソーティングネットワークの実装を通して、特化したソートアルゴリズムが典型的なソートアルゴリズムより高速になる例を見てみました
- ちょっとした実装のチューニングで大きな効果を得られるケースにも遭遇しました

## Appendix

### 余談

筆者がソーティングネットワークを知ったのは[Let Over Lambda](https://letoverlambda.com/)を読んだときです。難解な記述が多かった（理解できない箇所がたくさんあった）のですが、面白い本でした。

### ソースコード
- <https://gist.github.com/saka1/beedcd65b7f3bfa209b3776d49ade008>

### 参考リンク

- <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>