Algoritme de la raïl quadrada inversa ràpida

Algoritme de la raïl quadrada inversa ràpida

Els càlculs de l'allumenament i reflexió (mostrats ací en el shooter en primera persona OpenArena) usen el còdic de la raïl quadrada inversa ràpida per a calcular els ànguls d'incidència i reflexos.

La raïl quadrada inversa ràpida, a voltes coneguda com Fast InvSqrt() o per la constant hexadecimal 0x5F3759DF, és un algoritme que estima $\frac{1}{\sqrt{x}}$, el recíproc (o invers multiplicativo) de la raïl quadrada d'un número en punt flotant de 32 bits. Esta operació s'utilisa en el processament digital de senyals per a normalisar un vector, és dir, convertir-ho en un vector de mòdul 1. Per eixemple, els programes de gràfics per ordenador utilisen les raïls quadrades inverses per a calcular els ànguls d'incidència i reflexió per a l'allumenament i el sombrejat. L'algoritme és més conegut per la seua implementació en 1999 en el còdic font del videojoc de dispars en primera persona Quake III Arena,  que feya un gran us dels gràfics en 3D. L'algoritme no va escomençar a aparéixer en fòrums públics com Usenet fins a 2002 o 2003. El càlcul de les raïls quadrades sol dependre de moltes operacions de divisió, que per als números en menge flotant són computacionalment costoses. El quadrat invers ràpit genera una bona aproximació en un sol pas de divisió. S'han descobert atres videojocs anteriors a Quake 3 que utilisen un algoritme similar, encara que l'implementació de Quake seguix sent l'eixemple més conegut.

L'algoritme accepta un número en menge flotant de 32 bits com a entrada i almagasena un valor dividit a la mitat per al seu us posterior. A continuació, tractant els bits que representen el número en punt flotant com un sancer de 32 bits, es realisa un desplaçament llògic a la dreta d'un bit i el resultat es resta del número 0x5F3759DF (en notació decimal: 1.597.463.007), que és una representació en punt flotant d'una aproximació de $\sqrt{2^{127}}$ . Açò resulta en la primera aproximació de la raïl quadrada inversa de l'entrada. Tractant els bits de nou com un número en punt flotant, s'eixecuta una iteración del método de Newton, produint una aproximació més precisa.

L'algoritme es va atribuir originalment a John Carmack, pero una investigació va demostrar que el còdic tenia raïls més profundes tant en el hardware com en el software dels gràfics per ordenador. Els ajusts i alteracions varen passar per Silicon Graphics i 3dfx Interactive, i la constant original es va derivar d'una colaboració entre Cleve Moldre i Gregory Walsh, mentres Gregory treballava per a Ardent Computing a finals de la década dels 80. Walsh i Moldre varen adaptar la seua versió a partir d'un document inèdit de William Kahan i K.C. Ng difòs en maig de 1986.

En els posteriors alvanços en el hardware, especialment els rsqrtss en instruccions SSE per a x86, este método no és aplicable en general a l'informàtica de propòsit general, encara que seguix sent un eixemple interessant històricament, aixina com per a màquines més llimitades, com els sistemes de baix cost. No obstant, cada volta són més els fabricants que inclouen acceleradors trigonométricos i atres acceleradors matemàtics com CORDIC, obviant la necessitat d'estos algoritmes.

Motivació

Les superfícies normals són molt utilisades en els càlculs de l'allumenament i el sombrejat, que requerixen els càlculs de les normes dels vectores. Ací es mostra un camp de vectores normals a una superfície.

Un eixemple bidimensional de cóm utilisar la normal $C$ per a trobar l'àngul de reflexió i l'àngul d'incidència; en este cas, llum reflectida en un espill curve. L'algoritme de la raïl quadrada inversa ràpida generalisa este càlcul a espai tridimensional.

La raïl quadrada d'un número en menge flotant s'utilisa per a calcular un vector unitari. Els programes poden utilisar vectores unitaris per a determinar els ànguls d'incidència i reflexió. Els programes de gràfics 3D deuen portar a terme millons d'estos càlculs cada segon per a simular l'allumenament. Quan el còdic va ser desenrollat a principis del 1990, la majoria de la potència de processament del punt flotant es va quedar arrere sobre la velocitat de processament dels número entero. Açò era difícil per als programa de gràfics 3D abans de l'arribada d'hardware específic que s'encarregara de la transformació i allumenament.

La magnitut del vector es determina calculant la seua norma euclidiana: la raïl quadrada de la suma de quadrats dels components del vector. Quan cada component del vector es dividix per esta distància el nou vector serà un vector unitari que senyalarà cap a la mateixa direcció i sentit. En un programa de gràfics 3D tots els vectores estan en l'espai tridimensional aixina que v seria un vector $𝒗$$(v_1,v_2,v_3)$

$\Vert 𝒗\Vert =\sqrt{v_1^2+v_2^2+v_3^2}$

és la norma euclidiana o mòdul del vector.

$\widehat{𝒗}=\frac{𝒗}{\Vert 𝒗\Vert }$

és el vector unitari. Utilisant $\Vert 𝒗{\Vert }^2$ per a representar $v_1^2+v_2^2+v_3^2$

$\widehat{𝒗}=\frac{𝒗}{\sqrt{{\Vert 𝒗\Vert }^2}}$

el qual relaciona el vector d'unitat en l'invers de la raïl quadrada dels components de la distància distancia. L'invers de la raïl quadrada es pot usar per a calcular $\widehat{𝒗}$

$\widehat{𝒗}=𝒗\frac{1}{\sqrt{{\Vert 𝒗\Vert }^2}}$

on el terme fracció és l'invers de la raïl quadrada de $\Vert 𝒗{\Vert }^2$.

En aquella l'época, la divisió en punt flotant era més llarga comparada en la multiplicació; l'algoritme de la raïl quadrada inversa ràpida es va botar el pas de la divisió, proporcionant  ventaja en el rendiment. Quake III Arena, un videojoc de dispars en primera persona, va utilisar este algoritme per a accelerar la computació dels gràfics. L'algoritme des de llavors ha segut implementat en hardware de sombrejatres (shaders) de vèrtiços utilisant una matriu de portes llògiques programable (FPGA).

Visió general del còdic

El següent còdic és l'implementació de la raïl quadrada inversa ràpida de Quake III Arena, tret de les directives del preprocesador de C, pero incloent el comentari de text original:

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

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

	return i;
}

En aquell moment, el método general per a calcular la raïl quadrada inversa era calcular una aproximació per a

$\frac{1}{\sqrt{x}}$

i després comprovar l'aproximació per mig d'un atre método fins que s'obtenia en un marge d'error acceptable respecte al resultat correcte. El càlcul de raïls quadrades a principi de la década de 1990 obtenia estes aproximacions a partir d'una taula de consulta. La clau de la raïl quadrada inversa ràpida era calcular directament una aproximació utilisant l'estructura dels números en punt flotant, demostrant ser més ràpida que consultar estes taules. L'algoritme era aproximadament quatre voltes més ràpit que calcular la raïl quadrada en un atre método i obtindre l'inversa per mig de la divisió en punt flotant. L'algoritme va ser dissenyat pensant en l'especificació en punt flotant de 32 bits d'IEEE 754-1985, pero l'investigació de Chris Lomont va demostrar que podia ser implementat en atres especificacions en punt flotant.

Les ventages en velocitat oferides pel truc de la raïl quadrada inversa ràpida procedixen de tractar la paraula en punt flotant de 32 bits com un número entero i després restar-li-ho a una constant "màgica", 0x5F3759DF. Este canvi de bits i resta de sancers proporciona un patró que, quan és redefinit com un número en punt flotant, és una aproximació de la raïl quadrada inversa del número. Es realisa una iteración del método de Newton per a guanyar una miqueta més d'exactitut, i el còdic està terminat. L'algoritme genera resultats raonablement precisos utilisant una única primera aproximació del método de Newton; no obstant, és molt més llent i menys precís que usar l'instrucció de SSE rsqrtss en els processadors x86 també llançats en 1999.

Eixemple pràctic

Per eixemple, el número .$x=0.15625$ es pot usar per a calcular$\frac{1}{\sqrt{x}}\approx 2.52982$ Els primers passos de l'algoritme estan ilustrats avall:

0011_1110_0010_0000_0000_0000_0000_0000  Patró de bits de "x" i "i"
0001_1111_0001_0000_0000_0000_0000_0000  moure bits cap a la dreta una posició: (i >> 1)
0101_1111_0011_0111_0101_1001_1101_1111  L'número mágico 0x5F3759DF
0100_0000_0010_0111_0101_1001_1101_1111  El resultat de 0x5F3759DF - (i >> 1)

Interpretació com una representació 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 últim patró de bits com a números en punt flotant dona l'aproximació $y=2.61486$, que té un error d'al voltant del 3,4%. Despuix d'una sola iteración del método de Newton, el resultat final és $y=2.52549$ , en un error tan sol de 0,17%.

Evitar comportaments inesperats

Segons la biblioteca estàndar de C, reinterpretar un valor en punt flotant com un número entero eliminant-li la busca pot aplegar a causar comportament inesperat (comportament indefinit). Una atra manera de fer-ho seria situar el valor en punt flotant en una unió de senyes anònima que continga un membre sancer sense signe de 32 bits adicional, i els accessos a est proporcionarien una visió a nivell de bits del contingut del valor de punt flotant.

#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 les versions modernes de C++ (C++20 en avant), el método més comú per a implementar el cast d'esta funció és a través de std::bit_cast. Açò també permet a la funció treballar en un context 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 i = std::bit_cast<float>(
        0x5f3759df - (std::bit_cast<std::uint32_t>(number) >> 1));
	return i * (1.5f - (number * 0.5f * i * i));
}