アセンブリ解読 その4

その2で宣言した通り、今回の対象はループだ。その2に続きが発生したので番号的には飛んでしまっているが、続きという点は変わらない。 しかしループとなると絶対にSIMDとループアンロールが絡んで面倒なことになる、と思っているので少し腰が重いが、最適化の勉強というならむしろこれこそやらねばならない。

今回とりあえずこんなコードを書いた。

#include <stdint.h>

int32_t accum(const int32_t* begin, const int32_t* end)
{
    int32_t sum = 0;
    for(; begin != end; ++begin)
        sum += *begin;
    return sum;
}

beginからendが到達可能ならそこまでの範囲を総和する。到達不能なら知らない。

最適化なし

では最適化なしでできたバイナリを逆アセンブルしよう。

$ gcc -std=c99 -c -O0 accumulate.c && objdump -d accumulate.o > accumulate.asm

すると、

accumulate.o:     file format elf64-x86-64


Disassembly of section .text:

0000000000000000 <accum>:
   0: push   %rbp            ; いつもの
   1: mov    %rsp,%rbp       ; いつもの
   4: mov    %rdi,-0x18(%rbp); 第一引数をスタックへ
   8: mov    %rsi,-0x20(%rbp); 第二引数をスタックへ
   c: movl   $0x0,-0x4(%rbp) ; sumを即値0で初期化
  13: jmp    23 <accum+0x23> ; 23へ無条件ジャンプ
  15: mov    -0x18(%rbp),%rax; beginの持つアドレス値をraxへロード
  19: mov    (%rax),%eax     ; raxの指すメモリの内容をeaxへロード
  1b: add    %eax,-0x4(%rbp) ; eaxをsumへ足す
  1e: addq   $0x4,-0x18(%rbp); beginに即値4(sizeof(int32_t))を足す
  23: mov    -0x18(%rbp),%rax; 現在のポインタbeginをraxレジスタへロード
  27: cmp    -0x20(%rbp),%rax; endとbeginを比較。同一ならZF==1
  2b: jne    15 <accum+0x15> ; ZF==0なら15まで戻る
  2d: mov    -0x4(%rbp),%eax ; sumを返却値レジスタへ
  30: pop    %rbp            ; いつもの
  31: retq   

既にかなり短い。そして、ここまでやった内容で十分わかる。で、ここまでの経験からこれを最適化するとスタックメモリを使わないように頑張ると予想され、以下のようなコードが想定される(動作確認とかしてないのでミスや誤解があるかも)。

accumulate.o:     file format elf64-x86-64


Disassembly of section .text:

0000000000000000 <accum>:
  movl   $0x0,%eax
  jmp    judge
inside_loop:
  add    (%rdi),%eax
  addq   $0x4,%rdi
judge:
  cmp    %rsi,%rdi
  jne    inside_loop <accum+0x15>
  retq

だが、現実は甘くはなかった。

-O3

では最適化しよう。

$ gcc -std=c99 -c -O3 accumulate.c && objdump -d accumulate.o > accumulate_O3.asm

結果、目を疑うようなコードが出てきた。

accumulate.o:     file format elf64-x86-64


Disassembly of section .text:

0000000000000000 <accum>:
   0: cmp    %rsi,%rdi
   3: je     128 <accum+0x128>
   9: lea    0x4(%rdi),%rcx
   d: mov    %rsi,%r9
  10: mov    %rdi,%rdx
  13: and    $0xf,%edx
  16: sub    %rcx,%r9
  19: shr    $0x2,%rdx
  1d: shr    $0x2,%r9
  21: neg    %rdx
  24: lea    0x1(%r9),%r8
  28: and    $0x3,%edx
  2b: cmp    %r8,%rdx
  2e: cmova  %r8,%rdx
  32: cmp    $0x8,%r8
  36: ja     138 <accum+0x138>
  3c: mov    %r8,%rdx
  3f: cmp    $0x1,%rdx
  43: mov    (%rdi),%eax
  45: je     9c <accum+0x9c>
  47: add    0x4(%rdi),%eax
  4a: cmp    $0x2,%rdx
  4e: lea    0x8(%rdi),%rcx
  52: je     9c <accum+0x9c>
  54: add    0x8(%rdi),%eax
  57: cmp    $0x3,%rdx
  5b: lea    0xc(%rdi),%rcx
  5f: je     9c <accum+0x9c>
  61: add    0xc(%rdi),%eax
  64: cmp    $0x4,%rdx
  68: lea    0x10(%rdi),%rcx
  6c: je     9c <accum+0x9c>
  6e: add    0x10(%rdi),%eax
  71: cmp    $0x5,%rdx
  75: lea    0x14(%rdi),%rcx
  79: je     9c <accum+0x9c>
  7b: add    0x14(%rdi),%eax
  7e: cmp    $0x6,%rdx
  82: lea    0x18(%rdi),%rcx
  86: je     9c <accum+0x9c>
  88: add    0x18(%rdi),%eax
  8b: cmp    $0x8,%rdx
  8f: lea    0x1c(%rdi),%rcx
  93: jne    9c <accum+0x9c>
  95: add    0x1c(%rdi),%eax
  98: lea    0x20(%rdi),%rcx
  9c: cmp    %rdx,%r8
  9f: je     130 <accum+0x130>
  a5: sub    %rdx,%r8
  a8: sub    %rdx,%r9
  ab: lea    -0x4(%r8),%r10
  af: shr    $0x2,%r10
  b3: add    $0x1,%r10
  b7: cmp    $0x2,%r9
  bb: lea    0x0(,%r10,4),%r11
  c2:
  c3: jbe    109 <accum+0x109>
  c5: pxor   %xmm0,%xmm0
  c9: lea    (%rdi,%rdx,4),%rdi
  cd: xor    %edx,%edx
  cf: add    $0x1,%rdx
  d3: paddd  (%rdi),%xmm0
  d7: add    $0x10,%rdi
  db: cmp    %rdx,%r10
  de: ja     cf <accum+0xcf>
  e0: movdqa %xmm0,%xmm1
  e4: lea    (%rcx,%r11,4),%rcx
  e8: psrldq $0x8,%xmm1
  ed: paddd  %xmm1,%xmm0
  f1: movdqa %xmm0,%xmm1
  f5: psrldq $0x4,%xmm1
  fa: paddd  %xmm1,%xmm0
  fe: movd   %xmm0,%edx
 102: add    %edx,%eax
 104: cmp    %r11,%r8
 107: je     12a <accum+0x12a>
 109: lea    0x4(%rcx),%rdx
 10d: add    (%rcx),%eax
 10f: cmp    %rdx,%rsi
 112: je     12a <accum+0x12a>
 114: lea    0x8(%rcx),%rdx
 118: add    0x4(%rcx),%eax
 11b: cmp    %rdx,%rsi
 11e: je     150 <accum+0x150>
 120: add    0x8(%rcx),%eax
 123: retq   
 124: nopl   0x0(%rax)
 128: xor    %eax,%eax
 12a: repz retq 
 12c: nopl   0x0(%rax)
 130: repz retq 
 132: nopw   0x0(%rax,%rax,1)
 138: test   %rdx,%rdx
 13b: jne    3f <accum+0x3f>
 141: mov    %rdi,%rcx
 144: xor    %eax,%eax
 146: jmpq   a5 <accum+0xa5>
 14b: nopl   0x0(%rax,%rax,1)
 150: repz retq 

な ん だ こ れ は

基本的に命令の数が少ないほど速度は向上すると思っていたので、これは以外だ。確かにループを一回一回回るよりも効率のよい(と思われる)方法はたくさんある。冒頭でも言及したループアンローリングSIMD命令などだ。だがここまでくると本当にこの方が速いのか少し疑わしくなる。

ついでに、くやしいので私が予想したコードを出しそうなオプションを使ってみる。-Osだ。これはバイナリのサイズを小さくするように働く。

$ gcc -std=c99 -c -Os accumulate.c && objdump -d accumulate.o > accumulate_Os.asm
accumulate.o:     file format elf64-x86-64


Disassembly of section .text:

0000000000000000 <accum>:
   0: xor    %eax,%eax
   2: cmp    %rsi,%rdi
   5: je     f <accum+0xf>
   7: add    (%rdi),%eax
   9: add    $0x4,%rdi
   d: jmp    2 <accum+0x2>
   f: retq   

私の書いたコードと異なる点がいくつかある。前回も見たが、ゼロクリアは即値の代入でなくxorになっている。初っ端からジャンプするのではなく、先にcmpjeを発行している。私が書いたコードと比べて、渡された領域のサイズが0だった場合にジャンプが一回少なくなる。だが基本的には同じ発想だ。

また、O2も試した。

accumulate.o:     file format elf64-x86-64


Disassembly of section .text:

0000000000000000 <accum>:
 0: xor    %eax,%eax
 2: cmp    %rsi,%rdi
 5: je     1d <accum+0x1d>
 7: nopw   0x0(%rax,%rax,1)
 e:
10: add    (%rdi),%eax
12: add    $0x4,%rdi
16: cmp    %rdi,%rsi
19: jne    10 <accum+0x10>
1b: repz retq
1d: repz retq

こちらも同様だ。最低限の命令数で素直にループを回っている。

マイクロベンチマーク

O3のあのとても長いバイナリで本当に高速になるとは俄には信じられないのでちょっとしたマイクロベンチマークをする。

#include "accumulate.h"
#include <stdio.h>
#include <stddef.h>

int main()
{
    const size_t N = 100000;
    int32_t a[N];
    for(int32_t i=0; i<N; ++i)
        a[i] = i;

    for(size_t i = 0; i < N; ++i)
    {
        int32_t b = accum(a, a+N);
    }
    return 0;
}

前回見たとおりLTOなしで別々にコンパイルした後リンクすればインライン化の恐れはない(と思いたかったが、一応確認した)。

ここで、-march=nativeも戯れにコンパイルし、逆アセンブルした。余裕があれば後で持ち出すが、SIMDまわりの違いだった。結果は、

$ time ./opts # -Osでコンパイルした
3.23user 0.00system 0:03.23elapsed 99%CPU
$ time ./opt3 # -O3でコンパイルした
1.66user 0.00system 0:01.66elapsed 99%CPU
$ time ./optnative # -O3 -march=nativeでコンパイルした
1.18user 0.00system 0:01.18elapsed 99%CPU

となった。実際ハヤイ。ちなみに3回くらい試したが、ほぼこの関係は変わらなかった。

解読・前半戦

こうなるとコンパイラの出した答えに納得するしかないので最適化済みのアセンブリに戻ろう。とりあえずもう一度ざっと見て、何が起きているか把握したい。そして、できれば分割して解読したい。長くなると自分の考えを挟みづらいからだ。 しかし長いので、予想を立てる。多分、ここまで何度か言及しているループアンロールとSIMDを使いこなしているに違いない。というかそれ以外では並列化しかこの手のループが早くなる方法を私は知らない。そして多分並列化はされていない。

この予想を念頭に置いてもう一度前からざっと見ていこう。 まず、最初にポインタに対する演算をいくつかしている。その後、似たような命令のパターンが7回ほど続き、それに続く少し複雑そうな命令の中でpadddが使われている。確定だ。SIMD命令が使われている。 後ろでSIMD命令を使うことを考え、先に何回SIMD演算をすることになるかの計算と、あまりの分の処理をしておいているのだろう。

大まかな見当がついたので、分割して順に見ていこう。最初の部分は最後尾へのジャンプをいくつか含むので最後尾も一緒に見る。

0000000000000000 <accum>:
   0: cmp    %rsi,%rdi         ; 引数同士を比較
   3: je     128 <accum+0x128> ; 同じなら範囲が空なので0を返す
                               ; ここ以降では要素数が1以上と確定する
   9: lea    0x4(%rdi),%rcx    ; rdi+4をrcxへ代入
   d: mov    %rsi,%r9          ; endをr9レジスタへコピー
  10: mov    %rdi,%rdx         ; beginをrdxレジスタへコピー
  13: and    $0xf,%edx         ; edxの下4bitだけを残す
  16: sub    %rcx,%r9          ; r9(end)からrcx(begin+1)を引く
  19: shr    $0x2,%rdx         ; rdxを2回右シフト
  1d: shr    $0x2,%r9          ; r9を2回右シフト
  21: neg    %rdx              ; rdxの符号を反転
  24: lea    0x1(%r9),%r8      ; r9+0x1をr8へ
  28: and    $0x3,%edx         ; edxの下位2bitだけを残す
  2b: cmp    %r8,%rdx          ; r8とrdxを比較
  2e: cmova  %r8,%rdx          ; if(r8 < rdx) rdx = r8
  32: cmp    $0x8,%r8          ; 8とr8を比較
  36: ja     138 <accum+0x138> ; if(8 < r8) goto 138
;...
 128: xor    %eax,%eax         ; eaxをゼロクリア
 12a: repz retq                ; return
 12c: nopl   0x0(%rax)         ;
 130: repz retq                ; return
 132: nopw   0x0(%rax,%rax,1)  ;
 138: test   %rdx,%rdx         ; rdxが0ならZF=1
 13b: jne    3f <accum+0x3f>   ; ZF == 0なら3fへジャンプ
 141: mov    %rdi,%rcx         ; rdiをrcxへ代入
 144: xor    %eax,%eax         ; eaxをゼロクリア
 146: jmpq   a5 <accum+0xa5>   ; a5へジャンプ
 14b: nopl   0x0(%rax,%rax,1)  ;
 150: repz retq

最初の2行はわかりやすい。範囲が空だった場合に即リターンするためのものだ。そして、続く部分では要素が少なくとも1つ以上あるという前提で話を進めている。

先頭では2つのことが同時に行われているようで、意識があっちこっちに振り回されて読みにくい。とりあえず片方を無視してわかりやすい方だけ追っていこう。そちら側だけ抽出するとこうなる。

   9: lea    0x4(%rdi),%rcx    ; begin+1をrcxへ代入
   d: mov    %rsi,%r9          ; endをr9レジスタへコピー
  16: sub    %rcx,%r9          ; r9(end)からrcx(begin+1)を引く
  1d: shr    $0x2,%r9          ; r9((要素数-1)*4)を2回右シフト
  24: lea    0x1(%r9),%r8      ; r9+0x1をr8へ

まず9:begin+1を用意している(ここでsizeof(int32_t) == 4なので、0x4(%rdi)begin+1を意味していることになる)。 続いて、d:endr9へコピーし、16:でそこからbegin+1を引いている。これで、(要素数-1)*4の値が得られた(この最後の4倍はsizeof(int32_t)だ)。 これを1d:で2回右シフトすることで、要素数-1を得ている。右シフト2回は4での除算になるからだ。そして24:で要素数そのものがr8へ代入されている。 得られた値がわかりやすいのでこちらが何をしているかはだいたいわかったと思っていいだろう。反対側は何をしているんだろう。

  10: mov    %rdi,%rdx         ; beginをrdxレジスタへコピー
  13: and    $0xf,%edx         ; edxの下4bitだけを残す
  19: shr    $0x2,%rdx         ; rdxを2回右シフト
  21: neg    %rdx              ; rdxの符号を反転
  28: and    $0x3,%edx         ; edxの下位2bitだけを残す

じっくり何が起きるか考えてみないことには何が得られるかすらわからない。使われ方を追いかけてみたが、色々なところにジャンプするばかりでそもそもこの値が何かわからないと理解の足しにはならなかった。 そもそも、これはbeginのアドレス値に対してなにかしらを行っているので、要素の数とかそういうことはわからないはずだ(endの値を使っていないので)。だとするとなんだろう。 アドレスそのものの値に関係し、ループするときに、SIMDの前にしておく必要があり、その値によっては先に要素を処理しておいたりしなければならない……。ここで思いついた。メモリアライメントだ。 アライメントだとすれば、下位4bitだけ残した後に2回右シフトしてその後下位2bitだけ残した理由がわかる(符号反転は直感ではわからないが)。16バイト境界にアラインしているか確認しているのだ。

本当にそうだろうか。それは、今後のメモリの使われ方を確認すればわかる。SIMD命令の初っ端を確認すればいい。メモリからロードしているなら、整列していることが要求される。

  d3: paddd  (%rdi),%xmm0

rdiのアドレスが直接使われている。paddd命令は32bit整数4つを一度に加算する命令だ。ソースオペランドがメモリアドレスなら、それは16バイト境界に整列していなければならない。 ここでは要素数の余りの他にも、メモリアライメントの余りも考慮したバイナリが生成されているのだ。後でのSIMD命令を可能な限り効率よくするために! それによく見てみれば、SIMD演算が終了した後にひっそりと、通常のadd命令が含まれている。これは、アライメントされた領域を脱した後に要素が残っているならそれを足す、という操作をしているのだ。

 109: lea    0x4(%rcx),%rdx    ; rcxから32bit進んだアドレスをrdxへ
 10d: add    (%rcx),%eax       ; 返却値へ余った要素を足す
 10f: cmp    %rdx,%rsi         ; 先のrdxとendを比較
 112: je     12a <accum+0x12a> ; 一致していればreturn
 114: lea    0x8(%rcx),%rdx    ; rcxから2つ分進んだアドレスをrdxへ
 118: add    0x4(%rcx),%eax    ; 返却値へ余った要素を足す
 11b: cmp    %rdx,%rsi         ; rdxとendと比較
 11e: je     150 <accum+0x150> ; 一致していればreturn
 120: add    0x8(%rcx),%eax    ; 最大でも3つしか余らないので
 123: retq                     ; 足してreturn

ちなみに、中盤の同じ命令のパターンが続く部分も、基本的にはこれと同じことをしている。16バイト境界に整列したとみなせる部分まで、手動で足しているのだ。 しかし、それなら最大でも3つのはずで(int32_tは32bit、つまり4バイトなので、16バイトでアラインするまで4つは超えない)、7つある理由がわからない……と思っていたが、そんなことはない。 SIMDで4つずつ足すとしよう。手元に4つの要素がある。これをSIMDにロードして……それでどうする? 配列の要素は既に全てSIMDレジスタにロードされた。これ以上することはない。 あとは総和を取るだけだが、これは素でやってもそう変わらない。なので、SIMDにして高速化を期待するためには、8個以上の要素数が必要なのだ。だから最初に7個分の処理があるのだろう。

解読・後半戦

最後に、SIMDでループしている部分、先のアセンブリコードの直前の部分を読もう。

  c5: pxor   %xmm0,%xmm0        ; xmm0をゼロクリア
  c9: lea    (%rdi,%rdx,4),%rdi ; rdi = rdi + rdx * 4
  cd: xor    %edx,%edx          ; edxをゼロクリア
  cf: add    $0x1,%rdx          ; rdx += 1(カウンタインクリメント)
  d3: paddd  (%rdi),%xmm0       ; xmm0へrdiから128bit分足す
  d7: add    $0x10,%rdi         ; rdiを16byte前進
  db: cmp    %rdx,%r10          ; rdxをr10(多分4で割った余り)と比較
  de: ja     cf <accum+0xcf>    ; 終了していなければcfへ
  e0: movdqa %xmm0,%xmm1        ; xmm0をxmm1へコピー
  e4: lea    (%rcx,%r11,4),%rcx ; rcx += r11*4
  e8: psrldq $0x8,%xmm1         ; xmm1を8byteシフト
  ed: paddd  %xmm1,%xmm0        ; xmm1をxmm0に足す
  f1: movdqa %xmm0,%xmm1        ; xmm0をxmm1へコピー
  f5: psrldq $0x4,%xmm1         ; xmm1を4byteシフト
  fa: paddd  %xmm1,%xmm0        ; xmm1をxmm0に足す
  fe: movd   %xmm0,%edx         ; xmm0の下位ビットをedxへ
 102: add    %edx,%eax          ; edxをeaxへ足す
 104: cmp    %r11,%r8           ; r11とr8を比較
 107: je     12a <accum+0x12a>  ; 一致すればリターン

ここでは、メモリアラインが済んでいると考えてSIMDループを回した後、総和を取って返却値へ足している。総和はシフトしつつ足すという方法によって行われている。 8byteシフトすると、上2つのint32_tと下2つのint32_tがそれぞれ足されることになり、その後その値をまたシフトして足せば総和が求まる。 これを、余った部分を既に足していた返却値レジスタへ足し、この領域を抜けている。この後に続くのは最後に余った部分の総和だ。

bit単位で完全に理解しながら読んできたわけではないが、コンパイラが異様に賢いことはよくわかった。一応、-march=nativeや手動ループアンロールのコードも用意していたのだが、結構長くなったし、それ以上に私が疲れた。 今回はここまでで、今手持ちのストックは日を改めて解読していくことにする。

では。

(追記)

しかし、このコードだと要素数が非常に小さかったときはかえって遅くなるのでは、と考えた私は、7要素(SIMDにできない限界だ)の配列を作り、-O0-O3-O3 -march=nativeで比較してみた。すると、

$ gcc -std=c99 -O0 main.c accumulate.c
$ time ./a.out
0.22user 0.00system 0:00.22elapsed 98%CPU
$ gcc -std=c99 -O2 main.c accumulate.c
$ time ./a.out
0.08user 0.00system 0:00.08elapsed 98%CPU
$ gcc -std=c99 -O3 main.c accumulate.c
$ time ./a.out
0.10user 0.00system 0:00.10elapsed 99%CPU
$ gcc -std=c99 -O3 -march=native main.c accumulate.c
$ time ./a.out
0.07user 0.00system 0:00.07elapsed 98%CPU

妥当な結果が出た。O0だとスタック読み書きの効果が出て遅いが、-O2だとその分がなくなるので、SIMDの恩恵が受けられない場合はO3march=nativeとも十分張り合える速度になる。

もう少しよく見てみたいと思ったので、ループ回数を10倍にしたところ、

$ gcc -std=c99 -O2 main.c accumulate.c
$ time ./a.out
0.46user 0.00system 0:00.47elapsed 99%CPU
0.45user 0.00system 0:00.45elapsed 99%CPU
0.46user 0.00system 0:00.46elapsed 100%CPU
$ gcc -std=c99 -O3 main.c accumulate.c
$ time ./a.out
0.55user 0.00system 0:00.55elapsed 99%CPU
0.57user 0.00system 0:00.57elapsed 100%CPU
0.57user 0.00system 0:00.57elapsed 99%CPU
$ gcc -std=c99 -O3 -march=native main.c accumulate.c
$ time ./a.out
0.50user 0.00system 0:00.51elapsed 99%CPU
0.51user 0.00system 0:00.51elapsed 99%CPU
0.48user 0.00system 0:00.49elapsed 99%CPU

となった。若干ではあるが(そして3回ずつしか実行していないので統計的にはあやしいが)、O2が一番早くなっている。SIMDが使えない条件だと、メモリ配置や要素数の計算はオーバーヘッドでしかない。