hnwの日記

ぼくのかんがえたさいきょうのround関数

浮動小数点数の丸めにおいて丸め桁数を指定でき、それでいて精度を失わないようなround関数をCで実装してみました。

実装としては、受け取った浮動小数点数から最短になる10進表記に変換し、浮動小数点をズラすことなく10進表記のまま四捨五入を行うものです。これを元に偶数丸めを実装するのも容易でしょう。

実際、前回記事「RubyとPythonとC#のround関数のバグっぽい挙動について」で指摘した5.015の例についても期待通りに丸めることができます。

#include <stdio.h>

extern double precise_round(double x, int digits);

int main()
{
  printf("%f\n", precise_round(5.015, 2)); // 5.02
  printf("%f\n", precise_round(5.0149999999999987921, 2)); // 5.01
}

浮動小数点数とその10進表記の正確な相互変換は1990年の論文で決着済み

浮動小数点数とその10進表記との間で最も近い数に変換する方法は「How to Read Floating Point Numbers Accurately *1」「How to Print Floating-Point Numbers Accurately *2」という論文で示されています。また、この論文を元にしてDavid M. Gayが実装したdtoa.cは、多くのオープンソースプロジェクトで利用されています。

前回記事で紹介したMySQLの処理でも浮動小数点数から最短の10進表記を作り出していますが、この処理には上記dtoa.cの改変版が利用されています。ちなみに僕の今回の実装はオリジナルのdtoa.cをそのまま利用しました。

浮動小数点数リテラル浮動小数点数に直す処理でも似た議論があった

実は以前の記事「PHP以外全員不正解」「Rubyの浮動小数点数リテラルの扱いは正しいのか」「MySQLの自前strtod実装がタコすぎる」においても似たような話題がありました。これらの記事で、一部の言語やミドルウェアでは浮動小数リテラルを最近接の浮動小数点数に変換できていない、という指摘をしました。

これらはどれも自前で小数点をずらす処理を行っていたり10.0倍を繰り返したりと誤差が蓄積しかねないナイーブな実装になっていました。最新版では修正されて上記dtoa.cの成果を利用しているものがほとんどだと思います。

僕は、今回の各言語のround関数の実装も同様にナイーブな実装だと感じました。また、ナイーブな実装が生み出す誤差を避けられない誤差として受け止めている人が一定数いる気もしており、その意味でも非常に似ていると思います。

今回のround関数を各言語が実装すべきかどうかについて

結論から言うと、各言語のround関数が今回作ったような正確な処理になっている必要はないと思います。というのも、round関数は精度を落とす演算なので、ここで多少の誤差が混入したところで致命的な誤差の蓄積にはつながらないと考えられるからです。また、浮動小数点数リテラル浮動小数点数にする処理と異なり、何かの仕様で正確さについての要請があるわけでもありません。

つまり、僕個人としては、マニュアルに実装が明記してあればどんな実装でもいいんじゃないかと思っています。万一困る人がいたら自前実装したり、10進型を採用すればいいわけですし。

とはいえ、無用な誤差を生みがちなので、浮動小数点をずらして戻すような処理を見たら怖いと感じることは重要だと思います。また、言語実装をするような人も同じ直感を持っているはずだと思うので、各言語で不正確な実装になっているのを意外に感じたというのが前回記事を書いたきっかけでした。

*1:William D Clinger, "How to Read Floating Point Numbers Accurately", ACM SIGPLAN '90, pp.92--101, 1990. (同じ論文のretrospective版のPDF版

*2:Guy L. Steele, Jr. and Jon L. White, "How to Print Floating-Point Numbers Accurately", ACM SIGPLAN '90, pp.112--126, 1990. (同じ論文のretrospective版のPDF版