c++ - raiz - ¿La forma más rápida de obtener la parte entera de sqrt(n)?

raiz enesima en c++

Question

Como sabemos que si n no es un cuadrado perfecto, entonces sqrt(n) no sería un número entero. Como solo necesito la parte entera, creo que llamar a sqrt(n) no sería tan rápido, ya que también lleva tiempo calcular la parte fraccionaria.

Entonces mi pregunta es

¿Podemos obtener solo la parte entera de sqrt (n) sin calcular el valor real de sqrt(n) ? El algoritmo debe ser más rápido que sqrt(n) (definido en <math.h> o <cmath> )?

Si es posible, también puede escribir el código en el bloque asm .

Answer 1

Editar: esta respuesta es tonta - use (int) sqrt(i)

Después de -march=native -m64 -O3 perfiles con la configuración adecuada ( -march=native -m64 -O3 ), lo anterior fue mucho más rápido.

Muy bien, es una pregunta un poco vieja, pero la respuesta "más rápida" aún no se ha dado. El más rápido (creo) es el algoritmo Binary Square Root, explicado completamente en este artículo de Embedded.com .

Básicamente se trata de esto:

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

En mi máquina (Q6600, Ubuntu 10.10) hice un perfil tomando la raíz cuadrada de los números 1-100000000. Usar iqsrt(i) tomó 2750 ms. Usar (unsigned short) sqrt((float) i) tomó 3600ms. Esto se hizo usando g++ -O3 . Usando la opción de compilación -ffast-math los tiempos fueron 2100ms y 3100ms respectivamente. Tenga en cuenta que esto es sin usar ni siquiera una sola línea de ensamblador por lo que probablemente aún podría ser mucho más rápido.

El código anterior funciona tanto para C como para C ++ y con cambios de sintaxis menores también para Java.

Lo que funciona aún mejor para un rango limitado es una búsqueda binaria. En mi máquina, esto hace que la versión anterior salga del agua por un factor de 4. Por desgracia, tiene un alcance muy limitado:

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

Una versión de 32 bits se puede descargar aquí: https://gist.github.com/3481770

Answer 2

¿Por qué nadie sugiere el método más rápido?

Si:

1. el rango de números es limitado
2. el consumo de memoria no es crucial
3. el tiempo de lanzamiento de la aplicación no es crítico

luego cree int[MAX_X] lleno (en el inicio) con sqrt(x) (no necesita usar la función sqrt() para ello).

Todas estas condiciones se ajustan bastante bien a mi programa. Particularmente, una matriz int[10000000] va a consumir 40MB .

¿Qué piensas sobre esto?

Answer 3

Creo que la Google search proporciona buenos artículos como Calculate an integer square root que discutió sobre muchas formas posibles de cálculo rápido y hay buenos artículos de referencia, creo que nadie aquí puede proporcionar mejores que ellos (y si alguien primero puede producir un artículo sobre ), pero si los lees y hay ambigüedad con ellos, entonces podemos ser capaces de ayudarte bien.

Answer 4

En mi computadora con gcc, con -ffast-math, convertir un entero de 32 bits para flotar y usar sqrtf toma 1.2 s por 10 ^ 9 ops (sin -ffast-math toma 3.54 s).

El siguiente algoritmo usa 0.87 s por 10 ^ 9 a expensas de cierta precisión: los errores pueden ser hasta -7 o +1 aunque el error RMS es solo 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;
  }
}

La tabla se construye usando:

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

Descubrí que refinar la bisección usando declaraciones if si mejora la precisión, pero también ralentiza las cosas hasta el punto de que sqrtf es más rápido, al menos con -ffast-math.

Answer 5

En muchos casos, incluso el valor exacto de enteros enteros no es necesario, ya que tiene una buena aproximación. (Por ejemplo, a menudo ocurre en la optimización DSP, cuando la señal de 32 bits debe comprimirse a 16 bits, o de 16 bits a 8 bits, sin perder mucha precisión alrededor de cero).

Encontré esta ecuación útil:

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.

Esta ecuación genera una curva suave (n, sqrt (n)), sus valores no son muy diferentes de los reales sqrt (n) y, por lo tanto, pueden ser útiles cuando la precisión aproximada es suficiente.

Answer 6

Esto es tan corto que tiene un 99% de líneas.

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

¿Por qué limpiar xmm0 ? Documentación de cvtsi2ss

El operando de destino es un registro XMM. El resultado se almacena en la doble palabra baja del operando de destino, y las tres palabras dobles superiores se dejan sin cambios.

Versión intrínseca de GCC (solo se ejecuta en 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);
}

Versión intrínseca de Intel (probada en 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_cvt_ss2si(xmm0);
}

^^^^ Todos ellos requieren SSE 1. (ni siquiera SSE 2)

Answer 7

Intentaría el truco Fast Inverse Square Root .

Es una forma de obtener una muy buena aproximación de 1/sqrt(n) sin ninguna rama, basada en algunos giros de bits, por lo que no es portátil (especialmente entre las plataformas de 32 y 64 bits).

Una vez que lo obtiene, solo necesita invertir el resultado y toma la parte entera.

Puede haber trucos más rápidos, por supuesto, ya que este es un poco redondo.

EDIT : ¡hagámoslo!

Primero un pequeño ayudante:

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

Entonces el cuerpo principal:

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

Y los resultados:

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

Donde como se espera, el cálculo rápido funciona mucho mejor que el cálculo Int .

Ah, y por cierto, sqrt es más rápido :)

Answer 8

Para hacer enteros sqrt puedes usar esta especialización del método newtons:

Def isqrt(N):

    a = 1
    b = N

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

    return a

Básicamente para cualquier x el sqrt se encuentra en el rango (x ... N / x), por lo que solo bisectamos ese intervalo en cada ciclo para la nueva suposición. Algo así como la búsqueda binaria pero converge debe ser más rápida.

Esto converge en O (loglog (N)) que es muy rápido. Tampoco usa punto flotante, y también funcionará bien para enteros de precisión arbitrarios.

Answer 9

Si bien sospecho que puede encontrar muchas opciones buscando "raíz cuadrada entera rápida", aquí hay algunas ideas potencialmente nuevas que pueden funcionar bien (cada una independiente, o tal vez puede combinarlas):

1. Crea una matriz de static const de todos los cuadrados perfectos en el dominio que deseas admitir, y realiza una búsqueda binaria rápida sin ramas en él. El índice resultante en la matriz es la raíz cuadrada.
2. Convierta el número a punto flotante y divídalo en mantisa y exponente. Reduce a la mitad el exponente y multiplica la mantisa por algún factor mágico (tu trabajo para encontrarlo). Esto debería poder darle una aproximación muy cercana. Incluya un paso final para ajustarlo si no es exacto (o úselo como punto de partida para la búsqueda binaria anterior).

Answer 10

Si necesita rendimiento en la computación de la raíz cuadrada, supongo que calculará muchos de ellos. Entonces, ¿por qué no almacenar en caché la respuesta? No conozco el rango para N en su caso, ni si calculará muchas veces la raíz cuadrada del mismo entero, pero si es así, entonces puede almacenar en caché el resultado cada vez que se llame a su método (en una matriz sería el más eficiente si no demasiado grande).

Answer 11

Si no te importa una aproximación, ¿qué hay de esta función entera de sqrt que improvisé?

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

Utiliza el algoritmo descrito en este artículo de Wikipedia . En mi máquina es casi dos veces más rápido que sqrt :)