Todo lo que debes saber sobre el algoritmo de la inversa rápida de la raíz cuadrada

 (github.com/francisrstokes)
5 puntos por GN⁺ 2024-06-03 | 2 comentarios | Compartir por WhatsApp
  • El fast inverse square root, famoso por Quake 3, fue una solución de rendimiento de su época para aproximar rápidamente 1 / sqrt(x) mediante reinterpretación de bits de float y corrección con Newton-Raphson
  • La clave es que el patrón de bits entero de un float IEEE-754 de 32 bits puede tratarse como una aproximación escalada y desplazada de log2(x)
  • 0x5f3759df - (i >> 1) es una forma de trasladar log2(x^-0.5) = -0.5 * log2(x) a un desplazamiento de enteros y una resta, y la constante mágica proviene de 3/2 * 2^23 * (127 - σ)
  • Después se aplica una corrección de Newton-Raphson con una sola iteración de y = y * (1.5 - 0.5x * y * y), y la segunda iteración en el código de Quake está comentada
  • En 1999 se necesitaba la inversa de la raíz cuadrada cientos o miles de veces por segundo para iluminación y normalización de vectores 3D, pero en el hardware moderno la utilidad práctica de este truco ha disminuido gracias al procesamiento dedicado de punto flotante

Qué hace el código de Quake

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

  return y;
}
  • Esta función calcula un valor aproximado de la inversa de la raíz cuadrada 1 / sqrt(number) para number
  • La parte más famosa es la manipulación de bits que interpreta el valor float como si fuera un long y luego ejecuta 0x5f3759df - (i >> 1)
  • Cuando salió Quake 3 en 1999, la inversa de la raíz cuadrada era una operación lenta y costosa, y se necesitaba cientos o miles de veces por segundo en ecuaciones de iluminación y cálculos de vectores 3D que requerían normalización
  • En el hardware moderno, estos cálculos ya no se ejecutan en el CPU o, incluso cuando sí se ejecutan ahí, son rápidos gracias al hardware especializado de punto flotante que ha mejorado mucho

Representación float IEEE-754 de 32 bits

  • Un float de 32 bits está compuesto por tres partes
    • Sign: 1 bit, indica si es positivo o negativo
    • Exponent: 8 bits, define el rango al que pertenece el valor
    • Mantissa: 23 bits, representa linealmente la posición dentro de ese rango
  • Un valor normal se interpreta de la siguiente forma
N = (-1)^S * 2^(E - 127) * (1 + M / 2^23)
  • B = 127 es el valor de sesgo usado en el biased exponent, y el exponente real es e = E - B
  • La mantissa no se usa simplemente como multiplicación por m, sino en la forma 1 + m
    • Si m = 0, entonces es 2^e
    • Cuando m se acerca a 1, representa valores hasta justo antes del siguiente rango de exponente, 2^(e+1)
  • Si todos los bits del exponente son 0, se trata de un número sub-normal, y la fórmula cambia
N = (-1)^S * 2^-126 * m
  • Los sub-normal son necesarios para representar 0 y números muy pequeños, cercanos a 0
  • Si todos los bits del exponente son 1, se trata como un valor especial
    • Si E = 255 y M = 0, entonces es Infinity o -Infinity
    • Si M != 0, entonces es NaN

La relación logarítmica que aparece al ver los bits del float como entero

  • Si la representación interna del float se mira como un entero de 32 bits, puede expresarse con la siguiente fórmula
I_x = 2^31 S + 2^23 E + M
  • Como la inversa de la raíz cuadrada trabaja con entradas positivas, se toma S = 0, lo que simplifica la expresión
L = 2^23
I_x = L E + M
  • Dentro del mismo rango de exponente, la mantissa indica linealmente la posición, pero a medida que el exponente crece, la misma cantidad de pasos de mantissa cubre un tramo más amplio de la recta numérica
    • E = 127, es decir e = 0, corresponde aproximadamente al rango [1, 2)
    • E = 128, es decir e = 1, corresponde aproximadamente al rango [2, 4)
    • Ambos rangos tienen la misma cantidad de pasos de mantissa, pero el segundo es el doble de ancho
  • Por esta estructura, al ver el patrón bruto de bits del float como entero aparece una relación logarítmica

Los bits en bruto son una aproximación de log2(x)

  • Si el patrón de bits del float se interpreta como el entero I_x, puede verse como una aproximación lineal por tramos de log2(x)
  • Esta relación se expresa con la siguiente aproximación
log2(x) ≈ I_x / L - B
  • Si al entero bruto de bits se le divide por el tamaño de la mantissa L = 2^23 y se le resta el sesgo del exponente B = 127, se obtiene un valor cercano a log2(x)
  • El logaritmo dentro del intervalo de mantissa se trata como una aproximación lineal
log2(1 + x) ≈ x + σ
  • σ es un parámetro de ajuste que corrige la aproximación, y x representa la posición dentro del rango del exponente, en el intervalo [0, 1]

Convertir la inversa de la raíz cuadrada en una identidad logarítmica

  • El objetivo es obtener el siguiente valor
y = 1 / sqrt(x)
  • Esto puede reescribirse en forma exponencial como sigue
y = x^-0.5
  • Si se aplica la identidad logarítmica, el cálculo de la inversa de la raíz cuadrada queda en la siguiente relación
log2(1 / sqrt(x)) = log2(x^-0.5) = -0.5 * log2(x)
  • Aprovechando que los bits del float se comportan como una aproximación de log2(x), es posible aproximar directamente la representación entera de bits I_y para y a partir de la representación entera de bits I_x de x
I_y ≈ -0.5 I_x + 1.5 L (B - σ)
  • Esta fórmula lleva directamente a la línea clave del código de Quake
i = 0x5f3759df - ( i >> 1 );
  • i >> 1 desplaza los bits enteros una posición a la derecha y cumple el papel de multiplicar por 1/2
  • La constante de adelante, 0x5f3759df, corresponde a 1.5 * L * (B - σ)

Qué es la constante 0x5f3759df

  • Si se toma σ = 0, la constante se calcula así
1.5 * 2^23 * 127 = 1598029824
  • La representación hexadecimal de ese valor es 0x5f400000
  • Hay una diferencia de 566817 con respecto a la constante real de Quake, 0x5f3759df
  • A partir de esa diferencia, el valor de σ que corresponde al código de Quake se calcula de la siguiente manera
σ = 377878 / 2^23
σ = 0.04504656
  • En código C, puede calcularse la misma constante así
int32_t compute_magic(void) {
  double sigma = 0.0450465;
  double expression = 1.5 * pow(2.0, 23.0) * (127.0 - sigma);
  int32_t i = expression;
  return i;
}

// -> 0x5f3759df
  • Aquí se usa double, y la conversión a entero es un casting normal, no una reinterpretación de bits
  • Este valor de σ se eligió para optimizar la aproximación, pero no es el valor óptimo real, ni está claro quién lo creó exactamente

Por qué no es solo un hack simple

  • 0x5f3759df - (i >> 1) es una fórmula para construir un valor inicial de la inversa de la raíz cuadrada aprovechando que los bits brutos del float son una aproximación logarítmica
  • Aunque se basa en relaciones matemáticas complejas, en tiempo de ejecución solo usa operaciones rápidas como corrimientos y restas
  • En esa época había que procesar miles de operaciones costosas por segundo, así que este método fue un diseño de ingeniería ajustado a las limitaciones del hardware
  • Eso sí, este algoritmo solo funciona con normal float
    • En valores sub-normal no se cumple la premisa de la aproximación log2(1 + x) ≈ x + σ
    • En sub-normal, en la práctica aparece algo más cercano a 0 + x, así que la aproximación se rompe

Reducir el error con la corrección de Newton-Raphson

  • El valor inicial obtenido por manipulación de bits es bastante bueno, pero todavía deja un error medible
  • La siguiente línea mejora mucho la aproximación
y = y * ( threehalfs - ( x2 * y * y ) );
  • Esta línea aplica el método de Newton-Raphson
  • Para adaptar el problema de la inversa de la raíz cuadrada al método de Newton, se convierte en un problema de encontrar la raíz de la siguiente función
f(y) = 1 / y^2 - x = 0
  • El método de Newton genera una mejor aproximación y_(n+1) a partir de la aproximación actual y_n de la siguiente forma
y_(n+1) = y_n - f(y_n) / f'(y_n)
  • La derivada de f(y) = y^-2 - x es la siguiente
f'(y) = -2y^-3 = -2 / y^3

Fórmula de corrección de Newton sin divisiones

  • Si se usa la fórmula de Newton tal cual, aparecen varias divisiones en punto flotante
  • Una de las razones por las que este algoritmo es rápido es que evita la división en punto flotante
  • Si se reordena algebraicamente, queda una forma que usa solo multiplicaciones y no divisiones
y_(n+1) = y_n * (1.5 - 0.5x * y_n^2)
  • En el código de Quake, x2 = number * 0.5F precalcula 0.5x, y luego se usa en la siguiente línea
y = y * ( threehalfs - ( x2 * y * y ) );
  • Después de esta sola iteración, el error absoluto máximo es de 0.175%, y en muchos casos el error es menor
  • El código original incluye una segunda iteración de Newton, pero está comentada
// y  = y * ( threehalfs - ( x2 * y * y ) );   // 2nd iteration, this can be removed

Origen y algoritmos relacionados

  • Este algoritmo no fue inventado por John Carmack, y su origen exacto no es 100% seguro
  • Como referencia relacionada, se enlaza este artículo de Beyond3D: The truth is the exact origin is not 100% certain
  • Chris Lomont escribió un paper intentando encontrar el valor óptimo de sigma en la etapa de aproximación logarítmica: InvSqrt.pdf
  • CORDIC es un algoritmo que calcula seno y coseno sin punto flotante, usando solo sumas y desplazamientos de bits, y su funcionamiento detallado es bastante distinto al de fast inverse square root
  • Ambos algoritmos comparten la idea de aplicar observaciones matemáticas de manera eficiente según las limitaciones del hardware de su época

Lecturas relacionadas

2 comentarios

 
joyfui 2024-06-03

Ese código curioso que vuelve a aparecer de vez en cuando... jaja
 
GN⁺ 2024-06-03
Comentarios en Hacker News

  • Si es una computadora fabricada después de 1999, por lo general soporta el conjunto de instrucciones SSE, que incluye _mm_rsqrt_ps, capaz de calcular 4 raíces cuadradas inversas más rápido de una sola vez: https://www.intel.com/content/www/us/en/docs/intrinsics-guid...
    Aun así, la técnica tratada aquí no es completamente irrelevante todavía. La conversión float/int es rápida, pero todavía existe hardware sin instrucciones rsqrt, sqrt, pow ni log, y esas operaciones pueden aproximarse con este truco

    • La instrucción SSE de recíproco en punto flotante puede dar resultados ligeramente distintos entre Intel y AMD, así que si esperas resultados deterministas entre PCs puede volverse un dolor de cabeza: https://robert.ocallahan.org/2021/09/rr-trace-portability-di...
    • Curiosamente, SSE también tiene una instrucción de raíz cuadrada normal, pero es mucho más lenta que la de raíz cuadrada inversa, así que si puedes tolerar la pérdida de precisión, es más rápido calcular sqrt(x) como x * 1/sqrt(x)
    • En realidad, la inmensa mayoría de las computadoras ni siquiera soporta SSE, y mucho menos el conjunto de instrucciones i386/amd64; la proporción que lo soporta de forma no emulada sigue bajando
      En conjuntos de instrucciones de GPU, ARM, RISC-V, AVR, PIC, 8051, FPGA, etc., a menudo también viene integrada una operación aproximada de raíz cuadrada inversa, pero probablemente esté implementada con un algoritmo como este

  • Para ponerle un pequeño pero al artículo, no es correcto dar a entender que este tipo de cálculo no ocurre en los CPU modernos. Es un malentendido común pensar que los juegos o las apps con mucha carga de punto flotante quieren mandar todas las operaciones de punto flotante al GPU
    En realidad, solo tiene sentido mandarlo al GPU cuando se trata de trabajos grandes y uniformes. Si haces una normalización de vectores de una sola vez, como al construir una matriz de rotación para que un objeto mire hacia otro, normalmente es más rápido dejarlo en el CPU. Incluso ignorando el tiempo de transferencia al GPU, una sola operación de punto flotante suele ser más rápida en CPU, porque el GPU normalmente tiene menor reloj y alcanza altos FLOP gracias al paralelismo

    • Parece que aquí se referían no al GPU sino a la FPU. Antes la FPU hacía cálculos asíncronos, y ahora se considera una parte integrada del CPU

  • Escribí una implementación para MMIX, asumiendo originalmente que el valor de entrada es mayor que 2^-1021

  • Si te interesa, Wikipedia también tiene una buena explicación de esta función y su historia: https://en.wikipedia.org/wiki/Fast_inverse_square_root

  • Reuní varios de estos aquí: https://github.com/ncruces/fastmath/blob/main/fast.go
    También hay una publicación relacionada en StackOverflow: https://stackoverflow.com/questions/32042673/optimized-low-a...

    • Me resulta útil porque justo estaba pensando en empezar una recopilación de técnicas como estas para reutilizar un motor 3D estilo finales de los 80 que hice hace tiempo
    • También me gustaría ver los benchmarks del paquete fastmath

  • Ya es hora de ponerse quisquillosos. Hay un error tipográfico en la fórmula de float: debería ser (-1)^S y no -1^S. La segunda siempre da -1
    Tampoco es exacta la explicación de que interpretar el patrón de bits crudo sea una aproximación lineal por tramos de un logaritmo. Las líneas entre los puntos de datos de la gráfica azul no existen realmente, y no puedes tener solo la mitad de un bit en 1. Más bien se parece a una versión discreta del logaritmo, y los puntos de datos que sí existen, es decir, donde la línea roja y la azul se cruzan, son literalmente un log escalado y desplazado. Fuera de eso, es un buen artículo

    • No lo termino de entender. Si piensas en un float diminuto de 6 bits, con 1 bit de signo, 2 de exponente y 3 de mantisa, el rango [010000, 010111] contiene 2, 2.25, 2.5, 2.75, 3, 3.25, 3.5 y 3.75
      Pero las mantisas implícitas por el logaritmo base 2 de esos números son .0000000, .0010101, .0101001, .0111010, .1001010, .1011001, .1100111, .1110100, y salvo la primera no coinciden con 001, 010, etc. del float. Como los floats en el intervalo [2,4) están espaciados linealmente pero los logaritmos correspondientes no, sí puede verse, como dice el artículo, que el float es una aproximación lineal por tramos de un logaritmo
    • No es una aproximación lineal por tramos continua sino una aproximación lineal por tramos discreta. Es cierto que la línea azul no es continua, pero la interpretación es incorrecta. La gráfica azul no consiste solo en unos cuantos puntos de cruce, sino en 256 puntos individuales distribuidos uniformemente a lo largo del eje x
      Si fuera la gráfica completa habría 2^32 posibilidades dentro del patrón lineal por tramos, pero eso no es lo que dibuja el texto original. Dado que el artículo trata operaciones con enteros de 32 bits y floats IEEE-754 de 32 bits, me parece aceptable omitir lo de “discreta” en la explicación

  • Es un buen artículo que explica muchos conceptos interesantes, pero la manipulación algebraica de una sección es sorprendentemente mala
    Después de “hay muchos pasos exactos para pasar de la primera forma a esta, pero los incluí todos por completitud”, el desarrollo tiene pasos innecesarios y varios errores de signo que se cancelan entre sí. En particular, al pasar de la segunda a la tercera línea no se distribuye correctamente el signo negativo. A partir de la segunda línea, se podría ir mucho más corto y de forma correcta, empezando con y_n+1 = y_n + (1 - x * y_n^2) / y_n^2 * (y_n^3 / 2) y llegando a y_n+1 = y_n (1.5 * y_n - 0.5 * x * y_n * y_n). Me parece que los pasos intermedios son obvios para cualquiera que entienda álgebra

  • El número mágico de ese famoso fragmento de código no es la constante óptima. Usando otra constante probablemente se puede reducir el error relativo como en un 0.5% más
    En ese tiempo quizá era difícil encontrar el valor absolutamente óptimo, pero ahora es relativamente fácil. Yo también me metí alguna vez en esa madriguera y tengo un cuaderno de Jupyter para encontrar los números mágicos óptimos de (1/x^2) y (1/x)

    • El enlace a un artículo que explora esa pregunta está al final del texto

  • Lo más interesante de este artículo para mí fue el enlace a “How Java's Floating-Point Hurts Everyone Everywhere”: https://people.eecs.berkeley.edu/~wkahan/JAVAhurt.pdf
    El autor es William Kahan, también conocido como “Old Man of Floating-Point”: https://news.ycombinator.com/item?id=29042853 - An Interview with the Old Man of Floating-Point (1998)

    • No tiene que ver con el tema, pero empecé a leer el PDF de JAVAhurt y la tipografía es espantosa. Parece que usaron un paquete de TeX que abre demasiado y de forma muy irregular el espacio entre palabras, y casi da la impresión de que hicieron OCR a otro documento y se colaron espacios de más
      Incluso en las partes con fuente monoespaciada hay espacios extra raros. Me costó muchísimo concentrarme al leerlo y, aunque sé que no lo es, casi se siente como un manifiesto de científico loco

  • Este video que vi hace tiempo me gustó muchísimo: https://www.youtube.com/watch?v=p8u_k2LIZyo