IEEE 754 浮動小数点入門 - ushiostarfish’s diary

IEEE 754 浮動小数点

IEEE 754 浮動小数点は、本当は完全に扱えるはずもない実数をコンピューター上であたかもそれが実数を扱えるかのような錯覚を感じさせるのに成功している浮動小数点の標準規格です。過去にはたくさんの表現方法があったようですし、新しい表現方法が研究されていたりもするようですが、本記事では浮動小数点はすべて IEEE Standard for Floating-Point Arithmetic のことを指します。

レイトレーシングに限らずコンピューターで浮動小数点を使った数値計算は極めて広い応用がありますが、それにもかかわらず浮動小数点の負の側面、誤差やそれに対する傾向と対策についてはそれほどスポットが当たっていないように思います。そこで今回はゼロから浮動小数点の扱いについて考えてみようと思います。

内部表現

浮動小数点を理解するには、やはりまずその表現について把握することから始めるのが良いでしょう。浮動小数点の内部表現といっても、何ビットで表現するか？という問題があります。最も一般的なのは、32bitのbinary32, 64bitのbinary64でしょうか。C++ならそれに対応するのがfloat, doubleというのは極めて身近なデータ型でしょう。しかしC++はIEEE 754 浮動小数点が使われているかは処理系定義であることは注意が必要ですが、実際はそのように仮定しても問題ない環境は多いのではないでしょうか。ここではわかりやすさのため、binary32を仮定して話を進めます。テストコードについてはC++で行います。

binary32では浮動小数点は以下のように定義されます。

$\displaystyle { { \left( -1 \right) }^{ S }\ast 1.M\ast { 2 }^{ \left( E-127 \right) } }$

ここで $S$ は符号部で 0 か 1 の1bit 、 $M$ は二進数の仮数部で、Mが23bitつまり23桁あります。 Eは指数部で0から255、つまり8bitです。

wikipediaの図にするとビットの格納方法が一目でわかって便利です。

f:id:ushiostarfish:20190810214723p:plain ※https://ja.wikipedia.org/wiki/IEEE_754 より

...これで浮動小数点を完全に理解できるならばどんなに幸せでしょうか。しかしあきらめずに少しずつ紐解いていくことにします。

仮数部

先立って仮数部については簡単に補足します。仮数部は知っての通り、

f:id:ushiostarfish:20190810212226p:plain

のように10進数と対応します。ただしこの図では整数の部分が11の例ですが、浮動小数点では一番左は 1. となる決まりであり、0.とか11.とかはありえません。１つの理由としては、例えば 0.000001 を表現するのに左側のゼロ6個に6ビットつかうのはやはり無駄だからです。せっかく指数部があるのですから、小数点の位置を右に6回動かしてしまえば一番左は 1. になり、もしさらに右に桁があった場合には無駄なく使うことができます。

$\displaystyle { 0.000001=1.0\ast 2^{ -6 } }$

さらに1. が固定ならば、そこを除いて格納すれば節約になります。というわけで $1.M$ という表現になっています。ここで $1.M$ はどんな数かを少し考えてみます。最初に仮にMが4桁で[M=1010]だったとすると、

$\displaystyle { 1.M=1.1010=1+\frac { 0.1010\ast { 2 }^{ 4 } }{ { 2 }^{ 4 } } \\ =1+\frac { 1010 }{ { 2 }^{ 4 } } }$

ということは、binary32ではMは23桁であるので、

$\displaystyle { 1.M=1+\frac { M }{ { 2 }^{ 23 } } }$

という数だということがわかります。ということは、仮数部の項というのは、どうやらMをいじくると1以上2未満の間を行ったり来たりできることがわかります。

$\displaystyle { 1\le 1.M<2 }$

指数部

指数部は二進小数点の点の位置を移動させていると先に述べました。しかしこれだけではあまり脳内にイメージできないかもしれません。そこで数値が与えられたときどのような指数を持つのかを簡単に実験してみます。実験に先立ってfloatとuint32を相互変換できると便利ですので、以下の関数を準備します。

    inline uint32_t as_uint32(float f) {
        return *reinterpret_cast<uint32_t *>(&f);
    }
    inline float as_float(uint32_t u) {
        return *reinterpret_cast<float *>(&u);
    }

※ unionを使った再解釈は、C++では未定義動作とのことです

そして以下の関数で指数部を取り出してみます。

    inline uint8_t get_exponent(float f) {
        uint32_t e = as_uint32(f) & 0x7F800000;
        return e >> 23;
    }

-10から10くらいの間で浮動小数点を生成して、値をx、 $E-127$ すなわち指数部が2の何乗であるか？を y にしてプロットしたものが以下の図になります。格子の大きさは1です。

f:id:ushiostarfish:20190810221230p:plain

xが 0から1のエリアを見てみると、指数は0になっています。つまり、

$\displaystyle { 1\le 1.M\ast 2^{ 0 }<2\\ 1\le 1.M<2 }$

右隣はどうでしょうか？2から4のエリアです。

$\displaystyle { 1\ast 2^{ 1 }\le 1.M\ast 2^{ 1 }<2\ast 2^{ 1 }\\ 2\le 1.M\ast 2^{ 1 }<4 }$

なるほど、確かに式の通りです。左隣の0.5から1のエリアはどうでしょうか？

$\displaystyle { 1\ast 2^{ -1 }\le 1.M\ast 2^{ -1 }<2\ast 2^{ -1 }\\ \frac { 1 }{ 2 } \le 1.M\ast 2^{ -1 }<1 }$

元の範囲が1から2であるため、2倍もしくは1/2倍にしていけば確かに領域を重なることなく埋め尽くすことができます。なるほどよく考えられているものです。

f:id:ushiostarfish:20190810222808p:plain

このことから、指数部というのは担当領域、ある種窓のようなものだと捉えることができるわけです。

再び仮数部へ

指数部によって担当領域が決まること、そしてその幅が可変であることを確認しました。一方で仮数部の項は範囲が固定されているのも確認しました。

$\displaystyle { 1\le 1.M<2 }$

ここで表現できる数の精度限界について考えてみます。指数部と符号はあくまでも範囲を決めているわけなので、精度について重要なのはやはり仮数部になります。浮動小数点の定義から考えると、指数部で決まった範囲を仮数部が均等に刻んでいる、と捉えることができます。

例えばMが2ビットだった場合は、以下のように範囲が刻まれるでしょう。

f:id:ushiostarfish:20190810225114p:plain

実際にはbinary32では23ビットであるため、区間は $2^{23}$ 個に分割されることになります。ということは、例えば1から2の区間の長さは1であるため、最小幅は、

$\displaystyle { { 2 }^{ -23 } }$

となることがわかります。この数は特別に、”1より大きい最小の数”と1との差、として計算機イプシロン、もしくはマシンイプシロンと呼ばれています。ここでは https://github.com/catchorg/Catch2 を使って簡単なテストとして書いてみると、

TEST_CASE("machine_epsilon") {
    REQUIRE(std::ldexp(1.0f, -23) == std::numeric_limits<float>::epsilon());
}

と、このように実験で確認することができます。

しかし実際には指数部によって定義される区間の幅は可変長であるため、この最小幅も可変であることに気づきます。そこである数を絶対値が大きくなる方向へ1歩だけ動かす最小幅はどんな数であるか考えると、ある数のMを1インクリメントしたときの差分を考えればよさそうです。

$\displaystyle { 1.M_{ \left( +1 \right) }\ast { 2 }^{ \left( E-127 \right) }-1.M\ast { 2 }^{ \left( E-127 \right) }\\ =\left( 1.M_{ \left( +1 \right) }-1.M \right) \ast { 2 }^{ \left( E-127 \right) } }$

ここで、 $1.M$ が、

$\displaystyle { 1.M=1+\frac { M }{ { 2 }^{ 23 } } }$

であったのを思い出して、

$\displaystyle { \left( 1.M_{ \left( +1 \right) }-1.M \right) \ast { 2 }^{ \left( E-127 \right) }=\left( \frac { M_{ \left( +1 \right) } }{ { 2 }^{ 23 } } -\frac { M }{ { 2 }^{ 23 } } \right) \ast { 2 }^{ \left( E-127 \right) }\\ =\left( \frac { 1 }{ { 2 }^{ 23 } } \right) \ast { 2 }^{ \left( E-127 \right) }\\ ulp={ 2 }^{ \left( E-127-23 \right) }\\ ulp={ 2 }^{ \left( E-23-127 \right) }\\ E_{ ulp }=E-23\\ ulp={ 2 }^{ \left( E_{ ulp }-127 \right) } }$

とてもシンプルになりました。そしてこの数をulp(unit in the last place) と呼びます。仮数部の最小の桁の単位という気持ちですね。しかしこの式は浮動小数点の内部表現の $E$ に強く依存しており、若干泥臭いです。そこで少し違うアプローチで表現することを考えます。

ulpと似た概念でufp(unit in the first place)というのがあります。それは、

$\displaystyle { ufp\left( x \right) =2^{ \left\lfloor \log _{ 2 } \left| x \right| \right\rfloor } }$

ここで、 $\left\lfloor x \right\rfloor$ は床関数になります。これは何かというと、xの値の、指数部が定義する範囲の幅を表しています。

f:id:ushiostarfish:20190810232354p:plain

ということはulpは、その範囲の幅が仮数部によって分解された１つであることから、以下の式で定義できます。

$\displaystyle { ulp\left( E \right) ={ 2 }^{ \left( -23 \right) }ufp\left( x \right) }$

※ aは0でないとします。

このようにEの代わりに数値そのものと結びつけた表現として扱うこともできます。ついでにいうと左の項は計算機イプシロンになっていますね。

ここである実数 $a\in\mathbb{R}$ を最も近い浮動小数点に丸める関数を

$\displaystyle { RN\left( a \right) }$

とすると、丸めたときの誤差は、

$\displaystyle { \left| RN\left( a \right) -a \right| \le \frac { 1 }{ 2 } ulp\left( RN\left( a \right) \right) }$

が満たされることになります。ulpを使わない表現としては、

$\displaystyle { \left| RN\left( a \right) -a \right| \le \frac { 1 }{ 2 } { 2 }^{ \left( -23 \right) }ufp\left( a \right) }$

があり、このとき、

$\displaystyle { \frac { 1 }{ 2 } { 2 }^{ \left( -23 \right) }={ 2 }^{ \left( -24 \right) }=u }$

のuを単位相対丸め、という名前で使われることがあります。24は仮数部の桁数でもあるので、これをpとして、

$\displaystyle { { 2 }^{ \left( -p \right) }=u }$

と表現されることもあります。これを使うと先の式は、

$\displaystyle { \left| RN\left( a \right) -a \right| \le u\cdot ufp\left( a \right) }$

となります。ただし、計算機イプシロンとは違う値になるので、注意が必要です。例えば、GDC Vault - Math for Game Developers: Understanding and Tracing Numerical Errors in C++ ではその点に間違いがあったりします。※もちろんこれはこれで良い資料です。

以下は範囲とulpを表にしたものです。特筆すべき点としては、16777216以上ではもうulpが2以上になるため、もはや整数であっても表すことができない数字がここから出始めています。16777216は32bit整数型ではまだ全然余裕があることを考えると、これはかなりのデメリットです。

また、この表では指数関数的にulpが拡大するようにも見えますが、範囲も指数関数的に広がっているので、実際には表す数字とulpは比例関係（もちろん飛び飛びですが）になっていることには注意しましょう。

ゼロにならない？

浮動小数点の定義より、0に近い数字を表すためにはEにできるだけ小さな値を設定すれば良さそうです。そんなわけで最小値は、

$\displaystyle { 1.0\ast { 2 }^{ \left( 0-127 \right) }={ 2 }^{ \left( -127 \right) } }$

この値は大体、

$\displaystyle { 0.0000000000000000000000000000000000000058774717541114375398436826861112283890933277838604376075437585313920862972736358642578125 }$

なので、かなり0に近いは近いのですが、0ではありません。さらに問題なのは0までの距離です。E=0付近を図にすると、

f:id:ushiostarfish:20190811135706p:plain

$2^{-127}$ と $2^{-126}$ の間は仮数部により細かく刻まれていることを思い出すと、この0までの距離は ${ 2 }^{ 23 }=8388608$ 倍もの大きさの隙間であり、マイナスも考えると、その２倍の隙間が0の場所で存在することになります。そう考えると、これを無視するわけにはいきません。

非正規化数

そんなわけで、浮動小数点では、E=0のときだけ定義を変えてしまおう、ということになっています。それは以下の定義です。

$\displaystyle { { \left( -1 \right) }^{ S }\ast 0.M\ast { 2 }^{ \left( -126 \right) } }$

違いはわずかですが、最初の定義では最小の区間は、

$\displaystyle { 1\le 1.M<2\\ 1*{ 2 }^{ \left( -127 \right) }\le 1.M*{ 2 }^{ \left( -127 \right) }<2*{ 2 }^{ \left( -127 \right) }\\ { 2 }^{ \left( -127 \right) }\le 1.M*{ 2 }^{ \left( -127 \right) }<{ 2 }^{ \left( -126 \right) } }$

という区間でしたが、新しい定義では、

$\displaystyle { 0\le 0.M<1\\ 0*{ 2 }^{ \left( -126 \right) }\le 0.M*{ 2 }^{ \left( -126 \right) }<1*{ 2 }^{ \left( -126 \right) }\\ 0\le 0.M*{ 2 }^{ \left( -126 \right) }<{ 2 }^{ \left( -126 \right) } }$

そう、なんてことはない、０にもっとも近い部分の範囲だけ特別に０まで引き伸ばしているというわけです。

f:id:ushiostarfish:20190811140801p:plain

この特殊化された状態にある浮動小数点を、非正規化数と呼び、最初の定義のほうを正規化数と呼びます。これなら、0は、 $E=0, M=0$ で表すことができます。

$\displaystyle { 0.M*{ 2 }^{ \left( E-126 \right) }\\ =0.0*{ 2 }^{ \left( -126 \right) } }$

0付近でいきなり"穴"ができるのも防ぐこともできています。しかしいいことばかりではありません。代償は特殊化による演算パフォーマンスの低下であり、もしも非正規化数の取り扱いを辞めてもそれほど最終結果に影響が出ないと判断できる場合は、非正規化数になったとき0にフラッシュするようなモードを使ってこのパフォーマンス低下を回避することも考えられます。例えばembreeではそれをstrongly recommendedとしています。

https://www.embree.org/api.html#performance-recommendations

FLT_MIN, std::numeric_limits::min()

ところで FLT_MIN, std::numeric_limits::min() は負の数ではありません。これは正規化数で表現できる最小の正の数です。これは今までの議論によれば、

$\displaystyle { 1.0*{ 2 }^{ \left( -126 \right) } }$

ですので、

TEST_CASE("FLT_MIN") {
    REQUIRE(FLT_MIN == encode(PLUS_SIGN_BIT, 1, 0));
    REQUIRE(std::numeric_limits<float>::min() == encode(PLUS_SIGN_BIT, 1, 0));
    REQUIRE(std::numeric_limits<float>::min() == std::ldexp(1.0f, -126));
}

といったコードで確かめることができます。これは名前から勘違いしやすいところです。

-0 と +0

0近傍の話のついでに、0の符号についてついでに確認します。先の話では、0は、 $E=0, M=0$ で表すことができると述べました。この時実は符号は 0のときと1のときが存在します。これが-0 と + 0です。例えばcopysignといった関数ではちゃんと符号が反映されます。

TEST_CASE("plus_zero, minus_zero") {
    REQUIRE(-1.0f == std::copysign(1.0f, -0.0f));
    REQUIRE(+1.0f == std::copysign(1.0f, +0.0f));
}

注意しなければならない点としては、数としてはやはり0であることから、-0と+0はバイナリとして一致しないにもかかわらず、==で比較するとtureになるということです。

TEST_CASE("special_cases") {
    REQUIRE(-0.0f == +0.0f);
}

そしてこのことは大なり小なりでもそのように動作します。


TEST_CASE("special_cases") {
    REQUIRE(+0.0f >= -0.0f);
    REQUIRE(+0.0f <= -0.0f);
    REQUIRE((+0.0f < -0.0f) == false);
    REQUIRE((+0.0f > -0.0f) == false);
}

これは例えば、


    float x = -0.0f;
    if (0.0f <= x) {
        printf("%f is positive number?\n", x);
    }
    if (x <= -0.0f) {
        printf("%f is negative number?\n", x);
    }

というコードを試してみると、以下のような出力になってしまいます。

-0.000000 is positive number?
-0.000000 is negative number?

上のif文は排他的ではなかったのです。これはややエッジケースすぎる例かもしれませんが、意外と気づかれていないことかもしれません。

memset(&x, 0, sizeof(x));

ところでfloatやdoubleをmemsetなどでbitを全部0で埋めたらどうなるでしょうか？その場合は当然 $S=0, E=0, M=0$ ですから、結局値は+0です。

    float x;
    memset(&x, 0, sizeof(x));
    printf("%f\n", x); // 0.000000

このあたりも良く考えられていて、驚かされます。

例外的な数

$E=0$ のときの特殊化の話をしましたが、 $E=255$ のときも特殊化があり、それは-INF, +INF, NANです。

S	E	M
0	255	0	+INF
1	255	0	-INF
0~1	255	non 0	NAN

ただ、これについては今回はこれ以上踏み込みません

ulpを手に入れる

0付近の話が整理できたので、一度ulpの話に戻ることにします。floatのある数が与えられたとき、そのulpが欲しいときにはどうすればよいでしょうか。といってもさっきの議論によって、

$\displaystyle { E_{ ulp }=E-23 }$

であることがわかっているので、これを実装すると、

 const uint32_t MAX_FRACTION = 0x7FFFFF;
    const bool PLUS_SIGN_BIT  = false;
    const bool MINUS_SIGN_BIT = true;

    inline float encode(bool signbit, uint8_t expornent, uint32_t significand) {
        FT_ASSERT(significand < 0x800000);
        return as_float((signbit ? 0x80000000 : 0) | (static_cast<uint32_t>(expornent) << 23) | significand);
    }

    inline float ulp_only_normal_number(float f) {
        int32_t e = get_exponent(f);
        return encode(PLUS_SIGN_BIT, std::max(e - 23, 1), 0);
    }

ここで問題なのは、正規化数の範囲では、1が最小であり、

$\displaystyle { 1\le E_{ ulp }\\ 1\le E-23 }$

を満たさなければなりません。これは結局、

$\displaystyle { 24\le E }$

の数のulpだけが正規化数で表現できることを意味します。しかし非正規化数が使えるなら話は別です。非正規化数での表現を見つけるには、

$\displaystyle { 0.M\ast { 2 }^{ \left( -126 \right) }={ 2 }^{ \left( E-23-127 \right) } }$

を解いてMを求めれば良さそうです。

$\displaystyle { 0.M\ast { 2 }^{ \left( -126 \right) }={ 2 }^{ \left( E-23-127 \right) }\\ 0.M\ast { 2 }^{ \left( -126 \right) }=2^{ \left( E-24 \right) }\ast { 2 }^{ \left( -126 \right) }\\ 0.M={ 2 }^{ \left( E-24 \right) }\\ 0.M={ 2 }^{ \left( E-1 \right) }{ 2 }^{ \left( E-23 \right) }\\ \frac { M }{ { 2 }^{ 23 } } ={ 2 }^{ \left( E-1 \right) }\frac { 1 }{ { 2 }^{ 23 } } \\ M={ 2 }^{ \left( E-1 \right) } }$

求まりました。というわけで、実装は以下のようになります。

    inline float ulp(float f) {
        int32_t e = get_exponent(f);
        if (24 <= e) {
            return encode(PLUS_SIGN_BIT, e - 23, 0);
        }
        return encode(PLUS_SIGN_BIT, 0, 1 << std::max(e - 1, 0));
    }

e-1にmaxが必要なのは、非正規化数のulpを考えるときです。しかし非正規化数のulpは $E=1$ のときのulpと同じであることは、

f:id:ushiostarfish:20190811155338p:plain

図から明らかです。したがって、maxでガードすれば良いということになります。"floatのulpはfloatで表現することができる"、というわけですね。

floatの隣の数

実数には隣の数はありませんが、floatには隣の数があります。今までの議論から、あるfloatの数の隣の数が欲しいとき、仮数部の数値を増やしたり減らしたりすれば手にはいりそうですが、仮数部がオーバーフローやアンダーフローすることがあります。ところがここでfloatの内部表現を再確認すると、

f:id:ushiostarfish:20190810214723p:plain ※https://ja.wikipedia.org/wiki/IEEE_754 より

指数部が仮数部の上位ビットとして存在しています。ということは全体を整数として解釈するならば、仮数部がオーバーフローやアンダーフローしても、自動的に指数部がインクリメント・デクリメントされるだけであることに気づきます。つまり隣の数が欲しければ整数として解釈してインクリメント・デクリメントしても良いということです。

ただ残念なことに、正のfloatは良いのですが、負のfloatの場合はインクリメントすると絶対値が大きくなる方向に進みますが、一方intをインクリメントすると、絶対値は小さくなります。この対応が逆であることは注意しなければなりません。

2の補数 - Wikipedia

※ただ、正の値のみを考えれば良い問題も多いので、状況によってこのことは無視できます。

この点に注意しながら、floatを順序づけされたintと相互変換する関数を考えてみます。

    inline int32_t to_ordered(float f) {
        uint32_t b = as_uint32(f);

        uint32_t s = b & 0x80000000; // sign bit
        int32_t  x = b & 0x7FFFFFFF; // expornent and significand

        return s ? -x : x;
    }

    inline float from_ordered(int32_t ordered) {
        if (ordered < 0) {
            uint32_t x = -ordered;
            return as_float(x | 0x80000000);
        }
        return as_float(ordered);
    }

これを使うと、floatを任意の歩数だけ移動させる関数は、

    inline float next_float(float x, int move) {
        return from_ordered(to_ordered(x) + move);
    }

と簡単に実装できます。ただこの実装では-0.0と+0.0は両方とも 0にマッピングされてしまうところは注意が必要です。ただこれはメリットでもあります。なぜなら-0の次の数は+0ではなく、 $S=0, E=0, M=1$ の数であるべきだからです。これは標準の関数である std::nextafter の挙動もそうなっています。もし-0の次が+0だと、数として全く進んでいないことになります。

TEST_CASE("nextafter") {
    REQUIRE(std::nextafter(-0.0f, +std::numeric_limits<float>::max()) != +0.0f);
    REQUIRE(std::nextafter(-0.0f, +std::numeric_limits<float>::max()) == (-0.0f + ulp(-0.0f)));
}

floatの順序性の応用

それって本当に便利なの？という人もいるかもしれません。なので一つここで応用例を挙げます。例えばOpenCLではfloatに対してatomic な max関数がありません。しかし一方で、uintやintに対してはatomic maxが用意されています。

atomic_max

したがって、先のto_ordered関数を使ったり、もし扱う数がすべて正の数に限定できればもっと単純にuintとしてatomic_maxを使っても結果的に正しく処理できるはずです。私独自のラッパーを使っており、少し読みにくいかもしれませんが、例えば以下のようになります。

__kernel void find_max(__global float *values, __global uint *max_value) {
    __local uint local_max;
    if(get_local_id(0) == 0) { 
        local_max = 0;
    }

    barrier(CLK_LOCAL_MEM_FENCE);

    float value = values[get_global_id(0)];
    atomic_max(&local_max, as_uint(value));
    
    barrier(CLK_LOCAL_MEM_FENCE);

    if(get_local_id(0) == 0) { 
        atomic_max(max_value, local_max);
    }
}

TEST_CASE("max_value") {
    using namespace rt;
    auto &env = OpenCLProgramEnvioronment::instance();
    env.setSourceDirectory("../kernels");

    rt::Xoshiro128StarStar random;

    OpenCLContext context;

    int deviceCount = context.deviceCount();
    for (int device_index = 0; device_index < deviceCount; ++device_index) {
        // INFO("device name : " << context.device_info(device_index).name);

        auto lane = context.lane(device_index);

        OpenCLProgram program("find_max.cl", lane.context, lane.device_id);
        OpenCLKernel kernel("find_max", program.program());

        for (int j = 0; j < 100; ++j) {
            std::vector<float> values;
            for (int i = 0; i < 1000; ++i) {
                values.push_back(random.uniform(0.0f, 1000.0f));
            }
            float min_value_truth = 0.0f;
            for (auto f : values) {
                min_value_truth = std::max(min_value_truth, f);
            }
            OpenCLBuffer<float> values_gpu(lane.context, values.data(), values.size(), rt::OpenCLKernelBufferMode::ReadWrite);

            float min_value = 0.0f;
            OpenCLBuffer<uint32_t> min_value_gpu(lane.context, (uint32_t *)(&min_value), 1, rt::OpenCLKernelBufferMode::ReadWrite);

            kernel.setArguments(values_gpu.memory(), min_value_gpu.memory());
            kernel.launch(lane.queue, 0, values.size());

            min_value_gpu.read_immediately((uint32_t *)(&min_value), lane.queue);

            REQUIRE(min_value == min_value_truth);
        }
    }
}

※0以上を仮定しています

このようなハックが可能なのは内部表現の規則のおかげです。

整数

浮動小数点は誤差がある、誤差に気をつけなさい、と耳にタコができるくらいに聞きすぎると、整数の演算にも丸め誤差などが入り込むような錯覚を覚えることがあるかもしれません。しかし仮数部はbinary32の場合二進数で24桁あり、この範囲では丸め誤差なく完全に表現することができます。

※仮数部は、 $1.M$ となっているためビット数よりも桁が一つ多いです

1100110......10.
[       24桁      ]

せいぜい23歩程度の小数点の移動であり、指数部のビットが足りなくなることもありませんね。24桁ですべて1なら、16777215になります。しかし、 $16777216=2^{24}$ であるため、-16777216から16777216までが完全に表現できる最大ということになります。16777217は残念ながら表現できません。ulpが2だからです。この範囲であれば例えば、

    int x = -16777216; /* 2^24 */
    for (float i = x; i <= 16777215.0f /* 2^24 - 1 */; i += 1.0f) {
        REQUIRE(i == x);
        ++x;
    }
    REQUIRE(x == 16777216);

は意図通りに動作します。このように整数をちゃんと表せることが、JavaScriptの数値は64ビットの浮動小数点数だけ、というのを正当化する一つの理由となっています。

誤差

誤差については今まで真の値を丸めたときの誤差しか触れていませんでした。浮動小数点演算では誤差が発生する、と多くの人は認識しており、floatの比較には気を付けましょう、といったことや、答えが0に近づく引き算をしない、などのTipを聞くことはよくあります。有名な既存手法としては、以下のものがあります。

区間演算

floatやdoubleや、より高価な演算を使っても誤差が出るものは出る。ならば正確な答えを出すことよりも、正確な答えのある範囲を確実に求めよう、というのが区間演算のアイディアです。IEEE 754 浮動小数点では四則演算とsqrtは、無限の精度で計算され、それが丸められて答えになるように定められている事実を利用して、精度の保証を行います。

区間演算と交差判定

精度保証付き数値計算の必要性

エラーフリー変換

$\displaystyle { fl\left( a\circ b \right) }$

を、なにかの演算を無限の精度で行ってそれを浮動小数点に丸める計算だとします。多くのシチュエーションでは、

$\displaystyle { fl\left( a+b \right) \neq a+b }$

のように丸め誤差のために浮動小数点で真の値を得ることができません。しかし、

$\displaystyle { fl\left( a+b \right) =x+y }$

ここでxはなるべく絶対値が大きな数、yは絶対値がなるべく小さい数です。このような変換ができることが知られています。これを利用してあたかも倍の精度で演算が行えるような計算手法が提案されています。

Stef Graillat and Valérie Ménissier-Morain, "Error-free transformations in real and complex floating point arithmetic"

あたかも倍の精度の演算が行えることから double-double演算とも呼ばれるそうです。カハンの加算アルゴリズムもこの手法にアイディアは似ています。

しかし残念ながらこれらの手法はどのように計算をどのように改善すれば良いか？に直接答えてくれるわけではありません。なのでここでは紹介にとどめます。

Practically Accurate Floating-Point Math

そこで、既存の計算アルゴリズムをどのように改善したらよいのか？という視点では、

Neil Toronto and Jay McCarthy | Brigham Young University, "Practically Accurate Floating-Point Math"

が大変よくまとまっており、有益だと感じました。そこで今回はこちらを紹介したいと思います。

エラーの定義

まずはエラーとはどういうものかを定義します。

$\displaystyle { error\left( x,r \right) =\frac { \left| x-r \right| }{ ulp\left( fl\left( r \right) \right) } }$

ここでxは誤差があるかもしれない浮動小数点数、rは真の値です。真の値から何ulpずれているのか？というのがエラーである、というわけです。ulpは相対的な値ですから、相対誤差にも近い概念です。

例えば以下のような関数でいつでも評価できるようにすると便利です。

    template <class T>
    float flulp_error(float x, T r) {
        T e = (T(x) - r) / T(ulp(float(r)));
        return std::fabs(float(e));
    }

Racketではリファレンスとなる値の計算のためにMPFRを使っているそうです。今回私はwindows, C++の環境のため、vcpkg から、

vcpkg install mpfr:x64-windows
vcpkg install mpir:x64-windows

を入れたあと、boost::multiprecision を使います。

    using namespace boost::multiprecision;
    using Real = mpfr_float_1000;

    printf("%f\n", flulp_error(next_float(1.0f, 4), Real(1.0f))); // 4.000000

このように使うことができます。ほかにも例えばfloatのexpについて調べたい場合、

    using namespace std;
    using namespace boost::multiprecision;

    printf("%f\n", flulp_error<mpfr_float_1000>(exp(10.0f), exp(mpfr_float_1000(10.0f)))); //

この例では実行すると0.486941が得られます。これは0.5を超えないため、なかなかいい感じではないでしょうか。無限の精度で演算した後に丸められたとしても、0.5ulpはずれるのは回避できないからです。

floatを調べるのにdoubleはリファレンスとして使えるか？

floatよりもdouble, doubleよりもlong doubleで計算したものをリファレンスとするのは、おそらく多くのケースでうまくいくでしょう。しかし反例もあります。

以下の関数があり、 $x=77617,y=33096$ の時の値を求めます。

$\displaystyle { f\left( x,y \right) =333.75{ y }^{ 6 }-x^{ 4 }{ y }^{ 6 }+x^{ 2 }\left( 11x^{ 2 }{ y }^{ 2 }-121{ y }^{ 4 }-2 \right) +5.5{ y }^{ 8 }+\frac { x }{ 2y } }$

これをC++で組むと、

template <class T>
T f(T x, T y) {
    auto x2 = x * x;
    auto y2 = y * y;
    auto y4 = y2 * y2;
    auto y6 = y4 * y2;
    auto y8 = y4 * y4;
    return (T(333.75) * y6 - x2 * y6) + x2 * (T(11.0) * x2 * y2 - T(121.0) * y4 - T(2.0)) + T(5.5) * y8 + x / (T(2.0) * y);
}

そして、型を変えながら実行してみます。

    float a = 77617.0f;
    float b = 33096.0f;
    printf("%.30f, float  (8)\n", (float)f<float>(a, b));
    printf("%.30f, double (15)\n", (float)f<double>(a, b));
    printf("%.30f, mpfr   (30)\n", (float)f<number<mpfr_float_backend<30> >>(a, b));
    printf("%.30f, mpfr   (50)\n", (float)f<mpfr_float_50>(a, b));
    printf("%.30f, mpfr   (100)\n", (float)f<mpfr_float_100>(a, b));
    printf("%.30f, mpfr   (1000)\n", (float)f<mpfr_float_1000>(a, b));

すると、以下の出力が得られました。

1.172603964805603027343750000000, float  (8)
1.172603964805603027343750000000, double (15)
1.172603964805603027343750000000, mpfr   (30)
-0.827396035194396972656250000000, mpfr   (50)
-0.827396035194396972656250000000, mpfr   (100)
-0.827396035194396972656250000000, mpfr   (1000)

※ 括弧の中は10進数での桁数

驚くべきことに、long doubleまで上げても、変化がないので真実かと思いきや、さらに精度を上げていくと突然符号が変化します。ちょっと作為的な例ですが、浮動小数点の信用ならないふるまいを垣間見ることができます。これはRumpの例題というらしいです。"測定"を何度も行って真実を見極めるのは科学として自然なことですが、数値計算でも条件を変えて色々な角度から観察することは有効ではないでしょうか。

エラーの伝播

何か新しい関数を実装するには、多くの場合、浮動小数点の既存の演算を組み合わせて行われます。ある関数が、もし無限の精度で中間の演算が行われ、最後だけ丸められるとしたらその関数は最大でも0.5ulpの誤差しか生まない極めて理想的な関数であると言えます。しかしながら現実的にはパフォーマンス上の問題もありすべての関数をそのように組み立てるのは困難です。一方で演算→丸め→演算→丸め→演算→丸め→...を繰り返しているとどんどん誤差は蓄積していき、最終的には演算の意味消失まであり得ます。青天井です。ではどこで誤差は拡大するのでしょうか？

これを観察するために、もし入力値に1ulpの誤差があったとき、結果にどの程度影響するのか？を以下のようにテストしてみます。

printf("%f\n", flulp_error<mpfr_float_1000>(log(next_float(1.1f, 1)), log(mpfr_float_1000(1.1f))));

// 14.302301

もうちょっと0に近づいてみます。

printf("%f\n", flulp_error<mpfr_float_1000>(log(next_float(1.01f, 1)), log(mpfr_float_1000(1.01f))));

// 126.738930

さらに

printf("%f\n", flulp_error<mpfr_float_1000>(log(next_float(1.0001f, 1)), log(mpfr_float_1000(1.0001f))));

// 16382.376953

入力値が1に近づけば近づくほど、爆発的にエラーが増加しているのがわかります。これはちょうどx軸との交点、結果が0に近づく場所ですね。

f:id:ushiostarfish:20190812133803p:plain

これは結局、入力値に1ulpの誤差が含まれただけで、結果の下位ビットの多くの桁が無意味化するということです。これは計算を組み合わせるときに大いに問題になります。ただしこれはあくまでも入力値に誤差が含まれた時の話です。もし入力値が真実に十分近ければ問題にならないこともあります。例えば以下のように

    printf("%f\n", flulp_error<mpfr_float_1000>(log(1.1f), log(mpfr_float_1000(1.1f))));
    printf("%f\n", flulp_error<mpfr_float_1000>(log(1.01f), log(mpfr_float_1000(1.01f))));
    printf("%f\n", flulp_error<mpfr_float_1000>(log(1.0001f), log(mpfr_float_1000(1.0001f))));

        // 0.302301
        // 0.261070
        // 0.377019

全然問題ありません。伝播のしやすさともともとの関数の精度はまったく関係がないのです。また、時々桁落ちの誤差の説明で、以下のような説明を目にします。

1.3334-1.3333 を計算したいが、答えは 0.0001 になってしまい、有効数字が5桁から1桁に減少する。

なんだか狐につままれたような感じがします。はっきりしているのは、もしもこの二つの入力値が真実なら、答えの0.0001はやはり真実だということです。実際、問題は入力値が丸め誤差などによって信頼できない場合です。

例えば真実の下位の桁が [ ] のような数だったとすると、

1.3334 [ 439475 ] - 1.3332 [ 635621 ] = 0.0001803854

最後に丸めても、

0.00018039

このように最終結果に本来反映されるべきは [ ] の中だったわけで、真実にかなり近い数を保持できるポテンシャルがありながら、丸められた入力値で計算したことで 0.0001のさらに下の桁が失われてしまったわけです。

伝播は誰の仕業か？

ところでこの伝播はどこで拡大するのでしょうか？いったんエラーの定義を確認します。

$\displaystyle { error\left( x,r \right) =\frac { \left| x-r \right| }{ ulp\left( fr\left( r \right) \right) } }$

ここでulpは以下のように近似できると仮定します。

$\displaystyle { ulp\left( x \right) \simeq ulp\left( 1 \right) \left| x \right| }$

このとき、入力値の 1 ulpのずれを観察します。

$\displaystyle { x=f\left( x+ulp\left( x \right) \right) \\ r=f\left( x \right) \\ error\left( x,r \right) =\frac { \left| f\left( x+ulp\left( x \right) \right) -f\left( x \right) \right| }{ ulp\left( fr\left( f\left( x \right) \right) \right) } \\ =\frac { \left| f\left( x+ulp\left( x \right) \right) -f\left( x \right) \right| }{ ulp\left( 1 \right) \left| fr\left( f\left( x \right) \right) \right| } \\ =\frac { \left| x \right| \cdot \left| f\left( x+ulp\left( x \right) \right) -f\left( x \right) \right| }{ \left| x \right| \cdot ulp\left( 1 \right) \left| fr\left( f\left( x \right) \right) \right| } \\ =\frac { \left| x \right| \cdot \left| f\left( x+ulp\left( x \right) \right) -f\left( x \right) \right| }{ ulp\left( x \right) \left| fr\left( f\left( x \right) \right) \right| } }$

ここで、

$\displaystyle { \Delta x=ulp\left( x \right) }$

と捉えると、

$\displaystyle { =\frac { \left| x \right| \cdot \left| f\left( x+\Delta x \right) -f\left( x \right) \right| }{ \left| fr\left( f\left( x \right) \right) \right| \Delta x } }$

これを微分として解釈してしまうと、

$\displaystyle { error\left( x,r \right) =\frac { \left| x \right| \cdot \left| f'\left( x \right) \right| }{ \left| fr\left( f\left( x \right) \right) \right| } \\ =\left| \frac { x\cdot f'\left( x \right) }{ fr\left( f\left( x \right) \right) } \right| }$

最後に浮動小数点丸めを見なかったことにしてもさほど影響がないので無視すると、

$\displaystyle { error\left( x,r \right) =\left| \frac { x\cdot f'\left( x \right) }{ f\left( x \right) } \right| }$

ここからわかることは、このエラーの伝播の性質というのは浮動小数点の問題ではなく、本来その関数が持つ性質であったということです。逆に考えると、エラーの伝播のしやすさは、計算アルゴリズムを工夫しても改善できるものではないということです。我々ができるのは、エラーの伝播＆拡大しやすいところになるべく立ち入らず、立ち入る場合は誤差をなるべく少なくした状態で入るようにする、といったことなのです。

バッドランド

エラーの伝播＆拡大しやすいところを把握することは重要です。Neil Toronto and Jay McCarthy | Brigham Young University, "Practically Accurate Floating-Point Math" では、

$\displaystyle { 4\le \left| \frac { x\cdot f'\left( x \right) }{ f\left( x \right) } \right| }$

の領域を、エラーの伝播＆拡大のしやすいところとし、バッドランド(badlands)と呼んでいます。荒地、などと呼んでもいいかもしれませんが、逆にわかりにくい気がするのでカタカナで呼びます。

例えば先のlogのバッドランドを見てみます。gnuplotは割と使いやすいです。

class GNUPlot2D {
public:
    GNUPlot2D() {
        _fp = std::shared_ptr<FILE>(_popen("gnuplot", "w"), _pclose);
        gnu_printf("$data << EOD\n");
    }
    void add(double x, double y) {
        gnu_printf("%f %f\n", x, y);
    }
    void show(const char *label, double zMin, double zMax) {
        gnu_printf("EOD\n");
    
        gnu_printf("set grid\n");
        gnu_printf("set yrange[%e:%e]\n", zMin, zMax);
        gnu_printf("plot '$data' using 1:2 with lines title '%s'\n", label);

        flush();
    }
private:
    void gnu_printf(const char *format, ...) {
        va_list ap;

        va_start(ap, format);
        vfprintf(_fp.get(), format, ap);
        va_end(ap);
    }
    void flush() {
        fflush(_fp.get());
    }
    std::shared_ptr<FILE> _fp;
};

    GNUPlot2D p;
    int N = 1000;
    for (int i = 0; i < N; ++i) {
        float x = glm::mix(0.0f, 2.0f, float(i) / N);
        float e = flulp_error<mpfr_float_1000>(log(next_float(x, 1)), log(mpfr_float_1000(x)));
        p.add(x, e);
    }
    p.show("log", 0.0, 100.0);

f:id:ushiostarfish:20190812172400p:plain

x が1に近い所でバッドランドが観測できます。また、desmosのようなツールを使ってバッドランドを確認するのも良いでしょう。

f:id:ushiostarfish:20190812184639p:plain

ちなみにlogの場合は、バッドランドの定義より、解析的に範囲を求めることができます。

$\displaystyle { \left| \frac { x\cdot f'\left( x \right) }{ f\left( x \right) } \right| \le 4\\ \left| \frac { x\cdot \frac { 1 }{ x } }{ \log { \left( x \right) } } \right| \le 4\\ \left| \frac { 1 }{ \log { \left( x \right) } } \right| \le 4\\ \left| \log { \left( x \right) } \right| \le \frac { 1 }{ 4 } \\ -\log { \left( x \right) } \le \frac { 1 }{ 4 } ,\log { \left( x \right) } \le \frac { 1 }{ 4 } \\ x\ge exp\left( -\frac { 1 }{ 4 } \right) ,x\le exp\left( \frac { 1 }{ 4 } \right) \\ exp\left( -\frac { 1 }{ 4 } \right) \le x\le exp\left( \frac { 1 }{ 4 } \right) }$

続いて、桁落ちの原因とされる x - y のバッドランドも見てみましょう。


class GNUPlot3D {
public:
    GNUPlot3D() {
        _fp = std::shared_ptr<FILE>(_popen("gnuplot", "w"), _pclose);
        gnu_printf("$data << EOD\n");
    }
    void add(double x, double y, double z) {
        gnu_printf("%e %e %e\n", x, y, z);
    }
    void next_row() {
        gnu_printf("\n");
    }
    void show(const char *label, double zMin, double zMax) {
        gnu_printf("EOD\n");
        
        gnu_printf("set pm3d\n");
        gnu_printf("set cbrange[0:8]\n");
        gnu_printf("set palette defined(-1 \"blue\", 0 \"white\", 1 \"red\")\n");
        
        gnu_printf("set ticslevel 0\n");
        gnu_printf("set zrange[%e:%e]\n", zMin, zMax);
        gnu_printf("set xlabel 'x'\n");
        gnu_printf("set ylabel 'y'\n");
        gnu_printf("splot '$data' using 1:2:3 w pm3d title '%s'\n", label);
        
        flush();
    }
private:
    void gnu_printf(const char *format, ...) {
        va_list ap;

        va_start(ap, format);
        vfprintf(_fp.get(), format, ap);
        va_end(ap);
    }
    void flush() {
        fflush(_fp.get());
    }
private:
    std::shared_ptr<FILE> _fp;
};

    GNUPlot3D p;
    int N = 64;
    for (int i = 0; i < N; ++i) {
        float y = glm::mix(-1.0f, 1.0f, float(i) / N);

        for (int j = 0; j < N; ++j) {
            float x = glm::mix(-1.0f, 1.0f, float(j) / N);

            float value = next_float(x, 1) - next_float(y, -1);
            mpfr_float_1000 truth = mpfr_float_1000(x) - mpfr_float_1000(y);
            float e = flulp_error<mpfr_float_1000>(value, truth);
            p.add(x, y, e);
        }
        p.next_row();
    }
    p.show("x - y", 0.0, 16.0);

f:id:ushiostarfish:20190812172540p:plain

0付近は値が大きいので飛んでしまっていますが、xとyが近い値のとき、バッドランドになっていることがわかります。他にも Neil Toronto and Jay McCarthy | Brigham Young University, "Practically Accurate Floating-Point Math" で既にプロットされたものがあるので、覗いてみるのも良いでしょう。

また、便利なバッドランドリストもあります。

f:id:ushiostarfish:20190812183148p:plain ※ Figure 6 Neil Toronto and Jay McCarthy | Brigham Young University, "Practically Accurate Floating-Point Math"

関数を改善する

今までの知識を用いて、以下の関数を考えます。ただし簡単のため、x, yは正の数とします。

$\displaystyle { f\left( x,y \right) =\log { \left( \frac { x }{ y } \right) } }$

まずはこれがどれだけの精度があるのかを確認してみます。

float f(float x, float y) {
    return log(x / y);
}

mpfr_float_1000 f_truth(mpfr_float_1000 x, mpfr_float_1000 y) {
    return log(x / y);
}

    GNUPlot3D p;
    int N = 64;
    for (int i = 0; i < N; ++i) {
        float y = float(i) / N;

        for (int j = 0; j < N; ++j) {
            float x = float(j) / N;

            float value = f(x, y);
            mpfr_float_1000 truth = f_truth(x, y);
            float e = flulp_error<mpfr_float_1000>(value, truth);
            p.add(x, y, e);
        }
        p.next_row();
    }
    p.show("log(x/y)", 0.0, 16.0);

f:id:ushiostarfish:20190812183514p:plain

$\displaystyle { f\left( x,y \right) =\log { \left( \frac { x }{ y } \right) } }$

なかなか誤差が大きいエリアがあることが見て取れます。どうやら $\frac { x }{ y }$ が1に近いエリアでどうやら数値が暴れています。これは先に述べたlogのバッドランドそのものです。割り算の結果が丸められ、その誤差がバッドランドで伝播・拡大したことが原因と考えられます。バッドランドリストによれば、logのバッドランドは、

$\displaystyle { \exp\left(-\frac{1}{4}\right)\le x\le\exp\left(\frac{1}{4}\right) }$

のエリアがバッドランドであり、このエリアだけ回避するのが良いと思われます。そこで、ここではlogを直接使うのではなく、以下の関数を活用することを考えます。

$\displaystyle { \log _{ 1p }{ \left( z \right) } =\log { \left( z+1 \right) } }$

どうやらいろんな環境でも提供されている関数のようです。

std::log1p, std::log1pf, std::log1pl - cppreference.com

これのバッドランドは、オリジナルとは少しずれています。

f:id:ushiostarfish:20190812190553p:plain

ここで、これを使うように式を変形します。

$\displaystyle { f\left( x,y \right) =\log { \left( \frac { x }{ y } +1-1 \right) } \\ =\log _{ 1p }{ \left( \frac { x }{ y } -1 \right) } \\ =\log _{ 1p }{ \left( \frac { x-y }{ y } \right) } }$

このとき、 $\frac { x-y }{ y }$ は、 $exp\left( -\frac { 1 }{ 4 } \right) \le \frac { x }{ y } \le exp\left( \frac { 1 }{ 4 } \right)$ であるとき、

$\displaystyle { exp\left( -\frac { 1 }{ 4 } \right) \le \frac { x }{ y } \le exp\left( \frac { 1 }{ 4 } \right) \\ exp\left( -\frac { 1 }{ 4 } \right) -1\le \frac { x }{ y } -1\le exp\left( \frac { 1 }{ 4 } \right) -1\\ exp\left( -\frac { 1 }{ 4 } \right) -1\le \frac { x-y }{ y } \le exp\left( \frac { 1 }{ 4 } \right) -1 }$

となりますので、これが $\log _{ 1p }$ のバッドランドを踏んでいなければ良さそうです。確認してみます。

f:id:ushiostarfish:20190812191527p:plain

エクセレント！大丈夫のようです。これを試してみましょう。

float f_enhanced(float x, float y) {
    static const float lower = std::exp(-0.25f);
    static const float upper = std::exp(+0.25f);

    float x_div_y = x / y;
    if (lower <= x_div_y && x_div_y < upper) {
        return log1p((x - y) / y);
    }
    return log(x_div_y);
}

f:id:ushiostarfish:20190812192151p:plain

2ulpを超える所がすべて排除されています！

ところでじゃあ最初から、

$\displaystyle { f\left( x,y \right) =\log _{ 1p }{ \left( \frac { x-y }{ y } \right) } }$

を使うのはダメなんだろうか？というのが頭をよぎります。試してみます。

f:id:ushiostarfish:20190812192509p:plain

確かにもともと存在していた山は改善しています。しかし、端の $x=0$ 付近では誤差が大きいことがわかります。

$\displaystyle { \frac { x-y }{ y } }$

は x が0なら、 $\frac { -y }{ y } =-1$ になるわけで、おっと、この場所はlogp1のバッドランドではありませんか！というわけで結局今回のアプローチでは片方だけではうまくいかないということになります。

ちなみにこの話にはまだ続きがあり、Neil Toronto and Jay McCarthy | Brigham Young University, "Practically Accurate Floating-Point Math" にはさらにロバストにするには、という話もありますし、他の例もありますので、興味のある方は覗いてみると良いと思います。

まとめ

IEEE 754 浮動小数点の内部表現から、誤差の定式化、そしてどのような回避方法があるのかを調べ、そして簡単な改善アプローチを紹介しました。しかし実際に活用する時には、対処したい問題のドメインの知識もさらに加わりますし、６次元、７次元の関数を組み立てるときの対処方法についてもまだあまり良いアプローチを確立できていません。誤差の定式化も、これが最善かというと、確信は持てません。非正規化数のパフォーマンスの低下、フラッシュによる誤差の問題もあります。このように結局銀の弾丸はまだありませんが、浮動小数点がブラックボックスである状況に比べたら余程対策の立てようがあるのではないでしょうか。また、誤差を改善しました、という系の論文を読むときにもこの知識は大いに役に立つのではないかと思います。