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

C言語でのやり方

まず、C言語での丸めモードの変更方法を確認しておく。

標準Cで定められた方法としては、 fenv.h の fesetround/fegetround 関数を使う。

ただ、コンパイラーの最適化等はデフォルトでは丸めモードの変更を考慮していないかもしれない。コンパイラーに対しては #pragma STDC FENV_ACCESS ON という命令によって、プログラムが浮動小数点数の丸めモードを変更したり、その他FPUの状態にアクセスしたりする可能性があることを通知できる。

試しに、丸めモードを正の無限大方向にして足し算をする関数 roundeAddを書いてみよう。

#include <stdio.h>
#include <fenv.h>
#pragma STDC FENV_ACCESS ON
double roundedAdd(double x, double y)
{
    int r = fegetround();
    fesetround(FE_UPWARD);
    double z = x + y;
    fesetround(r);
    return z;
}
int main()
{
    printf("FE_UPWARD = %08x\n", FE_UPWARD);
    printf("FE_DOWNWARD = %08x\n", FE_DOWNWARD);
    printf("%a\n", roundedAdd(1, 0x1p-1000));
}

実行結果の例：

FE_UPWARD = 00000800
FE_DOWNWARD = 00000400
0x1.0000000000001p+0

main関数からは roundedAdd(1, 0x1p-1000) という呼び出しによって、 1+2^(-1000) を正方向へ丸めたものを計算させている。 1+2^(-1000) は倍精度浮動小数点数では正確に表現できず、一番近い倍精度浮動小数点数は 1 だが、正方向への丸めを指定しているので「1の次に大きい倍精度浮動小数点数」が返ってくる。

なお、Clangで上記のプログラムをコンパイルすると

test.c:3:14: warning: pragma STDC FENV_ACCESS ON is not supported, ignoring
      pragma [-Wunknown-pragmas]
#pragma STDC FENV_ACCESS ON
             ^
1 warning generated.

という警告が出る。そう、Clangは #pragma STDC FENV_ACCESS に対応していないのだ。状況はGCCでも同じで、同じような警告が出る。

つまりClangやGCCで上記のプログラムが期待通りに動作したのは偶然に過ぎない。

実際、筆者の環境では -O2 をつけて最適化を有効にしてコンパイルしたところ、 0x1p+0 という正しくない結果が得られた（roundedAdd 関数がインライン化されて演算が定数畳み込みされたことが原因だと考えられるので、 roundedAdd関数に __attribute__((noinline)) を指定したところ期待通りに動作した）。

ちなみに、GCCには -frounding-math というオプションでファイル単位で #pragma STDC FENV_ACCESS と同様の指定をできるらしいが、筆者の環境では上記のプログラムを -frounding-math 付きでコンパイルしても不正確な結果が返ってきた。

なお、この記事では簡単のために丸めモードの変更には fesetround を使うが、処理系によっては正しく実装されていなかったり、compiler intrinsicsを使った方が高速だったりする。この辺の環境依存な話については、柏木氏のページが参考になる： http://verifiedby.me/kv/rounding/

C FFIでfesetroundを呼ぶ

さっきのコードをHaskellに移植してみよう。fesetround はC FFIで呼べる。FE_UPWARD 等の定数は、移植性のある方法を考えるのが面倒なので筆者の環境での値を直接書き込んでしまう（ポータブルにしたかったら c2hs とかを使えば良いのだろうか）。

{-# LANGUAGE HexFloatLiterals #-}
import Numeric
import Foreign.C.Types

foreign import ccall unsafe "fesetround" c_fesetround :: CInt -> IO CInt
foreign import ccall unsafe "fegetround" c_fegetround :: IO CInt

roundedAdd :: Double -> Double -> IO Double
roundedAdd x y = do
  r <- c_fegetround
  c_fesetround 0x00000800 -- FE_UPWARD
  let z = x+y
  c_fesetround r
  return z

main = do
  x <- roundedAdd 1 0x1p-1000
  putStrLn $ showHFloat x ""

だが、これはもちろん期待通りには動作しない。Haskellに詳しい諸賢はどこが悪いか、全部指摘できるだろうか？

まず、遅延評価の問題がある。roundedAdd関数の中で let z しているが、 z の値はこの関数の中では必要ないので、roundedAdd関数の中では x+y は計算されない。roundedAdd関数を抜けてしまうと、浮動小数点数の丸めモードは元に戻っているので、デフォルトの最近接丸めで計算されてしまう。

これを解決するにはBangPatterns拡張を使って

let !z = x + y

と書く、という手がある。

遅延評価がらみでもう一つ、引数の方も正格評価しておく必要がある。roundedAdd (1+0x1p-1000) 0x1p-1000 と書いた時に 1+0x1p-1000 の方はHaskell標準のやり方（最近接丸め）で評価されて欲しいからだ。

これで、一応期待通りに動作するようになった。とはいえ、GHCの最適化は浮動小数点数の丸めモード変更を考慮していないだろうし、マルチスレッドの場合は丸めモードの変更が他のスレッドに波及する可能性がある。

詳しくはGHC User’s GuideのFFIの章を読んで頂きたい：https://downloads.haskell.org/~ghc/latest/docs/html/users_guide/ffi-chap.html#floating-point-and-the-ffi

そうすると丸めモードを変更する部分をHaskellで書くのは避けて、丸めモードを変更して戻す部分までをCで書き、それをFFIで呼び出す形にした方が良さそうである。

Cで実装した丸めモード指定付き演算をFFIで呼び出す

要するに、さっきのCで書いた roundedAdd 関数をFFIで呼び出せば、Haskellから「正方向に丸める足し算」を使えるということになる。

そういう風にして書いたCの関数をラップすれば、四則演算や平方根について、

data RoundingMode = TowardNearest | TowardInf | TowardNegInf | TowardZero
add :: RoundingMode -> Double -> Double -> Double
mul :: RoundingMode -> Double -> Double -> Double

という風なHaskellの関数を定義できる。

だが、しかし。話はこれでは終わらない。

丸めモードの変更というのは、プラットフォームにもよるかもしれないが、コストがかかる。変更回数（fesetroundの呼び出し回数）はなるべく少ない方が良い。

丸めモードの変更回数を最小化する、という観点からすると、演算ごとに「丸めモード変更→戻す」をするのは不利である。

例えば、区間演算で使う「正方向の丸めと負方向の丸めを計算する」というコードはC言語なら

#pragma STDC FENV_ACCESS ON
struct interval { double lo, hi; };
struct interval interval_add(struct interval a, struct interval b) {
    int r = fegetround();
    fesetround(FE_DOWNWARD);
    double lo = a.lo + b.lo;
    fesetround(FE_UPWARD);
    double hi = a.hi + b.hi;
    fesetround(r);
    return (struct interval){.lo = lo, .hi = hi};
}

と書ける（GCCやClangで期待通りに動作するとは言ってない）。このコードはfesetroundの呼び出しが3回である。一方、C言語で書いた「丸めモード指定付き加算」関数をそれぞれFFIで呼び出す方式だと、fesetroundが4回呼ばれることになる。

また、区間演算の乗算はもっと顕著で、C言語なら

#pragma STDC FENV_ACCESS ON
struct interval interval_mul(struct interval a, struct interval b) {
    int r = fegetround();
    fesetround(FE_DOWNWARD);
    double lo = fmin(a.lo * b.lo, fmin(a.lo * b.hi, fmin(a.hi * b.lo, a.hi * b.hi)));
    fesetround(FE_UPWARD);
    double hi = fmax(a.lo * b.lo, fmax(a.lo * b.hi, fmax(a.hi * b.lo, a.hi * b.hi)));
    fesetround(r);
    return (struct interval){.lo = lo, .hi = hi};
}

と書けるのに対し、FFIで「丸めモード指定付き乗算」をそれぞれ呼び出していたら、fesetroundを16回も呼び出す羽目になる。

もっと極端なことを言うと、配列の総和はC言語なら

#pragma STDC FENV_ACCESS ON
struct interval interval_sum(const struct interval *a, size_t n) {
    int r = fegetround();
    fesetround(FE_DOWNWARD);
    double lo = 0;
    for (size_t i = 0; i < n; ++i) {
        lo += a[i].lo;
    }
    fesetround(FE_UPWARD);
    double hi = 0;
    for (size_t i = 0; i < n; ++i) {
        hi += a[i].hi;
    }
    fesetround(r);
    return (struct interval){.lo = lo, .hi = hi};
}

と書けるが、FFIで「丸めモード指定付き加算」を呼び出す場合は、fesetroundを4n回も呼び出す羽目になる。

というわけで、Haskellで方向付き丸めのライブラリーを作りたかったら、C言語側で用意する関数は四則演算や平方根などのプリミティブな関数だけではなく、区間演算や総和などのもう少しハイレベルな操作も提供するべきである。

要は、ここで例示した interval_add 等の関数をHaskellからFFIで呼べばいいのだが、これらの関数はいくつかの理由から、C言語で実装するのには向いていない：

C言語のコードでは構造体返しを使っているが、GHCのC FFIはCの関数から返された構造体を受け取れない。ポインターを介さずに値をレジスターに乗せたまま返すには、アセンブリーでGHCの呼び出し規約に沿った関数を書く必要がある。
（こっちの方が致命的）ClangやGCC等のポピュラーなCの処理系は「関数の中で複数の異なる丸めモードを使うことを想定していない」ので間違ったコードを生成する。

まあ、流石にあれもこれもアセンブリー言語で書いていたら大変なので、区間の乗算や総和はそれぞれの丸めモードごとの関数をC言語で書くという形でも十分だろうとは思う。

ちなみに、GHCのFFIで構造体返しみたいなこと（ポインターを介さずに値をレジスターに乗せたまま返す）をするには、mod_poppo氏の次の記事が参考になる：

【低レベルHaskell】Haskell (GHC) でもインラインアセンブリに肉薄したい！

作っているもの

そういう感じで、最近は「ハードウェア組み込みの浮動小数点数に関する方向付き丸めをHaskellから使うためのライブラリー」を作っている。いつリリースできるかは不明。

https://github.com/minoki/rounded-hw

雑記帳

人生やっていき

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

C言語でのやり方

C FFIでfesetroundを呼ぶ

Cで実装した丸めモード指定付き演算をFFIで呼び出す

作っているもの

“Haskellでの浮動小数点数の方向付き丸めを考える” に1件のフィードバックがあります