整数 sqrt 的精度

我有一个这样的循环：

for(uint64_t i=0; i*i<n; i++) {

这需要每次迭代都进行乘法。如果我可以在循环之前计算 sqrt 那么我可以避免这种情况。

unsigned cut = sqrt(n)
for(uint64_t i=0; i<cut; i++) {

就我而言，如果 sqrt 函数向上舍入到下一个整数是可以的，但如果向下舍入则不行。

我的问题是： sqrt 函数是否足够准确，可以在所有情况下执行此操作？

编辑：让我列出一些案例。如果 n 是完全平方数，那么n = y^2我的问题是 - 是cut=sqrt(n)>=y对于所有 n?如果 cut=y-1 那么就有问题了。例如。如果 n = 120 并且 cut = 10 没问题，但是如果 n=121 (11^2) 并且 cut 仍然是 10 那么它就不起作用。

我首先担心的是 float 的小数部分只有 23 位和 double 52，因此它们无法存储某些 32 位或 64 位整数的所有数字。不过，我认为这不是问题。假设我们想要某个数字 y 的 sqrt，但我们无法存储 y 的所有数字。如果我们让可以存储的 y 分数为 x，我们可以写成 y = x + dx，那么我们要确保我们选择的任何 dx 都不会移动到下一个整数。

sqrt(x+dx) < sqrt(x) + 1  //solve
dx < 2*sqrt(x) + 1 
// e.g for x = 100 dx < 21
// sqrt(100+20) < sqrt(100) + 1

浮点数可以存储 23 位，因此我们让 y = 2^23 + 2^9。这已经足够了，因为 2^9

unsigned xi = -1-1;
printf("%u %u\n", xi, (unsigned)(float)xi);  //4294967294 4294967295
printf("%u %u\n", (unsigned)sqrt(xi), (unsigned)sqrtf(xi));  //65535 65536

由于 float 无法存储 2^31-2 的所有数字，因此 double 可以得到不同的 sqrt 结果。但 sqrt 的浮点版本要大一个整数。这就是我要的。对于 64 位整数，只要 double 的 sqrt 始终向上舍入就可以了。

首先，整数乘法确实非常便宜。只要每个循环迭代有多个工作周期和一个备用执行槽，它就应该通过在大多数非小型处理器上重新排序而完全隐藏。

如果您确实有一个整数乘法速度非常慢的处理器，那么一个真正聪明的编译器可能会将您的循环转换为：

for (uint64_t i = 0, j = 0; j < cut; j += 2*i+1, i++)

将乘法替换为lea或一个班次和两个加法。

撇开这些注释不谈，让我们看看你所说的问题。不，你不能只使用i < sqrt(n)。反例：n = 0x20000000000000。假设遵守 IEEE-754，您将拥有cut = 0x5a82799, and cut*cut is 0x1ffffff8eff971.

然而，基本的浮点误差分析表明计算中的误差sqrt(n)（转换为整数之前）以 ULP 的 3/4 为界。所以您可以安全地使用：

uint32_t cut = sqrt(n) + 1;

并且您最多将执行一次额外的循环迭代，这可能是可以接受的。如果您想完全精确，请改用：

uint32_t cut = sqrt(n);
cut += (uint64_t)cut*cut < n;

Edit: z boson澄清说，就他的目的而言，这仅在以下情况下才重要：n是一个精确的平方（否则，得到的值为cut“太小了”是可以接受的）。在这种情况下，不需要进行调整，可以安全地使用：

uint32_t cut = sqrt(n);

为什么这是真的？实际上，这很简单。转换n to double引入扰动：

double_n = n*(1 + e)

这满足|e| < 2^-53。该值的数学平方根可以展开如下：

square_root(double_n) = square_root(n)*square_root(1+e)

现在，自从n假设是最多 64 位的完美平方，square_root(n)是一个最多 32 位的精确整数，并且是我们希望计算的数学精确值。来分析square_root(1+e)术语，使用泰勒级数1:

square_root(1+e) = 1 + e/2 + O(e^2)
                 = 1 + d with |d| <~ 2^-54

因此，数学上的精确值square_root(double_n)与[1]所需的精确答案的距离不到 ULP 的一半，并且必须舍入到该值。

[1] 我在这里滥用相对误差估计，其中 ULP 的相对大小实际上在二进制文件中有所不同 - 我试图在不陷入太多困境的情况下提供一点证明的味道往下看细节。这一切都可以做得非常严格，只是对于 Stack Overflow 来说有点啰嗦。