整数 sqrt 的精度

Question

我有一个这样的循环：

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 我的问题是 - 所有 n 都是 cut=sqrt(n)>=y 吗？ 如果 cut=y-1 则有问题。例如。如果 n = 120 并且 cut = 10 没关系，但如果 n=121 (11^2) 并且 cut 仍然是 10，那么它将不起作用。

我首先担心的是浮点数的小数部分只有 23 位和双 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 总是向上取整就可以了。

Answer 1

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

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

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

将乘法替换为 lea 或移位和两个加法。

---

抛开这些注释，让我们看看您的问题。不，您不能只使用i < sqrt(n)。反例：n = 0x20000000000000。假设遵守 IEEE-754，您将拥有cut = 0x5a82799，而cut*cut 是0x1ffffff8eff971。

但是，基本的浮点错误分析表明，计算 sqrt(n)（转换为整数之前）的错误以 ULP 的 3/4 为界。所以你可以放心使用：

uint32_t cut = sqrt(n) + 1;

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

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

---

编辑：z boson 澄清说，就他的目的而言，这仅在 n 是一个精确的正方形时才重要（否则，得到一个“太小”的 cut 值是可以接受）。在这种情况下，无需调整，直接使用即可：

uint32_t cut = sqrt(n);

为什么这是真的？实际上，这很简单。将n 转换为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 的相对大小实际上会随着 binade 的变化而变化——我试图在没有的情况下给出一些证明的味道太沉迷于细节。这一切都可以做得非常严格，只是对于 Stack Overflow 来说有点罗嗦。

Answer 2

如果您可以访问 IEEE 754 双精度浮点，我所有的答案都是无用的，因为斯蒂芬佳能证明了两者

• 一种避免 imul in loop 的简单方法
• 一种计算天花板平方的简单方法

否则，如果由于某种原因您的平台不符合 IEEE 754 标准，或者只有单精度，您可以使用简单的 Newton-Raphson 循环获得平方根的整数部分。例如在 Squeak Smalltalk 中，我们在 Integer 中有这个方法：

sqrtFloor
    "Return the integer part of the square root of self"

    | guess delta |
    guess := 1 bitShift: (self highBit + 1) // 2.
    [
        delta := (guess squared - self) // (guess + guess).
        delta = 0 ] whileFalse: [
            guess := guess - delta ].
    ^guess - 1

其中 // 是整数除法商的运算符。
如果最初的猜测超过了精确解，就像这里的情况一样，可以避免最终保护guess*guess <= self ifTrue: [^guess].。
使用近似浮点 sqrt 初始化不是一个选项，因为整数任意大并且可能溢出

但是在这里，您可以使用浮点 sqrt 近似作为初始猜测的种子，我敢打赌，将在很少的循环中找到确切的解决方案。在 C 语言中是：

uint32_t sqrtFloor(uint64_t n)
{
    int64_t diff;
    int64_t delta;
    uint64_t guess=sqrt(n); /* implicit conversions here... */
    while( (delta = (diff=guess*guess-n) / (guess+guess)) != 0 )
        guess -= delta;
    return guess-(diff>0);
}

这是一些整数乘法和除法，但在主循环之外。

Answer 3

您正在寻找的是一种计算自然数平方根有理上限的方法。连分数是你需要的见wikipedia。

对于 x>0，有 .

为了使符号更简洁，将上面的公式改写为

通过在每个递归深度去除尾项 (x-1)/2's 来截断连分数，得到 sqrt(x) 的一系列近似值，如下所示：

上限出现在具有奇数行号的行上，并且越来越紧。当上限与其相邻下限之间的距离小于 1 时，您需要该近似值。用那个值作为cut的值，这里cut必须是浮点数，解决问题。

对于非常大的数字，应该使用有理数，所以在整数和浮点数之间的转换过程中不会丢失精度。