ukázka [SNIPPET] Fast Inverse Square Root

[ffredyk 167]

Zdravím přátelé, zase jsem narazil na zajímavý kousek kódu (evidentně již celkem profláknutej kousek kódu), který byl proslaven hrou Quake 3 Arena. 

jedná se o extrémně rychlou funkci pro výpočet 

Která v klasické funkci není přímo lehká pro výpočet (škodlivá pro time budget procesoru). Nicméně kdysi dávno ještě před Johnem Carmackem a Quakem se objevil kdosi jménem Greg Walsh,  který s touto magií přišel nejspíše první a trochu zachránil svět 3D grafiky v dobách, kdy zobrazení beztexturové krychle na obrazovce trvalo dvě a půl vteřiny.

K čemu se tato funkce vůbec využívá a proč to byl takový hit 3D světa? Protože zrovna tato funkce je zásadní pro výpočet kolmic povrchů 3D objektů, zastoupené 3D vektory o délce 1, pro vyjádření osvícení a odrazů. Těchto výpočtů probíhá spousta najednou a jejich výpočet zahrnuje normalizaci velkého množství vektorů. A jak se normalizuje vektor? Zjistíš délku vektoru a vydělíš s ní každou jeho část (x,y,z) - tedy vynásobíš každou část formulí:

Výpočet samotného x2 + y2 + z2 je relativně rychlý a jednoduchý. Hledání odmocniny a její využití v dělení je náročné.

Originální funkce z Quake III Arena:

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

Jedná se o malou a velmi efektivní rychlou super funkci. Na první pohled nás nejspíš všechny zaujal řádek "what the fuck?" - jedná se právě o to kouzlo celé funkce. Greg totiž přišel na způsob jakým efektivně napodobit matematickou funkci log2 , která se využívá při výpočtu "náročným" způsobem a kompletně ji obešel (tím odlehčil náročnost na výpočet). Kompletní matematické rovnice vám nevysvětlím já, nýbrž je najdete v přiložených zdrojových dokumentech, ve kterých jsem čerpal (ale je to fakt advanced shit, ne-matikářům to zruší hlavu).

Co se ve funkci vlastně tedy děje?

• Do proměnné x2 si uložíme polovinu z čísla
• Do proměnné y zkopírujeme naše číslo
• Změna typu čísla z float na integer (odstranění desetinné čárky na bit-levelu bez ztráty mantisy)
• Posunutí bitů integeru o 1 doprava (bit exponentu se přesune do mantisy) a výsledek se odečte od magic konstanty
• Přetypování zpět na float
• Dovýpočet přesnosti pomocí Newtonovi metody (1 iterace)
• (případná druhá iterace pro ještě přesnější výsledek)

Jak přesně funguje float na bit-levelu?

• Float na bit-levelu má tři části - sign, exponent a mantissa.
• Sign určuje znaménko čísla (pozitivní nebo negativní číslo).
• Exponent může být pozitivní i negativní, na bit-levelu je ale uložen jako unsigned (to je určeno pomocí sign části floatu). Hodnoty všech nul jsou rezervované pro nulové nebo neexistující výsledky a hodnoty všech jedniček jsou rezervované jako reprezentace nekonečna nebo neznámého.
• Mantisa je signifikantní číslo, které po umocnění exponentem, představuje uloženou hodnotu (první bit se neukládá - předpokládá se, že je vždy 1)

Bit shift operátor v tomto případě posunul jeden bit exponenta do mantisy a tím efektivně celé číslo vydělil 2. Poté odečítáme celé číslo od konstanty, která jako zázrakem vytvoří z našeho čísla relativně přesný odhad výpočtu odmocniny. Nikdo už se neshodl jakým způsobem bylo původní číslo nalezeno, nicméně existuje několik metod k nalezení toho nejlepšího (podrobnosti v odkazech na zdroje).

I když je tato funkce v moderní době irelevantní (moderní procesory a mikročipy se naučili pracovat s floating-point čísly efektivněji), její rychlosti se dalo využívat až do začátku 21. století, kdy už výkon floating-point logiky dohnal výkon integer výpočtů. Ale do té doby to byl nejlepší algoritmus v 3D grafice.

 

Zdroje:

• https://en.wikipedia.org/wiki/Floating-point_arithmetic
• https://en.wikipedia.org/wiki/Fast_inverse_square_root
• http://h14s.p5r.org/2012/09/0x5f3759df.html (Podrobnější popis celé funkce)
• http://www.lomont.org/papers/2003/InvSqrt.pdf (Detailní výpočet magic konstanty)

[bug 131]

Citace

Mantisa je signifikantní číslo, které po umocnění exponentem, představuje uloženou hodnotu (první bit se neukládá - předpokládá se, že je vždy 1)

tu by som chcel len doplniť že prvý bit sa síce neukladá (vynecháva) pri float (32bit) a double (64bit) kvôli šetreniu miesta v pamäti ale napríklad pri long double (80bit) už sa nevynecháva