¿Por qué SSE scalar sqrt (x) es más lento que rsqrt (x) * x?

He estado perfilando algunas de nuestras matemáticas básicas en un Intel Core Duo, y mientras observaba varios enfoques de raíz cuadrada, noté algo extraño: al usar las operaciones escalares SSE, es más rápido tomar una raíz cuadrada recíproca y multiplicarla para obtener el sqrt, que usar el código de operación sqrt nativo.

Lo estoy probando con un bucle algo como:

inline float TestSqrtFunction( float in );

void TestFunc()
{
  #define ARRAYSIZE 4096
  #define NUMITERS 16386
  float flIn[ ARRAYSIZE ]; // filled with random numbers ( 0 .. 2^22 )
  float flOut [ ARRAYSIZE ]; // filled with 0 to force fetch into L1 cache

  cyclecounter.Start();
  for ( int i = 0 ; i < NUMITERS ; ++i )
    for ( int j = 0 ; j < ARRAYSIZE ; ++j )
    {
       flOut[j] = TestSqrtFunction( flIn[j] );
       // unrolling this loop makes no difference -- I tested it.
    }
  cyclecounter.Stop();
  printf( "%d loops over %d floats took %.3f milliseconds",
          NUMITERS, ARRAYSIZE, cyclecounter.Milliseconds() );
}

Probé esto con algunos cuerpos diferentes para TestSqrtFunction, y tengo algunos tiempos que realmente me están rascando la cabeza. Lo peor de todo fue usar la función sqrt () nativa y dejar que el compilador "inteligente" se "optimice". A 24ns / float, usar el x87 FPU esto fue patéticamente malo:

inline float TestSqrtFunction( float in )
{  return sqrt(in); }

Lo siguiente que intenté fue usar un intrínseco para forzar al compilador a usar el código de operación sqrt escalar de SSE:

inline void SSESqrt( float * restrict pOut, float * restrict pIn )
{
   _mm_store_ss( pOut, _mm_sqrt_ss( _mm_load_ss( pIn ) ) );
   // compiles to movss, sqrtss, movss
}

Esto fue mejor, a 11,9 ns / flotación. También probé la extraña técnica de aproximación Newton-Raphson de Carmack , que funcionó incluso mejor que el hardware, a 4.3ns / float, aunque con un error de 1 en 2 ¹⁰ (que es demasiado para mis propósitos).

La sorpresa fue cuando probé la operación SSE para la raíz cuadrada recíproca , y luego usé una multiplicación para obtener la raíz cuadrada (x * 1 / √x = √x). Aunque esto elimina dos operaciones dependientes, que era la solución más rápida, con mucho, en 1.24ns / flotador y una precisión de 2 ^-14 :

inline void SSESqrt_Recip_Times_X( float * restrict pOut, float * restrict pIn )
{
   __m128 in = _mm_load_ss( pIn );
   _mm_store_ss( pOut, _mm_mul_ss( in, _mm_rsqrt_ss( in ) ) );
   // compiles to movss, movaps, rsqrtss, mulss, movss
}

Mi pregunta es básicamente ¿qué da ? ¿Por qué el código de operación de raíz cuadrada integrado en el hardware de SSE es más lento que sintetizarlo a partir de otras dos operaciones matemáticas?

Estoy seguro de que este es realmente el costo de la operación en sí, porque he verificado:

• Todos los datos caben en la caché y los accesos son secuenciales
• las funciones están en línea
• desenrollar el bucle no hace ninguna diferencia
• Los indicadores del compilador están configurados para la optimización completa (y el ensamblaje es bueno, lo verifiqué)

( editar : stephentyrone señala correctamente que las operaciones en cadenas largas de números deben usar las operaciones empaquetadas SIMD de vectorización, como rsqrtps, pero la estructura de datos de la matriz aquí es solo para fines de prueba: lo que realmente estoy tratando de medir es el rendimiento escalar para usar en el código que no se puede vectorizar).

sqrtssda un resultado correctamente redondeado. rsqrtssda una aproximación al recíproco, con una precisión de aproximadamente 11 bits.

sqrtssestá generando un resultado mucho más preciso, para cuando se requiere precisión. rsqrtssexiste para los casos en que una aproximación es suficiente, pero se requiere velocidad. Si lee la documentación de Intel, también encontrará una secuencia de instrucciones (aproximación recíproca de raíz cuadrada seguida de un solo paso de Newton-Raphson) que brinda una precisión casi total (~ 23 bits de precisión, si mal no recuerdo), y todavía es algo más rápido que sqrtss.

editar: Si la velocidad es crítica, y realmente está llamando a esto en un bucle para muchos valores, debería usar las versiones vectorizadas de estas instrucciones, rsqrtpso sqrtpsambas procesan cuatro flotantes por instrucción.

Esto también es válido para la división. MULSS (a, RCPSS (b)) es mucho más rápido que DIVSS (a, b). De hecho, es aún más rápido incluso cuando aumenta su precisión con una iteración Newton-Raphson.

Tanto Intel como AMD recomiendan esta técnica en sus manuales de optimización. En aplicaciones que no requieren el cumplimiento de IEEE-754, la única razón para usar div / sqrt es la legibilidad del código.

En lugar de proporcionar una respuesta, eso en realidad podría ser incorrecto (tampoco voy a verificar ni discutir sobre la caché y otras cosas, digamos que son idénticas), intentaré señalarle la fuente que puede responder a su pregunta.
La diferencia podría estar en cómo se calculan sqrt y rsqrt. Puede leer más aquí http://www.intel.com/products/processor/manuals/ . Sugeriría comenzar leyendo sobre las funciones del procesador que está usando, hay algo de información, especialmente sobre rsqrt (la CPU está usando una tabla de búsqueda interna con una gran aproximación, lo que hace que sea mucho más simple obtener el resultado). Puede parecer que rsqrt es mucho más rápido que sqrt, que 1 operación mul adicional (que no es demasiado costosa) podría no cambiar la situación aquí.

Editar: Algunos hechos que valdría la pena mencionar:
1. Una vez estaba haciendo algunas microoptimizaciones para mi biblioteca de gráficos y usé rsqrt para calcular la longitud de los vectores. (en lugar de sqrt, he multiplicado mi suma de cuadrados por rsqrt, que es exactamente lo que has hecho en tus pruebas), y funcionó mejor.
2. Calcular rsqrt usando una tabla de búsqueda simple podría ser más fácil, como para rsqrt, cuando x va al infinito, 1 / sqrt (x) va a 0, por lo que para las x pequeñas los valores de la función no cambian (mucho), mientras que para sqrt: va al infinito, por lo que es un caso simple;).

Además, una aclaración: no estoy seguro de dónde lo encontré en los libros que he vinculado, pero estoy bastante seguro de haber leído que rsqrt está usando una tabla de búsqueda, y solo debe usarse cuando el resultado no es necesario que sea exacto, aunque yo también podría estar equivocado, como lo fue hace algún tiempo :)

Newton-Raphson converge al cero de f(x)usar incrementos iguales a -f/f' donde f'está la derivada.

Para x=sqrt(y), se puede tratar de resolver f(x) = 0por xel uso f(x) = x^2 - y;

Entonces el incremento es: dx = -f/f' = 1/2 (x - y/x) = 1/2 (x^2 - y) / x que tiene una división lenta.

Puedes probar otras funciones (como f(x) = 1/y - 1/x^2) pero serán igualmente complicadas.

Veamos 1/sqrt(y)ahora. Puedes intentarlo f(x) = x^2 - 1/y, pero será igualmente complicado: dx = 2xy / (y*x^2 - 1)por ejemplo. Una opción alternativa no obvia para f(x)es:f(x) = y - 1/x^2

Luego: dx = -f/f' = (y - 1/x^2) / (2/x^3) = 1/2 * x * (1 - y * x^2)

¡Ah! No es una expresión trivial, pero solo tiene multiplicaciones, no divide. => ¡Más rápido!

Y: el paso de actualización completo new_x = x + dxluego dice:

x *= 3/2 - y/2 * x * x que también es fácil.

Hay varias otras respuestas a esto ya desde hace unos años. Esto es lo que acertó el consenso:

• Las instrucciones rsqrt * calculan una aproximación a la raíz cuadrada recíproca, buena para aproximadamente 11-12 bits.
• Se implementa con una tabla de búsqueda (es decir, una ROM) indexada por la mantisa. (De hecho, es una tabla de búsqueda comprimida, similar a las tablas matemáticas de antaño, que utiliza ajustes en los bits de orden inferior para ahorrar en transistores).
• La razón por la que está disponible es que es la estimación inicial utilizada por la FPU para el algoritmo de raíz cuadrada "real".
• También hay una instrucción recíproca aproximada, rcp. Ambas instrucciones son una pista de cómo la FPU implementa la raíz cuadrada y la división.

Esto es lo que se equivocó en el consenso:

• Las FPU de la era SSE no utilizan Newton-Raphson para calcular raíces cuadradas. Es un gran método en software, pero sería un error implementarlo de esa manera en hardware.

El algoritmo NR para calcular la raíz cuadrada recíproca tiene este paso de actualización, como otros han señalado:

x' = 0.5 * x * (3 - n*x*x);

Eso es un montón de multiplicaciones dependientes de datos y una resta.

Lo que sigue es el algoritmo que utilizan realmente las FPU modernas.

Dado b[0] = n, suponga que podemos encontrar una serie de números Y[i]tal que se b[n] = b[0] * Y[0]^2 * Y[1]^2 * ... * Y[n]^2acerque a 1. Luego considere:

x[n] = b[0] * Y[0] * Y[1] * ... * Y[n]
y[n] = Y[0] * Y[1] * ... * Y[n]

Es evidente que x[n]se acerca sqrt(n)y y[n]se acerca 1/sqrt(n).

Podemos usar el paso de actualización de Newton-Raphson para la raíz cuadrada recíproca para obtener un buen Y[i]:

b[i] = b[i-1] * Y[i-1]^2
Y[i] = 0.5 * (3 - b[i])

Luego:

x[0] = n Y[0]
x[i] = x[i-1] * Y[i]

y:

y[0] = Y[0]
y[i] = y[i-1] * Y[i]

La siguiente observación clave es esa b[i] = x[i-1] * y[i-1]. Entonces:

Y[i] = 0.5 * (3 - x[i-1] * y[i-1])
     = 1 + 0.5 * (1 - x[i-1] * y[i-1])

Luego:

x[i] = x[i-1] * (1 + 0.5 * (1 - x[i-1] * y[i-1]))
     = x[i-1] + x[i-1] * 0.5 * (1 - x[i-1] * y[i-1]))
y[i] = y[i-1] * (1 + 0.5 * (1 - x[i-1] * y[i-1]))
     = y[i-1] + y[i-1] * 0.5 * (1 - x[i-1] * y[i-1]))

Es decir, dados xey iniciales, podemos usar el siguiente paso de actualización:

r = 0.5 * (1 - x * y)
x' = x + x * r
y' = y + y * r

O, incluso más elegante, podemos configurar h = 0.5 * y. Esta es la inicialización:

Y = approx_rsqrt(n)
x = Y * n
h = Y * 0.5

Y este es el paso de actualización:

r = 0.5 - x * h
x' = x + x * r
h' = h + h * r

Este es el algoritmo de Goldschmidt, y tiene una gran ventaja si lo está implementando en hardware: el "bucle interno" es tres adiciones múltiples y nada más, y dos de ellos son independientes y se pueden canalizar.

En 1999, las FPU ya necesitaban un circuito de suma / resta en canalización y un circuito de multiplicación en canalización; de lo contrario, SSE no sería muy "fluido". En 1999, solo se necesitaba uno de cada circuito para implementar este bucle interno de una manera completamente canalizada sin desperdiciar mucho hardware solo en la raíz cuadrada.

Hoy, por supuesto, hemos fusionado la adición múltiple expuesta al programador. Una vez más, el bucle interno son tres FMA canalizados, que son (nuevamente) generalmente útiles incluso si no está calculando raíces cuadradas.

Es más rápido porque estas instrucciones ignoran los modos de redondeo y no manejan excepciones de punto flotante o números desnormalizados. Por estas razones, es mucho más fácil canalizar, especular y ejecutar otras instrucciones fp fuera de servicio.