Raíz cuadrada inversa inusual de John Carmack (Quake III)

Question

John Carmack tiene una función especial en el código fuente de Quake III que calcula la raíz cuadrada inversa de un flotador, 4 veces más rápido que (float)(1.0/sqrt(x)) normal, incluida una constante 0x5f3759df constante. Vea el código a continuación. ¿Puede alguien explicar línea por línea qué está sucediendo exactamente aquí y por qué esto funciona mucho más rápido que la implementación normal?

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

    x2 = number * 0.5F; 
    y = number; 
    i = * (long *) &y; 
    i = 0x5f3759df - (i >> 1); 
    y = * (float *) &i; 
    y = y * (threehalfs - (x2 * y * y)); 

    #ifndef Q3_VM 
    #ifdef __linux__ 
    assert(!isnan(y)); 
    #endif 
    #endif 
    return y; 
} 

Answer 1

FYI. Carmack no lo escribió. Terje Mathisen y Gary Tarolli tienen un crédito parcial (y muy modesto) por ello, y acreditan otras fuentes.

Cómo se derivó la constante mítica es algo así como un misterio.

Para citar Gary Tarolli:

Que en realidad está haciendo un cálculo de coma flotante en número entero - tomó un largo tiempo para averiguar cómo y por qué esto funciona, y no puedo recordar el detalles más.

Un poco mejor constante, desarrollado por un experto matemático (Chris Lomont) tratando de averiguar cómo el algoritmo original trabajado es:

float InvSqrt(float x) 
{ 
    float xhalf = 0.5f * x; 
    int i = *(int*)&x;    // get bits for floating value 
    i = 0x5f375a86 - (i >> 1);  // gives initial guess y0 
    x = *(float*)&i;    // convert bits back to float 
    x = x * (1.5f - xhalf * x * x); // Newton step, repeating increases accuracy 
    return x; 
} 

A pesar de esto, su primer intento de una forma matemática 'superior La versión de sqrt de id (que llegó a ser casi la misma constante) demostró ser inferior a la desarrollada inicialmente por Gary a pesar de ser matemáticamente mucho más "pura". No podía explicar por qué Id's era tan excelente iirc.

Answer 2

Según to this nice article escrito hace un tiempo ...

La magia del código, incluso si no puede seguirlo, se destaca como el i = 0x5f3759df - (i >> 1) ; línea. Simplificado, Newton-Raphson es una aproximación que comienza con una conjetura y lo refina con iteración. Tomando ventaja de la naturaleza de x86 procesadores de 32 bits, i, un entero, está establece inicialmente en el valor del número de punto flotante desea tomar la inversa del cuadrado de, utilizando un molde entero. I se establece en 0x5f3759df, menos se desplazó un poco a la derecha. El cambio a la derecha deja caer el bit menos significativo de i, esencialmente reduciéndolo a la mitad.

Es una muy buena lectura. Esto es solo una pequeña parte de eso.

Answer 3

Por supuesto, actualmente es mucho más lento que usar un sqrt de FPU (especialmente en 360/PS3), porque el intercambio entre float e int induce una carga-hit-store, mientras que la unidad de coma flotante puede hacer raíz cuadrada recíproca en el hardware.

Muestra cómo deben evolucionar las optimizaciones a medida que cambia la naturaleza del hardware subyacente.

Answer 4

Greg Hewgill y IllidanS4 dio un enlace con una excelente explicación matemática. Trataré de resumirlo aquí para aquellos que no quieren profundizar en los detalles.

Cualquier función matemática, con algunas excepciones, puede ser representado por una suma polinómica:

y = f(x) 

puede ser exactamente transforma en:

y = a0 + a1*x + a2*(x^2) + a3*(x^3) + a4*(x^4) + ... 

Dónde a0, a1, a2 ,. .. son constantes. El problema es que para muchas funciones, como la raíz cuadrada, para el valor exacto, esta suma tiene un número infinito de miembros, no termina en unos x^n. Pero, si nos detenemos en x^n todavía tendríamos un resultado hasta cierta precisión.

lo tanto, si tenemos:

y = 1/sqrt(x) 

En este caso particular se decidió descartar todos los miembros de polinomios anteriores segundos, probablemente debido a la velocidad de cálculo:

y = a0 + a1*x + [...discarded...] 

Y la tarea ahora tiene vino hacia abajo para calcular a0 y a1 para que y tenga la menor diferencia con respecto al valor exacto. Ellos han calculado que los valores más adecuados son:

a0 = 0x5f375a86 
a1 = -0.5 

Así que cuando se pone esto en la ecuación que se obtiene:

y = 0x5f375a86 - 0.5*x 

¿Qué es la misma que la línea que se ve en el código:

i = 0x5f375a86 - (i >> 1); 

Editar: en realidad aquí y = 0x5f375a86 - 0.5*x no es lo mismo que i = 0x5f375a86 - (i >> 1); ya que el desplazamiento de flotación como un entero no solo se divide por dos sino que también divide el exponente por dos y causa algunos otros artefactos, pero aún se reduce al cálculo de algunos coeficientes a0, a1, a2 ....

En este punto han descubierto que la precisión de este resultado no es suficiente para este propósito. Así que, además, lo hicieron sólo un paso de la iteración de Newton para mejorar la precisión de los resultados:

x = x * (1.5f - xhalf * x * x) 

podrían haber hecho un poco más iteraciones de un bucle, cada una mejora resultado, hasta que se cumpla la precisión requerida. ¡Así funciona exactamente en CPU/FPU! Pero parece que solo una iteración fue suficiente, lo que también fue una bendición para la velocidad. CPU/FPU realiza tantas iteraciones como sea necesario para alcanzar la precisión para el número de punto flotante en el que se almacena el resultado y tiene un algoritmo más general que funciona para todos los casos.

---

Así que en resumen, lo que hicieron es:

uso (casi) el mismo algoritmo que la CPU/FPU, explotar la mejora de las condiciones iniciales para el caso especial de 1/sqrt (x) y no calcule todo el camino hasta que la CPU/FPU de precisión vaya a detenerse antes, ganando velocidad de cálculo.

Answer 5

Tenía curiosidad por ver cuál era la constante como un flotador, así que simplemente escribí este fragmento de código y busqué en Google el número entero que apareció.

long i = 0x5F3759DF; 
    float* fp = (float*)&i; 
    printf("(2^127)^(1/2) = %f\n", *fp); 
    //Output 
    //(2^127)^(1/2) = 13211836172961054720.000000 

Parece que la constante es "Una aproximación número entero a la raíz cuadrada de 2^127 más conocido por el formato hexadecimal de su representación de coma flotante, 0x5f3759df" https://mrob.com/pub/math/numbers-18.html

En el mismo sitio se explica todo. https://mrob.com/pub/math/numbers-16.html#le009_16