double型同士の比較演算・内部表現

Question

前提・実現したいこと

double型の内部表現を明らかにし、誤差に対する理解を深める目的でdouble型をビット単位で表示する関数を
ウェブサイトを参考にしながら作りました。(https://www.k-cube.co.jp/wakaba/server/floating_point.html)
double型の10.1, 10.2, 20.3をそれぞれビット表示し、誤差が原因で10.1 + 10.2 != 20.3であることを確認しました。
そして、計算機イプシロンを左辺に加えることで両辺のビット表示が同じになることを確認し、それらを比較演算したところ、予想に反し、両辺の値が違うという結果が出ました。
どのようなことが起こっているのか教えていただけるとありがたいです。

該当のソースコード

#include <stdio.h>
#include <float.h>

void dump(unsigned char *p);

void main() {

    double x, y, z;
    x = 10.1;
    y = 10.2;
    z = 20.3;
	union {
		double d;
		unsigned char c[8];
	} uni;

    uni.d = x + y;
    printf("      x + y: ");
    dump(uni.c);

    uni.d = DBL_EPSILON + x + y;
    printf("eps + x + y: ");
    dump(uni.c);

    uni.d = z;
    printf("          z: ");
    dump(uni.c);

    if (DBL_EPSILON + x + y != z)
        puts("eps + x + y != z");
}        
void dump(unsigned char *p) {
	int i, j;
	
    //p[7]から下ってアクセスしていく
	for (i = 7; i >= 0; i--) {

        //1バイトを2進数表示する
        for (j = 7; j >= 0; j--) {
            printf("%d", (p[i] >> j) % 2);
        }
        putchar(' ');
    }
	putchar('\n');
}

実行結果

      x + y: 01000000 00110100 01001100 11001100 11001100 11001100 11001100 11001100
eps + x + y: 01000000 00110100 01001100 11001100 11001100 11001100 11001100 11001101
          z: 01000000 00110100 01001100 11001100 11001100 11001100 11001100 11001101
eps + x + y != z

Answer 1

Windows10でやってみましたが、コンパイラによって違いますね。

gcc version 10.1.0 (Rev3, Built by MSYS2 project)
      x + y: 01000000 00110100 01001100 11001100 11001100 11001100 11001100 1100 1100
eps + x + y: 01000000 00110100 01001100 11001100 11001100 11001100 11001100 1100 1100
          z: 01000000 00110100 01001100 11001100 11001100 11001100 11001100 1100 1101
eps + x + y != z

で、見たとおりです。
Borland C++ 5.5.1 for Win32 Copyright (c) 1993, 2000 Borlandだと、質問文と同じ結果になります。

x86系の浮動小数点レジスタは80bitで仮数部64bitなので、doubleへの格納で丸めが発生しますが、その丸めの仕方がコンパイラによって違うのでしょう。
Cでの!=は、浮動小数点同士を比べているということでしょう。

Comments

[otn]

DBL_EPSILON + x + y を一旦変数に入れて、その変数と z を比べれば 等しくなると思います。

[dotbot]

なるほど...計算した時点では仮数部が64bitもあるんですね...勉強になります。

Answer 2

じつは、CPU内部の浮動小数点演算ユニットの演算幅は６４ビットではなかったりします
doubleとして表現する場合は、６４ビットに丸められますが、内部演算では、６４ビットとして行われるとは限りません。

ということで、この演算が真となるか偽となるかは、CPUによって変わったりなんかします。

結論、浮動小数点数の==あるいは!=の比較はしてはいけません

インテル(R) アーキテクチャ (IA) 浮動小数点ユニット (FPU)、ストリーミング SIMD 拡張命令 (SSE)、ストリーミング SIMD 拡張命令2 (SSE2) を使用した浮動小数点算術演算 - w_fp_precision_j.pdf

Comments

[Zuishin]

double の精度が CPU によって変わるんですか？

[y_waiwai]

doubleじゃなく、FPUの搭載有無で内部演算結果が変わります

[Zuishin]

私の PC では eps + x + y の表示が質問とは異なるものになりました。なるほどありそうです。

[dotbot]

そうなんですね...浮動小数点の比較の際は気を付けます。

Answer 3

ちょっと面白そうだったので詳しく調べてみました。
otnさんの回答を参考にBorland C++ 5.5.1 for Win32で(unionの宣言でエラーが出たので位置だけ変えて)コンパイルしたところ、以下のようになりました。(抜粋)

  0040115C: C7 45 F8 33 33 33  mov         dword ptr [ebp-8],33333333h
            33
  00401163: C7 45 FC 33 33 24  mov         dword ptr [ebp-4],40243333h
            40
  0040116A: C7 45 F0 66 66 66  mov         dword ptr [ebp-10h],66666666h
            66
  00401171: C7 45 F4 66 66 24  mov         dword ptr [ebp-0Ch],40246666h
            40
  00401178: C7 45 E8 CD CC CC  mov         dword ptr [ebp-18h],0CCCCCCCDh
            CC
  0040117F: C7 45 EC CC 4C 34  mov         dword ptr [ebp-14h],40344CCCh
            40
…
  004011A0: D9 05 04 12 40 00  fld         dword ptr ds:[00401204h]
  004011A6: DC 45 F8           fadd        qword ptr [ebp-8]
  004011A9: DC 45 F0           fadd        qword ptr [ebp-10h]
  004011AC: DD 1B              fstp        qword ptr [ebx]
…
  004011DD: D9 05 04 12 40 00  fld         dword ptr ds:[00401204h]
  004011E3: DC 45 F8           fadd        qword ptr [ebp-8]
  004011E6: DC 45 F0           fadd        qword ptr [ebp-10h]
  004011E9: DC 5D E8           fcomp       qword ptr [ebp-18h]
…
  00401203: 00 00              add         byte ptr [eax],al
  00401205: 00 80 25 55 8B EC  add         byte ptr [eax+EC8B5525h],al

動作を説明すると、
10.1(をdouble精度で表したもの; 以下同)と10.2と10.3をメモリに用意、
x87命令で
DBL_EPSILONと10.1と10.2をそれぞれ足し、doubleとしてメモリに書く
DBL_EPSILONと10.1と10.2をそれぞれ足し、20.3と比較
という動作をしています。
(ところでDBL_EPSILONはfloatで持ってるんですね。ちょっと意外)

一方そのへんのgccでコンパイルしたところ、x87命令は使われずSSE命令が使われていました。これはdouble精度の演算です。
つまりotnさんの言うように「丸めの仕方がコンパイラによって違う」のではなく、演算の精度がそもそも異なる可能性が高いです。丸めの仕方(切り上げか切り下げかなど)がコンパイラによって違うというのは考えづらいと思います。

x87精度の演算でDBL_EPSILON+10.1＋10.2は20.3とは異なりますが、その値をdoubleへ丸めると20.3に丸まります。
一方double精度でDBL_EPSILON+10.1=10.1であり、10.1+10.2は丸めた際に20.3とは異なる値です。