快速整数平方根近似

Question

我目前正在寻找一个非常快的整数平方根近似值，其中 floor(sqrt(x)) <= veryFastIntegerSquareRoot(x) <= x

平方根例程用于计算素数，如果仅sqrt(x)检查小于或等于 的值是否为 的除数，则计算速度会快得多x。

我目前拥有的是来自 Wikipedia 的这个函数，稍微调整了一下以使用 64 位整数。

因为我没有其他函数可以比较（或者更准确地说，该函数对于我的目的来说太精确了，而且它可能需要更多的时间，而不是高于实际结果。）

Answer 1

无循环/无跳跃（好吧：几乎；-）牛顿-拉夫森：

/* static will allow inlining */
static unsigned usqrt4(unsigned val) {
    unsigned a, b;

    if (val < 2) return val; /* avoid div/0 */

    a = 1255;       /* starting point is relatively unimportant */

    b = val / a; a = (a+b) /2;
    b = val / a; a = (a+b) /2;
    b = val / a; a = (a+b) /2;
    b = val / a; a = (a+b) /2;

    return a;
}

对于 64 位整数，您将需要更多步骤（我的猜测：6）

Answer 2

floor(sqrt(x))精确计算

这是我的解决方案，它基于维基百科上提出的位猜测方法。不幸的是，维基百科上提供的伪代码有一些错误，所以我不得不做出一些调整：

unsigned char bit_width(unsigned long long x) {
    return x == 0 ? 1 : 64 - __builtin_clzll(x);
}

// implementation for all unsigned integer types
unsigned sqrt(const unsigned n) {
    unsigned char shift = bit_width(n);
    shift += shift & 1; // round up to next multiple of 2

    unsigned result = 0;

    do {
        shift -= 2;
        result <<= 1; // leftshift the result to make the next guess
        result |= 1;  // guess that the next bit is 1
        result ^= result * result > (n >> shift); // revert if guess too high
    } while (shift != 0);

    return result;
}

bit_width可以在常数时间内求值，并且循环最多会迭代ceil(bit_width / 2)。因此，即使对于 64 位整数，这最多也需要 32 次基本算术和按位运算迭代。

与迄今为止提出的所有其他答案不同，这实际上为您提供了最佳的近似值，即floor(sqrt(x))。对于任何 x ²，这将准确返回 x。

_{使用 log 2} (x)进行猜测

如果这对您来说仍然太慢，您可以仅根据二进制对数进行猜测。基本思想是我们可以使用 2 ^x/2sqrt计算任意数字 2 ^x的。可能有余数，所以我们不能总是精确地计算它，但我们可以计算上限和下限。^x/2

例如：

1. 我们有 25
2. 计算floor(log_2(25)) = 4
3. 计算ceil(log_2(25)) = 5
4. 下限：pow(2, floor(4 / 2)) = 4
5. 上限：pow(2, ceil(5 / 2)) = 8

事实上，实际的sqrt(25) = 5. 我们发现sqrt(16) >= 4和sqrt(32) <= 8。这意味着：

4 <= sqrt(16) <= sqrt(25) <= sqrt(32) <= 8
            4 <= sqrt(25) <= 8

这就是我们如何实现这些猜测，我们将其称为sqrt_lo和sqrt_hi。

// this function computes a lower bound
unsigned sqrt_lo(const unsigned n) noexcept
{
    unsigned log2floor = bit_width(n) - 1;
    return (unsigned) (n != 0) << (log2floor >> 1);
}

// this function computes an upper bound
unsigned sqrt_hi(const unsigned n)
{
    bool isnt_pow2 = ((n - 1) & n) != 0; // check if n is a power of 2
    unsigned log2ceil = bit_width(n) - 1 + isnt_pow2;
    log2ceil += log2ceil & 1; // round up to multiple of 2
    return (unsigned) (n != 0) << (log2ceil >> 1);
}

对于这两个函数，以下语句始终正确：

sqrt_lo(x) <= floor(sqrt(x)) <= sqrt(x) <= sqrt_hi(x) <= x

请注意，如果我们假设输入永远不为零，则(unsigned) (n != 0)可以简化为1并且该语句仍然为真。

这些功能可以O(1)通过硬件__builtin_clzll指令进行评估。他们只给出数字 2 ^2x、所以256、64、16等的精确结果。