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

以前、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(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(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向けのコードは試していない。

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

動的な丸めモード指定は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の関数が丸めモードを引数として受け取るのに近い)

静的な丸めモード指定

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

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の丸めモード指定に対応しそうなものは見当たらない。

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で浮動小数点数の丸めモードを制御するライブラリー」だが、数ヶ月放置した後、最近またいじり始めている。

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

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

  • C99
  • SSE2

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

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

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

TeXConf 2019向けのClutTeX発表資料

これは TeX & LaTeX Advent Calendar 2019 の10日目の記事です。9日目は munepi さんの ぼくのかんがえたさいきょうのLaTeX索引スタイルファイル でした。11日目は doraTeX さんです。

10月に予定されていたTeXConf 2019 が残念ながら台風によって中止になってしまいました。そこで、この記事では TeX 言語で台風の進路を制御する方法……ではなく、 TeXConf 2019 の発表で使う予定だった資料の供養をします。

続きを読む

Haskellのscan系関数を使いこなす

Haskellはリストを操作する関数を多数提供しています。map, filter, foldあたりが代表的で、これらは他の言語でもおなじみかと思います。

一方で、scan系関数(scanl, scanr)は他の言語ではあまり見かけない気がします。同じ関数型言語のSMLやOCamlにも標準では入っていないようです。

この記事では、scan系関数がどういう場合に利用できるかを紹介します。

続きを読む

技術書典7に、Haskellで競技プログラミングをやる本を出します

今週末の9月22日(日)に池袋で開催される技術書典7に、新刊「Haskellで戦う競技プログラミング」を出します。既刊「LaTeX処理自動化ツールClutTeX 使い方とその仕組み」も改訂して紙の本を頒布します。

技術書典7について、公式サイトより引用:

日時 2019/09/22 (日) 11:00〜17:00
場所 池袋サンシャインシティ 展示ホールC/D(文化会館ビル2/3F)
主催 TechBooster/達人出版会
一般入場は11:00~13:00のみ有料です。詳細はこちら

技術書典7

気になる方はサークル詳細からチェックリストに入れてください:

当サークルが配置されたのは「し03D」です。今回の技術書典は2Fの展示ホールDと3階の展示ホールCで行われますが、当サークルは2Fの展示ホールDです。会場が1箇所だった前回と同じ場所です。

続きを読む

TeXにとってやばい入力ファイル名

空白や記号類を名前に含むファイル名をTeXで処理させようとすると、うまく処理できないことがある。

例えば foo%bar.tex というファイルを latex foo%bar.tex のような感じでLaTeXに処理させようとすると、 % 以下の部分がコメント扱いされ、LaTeXは代わりに foo(.tex) というファイルを処理しようとする。

やばい名前のファイルを難しいことを考えずに処理できればいい方へ

8月中にリリースされる ClutTeX v0.4 を使って cluttex -e latex foo%bar.tex とすれば良い。

どういうファイル名がやばいのか、ClutTeX v0.4はどのように問題を解決するのか知りたい方へ

例によって結論に至るまでの経過をくどくど書いているので、結論だけ知りたい方は適当に読み飛ばすことをお勧めする。

続きを読む

AtCoderで青色になった

7月7日のAtCoder Beginner Contest 133で青色になりました。青色というのは上位7%に相当するらしいです。

「水色に上がった」記事からまだ1ヶ月ちょっとしか経っていないので、新しいHaskellネタも「青色に上がるためにこれを勉強した」みたいなネタもないのですが、適当に自分語りでもしてお茶を濁そうと思います。

続きを読む

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

週刊 代数的実数を作る」の #5 で、区間演算と方向付き丸めの話を書いた。浮動小数点数の計算は不正確だと思われがちだが、方向付き丸め等をうまく使えばある種の「正しい結果」(この数は確実に1.0より大きい、等)を得ることができる、という話だ。

MPFRのようなソフトウェア実装の浮動小数点数だと引数で丸めモードを指定できる。しかし、ハードウェア組み込み(注)の浮動小数点数(floatやdouble)の丸めモードを指定する方法は、言語や実行環境に依存する。

(注:環境によってはCPUに浮動小数点演算器が組み込まれておらず、floatやdoubleの演算もソフトウェア実装だったりするが、我々が普段使うPCではfloatやdoubleはほぼ確実にハードウェア実装されているため、以下「ハードウェア実装」で通す)

CやFortranのような低レベルかつ数値計算のニーズがあるようなプログラミング言語だと、丸めモードを変更する手段が用意されているが、JavaScriptやHaskellなど、そういうニーズが薄い言語では丸めモードの変更には対応していないことが多い。せいぜい、「浮動小数点数演算には最近接丸めを使用する」と規定されているのが関の山である。

しかし、「Haskellではできない処理がある」というのはどうにも気にくわない。どうにかして、ハードウェア組み込みの浮動小数点数の方向付き丸めをHaskellで扱う方法を考えたい。

ちなみに、ハードウェア組み込みの浮動小数点数にこだわらないのであれば、MPFRの丸めモード指定を使うHaskellライブラリーがすでに存在する:

http://hackage.haskell.org/package/rounded

続きを読む