Самый быстрый способ получить целую часть 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

Где, как и ожидалось, вычисления Fast работают намного лучше, чем вычисления Int .

Да, и, кстати, sqrt быстрее:)

Answer 2

Изменить: этот ответ глупый - используйте (int) sqrt(i)

После профилирования с правильными настройками (-march=native -m64 -O3) вышеописанное было лотом быстрее.

---

Хорошо, немного старый вопрос, но самый быстрый ответ еще не дан.Самым быстрым (я думаю) является алгоритм Binary Square Root, полностью объясненный в этой статье 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, так и для 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

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

Answer 4

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

Answer 5

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

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

Answer 6

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

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)), который очень быстро.Он также вообще не использует числа с плавающей запятой и будет хорошо работать для целых чисел произвольной точности.

Answer 7

Это так коротко, что в нем 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).

Answer 8

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

Если:

1. диапазон номеров ограничен
2. потребление памяти не критично
3. время запуска приложения не критично

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

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

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

Answer 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) и, следовательно, может быть полезным, когда приблизительная точность достаточна.

Answer 10

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

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

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