分岐予測ミスペナルティを測ってみた

近年のほとんどのプロセッサには分岐予測機構が搭載されています。 分岐予測機構は、この命令は分岐命令か否か、分岐命令だとしたらとび先はどこか、条件分岐ならそれが成立するかしないか、などを予測します。 予測が外れていた場合、正しい命令のフェッチからやり直しになり、パイプラインが埋まるまでには時間がかかります。 このパイプラインが埋まるまでにかかる時間を分岐予測ミスペナルティとします。

実験環境

実験1: リターンアドレスを破壊してみる

以下に示した実験コードでは、スタックに保存されているリターンアドレスを勝手に書き替えています。 リターン命令で戻る先を関数コール命令の次ではないところにすることで、分岐予測を外します。

main:
        push rax
        mov  rcx, 1000000000
        call f
        pop  rax
        ret

f:
        sub rcx, 1
        jz L1
        sub QWORD PTR[rsp], 5
L1:
        ret

これを実行するのにかかった時間は10.8秒ほどでした。 3.1GHzで動いていることから、33G cycleほどかかったことがわかります。 1G回ループを実行していることから、分岐予測ミスペナルティは最大で33 cycle程度ということになります。

実際の分岐予測ミスペナルティは、ループの中にcall, sub, retといったメモリに触る命令があるため解析が難しいです。 メモリアクセスには4cycleかかるとし、メモリを介した依存が正しく予測されているとすれば、callのストアに4cycle、subのロードに4cycle、subの演算に1cycle、subのストアに4cycle、retのロードに4cycle、で分岐予測ミスペナルティ以外に17cycleかかることになります。よって、分岐予測ミスペナルティ自体は16cycleということになります。

実験2: 疑似乱数を用いて方向予測器を外してみる

以下に示した実験コードでは、疑似乱数(64bitのxorshift)の下1bitを使った条件分岐があります。 成立しても成立しなくても結局同じとび先なのですが、正しく予測できなければ分岐予測ミス扱いとなり、ペナルティが発生します。 疑似乱数の下1bitを予測することは難しいので、50%の確率で分岐予測ミスが発生します。

main:
        mov     rcx, 2000000000
        mov     rax, 1

L0:
        test    eax, 1
        je      L1
L1:

        mov     rdx, rax
        shl     rdx, 7
        xor     rdx, rax
        mov     rax, rdx
        shr     rax, 9
        xor     rax, rdx

        sub     rcx, 1
        jne     L0
        ret

このコードを実行すると、8.45秒ほどかかりました。 3.1GHzで動いていることから、プログラムの実行に26G cycleかかったことになります。

mov rax, 1の部分をmov rax, 0に変更すると、xorshiftの計算自体は行われるものの常に0が出力され、分岐予測が当たるようになります。 この場合には2.75秒かかりました。 3.1GHzで動いていることから、プログラムの実行に8G cycleかかったことになります。 ループ一周当たりのレイテンシが4cycleであり2G回のループであるため、予想通りの結果です。

これらの差分から、分岐予測ミスペナルティは合計で18G cycle程度であることがわかります。 分岐予測ミスが50%の確率で発生するので、分岐予測ミス回数は1G回です。 よって、分岐予測ミスペナルティは18 cycleということになります。

まとめ

Intel Core i5 7200Uの分岐予測ミスペナルティは16~18cycle程度のようです。

量子コンピュータの勉強(その1)

[1204.4221] Magic-state distillation with the four-qubit codeという論文のFig. 1(c)でやっていることを手計算で試して理解しました。

マジック状態とは

量子コンピュータの勉強をやると、Xゲート・Zゲート・Hゲート・CNOTゲート、などが出てきますが、これだけでは万能な量子計算はできません。

万能な量子計算をするためには、例えばブロッホ球上で45°回すみたいなゲートが追加であれば十分です。

しかし、そのようなゲートは、量子誤り訂正符号上で実現することが困難です。

この問題は、量子テレポーテーションの時と同様のうまい手順で解決することができます。

具体的には、「適切な状態」を事前に準備しておき、(簡単なゲート操作を行った後)適切な観測を行い、その観測結果によってゲートを通したり通さなかったりする、という手順です。

この時に必要となる「適切な状態」のことをマジック状態と呼びます。

以下の図では、マジック状態である \left|H\right\rangle = \cos\left(\frac\pi8\right)\left|0\right\rangle+\sin\left(\frac\pi8\right)\left|1\right\rangleを用いたY軸周り45°回転を行う計算(以下図)をやってみます。

f:id:lpha_z:20201101224330p:plain
マジック状態|H>を用いたY軸周り45°回転手順(arXiv:1204.4221 の図1より一部抜粋)

落とし穴

論文中に書いてありますが、単独のYゲート操作と、角度パラメータ付きの Y(\theta)ゲート操作は別物です。

 Y = \left( \begin{array}{cc} 0 & -i \\
 i & 0 \end{array} \right)
 Y(\theta) = e^{-iY\theta/2} = \left( \begin{array}{cc} \cos(\theta/2) & -\sin(\theta/2) \\
\sin(\theta/2) & \cos(\theta/2) \end{array} \right)

計算

上の量子ビットである \left|H\right\rangleを、一般的に \left|H\right\rangle = c\left|0\right\rangle+s\left|1\right\rangleと置いてみます。

下の量子ビット \left|H\right\rangle = a\left|0\right\rangle+b\left|1\right\rangleと置きます。

すると、初期状態は ca\left|0\right\rangle\left|0\right\rangle+sa\left|1\right\rangle\left|0\right\rangle+cb\left|0\right\rangle\left|1\right\rangle+sb\left|1\right\rangle\left|1\right\rangleとなります。

Control-Yゲートを通った後の状態は、 ca\left|0\right\rangle\left|0\right\rangle-isb\left|1\right\rangle\left|0\right\rangle+cb\left|0\right\rangle\left|1\right\rangle+isa\left|1\right\rangle\left|1\right\rangleになっています。

ここで、基底を計算基底( \left|0\right\rangle \left|1\right\rangleを用いた表現)からYの固有状態( \left|+y\right\rangle = \frac{ 1}{ \sqrt 2} \left|0\right\rangle + \frac{ i}{ \sqrt 2 }\left|1\right\rangle \left|-y\right\rangle = \frac{ 1}{ \sqrt 2} \left|0\right\rangle - \frac{ i}{ \sqrt 2 }\left|1\right\rangle)を使った基底に変換します。これは、Yで測定するとYの固有状態のどちらかがその係数の絶対値の二乗の確率で出てくるため、その計算をやりやすくするためです。

すると、 ca\left|0\right\rangle\left|0\right\rangle-isb\left|1\right\rangle\left|0\right\rangle+cb\left|0\right\rangle\left|1\right\rangle+isa\left|1\right\rangle\left|1\right\rangle  = \frac{ca - sb}{\sqrt 2}\left|+y\right\rangle\left|0\right\rangle+\frac{ca + sb}{\sqrt 2}\left|-y\right\rangle\left|0\right\rangle+\frac{cb + sa}{\sqrt 2}\left|+y\right\rangle\left|1\right\rangle+\frac{cb - sa}{\sqrt 2}\left|-y\right\rangle\left|1\right\rangleです。

これをYで測定して固有状態 \left|+y\right\rangleが出た場合(確率1/2)、下の量子ビット (ca - sb)\left|0\right\rangle+(cb + sa)\left|1\right\rangleになっています。これは、上の量子ビット \left|H\right\rangleだった場合、Y軸周りに45°回転したことになります。したがって、特に追加のゲート操作は不要です*1

一方、固有状態 \left|-y\right\rangleが出た場合(確率1/2)、下の量子ビット (ca + sb)\left|0\right\rangle+(cb - sa)\left|1\right\rangleになっています。これは、上の量子ビット \left|H\right\rangleだった場合、-Y軸周りに45°回転したことになります。つまり、Y軸周りに-45°回転したことになります。したがって、追加でY軸周りに90°のゲート操作を行えば、合計でY軸周りに45°回転したことになる、という具合です。

結局、どちらの固有状態が出た場合でも、その観測結果に応じて適切にゲートを通したり通さなかったりすることで、ブロッホ球上でのY軸周り45°回転を行うことができました。

ちなみに、二つの測定結果はどちらも半々の確率で出ることから、測定結果から何らかの情報を得ることはできず、下の量子ビットに関する情報は漏れていない(=下の量子ビットを破壊していない)ということも確認できます。

*1:論文ではpositive mesurementの場合にY(π/2)ゲートをかけるとしていましたが、逆だと思います。

Intel MKLのSIMD exp関数を読んだ

Intel CPUで高速に動く、数学関数のSIMD実装に Intel Math Kernel Library というものがあります。

それに含まれる指数関数実装を読んでみました。

AVX-512(512bitレジスタがあるので、doubleを最大で八並列で演算できる)用の実装は、以下のようになっていました*1

00000000004013d0 <__svml_exp8_z0>:
  4013d0:  f3 0f 1e fa                     repz nop     edx
  4013d4:  62 f1 7c 48 10 25 62 29 00 00   vmovups      zmm4 , [rip+0x2962]
  4013de:  62 f1 7c 48 10 0d 98 29 00 00   vmovups      zmm1 , [rip+0x2998]
  4013e8:  62 f1 7c 48 10 15 ce 29 00 00   vmovups      zmm2 , [rip+0x29ce]
  4013f2:  62 f1 7c 48 10 1d 04 2a 00 00   vmovups      zmm3 , [rip+0x2a04]
  4013fc:  62 f1 7c 48 10 2d 7a 2a 00 00   vmovups      zmm5 , [rip+0x2a7a]
  401406:  62 71 7c 48 10 2d b0 2a 00 00   vmovups      zmm13, [rip+0x2ab0]
  401410:  62 f2 fd 78 a8 e1               vfmadd213pd  zmm4 , zmm0 , zmm1 , {rz-sae}
  401416:  62 f1 7c 48 10 35 e0 2a 00 00   vmovups      zmm6 , [rip+0x2ae0]
  401420:  62 71 7c 48 10 0d 16 2b 00 00   vmovups      zmm9 , [rip+0x2b16]
  40142a:  62 71 7c 48 10 1d 4c 2b 00 00   vmovups      zmm11, [rip+0x2b4c]
  401434:  62 71 7c 48 10 25 02 28 00 00   vmovups      zmm12, [rip+0x2802]
  40143e:  62 71 dd 18 5c f1               vsubpd       zmm14, zmm4 , zmm1 , {rn-sae}
  401444:  62 72 dd 48 7f 25 32 28 00 00   vpermt2pd    zmm12, zmm4 , [rip+0x2832]
  40144e:  62 d2 ed 18 bc c6               vfnmadd231pd zmm0 , zmm2 , zmm14, {rn-sae}
  401454:  62 d2 e5 18 bc c6               vfnmadd231pd zmm0 , zmm3 , zmm14, {rn-sae}
  40145a:  62 f1 fd 48 54 3d dc 29 00 00   vandpd       zmm7 , zmm0 , [rip+0x29dc]
  401464:  62 71 c5 18 59 c7               vmulpd       zmm8 , zmm7 , zmm7 , {rn-sae}
  40146a:  62 72 d5 18 b8 ef               vfmadd231pd  zmm13, zmm5 , zmm7 , {rn-sae}
  401470:  62 72 cd 18 b8 cf               vfmadd231pd  zmm9 , zmm6 , zmm7 , {rn-sae}
  401476:  62 71 bd 18 59 d7               vmulpd       zmm10, zmm8 , zmm7 , {rn-sae}
  40147c:  62 72 bd 18 a8 df               vfmadd213pd  zmm11, zmm8 , zmm7 , {rn-sae}
  401482:  62 52 bd 18 a8 e9               vfmadd213pd  zmm13, zmm8 , zmm9 , {rn-sae}
  401488:  62 52 ad 18 a8 eb               vfmadd213pd  zmm13, zmm10, zmm11, {rn-sae}
  40148e:  62 52 9d 18 a8 ec               vfmadd213pd  zmm13, zmm12, zmm12, {rn-sae}
  401494:  62 d2 95 18 2c c6               vscalefpd    zmm0 , zmm13, zmm14, {rn-sae}
  40149a:  c3                              ret

この実装では、 \exp(x)を求めるとき、まず \frac{x}{\log 2}を1/16の位まで求めます。この値の整数部分は、最後のvscalefpdで使用されます。小数部分は、vpermt2pdを使った16エントリの表引きに使われます。

 x \frac{1}{\log 2}で割ったあまりを tとすると、絶対値は \frac{\log 2}{16}未満のはずなので、 \exp(t)を六次の多項式を使って近似します。

vandpdは、 xの絶対値が大きすぎて \frac{x}{\log 2}を1/16の位まで求める計算がおかしくなった時に tがおかしな値となり、最終結果がおかしくなることを防ぐ役割があるようです。 vandpdで使われるマスク定数は0xbfffffffffffffffであり、絶対値が2以上の数値やNaN±Infを絶対値が2未満の数に変更する効果があります。 xの絶対値が大きすぎる場合はどうせvscalefpd0+Infになるので、適当に小さくすればよいということのようです。

こうして見てみると、vperfmt2命令やvscalef命令は指数関数を実装するためにあるという感じの命令ですね……。

*1:マスクレジスタを使わないことが確定している場合の実装です。

Intel CPU の分岐予測器を調べてみた(ループの場合だけ)

なんだか分岐予測器で当てれそうなところを当ててくれなかったことがあったので、Intel CPUの分岐予測器はどんな時に分岐予測を当てることができるのかを調べてみました。

実際に調べてみると、ループのような単純な場合であっても非常に複雑怪奇な挙動を示すため、内部構造がどうなっているか全く推測できませんでした。

とりあえず、結果だけまとめておきます。

実験方法

単純なループを10億回くりかえすだけの以下のアセンブリコードを実行したときの分岐予測ミス回数を計測しました。

        .globl  main
main:
        xorl      %eax, %eax
L6:
        xorl      %edx, %edx
L7:
        incl      %edx
        cmpq      $__, %rdx
        jb        L7
        incl      %eax
        cmpl      $1000000000, %eax
        jb        L6
        ret

ループ回数は、以下のスクリプトで設定しました。

for i in $(seq 1 100); do
 sed s/__/$i/ main.s > m.s
 gcc m.s
 perf stat ./a.out
done

実験に使用したCPU

ここ10年くらいのXeonで実験することができました。

  • Xeon X5550(Neharem EP、2009年)
  • Xeon E5-2690(Sandy Bridge EP、2012年)
  • Xeon E5-2643 v3(Haswell EP、2014年)
  • Xeon E5-1620 v4(Broadwell EP、2016年)
  • Xeon Gold 5115(Skylake SP、2017年)
  • Xeon Gold 6126(Skylake SP、2017年)
  • Xeon platinum 8280(Canscadelake SP、2019年)

結果

Xeon X5550(Neharem EP、2009年)

  • ループが64周以下ならば、分岐予測ミスは発生しない。
  • ループが65周以上であれば、(おそらくループ脱出時に)分岐予測ミスが1回発生する。

分岐履歴長が63であることを示唆しているように思えます。

非常に素直な結果になりました。

Xeon E5-2690(Sandy Bridge EP、2012年)

  • ループが33周以下ならば、分岐予測ミスは発生しない。
  • ループが34周以上であれば、(おそらくループ脱出時に)分岐予測ミスが1回発生する。

分岐履歴長が32であることを示唆しているように思えます。

非常に素直な結果ですが、なぜ3年前のプロセッサよりも劣化しているのでしょうか。 ループ以外の分岐予測の精度が上がったと信じたいです。

Xeon E5-2643 v3(Haswell EP、2014年)

  • ループが35周、40周、42周、44周、の場合にはループ一周につき1回以上の分岐予測ミスが発生する。
  • ループが53周、54周、55周、60周~74周、91周~94周の場合はループ一周につき1回未満だが0.1回以上の分岐予測ミスが発生する。
  • それ以外の場合でループが59周以下であれば、分岐予測ミスはほとんど発生しない。
  • ループが75周~90周の場合、ループ一周につき1.2回程度の分岐予測ミスが発生する。
  • ループが95周以上であれば、(おそらくループ脱出時に)分岐予測ミスが1回発生する。

わけの分からない結果になりました。

Xeon E5-1620 v4(Broadwell EP、2016年)

やるたびに結果が変わってよくわかりませんが、

  • ループが35周、40周、42周、44周の場合はループ一周につき1回程度の分岐予測ミスが発生しやすい
  • ループが40周以上だと当てられなくなる時と、60周以上だと当てられなくなる時がある

ループ35周、40周、42周、44周の時に当てられないという挙動が、Xeon E5-2643 v3(Haswell EP、2014年)と酷似しています。

Xeon Gold 5115(Skylake SP、2017年)と Xeon Gold 6126(Skylake SP、2017年)

やはりやるたびに結果が変わってよくわかりませんが、

  • ループが35周未満であれば、分岐予測ミスはほとんど発生しない
  • ループが35周のとき、分岐予測ミスが発生しやすい

やはり35周のところに何かがあるようです。

Xeon platinum 8280(Canscadelake SP、2019年)

  • ループが23周以下ならば、分岐予測ミスは発生しない。
  • ループが24周~28周の場合はループ一周につき1回未満だが0.05回以上の分岐予測ミスが発生する。
  • ループが29周以上であれば、(おそらくループ脱出時に)分岐予測ミスが1回発生する。

素直な結果に戻りましたが、7年前のプロセッサよりもさらに劣化しています。 "improved branch predictor" とはいったい何だったのでしょう……

まとめ

  • Intel CPU にはループ予測器はついていなさそう
  • ループ脱出の予測に関しては、Intel CPUはここ10年で全く進化していない(どころか当てられるループ周回数が減少している)
  • 30周程度のループでも当てられないことがある


  • 35周の時に当てられない現象は、分岐命令のアドレスなどを変えて再検証する必要あり

「単純和」プログラムの速度をJuliaとC言語で比較した

先週、浮動小数点数を足していくプログラムについて、Juliaで書いた場合とC言語で書いた場合とで速度を比較しました(浮動小数点数を足すだけのプログラムがC言語よりJuliaのほうが速い、について - よーる)。

この記事について、黒木玄さん(黒木玄 Gen Kuroki (@genkuroki) | Twitter)から「Juliaでは単純な和には普通にsimdを使うので、JuliaとCの両方でsimdを使って比較した方が良かったと思います。」というコメントをいただきました。 先週の記事で比較に用いたプログラムには並列性がないため、SIMD (single instruction multiple data) 命令を使った高速化は不可能なはずですが、Juliaでsimdを使うと10倍高速化するそうです。 一体どういうことなのだろうと思ってJuliaのマニュアル(Base.SimdLoop.@simd (macro))を見てみると、@simdというアノテーション浮動小数点数演算の順序を変える効果があるそうです。

浮動小数点数演算は一般に結合則が成り立ちません。 したがって、浮動小数点演算の順序を変えてしまうと、プログラムの意味が変わってしまいます(一般に計算結果も変わってしまいます)。 先週の記事で比較したプログラムはあくまで「浮動小数点数を前から順番に足していく」プログラムであって、「いい感じに浮動小数点数の和を求める」プログラムではありません。

この辺は、多くの読者が勘違いされているところだと感じました。

とはいえ、「浮動小数点数を前から順番に足していった値を正確に求めたい」というワークロードが常識的ではないのも事実です。 そこで、「浮動小数点数を順番を気にせずに足した値を求める」というワークロードで、JuliaとC言語の速度比較をしてみます。

使ったソースコード

ちょうど黒木玄さんがJulia側のソースコードを示してくださいました(simd example · GitHub)ので、それを使うことにします*1。 ただし、実行時間が短すぎて正確な評価ができなかったので、外側ループの回数を100倍にしました。

C言語側のソースコードは、私が素直に移植したものになります。

function test_inbounds_simd()
    a = fill(1e-10, 1000)
    s = 0.0
    for j = 1:1000000000
        @inbounds @simd for i = 1:1000
            s += a[i]
        end 
    end
    s
end
test_inbounds_simd()
double test() {
  unsigned long long i, j;
  double a[1000];
  double s = 0.0;
  for( i = 0; i < 1000; ++i ) { a[i] = 1e-10; }

  for( j = 0; j < 1000000000; ++j ) {
    for( i = 0; i < 1000; ++i ) {
      s += a[i];
    }
  }

  return s;
}

int main() {
  return test();
}

使用したハードウェア

  • Intel(R) Xeon(R) Platinum 8280 CPU @ 2.70GHz
    • 実際には3.6GHzくらいで動いているっぽい(?)
    • AVX-512使用可能(一命令で倍精度浮動小数点数8個を同時に演算可能)
    • FMA実行ユニット2基(1cycleに浮動小数点数加算命令が2回実行できる)
    • 浮動小数点数加算命令のレイテンシは4cycle

使用したソフトウェア

  • julia version 1.5.2(内部で使用しているLLVMのバージョンはLLVM version 9.0.1jl
  • clang version 10.0.0
  • icc (ICC) 19.0.5.281 20190815
  • perf
  • gdb

理論限界

ベンチマークプログラムは1000G個の浮動小数点数を加算するものです。 実験に使用したプロセッサは1サイクルに16個の浮動小数点数を加算することができますので、理論限界は62.5Gサイクルです。 3.6GHzで動いているとすると17秒程度かかることになりますが、周波数が安定しないので実時間は参考程度ということにします。

Juliaの結果

  • 199G cycles
  • 278G insns
  • 43G branches
  • 55秒くらい

JITコンパイルされた機械語gdbで確認してみると、以下のようになっていました。

   0x00007fffc565f0d0:  cmp    rax,0x3b9aca00
   0x00007fffc565f0d6:  lea    rax,[rax+0x1]
   0x00007fffc565f0da:  je     0x7fffc565f187
   0x00007fffc565f0e0:  vmovq  xmm0,xmm0
   0x00007fffc565f0e4:  vxorpd xmm1,xmm1,xmm1
   0x00007fffc565f0e8:  mov    esi,0x18
   0x00007fffc565f0ed:  vxorpd xmm2,xmm2,xmm2
   0x00007fffc565f0f1:  vxorpd xmm3,xmm3,xmm3
   0x00007fffc565f0f5:  nop    WORD PTR cs:[rax+rax*1+0x0]
   0x00007fffc565f0ff:  nop
   0x00007fffc565f100:  vaddpd zmm0,zmm0,ZMMWORD PTR [rcx+rsi*8-0xc0]
   0x00007fffc565f108:  vaddpd zmm1,zmm1,ZMMWORD PTR [rcx+rsi*8-0x80]
   0x00007fffc565f110:  vaddpd zmm2,zmm2,ZMMWORD PTR [rcx+rsi*8-0x40]
   0x00007fffc565f118:  vaddpd zmm3,zmm3,ZMMWORD PTR [rcx+rsi*8]
   0x00007fffc565f11f:  add    rsi,0x20
   0x00007fffc565f123:  cmp    rsi,0x3f8
   0x00007fffc565f12a:  jne    0x7fffc565f100
   0x00007fffc565f12c:  vaddpd zmm0,zmm1,zmm0
   0x00007fffc565f132:  vaddpd zmm0,zmm2,zmm0
   0x00007fffc565f138:  vaddpd zmm0,zmm3,zmm0
   0x00007fffc565f13e:  vextractf64x4 ymm1,zmm0,0x1
   0x00007fffc565f145:  vaddpd zmm0,zmm0,zmm1
   0x00007fffc565f14b:  vextractf128 xmm1,ymm0,0x1
   0x00007fffc565f151:  vaddpd xmm0,xmm0,xmm1
   0x00007fffc565f155:  vpermilpd xmm1,xmm0,0x1
   0x00007fffc565f15b:  vaddsd xmm0,xmm0,xmm1
   0x00007fffc565f15f:  test   dl,dl
   0x00007fffc565f161:  jne    0x7fffc565f0d0
   0x00007fffc565f167:  xor    esi,esi
   0x00007fffc565f169:  nop    DWORD PTR [rax+0x0]
   0x00007fffc565f170:  vaddsd xmm0,xmm0,QWORD PTR [rcx+rsi*8+0x1f00]
   0x00007fffc565f179:  inc    rsi
   0x00007fffc565f17c:  cmp    rsi,0x8
   0x00007fffc565f180:  jne    0x7fffc565f170
   0x00007fffc565f182:  jmp    0x7fffc565f0d0

この機械語コードから予測すると、1000個の浮動小数点数を足すのに

  • 192cycle(zmmレジスタを使った主ループが4cycle×31回、マージしていく部分が36cycle、端数処理ループが4cycle×8回)
  • 273命令(準備に11命令、主ループが7命令×31回、マージしていく部分が13命令、端数処理ループが4命令×8回)
  • 42分岐命令

が必要です。

これはperfコマンドの結果と整合します。

理論限界と大きく差がついている理由は、以下の通りです。

  • 主ループのアンローリングが不足している
    • 4cycleレイテンシを隠蔽しつつFMAユニット2基をフル稼働させるためには、依存関係にない命令が最低でも8つ必要
  • 端数処理が直列化されている

C言語の結果(clangでコンパイルした場合)

コンパイルオプションは、-O2 -mavx512f -funsafe-math-optimizationsです*2C言語の範囲では「浮動小数点数演算の順序を自由に変更してかまわない」という記述をすることができないので、-funsafe-math-optimizationsをつけます*3*4

  • 194G cycles
  • 213G insns
  • 32G branches
  • 55秒くらい

Juliaで書いた場合と同じくらいの速度になりました。それもそのはずで、JuliaのJITエンジンはLLVM、つまりclangのバックエンドと同じです。 というわけでコード生成結果や実行速度が大きく変わるはずはないのです。

命令数が少なくなっているのは、JITコンパイルより時間をかけてコンパイルすることできれいなコードを生成できたことによるものと思われます。

.LBB0_5:
        vaddpd  zmm0, zmm2, zmm0
        vaddpd  zmm0, zmm3, zmm0
        vaddpd  zmm0, zmm4, zmm0
        vextractf64x4   ymm2, zmm0, 1
        vaddpd  zmm0, zmm0, zmm2
        vextractf128    xmm2, ymm0, 1
        vaddpd  xmm0, xmm0, xmm2
        vpermilpd       xmm2, xmm0, 1
        vaddsd  xmm0, xmm0, xmm2
        vaddsd  xmm0, xmm0, xmm1
        vaddsd  xmm0, xmm0, xmm1
        vaddsd  xmm0, xmm0, xmm1
        vaddsd  xmm0, xmm0, xmm1
        vaddsd  xmm0, xmm0, xmm1
        vaddsd  xmm0, xmm0, xmm1
        vaddsd  xmm0, xmm0, xmm1
        vaddsd  xmm0, xmm0, xmm1
        add     rax, 1
        cmp     rax, 1000000000
        je      .LBB0_6
.LBB0_3:
        vmovq   xmm0, xmm0
        vxorpd  xmm2, xmm2, xmm2
        mov     ecx, 56
        vxorpd  xmm3, xmm3, xmm3
        vxorpd  xmm4, xmm4, xmm4
        .p2align        4, 0x90
.LBB0_4:
        vaddpd  zmm0, zmm0, zmmword ptr [rsp + 8*rcx - 576]
        vaddpd  zmm2, zmm2, zmmword ptr [rsp + 8*rcx - 512]
        vaddpd  zmm3, zmm3, zmmword ptr [rsp + 8*rcx - 448]
        vaddpd  zmm4, zmm4, zmmword ptr [rsp + 8*rcx - 384]
        cmp     rcx, 1016
        je      .LBB0_5
        vaddpd  zmm0, zmm0, zmmword ptr [rsp + 8*rcx - 320]
        vaddpd  zmm2, zmm2, zmmword ptr [rsp + 8*rcx - 256]
        vaddpd  zmm3, zmm3, zmmword ptr [rsp + 8*rcx - 192]
        vaddpd  zmm4, zmm4, zmmword ptr [rsp + 8*rcx - 128]
        add     rcx, 64
        jmp     .LBB0_4

zmmレジスタを使った主ループが二倍アンローリングされたほか、端数ループが八倍(完全)アンローリングされています。これらが命令数削減に役立っているようです。

この機械語コードから予測すると、1000個の浮動小数点数を足すのに

  • 192cycle(zmmレジスタを使った主ループが8cycle×15.5回、マージしていく部分が36cycle、端数処理が4cycle×8回)
  • 213命令(ループの外が27命令、ループの中が12命令×15.5回)
  • 32分岐命令

が必要です。

これはperfコマンドの結果と整合します。

理論限界と大きく差がついている理由は、Juliaの場合と全く同様です。

どうでもいいですが、配列の最後のほうがループ内定数とみなされてxmmレジスタに乗っていますね……。これは微妙にずるな気がしますが計算時間や命令数に変わりはないのでよしとしましょう。

C言語の本気(iccコンパイルした場合)

C言語の強みは何といっても ハードウェアを作っている会社が専用コンパイラを提供してくれる ことにあります。

Intel社が作ったコンパイラiccを使ってコンパイルしてみましょう。 コンパイルオプションは、-axCORE-AVX512 -qopt-zmm-usage=highです。

結果は以下の通りです。

  • 62.7G cycles
  • 172G insns
  • 15.6G branches
  • 18秒くらい

理論限界の62.5Gサイクルにかなり近い速度が達成できていることがわかります*5

最内ループを確認してみると、お互いに依存しない8個のvaddpd命令になっていました。 さすがにこの辺の最適化はちゃんとしていますね。

まとめ

*1:この手の言語間の速度比較では、不適切なソースコードが使われることにより特定の言語の性能が不当に低く出てしまうということがありがちですが、そういうことにならない有意義な比較ができることを感謝いたします。

*2:-O3や-Ofastだと完全にアンローリングされてしまって見づらかったので-O2にしました(実行にかかるサイクル数は同じでした)。

*3:-funsafe-math-optimizationsは-ffast-mathに含まれるので、-ffast-mathでもよいです。

*4:-funsafe-math-optimizationsは複数のオプションを有効にしますが、実際に役立っているのは-fassociative-math -fno-signed-zerosのようです。

*5:ループ交換が行われているので本来とは全く異なる順序で計算されていますが、加算は1000G回行われています。

浮動小数点数を足すだけのプログラムがC言語よりJuliaのほうが速い、について

なんかすごく前に話題になっていたので、確かめてみました。

数値計算に強いプログラミング言語と言えば従来FortranC言語でしたが、近年Pythonのように書けて実用上十分な速度を達成できるJuliaが人気を集めているようです*1

Juliaは動的言語ですが、(1) 型推論を行うためfor文などでいちいち型検査するPythonより速い (2) JITコンパイルされるので関数に切り出した場合に高速、といった点が高速化に寄与しています。

C言語で書いたほうが速いに決まっている」などの決めつけはよくないので、実際に計測してみます。

そもそもベンチマークソースコードが提示されていないので議論のしようがないのですが、「浮動小数点数を足すだけ」という文章を私が解釈したコードは以下の通りになります。

function test()
  a = zeros(Float64, 1000)
  a[1] = 1e-16

  sum = 1.0

  for j = 1:10000000, i = 1:1000
    sum += a[i]
  end

  return sum
end

test()
double test() {
  double a[1000] = {};
  a[0] = 1e-16;

  double sum = 1.0;

  for( unsigned long long j = 0; j < 10000000; ++j ) {
    for( unsigned long long i = 0; i < 1000; ++i ) {
      sum += a[i];
    }
  }

  return sum;
}

int main() {
  return test();
}

メモリアクセスの速度が測りたいわけではないので、配列のサイズはL1キャッシュに乗る程度の大きさとしました。

環境

  • gcc 9.1.0、-O3をつけてコンパイル
  • julia version 1.5.2
  • Intel(R) Xeon(R) CPU E5-2603 v3 @ 1.60GHz(addsd命令は3cycleレイテンシ)

結果

C言語

  • 30.1G cycles
  • 35.1G instructions
  • 5.0G branches
  • 18.8 seconds

Julia

  • 31.9G cycles
  • 61.5G instructions
  • 20.3G branches
  • 19.1 seconds

考察

数値上はJuliaのほうがほんのわずかに遅くなっていますが、C言語は事前コンパイルが必要でそれに0.2秒ほどかかることを考慮に入れると速度は同じと言ってよいでしょう。

そもそもこのベンチマークの速度はレイテンシに支配されているため、それなりの最適化がかかっていれば速度が変わるわけはないのです。

C言語コードをコンパイルした結果は以下のようになっており、ループ一周は7命令です。アセンブリ上でのループ一周は元のソースコードの二周分に相当しています(アンローリングされています)。 この中には直列に依存したaddsd命令が2つあるので、ループ一周当たり6cycleかかることになります。 C言語コンパイルした結果のperfのサイクル数、命令数、分岐命令数は、これらと整合します。

.L3:
        movsd   xmm1, QWORD PTR [rax]
        add     rax, 16
        addsd   xmm1, xmm0
        movsd   xmm0, QWORD PTR [rax-8]
        addsd   xmm0, xmm1
        cmp     rdx, rax
        jne     .L3

一方Juliaの方はサイクル数こそほぼ一緒ですが、命令数がかなり多くなっています。 これは配列の境界検査を毎回行っているためです。 C言語と条件をそろえ、配列境界検査を行わないようにするため、@inboundsをつけてみると、以下のように命令数と分岐命令数が減ります。

  • 31.9G cycles
  • 41.5G instructions
  • 10.3G branches
  • 19.1 seconds

命令数は一ループ当たり2命令、分岐命令数は一ループ当たり1命令減りました。この減った分は、比較命令と分岐命令ということになります。

分岐命令数から察するにループアンローリングは行われていないようですが、Juliaで書いたコードをJITコンパイルした結果はほぼ最適です。

まとめ

浮動小数点数を足すだけのプログラム」の速度がC言語で書いた場合とJuliaで書いた場合で変わることはありませんでした。

*1:ほかの言語では数値計算ができないとか遅いとか言っているわけではありませんが、あまり見かけないですね。最近はC++で書かれたソフトウェアもちらほら見かけ始めました(知っているところではALPSとか)。

パソコンの修理がおわった(二日連続・二回目)

同時に修理に出していたもう一つのパソコンの修理が終わったようです。戻ってくるときも同じとはちょっと不思議です。