Algoritmo de la raíz cuadrada inversa rápida

De Wikipedia, la enciclopedia libre
Ir a la navegación Ir a la búsqueda
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 , 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 . 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[editar]

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

es la norma euclidiana o módulo del vector.

es el vector unitario. Utilizando para representar

el cuál 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

dónde el término fracción es el inverso de la raíz cuadrada de .

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[editar]

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 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[editar]

Por ejemplo, el número . se puede usar para calcular Los primeros pasos del algoritmo están ilustrados abajo:

0011_1110_0010_0000_0000_0000_0000_0000  patrón de ambos x y i
0001_1111_0001_0000_0000_0000_0000_0000  Cambio derecho 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 , 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 , con un error tan solo de 0,17%.

Evitar comportamientos inesperados[editar]

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[editar]

El algoritmo calcula mediante los pasos siguientes:

  1. Asigna el argumento a un entero como manera de calcular una aproximación del logaritmo en binario
  2. Usa esta aproximación para calcular una aproximación de
  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[editar]

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 única, se proporciona aquí una breve descripción de esta representación. Para codificar un número real distinto de cero x, el primer paso es para escribir como número binario normalizado:

dónde el exponente es un entero , y es la representación binaria del "significante" . Ya que el bit solo antes del punto en el significante es siempre 1, no necesita ser almacenado. De esta forma, se calculan tres enteros sin signo.

  • ,el bit de signo, es si es positivo y si es negativo o cero (1 bit)
  • es el "exponente desviado", dónde es el "sesgo de exponente" (8 bits)
  • , dónde (23 bits)

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

Cuando un ejemplo, consideramos otra vez el número Normalizando se obtiene

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

Estos campos se codifican como son mostrados en la figura abajo:

Float w significand 2.svg

Con inverso de la raíz cuadrada,.

El método de Newton[editar]

La aproximación proporcionada por los primeros pasos puede ser refinada utilizando El método del Newton: si hay una aproximación , para , entonces una aproximación mejor puede ser calculada para tomando , dónde es la derivada de en .[1]​ Para ese

dónde y .

Tratando como un número en punto flotante, y = y* (threehalfs - x/2*y*y); es equivalente a

Repitiendo este paso, utilizando el resultado de la función () como la entrada de la iteración siguiente, provoca la convergencia de la sucesión al inverso raíz cuadrada. Para los propósitos del motor de Quake III motor, sólo se utilizó una iteración. Una segunda iteración se mantiene en el código pero marcada como comentario.[2]


Historia[editar]

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. {{Listaref|refs=[3][4][5][6]

Referencias[editar]

  1. Hardy, 1908, p. 323.
  2. Eberly, 2001, p. 2.
  3. «sqrt implementation in fdlibm». 
  4. Moler, Cleve. «Symplectic Spacewar». MATLAB Central - Cleve's Corner. MATLAB. Consultado el 21 de julio de 2014. 
  5. Sommefeldt, Rys (29 de noviembre de 2006). «Origin of Quake3's Fast InvSqrt()». Beyond3D. Consultado el 12 de febrero de 2009. 
  6. «quake3-1.32b/code/game/q_math.c». Quake III Arena. id Software. Archivado desde el original el 2 de agosto de 2020. Consultado el 21 de enero de 2017. 

Bibliografía[editar]