Help us understand the problem. What is going on with this article?

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

典型的な条件のもとでは、ソートに時間計算量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::maxcmov系の命令で実現することで、分岐命令を使わないコード生成をしてくれたようです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を読んだときです。難解な記述が多かった(理解できない箇所がたくさんあった)のですが、面白い本でした。

ソースコード

参考リンク


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

  2. mega cyclesつまり100万クロックサイクルです。 

  3. Intelのアーキテクチャ最適化マニュアルによると、cmovを使った分岐命令の削除は、分岐が予測不能である場合に検討すべきとされているようです。まさに今回のような状況には当てはまるように思います。 

Why do not you register as a user and use Qiita more conveniently?
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away
Comments
Sign up for free and join this conversation.
If you already have a Qiita account
Why do not you register as a user and use Qiita more conveniently?
You need to log in to use this function. Qiita can be used more conveniently after logging in.
You seem to be reading articles frequently this month. Qiita can be used more conveniently after logging in.
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away