Libc hypot 函数似乎 return double 类型的错误结果...为什么？

Question

#include <tgmath.h>
#include <iostream>
int main(int argc, char** argv) {

        #define NUM1 -0.031679909079365576
        #define NUM2 -0.11491794452567111

        std::cout << "double precision :"<< std::endl;
        typedef std::numeric_limits< double > dbl;
        std::cout.precision(dbl::max_digits10);
        std::cout << std::hypot((double)NUM1, (double)NUM2);
        std::cout << " VS sqrt :" << sqrt((double )NUM1*(double )NUM1 
                                  + (double )NUM2*(double )NUM2) << std::endl;

        std::cout << "long double precision :"<< std::endl;
        typedef std::numeric_limits<long double > ldbl;
        std::cout.precision(ldbl::max_digits10);
        std::cout << std::hypot((long double)NUM1, (long double)NUM2);
        std::cout << " VS sqrt :" << sqrt((long double )NUM1*(long double )NUM1 + (long double )NUM2*(long double )NUM2);
}

Returns 在 Linux 下（Ubuntu 18.04 clang 或 gcc，无论优化如何，glic 2.25）：

双精度： 0.1192046585217293 VS 平方：0.11920465852172932

长双精度： 0.119204658521729311251 VS 开方：0.119204658521729311251

根据 cppreference :

实现通常保证精度小于1 ulp（最后一位的单位）：GNU、BSD、Open64 std::hypot(x, y) 等同于 std::abs(std::complex(x,y)) POSIX 指定仅当两个参数都低于正规且正确结果也低于正规时才可能发生下溢（这禁止天真的实现）

所以，hypot((double)NUM1, (double)NUM2) 应该 return 0.11920465852172932，我想（作为天真的 sqrt 实现）。在windows上，使用MSVC 64位，是这样的。

为什么我们在使用 glibc 时会看到这种差异？如何解决这种不一致？

Answer 1

0.11920465852172932 表示为 0x1.e84324de1b576p-4（双精度）
0.11920465852172930 表示为 0x1.e84324de1b575p-4（双精度）
0.119204658521729311251 是 long-double 结果，我们可以假设它在小数点后几位是正确的。即准确结果更接近四舍五入的结果。

这些 FP 位模式仅在尾数的低位（也称为有效数）不同，确切的结果在它们之间。 所以它们每个都有不到 1 ulp 的舍入误差， 实现了典型实现（包括 glibc）的目标。

与 IEEE-754“基本”操作 (add/sub/mul/div/sqrt) 不同，hypot 不需要“正确舍入”。这意味着 <= 0.5 ulp 的误差。对于硬件不直接提供的操作，实现这一点会慢得多。（例如，使用至少几个额外的绝对正确位进行扩展精度计算，因此您可以四舍五入到最接近的双精度，就像硬件对基本操作所做的那样）

碰巧在这种情况下，朴素的计算方法产生了正确舍入的结果，而 glibc 的 std::hypot 的“安全”实现（在加法前对小数进行平方时必须避免下溢）产生了结果误差 >0.5 但 <1 ulp。

您没有指定您是否在 32 位模式下使用 MSVC。

大概 32 位模式将使用 x87 进行 FP 数学运算，从而提供额外的临时精度。尽管某些 MSVC 版本的 CRT 代码将 x87 FPU 的内部精度设置为在每次操作后四舍五入为 53 位尾数，因此它的行为类似于使用实际 double 的 SSE2，只是指数范围更宽。参见 Bruce Dawson's blog post。

所以我不知道除了运气之外是否还有其他原因让 MSVC std::hypot 得到了正确的舍入结果。

注意MSVC中的long double与64位的double是同一类型；该 C++ 实现不公开 x86 / x86-64 的 80 位硬件扩展精度类型。（64 位尾数）。

Libc hypot 函数似乎 return double 类型的错误结果...为什么？

Libc hypot function seems to return incorrect results for double type... why?

c++

floating-point

sse

glibc

hypotenuse