hnwの日記

Pythonのround関数で奇数を丸めたら偶数になった

浮動小数点数

Pythonのround関数にバグらしきものを見つけたよ、という報告です。下記は僕のMacBookでの実行結果です。

$ python -c 'x=9007199254740991.0; print "%.19f\n%.19f" % (x, round(x))'
9007199254740991.0000000000000000000
9007199254740992.0000000000000000000


xは9000兆より少し大きい整数で、IEEE754倍精度浮動小数点数で誤差無く表現できる数です。ところが、これをround関数で丸めたら1大きい数になってしまいました。

再現方法

今回も微妙な話題なので、環境によって起きたり起こらなかったりします。私の手元の環境で言うと、MacOSXとFreeBSDで起こり、Linux環境では起こりませんでした。コンパイルオプションによってはLinux環境でも同じ現象が観察できるかもしれません。

原因

整数を丸めて別の整数になるというのは奇妙な気がしますが、この原因はPythonのround関数の実装にあります。Pythonの実装ではこの場合のroundをfloor(x+0.5)で求めるのですが、今回のxについて0.5を足してしまうと9007199254740991.5となり、この数はdoubleでピッタリ表現できません。そこで一番近い数に丸めようとして偶数丸めが起こり、9007199254740992.0になるというわけです。


4503599627370497.0から9007199254740991.0までの奇数は全て同様の条件を満たすので、これらの数はPythonのround関数で丸めると1大きくなります。


実は今回の例は、直前のエントリ「Rubyのround関数の実装が変わってた」でも紹介しているのですが、Rubyの1.8.6-p114以前やPHPの5.2.6以前でも同じ現象が見られます。原因も全く同じで、これらのバージョンでは引数に0.5を足す実装だったためです。

libmのround関数

C99からはround関数が利用できるわけですから、このround関数がどんな挙動なのか調べてみましょう。全環境で中身が同一である保証はありませんが、FreeBSDのroundの実装が参考になりそうです。この実装と、Pythonの実装と、両者を実際に動かしてみることにしましょう。

#include<stdio.h>
#include<math.h>

double bsd_round(double x) {
  double ret;
  if (x >= 0) {
    ret = ceil(x);
    if (ret - x > 0.5) {
      ret -= 1.0;
    }
  } else {
    ret = ceil(-x);
    if (ret + x > 0.5) {
      ret -= 1.0;
    }
    ret = -ret;
  }
  return ret;
}

double python_round(double x) {
  double ret;
  if (x >= 0) {
    ret = floor(x+0.5);
  } else {
    ret = ceil(x-0.5);
  }
  return ret;
}

int main(void) {
  double x = 9007199254740991.0;
  printf("%f\n", bsd_round(x));
  printf("%f\n", python_round(x));
}

$ gcc -O3 round-test.c -lm ; ./a.out
9007199254740991.000000
9007199254740992.000000


FreeBSDの実装では奇数は奇数のままですが、Pythonの実装では偶数になることがわかりました。


両者の実装は計算順序が変わっただけで同一に見えますが、浮動小数点数の演算では計算順序が計算結果に影響することもあるわけです。


xが正の場合について言うと、FreeBSDの実装ではceil(x)-xが0.5より大きいかどうかでどちらに丸めるかを決定します。この演算の結果は必ずdoubleでピッタリの数になるため、Pythonの実装のように計算途中で丸めが入ることがありません。

これがバグかどうかについて

まあバグですよね。Pythonのリファレンスマニュアルによれば「rounded to the closest multiple of 10 to the power minus n;」だそうですから、第二引数が省略された場合には最も近い整数に丸められるはずです。


実装のバグなのか、ドキュメントのバグなのかは僕にはわかりません。別にどっちでもいいと思います。

おまけ：Javaのround関数

Javaのマニュアルが素晴らしいというのは以前も書きましたが、再度紹介します。


今回の不思議な挙動は、Javaのround関数でも起こると思います。でも、JavaのAPIリファレンスに下記の通り書いてありますから、こんな挙動を見つけても仕様であることが一目でわかります。

public static long round(double a)


Returns the closest long 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 long. In other words, the result is equal to the value of the expression:


(long)Math.floor(a + 0.5d)

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


やはりマニュアルは重要ですよね。仕様の善し悪しよりドキュメントの整備の方が断然重要です。