浮動小数点数の演算の丸めモードを制御する

以前、Haskellで浮動小数点数の丸めモードを制御したいという話を書いた。

Haskellでの浮動小数点数の方向付き丸めを考える

以前の記事では、丸めモードを変える方法としてはC言語のfesetroundだけを扱ったが、高速化するなら各アーキテクチャ固有の方法を使えた方が有利である。また、動的な丸めモード指定ではない静的な丸めモードについてもこの記事では触れる。

動的な丸めモード指定

浮動小数点数演算のできるプロセッサには大抵浮動小数点数演算の挙動を指定するレジスタがあり、そこで丸めモードを指定できることが多い^[要出典]。

C言語などの一部のプログラミング言語でも、そういうプロセッサの丸めモードを扱えるような機能が用意されている。

C99

<fenv.h> の関数と #pragma を併用することにより丸めモードを変更できる。

過去記事でも触れたので、詳しくは省略する。

実行環境によっては浮動小数点数の丸めモードの変更がサポートされていないこともありうる（WebAssemblyとかはいかにも非対応そうだ）。そういう場合、fesetroundは0以外の値を返す。

欠点として、環境によっては遅い場合がある。丸めモードを変更する関数 fesetround はCの関数として実装されている以上、コンパイラーがインライン化でもしてくれない限り関数呼び出しのコストが発生する。

なので、実行環境ごとに #ifdef する覚悟があるのであれば、C99の方法ではなく実行環境ごとの組み込み関数・インラインアセンブリを使った方が早い。

x86系

x86系には歴史的事情により、x87 FPUとSSE系（という呼び方でいいのか）の2つの浮動小数点数演算器がある。

x86系におけるfesetroundはx87 FPUとSSE系の両方の丸めモードを変更する。もしもx87 FPUかSSE系のいずれかしか使わないのであれば、fesetroundの代わりにここで述べるプラットフォーム依存な方法を使うことで丸めモードの切り替えを高速化できる。

x87 FPU

x87 FPUは16ビットのControl Wordと呼ばれる状態を持ち、そのうち2ビットが丸めモードを指定する。

（x87 FPUのデータレジスタは内部的に80ビット（精度64ビット）の幅があり、デフォルトではIEEE754のbinary64やbinary32準拠とならないことで悪名高い。Control WordのうちPrecision Control fieldと呼ばれるビットをいじればこの辺の挙動を変えて演算を精度53ビットや精度24ビットで行うことが可能となる。）

x87 FPUのControl Wordを変更するには、FSTCW (Store x87 FPU Control Word), FLDCW (Load x87 FPU Control Word)の命令を使う。「Store」はControl Wordの値をメモリにストアする、「Load」はメモリの値をControl Wordにロードする方である。

GCC系のインラインアセンブリを使って、丸めモードを指定した加算を行うコード例は次のようになる：

/* mode は 0: To Nearest, 1: Toward -∞, 2: Toward +∞, 3: Toward Zero */
long double rounded_add(unsigned int mode, long double x, long double y) {
    uint16_t oldcword;
    asm("fstcw %0" : "=m"(oldcword));
    uint16_t newcword = (oldcword & ~(3u << 10)) | ((uint16_t)mode << 10);
    asm("fldcw %0" : : "m"(newcword));
    long double result = x + y; // volatileが要るかも
    asm("fldcw %0" : : "m"(oldcword));
    return result;
}

x87 FPUがもともとコプロセッサだったことが関係するのかは知らないが、これらの命令に指定できるのはメモリ上のアドレスであり、汎用レジスタは指定できない。

SSE2

最近のx86向け（特にx86_64）にコンパイルしたコードでは、単精度・倍精度の浮動小数点数演算にはSSE2が使われることが多いと思う。

SSE2ではMXCSRレジスタを使って浮動小数点数の諸々を設定できる。

MXCSRレジスタをC言語からいじるには ﾅﾝﾄｶintrin.h を #include することで使えるようになる組み込み関数 _mm_getcsr, _mm_setcsr を使う。これらを使うと、丸めモードを指定して浮動小数点数の演算を行うコードは次のようになる：

double rounded_add(unsigned int mode, double x, double y) {
    unsigned int reg = _mm_getcsr();
    _mm_setcsr((reg & ~(3u << 13)) | (mode << 13));
    double result = x + y; // volatileが要るかもしれない
    _mm_setcsr(reg);
    return result;
}

余談だが、MXCSRレジスタを読み書きする命令（STMXCSR, LDMXCSR）もx87 FPUのそれと同じように、汎用レジスタではなく、メモリ上のアドレスを指定する必要がある。なんで？

ARM

ARMではFPSCRレジスタを使って浮動小数点数に関する諸々を設定できる。

GCCのビルトイン関数を使ってFPSCRレジスタを操作するコードの例は次のようになる：

double rounded_add(unsigned int mode, double x, double y) {
    unsigned int reg = __builtin_arm_get_fpscr();
    __builtin_arm_set_fpscr((reg & ~(3u << 22)) | (mode << 22));
    double result = x + y; // volatileが要るかも（ラズパイ上のGCCで試した時は必要だった）
    __builtin_arm_set_fpscr(reg);
    return result;
}

AArch64ではFPCRレジスタを使う。

double rounded_add(unsigned int mode, double x, double y) {
    unsigned int reg = __builtin_aarch64_get_fpcr();
    __builtin_aarch64_set_fpcr((reg & ~(3u << 22)) | (mode << 22));
    double result = x + y; // volatileが要るかも
    __builtin_aarch64_set_fpcr(reg);
    return result;
}

普段使いのパソコンがARM系という人はまだ少ないだろう。これらのコードを実機で試すには、ラズパイを使っておくのが手軽だと思う。上記のARM32向けのコードは筆者の家に転がっていた古いラズパイで動作を確かめた。筆者の家にはラズパイ3も多分転がっているので、それを使えばAArch64向けのコードも検証できるはずだが、発掘が面倒なのでAArch64向けのコードは試していない。

【2020年8月29日】コードを修正した。

動的な丸めモード指定の問題点

動的な丸めモード指定はC言語で標準化されていたり、多くの（？）プロセッサで実装されていたりするが、問題点もある。

まず、コンパイラの最適化との相性が悪い。

上述のコード例のコメントに「volatileが要るかもしれない」と書いたが、GCCやClang等のコンパイラは、volatileを指定しない場合実際の浮動小数点数演算 (x + y) と丸めモードの変更の順番を入れ替える可能性がある。

一応C言語では #pragma STDC FENV_ACCESS ON みたいなプラグマでコンパイラに注意を促すことができるようになってはいるが、GCCやClangはこのプラグマを実装していない。

性能の問題もある。

頻繁に丸めモードを変更する場合、実際の演算ではなく丸めモードを変更する命令の実行だけでそれなりに時間がかかってしまいそうだ（筆者は実測したわけではないので具体的な数値は出せない）。

これらの理由から、動的に丸めモードを指定するのではなく、静的に、命令自体で丸めモードを指定できることが望ましい。

擬似コードで書けば、

set_rounding_mode(UPWARD);
result = x + y;
set_rounding_mode(NEAREST);

ではなく

z = rounded_add(UPWARD, x, y); // あるいは add_up(x, y) みたいに関数自体を分ける

と書きたい。擬似コード中の rounded_add はCPUの命令に直接対応するような組み込み関数を想定している。

（MPFRの関数が丸めモードを引数として受け取るのに近い）

静的な丸めモード指定

最近のアーキテクチャでは、静的な丸めモード指定（命令自体に丸めモードを埋め込む）ができるものがある。また、プログラミング言語やコンパイラ基盤でも、静的な丸めモード指定に対応しようという動きがあるようだ（たとえアーキテクチャが静的な丸めモードに対応していなくても、処理系が正しい最適化をできるようになるという点で、プログラミング言語や中間言語で静的な丸めモード指定をできるようにする意義はある）。

【2020年8月29日追記】IEEE 754-2008以降では、”mode” という語はもっぱら動的な指定に使われ、「丸め方法」を含む浮動小数点演算に関する属性には “attribute” を使っている。「rounding-direction attributesを指定する」という感じである。ただ、この記事では執筆当時のまま「丸めモード指定」という言い方を使う。

AVX512 (EVEX prefix)

SSE4.1のROUNDSD等の命令（浮動小数点数を整数に変換する）では、（MXCSR等のレジスタを参照するのではなく）命令自体に丸めモードを埋め込めるようになっている。

AVX512のEVEX encodingではこれを発展させて、他の算術命令についても命令自体に丸めモードを埋め込めるようになった。

C言語でAVX512の機能を使って丸めモードを指定するには、 _mm_add_round_sd 等の組み込み関数を使う。第3引数で丸めモードを指定するのだが、ここはコンパイル時定数である必要がある。

double rounded_add_up(double a, double b) {
    __m128d av = _mm_set_sd(a);
    __m128d bv = _mm_set_sd(b);
    __m128d resultv = _mm_add_round_sd(av, bv, _MM_FROUND_TO_POS_INF | _MM_FROUND_NO_EXC);
    double result;
    _mm_store_sd(&result, resultv);
    return result;
}

長ったらしいコードに見えるが、この関数のコンパイル結果は実質 vaddsd {ru-sae}, %xmm1, %xmm0, %xmm0 の1命令となる。{ru-sae} が丸めモード指定の部分で、 ru は round upward, sae は suppress all exceptions の略だと思われる。

他の丸めモードを指定するには、C言語では _MM_FROUND_TO_NEAREST_INT, _MM_FROUND_TO_NEG_INF, _MM_FROUND_TO_POS_INF, _MM_FROUND_TO_ZERO 等を指定する。最近接丸めの指定が *_TO_NEAREST ではなく *_TO_NEAREST_INT なのが、ROUNDSD（浮動小数点数を整数に変換する）用に定義した定数を流用した感を醸し出している。

筆者を含めて、AVX512用のコードを動かせるマシンが手元にないという人も多いだろう。その場合は、Intelが出しているエミュレーター（sdeとかいうやつ）を使えば良い。あるいは、クラウドでAVX512が使えるマシンを借りるという手もあるだろう。

CUDA

CUDAには丸めモードを指定して演算できる組み込み関数／命令があるらしい。

GPGPUはあまりやったことがないのでここでは詳しくは触れない。

RISC-V

RISC-Vは最近流行りの、新しいプロセッサアーキテクチャである。後発なだけあって、いろいろ筋の良い設計をしている、らしい。

RISC-Vの仕様はここで見れる。RISC-Vの仕様はいくつかの細かい仕様に分かれていて、実装者が必要な仕様を選んで実装できる（不要な機能は実装しなくても良い）。整数の乗算もオプショナルである。

組み込み向けのマイコンでは浮動小数点数は省略されがちだが、RISC-Vの浮動小数点数も当然オプショナルである。実装する場合の仕様はF（単精度）, D（倍精度）, Q（四倍精度）のそれぞれがある。

RISC-VのF, D, Qの各仕様では、命令のうち3ビットを使って丸めモードを埋め込むことができる。

x86系やARMの制御レジスタでは丸めモードは4通りだったので2ビットで済んでいたが、4通りの丸めモードの他に「動的な丸めモードを使う」指定と、IEEE754-2008で規定された「roundTiesToAway（最近接丸め、ただし最も近い数が2通りある場合は絶対値が大きい方を採用）」を指定するために、3ビット使うようになっている。

（roundTiesToAwayに対応しているのがいかにも後発のアーキテクチャという感じがする。）

ただ、浮動小数点数の使えるRISC-Vマシンを個人が気軽に入手できる未来はもうちょっと先な気がする。HiFive Unleashedというやつは浮動小数点数も使えるみたいだが、個人にはちょっとハードルが高い。

また、プログラミング環境的な意味でも、丸めモード指定で遊べるのはもうちょっと先になりそうだ。GCCのビルトイン関数を見てもRISC-Vの丸めモード指定に対応しそうなものは見当たらない。

【2020年8月29日追記】GCCのインラインアセンブリーでRISC-Vの静的な丸め指定を使うと次のようになる。4番目のオペランドとして、 rne (roundTiesToEven), rtz (roundTowardZero), rdn (round down; roundTowardNegative), rup (round up; roundTowardPositive), rmm (round to nearest, ties to max magnitude; roundTiesToAway), dyn (dynamic rounding mode) のいずれかを指定する。

double add_up(double x, double y)
{
    double result;
    asm("fadd.d %0, %1, %2, rup" : "=f"(result) : "f"(x), "f"(y));
    return result;
}

double sub_up(double x, double y)
{
    double result;
    asm("fsub.d %0, %1, %2, rup" : "=f"(result) : "f"(x), "f"(y));
    return result;
}

テストするにはLinux上のQEMU User space emulatorを使うのが楽で、必要な諸々をインストールした後

$ riscv64-linux-gnu-gcc foo.c
$ qemu-riscv64 -L /usr/riscv64-linux-gnu/ ./a.out

という風にコンパイル・実行する。実機としては中国産のKendryte K210というSoCを搭載したマイコンボードやモジュールがいくつか出ているので、それを試してみるのも良いかもしれない。

C2x

我らがC言語でも、静的な丸めモード指定を標準化しようという動きがあるようだ（単体の技術仕様としてはTS 18661-1と2で、C2xにマージしたドラフトはn2314/n2315あたりを見れば良さそう）。

丸めモードの指定には（先述の擬似コードのような組み込み関数ではなく） #pragma を使用する。使用例は多分次のコードの感じ：

double add_up(double x, double y) {
    #pragma STDC FENV_ROUND FE_UPWARDS
    return x + y;
}

仕様にちゃんと目を通したわけではないのでアレだけど、Cコンパイラーは #pragma STDC FENV_ROUND に遭遇したら適切に fesetround 相当の呼び出しを挿入したり対応する機械語の命令を挿入したり、あるいはスコープ内の演算を静的な丸めモード指定付きの命令に置き換える、みたいなことをすることになるのだと思う。

四則演算とsqrtだけでなく、他の数学関数もこの #pragma の影響を受けるようだ。

LLVM intrinsics

LLVMの命令にも静的な丸めモードを指定できるものがあるようだ（ただし現段階では実験的なもので、今後どうなるかはわからない。そもそも現段階でどういう実装がされているのかもry）。

http://llvm.org/docs/LangRef.html#constrained-floating-point-intrinsics

Constrained Floating-Point Intrinsicsというやつに丸めモードを指定できるらしい。

これはC2xのやつとは違い、あくまで最適化のためのヒントであって、実際に実行時に丸めモードを変えるにはプログラマーは別途 fesetround かそれに相当する命令を実行する必要がある。…という風にドキュメントからは読める。違ってたらごめん。

Haskellライブラリーの進捗状況

以前の記事で最後に紹介した「Haskellで浮動小数点数の丸めモードを制御するライブラリー」だが、数ヶ月放置した後、最近またいじり始めている。

https://github.com/minoki/rounded-hw (wipブランチ)

（rounded-hw の関数は丸めモードを（グローバル変数ではなく）引数として指定できるようになっている。なので、Haskellのユーザーに対しては「静的な丸めモード指定」のできるインターフェースを提供していることになる。実行時の効率はともかくとして。）

現時点では、丸めモードの変更方法としては

C99
SSE2

に対応しており、余力があればAVX512やx87 FPU (long double用) にも対応したいと考えている。この記事に書いた方法でARMにも対応できるはずだが、筆者がARMでのHaskell環境を構築できていないので、対応するのは当分先のことだろう。

また、JavaScriptやJavaやWebAssemblyなど、非ネイティブなGHCバックエンドのことも想定して、FFIを使わない実装も用意する。現時点でのそれは「全部 Rational で計算してそれを自作関数により Double に丸める」実装なので劇遅だが、twoSum等のアルゴリズムを使えばもう少し高速化の余地はあるだろう。

とはいえ、初期リリースに関してはあまり欲張らずに、なるべく早く出したい。