Raíz cuadrada inversa rápida

La versión actual de la página aún no ha sido revisada por colaboradores experimentados y puede diferir significativamente de la versión revisada el 17 de enero de 2022; las comprobaciones requieren 76 ediciones .

Fast Reverse Square Root (también Fast InvSqrt() o 0x5F3759DF por la constante "mágica" utilizada , en decimal 1 597 463 007) es un algoritmo aproximado rápido para calcular la raíz cuadrada inversa para números de coma flotante positivos de 32 bits . El algoritmo utiliza operaciones de enteros "restar" y " desplazamiento de bits ", así como "restar" fraccionario y "multiplicar" - sin operaciones lentas "dividir" y "raíz cuadrada". A pesar de ser "hacky" a nivel de bits, la aproximación es monótona y continua : los argumentos cercanos dan resultados cercanos. La precisión (menos del 0,2% hacia abajo y nunca hacia arriba) [1] [2] no es suficiente para cálculos numéricos reales, pero es suficiente para gráficos tridimensionales . $y={\frac{1}{{\sqrt {x))))$

Algoritmo

El algoritmo toma un número de punto flotante de 32 bits ( precisión simple en formato IEEE 754 ) como entrada y realiza las siguientes operaciones en él:

Tratando un número fraccionario de 32 bits como un entero, realice la operación y 0 = 5F3759DF 16 − (x >> 1) , donde >> es un desplazamiento de bits hacia la derecha. El resultado se trata nuevamente como un número fraccionario de 32 bits.
Para aclarar, se puede realizar una iteración del método de Newton : y 1 \u003d y 0 (1.5 - 0.5xy 0 ²) .

Algoritmo original:

float Q_rsqrt ( número flotante ) { yo largo ; flotar x2 , y ; const float tres mitades = 1.5F ; x2 = número * 0.5F ; y = número ; i = * ( largo * ) & y ; // malvada piratería de nivel de bit de punto flotante yo = 0x5f3759df - ( yo >> 1 ); // ¿Qué diablos? y = * ( flotante * ) & yo ; y = y * ( tres mitades - ( x2 * y * y ) ); // 1ra iteración // y = y * ( tres mitades - ( x2 * y * y ) ); // 2da iteración, esto se puede eliminar devolver y ; }

Correcto según los estándares de la implementación C moderna, teniendo en cuenta posibles optimizaciones y multiplataforma :

float Q_rsqrt ( número flotante ) { const flotante x2 = número * 0.5F ; const float tres mitades = 1.5F ; unión { flotante f ; uint32_ti ; _ } conversión = { número }; // miembro 'f' establecido en el valor de 'número'. conv . i = 0x5f3759df - ( conv . i >> 1 ); conv . f *= tres mitades - x2 * conv . f * conv . f ; conversión de retorno f ; }

La implementación de Quake III: Arena [3] considera que la longitud es igual a , y utiliza punteros para la conversión (la optimización "si , ninguno ha cambiado" puede funcionar de forma errónea; en GCC, al compilar para "liberar", aparece una advertencia motivado). Y también contiene una palabra obscena : John Carmack , al exponer el juego en el dominio público, no entendió lo que se estaba haciendo allí. floatlongfloatlong

En C++20 , puede usar la nueva función . bit_cast

#incluir <bit> #incluir <límites> #incluir <cstdint> constexpr float Q_rsqrt ( número flotante ) noexcept { static_assert ( std :: numeric_limits < float >:: is_iec559 ); // Comprobar si la máquina de destino es compatible float const y = std :: bit_cast < float > ( 0x5f3759df - ( estándar :: bit_cast < estándar :: uint32_t > ( número ) >> 1 )); devuelve y * ( 1.5f - ( número * 0.5f * y * y )); }

Historia

El algoritmo probablemente fue desarrollado en Silicon Graphics en la década de 1990, y apareció una implementación en 1999 en el código fuente del juego de computadora Quake III Arena , pero el método no apareció en foros públicos como Usenet hasta 2002-2003. El algoritmo genera resultados razonablemente precisos utilizando la única primera aproximación del método de Newton . En ese momento, la principal ventaja del algoritmo era el rechazo de los costosos cálculos de punto flotante en favor de las operaciones con enteros. Las raíces cuadradas inversas se utilizan para calcular los ángulos de incidencia y reflexión para la iluminación y el sombreado en los gráficos por computadora.

El algoritmo se atribuyó originalmente a John Carmack , pero sugirió que Michael Abrash , un especialista en gráficos, o Terje Mathisen, un especialista en lenguaje ensamblador, lo trajeron a id Software . El estudio del problema mostró que el código tenía raíces más profundas tanto en las áreas de hardware como de software de los gráficos por computadora. Tanto Silicon Graphics como 3dfx Interactive han realizado correcciones y cambios , y la primera versión conocida fue escrita por Gary Tarolli para SGI Indigo . Quizás el algoritmo fue inventado por Greg Walsh y Clive Moler , los colegas de Gary en Ardent Computer [5] .

Jim Blinn, un especialista en gráficos 3D, ha propuesto un método de raíz cuadrada inversa tabular similar [6] que cuenta hasta 4 dígitos (0,01%) y se encontró al desmontar el juego Interstate '76 (1997) [7] .

Con el lanzamiento de 3DNow! Los procesadores AMD introdujeron la instrucción del ensamblador PFRSXRT [8] para un cálculo aproximado rápido de la raíz cuadrada inversa. La versión para el doble no tiene sentido: la precisión de los cálculos no aumentará [2] , por lo tanto, no se agregó. En 2000, se añadió a SSE2 la función RSQRTSS [9] , que es más precisa que este algoritmo (0,04 % frente a 0,2 %).

Análisis y error

La representación de bits de un número fraccionario de 4 bytes en formato IEEE 754 se ve así:

Señal
	Ordenar								mantisa
0	0	una	una	una	una	una	0	0	0	una	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	$=(1+2^{-2})\cdot 2^{-3}=0{,}15625$
31				24				23				dieciséis				quince				ocho				7				0

Tratamos solo con números positivos (el bit de signo es cero), no desnormalizados , no ∞ y no NaN . Dichos números se escriben en forma estándar como 1,mmmm 2 ·2 e . La parte 1, mmmm se llama mantisa , e es el orden . La unidad principal no se almacena ( unidad implícita ), así que llamemos al valor 0, mmmm la parte explícita de la mantisa. Además, los números fraccionarios de máquina tienen un orden desplazado : 2 0 se escribe como 011.1111.1 2 .

En números positivos , la biyección "fraccional ↔ entero" (indicada a continuación como ) es continua como una función lineal por partes y monótona . De esto podemos afirmar inmediatamente que la raíz inversa rápida, como combinación de funciones continuas, es continua. Y su primera parte, desplazamiento-resta, también es monótona y lineal por partes. La biyección es complicada, pero casi "gratuita": según la arquitectura del procesador y las convenciones de llamada , no necesita hacer nada o mover el número del registro fraccionario al entero. ${\ Displaystyle I_ {x}}$

Por ejemplo, la representación binaria del entero hexadecimal 0x5F3759DF es 0|101.1111.0|011.0111.0101.1001.1101.1111 2 (los puntos son límites de nibble, las líneas verticales son límites de campo fraccional de computadora). El orden de 101 1111 0 2 es 190 10 , luego de restar el desplazamiento 127 10 obtenemos el exponente 63 10 . La parte explícita de la mantisa 01 101 110 101 100 111 011 111 2 después de agregar la unidad principal implícita se convierte en 1.011 011 101 011 001 110 111 11 2 = 1.432 430 148… 10 . Teniendo en cuenta la precisión real de la computadora fraccional 0x5F3759DF ↔ 1.4324301 10 2 63 .

Denotamos la parte explícita de la mantisa del número , - orden no desplazado, - longitud en bits de la mantisa, - desplazamiento del orden. El número , escrito en una cuadrícula de bits logarítmica lineal de fracciones de computadora, se puede [10] [3] aproximar mediante una cuadrícula logarítmica como , donde es un parámetro utilizado para ajustar la precisión de la aproximación. Este parámetro va de 0 (exacto en y ) a 0,086 (exacto en un punto, ) ${\ estilo de visualización m_ {x} \ en [0,1)}$ $X$ $e_{x}\en \mathbb {Z}$ ${\ estilo de visualización L = 2 ^ {23}}$ ${\ estilo de visualización B = 127}$ $x\equiv 2^{e_{x}}(1+m_{x})$ $\log _{2}x\equiv e_{x}+\log _{2}(1+m_{x})\approx e_{x}+m_{x}+\sigma$ $\sigma$ ${\ estilo de visualización m_ {x} = 0}$ $una$ $m_{x}=0{,}443$

Usando esta aproximación, la representación entera de un número se puede aproximar como $X$

I_{x}\equiv L(e_{x}+B+m_{x})\approx L\log _{2}x+L(B-\sigma )

En consecuencia, . $\log _{2}x\approx {\frac {I_{x}}{L}}-(B-\sigma )$

Hagamos lo mismo [3] para (respectivamente ) y obtengamos ${\ Displaystyle y = {\ tfrac {1} {\ sqrt {x}}}}$ $\log _{2}y=-{\tfrac {1}{2}}\log _{2}x$

I_{y}\approx {\tfrac {3}{2}}L(B-\sigma )-{\tfrac {1}{2}}I_{x}

y\approx I^{-1}\left[{\tfrac {3}{2}}L(B-\sigma )-{\tfrac {1}{2}}I_{x}\right]

La constante mágica , teniendo en cuenta los límites , en aritmética fraccionaria tendrá un orden imparcial y una mantisa ), y en notación binaria - 0|101.1111.0|01 1 ... (1 es una unidad implícita; 0.5 proviene del orden ; una pequeña unidad corresponde al rango [1.375; 1.5) y por lo tanto es altamente probable, pero no está garantizado por nuestras estimaciones). ${\tfrac {3}{2}}L(B-\sigma )$ $\sigma$ $\left\lpiso {\tfrac {3L}{2L}}(B-\sigma)\right\rpiso -B=63$ $c=1+0{,}5-{\tfrac {3}{2}}\sigma \in (1{,}37;1{,}5$

Es posible calcular a qué equivale la primera aproximación lineal por partes [11] (en la fuente, no se usa la mantisa en sí, sino su parte explícita ): ${\ estilo de visualización t = c-1}$

para : ; $x\en [0{,}5;\;c-0{,}5)$ $y_{01}=-x+t+{\tfrac {3}{2}}=-x+c+{\tfrac {1}{2}}$
para : ; $x\en [c-0{,}5;\;1)$ $y_{02}=-{\tfrac {1}{2}}x+{\tfrac {1}{2}}t+{\tfrac {5}{4}}=-{\tfrac {1}{ 2}}x+{\tfrac{1}{2}}c+{\tfrac {3}{4}}$
para : . ${\ estilo de visualización x \ en [1; \; 2)}$ $y_{03}=-{\tfrac {1}{4}}x+{\tfrac {1}{2}}t+1=-{\tfrac {1}{4}}x+{\tfrac { 1}{2}}c+{\tfrac {1}{2}}$

En mayor o menor , el resultado cambia proporcionalmente: cuando se cuadruplica , el resultado se reduce exactamente a la mitad. $X$ $X$

El método de Newton da [11] , y . La función es decreciente y convexa hacia abajo; en tales funciones, el método de Newton se acerca al valor verdadero desde la izquierda; por lo tanto, el algoritmo siempre subestima la respuesta. $f(y)={\frac {1}{y^{2))}-x$ ${\displaystyle f'(y)=-{\frac {2}{y^{3))))$ $y_{n+1}=y_{n}-{\frac {f(y_{n})}{f'(y_{n))}}={\frac {y_{n}(3- xy_{n}^{2})}{2}}=y_{n}(1{,}5-0{,}5\,xy_{n}^{2})$ $f(y)$

No se sabe de dónde vino la constante 0x5F3759DF ↔ [a] 1.4324301 2 63 . Por fuerza bruta, Chris Lomont y Matthew Robertson descubrieron [1] [2] que la mejor constante [b] en términos de error relativo marginal para float es 0x5F375A86 ↔ 1.4324500 2 63 , para double es 0x5FE6EB50C7B537A9. Cierto, para el doble el algoritmo no tiene sentido (no proporciona una ganancia en precisión en comparación con la flotación) [2] . La constante de Lomont también se obtuvo analíticamente ( c = 1.432450084790142642179 ) [b] , pero los cálculos son bastante complicados [11] [2] .

Después de un paso del método de Newton, el resultado es bastante preciso ( +0 % -0,18 % ) [1] [2] , lo cual es más que adecuado para propósitos de gráficos por computadora ( 1 ⁄ 256 ≈ 0,39 % ). Tal error se conserva en todo el rango de números fraccionarios normalizados. Dos pasos dan una precisión de 5 dígitos [1] , después de cuatro pasos se alcanza el error del doble . Si lo desea, puede reequilibrar el error multiplicando los coeficientes 1,5 y 0,5 por 1,0009 para que el método dé simétricamente ±0,09%; esto es lo que hicieron en el juego Interstate '76 y el método Blinn, que también iteran el método de Newton [7 ] .

El método de Newton no garantiza la monotonicidad, pero la enumeración por computadora muestra que la monotonicidad existe.

Código fuente (C++) #incluir <iostream> unión FlotanteInt { flotar comoFlotar ; int32_t comoInt ; }; int floatToInt ( flotante x ) { FlotanteInt r ; r _ comoFlotante = x ; volver r . asInt ; } flotante intToFloat ( int x ) { FlotanteInt r ; r _ comoInt = x ; volver r . comoFlotador ; } float Q_rsqrt ( número flotante ) { yo largo ; flotar x2 , y ; const float tres mitades = 1.5F ; x2 = número * 0.5F ; y = número ; i = * ( largo * ) & y ; // piratería malvada de nivel de bit de punto flotante i = 0x5f3759df - ( i >> 1 ); // ¿Qué diablos? y = * ( flotante * ) & yo ; // ¡No sé qué diablos! y = y * ( tres mitades - ( x2 * y * y ) ); // 1ra iteración devolver y ; } int principal () { int iStart = floatToInt ( 1.0 ); int iEnd = floatToInt ( 4.0 ); std :: cout << "Números por recorrer: " << iEnd - iStart << std :: endl ; int nProblemas = 0 ; float oldResult = std :: numeric_limits < float >:: infinity (); for ( int i = iStart ; i <= iEnd ; ++ i ) { float x = intToFloat ( i ); resultado flotante = Q_rsqrt ( x ); if ( resultado > resultado anterior ) { std :: cout << "Se encontró un problema en " << x << std :: endl ; ++ nProblemas ; } } std :: cout << "Problemas totales: " << nProblems << std :: endl ; devolver 0 ; }

Existen algoritmos similares para otras potencias, como la raíz cuadrada o cúbica [3] .

Motivación

La superposición "directa" de la iluminación en un modelo tridimensional , incluso uno de gran cantidad de polígonos, incluso teniendo en cuenta la ley de Lambert y otras fórmulas de reflexión y dispersión, dará inmediatamente un aspecto poligonal: el espectador verá la diferencia en la iluminación a lo largo del aristas del poliedro [12] . A veces esto es necesario, si el objeto es realmente angular. Y para los objetos curvilíneos, hacen esto: en un programa tridimensional, indican si un borde afilado o uno liso [12] . Dependiendo de esto, incluso al exportar el modelo al juego, las esquinas de los triángulos calculan la normal de la unidad de longitud a la superficie curva. Durante la animación y la rotación, el juego transforma estas normales junto con el resto de los datos 3D; al aplicar iluminación, se interpola sobre todo el triángulo y se normaliza (lo lleva a una unidad de longitud).

Para normalizar un vector, debemos dividir sus tres componentes por la longitud. O, mejor, multiplícalos por el recíproco de la longitud: . Se deben calcular millones de estas raíces por segundo. Antes de que se creara el hardware de procesamiento de iluminación y transformación dedicado , el software informático podía ser lento. En particular, a principios de la década de 1990, cuando se desarrolló el código, la mayoría de los cálculos de punto flotante se quedaron atrás en rendimiento con respecto a las operaciones con números enteros. $(x',y',z')=(x,y,z){\frac {1}{\sqrt {x^{2}+y^{2}+z^{2))} }$

Quake III Arena utiliza un algoritmo rápido de raíz cuadrada inversa para acelerar el procesamiento de gráficos por parte de la CPU, pero desde entonces el algoritmo se ha implementado en algunos sombreadores de vértices de hardware especializados que utilizan matrices programables dedicadas (FPGA).

Incluso en computadoras de la década de 2010, dependiendo de la carga del coprocesador fraccional , la velocidad puede ser de tres a cuatro veces mayor que usando funciones estándar [11] .

Comentarios

↑ Aquí, la flecha ↔ significa la biyección de la representación binaria de un número entero y la representación binaria de un número de coma flotante en el formato IEEE 754 explicado anteriormente .
↑ 1 2 Si pones 127 en el campo de pedido, obtienes 0x3FB75A86. La biblioteca GRISU2, que es completamente entera e independiente de las sutilezas del coprocesador, dice que 0x3FB75A86 ↔ 1.43245 es el número decimal más corto que se convierte en un flotante dado . Sin embargo, la unidad menos significativa sigue siendo 1,19 10 −7 y 0x3FB75A86 = 1,432450056 ≈ 1,4324501. El siguiente fraccionario es 0x3FB75A87 ↔ 1.4324502 sin sutilezas. De ahí el redondeo no intuitivo de 1,43245008 a 1,4324500.

Notas

↑ 1 2 3 4 Copia archivada . Consultado el 25 de agosto de 2019. Archivado desde el original el 6 de febrero de 2009. (indefinido)
↑ 1 2 3 4 5 6 https://web.archive.org/web/20140202234227/http://shelfflag.com/rsqrt.pdf
↑ 1 2 3 4 Hummus e imanes . Consultado el 1 de febrero de 2017. Archivado desde el original el 13 de enero de 2017. (indefinido)
↑ Beyond3D - Origen de Fast InvSqrt() de Quake3 . Consultado el 4 de octubre de 2019. Archivado desde el original el 10 de abril de 2017. (indefinido)
↑ Beyond3D - Origen de Fast InvSqrt() de Quake3 - Segunda parte . Consultado el 25 de agosto de 2019. Archivado desde el original el 25 de agosto de 2019. (indefinido)
↑ Trucos de punto flotante | Diarios y revistas IEEE | Explorar IEEE
↑ 1 2 Raíz cuadrada recíproca rápida... ¿en 1997? — Blog de Shane Peelar
↑ PFRSQRT - Calcular el valor aproximado del recíproco de la raíz cuadrada de un valor real corto - Club155.ru . Consultado el 4 de octubre de 2019. Archivado desde el original el 16 de octubre de 2019. (indefinido)
↑ RSQRTSS - Calcular el recíproco de la raíz cuadrada del valor escalar de punto flotante de precisión simple . Consultado el 6 de octubre de 2019. Archivado desde el original el 12 de agosto de 2019. (indefinido)
↑ https://web.archive.org/web/20150511044204/http://www.daxia.com/bibis/upload/406Fast_Inverse_Square_Root.pdf
↑ 1 2 3 4 Cálculo de Schwidke de la raíz cuadrada envuelta con las constantes mágicas: pidkhid analítico . Consultado el 12 de junio de 2022. Archivado desde el original el 17 de abril de 2022. (indefinido)
↑ 1 2 3 Esta es la norma: qué son los mapas normales y cómo funcionan / Sudo Null IT News

Enlaces

C. Lomont, Raíz cuadrada inversa rápida , Informe técnico, 2003.
Una breve historia de InvSqrt por Matthew Robertson
0x5f3759df , investigaciones adicionales sobre la precisión y la generalización del algoritmo por Christian Plesner Hansen