获取sqrt(n)整数部分的最快方法？

Question

我们知道如果n不是一个完美的正方形,那么sqrt(n)就不会是一个整数.由于我只需要整数部分,我觉得调用sqrt(n)不会那么快,因为计算小数部分也需要时间.

所以我的问题是,

我们是否只能获得sqrt(n)的整数部分而不计算实际值sqrt(n)？算法应该比sqrt(n)(在<math.h>或中定义<cmath>)更快？

如果可能,您也可以在asm块中编写代码.

Answer 1

我会尝试Fast Inverse Square Root技巧.

这是一种在1/sqrt(n)没有任何分支的情况下获得非常好的近似的方法,基于一些比特错误,因此不可移植(特别是在32位和64位平台之间).

一旦得到它,你只需要反转结果,并取整数部分.

当然,可能有更快的技巧,因为这个有点圆.

编辑:让我们做!

首先是一个小帮手:

// benchmark.h
#include <sys/time.h>

template <typename Func>
double benchmark(Func f, size_t iterations)
{
  f();

  timeval a, b;
  gettimeofday(&a, 0);
  for (; iterations --> 0;)
  {
    f();
  }
  gettimeofday(&b, 0);
  return (b.tv_sec * (unsigned int)1e6 + b.tv_usec) -
         (a.tv_sec * (unsigned int)1e6 + a.tv_usec);
}

然后主体:

#include <iostream>

#include <cmath>

#include "benchmark.h"

class Sqrt
{
public:
  Sqrt(int n): _number(n) {}

  int operator()() const
  {
    double d = _number;
    return static_cast<int>(std::sqrt(d) + 0.5);
  }

private:
  int _number;
};

// http://www.codecodex.com/wiki/Calculate_an_integer_square_root
class IntSqrt
{
public:
  IntSqrt(int n): _number(n) {}

  int operator()() const 
  {
    int remainder = _number;
    if (remainder < 0) { return 0; }

    int place = 1 <<(sizeof(int)*8 -2);

    while (place > remainder) { place /= 4; }

    int root = 0;
    while (place)
    {
      if (remainder >= root + place)
      {
        remainder -= root + place;
        root += place*2;
      }
      root /= 2;
      place /= 4;
    }
    return root;
  }

private:
  int _number;
};

// http://en.wikipedia.org/wiki/Fast_inverse_square_root
class FastSqrt
{
public:
  FastSqrt(int n): _number(n) {}

  int operator()() const
  {
    float number = _number;

    float x2 = number * 0.5F;
    float y = number;
    long i = *(long*)&y;
    //i = (long)0x5fe6ec85e7de30da - (i >> 1);
    i = 0x5f3759df - (i >> 1);
    y = *(float*)&i;

    y = y * (1.5F - (x2*y*y));
    y = y * (1.5F - (x2*y*y)); // let's be precise

    return static_cast<int>(1/y + 0.5f);
  }

private:
  int _number;
};


int main(int argc, char* argv[])
{
  if (argc != 3) {
    std::cerr << "Usage: %prog integer iterations\n";
    return 1;
  }

  int n = atoi(argv[1]);
  int it = atoi(argv[2]);

  assert(Sqrt(n)() == IntSqrt(n)() &&
          Sqrt(n)() == FastSqrt(n)() && "Different Roots!");
  std::cout << "sqrt(" << n << ") = " << Sqrt(n)() << "\n";

  double time = benchmark(Sqrt(n), it);
  double intTime = benchmark(IntSqrt(n), it);
  double fastTime = benchmark(FastSqrt(n), it);

  std::cout << "Number iterations: " << it << "\n"
               "Sqrt computation : " << time << "\n"
               "Int computation  : " << intTime << "\n"
               "Fast computation : " << fastTime << "\n";

  return 0;
}

结果如下:

sqrt(82) = 9
Number iterations: 4096
Sqrt computation : 56
Int computation  : 217
Fast computation : 119

// Note had to tweak the program here as Int here returns -1 :/
sqrt(2147483647) = 46341 // real answer sqrt(2 147 483 647) = 46 340.95
Number iterations: 4096
Sqrt computation : 57
Int computation  : 313
Fast computation : 119

在预期的情况下,快速计算比Int计算执行得更好.

哦,顺便说一句,sqrt更快:)

Answer 2

编辑:这个答案是愚蠢的 - 使用 (int) sqrt(i)

与剖析后适当设置(-march=native -m64 -O3)上面是一个很大更快.

好吧,有点老问题,但尚未给出"最快"的答案.最快(我认为)是二进制平方根算法,在此Embedded.com文章中有详细解释.

它基本上归结为:

unsigned short isqrt(unsigned long a) {
    unsigned long rem = 0;
    int root = 0;
    int i;

    for (i = 0; i < 16; i++) {
        root <<= 1;
        rem <<= 2;
        rem += a >> 30;
        a <<= 2;

        if (root < rem) {
            root++;
            rem -= root;
            root++;
        }
    }

    return (unsigned short) (root >> 1);
}

在我的机器上(Q6600,Ubuntu 10.10)我通过取数字1-100000000的平方根进行分析.使用iqsrt(i)耗时2750毫秒.用(unsigned short) sqrt((float) i)了3600ms.这是使用完成的g++ -O3.使用-ffast-math编译选项的时间分别为2100ms和3100ms.请注意,这甚至不使用单行汇编程序,因此它可能仍然会更快.

上面的代码适用于C和C++,并且对Java也有一些小的语法更改.

对于有限范围更有效的是二分搜索.在我的机器上,它将上面的版本从水中吹出4倍.可悲的是它的范围非常有限:

#include <stdint.h>

const uint16_t squares[] = {
    0, 1, 4, 9,
    16, 25, 36, 49,
    64, 81, 100, 121,
    144, 169, 196, 225,
    256, 289, 324, 361,
    400, 441, 484, 529,
    576, 625, 676, 729,
    784, 841, 900, 961,
    1024, 1089, 1156, 1225,
    1296, 1369, 1444, 1521,
    1600, 1681, 1764, 1849,
    1936, 2025, 2116, 2209,
    2304, 2401, 2500, 2601,
    2704, 2809, 2916, 3025,
    3136, 3249, 3364, 3481,
    3600, 3721, 3844, 3969,
    4096, 4225, 4356, 4489,
    4624, 4761, 4900, 5041,
    5184, 5329, 5476, 5625,
    5776, 5929, 6084, 6241,
    6400, 6561, 6724, 6889,
    7056, 7225, 7396, 7569,
    7744, 7921, 8100, 8281,
    8464, 8649, 8836, 9025,
    9216, 9409, 9604, 9801,
    10000, 10201, 10404, 10609,
    10816, 11025, 11236, 11449,
    11664, 11881, 12100, 12321,
    12544, 12769, 12996, 13225,
    13456, 13689, 13924, 14161,
    14400, 14641, 14884, 15129,
    15376, 15625, 15876, 16129,
    16384, 16641, 16900, 17161,
    17424, 17689, 17956, 18225,
    18496, 18769, 19044, 19321,
    19600, 19881, 20164, 20449,
    20736, 21025, 21316, 21609,
    21904, 22201, 22500, 22801,
    23104, 23409, 23716, 24025,
    24336, 24649, 24964, 25281,
    25600, 25921, 26244, 26569,
    26896, 27225, 27556, 27889,
    28224, 28561, 28900, 29241,
    29584, 29929, 30276, 30625,
    30976, 31329, 31684, 32041,
    32400, 32761, 33124, 33489,
    33856, 34225, 34596, 34969,
    35344, 35721, 36100, 36481,
    36864, 37249, 37636, 38025,
    38416, 38809, 39204, 39601,
    40000, 40401, 40804, 41209,
    41616, 42025, 42436, 42849,
    43264, 43681, 44100, 44521,
    44944, 45369, 45796, 46225,
    46656, 47089, 47524, 47961,
    48400, 48841, 49284, 49729,
    50176, 50625, 51076, 51529,
    51984, 52441, 52900, 53361,
    53824, 54289, 54756, 55225,
    55696, 56169, 56644, 57121,
    57600, 58081, 58564, 59049,
    59536, 60025, 60516, 61009,
    61504, 62001, 62500, 63001,
    63504, 64009, 64516, 65025
};

inline int isqrt(uint16_t x) {
    const uint16_t *p = squares;

    if (p[128] <= x) p += 128;
    if (p[ 64] <= x) p +=  64;
    if (p[ 32] <= x) p +=  32;
    if (p[ 16] <= x) p +=  16;
    if (p[  8] <= x) p +=   8;
    if (p[  4] <= x) p +=   4;
    if (p[  2] <= x) p +=   2;
    if (p[  1] <= x) p +=   1;

    return p - squares;
}

可以在此处下载32位版本:https://gist.github.com/3481770

Answer 3

以下解决方案计算整数部分，含义floor(sqrt(x))准确，没有舍入误差。

其他方法的问题

• 使用floatordouble既不便携也不足够精确
• @orlpisqrt给出了疯狂的结果，例如isqrt(100) = 15
• 基于巨大查找表的方法在 32 位以上不实用
• 使用快速反平方根非常不精确，你最好使用sqrtf
• 牛顿方法需要昂贵的整数除法和良好的初始猜测

我的方法

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

// C++20 also provides std::bit_width in its <bit> header
unsigned char bit_width(unsigned long long x) {
    return x == 0 ? 1 : 64 - __builtin_clzll(x);
}

template <typename Int, std::enable_if_t<std::is_unsigned<Int, int = 0>>
Int sqrt(const Int n) {
    unsigned char shift = bit_width(n);
    shift += shift & 1; // round up to next multiple of 2

    Int result = 0;

    do {
        shift -= 2;
        result <<= 1; // make space for the next guessed bit
        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 次基本算术和按位运算迭代。

编译输出只有大约 20 条指令。

表现

我float通过统一生成输入来将我的方法与 -bases 方法进行基准测试。请注意，在现实世界中，大多数输入将更接近于 0，而不是std::numeric_limits<...>::max()。

• 因为uint32_t这25x比使用性能差std::sqrt(float)
• 因为uint64_t这30x比使用性能差std::sqrt(double)

准确性

与使用浮点数学的方法不同，此方法始终非常准确。

• 使用sqrtf可能会在 [2 ²⁸ , 2 ³² ) 范围内提供不正确的舍入。例如，sqrtf(0xffffffff) = 65536当平方根实际上是 时65535.99999。
• ^{双精度在 [2 60} , 2 ⁶⁴ ) 范围内不一致。例如，sqrt(0x3fff...) = 2147483648当平方根实际上是 时2147483647.999999。

唯一涵盖所有 64 位整数的是 x86 扩展精度long double，因为它可以容纳整个 64 位整数。

结论

正如我所说，这是正确处理所有输入、避免整数除法并且不需要查找表的唯一解决方案。总之，如果您需要一种独立于精度并且不需要庞大的查找表的方法，那么这是您唯一的选择。它可能特别有用constexpr在性能并不重要且获得 100% 准确结果可能更为重要的情况下，它

使用牛顿法的替代方法

当从一个好的猜测开始时，牛顿法可以相当快。对于我们的猜测，我们将向下舍入到 2 的下一个幂并在常数时间内计算平方根。对于任何数字 2 ^{x ，我们可以使用 2 x/2}获得平方根。

template <typename Int, std::enable_if_t<std::is_unsigned_v<Int>, int> = 0>
Int sqrt_guess(const Int n)
{
    Int log2floor = bit_width(n) - 1;
    // sqrt(x) is equivalent to pow(2, x / 2 = x >> 1)
    // pow(2, x) is equivalent to 1 << x
    return 1 << (log2floor >> 1);
}

请注意，这并不完全是 2 ^x/2，因为我们在右移过程中丢失了一些精度。相反，它是 2 ^Floor(x/2)。另请注意，sqrt_guess(0) = 1这实际上是避免在第一次迭代中除以零所必需的：

template <typename Int, std::enable_if_t<std::is_unsigned_v<Int>, int> = 0>
Int sqrt_newton(const Int n)
{
    Int a = sqrt_guess(n);
    Int b = n;
    
    // compute unsigned difference
    while (std::max(a, b) - std::min(a, b) > 1) {
        b = n / a;
        a = (a + b) / 2;
    }

    // a is now either floor(sqrt(n)) or ceil(sqrt(n))
    // we decrement in the latter case
    // this is overflow-safe as long as we start with a lower bound guess
    return a - (a * a > n);
}

这种替代方法的性能大致相当于第一个提案，但通常要快几个百分点。然而，它严重依赖于高效的硬件划分，结果可能会有很大差异。

使用会sqrt_guess产生巨大的差异。1它比用作初始猜测大约快五倍。

Answer 4

虽然我怀疑你可以通过搜索"快速整数平方根"找到很多选项,但这里有一些可能很好用的新想法(每个独立,或者你可以将它们组合起来):

1. static const在您想要支持的域中创建所有完美正方形的数组,并对其执行快速无分支二进制搜索.数组中生成的索引是平方根.
2. 将数字转换为浮点数并将其分解为尾数和指数.将指数减半并将尾数乘以某个神奇因子(找到它的工作).这应该能够给你一个非常接近的近似值.如果它不准确,则包括调整它的最后一步(或者将其用作上面二进制搜索的起点).

Answer 5

我认为Google search提供了很好的文章Calculate an integer square root,讨论了很多可能的快速计算方法,并且有很好的参考文章,我认为这里没有人可以提供比他们更好的(如果有人可以先生成关于它的论文),但如果你读了他们和他们有歧义,那么我们可以帮助你.

Answer 6

如果你不介意近似,那么这个整数sqrt函数如何拼凑在一起.

int sqrti(int x)
{
    union { float f; int x; } v; 

    // convert to float
    v.f = (float)x;

    // fast aprox sqrt
    //  assumes float is in IEEE 754 single precision format 
    //  assumes int is 32 bits
    //  b = exponent bias
    //  m = number of mantissa bits
    v.x  -= 1 << 23; // subtract 2^m 
    v.x >>= 1;       // divide by 2
    v.x  += 1 << 29; // add ((b + 1) / 2) * 2^m

    // convert to int
    return (int)v.f;
}

它使用本维基百科文章中描述的算法.在我的机器上,它几乎是sqrt的两倍:)