リオーダーバッファのサイズを測ってみた

リオーダーバッファは、アウトオブオーダー実行をするプロセッサにおいて、その投機的な状態を保存するバッファです。 リオーダーバッファのサイズは、命令ウィンドウ(実行順序を並べ替える範囲)の大きさの上限を決めます。

命令ウィンドウに含まれる、依存関係のない複数の命令列は並列に実行できますが、命令ウィンドウに含まれない場合は並列に実行できません。 そこで、依存関係のない複数の命令列の間の命令数を徐々に増やしていくと、並列に実行できる命令列から並列に実行できない命令列に変わっていきます。 これを利用して、リオーダーバッファのサイズを測ってみます。

使用するアセンブリプログラム

.text
.intel_syntax noprefix
.globl main
main:
.cfi_startproc
        mov eax, 10000000
.p2align 4, 0x90
L:
        xorpd xmm0, xmm0
        addsd xmm0, xmm0
        addsd xmm0, xmm0
        (略:addsdが30個)
        addsd xmm0, xmm0
        addsd xmm0, xmm0
        xor edx, edx
        xor edx, edx
        (略:xorがN個)
        xor edx, edx
        xor edx, edx
        dec eax
        jne L
.cfi_endproc

addsd命令を直列につなげることで、なかなか実行が終わらない命令を作り出します。 xor命令は、次のループに含まれるaddsd命令列との距離を離すために使用する、ただの無意味な命令です。

測定結果

HaswellアーキテクチャのCPUで測定した結果が以下の図になります。横軸がxor命令の数、縦軸が実行にかかったサイクル数です。

f:id:lpha_z:20210221231421p:plain

Haswellアーキテクチャのリオーダーバッファのサイズは192なので、xor命令が188個以上の場合、直列に並んだ最後のaddsd命令が終わるまで次のループのaddsd命令を命令ウィンドウに取り込めません。 測定結果は、これと整合します。

また、xor命令が78個以上の場合、直列に並んだ最後のaddsd命令が終わるまで次の次のループのaddsd命令を命令ウィンドウに取り込めません。 これも測定結果と整合しています。

さらに、xor命令が41個以上の場合、直列に並んだ最後のaddsd命令が終わるまで3つ次のループのaddsd命令を命令ウィンドウに取り込めません。 これも測定結果と整合しています。

Haswellアーキテクチャではaddsd命令は3並列でしか実行できないので、それ以上xor命令を減らしても実行にかかる時間はほぼ一定です。

なお、n=1, 2, 3以外の部分で少しずつ性能が下がるのは、無意味な命令がフェッチ幅を浪費していることによるものです。 実際、その部分での傾きは0.25であり、Haswellアーキテクチャのフェッチ幅4 instructions/cycleと整合します。

Apple M1 でもやってみた

同様のプログラムをARMでも作成し、Apple M1で動かしてみました。 長いレイテンシを作るための命令は、50個配置しました。 横軸は無意味な命令として使ったmov x9, #0命令の数、縦軸は実行にかかった秒数です。

f:id:lpha_z:20210221230906p:plain

無意味な命令の数が622個のところで傾きが変わっていることから、Apple M1(の高性能コア)のリオーダーバッファのサイズは626であるとわかります。

これは、他の方法(非常に長いレイテンシを持つメモリアクセス命令を用いた方法)により測定された結果(Apple's Humongous CPU Microarchitecture - Apple Announces The Apple Silicon M1: Ditching x86 - What to Expect, Based on A14)とも一致しています。

RISC-Vで関数ポインタ呼び出しにjalr t0を使ってはいけない

概要

RISC-Vでは、jalr t0という命令には特別な意味が割り当てられているので、関数ポインタを用いた関数呼び出しのために使うと一部のプロセッサで性能低下を引き起こします。 t0以外のレジスタを使う場合は問題なく動作するので、関数ポインタの格納にはt0以外のレジスタを使いましょう。

詳細

高性能プロセッサには、分岐命令を実行せずに分岐先を予測する、分岐予測といった仕組みが実装されています。 分岐命令のうち、関数からの復帰命令は分岐先が毎回同じとは限らないので、「前回の分岐先を記憶しておく」などの方法では予測ができないことがあります。 これに対しては、関数呼び出し命令の次の番地を記録しておくスタックを用意することで簡単に解決が可能です。 このような分岐予測に使うスタックをreturn address stack (RAS)と呼びます。

ところで、「これは関数呼び出し命令」「これは関数からの復帰命令」というのは高級言語からコンパイルするときの約束(呼び出し規約)です。 これを守らずに書かれた手書きアセンブリコードは、RASが搭載されたプロセッサでは正しく動作しないのでしょうか。

もちろん、そんな風にはなっていません。 RASが搭載されたプロセッサであっても、プログラムはRASが搭載されていないプロセッサと同じように動作します。 これは、実行ステージで必ず分岐予測が正しいかを確かめるためです。 したがって、分岐予測が間違っていたとしても、プログラムが誤動作することはありません。 とはいえ、分岐予測の失敗は一般に性能低下につながるため、なるべく避けたいものです。

分岐予測の失敗を防ぐためには、RASへのpush, popを正しく行う(RASを適切に動作させる)必要があります。 そのためには、pushすべき命令(関数呼び出し命令)やpopすべき命令(関数からの復帰命令)を正しく見分ける必要があります。 従来の命令セット設計では、専用の関数呼び出し命令や関数からの復帰命令を設けることが普通でした。 こうすることで、RASの適切な動作を自明に判断することができます。

これに対して、RISC-Vでは、専用の命令を設けることはせず、ソースレジスタ・デスティネーションレジスタとして何が指定されているかを根拠にRASの動作を決定します。 RASの動かし方の指南はRISC-V Unprivileged ISA V20191213のp. 22に書いてあります。

このマニュアルを見ると、関数ポインタ経由の関数呼び出しの時に使うjalr命令で、ソースレジスタx5、デスティネーションレジスタx1としたjalr x1, x5はpop, then pushと書かれています。 関数呼び出し命令ではpushのみすべきですから、関数呼び出しのつもりでjalr x1, x5を使うと分岐予測の失敗につながります。 なお、jalr x1, x5jalr t0と略記されることがあるので、本記事のタイトルということになります。

なぜRASの動かし方が変則的なのか

RISC-Vで、関数呼び出し専用の命令、などを設けていない理由は、命令数の削減にあります。 ソースレジスタやデスティネーションレジスタに何が指定されているかを根拠にRASの動作を決める、というのは実際可能です。

しかし、RISC-Vにalternate link registerが存在することを加味すると、話が一気にややこしくなります。 alternate link registerというのは「第二のリンクレジスタ」くらいの意味で、millicode呼び出しというコード圧縮のテクを実現するためのもので、それ自体は問題の原因ではありません。 問題の原因は、RISC-Vではこれを「二つのコルーチンを行き来するという状況は、リンクレジスタが二つあればリンクレジスタの退避が不要である」という思い付きを実現するために転用していることです。

確かに二つのコルーチンを行き来する状況はRASを適切に動作させる(popしてpushする)ことで分岐予測を成功させることができますが、無理やり感が否めません。 これのせいでjalr x1, x5はpop, then pushという動作をすると決められています。

落とし穴:アセンブリ言語で書かないから関係ないよ

コンパイラが間違ってjalr t0を出力することがあります。 少なくとも、LLVM9のRISC-Vコンパイラはこのコードを出力します。 gccはこのコードを避ける工夫が行われています。

LLVMRISC-Vコンパイラ(LLVM9で確認)は、引数が0個の関数ポインタ呼び出しの場合は関数ポインタの格納にa0を使います。 以下、引数が1個、2個、……、7個の場合、関数ポインタの格納にa1a2、……、a7を使います(要するに、使われていない最も若い引数用レジスタを使うという戦略のようです)。 しかし、引数が8個以上の場合、関数ポインタの格納にt0を使うため、jalr t0となって性能低下を引き起こすコードを出力します。

gccRISC-Vコンパイラ(gcc7で確認)の場合、引数が5個以下であればa5、6個であればa6、7個であればa7、8個以上であればt1を使います。 gccのコードを確認したところ、明示的にjalr t0を避ける工夫が行われていることがわかりました。 gcc/constraints.md at master · gcc-mirror/gcc · GitHub

一般的な場合を高速化?

二つのコルーチンを行き来するという限定的な状況でしか役に立たないこの機能は、「一般的な場合を高速化」という原則に反しているような気がします。 現状、この機能を活用するコンパイラはおそらく無いので、ハードウェアが複雑化するだけの仕様です。 というか誤爆するコンパイラがあるせいで、性能低下の原因にしかなっていません。 コンパイラの出すコードが正しいと思ってハードウェア側を直すプログラマもいそうです。

今後のコンパイラの進展に期待する場合、後方互換性のために事前に定義しておかないといけない、ということかもしれません。 あるいは、ハードウェアが提供している面白機能として見ると楽しいですが、具体的に役立つ例は思いつきません。

なお、三つ以上のコルーチンを行き来する状況では、この機能はほぼ役に立ちません。 対称コルーチンで三つ以上のコルーチンを巡回する場合は明らかに役に立ちません。 非対称コルーチンでMain→Croutine A→Main→Coroutine B→Mainのように行き来する場合も無駄です(Main側ではRAS pushする命令を、Coroutine側ではRAS popをする命令を使えば十分です)。

glibcのexp(x)をいじめる

glibcのexp関数は、倍精度浮動小数点数の指数関数をソフトウェアで計算するものです。 glibcの実装(IBM Accurate Mathematical Libraryの実装を用いたもの)は、おそらく完全精度を達成しています。 完全精度とは、入力がどんな値であっても真の値に最も近い浮動小数点数に丸めた結果を返すことを言います。

しかし、これを実現するのは非常に困難であることが知られています(テーブルメーカーのジレンマ)。 glibcの実装では、以下の三ステップを踏むことで、完全精度を達成しています。

  1. まず倍精度浮動小数点数演算とテーブル引きを駆使して近似値を求め、誤差を加味しても丸めた結果に影響がないと判断されれば、その値を返す
  2. 144bitの多倍長演算を使って近似値を求め、誤差を加味しても丸めた結果に影響がないと判断されれば、その値を返す
  3. 768bitの多倍長演算を使って近似値を求め、その値を返す

なお、倍精度指数関数については、triple-double(疑似六倍精度)で注意深く計算すれば、完全精度を達成できることが知られています(crlibmによる証明があります→https://hal-ens-lyon.archives-ouvertes.fr/ensl-01529804/file/crlibm.pdf)。

このような仕組みになっているため、多くの場合には高速であるものの、特定のワーストケース入力の場合には非常に遅くなります。

2.のステップに突入する数として、簡単な例に-0.1841があります。exp(-0.1841)を計算すると、普通の数に比べて90倍ほど低速です。

3.のステップに突入する数はそんなに簡単には見つかりません。 exp(x)の丸めがぎりぎりになる値を探索するという研究の結果(→https://hal.inria.fr/inria-00072594/file/RR2000-35.pdf)を参照すると、0x1.9e9cbbfd6080bp-31が最も丸めが難しいと記されていました。 これを入力してみると、普通の数に比べて1000倍ほど低速でした。

[Wandbox]三へ( へ՞ਊ ՞)へ ハッハッ

exp(0x1.9e9cbbfd6080bp-31)を正確に計算してみると、以下のようになり、丸める方向を正しく判断するために小数点以下112bit以上の計算が必要なことがわかります。

 1 0x1.0000000000000000000000000000000
 x 0x0.000000033D3977FAC10160000000000
 \frac 1 2 x^2 0x0.00000000000000053EFE9FFA55F3F7B...
 \frac 1 6 x^3 0x0.000000000000000000000005AA0EAD6...
 \exp(x) 0x1.000000033D397800000000000002A51...

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

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

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

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

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

  • 0x1.fffffffffffff8p+02.0に丸めた場合、丸め誤差は2^-53、2.0の最終桁の重みは2^-51なので、丸め誤差は0.25ULP(?)
  • 0x1.fffffffffffff8p+00x1.fffffffffffffp+0に丸めた場合、丸め誤差は2^-53、0x1.fffffffffffffの最終桁の重みは2^-52なので、丸め誤差は0.5ULP(?)

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

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

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

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

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

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

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

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

  • 真の値を挟む二つの浮動小数点数は、0x1.0000000000000p+00x1.0000000000001p+0なので、1.の定義では1ULP=2^-52
  • 真の値に最も近い二つの浮動小数点数は、0x1.fffffffffffffp-10x1.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でみた丸め誤差」という関数が不連続になって取り扱いが不便です。

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シミュレータは、その構造上、以下の本質的な問題が存在します。

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

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

  1. あるロード命令由来のアクセスが、L1Dキャッシュミスを起こした後L2キャッシュでヒットする。
  2. その次の命令が、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%も命令キャッシュミスが起こるとしています。

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を返すというアルゴリズムなので、この計算結果が選ばれることはありません。

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

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

実験環境

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