获得sqrt(n)的整数部分的最快方法是什么?

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

所以我的问题是,

我们可以得到唯一的sqrt(n)的整数部分,而不计算的实际值sqrt(n) 算法应该比sqrt(n)更快(在<math.h><cmath> )?

如果可能的话,你也可以在asm块中编写代码。


我会尝试快速反向平方根技巧。

这是一种非常好的近似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 iterationsn";
    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更快:)


编辑:这个答案是愚蠢的 - 使用(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


虽然我怀疑你可以通过搜索“快速整数平方根”找到很多选项,但下面是一些可能运作良好的新想法(每个独立的,或者你可以合并它们):

  • 在要支持的域中创建一个包含所有完美正方形的static const数组,并在其上执行快速无分叉二进制搜索。 数组中的结果索引是平方根。
  • 将数字转换为浮点数并将其分解为尾数和指数。 将指数减半并将尾数乘以一些魔法因子(您的工作找到它)。 这应该能够给你一个非常接近的近似值。 包括最后一步,如果它不精确(或将其用作上述二进制搜索的起点),则对其进行调整。
  • 链接地址: http://www.djcxy.com/p/86601.html

    上一篇: Fastest way to get the integer part of sqrt(n)?

    下一篇: Java `final` method: what does it promise?