VC6における浮動少数点数値について

質問

はじめて質問させて頂きます。

C言語の基礎演習レベルの題材ですが
想定外の結果が表示され解決できておりません。
内容としては2つのfloat型変数(以下ソースのb,c)に
全く同じ演算処理結果を代入し
コンソール表示させるというプログラムです。

#include <stdio.h>

void main()
{
float a = 0, k = 0;
float b = 0;
float c = 0;

a = 0.2;
k = 20.0;

k += a;
b = k - 20.0;
c = k - 20.0;

printf( "%f %f %f %f\n", a, b, c, k );

if( b == c )
{
printf( "*\n" );
}

while( getchar() != EOF ){}
}

予想実行結果：
0.200000 0.200000 0.200000 20.200000
*

実際の結果：
0.200000 0.200000 0.200001 20.200001

当方、浮動少数点数値の扱いにも不慣れな為、原因不明です。
なぜこのような結果になるのか、どなたかご教授頂けないでしょうか。

宜しくお願い致します。

回答No.10

Tacosan

ついでにいうと「20.0が２進化浮動小数点数の循環小数だった事」は意味不明＞#7. 「2進化浮動小数点数の循環小数」ってなんだろう?
まあ, いずれにしてもこの辺の細か～い誤差はあったりするので浮動小数同士を == や != で比較するのは割と危険.

回答No.9

nda23

そもそも、浮動小数点型を10進数にした時の信頼精度をご存知ですか？
単精度では6桁、倍精度では15桁です。これは少数部だけでなく、整数
でも発生する問題です。入口/出口だけではなく、途中の一時的な値
についても、同桁数以内になっていないと、誤差が拡大します。
それらの事情を考慮したうえで、ワークも含め型を決定してください。

演算誤差を抑えるため、精度の高い変数に代入してから演算する等の
方法を明示的に採用する方が分かり易いと思いますし、周囲の納得を
得られるのではないでしょうか。

尚、最適化を行ったため、不具合が生じる例は実数演算以外にも幾つか
見受けられます。

回答No.8

titokani

＞簡単に言えば、上記のプログラムは
＞b = (float)((double)20.2 - (double)20.0);
＞c = (float)((double)(float)20.2 - (double)20.0);
＞というプログラムと等価。
いやいや、いかに最適化といえど、勝手に型を変えてはいけません。
k+=aの型は、あくまでもfloatです。これは規格できっちり定められています。
今回は最適化によって、演算誤差が異なってしまうことの一例とだと思います。

回答No.7

chie65536（@chie65535）

ちょっと訂正と追加。

＞簡単に言えば、上記のプログラムは
＞b = (float)((double)20.2 - (double)0.2);
＞c = (float)((double)(float)20.2 - (double)0.2);
＞というプログラムと等価。

は

簡単に言えば、上記のプログラムは
b = (float)((double)20.2 - (double)20.0);
c = (float)((double)(float)20.2 - (double)20.0);
というプログラムと等価。

の間違い。

不運だったのは「20.0が２進化浮動小数点数の循環小数だった事」なのかも知れない。

回答No.6

chie65536（@chie65535）

最適化により
k += a;
b = k - 20.0;
c = k - 20.0;
は
b = (k += a) - 20.0;
c = k - 20.0;
となっています。

そして、暗黙的な型キャストにより、以下のようになります。

//部分式「(k += a) - 20.0」はすべてdouble型で計算されている
//floatへの変換により精度が落ちるのは、bへの代入とkへの代入のみ
//(k += a)の値そのものはdoubleの精度を保持している事に注意
b = (float)((k += a) - 20.0);
//kを取り出すと、floatからdoubleへの変換が起きる
//従って、この時点での部分式「(k)」は「(double)20.2」と等しくない
c = (double)(float)(k) - 20.0;

簡単に言えば、上記のプログラムは
b = (float)((double)20.2 - (double)0.2);
c = (float)((double)(float)20.2 - (double)0.2);
というプログラムと等価。

「誰がどう見ても、この２つの式は異なる」ので「結果が異なる」のは当たり前。

なお「こういう現象が起きた原因」は「最適化」にあり、言語仕様上「最適化は環境に依存し、コンパイラごとに異なる」ので、VC6のバグとは言えない。

コンパイル時に最適化を抑止しなかった事、最適化による実行結果の差異を意識しない書き方をした事を考えれば、こう書いた人間に責任があると言える。

双方で同じ結果を得たいのであれば
k += a;
c = 0.0;//上のkへの代入と下のkの参照をぶった切って最適化させない
b = k - 20.0;
c = k - 20.0;
と書くか
k += a;
b = (double)(float)(k) - 20.0;//明示的にfloatに精度を落とさせる
c = (double)(float)(k) - 20.0;//明示的にfloatに精度を落とさせる
と書くか
k += a;//これも最適化されて意味無くなるかも知れないが
b = (0.0,k) - 20.0;//カンマ演算子で「捨て演算」をして最適化させない
c = k - 20.0;
と書きましょう。

回答No.4

titokani

VC++2005で試したところ、
debug
0.200000 0.200001 0.200001 20.200001
*
release
0.200000 0.200001 0.200001 20.200001
*
VC6だと、
debug
0.200000 0.200000 0.200001 20.200001
release
0.200000 0.200000 0.200000 20.200000
*
でした。
VC6debugのアセンブリコードは、
8: float a = 0, k = 0;
00401028 mov dword ptr [ebp-4],0
0040102F mov dword ptr [ebp-8],0
9: float b = 0;
00401036 mov dword ptr [ebp-0Ch],0
10: float c = 0;
0040103D mov dword ptr [ebp-10h],0
11:
12: a = 0.2;
00401044 mov dword ptr [ebp-4],3E4CCCCDh
13: k = 20.0;
0040104B mov dword ptr [ebp-8],41A00000h
14:
15: k += a;
00401052 fld dword ptr [ebp-8]
00401055 fadd dword ptr [ebp-4]
00401058 fst dword ptr [ebp-8]
16: b = k - 20.0;
0040105B fsub qword ptr [__real@8@4003a000000000000000 (00426030)]
00401061 fstp dword ptr [ebp-0Ch]
17: c = k - 20.0;
00401064 fld dword ptr [ebp-8]
00401067 fsub qword ptr [__real@8@4003a000000000000000 (00426030)]
0040106D fstp dword ptr [ebp-10h]
18:
19: printf( "%f %f %f %f\n", a, b, c, k );
でした。
まさしく、＃２さんのおっしゃる通りですが、VC6のバグといっていいかもしれませんね。

お礼

回答ありがとうございます。

＃２さんへのお礼で記載させて頂いたソースの修正行った結果からも
納得のいく回答を頂きました。

ありがとうございます。

回答No.3

nda23

単精度実数では符号1、指数8、仮数24ビットで表現します。
コンピュータの中身は２進法なので、0.2の表現が難しい。
先ず、2の-3乗である0.125が入り、2の-4乗の0.0625が入り、・・・
という具合です。２進で示すと以下の通りです。
0-01111100-(1)10011001100110011001101
左端は符号で、正なので0です。次の8ビットは0x7Cですが、0x7Fが
0なので、0x7Cは-3を示します。(1)は実データ上は存在しませんが、
有効ビットの最上位を仮定するので、常に１があるものとして操作
します。その右は2の-4乗となりますね。ここに20を足すと、最上位は
2の4乗になるので、指数は0x83(指数の0点は0x7F)になります。
0-10000011-(1)01000011001100110011010 →1001101が失われる。
整数部を優先するため、最初の状態から右の7ビットが消えます。
失われたビットの最上位が1なので、切り上げられ、最後が1001では
なく1010になります。ここから20を引くと、上位5ビットが無くなり、
有効ビットが最上位になるようシフトされます。
0-10000011-(0)00000011001100110011010 ←上位5ビットが消える。
0-01111100-(1)10011001100110100000000 ←左シフトする。(0補充)
0-01111100-(1)10011001100110011001101 ←最初の0.2の時の内容
　　　　　　　　　　　　　　　　　 ↑
2の-18乗のところが切り上がっているのが分かりますか？
これが10進に変換した時の「誤差」になるのです。
因みに2の-1乗は0.5なので、右にシフトしてこぼれたビットが1なら
切り上げます。つまり、四捨五入というわけです。
ご理解いただけましたか？

回答No.2

Tacosan

そこではなく「同じ計算をしているはずの b と c でなぜ値が異なるのか」を問題にしているのだと思います＞#1.
「ひょっとして」ですが, ひょっとすると「b を計算するときには直前の計算で求めた k の値をそのまま使い, c を計算するときには k の値をメモリから読みだしている」可能性があります.
つまり, b を計算するときには直前の計算で求めた (double としての) k の値を使い, c を計算するときにはメモリに格納した (float としての) k の値を使っているという可能性です.
手元の VS2008 でこの状況をシミュレートすると, 確かに言われるような結果になります. また, printf のフォーマットをすべて %.16f にすると
0.2000000029802322 0.2000000029802322 0.2000007629394531 20.2000007629394530
という結果が得られています.

お礼

回答ありがとうございます。

titokaniさんも同意見のようです。
実際にソースに変更加えてみたらよく理解できました。
以下のような変更加えたらbとcは同じ結果が得られました。
大変参考になります。

k += a;
c = 0;　←　追加
b = k - 20.0;
c = k - 20.0;

結果：
0.200000 0.200001 0.200001 20.200001

回答No.1

SilverThaw

興味深い内容ですね。
VC6が手元にあるので簡単に確認してみました。

以下に要素のある個所からのメモリダンプ(一部)を付記します
>a = 0.2;
　a = cd,cc,4c,3e,
>k = 20.0;
　k = 00,00,a0,41,
>k += a;
　k = 9a,99,a1,41,
>b = k - 20.0;
　b = 00,cc,4c,3e,
a(0.2)とk(20.00)では指数部(4Byte目)の内容が違っています。
その為、計算前に同じ指数の条件変更した際に発生する誤差のようですね。

お礼

回答ありがとうございます。
丸め誤差のようなものだと理解させて頂きました。
大変参考になります。