¿Es posible rodar una versión mucho más rápido de sqrt

Question

En una aplicación que estoy perfiles, he encontrado que en algunos escenarios de esta función es capaz de hacerse cargo del 10% del tiempo total de ejecución.

he visto discusión sobre los años de implementaciones más rápido sqrt utilizando astuto engaño de punto flotante, pero no sé si estas cosas son anticuadas en las CPU modernas.

MSVC ++ 2008 se está utilizando el compilador, para referencia ... aunque yo supongo sqrt no va a añadir mucho sobrecarga sin embargo.

Ver también aquí para discusión similar en modf función.

EDIT: para referencia, este es un método ampliamente utilizado, pero es lo que realmente ¿mucho mas rapido? ¿Cuántos ciclos es SQRT todos modos en estos días?

Answer 1

Sí, es posible incluso sin engaños:

1) Precisión de sacrificio para la velocidad:. El algoritmo sqrt es iterativo, re-implementar con menos iteraciones

2) tablas de búsqueda:. O bien sólo para el punto de inicio de la iteración, o combinados con interpolación para conseguir que todo el camino

3) el almacenamiento en caché: están siempre sqrting el mismo conjunto limitado de valores? si es así, el almacenamiento en caché puede funcionar bien. He encontrado este útil en aplicaciones gráficas, donde se está calculando la misma cosa para un montón de formas del mismo tamaño, por lo que los resultados pueden almacenar en caché de manera útil.

Answer 2

Hay una gran tabla de comparación aquí: http://assemblyrequired.crashworks.org/timing-square-root/

Es una larga historia corta, ssqrts de SSE2 es aproximadamente 2 veces más rápido que la FPU fsqrt, y una aproximación + iteración es de aproximadamente 4 veces más rápido que eso (8x general).

Además, si usted está tratando de tomar una raíz cuadrada de precisión simple, asegúrese de que en realidad lo que está recibiendo. He oído hablar de por lo menos un compilador que convertiría el argumento de flotador a un doble, llama raíz cuadrada de doble precisión, a continuación, convertir de nuevo a flotador.

Answer 3

Usted es muy probable que ganar más mejoras en la velocidad cambiando algoritmos que cambiando su implementaciones : Trata de llamada sqrt() menos en lugar de hacer llamadas más rápidamente. (Y si usted piensa que esto no es posible - las mejoras para sqrt() que mencionas son sólo eso: la mejora de la algoritmo se utiliza para calcular una raíz cuadrada.)

Dado que se utiliza muy a menudo, lo más probable es que la implementación de la biblioteca estándar de sqrt() es óptima casi para el caso general. A menos que tenga un dominio restringido (por ejemplo, si necesita menos precisión) en el que el algoritmo puede tomar algunos atajos, es muy poco probable que alguien se le ocurre una aplicación que es más rápido.

Tenga en cuenta que, puesto que la función utiliza el 10% de su tiempo de ejecución, incluso si logra llegar a una aplicación que sólo toma 75% del tiempo de std::sqrt(), esto sigue siendo sólo hará que su tiempo de ejecución por 2,5% . Para la mayoría de las aplicaciones de los usuarios no se daría cuenta de esto, salvo que utilizan un reloj para medir.

Answer 4

¿Qué tan preciso que necesita su sqrt ser? Usted puede obtener una aproximación razonable muy rápidamente: ver de Quake3 excelente cuadrado inverso función raíz para la inspiración (nota que el código es GPL, por lo que no puede querer integrarlo directamente).

Answer 5

No sé si te fijas esto, pero he leído sobre él antes, y parece que la cosa más rápida de hacerlo es reemplazar la función sqrt con una versión en línea de montaje;

se puede ver una descripción de una carga de alternativas aquí .

Lo mejor es este fragmento de la magia:

double inline __declspec (naked) __fastcall sqrt(double n)
{
    _asm fld qword ptr [esp+4]
    _asm fsqrt
    _asm ret 8
} 

Se trata de 4.7x más rápido que la llamada sqrt estándar con la misma precisión.

Answer 6

Esta es una manera rápida con una tabla de consulta de tan sólo 8 KB. Mistake es de ~ 0,5% del resultado. Usted puede ampliar fácilmente la mesa, lo que reduce el error. Corre alrededor de 5 veces más rápido que la raíz cuadrada regular ()

// LUT for fast sqrt of floats. Table will be consist of 2 parts, half for sqrt(X) and half for sqrt(2X).
const int nBitsForSQRTprecision = 11;                       // Use only 11 most sagnificant bits from the 23 of float. We can use 15 bits instead. It will produce less error but take more place in a memory. 
const int nUnusedBits   = 23 - nBitsForSQRTprecision;       // Amount of bits we will disregard
const int tableSize     = (1 << (nBitsForSQRTprecision+1)); // 2^nBits*2 because we have 2 halves of the table.
static short sqrtTab[tableSize]; 
static unsigned char is_sqrttab_initialized = FALSE;        // Once initialized will be true

// Table of precalculated sqrt() for future fast calculation. Approximates the exact with an error of about 0.5%
// Note: To access the bits of a float in C quickly we must misuse pointers.
// More info in: http://en.wikipedia.org/wiki/Single_precision
void build_fsqrt_table(void){
    unsigned short i;
    float f;
    UINT32 *fi = (UINT32*)&f;

    if (is_sqrttab_initialized)
        return;

    const int halfTableSize = (tableSize>>1);
    for (i=0; i < halfTableSize; i++){
         *fi = 0;
         *fi = (i << nUnusedBits) | (127 << 23); // Build a float with the bit pattern i as mantissa, and an exponent of 0, stored as 127

         // Take the square root then strip the first 'nBitsForSQRTprecision' bits of the mantissa into the table
         f = sqrtf(f);
         sqrtTab[i] = (short)((*fi & 0x7fffff) >> nUnusedBits);

         // Repeat the process, this time with an exponent of 1, stored as 128
         *fi = 0;
         *fi = (i << nUnusedBits) | (128 << 23);
         f = sqrtf(f);
         sqrtTab[i+halfTableSize] = (short)((*fi & 0x7fffff) >> nUnusedBits);
    }
    is_sqrttab_initialized = TRUE;
}

// Calculation of a square root. Divide the exponent of float by 2 and sqrt() its mantissa using the precalculated table.
float fast_float_sqrt(float n){
    if (n <= 0.f) 
        return 0.f;                           // On 0 or negative return 0.
    UINT32 *num = (UINT32*)&n;
    short e;                                  // Exponent
    e = (*num >> 23) - 127;                   // In 'float' the exponent is stored with 127 added.
    *num &= 0x7fffff;                         // leave only the mantissa 

    // If the exponent is odd so we have to look it up in the second half of the lookup table, so we set the high bit.
    const int halfTableSize = (tableSize>>1);
    const int secondHalphTableIdBit = halfTableSize << nUnusedBits;
    if (e & 0x01) 
        *num |= secondHalphTableIdBit;  
    e >>= 1;                                  // Divide the exponent by two (note that in C the shift operators are sign preserving for signed operands

    // Do the table lookup, based on the quaternary mantissa, then reconstruct the result back into a float
    *num = ((sqrtTab[*num >> nUnusedBits]) << nUnusedBits) | ((e + 127) << 23);
    return n;
}