[Java,C] 浮動小数点での丸め誤差の例

floatやdoubleといった浮動小数点数で計算をしていると、丸め処理をしているつもりでも誤差が出ることがある。
これの簡単な例がないかなぁと思っていたんだけど、 0.35 * 3.0 でうまく表現できそう。

Javaで 0.35 * 0.3 を計算してみる

Javaのdoubleで 0.35 * 0.3 を計算してみると

System.out.println(0.35 * 3.0);
// 1.0499999999999998

1.05 になりそうなところが 1.0499999999999998 になっている。
これだと、計算結果を小数点２桁目で四捨五入なんてことをやろうとしても、 1.1 ではなく、1.0 になってしまう。

仮数部を取り出して、2進数表示してみると、1.05 とは 1 違うことが分かる。

System.out.println(
    Long.toBinaryString(
    0x000fffffffffffffL & Double.doubleToLongBits(0.35 * 3.0)));
System.out.println(
    Long.toBinaryString(
    0x000fffffffffffffL & Double.doubleToLongBits(1.05)));

// 0.35 * 3.0
// 110011001100110011001100110011001100110011001100
// 1.05
// 110011001100110011001100110011001100110011001101

なんで、こんなことになっちゃったんだろうか。

まず、0.35 の仮数部の2進表現を見てみよう。

System.out.println(
    Long.toBinaryString(
    0x000fffffffffffffL & Double.doubleToLongBits(0.35));

// 110011001100110011001100110011001100110011001100110

仮数部の正規化表現だと、先頭は1で始まるけど、先頭ビットは明示的には表現されていないので、53ビット目に1をつけてあげると以下のようになる。

10110011001100110011001100110011001100110011001100110

これを long で2倍、3倍してみると、以下のようになる。

long n = 0b10110011001100110011001100110011001100110011001100110L;
System.out.println(Long.toBinaryString(n));
System.out.println(Long.toBinaryString(2 * n));
System.out.println(Long.toBinaryString(3 * n));

// n
// 10110011001100110011001100110011001100110011001100110
// 2 * n
// 101100110011001100110011001100110011001100110011001100
// 3 * n
// 1000011001100110011001100110011001100110011001100110010

それぞれ、仮数部に収まる範囲を超えてしまっているので、端数を丸めることになるけど、この時の丸め方としては再近接偶数丸めとなる。
これは丸めた結果が丸める前の値に近い方に丸めるんだけど、丸めた結果が丸める前とちょうど等距離の場合は丸めた結果が偶数になるように丸めるやり方になる。

IEEE 754 最近接丸め
The Java® Language Specification 4.2.4. Floating-Point Operations

なので、3倍した値の端数の 10 が切り上げ方向、切り捨て方向のどちらに丸められるかは丸められる桁の1つ上の桁の値に依存する。
今回の場合は該当桁を見ると 0 になっているので、偶数にするためには切り捨て方向に丸めることになる。

最終的に丸めた結果は以下のようになる。

10000110011001100110011001100110011001100110011001100

これの先頭の 1 を除くと
110011001100110011001100110011001100110011001100

これで元の 0.35 * 3.0 の仮数部と一致した！

Cで丸め方向を指定してみる

Javaではfloat,doubleの計算は必ず再近接偶数丸めで計算されるけど、Cでは丸め方向を指定することができる。
ただし、これはプロセッサに依存するところが大きいので、どんな処理系でもできるとは限らない。

fenv.c

#include <fenv.h>
#include <stdio.h>

int main(void) {
    double dx = 0.35;
    double dy = 3.0;

    fesetround(FE_TONEAREST);
    double d1 = dx * dy;

    fesetround(FE_UPWARD);
    double d2 = dx * dy;

    fesetround(FE_TONEAREST);
    printf("%d\n", (int)(d1 * 100.0) % 10);
    printf("%d\n", (int)(d2 * 100.0) % 10);

    printf("%d\n", d1 == d2);
    printf("%d\n", d1 == 1.05);
    printf("%d\n", 1.05 == d2);

    return 0;
}

$ ./fenv
4
5
0
0
1
$

丸め方向は fesetround で設定する。
FE_TONEAREST で 0.35 * 3.0 を計算した結果はJava同様に 1.04 になってるけど、FE_UPWARD で切り上げ方向に丸めた場合は 1.05 になっている。

これを確認するために 100.0をかけて、intにキャストしてから、10 の剰余を取っている。