stringtranslate.com

Raíz cuadrada inversa rápida

Los cálculos de iluminación y reflexión, como en el videojuego OpenArena , utilizan el código de raíz cuadrada inversa rápida para calcular ángulos de incidencia y reflexión.

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.

Historia

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]

Motivación

Las normales de superficie se utilizan ampliamente en los cálculos de iluminación y sombreado, lo que requiere el cálculo de normas para vectores. Aquí se muestra un campo de vectores normales a una superficie.
Ejemplo bidimensional del uso de la normal para hallar el ángulo de reflexión a partir del ángulo de incidencia; en este caso, sobre la luz que se refleja en un espejo curvo. Se utiliza la raíz cuadrada inversa rápida para generalizar este cálculo al espacio tridimensional.

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.

Descripción general del código

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

Ejemplo resuelto

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

Evitar comportamientos indefinidos

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

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 )); }                          

Algoritmo

El algoritmo se calcula realizando los siguientes pasos:

  1. Alias: argumento de un número entero como forma de calcular una aproximación del logaritmo binario.
  2. Utilice esta aproximación para calcular una aproximación de
  3. Alias ​​de regreso a un flotante, como una forma de calcular una aproximación del exponencial de base 2
  4. Refinar la aproximación utilizando una sola iteración del método de Newton.

Representación en punto flotante

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:

Alias ​​a un entero como logaritmo aproximado

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.

El entero asociado a un número de punto flotante (en azul), comparado con un logaritmo escalado y desplazado (en gris).

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

Primera aproximación del resultado

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]

El método de Newton

Error relativo entre el cálculo directo y la raíz cuadrada inversa rápida al realizar 0, 1, 2, 3 y 4 iteraciones del método de búsqueda de raíces de Newton. Nótese que se adopta la precisión doble y que la diferencia representable más pequeña entre dos números de precisión doble se alcanza después de realizar 4 iteraciones.

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]

Exactitud

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.

Mejoras posteriores

Número mágico

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]

Hallazgo cero

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.

Obsolescencia

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 .

Véase también

Notas

  1. ^ El uso del tipo longreduce 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.
  2. ^ debe estar en el rango para que pueda representarse como un número normal .
  3. ^ Los únicos números reales que se pueden representar exactamente como punto flotante son aquellos para los que es un entero. Otros números solo se pueden representar de forma aproximada redondeándolos al número exactamente representable más cercano.
  4. ^ Dado que es positivo, .

Referencias

  1. ^ abc Ruskin, Elan (16 de octubre de 2009). "Raíz cuadrada de sincronización". Requiere ensamblaje . Archivado desde el original el 8 de febrero de 2021. Consultado el 7 de mayo de 2015 .
  2. ^ feilipu. "z88dk es una colección de herramientas de desarrollo de software destinadas a las computadoras 8080 y z80". GitHub .
  3. ^ Munafo, Robert. "Propiedades notables de números específicos". mrob.com . Archivado desde el original el 16 de noviembre de 2018.
  4. ^ "Implementación de sqrt en fdlibm - Consulte la discusión de W. Kahan y KC Ng en los comentarios en la mitad inferior de este código".
  5. ^ Moler, Cleve (19 de junio de 2012). "Symplectic Spacewar". MATLAB Central - El rincón de Cleve . MATLAB . Consultado el 21 de julio de 2014 .
  6. ^ abc Sommefeldt, Rys (19 de diciembre de 2006). "Origen de Fast InvSqrt() de Quake3 - Segunda parte". Beyond3D . Consultado el 19 de abril de 2008 .
  7. ^ abcd Sommefeldt, Rys (29 de noviembre de 2006). "Origen de Fast InvSqrt() de Quake3". Beyond3D . Consultado el 12 de febrero de 2009 .
  8. ^ Blinn 1997, págs. 80–84.
  9. ^ Peelar, Shane (1 de junio de 2021). "Raíz cuadrada recíproca rápida... ¿en 1997?".
  10. ^ "Discusión sobre CSDN". Archivado desde el original el 2 de julio de 2015.
  11. ^ Lomont 2003, pág. 1-2.
  12. ^ Zafar, Saad; Adapa, Raviteja (enero de 2014). "Diseño de arquitectura de hardware y mapeo del algoritmo 'Fast Inverse Square Root'". Conferencia internacional sobre avances en ingeniería eléctrica (ICAEE) de 2014. págs. 1–4. doi :10.1109/ICAEE.2014.6838433. ISBN 978-1-4799-3543-7.S2CID2005623  .​
  13. ^ Middendorf 2007, págs. 155-164.
  14. ^ Blinn 2003, pág. 130.
  15. ^ "quake3-1.32b/code/game/q_math.c". Quake III Arena . id Software . Archivado desde el original el 2017-07-29 . Consultado el 2017-01-21 .{{cite web}}: CS1 maint: bot: original URL status unknown (link)
  16. ^ Eberly 2001, pág. 504.
  17. ^ Lomont 2003, pág. 1.
  18. ^ Lomont 2003.
  19. ^ Lomont 2003, pág. 3.
  20. ^ McEniry 2007, pág. 2, 16.
  21. ^ desde Eberly 2001, pág. 2.
  22. ^ Fog, Agner. "Listas de latencias de instrucciones, rendimientos y desgloses de microoperaciones para CPU Intel, AMD y VIA" (PDF) . Consultado el 8 de septiembre de 2017 .
  23. ^ Goldberg 1991, pág. 7.
  24. ^ Goldberg 1991, págs. 15-20.
  25. ^ Goldberg 1991, pág. 16.
  26. ^ McEniry 2007, pág. 3.
  27. ^ Hennessey y Patterson 1998, pág. 305.
  28. ^ Hardy 1908, pág. 323.
  29. ^ McEniry 2007, pág. 6.
  30. ^ desde Lomont 2003, pág. 10.
  31. ^ Lomont 2003, págs. 10-11.
  32. ^ Matthew Robertson (24 de abril de 2012). "Una breve historia de InvSqrt" (PDF) . UNBSJ.
  33. ^ Kadlec, Jan (2010). "Řrřlog::Mejora de la raíz cuadrada inversa rápida" (blog personal). Archivado desde el original el 9 de julio de 2018. Consultado el 14 de diciembre de 2020 .
  34. ^ Moroz y otros. 2018.
  35. ^ Muller, Jean-Michel (diciembre de 2020). «Funciones elementales y computación aproximada». Actas del IEEE . 108 (12): 2146. doi :10.1109/JPROC.2020.2991885. ISSN  0018-9219. S2CID  219047769.

Bibliografía

Lectura adicional

Enlaces externos