hnwの日記

round関数その8:RubyとPythonのround関数は奇妙じゃないんですか?

さて、今回はRubyPythonの題名をつけてしまいましたけど、それらの言語だけに限らないような話題です*1。早速ですが、FreeBSD環境を用意してください(おそらくSPARC Solaris環境やMac OS X環境でも同じ結果だろうと思いますが、僕は確認していません)。

$ ruby -e 'x=0.49999999999999994; printf("%.19f\n%.19f\n", x, x.round());'
0.4999999999999999445
1.0000000000000000000
$ python -c 'x=0.49999999999999994; print "%.19f\n%.19f" % (x, round(x))'
0.4999999999999999445
1.0000000000000000000

どう見ても0.5未満の数を四捨五入して1になっています。環境は以下の通りです。

$ uname -mrs
FreeBSD 6.1-RELEASE-p10 i386  
$ ruby --version
ruby 1.8.5 (2006-08-25) [i386-freebsd6]
$ python -V
Python 2.4.3

これは何度か話題にしているIntel系CPUの浮動小数レジスタが80bitある件と関係していますので、Intel Linuxで本現象が発生する環境は稀だろうと思います。


さて、ではRubyのround関数のソースコードを見てみましょう。下記は1.8.6のものですが、1.8.5でも処理は同じです*2

/*
 *  call-seq:
 *     flt.round   => integer
 *
 *  Rounds <i>flt</i> to the nearest integer. Equivalent to:
 *
 *     def round
 *       return (self+0.5).floor if self > 0.0
 *       return (self-0.5).ceil  if self < 0.0
 *       return 0
 *     end
 *
 *     1.5.round      #=> 2
 *     (-1.5).round   #=> -2
 *
 */

static VALUE
flo_round(num)
    VALUE num;
{
    double f = RFLOAT(num)->value;
    long val;

    if (f > 0.0) f = floor(f+0.5);
    if (f < 0.0) f = ceil(f-0.5);

    if (!FIXABLE(f)) {
        return rb_dbl2big(f);
    }
    val = f;
    return LONG2FIX(val);
}

これがRubyのround関数の中身です。一見すると問題がなさそうなソースコードですが、なぜ上の結果になったのでしょうか。


ここで前回記事「round関数その7:偶数丸め」で話題にした偶数丸めが関係してきます。実はこの0.4999999999999999445というのは0.5-1/2^54、つまりIEEE64bit浮動小数点数で表現できる0.5未満の数の中で最大の数です。この数に上記ソースコードの通りに0.5を足すと1.0-1/2^54になってしまいますが、この数はIEEE64bit浮動小数点数では表現できません。IEEE64bit浮動小数点数で表現できる中で最も近い数は1.0と1.0-1/2^53で、このどちらも同じ距離です。等距離の場合には仮数部の最下位ビットが0となる方に丸めるのがIEEE754の偶数丸めですから、浮動小数点数が内部64bitの計算機では丸めの結果が1.0になり、floor(f+0.5)は1.0になります。


fが正の数の場合について言えば、例えばfloor(f)+0.5とfとを比較してfの方が大きかったらceil(f)を、小さかったらfloor(f)を返す実装にすればこんな問題は起きないわけですけど、fに0.5なんて数を不用意に足したせいで起こった問題だと言えます。


さて、これはバグなのでしょうか。各言語で似た実装がされているところを見ると、今まで多数の人が同様のソースコードをレビューし、もっと多くの人が自分のアプリケーションでこの実装のround関数を使っているはずですけど、誰も気づかなかったのでしょうか*3。それとも、round関数にふさわしい挙動として多くの人が受け入れているのでしょうか。


もしかすると参考になるかもしれませんので、JavaのMath.roundのマニュアルを紹介しておきます。

public static int round(float a)


Returns the closest int to the argument. The result is rounded to an integer by adding 1/2, taking the floor of the result, and casting the result to type int. In other words, the result is equal to the value of the expression:


(int)Math.floor(a + 0.5f)

http://java.sun.com/javase/6/docs/api/java/lang/Math.html#round(double)

これはround関数の動作を完璧に表現している、いいマニュアルですね。0.4999999999999999445を引数にして答えが1になったとしても*4、その挙動が仕様だと一目でわかります。他の言語のマニュアルもこういう書き方がしてあれば悩むことは無いんですけどね。

*1:今回の環境のPHPでも同じ挙動ですし、実はこの環境のlibmも同じ挙動です

*2:コメントが少しだけ変わっています

*3:少なくとも僕は先週まで気づいていませんでしたが…

*4:実際、僕の手元の環境では1になります