这是 sqrt 函数中的错误吗答案

【问题标题】：is this a bug in sqrt function这是 sqrt 函数中的错误吗
【发布时间】：2012-05-30 11:55:20
【问题描述】：

我创建了一个应用程序来计算 64 位范围内的素数，所以当我尝试使用来自 math.h 的 sqrt 函数计算 64 位数的平方根时，我发现答案不准确，例如当输入是~0ull 答案应该是~0u 但我得到的是0x100000000 这是不对的，所以我决定使用汇编x86 语言创建我自己的版本，看看这是否是一个错误，这是我的功能：

inline unsigned prime_isqrt(unsigned long long value)
{
    const unsigned one = 1;
    const unsigned two = 2;

    __asm
    {
        test dword ptr [value+4], 0x80000000
        jz ZERO
        mov eax, dword ptr [value]
        mov ecx, dword ptr [value + 4]

        shrd eax, ecx, 1
        shr  ecx, 1
        mov  dword ptr [value],eax 
        mov  dword ptr [value+4],ecx

        fild value
        fimul two
        fiadd one
        jmp REST
ZERO: 
        fild value
REST: 
        fsqrt
        fisttp value
        mov eax, dword ptr [value]
    }
}

输入是奇数，求平方根。当我用相同的输入测试我的函数时，结果是一样的。

我不明白为什么这些函数会对结果进行舍入，或者具体来说为什么sqrt 指令会舍入结果？

【问题讨论】：

标签： c++ assembly x86 standard-library x87

【解决方案1】：

sqrt 不四舍五入 - 当您将整数转换为双精度时您会这样做。双精度不能代表 64 位整数可以不损失精度的所有数字。具体从 2⁵³ 开始，有多个整数将表示为相同的双精度值。

因此，如果将大于 2⁵³ 的整数转换为 double，则会丢失一些最低有效位，这就是为什么 (double)(~0ull) 是 18446744073709552000.0，而不是 18446744073709551615.0（或者更准确地说是后者实际上等于前者，因为它们代表相同的双数）。

【讨论】：

谢谢，所以我可以用较低的 32 位计算它的平方根 - 称之为 (a) - 然后从 64 位数字中减去较低的 32 位并得到该数字的平方根 - 调用它 (b) - 现在我需要做的就是将 (a) 乘以 (b)，我们就完成了。
@Muhammadalaa：嗯，不。 64 位在数学上可以写成 (H * 2^32 + L)，其中 H 和 L 各为 32 位。您声明sqrt(H * 2^32 + L) = sqrt(L) * sqrt(H*2^32)。事实并非如此。但是，您可以将其计算为sqrt(H) * sqrt(2^32+L/H)。
舍入会给你一个大约 2^-53 的误差，但这里的误差是 2 倍：~0u == UINT_MAX 与 0x80000000 = UINT_MAX/2。我不认为这可以解释它。
@MSalters 什么？这里的错误是 1。他期望 ~0u = 0xFFFFFFFF。他得到了 0x100000000 = 0xFFFFFFFF + 1。
@MSalters 现在你提到它，他的代码不能返回 0x100000000 因为这不适合未签名的它......这就是他说他得到的（这就是他会得到的如果他返回一个 unsigned long long 并相应地设置寄存器）。可能他没有看返回值（应该是0），但是他用调试器发现value的值在函数返回之前是0x100000000。或者他确实查看了返回值，然后使用调试器找出为什么它是0，然后发现value是0x100000000。

【解决方案2】：

您对所调用的 C++ 函数不是很清楚。 sqrt 是重载名称。你可能想要sqrt(double(~0ull))。没有sqrt 重载需要unsigned long long。

【讨论】：

他知道这一点。他在他的代码中明确地将数字转换为双精度数。