先日、 ftoa (少数を文字列に変換する関数) を自作していたときに、不可解な現象に出くわしました。問題となったのは下記のようなプログラムです。


#include &ltstdio.h>

void main()
{
double f = 0.12;
int d;
f *= 10; // f = 1.2
d = (int)f; // d = 1
f -= d; // f = 0.2
f *= 10; // f = 2.0
d = (int)f; // d should be 2
f -= d; // f should be 0.0
printf("%lf\n", f);
}



この結果は当然ゼロとなるべきところですが。結果は、1.000000 となってしまいます。極論すると、

2.0 - 2 = 1.0

と同じということです。しかも、問題が発生するのは、0.12 という数字だけ。0.11 でも 0.13 でも 0.23 でも 0.34 でも問題は出ません。

あまりにも不可解なので、こんな問題に答えられるのは、世界中のプログラマが集まる掲示板 Stack Overflow くらいだろうということで、思い切って質問することにしました。ドキドキ。



http://stackoverflow.com/


日本の掲示板では答えに巡り合うまでに時間がかかることが多いですが、さすが世界の天才プログラマが跋扈する Stack Overflow、なんと数分後には回答が来ました。すげぇ・・・。アドバイスに従い、下記のプログラムを動かすと答えはすぐに見つかりました。



#include &ltstdio.h>

void main()
{
double f = 0.12;
int d;
printf("f :f: %.16f\n", f);

f *= 10;
printf("f *= 10 :f: %.16f\n", f);

d = (int)f;
printf("d = (int)f :d: %d\n", d);

f -= d;
printf("f -= d :f: %.16f\n", f);

f *= 10;
printf("f *= 10 :f: %.16f\n", f);

d = (int)f;
printf("d = (int)f :d: %d\n", d);

f -= d;
printf("f -= d :f: %.16f\n", f);
}



演算の結果はこちらです。



1 f :f: 0.1200000000000000
2 f *= 10 :f: 1.2000000000000000
3 d = (int)f:d: 1
4 f -= d :f: 0.2000000000000000
5 f *= 10 :f: 1.9999999999999996
6 d = (int)f:d: 1
7 f -= d :f: 0.9999999999999996



4行目から5行目でいきなり数値が豹変しています。ここで問題が発生していました。問題発生のメカニズムも教えてもらいました。

1.0 は、C言語の採用している binary64 では正しくは以下の数値となります。17桁以降は丸め誤差です。

1.0000000000000000444089...

上記の値を4行目の 1.2 -1.0 に適用すると下記のようになります。

1.2000000000000000 - 1.0000000000000000444089 = 0.1999999999999999555911

この段階では、小数点以下17桁以下の数は丸め込まれるので、0.2 と表示されます。

しかし、5行目で10をかけると、状況が一変し、この誤差がひょっこり顔をのぞかせます。末尾の一桁が、繰り上げられますので、f の値は、1.9999999999999996 となってしまいます。

あとはプログラムの通り、この微妙な誤差が致命的な違いとなって表れ、期待された値とまったく異なる結果となります。しかもこの条件にあてはまるのは、17桁以降の数値と丸め誤差が重なったときなので非常に稀です。

このようなプログラムが人命にかかわる車や飛行機、ロケットなどに混入していたらと思うとぞっとします。プログラマの皆さん、浮動小数点の扱いは十分に注意してくださいね。

せっかくなので、Stack Overflow で解決策ないの?と質問をしたら下記の対策を勧められました。



#include &ltstdio.h>
#include &ltfloat.h>

void main()
{
double f = 0.12;
int d;

printf("f :f: %.16f\n", f);

f *= 10;
printf("f *= 10 :f: %.16f\n", f);

d = (int)f;
printf("d = (int)f :d: %d\n", d);

f -= d;
printf("f -= d :f: %.16f\n", f);

f *= 10;
printf("f *= 10 :f: %.16f\n", f);

f *= (1 + DBL_EPSILON);
printf("DBL_EPSILON:f: %.16f\n", f);

d = (int)f;
printf("d = (int)f :d: %d\n", d);

f -= d;
printf("f -= d :f: %.16f\n", f);
}



EPSILON は、理系の人ならお馴染みのイプシロン・デルタのイプシロンです。float の場合は接頭子は FLTになりますので注意してください。詳細については下記のサイトが参考になります。


Is floating point math broken?
http://stackoverflow.com/questions/588004/is-floating-point-math-broken


普段、見てばかりいる Stack Overflow ですが、今回は自分で参加してみて、あらためて世界の叡智の凄さを体感しました。日本のプログラマの皆さんもどんどん世界に出ていきましょう!
(^_^)/





苦しんで覚えるC言語

  • 作者: MMGames
  • 出版社/メーカー: 秀和システム
  • 発売日: 2011/06/24
  • メディア: 単行本






絵で見てわかるシステムパフォーマンスの仕組み

  • 作者: 小田 圭二
  • 出版社/メーカー: 翔泳社
  • 発売日: 2014/06/21
  • メディア: 単行本(ソフトカバー)