这是sqrt函数的一个bug吗？

Question

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

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指令会四舍五入结果？

- Muhammad

Answer 1

sqrt不会对任何东西进行四舍五入——当你将整数转换为双精度浮点数时，你才会这样做。双精度浮点数无法在不丢失精度的情况下表示所有64位整数。特别是从2⁵³开始，有多个整数将被表示为相同的双精度浮点值。

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

- sepp2k

---

谢谢，那么我可以取低32位计算其平方根 - 称之为（a）- 然后从64位数字中减去低32位并获取该数字的平方根 - 称之为（b）- 现在我所需要做的就是将（a）乘以（b），我们就完成了。 - Muhammad

@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)。 - MSalters

@MSalters 什么？这里的错误是1。他期望 ~0u = 0xFFFFFFFF。但实际得到的是0x100000000 = 0xFFFFFFFF + 1。 - sepp2k

@sepp2k：抱歉，我又混淆了（不太熟悉x86汇编），所以正确的结果应该是 0x00000000FFFFFFFF？一开始我读取的结果是从 EAX 中获取的，不能是 0x100000000。 - MSalters

@MSalters现在你提到了，他的代码不能返回0x100000000，因为它不适合无符号整数...虽然这是他所说的（如果他返回一个unsigned long long并相应地设置寄存器，那么他将得到这个值）。可能他没有查看返回值（应该是0），但他使用调试器找出函数返回之前'value'的值为0x100000000。或者他确实查看了返回值，然后使用调试器找出为什么它是0，然后发现'value'是0x100000000。 - sepp2k

Answer 2

你对你调用的 C++ 函数不是很清楚。`sqrt` 是一个重载的名称。你可能想要使用 `sqrt(double(~0ull))`。没有接受 `unsigned long long` 的 `sqrt` 重载。

- MSalters

---

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