Самый быстрый способ получить целую часть sqrt (n)?

Ответ 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 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 быстрее:)

Ответ 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) заняло 3600 мс. Это было сделано с помощью g++ -O3. Используя параметр компиляции -ffast-math, время было 2100 мс и 3100 мс соответственно. Обратите внимание, что это не использует даже одну строку ассемблера, поэтому, возможно, она будет намного быстрее.

Вышеприведенный код работает как для 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

Ответ 3

Хотя я подозреваю, что вы можете найти множество опций, ища "быстрый целочисленный квадратный корень", вот некоторые потенциально новые идеи, которые могут работать хорошо (каждый независимый или, возможно, вы можете их комбинировать):

• Создайте массив static const всех совершенных квадратов в домене, который вы хотите поддерживать, и выполните быстрый разветвленный двоичный поиск на нем. Получившийся индекс в массиве - это квадратный корень.
• Преобразование числа в плавающую точку и разбиение его на мантиссы и экспонента. Уменьшите экспоненту и умножьте мантиссу на какой-то магический фактор (ваша работа, чтобы найти его). Это должно дать вам очень близкое приближение. Включите последний шаг, чтобы настроить его, если он не является точным (или использовать его как отправную точку для двоичного поиска выше).

Ответ 4

Я думаю, Google search содержит хорошие статьи, такие как Calculate an integer square root, в котором обсуждалось слишком много возможных способов быстрого расчета, и есть хорошие справочные статьи, я думаю, что никто здесь не может обеспечить их лучше (и если кто-то может сначала написать об этом бумагу), но если вы их прочтете и там есть двусмысленность с ними, тогда мы можем помочь вам хорошо.

Ответ 5

Если вы не против приближения, как насчет этой целочисленной функции 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;
}

Он использует алгоритм, описанный в этой статье Wikipedia. На моей машине это почти в два раза быстрее, чем sqrt:)

Ответ 6

Чтобы сделать integer sqrt, вы можете использовать эту специализацию метода newtons:

Def isqrt(N):

    a = 1
    b = N

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

    return a

В принципе для любого x sqrt лежит в диапазоне (x... N/x), поэтому мы просто делим пополам этот интервал в каждом цикле для нового предположения. Подобный бинарный поиск, но он сходится быстрее.

Это сходится в O (loglog (N)), что очень быстро. Он также не использует плавающие точки вообще, и он также будет хорошо работать для целых чисел точности.

Ответ 7

Почему никто не предлагает самый быстрый метод?

Если:

• диапазон чисел ограничен.
• Потребление памяти не имеет решающего значения.
• Время запуска приложения не критично.

затем создайте int[MAX_X], заполненный (при запуске) с помощью sqrt(x) (вам не нужно использовать для него функцию sqrt()).

Все эти условия вполне соответствуют моей программе. В частности, массив int[10000000] будет потреблять 40MB.

Что вы думаете об этом?

Ответ 8

Это так коротко, что в нем 99% строк:

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;
}

Зачем чистить xmm0? Документация cvtsi2ss

Операндом-адресатом является регистр XMM. Результат сохраняется в нижнем двойном слове операнда назначения, а три верхних двойных слова остаются без изменений.

Внутренняя версия GCC (работает только на 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 (протестирована на 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);
}

^^^^ Все они требуют SSE 1 (даже не SSE 2).

Ответ 9

Во многих случаях даже точное целочисленное значение sqrt не требуется, достаточно иметь хорошую аппроксимацию. (Например, это часто происходит при оптимизации DSP, когда 32-разрядный сигнал должен быть сжат до 16 бит или от 16 бит до 8 бит, без потери значительной точности вокруг нуля).

Я нашел это полезное уравнение:

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.

Это уравнение генерирует гладкую кривую (n, sqrt (n)), ее значения не сильно отличаются от реального sqrt (n) и поэтому могут быть полезны, когда достаточно приблизительной точности.

Ответ 10

Если вам нужна производительность при вычислении квадратного корня, я думаю, вы вычислите их много. Тогда почему бы не кешировать ответ? Я не знаю диапазон для N в вашем случае, и если вы будете много раз вычислять квадратный корень из того же целого числа, но если да, то вы можете кэшировать результат каждый раз, когда вы вызываете ваш метод (в массиве будет самый эффективный, если не слишком большой).

Ответ 11

На моем компьютере с gcc, с -ffast-math, преобразование 32-разрядного целого числа в float и использование sqrtf занимает 1,2 с на 10 ^ 9 операций (без -ffast-math требуется 3,54 с).

Следующий алгоритм использует 0,87 с на 10 ^ 9 за счет некоторой точности: ошибки могут достигать -7 или +1, хотя ошибка RMS составляет всего 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;
  }
}

Таблица построена с использованием:

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;
  }
}

Я обнаружил, что уточнение деления пополам с использованием дальнейшего, если утверждения повышают точность, но также замедляет работу до такой степени, что sqrtf быстрее, по крайней мере, с -ffast-math.