sqrt (n)의 정수 부분을 얻는 가장 빠른 방법은 무엇입니까?

sqrt (n)의 정수 부분을 얻는 가장 빠른 방법은 무엇입니까?

---

우리가 알다시피 n완전한 제곱 sqrt(n)이 아니라면 정수가 아닐 것입니다. 정수 부분 만 필요하기 때문에 sqrt(n)분수 부분도 계산하는 데 시간이 걸리기 때문에 호출 이 그렇게 빠르지 않을 것이라고 생각합니다 .

그래서 제 질문은

의 실제 값을 계산하지 않고 sqrt (n) 의 정수 부분 만 얻을 수 있습니까 sqrt(n)? 알고리즘은 sqrt(n)( <math.h>또는 에서 정의 됨 <cmath>) 보다 빠릅니다 .

가능하다면 asm블록에 코드를 작성할 수도 있습니다.

---

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

예상대로 Fast 계산이 Int 계산 보다 훨씬 더 잘 수행 됩니다.

아, 그건 그렇고, sqrt더 빠릅니다 :)

---

편집 :이 대답은 어리석은 것입니다. (int) sqrt(i)

적절한 설정 ( -march=native -m64 -O3)으로 프로파일 링 한 후 위의 작업이 훨씬 빨라졌습니다.

---

알겠습니다. 약간 오래된 질문이지만 "가장 빠른"답변은 아직 제공되지 않았습니다. 가장 빠른 방법은 (내 생각에) 이 Embedded.com 기사 에서 자세히 설명하는 Binary Square Root 알고리즘 입니다.

기본적으로 이것으로 귀결됩니다.

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)데 2750ms가 걸렸습니다. 사용 (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

---

"빠른 정수 제곱근"을 검색하여 많은 옵션을 찾을 수 있다고 생각하지만 다음은 잘 작동 할 수있는 잠재적으로 새로운 아이디어입니다 (각각 독립적이거나 조합 할 수 있음).

1. static const지원하려는 도메인의 모든 완벽한 사각형 배열을 만들고 그것에 대해 빠른 분기없는 이진 검색을 수행합니다. 배열의 결과 인덱스는 제곱근입니다.
2. 숫자를 부동 소수점으로 변환하고 가수와 지수로 나눕니다. 지수를 반으로 나누고 가수에 마법 요소를 곱하십시오 (찾는 작업). 이것은 당신에게 매우 가까운 근사치를 줄 수있을 것입니다. 정확하지 않은 경우 조정하기위한 마지막 단계를 포함합니다 (또는 위의 이진 검색의 시작점으로 사용).

---

나는 빠른 계산의 가능한 많은 방법에 대해 논의한 Google search좋은 기사를 제공 Calculate an integer square root하고 좋은 참조 기사가 있다고 생각합니다. 여기에있는 어느 누구도 그들보다 더 나은 것을 제공 할 수 없다고 생각합니다 (그리고 누군가가 그것에 대해 먼저 논문을 만들 수 있다면). 그들과 모호한 부분이 있다면 우리가 당신을 잘 도울 수있을 것입니다.

---

근사치에 신경 쓰지 않는다면 내가 함께 만든이 정수 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;
}

It uses the algorithm described in this Wikipedia article. On my machine it's almost twice as fast as sqrt :)

---

To do integer sqrt you can use this specialization of newtons method:

Def isqrt(N):

    a = 1
    b = N

    while |a-b| > 1
        b = N / a
        a = (a + b) / 2

    return a

Basically for any x the sqrt lies in the range (x ... N/x), so we just bisect that interval at every loop for the new guess. Sort of like binary search but it converges must faster.

This converges in O(loglog(N)) which is very fast. It also doesn't use floating point at all, and it will also work well for arbitrary precision integers.

---

Why nobody suggests the quickest method?

If:

1. the range of numbers is limited
2. memory consumption is not crucial
3. application launch time is not critical

then create int[MAX_X] filled (on launch) with sqrt(x) (you don't need to use the function sqrt() for it).

All these conditions fit my program quite well. Particularly, an int[10000000] array is going to consume 40MB.

What's your thoughts on this?

---

This is so short that it 99% inlines:

static inline int sqrtn(int num) {
    int i;
    __asm__ (
        "pxor %%xmm0, %%xmm0\n\t"   // clean xmm0 for cvtsi2ss
        "cvtsi2ss %1, %%xmm0\n\t"   // convert num to float, put it to xmm0
        "sqrtss %%xmm0, %%xmm0\n\t" // square root xmm0
        "cvttss2si %%xmm0, %0"      // float to int
        :"=r"(i):"r"(num):"%xmm0"); // i: result, num: input, xmm0: scratch register
    return i;
}

Why clean xmm0? Documentation of cvtsi2ss

The destination operand is an XMM register. The result is stored in the low doubleword of the destination operand, and the upper three doublewords are left unchanged.

GCC Intrinsic version (runs only on GCC):

#include <xmmintrin.h>
int sqrtn2(int num) {
    register __v4sf xmm0 = {0, 0, 0, 0};
    xmm0 = __builtin_ia32_cvtsi2ss(xmm0, num);
    xmm0 = __builtin_ia32_sqrtss(xmm0);
    return __builtin_ia32_cvttss2si(xmm0);
}

Intel Intrinsic version (tested on GCC, Clang, ICC):

#include <xmmintrin.h>
int sqrtn2(int num) {
    register __m128 xmm0 = _mm_setzero_ps();
    xmm0 = _mm_cvt_si2ss(xmm0, num);
    xmm0 = _mm_sqrt_ss(xmm0);
    return _mm_cvtt_ss2si(xmm0);
}

^^^^ All of them require SSE 1 (not even SSE 2).

---

In many cases, even exact integer sqrt value is not needed, enough having good approximation of it. (For example, it often happens in DSP optimization, when 32-bit signal should be compressed to 16-bit, or 16-bit to 8-bit, without loosing much precision around zero).

I've found this useful equation:

k = ceil(MSB(n)/2); - MSB(n) is the most significant bit of "n"

sqrt(n) ~= 2^(k-2)+(2^(k-1))*n/(2^(2*k))); - all multiplications and divisions here are very DSP-friendly, as they are only 2^k.

This equation generates smooth curve (n, sqrt(n)), its values are not very much different from real sqrt(n) and thus can be useful when approximate accuracy is enough.

---

If you need performance on computing square root, I guess you will compute a lot of them. Then why not caching the answer? I don't know the range for N in your case, nor if you will compute many times the square root of the same integer, but if yes, then you can cache the result each time your method is called (in an array would be the most efficient if not too large).

---

On my computer with gcc, with -ffast-math, converting a 32-bit integer to float and using sqrtf takes 1.2 s per 10^9 ops (without -ffast-math it takes 3.54 s).

The following algorithm uses 0.87 s per 10^9 at the expense of some accuracy: errors can be as much as -7 or +1 although the RMS error is only 0.79:

uint16_t SQRTTAB[65536];

inline uint16_t approxsqrt(uint32_t x) { 
  const uint32_t m1 = 0xff000000;
  const uint32_t m2 = 0x00ff0000;
  if (x&m1) {
    return SQRTTAB[x>>16];
  } else if (x&m2) {
    return SQRTTAB[x>>8]>>4;
  } else {
    return SQRTTAB[x]>>8;
  }
}

The table is constructed using:

void maketable() {
  for (int x=0; x<65536; x++) {
    double v = x/65535.0;
    v = sqrt(v);
    int y = int(v*65535.0+0.999);
    SQRTTAB[x] = y;
  }
}

I found that refining the bisection using further if statements does improve accuracy, but it also slows things down to the point that sqrtf is faster, at least with -ffast-math.

참고URL : https://stackoverflow.com/questions/4930307/fastest-way-to-get-the-integer-part-of-sqrtn