Radice quadrata inversa veloce

La radice quadrata inversa veloce (conosciuta in inglese come fast inverse square root o come Fast InvSqrt() o anche dalla costante esadecimale 0x5f3759df) è un metodo per calcolare x^−½, ovvero il reciproco (o inverso moltiplicativo) di una radice quadrata per un numero float a 32-bit nel formato "IEEE 754 single precision binary floating-point format: binary32". L'algoritmo è stato probabilmente sviluppato presso la Silicon Graphics all'inizio degli anni '90, l'implementazione apparve nel 1999 all'interno del codice sorgente del gioco Quake III Arena, ma il metodo non è apparso su forum pubblici come Usenet fino al 2002 o al 2003. All'epoca, il vantaggio principale dell'algoritmo proveniva dall'evitare operazioni su float, computazionalmente costose, a favore di operazioni su interi. La radice quadrata inversa è usata per calcolare l'angolo d'incidenza e i riflessi della luce e delle ombre nella computer grafica.

L'algoritmo accetta numeri float a 32-bit come input e ne salva il valore dimezzato per usarlo in seguito. Successivamente, trattando i bit che rappresentano il float come un integer a 32-bit, viene effettuato un logical shift a destra di un bit e il risultato è poi sottratto dalla "costante magica" 0x5f3759df. Questa è la prima approssimazione della radice quadrata inversa dell'input. Trattando i bit nuovamente come float si esegue una iterazione del metodo di Newton per ottenere una approssimazione più precisa. Questo calcola una approssimazione della radice quadrata inversa su un float con una velocità approssimativamente 4 volte più veloce della divisione su float.

Inizialmente l'algoritmo è stato attribuito a John Carmack, ma successive ricerche hanno mostrato che il codice ha radici più profonde che coinvolgono sia il software che l'hardware usato nella computer grafica. Modifiche ed alterazioni sono state introdotte nel tempo sia dalla Silicon Graphics che dalla 3dfx Interactive, fino ad arrivare al primo utilizzo noto con la implementazione di Gary Tarolli per il SGI Indigo.

Non è noto come la costante sia stata derivata originalmente anche se le ricerche hanno portato luce su alcuni possibili metodi di derivazione della stessa.

Il codice sorgente seguente è l'implementazione della radice quadrata inversa in Quake III Arena, senza le direttive al preprocessore C, ma con ancora i commenti originali^[1].

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

Per determinare la radice inversa, un programma dovrebbe calcolare un'approssimazione di ${\displaystyle x^{-1/2}}$ e poi dovrebbe applicare un metodo numerico per migliorare il risultato fino ad avere un errore accettabile. All'inizio degli anni 1990, una prima approssimazione del risultato veniva da una tabella di associazione^[2]. Questa nuova funzione si rivelò migliore delle tabelle di associazione e circa quattro volte più veloce di una divisione a virgola mobile classica^[3]. L'algoritmo è stato concepito secondo lo standard IEEE 754-1985 per numeri con virgola mobile 32-bit, ma Chris Lomont e Charles McEniry hanno mostrato come potesse essere implementato usando altre specifiche per la virgola mobile.

Il salto di velocità dell'algoritmo è dovuto al trattamento del double con il numero a virgola mobile, considerato con un intero che viene poi sottratto alla costante 0x5f3759df^[4][5][6]. Dopo la sottrazione fra interi e lo spostamento (shift) a destra, si ottiene un double che, se considerato a virgola mobile, diventa un'approssimazione grossolana della radice quadrata inversa del numero in ingresso. L'algoritmo termina con un'iterazione col metodo di Newton per guadagnare precisione sul risultato finale. L'algoritmo genera dei risultati ragionevolmente precisi usando un solo passaggio col metodo di Newton; rimane però più lento rispetto all'istruzione SSE rsqrtss del 1999 per microprocessori x86^[7].

Consideriamo il numero ${\displaystyle x=0{,}15625}$, per il quale si deve calcolare ${\displaystyle 1/{\sqrt {x}}\approx 2{,}52982}$. Seguono i primi passi dell'algoritmo:

0011_1110_0010_0000_0000_0000_0000_0000  Binario di x e i
0001_1111_0001_0000_0000_0000_0000_0000  Spostamento a destra di 1: (i >> 1)
0101_1111_0011_0111_0101_1001_1101_1111  Numero magico 0x5f3759df
0100_0000_0010_0111_0101_1001_1101_1111  Risultato della sottrazione: 0x5f3759df - (i >> 1)

Usando una rappresentazione IEEE 32-bit :

0_01111100_01000000000000000000000  1.25 * 2^-3
0_00111110_00100000000000000000000  1.125 * 2^-65
0_10111110_01101110101100111011111  1.432430... * 2^+63
0_10000000_01001110101100111011111  1.307430... * 2^+1

Interpretando l'ultimo numero binario come un numero a virgola mobile, si ottiene l'approssimazione ${\displaystyle y=2{,}61486}$ con un errore relativo di circa 3,4 %. Dopo l'applicazione di un'iterazione del metodo di Newton, il risultato finale è ${\displaystyle y=2{,}52549}$, con un errore di solo 0,17 %.

• Radice quadrata

• Wikimedia Commons contiene immagini o altri file sulla radice quadrata inversa veloce

Portale Informatica

Portale Matematica