raices - raiz cuadrada en c++

¿Es posible rodar una versión significativamente más rápida de sqrt?

Question

En una aplicación que estoy perfilando, encontré que en algunos escenarios esta función puede tomar más del 10% del tiempo total de ejecución.

He visto discusiones a lo largo de los años sobre implementaciones de sqrt más rápidas usando trucos astutos de coma flotante, pero no sé si tales cosas están desactualizadas en las CPU modernas.

El compilador de MSVC ++ 2008 está siendo utilizado, como referencia ... aunque supongo que sqrt no va a agregar mucha sobrecarga.

Ver también aquí para una discusión similar sobre la función modf .

EDITAR: como referencia, this es un método ampliamente utilizado, pero ¿es realmente mucho más rápido? ¿Cuántos ciclos tiene SQRT de todos modos en estos días?

Answer 1

¿Qué tan preciso necesita que sea su sqrt ? Puede obtener aproximaciones razonables muy rápidamente: consulte la excelente función de raíz cuadrada inversa de Quake3 para obtener inspiración (tenga en cuenta que el código tiene GPL, por lo que es posible que no desee integrarlo directamente).

Answer 2

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

Para abreviar, los ssqrts de SSE2 son aproximadamente 2 veces más rápidos que FPU fsqrt, y una aproximación + iteración es aproximadamente 4 veces más rápida que eso (8x en general).

Además, si intenta tomar un sqrt de precisión simple, asegúrese de que eso sea lo que realmente obtiene. He oído hablar de al menos un compilador que convertiría el argumento flotante en un doble, llamar a sqrt de precisión doble y luego volver a convertir en flotante.

Answer 3

Aquí hay una manera rápida con una tabla de búsqueda de solo 8 KB. El error es ~ 0.5% del resultado. Puede ampliar fácilmente la tabla, reduciendo así el error. Se ejecuta aproximadamente 5 veces más rápido que el sqrt 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;
}

Answer 4

Es muy probable que ganes más mejoras de velocidad cambiando tus algoritmos que cambiando sus implementaciones : intenta llamar a sqrt() menos en lugar de hacer llamadas más rápido. (Y si cree que esto no es posible, las mejoras para sqrt() que menciona son solo eso: mejoras del algoritmo utilizado para calcular una raíz cuadrada).

Dado que se usa con mucha frecuencia, es probable que la implementación de sqrt() su biblioteca estándar sea casi óptima para el caso general. A menos que tenga un dominio restringido (por ejemplo, si necesita menos precisión) donde el algoritmo puede tomar algunos accesos directos, es muy poco probable que a alguien se le ocurra una implementación que sea más rápida.

Tenga en cuenta que, dado que esa función usa el 10% de su tiempo de ejecución, incluso si logra crear una implementación que solo toma el 75% del tiempo de std::sqrt() , esto aún reducirá el tiempo de ejecución. 2,5% . Para la mayoría de las aplicaciones, los usuarios ni siquiera notarían esto, excepto si usan un reloj para medir.

Answer 5

No sé si solucionó esto, pero ya lo leí antes, y parece que lo más rápido es reemplazar la función sqrt por una versión de ensamblaje en línea;

puedes ver una descripción de una carga de alternativas here .

Lo mejor es este fragmento de magia:

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

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

Answer 6

Sí, es posible incluso sin engaños:

1) precisión de sacrificio para la velocidad: el algoritmo sqrt es iterativo, se vuelve a implementar con menos iteraciones.

2) tablas de búsqueda: ya sea para el punto de inicio de la iteración, o combinado con la interpolación para llegar hasta allí.

3) almacenamiento en caché: ¿siempre estás cuadrando el mismo conjunto limitado de valores? si es así, el almacenamiento en caché puede funcionar bien. Lo encontré útil en aplicaciones gráficas donde lo mismo se calcula para muchas formas del mismo tamaño, por lo que los resultados se pueden almacenar en caché.