Raíz cuadrada inversa rápida , a veces denominada Fast InvSqrt() o por la constante hexadecimal 0x5F3759DF , es un algoritmo que estima , el recíproco (o inverso multiplicativo) de la raíz cuadrada de un número de punto flotante de 32 bits en formato de punto flotante IEEE 754. El algoritmo es más conocido por su implementación en 1999 en Quake III Arena , un videojuego de disparos en primera persona basado en gran medida en gráficos 3D . Con los avances de hardware posteriores, especialmente la instrucción x86 SSE , este algoritmo no es generalmente la mejor opción para las computadoras modernas, [1] aunque sigue siendo un ejemplo histórico interesante. [2] rsqrtss
El algoritmo acepta un número de punto flotante de 32 bits como entrada y almacena un valor reducido a la mitad para su uso posterior. Luego, al tratar los bits que representan el número de punto flotante como un entero de 32 bits, se realiza un desplazamiento lógico a la derecha de un bit y el resultado se resta del número 0x 5F3759DF , que es una representación de punto flotante de una aproximación de . [3] Esto da como resultado la primera aproximación de la raíz cuadrada inversa de la entrada. Al tratar los bits nuevamente como un número de punto flotante, ejecuta una iteración del método de Newton , lo que produce una aproximación más precisa.
William Kahan y KC Ng, de Berkeley, escribieron un artículo inédito en mayo de 1986 en el que describían cómo calcular la raíz cuadrada utilizando técnicas de manipulación de bits seguidas de iteraciones de Newton. [4] A finales de los años 80, Cleve Moler, de Ardent Computer, se enteró de esta técnica [5] y se la transmitió a su compañero de trabajo Greg Walsh. Greg Walsh ideó el ahora famoso algoritmo de raíz cuadrada inversa constante y rápida. Gary Tarolli trabajaba como consultor para Kubota, la empresa que financiaba a Ardent en ese momento, y probablemente llevó el algoritmo a 3dfx Interactive alrededor de 1994. [6] [7]
Jim Blinn demostró una aproximación simple de la raíz cuadrada inversa en una columna de 1997 para IEEE Computer Graphics and Applications . [8] La ingeniería inversa de otros videojuegos 3D contemporáneos descubrió una variación del algoritmo en Interstate '76 de Activision de 1997. [9]
Quake III Arena , unvideojuego de disparos en primera persona , fue lanzado en 1999 por id Software y utilizó el algoritmo. Brian Hook pudo haber traído el algoritmo de 3dfx a id Software. [6] Una discusión del código apareció en el foro de desarrolladores chino CSDN en 2000, [10] y Usenet y el foro gamedev.net difundieron ampliamente el código en 2002 y 2003. [11] Surgió la especulación sobre quién escribió el algoritmo y cómo se derivó la constante; algunos supusieron que John Carmack . [7] El código fuente completo de Quake III fue lanzado en QuakeCon 2005 , pero no proporcionó respuestas. La cuestión de la autoría se resolvió en 2006 cuando Greg Walsh, el autor original, contactó a Beyond3D después de que su especulación ganara popularidad en Slashdot. [6]
En 2007, el algoritmo se implementó en algunos sombreadores de vértices de hardware dedicados utilizando matrices de puertas programables en campo (FPGA). [12] [13]
La raíz cuadrada inversa de un número de punto flotante se utiliza en el procesamiento de señales digitales para normalizar un vector, escalándolo a una longitud de 1 para producir un vector unitario . [14] Por ejemplo, los programas de gráficos de computadora utilizan raíces cuadradas inversas para calcular ángulos de incidencia y reflexión para iluminación y sombreado . Los programas de gráficos 3D deben realizar millones de estos cálculos cada segundo para simular la iluminación. Cuando se desarrolló el código a principios de la década de 1990, la mayor parte de la potencia de procesamiento de punto flotante estaba por detrás de la velocidad del procesamiento de números enteros. [7] Esto era problemático para los programas de gráficos 3D antes de la llegada del hardware especializado para manejar la transformación y la iluminación . El cálculo de raíces cuadradas generalmente depende de muchas operaciones de división, que para los números de punto flotante son computacionalmente costosas . El cuadrado inverso rápido genera una buena aproximación con solo un paso de división.
La longitud del vector se determina calculando su norma euclidiana : la raíz cuadrada de la suma de los cuadrados de los componentes del vector . Cuando cada componente del vector se divide por esa longitud, el nuevo vector será un vector unitario que apunta en la misma dirección. En un programa de gráficos 3D, todos los vectores están en el espacio tridimensional , por lo que sería un vector .
es la norma euclidiana del vector.
es el vector normalizado (unitario). Sustituyendo , la ecuación también se puede escribir como:
que relaciona el vector unitario con la raíz cuadrada inversa de los componentes de la distancia. La raíz cuadrada inversa se puede utilizar para calcular porque esta ecuación es equivalente a
donde el término fraccionario es la raíz cuadrada inversa de .
En ese momento, la división de punto flotante era generalmente costosa en comparación con la multiplicación; el algoritmo de raíz cuadrada inversa rápida omitió el paso de división, lo que le dio su ventaja de rendimiento.
El código siguiente es la implementación rápida de la raíz cuadrada inversa de Quake III Arena , despojada de las directivas del preprocesador C , pero incluyendo el texto exacto del comentario original: [15]
float Q_rsqrt ( float número ) { long i ; float x2 , y ; const float tres mitades = 1.5F ; x2 = número * 0.5F ; y = número ; i = * ( long * ) & y ; // hackeo malvado de nivel de bit de punto flotante i = 0x5f3759df - ( i >> 1 ); // ¿qué carajo? y = * ( float * ) & i ; y = y * ( threehalfs - ( x2 * y * y ) ); // 1.ª iteración // y = y * (threehalfs - (x2 * y * y) ); // 2.ª iteración, esto se puede eliminar devuelve y ; }
En ese momento, el método general para calcular la raíz cuadrada inversa era calcular una aproximación para , luego revisar esa aproximación mediante otro método hasta que estuviera dentro de un rango de error aceptable del resultado real. Los métodos de software comunes a principios de la década de 1990 extraían aproximaciones de una tabla de búsqueda . [16] La clave de la raíz cuadrada inversa rápida era calcular directamente una aproximación utilizando la estructura de números de punto flotante, lo que resultó más rápido que las búsquedas en tablas. El algoritmo era aproximadamente cuatro veces más rápido que calcular la raíz cuadrada con otro método y calcular el recíproco mediante una división de punto flotante. [17] El algoritmo fue diseñado teniendo en cuenta la especificación de punto flotante de 32 bits IEEE 754-1985 , pero la investigación de Chris Lomont mostró que podría implementarse en otras especificaciones de punto flotante. [18]
Las ventajas en velocidad ofrecidas por el truco de la raíz cuadrada inversa rápida vinieron de tratar la palabra de punto flotante de 32 bits [nota 1] como un entero , luego restarlo de una constante " mágica ", 0x 5F3759DF . [7] [19] [20] [21] Esta resta de enteros y desplazamiento de bits da como resultado un patrón de bits que, cuando se redefine como un número de punto flotante, es una aproximación aproximada para la raíz cuadrada inversa del número. Se realiza una iteración del método de Newton para ganar algo de precisión y el código está terminado. El algoritmo genera resultados razonablemente precisos utilizando una primera aproximación única para el método de Newton ; sin embargo, es mucho más lento y menos preciso que usar la instrucción SSE en procesadores x86 también lanzados en 1999. [1] [22]rsqrtss
Como ejemplo, se puede utilizar el número para calcular . Los primeros pasos del algoritmo se ilustran a continuación:
0011_1110_0010_0000_0000_0000_0000_0000 Patrón de bits de x e i0001_1111_0001_0000_0000_0000_0000_0000 Desplazar una posición a la derecha: (i >> 1)0101_1111_0011_0111_0101_1001_1101_1111 El número mágico 0x5F3759DF0100_0000_0010_0111_0101_1001_1101_1111 El resultado de 0x5F3759DF - (i >> 1)
Interpretando como representación IEEE de 32 bits:
0_01111100_010000000000000000000000 1,25 × 2 −3 0_00111110_001000000000000000000000 1,125 × 2 −65 0_10111110_01101110101100111011111 1.432430... × 2 63 0_10000000_01001110101100111011111 1.307430... × 2 1
Reinterpretando este último patrón de bits como un número de punto flotante se obtiene la aproximación , que tiene un error de aproximadamente el 3,4 %. Después de una única iteración del método de Newton , el resultado final es , un error de solo el 0,17 %.
Según el estándar C, no es válido reinterpretar un valor de punto flotante como un entero mediante la conversión y luego desreferenciando el puntero a él ( comportamiento indefinido ). Otra forma sería colocar el valor de punto flotante en una unión anónima que contenga un miembro entero sin signo adicional de 32 bits, y los accesos a ese entero proporcionan una vista a nivel de bits del contenido del valor de punto flotante. Sin embargo, el juego de palabras de tipos a través de una unión también es un comportamiento indefinido en C++.
# incluir <stdint.h> // uint32_t float Q_rsqrt ( float número ) { unión { float f ; uint32_t i ; } conv = { . f = número }; conv . i = 0x5f3759df - ( conv . i >> 1 ); conv . f *= 1.5F - ( número * 0.5F * conv . f * conv . f ); devolver conv . f ; }
En C++ moderno, el método recomendado para implementar las conversiones de esta función es mediante C++20 std::bit_cast
. Esto también permite que la función funcione en constexpr
contexto:
importar std ; constexpr std :: float32_t Q_rsqrt ( std :: float32_t número ) noexcept { const auto y = std :: bit_cast < std :: float32_t > ( 0x5f3759df - ( std :: bit_cast < std :: uint32_t > ( número ) >> 1 )); return y * ( 1.5f 32 - ( número * 0.5f 32 * y * y )); }
El algoritmo se calcula realizando los siguientes pasos:
Dado que este algoritmo depende en gran medida de la representación a nivel de bits de números de punto flotante de precisión simple, aquí se ofrece una breve descripción general de esta representación. Para codificar un número real distinto de cero como un número flotante de precisión simple, el primer paso es escribirlo como un número binario normalizado : [23]
donde el exponente es un entero y es la representación binaria de la mantisa . Dado que el bit único antes del punto en la mantisa es siempre 1, no es necesario almacenarlo. La ecuación se puede reescribir como:
donde significa , por lo tanto . A partir de esta forma, se calculan tres enteros sin signo: [24]
Estos campos se empaquetan luego, de izquierda a derecha, en un contenedor de 32 bits. [25]
A modo de ejemplo, consideremos nuevamente el número . Al normalizarlo obtenemos:
y por lo tanto, los tres campos enteros sin signo son:
Estos campos se empaquetan como se muestra en la siguiente figura:
Si se tuviera que calcular sin ordenador ni calculadora, sería útil una tabla de logaritmos , junto con la identidad , que es válida para cualquier base . La raíz cuadrada inversa rápida se basa en esta identidad y en el hecho de que al convertir un float32 en un entero se obtiene una aproximación aproximada de su logaritmo. Aquí se explica cómo:
Si es un número normal positivo :
entonces
y dado que , el logaritmo en el lado derecho se puede aproximar mediante [26]
donde es un parámetro libre que se utiliza para ajustar la aproximación. Por ejemplo, produce resultados exactos en ambos extremos del intervalo, mientras que produce la aproximación óptima (la mejor en el sentido de la norma uniforme del error). Sin embargo, el algoritmo no utiliza este valor ya que no tiene en cuenta los pasos posteriores.
Así pues, existe la aproximación
Al interpretar el patrón de bits de punto flotante de como un entero se obtiene [nota 4]
Entonces parece que es una aproximación lineal por partes escalada y desplazada de , como se ilustra en la figura de la derecha. En otras palabras, se aproxima mediante
El cálculo de se basa en la identidad
Utilizando la aproximación del logaritmo anterior, aplicada tanto a como , la ecuación anterior da:
Por lo tanto, una aproximación de es:
que está escrito en el código como
yo = 0x5f3759df - ( yo >> 1 );
El primer término anterior es el número mágico
de lo que se puede inferir que . El segundo término, , se calcula desplazando los bits de una posición hacia la derecha. [27]
Con como raíz cuadrada inversa, . La aproximación obtenida con los pasos anteriores se puede refinar utilizando un método de búsqueda de raíces , un método que encuentra el cero de una función . El algoritmo utiliza el método de Newton : si hay una aproximación, para , entonces se puede calcular una mejor aproximación tomando , donde es la derivada de en . [28] Para lo anterior ,
donde y .
Tratarlo como un número de punto flotante es equivalente ay = y*(threehalfs - x/2*y*y);
Al repetir este paso, utilizando la salida de la función ( ) como entrada de la siguiente iteración, el algoritmo hace que converja a la raíz cuadrada inversa. [29] Para los fines del motor Quake III , solo se utilizó una iteración. Una segunda iteración permaneció en el código, pero fue comentada . [21]
Como se señaló anteriormente, la aproximación es muy precisa. El gráfico único de la derecha representa el error de la función (es decir, el error de la aproximación después de que se ha mejorado ejecutando una iteración del método de Newton), para entradas que comienzan en 0,01, donde la biblioteca estándar da como resultado 10,0 e InvSqrt() da 9,982522, lo que hace que la diferencia relativa sea 0,0017478, o el 0,175 % del valor verdadero, 10. El error absoluto solo disminuye a partir de ese momento, y el error relativo se mantiene dentro de los mismos límites en todos los órdenes de magnitud.
No se sabe con precisión cómo se determinó el valor exacto del número mágico. Chris Lomont desarrolló una función para minimizar el error de aproximación eligiendo el número mágico en un rango. Primero calculó la constante óptima para el paso de aproximación lineal como 0x5F37642F , cerca de 0x5F3759DF , pero esta nueva constante dio un poco menos de precisión después de una iteración del método de Newton. [30] Luego, Lomont buscó una constante óptima incluso después de una y dos iteraciones de Newton y encontró 0x5F375A86 , que es más precisa que la original en cada etapa de iteración. [30] Concluyó preguntando si el valor exacto de la constante original se eligió a través de derivación o prueba y error . [31] Lomont dijo que el número mágico para el tipo de tamaño IEEE754 de 64 bits doble es 0x5FE6EC85E7DE30DA , pero más tarde Matthew Robertson demostró que era exactamente 0x5FE6EB50C7B537A9 . [32]
Jan Kadlec redujo el error relativo en un factor adicional de 2,7 ajustando también las constantes en la iteración única del método de Newton, [33] llegando después de una búsqueda exhaustiva a
conv.i = 0x5F1FFFF9 - ( conv.i >> 1 ) ; conv.f * = 0,703952253f * ( 2,38924456f - x * conv.f * conv.f ) ; devolver conv.f ;
Ahora está disponible un análisis matemático completo para determinar el número mágico para números de punto flotante de precisión simple. [34] [35]
El método intermedio entre una y dos iteraciones del método de Newton en términos de velocidad y precisión es una única iteración del método de Halley . En este caso, el método de Halley es equivalente a aplicar el método de Newton con la fórmula inicial . El paso de actualización es entonces
donde la implementación debe calcularse solo una vez, a través de una variable temporal.
Las incorporaciones posteriores de los fabricantes de hardware han hecho que este algoritmo sea redundante en su mayor parte. Por ejemplo, en x86 , Intel introdujo la instrucción SSErsqrtss
en 1999. En una prueba comparativa de 2009 en el Intel Core 2 , esta instrucción tomó 0,85 ns por punto flotante en comparación con los 3,54 ns del algoritmo rápido de raíz cuadrada inversa, y tuvo menos errores. [1]
Algunos sistemas embebidos de bajo costo no tienen instrucciones especializadas de raíz cuadrada. Sin embargo, los fabricantes de estos sistemas suelen proporcionar bibliotecas trigonométricas y de otras matemáticas basadas en algoritmos como CORDIC .
long
reduce la portabilidad de este código en los sistemas modernos. Para que el código se ejecute correctamente, sizeof(long)
debe tener 4 bytes, de lo contrario pueden producirse resultados negativos. En muchos sistemas modernos de 64 bits, sizeof(long)
tiene 8 bytes. El sustituto más portátil es int32_t
.{{cite web}}
: CS1 maint: bot: original URL status unknown (link)