Algoritmo de la raíz cuadrada inversa rápida

Los cálculos de la iluminación y reflexión (mostrados aquí en el shooter en primera persona OpenArena) usan el código de la raíz cuadrada inversa rápida para calcular los ángulos de incidencia y reflejos.

La raíz cuadrada inversa rápida, a veces conocida como Fast InvSqrt() o por la constante hexadecimal 0x5F3759DF, es un algoritmo que estima ${\frac {1}{\sqrt {x}}}$, el recíproco (o inverso multiplicativo) de la raíz cuadrada de un número en punto flotante de 32 bits. Esta operación se utiliza en el procesamiento digital de señales para normalizar un vector, es decir, convertirlo en un vector de módulo 1. Por ejemplo, los programas de gráficos por ordenador utilizan las raíces cuadradas inversas para calcular los ángulos de incidencia y reflexión para la iluminación y el sombreado. El algoritmo es más conocido por su implementación en 1999 en el código fuente del videojuego de disparos en primera persona Quake III Arena,  que hacía un gran uso de los gráficos en 3D. El algoritmo no empezó a aparecer en foros públicos como Usenet hasta 2002 o 2003. El cálculo de las raíces cuadradas suele depender de muchas operaciones de división, que para los números en coma flotante son computacionalmente costosas. El cuadrado inverso rápido genera una buena aproximación con un solo paso de división. Se han descubierto otros videojuegos anteriores a Quake 3 que utilizan un algoritmo similar, aunque la implementación de Quake sigue siendo el ejemplo más conocido.

El algoritmo acepta un número en coma flotante de 32 bits como entrada y almacena un valor dividido a la mitad para su uso posterior. A continuación, tratando los bits que representan el número en 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 0x5F3759DF (en notación decimal: 1.597.463.007), que es una representación en punto flotante de una aproximación de ${\sqrt {2^{127}}}$ . Esto resulta en la primera aproximación de la raíz cuadrada inversa de la entrada. Tratando los bits de nuevo como un número en punto flotante, se ejecuta una iteración del método de Newton, produciendo una aproximación más precisa.

El algoritmo se atribuyó originalmente a John Carmack, pero una investigación demostró que el código tenía raíces más profundas tanto en el hardware como en el software de los gráficos por ordenador. Los ajustes y alteraciones pasaron por Silicon Graphics y 3dfx Interactive, y la constante original se derivó de una colaboración entre Cleve Moler y Gregory Walsh, mientras Gregory trabajaba para Ardent Computing a finales de la década de los 80. Walsh y Moler adaptaron su versión a partir de un documento inédito de William Kahan y K.C. Ng difundido en mayo de 1986.

Con los posteriores avances en el hardware, especialmente los rsqrtss en instrucciones SSE para x86, este método no es aplicable en general a la informática de propósito general, aunque sigue siendo un ejemplo interesante históricamente, así como para máquinas más limitadas, como los sistemas de bajo coste. Sin embargo, cada vez son más los fabricantes que incluyen aceleradores trigonométricos y otros aceleradores matemáticos como CORDIC, obviando la necesidad de estos algoritmos.

Motivación

Las superficies normales son muy utilizadas en los cálculos de la iluminación y el sombreado, que requieren los cálculos de las normas de los vectores. Aquí se muestra un campo de vectores normales a una superficie.

Un ejemplo bidimensional de cómo utilizar la normal $C$ para encontrar el ángulo de reflexión y el ángulo de incidencia; en este caso, luz reflejada en un espejo curvo. El algoritmo de la raíz cuadrada inversa rápida generaliza este cálculo a espacio tridimensional.

La raíz cuadrada de un número en coma flotante se utiliza para calcular un vector unitario. Los programas pueden utilizar vectores unitarios para determinar los ángulos de incidencia y reflexión. Los programas de gráficos 3D deben llevar a cabo millones de estos cálculos cada segundo para simular la iluminación. Cuando el código fue desarrollado a principios del 1990, la mayoría de la potencia de procesamiento del punto flotante se quedó atrás con respecto a la velocidad de procesamiento de los números enteros. Esto era difícil para los programa de gráficos 3D antes de la llegada de hardware específico que se encargase de la transformación e iluminación.

La magnitud del vector se determina calculando su norma euclidiana: la raíz cuadrada de la suma de cuadrados de los componentes del vector. Cuando cada componente del vector se divide por esta distancia el nuevo vector será un vector unitario que señalará hacia la misma dirección y sentido. En un programa de gráficos 3D todos los vectores están en el espacio tridimensional así que v sería un vector ${\boldsymbol {v}}$$(v_{1},v_{2},v_{3})$

$\|{\boldsymbol {v}}\|={\sqrt {v_{1}^{2}+v_{2}^{2}+v_{3}^{2}}}$

es la norma euclidiana o módulo del vector.

${\boldsymbol {\hat {v}}}={\frac {\boldsymbol {v}}{\left\|{\boldsymbol {v}}\right\|}}$

es el vector unitario. Utilizando $\|{\boldsymbol {v}}\|^{2}$ para representar $v_{1}^{2}+v_{2}^{2}+v_{3}^{2}$

${\boldsymbol {\hat {v}}}={\frac {\boldsymbol {v}}{\sqrt {\left\|{\boldsymbol {v}}\right\|^{2}}}}$

el cual relaciona el vector de unidad con el inverso de la raíz cuadrada de los componentes de la distancia distancia. El inverso de la raíz cuadrada se puede usar para calcular ${\boldsymbol {\hat {v}}}$

${\boldsymbol {\hat {v}}}={\boldsymbol {v}}\,{\frac {1}{\sqrt {\left\|{\boldsymbol {v}}\right\|^{2}}}}$

dónde el término fracción es el inverso de la raíz cuadrada de $\|{\boldsymbol {v}}\|^{2}$.

En aquella la época, la división en punto flotante era más larga comparada con la multiplicación; el algoritmo de la raíz cuadrada inversa rápida se saltó el paso de la división, proporcionando  ventaja en el rendimiento. Quake III Arena, un videojuego de disparos en primera persona, utilizó este algoritmo para acelerar la computación de los gráficos. El algoritmo desde entonces ha sido implementado en hardware de sombreadores (shaders) de vértices utilizando una matriz de puertas lógicas programable (FPGA).

Visión general del código

El siguiente código es la implementación de la raíz cuadrada inversa rápida de Quake III Arena, sacado de las directivas del preprocesador de C, pero incluyendo el comentario de texto original:

float Q_rsqrt( float number )
{
	long i;
	float x2, y;
	const float threehalfs = 1.5F;

	x2 = number * 0.5F;
	y  = number;
	i  = * ( long * ) &y;                       // evil floating point bit level hacking
	i  = 0x5f3759df - ( i >> 1 );               // what the fuck? 
	y  = * ( float * ) &i;
	y  = y * ( threehalfs - ( x2 * y * y ) );   // 1st iteration
//	y  = y * ( threehalfs - ( x2 * y * y ) );   // 2nd iteration, this can be removed

	return y;
}

En aquel momento, el método general para calcular la raíz cuadrada inversa era calcular una aproximación para ${\frac {1}{\sqrt {x}}}$ y luego comprobar la aproximación mediante otro método hasta que se obtenía con un margen de error aceptable respecto al resultado correcto. El cálculo de raíces cuadradas a principio de la década de 1990 obtenía estas aproximaciones a partir de una tabla de consulta. La clave de la raíz cuadrada inversa rápida era calcular directamente una aproximación utilizando la estructura de los números en punto flotante, demostrando ser más rápida que consultar estas tablas. El algoritmo era aproximadamente cuatro veces más rápido que calcular la raíz cuadrada con otro método y obtener la inversa mediante la división en punto flotante. El algoritmo fue diseñado pensando en la especificación en punto flotante de 32 bits de IEEE 754-1985, pero la investigación de Chris Lomont demostró que podía ser implementado en otras especificaciones en punto flotante.

Las ventajas en velocidad ofrecidas por el truco de la raíz cuadrada inversa rápida proceden de tratar la palabra en punto flotante de 32 bits como un número entero y luego restárselo a una constante "mágica", 0x5F3759DF. Este cambio de bits y resta de enteros proporciona un patrón que, cuando es redefinido como un número en punto flotante, es una aproximación de la raíz cuadrada inversa del número. Se realiza una iteración del método de Newton para ganar algo más de exactitud, y el código está terminado. El algoritmo genera resultados razonablemente precisos utilizando una única primera aproximación del método de Newton; sin embargo, es mucho más lento y menos preciso que usar la instrucción de SSE rsqrtss en los procesadores x86 también lanzados en 1999.​​

Ejemplo práctico

Por ejemplo, el número .$x=0.15625$ se puede usar para calcular${\frac {1}{\sqrt {x}}}\approx 2.52982$ Los primeros pasos del algoritmo están ilustrados abajo:

0011_1110_0010_0000_0000_0000_0000_0000  Patrón de bits de "x" e "i"
0001_1111_0001_0000_0000_0000_0000_0000  mover bits hacia la derecha una posición: (i >> 1)
0101_1111_0011_0111_0101_1001_1101_1111  El número mágico 0x5F3759DF
0100_0000_0010_0111_0101_1001_1101_1111  El resultado de 0x5F3759DF - (i >> 1)

Interpretación como una representación en IEEE 32-bits:

0_01111100_01000000000000000000000  1.25 × 2−3
0_00111110_00100000000000000000000  1.125 × 2−65
0_10111110_01101110101100111011111  1.432430... × 263
0_10000000_01001110101100111011111  1.307430... × 21

Reinterpretar este último patrón de bits como números en punto flotante da la aproximación $y=2.61486$, que tiene un error de alrededor del 3,4%. Después de una sola iteración del método de Newton, el resultado final es $y=2.52549$ , con un error tan solo de 0,17%.

Evitar comportamientos inesperados

Según la biblioteca estándar de C, reinterpretar un valor en punto flotante como un número entero eliminándole el puntero puede llegar a causar comportamiento inesperado (comportamiento indefinido). Otra manera de hacerlo sería situar el valor en punto flotante en una unión de datos anónima que contenga un miembro entero sin signo de 32 bits adicional, y los accesos a este proporcionarían una visión a nivel de bits del contenido del valor de punto flotante.

#include <stdint.h> // uint32_t
	
float Q_rsqrt( float number )
{
	const float x2 = number * 0.5F;
	const float threehalfs = 1.5F;
	
	union {
		float    f;
		uint32_t i;
	} conv = { .f = number };
	conv.i  = 0x5f3759df - ( conv.i >> 1 );
	conv.f *= threehalfs - ( x2 * conv.f * conv.f );
	return conv.f;
}

En las versiones modernas de C++ (C++20 en adelante), el método más común para implementar el cast de esta función es a través de std::bit_cast. Esto también permite a la función trabajar en un contexto de constexpr.

//only works after C++20 standard
#include <bit>
#include <limits>
#include <cstdint>

constexpr float Q_rsqrt(float number) noexcept
{
	static_assert(std::numeric_limits<float>::is_iec559); // (enable only on IEEE 754)

	float const y = std::bit_cast<float>(
        0x5f3759df - (std::bit_cast<std::uint32_t>(number) >> 1));
	return y * (1.5f - (number * 0.5f * y * y));
}

Algoritmo

El algoritmo calcula ${\frac {1}{\sqrt {x}}}$ mediante los pasos siguientes:

1. Asigna el argumento $x$ a un entero como manera de calcular una aproximación del logaritmo en binario $\log _{2}(x)$
2. Usa esta aproximación para calcular una aproximación de $\log _{2}\left({\frac {1}{\sqrt {x}}}\right)=-{\frac {1}{2}}\log _{2}(x)$
3. Reasigna a un decimal como manera de calcular una aproximación en potencia de base 2
4. Refina esa aproximación mediante una única interacción del método del newton.

Representación en punto flotante de simple precisión

Dado que este algoritmo se basa en gran medida en la representación a nivel de bits de los números de punto flotante de precisión simple, se proporciona aquí una breve descripción de esta representación. Para codificar un número real distinto de cero, el primer paso es escribirlo como número binario normalizado:

${\begin{aligned}x&=\pm 1.b_{1}b_{2}b_{3}\ldots \times 2^{e_{x}}\\&=\pm 2^{e_{x}}(1+m_{x})\end{aligned}}$

dónde el exponente $e_{x}$ es un entero ,$m_{x}\in [0,1)$, y $1.b_{1}b_{2}b_{3}\ldots$ es la representación binaria del significando $(1+m_{x})$. Ya que el bit antes del punto en el significando es siempre 1, no necesita ser almacenado. De esta forma, se calculan tres enteros sin signo.

• $S_{x}$, el bit de signo, es $0$ si $x$ es positivo y $1$ si es negativo o cero (1 bit)
• $E_{x}=e_{x}+B$, es el "exponente sesgado", dónde $B=127$ es el "sesgo de exponente" (8 bits)
• $M_{x}=m_{x}\times L$, dónde $L=2^{23}$ (23 bits)

Estos campos son entonces empaquetados, de izquierda a derecha, en un contenedor de 32 bits.

Como un ejemplo, considere otra vez el número $x=0.15625=0.00101_{2}$ Normalizando $x$ se obtiene

$x=+2^{-3}(1+0.25)$

Y así, los tres campos de entero sin signo son:

• $S=0$
• $E=-3+127=124=0111\ 1100_{2}$
• $M=0.25\times 2^{23}=2\ 097\ 152=0010\ 0000\ 0000\ 0000\ 0000\ 0000_{2}$

Estos campos se codifican como se muestra en la figura siguiente:

El método de Newton

La aproximación proporcionada por los primeros pasos del código puede ser refinada utilizando el método del Newton, el cual permite encontrar aproximaciones de los ceros o raíces de una función real. La función $f(y)={\frac {1}{y^{2}}}-x$ vale cero cuando $y={\frac {1}{\sqrt {x}}}$ y su derivada es $f'(y)=-{\frac {2}{y^{3}}}$. Hallar el valor $y$, tal que $f(y)=0$ nos permite encontrar el inverso multiplicativo de la raíz de $x$. El método de Newton propone que si hay una aproximación $y_{n}$, para $y$, entonces una mejor aproximación $y_{n+1}$ puede ser calculada tomando $y_{n+1}=y_{n}-{\frac {f(y_{n})}{f'(y_{n})}}$ , dónde $f'(y_{n})$ es la derivada de $f(y)$ evaluada en $y_{n}$.​ Entonces para esta función $f(y)$

$y_{n+1}=y_{n}-{\frac {{\frac {1}{y_{n}^{2}}}-x}{-{\frac {2}{y_{n}^{3}}}}}=y_{n}+{\frac {y_{n}^{3}}{2}}\left({\frac {1}{y_{n}^{2}}}-x\right)=y_{n}+{\frac {y_{n}}{2}}-y_{n}^{3}{\frac {x}{2}}$

$y_{n+1}=y_{n}\left({\frac {3}{2}}-y_{n}^{2}{\frac {x}{2}}\right)$

Lo cual es equivalente a la línea de código y = y * ( threehalfs - ( x2 * y * y ) );.

Repitiendo este paso sucesivamente, utilizando el resultado $y_{n+1}$ como la entrada de la iteración siguiente, el valor converge al inverso multiplicativo de la raíz raíz cuadrada de $x$. Para los propósitos del motor de Quake III, sólo se utilizó una iteración. Una segunda iteración se mantiene en el código pero fue marcada como comentario.​

Historia

El código fuente de Quake III Arena no se publicó hasta la QuakeCon 2005, pero ya en 2002 o 2003 aparecieron copias del código de raíz cuadrada inversa rápida en Usenet y otros foros. La especulación inicial apuntaba a John Carmack como el probable autor del código, pero los autores originales demostraron estar mucho más atrás en la historia de los gráficos 3D por computadora. Rys Sommefeldt concluyó que el algoritmo original fue diseñado por Greg Walsh en Ardent Computer en consulta con Cleve Moler, el creador de MATLAB. Cleve Moler se enteró de este truco gracias al código escrito por William Kahan y K.C. Ng en Berkeley alrededor de 1986.

Jim Blinn también demostró una aproximación simple de la raíz cuadrada inversa en una columna de 1997 para IEEE Computer Graphics and Applications. La ingeniería inversa de otros videojuegos 3D contemporáneos descubrió una variación del algoritmo en Interstate '76 de Activision de 1997, dos años antes de que se publicara Quake III Arena.

Referencias

Bibliografía

• Hardy, Godfrey (1908). A Course of Pure Mathematics. Cambridge, UK: Cambridge University Press. ISBN 0-521-72055-9. 
• Eberly, David (2001). 3D Game Engine Design. Morgan Kaufmann. ISBN 978-1-55860-593-0. 

• Datos: Q32980
• Multimedia: Fast inverse square root / Q32980