2021-01-17

Unit of Last Place (ULP) の定義について

On the definition of ulp(x) という論文を読んだメモです。

浮動小数点数の誤差を議論するとき、最終桁の重み (Unit of Last Place, ULP) が単位として使われます。例えば、以下のような使い方をします。

最近接丸めをした場合、その丸め誤差は0.5ULP以下になる
方向丸め（切り捨てや切り上げ、零への丸めなど）を行った場合、その丸め誤差は1.0ULP未満になる

しかし、以下の倍精度二進浮動小数点数の例を考えてみると、丸め方によってULPの定義が変わるようでは困ることがわかります。

0x1.fffffffffffff8p+0を2.0に丸めた場合、丸め誤差は2^-53、2.0の最終桁の重みは2^-51なので、丸め誤差は0.25ULP（？）
0x1.fffffffffffff8p+0を0x1.fffffffffffffp+0に丸めた場合、丸め誤差は2^-53、0x1.fffffffffffffの最終桁の重みは2^-52なので、丸め誤差は0.5ULP（？）

したがって、「丸める前の値」のみに基づいてULPが定義されるべきであることがわかります。しかし、「丸める前の値」は浮動小数点数ではないので、「最終桁の重み」というのを正しく定義することができません。というわけで、やはり困ってしまいます。

そういう場合にも適応できるULPの定義は、下のようなものが挙げられます。

「丸める前の値」を挟む二つの浮動小数点数の距離を1ULPとする
「丸める前の値」に最も近い二つの浮動小数点数の距離を1ULPとする

1.の定義は、以下の欠点を抱えています。

「丸める前の値」がちょうど浮動小数点数で表せる場合に困る
- 実際に困るのは、基数のべき乗ぴったりの時
挟む二つの浮動小数点数の片方が±∞の時に困る

2.の定義は、実数の全領域で矛盾なく定義できるという利点がありますが、基数のべき乗に近いところで「最も近い二つの浮動小数点数」が「丸める前の値」を挟まない（両方とも「丸める前の値」よりも絶対値が小さい）ことがあり、やや直観に反します。

1.の定義と2.の定義は、浮動小数点数で表せるほとんどの領域上で同じ定義となりますが、オーバーフロー付近、アンダーフロー付近、基数のべき乗をわずかに超えた付近のみで異なる定義となります。浮動小数点数をオーバーフロー付近やアンダーフロー付近で使うべきではないですが、基数のべき乗をわずかに超えた付近は使わざるを得ないことがあります。よって、以下では基数のべき乗をわずかに超えた付近での問題を議論します。

基数のべき乗を超えた付近で定義が異なるというのは、例えば真の値が0x1.00000000000003p+0で倍精度二進浮動小数点数に丸める場合を考えてみるとわかります。

真の値を挟む二つの浮動小数点数は、0x1.0000000000000p+0と0x1.0000000000001p+0なので、1.の定義では1ULP=2^-52
真の値に最も近い二つの浮動小数点数は、0x1.fffffffffffffp-1と0x1.0000000000000p+0なので、2.の定義では1ULP=2^-53

1.の定義では、基数のべき乗のところでULPの大きさが変わりますが、2.の定義では、基数のべき乗を少し超えた付近でULPの大きさが変わります。

この切り替わる点は $r^N x_{cut}$ と表せます。倍精度二進浮動小数点数の場合、1.の定義では $r=2, x_{cut}=1$ 、2.の定義では $r=2, x_{cut}=1+2^{-54}=$ 0x1.0000000000004p+0となります。

これを一般化すると、浮動小数点数の基数をr、仮数部はp桁*1だとして、

1.の定義の場合、 $x_{cut} = 1$
2.の定義の場合、 $x_{cut} = 1 + \frac{r -1}2 r^{-p}$

となります。

ところで、ULPの定義にこだわる理由は、「丸め誤差がxxULP以内だから、○○な性質を持つ浮動小数点数に丸められているはず」「○○な性質を持つ浮動小数点数に丸めたから、丸め誤差がxxULP以内のはず」という議論をしたいからです。以下では、これらの性質が成り立つかを議論します。

丸め誤差が0.5ULP未満だから、唯一の最も近い浮動小数点数に丸めたはず

$1 + \left(\frac{r}{2}-1\right)r^{-p} \le x_{cut}$ を満たしていれば、上の性質が成り立ちます。

二進浮動小数点数であれば、1.の定義でも2.の定義でも条件を満たすため、上の性質が成り立ちます（1.の定義はぎりぎりです）。
三以上の基数の場合、2.の定義の場合のみ上の性質が成り立ちます。

最も近い浮動小数点数（のうちの一つ）に丸めたのだから、丸め誤差が0.5ULP以下のはず

$x_{cut} \le 1 + \frac 1 2 r^{-p}$ を満たしていれば、上の性質が成り立ちます。

二進浮動小数点数であれば、1.の定義でも2.の定義でも条件を満たすため、上の性質が成り立ちます（2.の定義はぎりぎりです）。
三以上の基数の場合、1.の定義の場合のみ上の性質が成り立ちます。

なお、三進浮動小数点数であっても、 $x_{cut} = \frac 1 2 3^{-p}$ とすれば、上の二つの性質を満たすULP関数を作ることができます。四進以上の場合、上の二つの性質を同時に満たすULP関数を作ることは不可能です。

切り上げで丸めても切り捨てで丸めても、丸め誤差は1.0ULP未満のはず

基数がいくつであっても、1.の定義の場合のみ上の性質が成り立ちます。

丸め誤差が1.0ULP未満なら、切り上げで丸めた数か切り捨てで丸めた数のはず

1.の定義の場合、成り立ちません。基数のべき乗をちょっと超えた数から1.0ULP未満の浮動小数点数は、r+1個あります。 2.の定義なら成り立つ気がしますが、記載されていませんでした。

結論

1.の定義がわかりやすく*2、二進浮動小数点数の場合には求められるすべての性質を満たすためこれを選ぶべきです。

ただし、無限大付近で問題になるため、絶対値が最大の浮動小数点数より大きい範囲では、2.の定義を流用するのが簡単です。 IEEE754で「こういう場合に無限大に丸めるべき」と定められたとおりに丸めれば、「最近接丸めを行った場合、その丸められた値が有限であれば、丸め誤差は0.5ULP以下」の性質は満たします。

*1:ケチ表現を使う二進浮動小数点数の場合、その1bitを含みます。倍精度二進浮動小数点数なら53桁と数えます。

*2:2.の定義の場合、「実数xを浮動小数点数に最近接丸めで丸めた場合のULPでみた丸め誤差」という関数が不連続になって取り扱いが不便です。

2021-01-03

CVPシミュレータを解読したメモ

新年になりました。今年もよろしくお願いします。

Value prediction（値予測）の世界的な大会、Second Championship Value Prediction @ HPCA'21というのが開催されているそうです。

値予測というのは投機的実行の一種で、レジスタに書き込まれる値を予測するものです。投機的実行の代表的な例として分岐予測がありますが、これは分岐命令のとび先を予測するものです。

どちらも、あくまで予測による投機的実行なので、実際に実行して予測が正しかったことを確認する必要があります。予測が間違っていた場合、その予測結果を（直接的にでも間接的にでも）使った命令を全てやり直す必要があります。特に、値予測の場合は分岐予測と違い、「なんでもいいから適当に予測しておく」のではなく「予測しない」という消極的選択肢を選んだ方が性能を高めることが多いです*1。

そういうわけで、値予測器の評価では（正答率やカバー率といった単純な指標ではなく、）実際にプログラムが高速に動作するかで評価を行う必要があります。そのために作られたシミュレータがCVP simulator（https://github.com/eric-rotenberg/CVP）です。

CVPシミュレータの特徴

CVPシミュレータの特徴は、以下の通りです。

Qualcomm社が提供するARM64のトレースを読み込んでシミュレーションする
- ARM64の命令は最大5つ*2のマイクロ命令に分解されてシミュレーションされる
（サイクル単位ではなく）命令単位でどのサイクルに何が起こるかをシミュレーションする
L1Iキャッシュをミスしない限り命令フェッチに時間がかからないとしている
分岐予測器にTAGE-SC-LとITTAGEを採用している
- BTBやRASは実装されていません

CVPシミュレータは、その構造上、以下の本質的な問題が存在します。

シミュレーションが正確ではない

サイクル順ではなくて命令順にシミュレーションしているため、（少なくとも）キャッシュの状態が正しくシミュレーションされません。具体的には、以下のようなパターンが問題になります。

あるロード命令由来のアクセスが、L1Dキャッシュミスを起こした後L2キャッシュでヒットする。
その次の命令が、L1Iキャッシュミスを起こした後L2キャッシュもミス。L3キャッシュのデータをL2キャッシュに挿入する。

このような場合、「L2キャッシュでヒットする」は「L2キャッシュに挿入する」の後に起こるかもしれず、その場合にはL2キャッシュの状態が変化しているので本当にヒットとは限りません。

命令キャッシュミスはそんなに発生しないので問題ないということなのでしょうか。しかし、配布されているトレースのうち一部のトレースでは23%も命令キャッシュミスが起きるらしい*3ので、「命令キャッシュミスはそんなに発生しない」という仮定は成り立っていないと思います。

とはいえ、サイクル順のシミュレーションをする場合はたくさんのコードを書かないといけないので、シミュレーションが多少不正確になっても問題ないとすれば、適当な妥協点なのかもしれません。

CVPシミュレータのおかしな点

CVPシミュレータの微妙な点（おそらく簡単に修正可能）は以下の通りです。

マイクロ命令への分解がナイーブ

ベースアドレス書き替え付きのロード命令がマイクロ命令に分解される場合、「ベースアドレス書き替えを実現する整数命令」と「ロード命令」に分けるべきですが、全てがロード命令に分解されます。意図していないアドレスの読み出しが行われるという問題もあります（aperaisさんも指摘しています→Aperais/cvp2v2 by aperais · Pull Request #3 · eric-rotenberg/CVP · GitHub）。

ベースアドレス更新のレイテンシに支配されているようにシミュレーションされているトレースでストライド予測を行ったらスコアが上がるのかもしれません。実際にはシミュレータの都合で不当に遅くなっているのを高速化しただけで、意味はないのですが……。

その他、フラグレジスタを介した依存関係を正しくシミュレーションできていません（トレースに記録されていないデータが必要なため、シミュレータで推測を行っているようですが、それが使われていません）。

ロード命令のアドレス変換には1cycleかかるのにストア命令のアドレス変換は0cycleで終わる

ロード命令の方はアドレス変換にかかる時間をシミュレーションしていますが、ストア命令の方はやっていません。

典型的にはストアアドレスの計算の方がデータ到着より早く終わるため、アドレス変換にかかる時間をカウントしていないのかもしれません。単に忘れていただけかもしれません。

ストアフォワーディングが0cycleで終わる

ストアを行ったサイクルにストアフォワーディングを待っているロード命令も終わります。

レジスタファイルの書き込みポート数をシミュレーションしていない

一サイクルに実行を開始できる命令の数（演算器の数）は正しくシミュレーションしていますが、一サイクルに実行を終了できる命令の数（レジスタファイルの書き込みポート数）はシミュレーションされていません。

固定レイテンシの命令だけであれば問題なさそうですが、キャッシュミスしたロード命令やストアフォワーディング待ちのロード命令は可変レイテンシであり、一サイクルに多数の命令が同時に終了します。

これは他のものと比べて実装難易度がやや高いですが、リプレイすべきサイクルを計算してリソース占有の予約を入れればシミュレーションできるはずです。

ここまで変な点をいくつか書きましたが、ChampSim（特にdevelopブランチは完全に壊れてしまっています）に比べたら全然ましです。

*1:実際、 https://www.microarch.org/cvp1/papers/Seznec_Unlimited_2020.pdfによれば、チューニングされた値予測器の正答率は99%を超えます。これは、プログラムの挙動が予測しやすいということではありません。値予測器には「予測しない」という選択肢が存在し、それを選ぶことが性能を高めるため、そうチューニングされているからです。

*2:SIMD レジスタを二つ読み込むload pair命令で、ベースレジスタ書き換え付きの場合に発生します。

*3:The 1st instruction prefetching championshipで使ったトレースはCVPで配布されているトレースを流用しているらしく、それを用いた評価（https://research.ece.ncsu.edu/ipc/wp-content/uploads/2020/05/eip.pdf）ではプリフェッチなしの場合23%も命令キャッシュミスが起こるとしています。

2020-12-20

CUDA libraryのexp(x)を読んだ

なんだか指数関数の実装ばかり見ているような気がしますが、GPGPU向けの実装はどうなっているんだろうなと思ったので見てみました。

Programming Guide :: CUDA Toolkit Documentationによれば、最大誤差1ULPを保証する実装になっているようです。 ULPの定義について詳しく語った論文を引用しているあたり、かなり丁寧なつくりになっていそうです。

中身を見てみたところ、11次の多項式近似をホーナー法で求め指数部をビット演算で操作する、というよくある感じの実装でした。

この方法でうまくいかない入力だった場合の例外処理が単精度浮動小数点数を使って実装されており面白かったです。これは、倍精度浮動小数点数演算器をなるべく使わないような実装をしているためだと考えられます。倍精度指数関数を計算する場合、倍精度浮動小数点数演算が大量に出てきますが、GPU上の倍精度浮動小数点数演算器は貴重なので、例外処理なんかに使いたくないということでしょう。

指数関数の入力は、大きく分けて以下の5つに分類されます*1。

小さすぎてexp(x)が0.0になるような入力（-Infを含む）
exp(x)が非正規化数になるような入力
exp(x)が正規化数になるような入力
大きすぎてexp(x)が+Infになるような入力（+Infを含む）
NaN（exp(NaN)=NaN）

この場合分けを、入力倍精度浮動小数点数の上位32bitを単精度浮動小数点数とみなして絶対値をとってから特定の定数と比較することで行っています。

「倍精度浮動小数点数の上位32bitを単精度浮動小数点数とみなして」というのがいきなり意味不明ですが、これはGPUで倍精度浮動小数点数を持つときのフォーマットを利用した技です。 GPU上で倍精度浮動小数点数を持つときには、隣り合う*232bitレジスタ二つに入れているようなので、これらのレジスタのうち片方を命令オペランドに指定することで、簡単に上位32bitを単精度浮動小数点数として扱えるということです。

倍精度浮動小数点数と単精度浮動小数点数だと指数部のbit数が違う点が気になるかもしれませんが、比較の両辺が正（+0.0含む、NaNを含まない）であればただの辞書順比較と同じなので大丈夫です。

そんなわけで、「元の倍精度浮動小数点数入力の上位32bitを単精度浮動小数点数とみなした時、その絶対値が特定の定数より小さければ、必ずexp(x)が正規化数になる」ように、「特定の定数」が選ばれているようです。実用上は大半の入力でこの条件が満たされるはずで、その場合にearly returnすれば追加の計算コストをほとんど払わずに済みます（この場合、単精度浮動小数点数の絶対値演算と比較演算がそれぞれ一回ずつで済みました）。

また、『元の倍精度浮動小数点数入力の上位32bitを単精度浮動小数点数とみなした時、その絶対値が別の特定の定数より大きい、またはその単精度浮動小数点数がNaNである*3場合、「小さすぎてexp(x)が0.0」「大きすぎてexp(x)が+Inf」「入力がNaNなのでexp(x)もNaN」のどれか』になるように「別の特定の定数」が選ばれているようです。このようなケースは例外ケースの中では頻出なので、ここも最適化テクニックで高速化されていました。具体的には、「exp(x)=x+(+Inf)を返す」という方法で「大きすぎてexp(x)が+Inf」と「入力がNaNなのでexp(x)がNaN」のケースをまとめて処理していました。どんな有限のxに対してもx+(+Inf)は+Infですし、xがNaNであればx+(+Inf)はNaNです*4。

それ以外のパターン、つまり「exp(x)が非正規化数になるような入力x」に加えて、単精度浮動小数点数を使ったがために厳密な場合分けができずに漏れてしまった「exp(x)がかなり大きな正規化数になるような入力x」「exp(x)がかなり小さな正規化数になるような入力x」「exp(x)が最も小さな非正規化数よりもやや小さくて0.0になるような入力x」の場合は、指数部のビット演算を工夫した後、浮動小数点数演算で最終結果を算出していました。これらのパターンが使われることはほとんどないですが、演算の量はそれなりに増えてしまうようです。

*1:実装アルゴリズムによっては、途中の計算結果がオーバーフローするが答えは正規化数になる、といったケースの例外処理が必要なこともあり得ます。

*2:論理レジスタ番号が2nと2n+1のレジスタという意味です。

*3:このような比較はunordered比較と言って一命令で実現可能です。

*4:xが-Infだとx+(+Inf)はNaNになりますが、まずxが負であれば0.0を返すというアルゴリズムなので、この計算結果が選ばれることはありません。

2020-11-15

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

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

実験環境

Intel Core i5 7200U（Turbo Boost時3.1GHz）

実験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程度のようです。

2020-11-01

量子コンピュータの勉強（その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)$ ゲート操作は別物です。

計算

上の量子ビットである $\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$ と置きます。

ここで、基底を計算基底（ $\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)ゲートをかけるとしていましたが、逆だと思います。

2020-10-25

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$ の絶対値が大きすぎる場合はどうせvscalefpdで0か+Infになるので、適当に小さくすればよいということのようです。

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

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

2020-10-18

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周の時に当てられない現象は、分岐命令のアドレスなどを変えて再検証する必要あり