202202701

sqrt(0x1.fffffffffffffp1) の値が Linux と UCRT64/MSYS2 で違う。他の数学関数と違い、 sqrt は exact な値を丸めた値に一致することが保証されているのではないか。

[[gnu::noinline]] double call_sqrt(double x) {
	return std::sqrt(x);
}

void dump(double x) {
    std::printf("%+.18e  %+.013a\n", x, x);
}

int main() {
    //fesetround(FE_TONEAREST); // no effects.
    double x = 0x1.fffffffffffffp1;
    dump(std::sqrt(x));
    //dump(__builtin_sqrt(x)); // always the same result as std::sqrt().
    dump(call_sqrt(x));
    dump(_mm_cvtsd_f64(_mm_sqrt_pd(_mm_set_sd(x))));
}

(clang|gcc)/Linux -O[03]
    +1.999999999999999778e+00  +0x1.fffffffffffffp+0
    +1.999999999999999778e+00  +0x1.fffffffffffffp+0
    +1.999999999999999778e+00  +0x1.fffffffffffffp+0
gcc/ucrt64 -O0:
    +2.000000000000000000e+00  +0x1.0000000000000p+1
    +2.000000000000000000e+00  +0x1.0000000000000p+1
    +1.999999999999999778e+00  +0x1.fffffffffffffp+0
gcc/ucrt64 -O3:
    +1.999999999999999778e+00  +0x1.fffffffffffffp+0
    +1.999999999999999778e+00  +0x1.fffffffffffffp+0
    +1.999999999999999778e+00  +0x1.fffffffffffffp+0
clang/ucrt64 -O0:
    +2.000000000000000000e+00  +0x1.0000000000000p+1
    +2.000000000000000000e+00  +0x1.0000000000000p+1
    +1.999999999999999778e+00  +0x1.fffffffffffffp+0
clang/ucrt64 -O3:
    +2.000000000000000000e+00  +0x1.0000000000000p+1
    +1.999999999999999778e+00  +0x1.fffffffffffffp+0
    +2.000000000000000000e+00  +0x1.0000000000000p+1

結果がバラバラに見えるが、何が起こっているのか。最適化しない場合、 Linux では glibc の sqrt を呼ぶが、 MSYS2 では埋め込みで…なんと x87 FPU を叩くコードを吐く。いやいやいや。つまり 80-bit 精度の sqrt の結果が丸められて 2 になっている。

一方、最適化した場合は SSE を使うか定数が埋め込まれる。後者はコンパイラのビルドオプションか何かの影響で、 64-bit 演算の結果に一致したり 80-bit の結果に一致したりしているのだと思う。