bcc32 double と unsigned 比較で

[Hermit]

お久しぶりです。
bcc32 で、以下のプログラムの動作がどうもよくわからないので、教えてください。

sqrtの精度はどのくらいなのかなと思って以下のようなプログラムでテストをしてみたのですが、

コード:

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

int main() {
	unsigned i = 1;
	union {
		double d;
		unsigned u[2];
	} d;
	for (;i<94906268;i++) {
		if (i != (d.d = sqrt((double)i*i))) {
			printf("%u err\n",i);
			break;
		}
	}
	printf("%u %f %x%x %d\n",i,d.d,d.u[0],d.u[1],(unsigned)d.d);
	return 0;
}

d:\usr>bcc32 foo.c
Borland C++ 5.5.1 for Win32 Copyright (c) 1993, 2000 Borland
foo.c:
Turbo Incremental Link 5.00 Copyright (c) 1997, 2000 Borland

d:\usr>foo
94906267 err
94906267 94906267.000000 6c0000004196a09e 94906267

の様になってしまいました。

いまひとつ計算が合わない理由がわからないのですが、
どうしてでしょう？

[beatle]

パッと見原因がわかりませんね．
gccでコンパイル＆実行したところ，

$ ./a.exe
94906268 94906267.000000 6c0000004196a09e 94906267

となりました．

[Hermit]

lcc_win32 では大丈夫でした。－１まで回ってくれます。
bcc32 は、なぜかエラーの部分を通ってしまいます。

ちょっと、コードが違ってたみたいです（やってるうちに頭がこんがらがってきて・・・）

コード:

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

int main() {
	unsigned i = 1;
	union {
		double d;
		unsigned u[2];
	} d,d1;
	for (;i<94906268;i++) {
		if (i != (d.d = sqrt((double)i*i))) {
			printf("%u err\n",i);
			break;
		}
	}
	d1.d=i;
	printf("%u %f %x%x %x%x\n",i,d.d,d.u[0],d.u[1],d1.u[0],d1.u[1]);
	return 0;
}

で、

d:\usr>foo
94906267 err
94906267 94906267.000000 6c0000004196a09e 6c0000004196a09e
の様になりました。うちでは。

[box]

質問者さんは、どういう結果を得ることを想定されているかを示す必要があるような気がします。

[Hermit]

はい、では
if (i != (d.d = sqrt((double)i*i)))
double のビット列が同じ時、同じと判断して、以後の
printf("%u err\n",i);
を通らないのではないかと思っているのです。

if (i != (unsigned)(d.d = sqrt((double)i*i)))
でも、
printf("%u err\n",i);
の方を通ってしまうので、
原因が知りたいのです。

[Hermit]

boxさん＞バグのないプログラムはない。

ｂｃｃ３２のバグってことでしょうか？
それならそれでいいんですが。

[beatle]

「バグのないプログラムはない。」はboxさんがHermitさんに書いたメッセージではありませんので，そこに返信するのは的外れだと思います．
まあ，bcc32のバグの可能性もありますけれど．

さて，僕もどうやったら問題を切り分けられるか悩んでいるのですが，まずif文の中で代入するのをやめて見ませんか．
それから，if文の中では，iがdoubleにキャストされてから!=によって比較されているはずですので，一度double型変数にiを代入してから，そのdouble型変数とd.dを比較してみてはいかがでしょうか．

[box]

ｂｃｃ３２のバグってことでしょうか？

あれは、メールの最後に付ける署名みたいなものです。
そこに反応されても、ちょっと…。

[Hermit]

＞あれは、メールの最後に付ける署名みたいなものです。
＞そこに反応されても、ちょっと…。

そういうのがあったんですね。なんか追伸で書かれたのかと思ってました。失礼しました。

で、とりあえずコードを以下のように変えてみました。

コード:

#include<stdio.h>
#include<math.h>
 
int main() {
    unsigned i = 1;
    union {
        double d;
        unsigned u[2];
    } d,d1;
    for (;i<94906268;i++) {
	d.d =sqrt((double)i*i); d1.d = i; 
        if (d1.d != d.d) {
            printf("%u err\n",i);
            break;
        }
        if (d1.d != sqrt((double)i*i)) {
            printf("%u sqrterr\n",i);
            break;
        }
    }
    printf("%u %f %x%x %x%x\n",i,d.d,d.u[0],d.u[1],d1.u[0],d1.u[1]);
    return 0;
}

結果、
c:\usr>foo
94906267 sqrterr
94906267 94906267.000000 6c0000004196a09e 6c0000004196a09e

if の比較は代入した d.d では同じになりましたが、sqrt() の返却値だと差異が有る様です。
とりあえず -S でアセンブラにしてみたけど・・・何をやってるのか良くわからなかった・・・
それ以前はいいのだから、コード自体は問題ないですよね。比較も qword みたいだし。

コード:

   ;	        if (d1.d != d.d) {
	fld       qword ptr [esi]
	fcomp     qword ptr [ebp-8]
	fnstsw ax
	sahf
	je        short @4
--------------------------------------
   ;	        if (d1.d != sqrt((double)i*i)) {
@4:前半省略
	call      _sqrt
	add       esp,8
	fcomp     qword ptr [esi]
	fnstsw ax
	sahf
	je        short @6

[Hermit]

ちょっとだけわかってきました。
単純に、sqrt() の返却値を そのまま使うと、誤差が出るようです。

コード:

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

int main() {
	unsigned i,k,l;
	double d,e;

	i = 94906267;
	e = (double)i*i;
	d = sqrt(e);
	k = (unsigned)sqrt(e);
	l = (unsigned)d;

	printf(	"i->%u\n"
		"e->%f\n"
		"d->%f\n"
		"k->%u\n"
		"l->%u\n"
		,i,e,d,k,l);
	return 0;
}

で、結果
C:\usr>bar
i->94906267
e->9007199515875288.000000
d->94906267.000000
k->94906266
l->94906267

で、直接キャストすると誤差が出てるようです。

で、アセンブラコードにすると、

コード:

   ;	
   ;		d = sqrt(e);
   ;	
	push      dword ptr [ebp-12]
	push      dword ptr [ebp-16]
	call      _sqrt
	add       esp,8
	fstp      qword ptr [ebp-8]
   ;	
   ;		k = (unsigned)sqrt(e);
   ;	
	push      dword ptr [ebp-12]
	push      dword ptr [ebp-16]
	call      _sqrt
	add       esp,8
	call      __ftol
	mov       esi,eax

の様になって、キャストすると、直接 __ftol しているから、一旦 double にする為に fstp している値と違うと見ました。
アセンブラは良くわからないんですが、sqrt の返却値は、直接 fpu命令の返却値を使い、
10byte で帰ってくる様に思うのですが、どうでしょう？

[かずま]

i->94906267
e->9007199515875288.000000

i * i は 9007199515875289 です。
これは、16進で 0x2000000f9097d9 ですから 54ビットの精度が必要です。
double の精度は 53ビットなので、sqrt() の引数として渡すとき
最下位のビットが落ちて 9007199515875288.0 となります。 なので、
sqrt() の結果は 94906267.0 にならず、94906266.99999999473... になる
のですが、これは内部の浮動小数点レジスタに持っている値ですから、
double の変数にストアすると、54ビット目を丸めて 53ビットにすると
94906267.0 になるのです。
sqrt() の結果を double の変数にストアせず、他の double の変数と
比較するとき、BCC は浮動小数点レジスタの値を丸めずに比較しています。

コード:

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

int main(void)
{
    int i;
    for (i = 94906260; i < 94906275; i++) {
        __int64 j = (__int64)i * i;
        long double e = sqrt((double)i * i);
        printf("%d %I64d %I64x %.15Lf\n", i, j, j, e);
    }
    return 0;
}

実行結果

コード:

94906260 9007198187187600 1fffffc05e6d90 94906260.000000000000000
94906261 9007198377000121 1fffffcbaebcb9 94906261.000000000000000
94906262 9007198566812644 1fffffd6ff0be4 94906262.000000000000000
94906263 9007198756625169 1fffffe24f5b11 94906263.000000000000000
94906264 9007198946437696 1fffffed9faa40 94906264.000000000000000
94906265 9007199136250225 1ffffff8eff971 94906265.000000000000000
94906266 9007199326062756 200000044048a4 94906266.000000000000000
94906267 9007199515875289 2000000f9097d9 94906266.999999994730000
94906268 9007199705687824 2000001ae0e710 94906268.000000000000000
94906269 9007199895500361 20000026313649 94906268.999999994730000
94906270 9007200085312900 20000031818584 94906270.000000000000000
94906271 9007200275125441 2000003cd1d4c1 94906270.999999994730000
94906272 9007200464937984 20000048222400 94906272.000000000000000
94906273 9007200654750529 20000053727341 94906272.999999994730000
94906274 9007200844563076 2000005ec2c284 94906274.000000000000000

i = 94906266
i * i = 9007199326062756 = 0x200000044048a4
これは下位ビットが 0100 なので、精度は 52ビットであり正しい結果が出ます。
94906268, 94906270, 94906272 などの偶数も同様です。

[かずま]

コード:

    printf("%u %f %x%x %x%x\n",i,d.d,d.u[0],d.u[1],d1.u[0],d1.u[1]);

94906267 94906267.000000 6c0000004196a09e 6c0000004196a09e

94906267.0 は、0x6c0000004196a09e ではなく、0x4196a09e6c000000 です。
x86 はリトルエンディアンですから。

コード:

#include <stdio.h>

int main(void)
{
	union {
		double d;
		unsigned u[2];
		__int64 x;
	} d;
	d.d = 94906267;
	printf("%.1f\n", d.d);
	printf("%08x%08x\n", d.u[1], d.u[0]);
	printf("%016I64x\n", d.x);
	return 0;
}

[Hermit]

>比較するとき、BCC は浮動小数点レジスタの値を丸めずに比較しています。

やはり、そういういことでしたか。理解できました。

ということは、bcc32 では、

コード:

#include <stdio.h>
#include <math.h>
int main() {
	double s2=sqrt(2);
	if (s2 == sqrt(2))
		printf("true");
	else
		printf("false");
	return 0;
}

の様なコードは、false になるってことですね。
多分、何かのオプションで回避できるんでしょうね。

＞54ビット目を丸めて 53ビットにすると
丸めが発生していたからですか。
オーバーフローは後から考えようと思ってたので、気にしていませんでした。
後回しにしていた事が早々とわかりました。ありがとうございます。

＞x86 はリトルエンディアンですから。
そうでしたね。忘れてました。

疑問が解けてすっきりしました。ありがとうございます。＞ かずまさん