浮動小数点数演算における保護桁について

浮動小数点数同士を足したり引いたりする場合、その数学的な結果は浮動小数点数で表せるとは限りません。そのため、結果を適当に浮動小数点数に戻す、丸めが行われます。丸めを行うためには、最終的な結果に使われる桁数より少し多く情報を保持しておく必要があります。そのような部分のことを、保護桁と言います。

正確に丸める最も簡単な方法は、丸める前の結果を正確に求める、というものです。しかしそうするには、かなり大きなビット幅の加算器が必要になります。幸運なことに、これを回避する方法があります。最終的な結果に使われる桁数よりもずっと離れた桁の計算はしないというものです。具体的には、最終的な結果に使われる桁数に加えて3bit余分にあれば、IEEE754で規定された丸めを行うことができます。

この記事では、なぜ余分に3bit必要なのかを説明します。この3bitには名前がついていて、それぞれ上の桁から「ガードビット」「ラウンドビット」「スティッキービット」となっています。ガードビットは保護桁の最上位のビット、ラウンドビットはその次のビットを意味しています。スティッキービットは特殊で、その次以降が全て0ではない（1がいくつか立っている）ことを意味するビットです。

以下では、説明が冗長になることを防ぐため、以下の太字用語を定義して使います。

加算は、同符号の浮動小数点数同士の加算、または異符号の浮動小数点数間の減算を意味します。
減算とは、異符号の浮動小数点数同士の加算、または同符号の浮動小数点数間の減算を意味します。

結論から書くと、以下のようになります。

加算の結果を偶数丸めする場合、ラウンドビットは不要で、ガードビットとスティッキービットのみが必要
減算の結果を偶数丸めする場合、ガードビットとスティッキービットだけでなくラウンドビットも必要
- 加算と違ってラウンドビットが必要なのは、繰り下がりが発生してガードビット部分が最終結果に使われる場合があるため

加算する場合

浮動小数点数では符号＋指数部＋絶対値表現なので、以下では符号のことは考えずに絶対値部分のみを考えます。

加算する場合、以下のように桁合わせをしてから加算器に入力することになります。ここで、桁合わせしてはみ出た部分は、最終結果には直接的には使われません。したがって、この部分が保護桁であるということになります。なお、最上位桁に繰上りが発生した場合は、最終結果に使われない部分が1bit増えますが、保護桁が不足するのではなく余るだけなので、問題にはなりません。

桁合わせしてはみ出た部分については、片側のオペランドが0なので、特に計算せずガードビット（G）とスティッキービット（S）が求まります。

ガードビットが必要なのは、最近接丸めを実現するためです。 0.5より少ないか多いかを判断するために必要になります。スティッキービットが必要なのは、偶数丸めを実現するためです。 0.5ちょうどなのか0.5より多いのかを判断するために必要になります。つまり、四捨五入丸め（round nearest, ties to away）をサポートすればよい場合、スティッキービットは不要です。また、ゼロへの丸めをサポートすればよい場合、ガードビットもスティッキービットも不要です。他にも丸めモードは多数ありますが、ガードビットとスティッキービットさえあればいずれにも対応できます。

丸めモードごとの各ビットの要不要は以下のようになります。ここで、「ガードビット不要・スティッキービット必要」とは、はみ出した部分に1があるかを判定するスティッキービットのみが必要であるという意味です。

丸めモード	G	S
ゼロへの丸め	不要	不要
正の無限大への丸め	不要	必要
負の無限大への丸め	不要	必要
偶数丸め	必要	必要
奇数丸め	不要	必要
四捨五入丸め	必要	不要

偶数丸めを行う場合、ガードビットとスティッキービットの使い方は以下のようになります。

G	S	動作
0	0	そのまま
0	1	そのまま
1	0	結果の仮数部の最下位が1ならば1を加える
1	1	結果の仮数部に1を加える

減算する場合

減算する場合も、以下のように桁合わせをしてから加算器に入力することになります。ここで、桁合わせしてはみ出た部分が保護桁なのですが、注意するべき点があります。それは、最上位桁で繰り下がりが発生した場合、最終結果に使われないとは言い切れない点です。ここでは、以下の二通りに分類して考えてみます。

指数部が近く、多数桁の桁落ちが発生しうる場合

両オペランドの指数部の差が0か1である場合が相当します。この場合、桁落ちが発生しますが、結果は浮動小数点数で正確に表せます。したがって、桁合わせしてはみ出た部分（高々1bit）を含めて加算器に入力し、出力を桁合わせすれば計算完了です。

指数部が近くないので、多数桁の桁落ちが発生しない場合

両オペランドの指数部の差が2以上である場合が相当します。多数桁の桁落ちが発生しないとは言っても、最上位桁の繰り下がりが発生する場合はあります。その場合、最終結果の有効数字が1bit足りなくなるので、保護桁から捻出する必要が発生します。よって、加算の時と比べてもう1bit必要で、ガードビット（G）・ラウンドビット（R）・スティッキービット（R）の3bitが必要です。

足し算の時はあふれたビットがそのまま保護桁となりましたが、引き算の場合はあふれたビットも含めて引き算をして保護桁の値を決定します。スティッキービットを決めるためには完全な引き算は必要なく、関係するビットを全てORしたものを減算器に入力すれば十分です。つまり、3bit分だけ幅広だけの減算器が必要ということになります。

減算の場合、四捨五入すればいい場合でもスティッキービットまで必要です（五捨五超入なら不要です）。

偶数丸めを行う場合、保護桁の使い方は以下のようになります。

最上位桁の繰り下がり	G	R	S	動作
なし	0	x	x	そのまま
なし	1	0	0	結果の仮数部の最下位が1ならば1を加える
なし	1	0	1	結果の仮数部に1を加える
なし	1	1	x	結果の仮数部に1を加える
あり	-	0	x	そのまま（Gが仮数部の最下位ビットとなる）
あり	-	1	0	Gが1ならば結果の仮数部に1を加える
あり	-	1	1	結果の仮数部に1を加える

なお、加算の時と回路を共有することを考えると、桁合わせの際に1bitずらして配置することで以下を共通化できるので、そういった構成法がよいかもしれません。

桁合わせに必要な1bit右シフト回路
- 1bitずらさないと、1bit左シフト回路も必要になる
ラウンドビットが不要
丸める際の補正が0か+1しか出てこない

参考文献

パターソン＆ヘネシーコンピュータの構成と設計上（第6版）日経BP，2021年
ディジタル数値演算回路の実用設計第6章浮動小数点演算の作法（https://www.cqpub.co.jp/hanbai/books/36/36171/36171_6syo.pdf）

2022-07-17

qsortの比較関数の書き方について

C言語標準ライブラリに含まれるqsortは、配列のソートを行ってくれる関数です。第一引数にソート対象の配列、第二引数に配列の要素数、第三引数一要素のサイズ、第四引数に比較関数、のように渡すことで、ソートを行ってくれます。ここで、第四引数に渡す比較関数をどのように書くかについてまとめてみます。

比較関数の仕様

int comp( const void* lhs, const void* rhs )のような関数を書く必要があります。ここで、返り値は、

lhsの方が小さい（より正確には、「配列の先頭に近い方に配置してほしい」）場合は、負の値
lhsとrhsが等しい場合は、0
lhsの方が大きい（より正確には、「配列の末尾に近い方に配置してほしい」）場合は、正の値

を返す必要があります。

C++におけるstd::sortの比較関数のように二値を返すのではなく、三値（以上）を返す必要がある点に注意します。

比較対象が`char`の場合

比較対象がcharの場合、うまい方法があります。例えば昇順にソートしたい場合、以下のように書けば、仕様を満たせます。

int comp( const void* lhs, const void* rhs ) {
    return *(char*)lhs - *(char*)rhs;
}

ここで、*(char*)lhsや*(char*)rhsはchar型ですが、算術演算のオペランドになっているため、int型に変換されてから演算が行われます。したがって、この計算でオーバーフローは発生せず、仕様通り「lhsの方が小さい時は負の値、等しければ0、lhsの方が大きい場合は正の値」を返すことになります。

比較対象が`int`の場合

比較対象がintの場合、先のような引き算を使った方法はうまくいきません。なぜなら、int同士の引き算の結果は、intに収まる保証がないためです。 long longなど、int型より大きい型を使えばよいと思うかもしれませんが、返り値をintに収める必要があるため、やはりうまくいきません。

※あまりにも絶対値の大きな値が来ないことが分かっている場合など、オーバーフローが発生しないことが確実な場合、引き算を使った方法を使っても安全です。

このような場合、以下のような分岐を用いた比較関数を書くことが紹介されることがあります。

int comp( const void* lhs, const void* rhs ) {
    int L = *(int*)lhs;
    int R = *(int*)rhs;
    if( L < R ) { return -1; }
    if( L > R ) { return 1; }
    return 0;
}

しかし、これよりも以下のように書いたほうが高速です。

int comp( const void* lhs, const void* rhs ) {
    int L = *(int*)lhs;
    int R = *(int*)rhs;
    int lt = L < R;
    int gt = L > R;
    return gt - lt;
}

この書き方が高速な理由は、分岐がなくなっているからです。ソート関数中の分岐を先取りしているだけなので分岐予測ミスの回数自体は変わらなさそうに思いますが、実際には後者の方が分岐予測ミス数が少なくなります。これはソート関数中の分岐は非常に予測困難であることが関係しているのかもしれません。

測定

分岐を使って比較結果を返すのではなく、return (L > R) - (L < R);と書いたほうが速いのかを実際に確かめてみます。測定環境は以下です。

Ubuntu 20.04.4 LTS (WSL2 on Windows 11)
gcc 9.4.0、最適化オプションは-O3 -march=native
Intel Core i9 12900K (4880MHzくらいで動いている）

実験に用いたソースコードは以下です。

#include <stdio.h>
#include <stdlib.h>
#include <time.h>

int comp( const void* lhs, const void* rhs ) {
  long L = *(long*)lhs;
  long R = *(long*)rhs;
  int gt = L > R;
  int lt = L < R;
  return gt - lt;
}

#define SIZE 100000000
long arr[SIZE];

int main() {
        for( int i = 0; i < SIZE; ++i ) {
                arr[i] = rand() + rand() * 0x100000000;
        }
        clock_t start = clock();

        qsort( arr, SIZE, sizeof(long), &comp );

        clock_t end = clock();

        printf( "Elapsed: %f\n", (double)(end - start) / CLOCKS_PER_SEC );
}

int comp( const void* lhs, const void* rhs ) {
  long L = *(long*)lhs;
  long R = *(long*)rhs;
  if (L > R) return 1;
  if (L < R) return -1;
  return 0;
}

測定の結果、前者は10.9±0.1秒程度でソートが終わりましたが、後者は11.3±0.1秒程度の時間がかかりました。このように、qsortに渡す比較関数中での分岐を減らすことは高速化につながることが分かります。

まとめ

charを比較するときには引き算を使うテクニックがある
int以上を比較したい場合には引き算のテクニックは使えない（オーバーフローするため）
オーバーフローを考慮に入れないといけない場合、(L > R) - (L < R)のように記述して分岐を減らす方が高速

2022-07-03

RISC-VのB拡張とB拡張対応gccのビルド方法

gcc12から、RISC-VのB拡張命令を使ったコンパイルが可能になりました。 B拡張には、ビット演算を取り扱う命令が含まれています。したがって、B拡張に対応したCPUでは、gcc12でコンパイルすることによりプログラムが高速化することが期待できます。

RISC-VのB拡張は、複数の拡張（Zba, Zbb, Zbc, Zbs, Zbp, Zbm, Zbf, Zbe, Zbr, Zbt）に分かれています。このうち、Zba, Zbb, Zbc, Zbs拡張のみがRISC-Vの公式に批准（Ratified）されており、残りの拡張はまだ仕様策定中です。 gcc12で対応しているのは、この批准されているZba, Zbb, Zbc, Zbs拡張のみです*1。

以下の図は、RISC-VのB拡張に含まれる命令が、どのZb*拡張に含まれているかを示したものです。例えば、sh1add命令はZba拡張に含まれており、andn命令はZbb拡張とZbp拡張に含まれています。 pack命令とその亜種は複数の拡張に含まれています。矢印の先にある命令（zext.h、rev8、orc.b、unzip8、unzip16）は、それぞれ矢印の元にある命令（pack[w]、grevi、gorci、unshfli、shfli）の特殊パターンであり、その命令ビットパターンは将来Zbp拡張を批准したときに問題とならないよう、互換性があるように定義されています。 zext.w命令はpack命令の特別なパターンとして定義される予定でしたが、批准されている拡張に含まれるadd.uw命令の特別なパターンの動作としても定義できるため、現在はそのように定義されているようです。

さて、B拡張に対応したgcc12をビルドする方法を以下に説明します。

まず、riscv-gnu-toolchainリポジトリをクローンします。ここで、--recursiveオプション付きでクローンすると以下の手順で失敗する*2ので気を付けてください。 riscv-gnu-toolchainリポジトリはsubmoduleをいくつか持っていますが、これらはビルド時に自動で取得されるため、--recursiveは不要です。クローン後、gcc-12ブランチに移動します。これは、現在のデフォルトブランチではgcc11がビルドされるためです。

git clone https://github.com/riscv/riscv-gnu-toolchain
cd riscv-gnu-toolchain/
git checkout gcc-12

./configureコマンドを実行することで、ビルド方法の設定を行います。

--prefixには、出来上がったgcc等をインストールする場所を指定します。

もし、RV64GB対応のgcc12を作りたい場合、--with-archオプションや--with-abiオプションは何も指定する必要がありません。 RV32GB対応のgcc12を作りたい場合、--with-arch=rv32g_zba_zbb_zbc_zbs --with-abi=ilp32を指定してください。

./configure --prefix ~/riscv64gb_gcc121
make linux -j$(nproc)

./configure --prefix ~/riscv32gb_gcc121 --with-arch=rv32g_zba_zbb_zbc_zbs --with-abi=ilp32
make newlib -j$(nproc)

あとは待てば、自動的にできあがったgccがインストールされます。

コンパイル時に-march=rv64g_zba_zbb_zbc_zbsなどとオプションをつけることで、B拡張が使われたバイナリが生成されます。

*1:clang12は、まだ批准されていないZbt拡張に含まれる命令も使うようです。

*2:asコマンドのビルドに失敗し、結果としてgcc等が使えない状態で終了します。

2022-05-08

Zen2の方向分岐予測器を調べてみる（その3）

前回に引き続き、Zen2の方向分岐予測器の仕組みを明らかにしていきます。

分岐履歴の仕組み

gselect(GAs)のような単純な分岐予測器は、分岐履歴をそのまま分岐予測テーブルのインデクスの一部として利用します。 gshareのような分岐予測器であっても、分岐履歴をそのままPCとXORして分岐予測テーブルのインデクスとして利用するのが一般的です。このようなことができるのは、（分岐履歴のビット数）≦（PHTのインデクスのビット幅）が成り立つときに限ります。

一方、Zen2の方向分岐予測器の分岐履歴は、一分岐命令当たり5bit×履歴長91で合計455bitもあります。したがって、分岐予測テーブルのインデクスに分岐履歴をそのまますべて使うことは不可能です。そのため、分岐履歴同士でXORするなどして情報を圧縮することが行われているはずです。この情報圧縮に使われる関数を特定してみます。今回はまず、分岐履歴の更新方法を調べてみます。

前回までは二つの分岐命令を使って実験していましたが、今回からは三つの分岐命令を使って実験していきます。実験に使うアセンブリコードは以下の通りです。

        .text
        .intel_syntax noprefix
        .globl main
main:
        mov rcx, 200000
.p2align 6, 0x90
L0:
        mov rax, （パターン長）
        （nopを何個か）
L1:
        sub rax, 1
        jne L1 # 分岐A
        （nopを何個か）
        je L2 # 分岐B
L2:
        （nopを何個か）
        sub rcx, 1
        jne L0 # 分岐C
        ret

以下では、nopの数を具体的に示すのではなく、分岐A, B, Cの最終バイトのアドレスをα, β, γとして記述します。

実験1：情報の追い出され方

β=0x401000, γ=0x402000に固定し、パターン長およびαを変化させました。その結果、分岐予測が200000回以上外れた組み合わせは以下のようになりました。

	82	83	84	85	86	87	88	89	90	91	92	93	備考
0x400001	×	×	×	×	×	×	×	×	×	×	×	×	Addr[0]は使われない
0x400002						×	×					×	Hash_0(Addr)=Addr[1]
0x400004					×	×					×	×	Hash_1(Addr)=Addr[2]^...
0x400008				×	×					×	×	×	Hash_2(Addr)=Addr[3]^...
0x400010					×	×					×	×	Hash_1(Addr)=Addr[4]^...
0x400020				×	×					×	×	×	Hash_2(Addr)=Addr[5]^...
0x400040			×	×					×	×	×	×	Hash_3(Addr)=Addr[6]^...
0x400080				×	×					×	×	×	Hash_2(Addr)=Addr[7]^...
0x400100			×	×					×	×	×	×	Hash_3(Addr)=Addr[8]^...
0x400200		×	×					×	×	×	×	×	Hash_4(Addr)=Addr[9]^...
0x400400			×	×					×	×	×	×	Hash_3(Addr)=Addr[10]^...
0x400800		×	×					×	×	×	×	×	Hash_4(Addr)=Addr[11]^...

この実験結果から、以下のことが分かります。

分岐履歴長を超えて古い情報が一度に無くなるのではなく、だんだんなくなっていく
- 普通のローリングハッシュではなく、左シフトで追い出されていくような感じ
なくなっていく順は、Hash_4、Hash_3、Hash_2、Hash_1、Hash_0の順
分岐履歴情報が完全になくなるまでの期間という意味において、分岐履歴長はやはり91

ところで、パターン長83～88のところでは、追い出されることとは無関係の分岐予測ミスが発生しています。以下の実験ではまず、このような現象の発生しないパターン長82で実験を行っていきます。

実験2：情報の重ね合わせ方

今度はパターン長を82、α=0x401000に固定し、βとγを変化させていきます。その結果、分岐予測が200000回以上外れた組み合わせは以下のようになりました。

これらのパターン（βγαα……α）で分岐予測に失敗するのは、αが連続しているパターン（αααα……α）と区別がつかないことが原因でしょう。つまりこれらのパターンでは、分岐履歴が「αの生成する分岐履歴情報＝0b00000が連続している場合の分岐履歴」と同一になっているはずです。

この結果から最低限分かることをまとめてみます。

分岐履歴中に1が立っているビットが1つ
- Hash_2(γ)のみが1の場合、区別がつかない
- Hash_1(γ)のみが1の場合、区別がつかない
- Hash_0(γ)のみが1の場合、区別がつかない
- Hash_1(β)のみが1の場合、区別がつかない
- Hash_0(β)のみが1の場合、区別がつかない
分岐履歴中に1が立っているビットが2つ
- Hash_3(β)とHash_4(γ)のみが1の場合、区別がつかない
- Hash_3(β)とHash_0(γ)のみが1の場合、区別がつかない
- Hash_2(β)とHash_3(γ)のみが1の場合、区別がつかない
- Hash_2(β)とHash_1(γ)のみが1の場合、区別がつかない
- Hash_1(β)とHash_2(γ)のみが1の場合、区別がつかない
- Hash_1(β)とHash_0(γ)のみが1の場合、区別がつかない
- Hash_0(β)とHash_3(γ)のみが1の場合、区別がつかない
- Hash_0(β)とHash_1(γ)のみが1の場合、区別がつかない

これらを観察すると、以下の非常に規則的な組み合わせでXORされていることが分かります。

Hash_3(β)とHash_4(γ)
Hash_2(β)とHash_3(γ)
Hash_1(β)とHash_2(γ)
Hash_0(β)とHash_1(γ)

その他の不規則な区別がつかない組み合わせは、分岐予測テーブルのインデクス関数由来だと考えられます。

また、パターン長75～81で実験しても同様の規則的な区別のつかない組み合わせが存在しました。この法則性は、分岐履歴を更新する際の操作が「1ビット左シフトしてから下位ビットに分岐履歴情報をXOR」という仕組みになっているのが由来だと考えられます。パターン長75～81で実験した際の不規則な区別のつかない組み合わせの出現パターンは、これと整合します。

wire Hash_0 = BrAddress[1];
wire Hash_1 = BrAddress[4] ^ BrAddress[2];
wire Hash_2 = BrAddress[7] ^ BrAddress[5] ^ BrAddress[3];
wire Hash_3 = BrAddress[10] ^ BrAddress[8] ^ BrAddress[6];
wire Hash_4 = BrAddress[11] ^ BrAddress[9];
wire [4:0] TargetHash = { Hash_4, Hash_3, Hash_2, Hash_1, Hash_0 };

TargetHistory <= (TargetHistory << 1) ^ TargetHash;

まとめ

分岐履歴は91bitのシフトレジスタになっている
古い分岐履歴情報は左シフトで追い出される（ローリングハッシュ的な追い出し方ではない）
分岐履歴の更新方法は、「1ビット左シフトしてから下位ビットに分岐履歴情報をXOR」になっている

2022-05-01

Zen2の方向分岐予測器を調べてみる（その2）

前回に引き続き、Zen2の方向分岐予測器の仕組みを明らかにしていきます。

分岐履歴に使われる情報

分岐元アドレス

前回と同様のアセンブリを用いて、分岐命令のアドレスのみを変化させていき、パターン長1～48の予測を当てることができるかを調べていきます。分岐命令のアドレスを変えるためには、NOPを詰める位置（分岐命令の上に置くか下に置くか）を変化させました。

その結果、以下のことが分かります。

分岐命令のアドレスが2n+1の時と2n+2の場合は正答パターンがほぼ一致
2n+1と2n+2の間で繰上りが発生することもある
- 上位ビットが同じで下位ビットが01と10の組でパターンが一致するということではなさそう

このことから、以下の仮説を立てます。

分岐命令のアドレス＋1が分岐履歴に使われる。この時、アドレスの下位1ビットは使われない。

分岐先アドレス

前回の調査の結果、パターン長が9以上で常に予測失敗する命令配置があることを知っています。このようなことが発生するのは、jne L1の生成する分岐履歴情報とjne L0の生成する分岐履歴情報が完全に一致しているからとみて間違いないでしょう。

これを利用すると、どのような分岐履歴情報を生成しているかを調べることができそうです。具体的には、jne L0の生成する分岐履歴情報を固定し、jne L1の生成する分岐履歴情報がそれと一致する場合を探索することで、法則性を見つけていきます。

以下では、手を抜いてパターン長が9の時だけ調べます。以下のアセンブリを実行したとき、分岐予測が十分多く外れるNOPを詰める量を調べてみます。

        .text
        .intel_syntax noprefix
        .globl main
main:
        mov rcx, 20000000
        mov rax, 1
.p2align 6, 0x90
L0:
        （nopを何個か）
        sub rax, 1
        jne L1
        mov rax, 9
        （nopを何個か）
L1:
        （nopを何個か）
        sub rcx, 1
        jne L0
        ret

すると、NOPを詰める量が例えば以下の値の時、分岐予測が十分多く外れることが分かります。

上	中	下	`jne L1`分岐元	`jne L1`分岐先	`jne L0`分岐元	`jne L0`分岐先
3	9	2	0x400547	0x400559	0x40055f	0x400540
3	9	3	0x400547	0x400559	0x400560	0x400540
3	8	4	0x400547	0x400558	0x400560	0x400540
3	7	5	0x400547	0x400557	0x400560	0x400540
3	6	6	0x400547	0x400556	0x400560	0x400540
3	5	7	0x400547	0x400555	0x400560	0x400540
3	4	8	0x400547	0x400554	0x400560	0x400540
3	3	9	0x400547	0x400553	0x400560	0x400540

この結果から、以下の結論を得ます。

分岐先アドレスの情報は分岐履歴に使われない

方向分岐しか存在しない場合、実行経路を一意に特定するためには分岐先アドレスは不要（間接分岐がある場合は分岐先アドレスも必要）であり、このような分岐履歴の持ち方になっていることは理解できます。

分岐元アドレスのハッシュ関数の特定

パターン長が9以上で常に予測失敗する命令配置をまとめてみました。

上	中	下	`jne L1`分岐元	`jne L0`分岐元	(A&~1)^(B&~1)	備考
1	0	0	0x400545	0x400552	0x14	前回発見
3	7	0	0x400547	0x40055b	0x14
3	8	0	0x400547	0x40055c	0x14
4	6	0	0x400548	0x40055b	0x14
4	7	0	0x400548	0x40055c	0x14
5	7	0	0x400549	0x40055d	0x14
5	8	0	0x400549	0x40055e	0x14
7	0	0	0x40054b	0x400558	0x14	前回発見
9	0	0	0x40054d	0x40055a	0x14	前回発見
3	9	2	0x400547	0x40055f	0x28
3	9	3	0x400547	0x400560	0x28
3	0	7	0x400547	0x40055b	0x14
3	0	11	0x400547	0x40055f	0x28
3	0	31	0x400547	0x400573	0x3c
3	0	123	0x400547	0x4005cf	0x9c	`jne L0`は6Byte命令
3	0	143	0x400547	0x4005e3	0xa0	`jne L0`は6Byte命令
3	0	163	0x400547	0x4005f7	0xb4	`jne L0`は6Byte命令
58	0	58	0x40057e	0x4005c5	0xb4	`jne L0`は6Byte命令
58	0	78	0x40057e	0x4005d9	0xa0	`jne L0`は6Byte命令
58	0	82	0x40057e	0x4005dd	0x9c	`jne L0`は6Byte命令
58	0	102	0x40057e	0x4005f1	0x88	`jne L0`は6Byte命令
60	0	6	0x400580	0x400593	0x14	`jne L0`は6Byte命令
60	0	26	0x400580	0x4005a7	0x28	`jne L0`は6Byte命令
60	0	46	0x400580	0x4005bb	0x3c	`jne L0`は6Byte命令
60	0	1142	0x400580	0x400a03	0xf88	`jne L0`は6Byte命令
60	0	1162	0x400580	0x400a17	0xf9c	`jne L0`は6Byte命令
60	1142	0	0x400580	0x400a07	0xf88	共に6Byte命令
60	1154	0	0x400580	0x400a13	0xf9c	共に6Byte命令
60	2542	0	0x400580	0x400f7f	0xa00	共に6Byte命令
60	2798	0	0x400580	0x40107f	0x1500	共に6Byte命令
60	3118	0	0x400580	0x4011bf	0x1440	共に6Byte命令
60	3886	0	0x400580	0x4014bf	0x1140	共に6Byte命令
60	3958	0	0x400580	0x401507	0x1088	共に6Byte命令
60	3982	0	0x400580	0x40151f	0x10a0	共に6Byte命令
60	4078	0	0x400580	0x40157f	0x1000	共に6Byte命令
60	4090	0	0x400580	0x40158b	0x1014	共に6Byte命令
60	4118	0	0x400580	0x4015a7	0x1028	共に6Byte命令
60	8174	0	0x400580	0x40257f	0x2000	共に6Byte命令
60	16366	0	0x400580	0x40457f	0x4000	共に6Byte命令
60	32750	0	0x400580	0x40857f	0x8000	共に6Byte命令

まず「分岐命令のアドレス＋1が分岐履歴に使われる」の謎が解明できました。先ほどまでの実験に使っている分岐命令は全て2Byteのものだったので＋1になっていましたが、正確には「分岐命令の最後のバイトのアドレスが使われる」ということのようです。

また、アドレスのXORで見た差分が0x14の時に分岐履歴情報が区別できないことから、アドレスの2bit目（Addr[2]）とアドレスの4bit目（Addr[4]）はXORされた形でのみ分岐履歴上に存在することが分かります。同様に、以下のような組み合わせでXORされていることが分かります。

(A&~1)^(B&~1)
0x14	Addr[2]	Addr[4]
0x28	Addr[3]	Addr[5]
0x88	Addr[3]	Addr[7]
0xa0	Addr[5]	Addr[7]
0x140	Addr[6]	Addr[8]
0x440	Addr[6]	Addr[10]
0x500	Addr[8]	Addr[10]
0xa00	Addr[9]	Addr[11]

したがって、ハッシュ関数は以下のようになっており、分岐元アドレスの情報を5bitに圧縮していることが分かります。

Hash_0(Addr) = Addr[1]
Hash_1(Addr) = Addr[2]^Addr[4]
Hash_2(Addr) = Addr[3]^Addr[5]^Addr[7]
Hash_3(Addr) = Addr[6]^Addr[8]^Addr[10]
Hash_4(Addr) = Addr[9]^Addr[11]

まとめ

分岐履歴の構築には、分岐命令の最後のバイトのアドレスが使用される
分岐履歴の構築には、分岐先のアドレスは使用されない
分岐命令の最後のバイトのアドレスの最下位ビットは使われない
- 全ての分岐命令は2Byte以上なので、最下位ビットのみ異なる分岐命令の組は通常存在しない
分岐命令の最後のバイトのアドレスの2の位～2048の位（合計11ビット）がハッシュ関数に通されて5ビットに圧縮されて使用される

2022-04-24

Zen2の方向分岐予測器を調べてみる（その1）

最初はGolden Coveの方向分岐予測器を調べていたのですが、規則的な分岐履歴系列を与えても正答率が100%に張り付かないなど、あまりにもアナログで複雑な挙動を示すのであきらめました。それと比べてZen2の方向分岐予測器は、正答率が0%や100%に張り付く非常にデジタルな動作をするようで、中身を推測しやすそうです。

測定環境

AMD EPYC 7702 64-Core Processor
gcc (GCC) 7.3.1 20180303 (Red Hat 7.3.1-5)
- スタティックリンクを用いた
perf version 3.10.0-1160.15.2.el7.x86_64.debug

使用される分岐履歴の長さの推定

単純なループの例

まずは、以下のアセンブリを実行した際の分岐予測ミス数を調べてみます。

        .text
        .intel_syntax noprefix
        .globl main
main:
        mov rcx, 20000000
        mov rax, 1
.p2align 6, 0x90
L0:
        sub rax, 1
        jne L1
        mov rax, （パターン長）
L1:
        sub rcx, 1
        jne L0
        ret

すると、パターン長が33回、39回、44回、45回、47回以上、の時に分岐予測ミスが（20000000÷パターン長）回程度発生します。この分岐予測ミスは、jne L1部分に由来するとみて間違いないでしょう。

しかし、ここから分岐履歴長は90程度である、と結論付けるのは危険です。もう少し調査してみましょう。

分岐命令のアドレスを変更してみると

以下のように、NOPをいくつか挿入することで、分岐命令のアドレスを変化させてみます。

        .text
        .intel_syntax noprefix
        .globl main
main:
        mov rcx, 20000000
        mov rax, 1
.p2align 6, 0x90
L0:
        （nopを何個か）
        sub rax, 1
        jne L1
        mov rax, （パターン長）
L1:
        sub rcx, 1
        jne L0
        ret

すると、分岐予測ミス数は以下の図のようになりました。 Excelの機能で、十分分岐予測ミス数が多かった部分を色付けしました。

この結果から、以下のことが言えます。

パターン長が47以上だと当てられない
分岐命令のアドレスによっては、パターン長が47未満でも予測を外すことがある
- 予測を外すパターン長はNOPの数により何グループかに分けられ、4NOPや8NOPを中心に対称
  - 分岐予測テーブルにアクセスする際のインデクス関数がXOR等の比較的単純なものであることを示唆

頑張れば分岐予測テーブルのアクセスに使うインデクス関数を予測出来そうです。しかしそのためにはまず、どういった種類の分岐履歴が使われているかを明らかにする必要があります。

分岐履歴の種類

分岐履歴に種類なんてあるのか、という感じですが、よく使われる分岐履歴には以下の二種類があります。

分岐方向履歴

一般的な教科書に書かれている、もっとも有名なものです。

分岐方向履歴はその名の通り、分岐方向（成立or不成立）の情報を一列に並べたものです。

成立（taken）の場合は1を、不成立（not taken）の場合は0を、ビット列の末尾に追加する、といった方法で作られます。

分岐先履歴

商用プロセッサでよく使われている[1]方法だそうです。

分岐先履歴は、分岐が成立したときのみ更新される履歴です。更新時には、その分岐命令のアドレスを分岐先のアドレスを組み合わせてハッシュ関数に通した値が使われます。

分岐先履歴の利点は以下のようにたくさんあり、商用プロセッサで使われるのも理解できます。

複数の分岐命令を並列に予測可能
- 分岐方向履歴だと、「直前の命令は分岐命令でない」「直前の命令は分岐命令で条件不成立」の二通りの可能性を考慮する必要がある
不成立の分岐を見逃しても履歴が壊れない
- 分岐方向履歴だと、「この命令は条件分岐命令か」も正しく予測する必要がある
（ハッシュ関数を通さずフルの情報が使える場合）実行経路情報を一意に特定できる
- 分岐方向履歴は、実行経路が異なっても同じになることがある

使用されている分岐履歴の種類と履歴長の推定

以下のアセンブリを実行した際の分岐予測ミス数を調べてみます。

        .text
        .intel_syntax noprefix
        .globl main
main:
        mov rcx, 20000000
        mov rax, 1
.p2align 6, 0x90
L0:
        （nopを何個か）
        sub rax, 1
        je L1
        sub rax, （パターン長）
L1:
        add rax, （パターン長）
        sub rcx, 1
        jne L0
        ret

すると前回とは違い、NOPの数によらず予測できなくなるパターン長は92以上の場合となります。

この違いは、分岐パターンの違いにあります。

最初に使ったアセンブリでは、分岐パターンは TTTTTTTT.....TTNT を繰り返すものでした（ここで、Tは条件成立（taken）、Nは条件不成立（not taken）を意味します）。一方、今回使ったアセンブリでは、分岐パターンはNTNTNTNT...NTTTを繰り返すものになっています。

つまり、分岐パターンにNが多い時は長いパターンを当てることができるようになっています。

したがって、Zen2の分岐予測器では、分岐先履歴が使われているといってよいでしょう。

また、その履歴長はおそらく91でしょう。その根拠は以下です。

最初に使ったアセンブリで、履歴長が92あればパターン長が47の時を当てられるはずだが当てられていない
今回使ったアセンブリで、履歴長が92あればパターン長が92の時を当てられるはずだが当てられていない

まとめ

Zen2の方向分岐予測器の仕組みを解析してみました。その結果、分岐先履歴が使われていること、履歴長は91であること、分岐予測テーブルにアクセスする際のインデクス関数は単純そうであること、が分かりました。

参考文献

[1] Y. Issii, et al.,"Rebasing Instruction Prefetching: An Industry Perspective", IEEE Computer Architecture Letters, 19(2), 2020.

2022-04-03

IntelのGracemontコアの浮動小数点数演算レイテンシの測定

第12世代インテル® Core™ プロセッサーのEコア（高効率コア）はGracemontというマイクロアーキテクチャでできています。このコアで命令を実行したときのレイテンシを調べたものにuops.info - Tableがありますが、測定方法がよくなくて正しい値になっていません。そこでnanoBenchを用いて独自に測定してみました。

測定方法

sudo taskset -c 23 ./nanoBench.sh -asm "測定コード" -config configs/cfg_AlderLakeE_common.txt

測定結果

`vaddsd` および `vshufpd`

vaddsd xmm0, xmm0, xmm1 → 3.00 cycle
vaddsd xmm0, xmm0, xmm1; vaddsd xmm0, xmm0, xmm1 → 6.00 cycle
vshufpd xmm0, xmm0, xmm1, 0 → 1.00 cycle
vshufpd xmm0, xmm0, xmm1, 0; vshufpd xmm0, xmm0, xmm1, 0 → 2.00 cycle
vaddsd xmm0, xmm0, xmm1; vshufpd xmm0, xmm0, xmm1, 0 → 7.00 cycle
vaddsd xmm0, xmm0, xmm1; vshufpd xmm0, xmm0, xmm1, 0; vaddsd xmm0, xmm0, xmm1; vshufpd xmm0, xmm0, xmm1, 0 → 14.00 cycle
vaddsd xmm0, xmm0, xmm1; vaddsd xmm0, xmm0, xmm1; vshufpd xmm0, xmm0, xmm1, 0; vshufpd xmm0, xmm0, xmm1, 0 → 11.00 cycle

以上の結果から、以下のことが分かります。

vaddsd命令のレイテンシは 3 cycle
vshufpd命令のレイテンシは 1 cycle
vaddsd命令とvshufpd命令を行き来すると、往復で 3 cycle のペナルティ
- uops.info の測定はこのペナルティを考慮できていないため間違い

`vmulsd` および `vfmadd132sd`

vmulsd xmm0, xmm1, xmm1 → 4.00 cycle
vmulsd xmm0, xmm0, xmm1; vmulsd xmm0, xmm0, xmm1 → 8.00 cycle
vmulsd xmm0, xmm0, xmm1; vshufpd xmm0, xmm0, xmm1, 0 → 7.00 cycle
vfmadd132sd xmm0, xmm1, xmm2 → 6.00 cycle
vfmadd132sd xmm0, xmm1, xmm2; vshufpd xmm0, xmm0, xmm1, 0 → 10.00 cycle
vmulsd xmm0, xmm0, xmm1; vaddsd xmm0, xmm0, xmm1 → 6.08 cycle
vmulsd xmm0, xmm0, xmm1; vaddsd xmm0, xmm0, xmm1; vaddsd xmm0, xmm0, xmm1 → 9.08 cycle
vmulsd xmm0, xmm0, xmm1; vaddsd xmm0, xmm0, xmm1; vmulsd xmm0, xmm0, xmm1 → 10.08 cycle
vmulsd xmm0, xmm0, xmm1; vaddsd xmm0, xmm0, xmm1; vmulsd xmm0, xmm0, xmm1; vaddsd xmm0, xmm0, xmm1 → 12.16 cycle
vmulsd xmm0, xmm0, xmm1; vmulsd xmm0, xmm0, xmm1; vaddsd xmm0, xmm0, xmm1; vaddsd xmm0, xmm0, xmm1 → 13.08 cycle
vfmadd132sd xmm0, xmm1, xmm2; vmulsd xmm0, xmm0, xmm1; vaddsd xmm0, xmm0, xmm1 → 12.08 cycle
vfmadd132sd xmm0, xmm1, xmm2; vaddsd xmm0, xmm0, xmm1; vmulsd xmm0, xmm0, xmm1 → 13.00 cycle

以上の結果から、以下のことが分かります。

vmulsd命令のレイテンシは 4 cycle
vfmadd132sd命令のレイテンシは 6 cycle
vmulsd命令とvshufpd命令と行き来すると往復で 2 cycle のペナルティ
vfmadd132sd命令とvshufpd命令と行き来すると往復で 3 cycle のペナルティ
- uops.info の測定はこのペナルティを考慮できていないため間違い
vmulsd命令からvaddsd命令にフォワーディングするとき、1 cycle 得することがある（92%ほど発生）

浮動小数点数演算命令と`vshufpd`の行き来はどこにペナルティがあるか

movq rax, xmm0; movq xmm0, rax → 10.00 cycle
vshufpd xmm0, xmm0, xmm1, 0; movq rax, xmm0; movq xmm0, rax→ 11.00 cycle
vaddsd xmm0, xmm0, xmm1; movq rax, xmm0; movq xmm0, rax → 13.00 cycle

以上の結果から、以下のことが分かります。

movq rax, xmm0; movq xmm0, rax命令列で移動させると合計で 10 cycle かかる
この命令列を使った場合、vshufpd命令やvaddsd命令と組み合わせてもペナルティは発生しない

これを使って以下のような実験を行うと、どこでペナルティが発生するかを特定することができます。

vaddsd xmm0, xmm0, xmm1; vshufpd xmm0, xmm0, xmm1, 0; movq rax, xmm0; movq xmm0, rax → 17.00 cycle
vshufpd xmm0, xmm0, xmm1, 0; vaddsd xmm0, xmm0, xmm1; movq rax, xmm0; movq xmm0, rax → 14.00 cycle

したがって、以下のことが分かります。

vaddsd命令からvshufpd命令にフォワーディングするとき、3 cycle のペナルティが発生する
vshufpd命令からvaddsd命令にフォワーディングするときは、ペナルティが発生しない

vmulsd命令、vfmadd132sd命令についても同じ実験を行うことで、どこでペナルティが発生するかを特定することができます。

vmulsd xmm0, xmm0, xmm1; vshufpd xmm0, xmm0, xmm1, 0; movq rax, xmm0; movq xmm0, rax → 17.00 cycle
vshufpd xmm0, xmm0, xmm1, 0; vmulsd xmm0, xmm0, xmm1; movq rax, xmm0; movq xmm0, rax → 15.00 cycle
vfmadd132sd xmm0, xmm1, xmm1; vshufpd xmm0, xmm0, xmm1, 0; movq rax, xmm0; movq xmm0, rax → 20.00 cycle
vshufpd xmm0, xmm0, xmm1, 0; vfmadd132sd xmm0, xmm1, xmm1; movq rax, xmm0; movq xmm0, rax → 17.00 cycle

したがって、以下のことが分かります。

vmulsd命令からvshufpd命令にフォワーディングするとき、2 cycle のペナルティが発生する
vfmadd132sd命令からvshufpd命令にフォワーディングするとき、3 cycle のペナルティが発生する

まとめ

vshufpd命令のレイテンシは 1 cycle
vaddsd命令のレイテンシは 3 cycle
vmulsd命令のレイテンシは 4 cycle
vfmadd132sd命令のレイテンシは 6 cycle
vmulsd命令からvaddsd命令にフォワーディングするとき、1 cycle 得することがある
vaddsd命令からvshufpd命令にフォワーディングするとき、3 cycle のペナルティが発生する
vmulsd命令からvshufpd命令にフォワーディングするとき、2 cycle のペナルティが発生する
vfmadd132sd命令からvshufpd命令にフォワーディングするとき、3 cycle のペナルティが発生する

加算する場合

減算する場合

指数部が近く、多数桁の桁落ちが発生しうる場合

指数部が近くないので、多数桁の桁落ちが発生しない場合

参考文献

比較関数の仕様

比較対象がcharの場合

比較対象がintの場合

測定

まとめ

分岐履歴の仕組み

実験1：情報の追い出され方

実験2：情報の重ね合わせ方

まとめ

分岐履歴に使われる情報

分岐元アドレス

分岐先アドレス

分岐元アドレスのハッシュ関数の特定

まとめ

測定環境

使用される分岐履歴の長さの推定

単純なループの例

分岐命令のアドレスを変更してみると

分岐履歴の種類

分岐方向履歴

分岐先履歴

使用されている分岐履歴の種類と履歴長の推定

まとめ

参考文献

測定方法

測定結果

vaddsd および vshufpd

vmulsd および vfmadd132sd

浮動小数点数演算命令とvshufpdの行き来はどこにペナルティがあるか

まとめ

比較対象が`char`の場合

比較対象が`int`の場合

`vaddsd` および `vshufpd`

`vmulsd` および `vfmadd132sd`

浮動小数点数演算命令と`vshufpd`の行き来はどこにペナルティがあるか