MSVC 在最大值及其附近为 std::uint64_t 产生不正确的平方根。错误还是预期？

Question

问题： 当我为整数配对函数编写一些单元测试时，我注意到 "bug" 不断出现。 std::sqrt 和 std::sqrtl 在最大 std::uint64_t 值附近产生相同的错误结果。这仅仅是因为浮点舍入错误或微软编译器中的错误而可以预料的吗？如果是预料之中的，那么有没有办法在不求助于 128 位数据或使用迭代算法来提高精度的情况下规避这个问题？也许我应该注意一个编译器标志？

示例： 使用这些配对函数，要取消配对输入，您必须获得输入平方根的底部，我无法绕过最大值。我把这个问题抽象成一个简单的例子来演示：

    auto max = std::numeric_limits<std::uint64_t>::max(); //18446744073709551615
    auto square_root = std::floorl(std::sqrtl(max));       //4294967296.0000000
    auto max2 = max - 1;
    auto square_root2 = std::floorl(std::sqrtl(max2));     //4294967296.0000000
    auto max3 = max - 2;
    auto square_root3 = std::floorl(std::sqrtl(max3));     //4294967296.0000000
    auto max4 = max - 5;
    auto square_root4 = std::floorl(std::sqrtl(max4));     //4294967296.0000000
    auto max5 = max / 10;
    auto square_root5 = std::floorl(std::sqrtl(max5));     //1358187913.0000000

根据wolframalpha，前4个平方根的正确值如下

//(respectively)
4294967295.99999999988358467817306518554529727818955797
4294967295.99999999976716935634613037108743911275823190
4294967295.99999999965075403451919555662642550370602178
4294967295.99999999930150806903839111322445201482408715

//To clarify, the floor of 4294967295.99999999930 should be 
//4294967295.0, where the code produces
//4294967296.0

备注：

我应该说，我预计这不会导致配对问题，这更像是一个概念性问题，如果这不是预期的行为，我会提交一个错误。
我正在使用 visual studio 2017 并在 x86 和 x64 版本中看到此行为
有关配对功能的更多信息，请访问：http://szudzik.com/ElegantPairing.pdf（虽然我不是在寻找不同的解决方案）

额外信息

感谢@Phil1970 的简洁回答和解释——除了他的回答，我还想分享一段我从 the article 中找到的有趣摘录，因为它与我发布的问题有关[=16] =]

Brown [1981] has proposed axioms for floating-point that include most of the existing floating-point hardware. However, proofs in this system cannot verify the algorithms of sections Cancellation and Exactly Rounded Operations, which require features not present on all hardware. Furthermore, Brown's axioms are more complex than simply defining operations to be performed exactly and then rounded. Thus proving theorems from Brown's axioms is usually more difficult than proving them assuming operations are exactly rounded.

There is not complete agreement on what operations a floating-point standard should cover. In addition to the basic operations +, -, × and /, the IEEE standard also specifies that square root, remainder, and conversion between integer and floating-point be correctly rounded. It also requires that conversion between internal formats and decimal be correctly rounded (except for very large numbers).

虽然 "correctly rounded (except for very large numbers)" 对我来说听起来有点可疑，这意味着大数的浮点规范 不需要 正确的舍入，我相信所经历的行为在这个例子中确实是预期的。这里要注意的是，浮点误差和舍入不限于使用小数点后小数的 mantessa 中的细微差别。意识到与整数数据类型之间的隐式转换的大数的（不）准确性也很重要。

Answer 1

一个std::uint64_t是64位

在 Visual Studio 的情况下，double 或 long double 也是 64 位。

IEE 754 格式的双精度数有 53 位尾数：https://en.wikipedia.org/wiki/Double-precision_floating-point_format。

64 - 53 位 = 11 位。由 1 组成的 11 位是 2047。由于数字四舍五入到最接近的双精度数，因此对于 max - 1023 到 max 之间的任何数字，我们得到 18446744073709551616。对于 max - 3070 和 max - 1024 之间的数字，我们得到 18446744073709549568。这两个数字之间的差是 2048 (2^11)。所以一个 double 值在 2^64 左右精确到 ±2048.

MSVC 使用 64 位长双精度。其他编译器可能使用 80 位甚至 128 位：https://en.wikipedia.org/wiki/Long_double.

我使用以下代码测试了 max

的一些偏移量

void test(std::uint64_t offset)
{
    auto max = std::numeric_limits<std::uint64_t>::max(); //18446744073709551615

    std::cout << "max - " << std::setw(10) << offset << " : " 
        << std::setprecision(20) << static_cast<double>(max - offset) << "\n";
}

这是一些数字的输出：

max -          0 : 18446744073709551616
max -       1023 : 18446744073709551616
max -       1024 : 18446744073709549568
max -       2047 : 18446744073709549568
max -       2048 : 18446744073709549568
max -       3070 : 18446744073709549568
max -       3071 : 18446744073709547520

补充信息

是的，如果转换整数值需要比浮点数的尾数所允许的位数更多的位数，则预计会降低一些精度。通常，使用 double 对于通常的应用来说这不是问题，因为一个大于 2^53 (9.007e15) 的数字在转换为 IEEE 754 double 时开始失去精度。

我不知道从整数到双精度的转换是否需要四舍五入，或者它是否由实现定义。相反方向，总是t运行cated。将文本转换为双精度值（在编译时或运行时），我相信它总是四舍五入到最接近的可表示值。

话虽如此，像 4294967295.99999999930150806903839111322445201482408715 这样的平方根四舍五入前的结果不能用足够精度的双精度表示。它需要 20 位有效数字来表示 4294967295.9999999993，而 64 位双精度数大约有 15 位有效数字。

那你能做什么？

如果您不需要处理那么高的精度，请忽略该问题。在单元测试中使用 2^50 左右的数字作为 max。
使用其他支持更大浮点类型的编译器，如 Intel 或 gcc。
使用专为大型精度浮点计算而制作的库。
使用一些公式和可能的迭代估计平方根。
将数字转换回整数。检查是否小于 2^32。如果是平方，如果高于原始数，则减一。

MSVC 在最大值及其附近为 std::uint64_t 产生不正确的平方根。错误还是预期？

MSVC yields incorrect square root for std::uint64_t at and around max. Bug or expected?

c++

msvcrt