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

floating-point square-root

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 el normal (float)(1.0/sqrt(x)) , incluida una extraña constante 0x5f3759df . 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

De acuerdo con este lindo artículo escrito hace un tiempo ...

La magia del código, incluso si no puedes seguirlo, se destaca como i = 0x5f3759df - (i >> 1); línea. Simplificado, Newton-Raphson es una aproximación que comienza con una suposición y la refina con iteración. Aprovechando la naturaleza de los procesadores x86 de 32 bits, i, un entero, se establece inicialmente en el valor del número de punto flotante del que desea tomar el cuadrado inverso, utilizando un molde de enteros. Luego se establece en 0x5f3759df, menos se desplazó un bit a la derecha. El cambio a la derecha arroja el bit menos significativo de i, esencialmente a la mitad.

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

Answer 2

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 a Gary Tarolli:

Que en realidad está haciendo un cálculo en coma flotante en número entero, tomó mucho tiempo descubrir cómo y por qué funciona, y ya no recuerdo los detalles.

Una constante un poco mejor, desarrollada por un matemático experto (Chris Lomont) tratando de averiguar cómo funcionaba el algoritmo original:

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 intento inicial de una versión matemáticamente "superior" de sqrt de ID (que llegó a ser casi la misma constante) resultó inferior al desarrollado inicialmente por Gary a pesar de ser matemáticamente mucho más "puro". No podía explicar por qué Id''s era tan excelente iirc.

Answer 3

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

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

Answer 4

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 entera a la raíz cuadrada de 2 ^ 127 más conocida por la forma hexadecimal de su representación en coma flotante, 0x5f3759df" https://mrob.com/pub/math/numbers-18.html

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

Answer 5

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

Cualquier función matemática, con algunas excepciones, puede representarse mediante una suma polinómica:

y = f(x)

puede transformarse exactamente en:

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

Donde 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 algún x ^ n . Pero, si nos detenemos en algo x ^ n , todavía tendríamos un resultado hasta cierta precisión.

Entonces, si tenemos:

y = 1/sqrt(x)

En este caso particular, decidieron descartar todos los miembros polinomiales por encima del segundo, probablemente debido a la velocidad de cálculo:

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

Y ahora la tarea se ha reducido para calcular a0 y a1 para que y tenga la menor diferencia con respecto al valor exacto. Han calculado que los valores más apropiados son:

a0 = 0x5f375a86
a1 = -0.5

Entonces cuando pones esto en ecuación obtienes:

y = 0x5f375a86 - 0.5*x

Que es lo mismo que la línea que ves 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 la coma 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, descubrieron que la precisión de este resultado no es suficiente para este propósito. Así que, además, hicieron solo un paso de la iteración de Newton para mejorar la precisión del resultado:

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

Podrían haber hecho más iteraciones en un bucle, cada una mejorando el 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.

En resumen, lo que hicieron fue:

Utilice (casi) el mismo algoritmo que CPU / FPU, explote 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, por lo tanto ganando velocidad de cálculo