Por que SSE escalar sqrt (x) é mais lento que rsqrt (x) * x?

Eu estive traçando um perfil de alguns de nossos cálculos matemáticos básicos em um Intel Core Duo e, ao examinar várias abordagens para a raiz quadrada, notei algo estranho: usando as operações escalares SSE, é mais rápido obter uma raiz quadrada recíproca e multiplicá-la para obter o sqrt, do que usar o opcode nativo sqrt!

Estou testando com um loop parecido com:

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

Eu tentei isso com alguns corpos diferentes para TestSqrtFunction, e eu tenho alguns tempos que estão realmente coçando minha cabeça. O pior de tudo, de longe, foi usar a função nativa sqrt () e deixar o compilador "inteligente" "otimizar". A 24 ns / flutuante, usar o FPU x87 era pateticamente ruim:

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

A próxima coisa que tentei foi usar um intrínseco para forçar o compilador a usar o opcode sqrt escalar do 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
}

Isso era melhor, a 11,9 ns / float. Eu também tentei a técnica de aproximação Newton-Raphson maluca do Carmack , que funcionou ainda melhor do que o hardware, a 4,3 ns / float, embora com um erro de 1 em 2 ¹⁰ (o que é demais para meus propósitos).

A surpresa foi quando tentei a operação de SSE para raiz quadrada recíproca e, em seguida, usei uma multiplicação para obter a raiz quadrada (x * 1 / √x = √x). Mesmo que isso leve duas operações dependentes, foi a solução mais rápida de longe, a 1,24 ns / flutuante e com precisão 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
}

Minha pergunta é basicamente o que dá ? Por que o opcode de raiz quadrada embutido no hardware do SSE é mais lento do que sintetizá-lo a partir de duas outras operações matemáticas?

Tenho certeza de que este é realmente o custo da operação em si, porque verifiquei:

• Todos os dados cabem no cache e os acessos são sequenciais
• as funções são embutidas
• desenrolar o loop não faz diferença
• os sinalizadores do compilador estão definidos para otimização total (e a montagem está boa, eu verifiquei)

( editar : stephentyrone aponta corretamente que as operações em longas sequências de números devem usar as operações empacotadas SIMD de vetorização, como rsqrtps- mas a estrutura de dados da matriz aqui é apenas para fins de teste: o que estou realmente tentando medir é o desempenho escalar para uso no código que não pode ser vetorizado.)

sqrtssdá um resultado arredondado corretamente. rsqrtssdá uma aproximação ao recíproco, com precisão de cerca de 11 bits.

sqrtssestá gerando um resultado muito mais preciso, para quando a precisão é necessária. rsqrtssexiste para os casos em que uma aproximação é suficiente, mas a velocidade é necessária. Se você ler a documentação da Intel, também encontrará uma sequência de instruções (aproximação de raiz quadrada recíproca seguida por uma única etapa de Newton-Raphson) que fornece precisão quase total (~ 23 bits de precisão, se bem me lembro), e ainda é um pouco mais rápido do que sqrtss.

editar: Se a velocidade é crítica, e você está realmente chamando isso em um loop para muitos valores, você deve usar as versões vetorizadas dessas instruções, rsqrtpsou sqrtpsambas processam quatro flutuadores por instrução.

Isso também se aplica à divisão. MULSS (a, RCPSS (b)) é muito mais rápido do que DIVSS (a, b). Na verdade, é ainda mais rápido, mesmo quando você aumenta sua precisão com uma iteração de Newton-Raphson.

A Intel e a AMD recomendam essa técnica em seus manuais de otimização. Em aplicativos que não exigem conformidade com IEEE-754, o único motivo para usar div / sqrt é a legibilidade do código.

Em vez de fornecer uma resposta, isso pode estar incorreto (também não vou verificar ou discutir sobre cache e outras coisas, digamos que sejam idênticos), tentarei apontar a fonte que pode responder à sua pergunta.
A diferença pode estar em como sqrt e rsqrt são calculados. Você pode ler mais aqui http://www.intel.com/products/processor/manuals/ . Eu sugiro começar lendo sobre as funções do processador que você está usando, há algumas informações, especialmente sobre rsqrt (cpu está usando a tabela de pesquisa interna com grande aproximação, o que torna muito mais simples obter o resultado). Pode parecer que rsqrt é muito mais rápido que sqrt, que 1 operação mul adicional (que não é tão cara) pode não mudar a situação aqui.

Edit: Alguns fatos que podem valer a pena mencionar:
1. Uma vez eu estava fazendo algumas microotimizações para minha biblioteca de gráficos e usei rsqrt para calcular o comprimento dos vetores. (em vez de sqrt, multipliquei minha soma do quadrado pelo rsqrt dele, que é exatamente o que você fez em seus testes), e o desempenho foi melhor.
2. Calcular rsqrt usando a tabela de pesquisa simples pode ser mais fácil, pois para rsqrt, quando x vai para o infinito, 1 / sqrt (x) vai para 0, então para x pequenos os valores da função não mudam (muito), enquanto para sqrt - vai ao infinito, então é um caso simples;).

Além disso, esclarecimento: não tenho certeza de onde encontrei nos livros que vinculei, mas tenho quase certeza de que li que rsqrt está usando alguma tabela de pesquisa e deve ser usada apenas quando o resultado não precisa ser exato, embora - eu também possa estar errado, como estava há algum tempo :).

Newton-Raphson converge para o zero de f(x)usar incrementos iguais a -f/f' onde f'está a derivada.

Pois x=sqrt(y), você pode tentar resolver f(x) = 0por xusar f(x) = x^2 - y;

Então o incremento é: dx = -f/f' = 1/2 (x - y/x) = 1/2 (x^2 - y) / x que tem uma divisão lenta nele.

Você pode tentar outras funções (como f(x) = 1/y - 1/x^2), mas elas serão igualmente complicadas.

Vamos dar uma olhada 1/sqrt(y)agora. Você pode tentar f(x) = x^2 - 1/y, mas será igualmente complicado: dx = 2xy / (y*x^2 - 1)por exemplo. Uma escolha alternativa não óbvia para f(x)é:f(x) = y - 1/x^2

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

Ah! Não é uma expressão trivial, mas você só tem multiplicadores nela, sem divisão. => Mais rápido!

E: a etapa de atualização completa new_x = x + dxentão diz:

x *= 3/2 - y/2 * x * x o que também é fácil.

Existem várias outras respostas para isso já de alguns anos atrás. Aqui está o que o consenso acertou:

• As instruções rsqrt * calculam uma aproximação da raiz quadrada recíproca, boa para cerca de 11-12 bits.
• É implementado com uma tabela de pesquisa (ou seja, uma ROM) indexada pela mantissa. (Na verdade, é uma tabela de pesquisa compactada, semelhante às tabelas matemáticas antigas, usando ajustes nos bits de ordem inferior para economizar nos transistores.)
• A razão pela qual está disponível é que é a estimativa inicial usada pela FPU para o algoritmo de raiz quadrada "real".
• Também há uma instrução recíproca aproximada, rcp. Ambas as instruções são uma pista de como a FPU implementa a raiz quadrada e a divisão.

Aqui está o que o consenso errou:

• FPUs da era SSE não usam Newton-Raphson para calcular raízes quadradas. É um ótimo método em software, mas seria um erro implementá-lo dessa forma em hardware.

O algoritmo NR para calcular a raiz quadrada recíproca tem esta etapa de atualização, como outros observaram:

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

São muitas multiplicações dependentes de dados e uma subtração.

O que se segue é o algoritmo que as FPUs modernas realmente usam.

Dado b[0] = n, suponha que possamos encontrar uma série de números Y[i]que se b[n] = b[0] * Y[0]^2 * Y[1]^2 * ... * Y[n]^2aproxime de 1. Em seguida, considere:

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

Claramente x[n]se aproxima sqrt(n)e y[n]se aproxima 1/sqrt(n).

Podemos usar a etapa de atualização Newton-Raphson para raiz quadrada recíproca para obter um bom Y[i]:

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

Então:

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

e:

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

A próxima observação chave é esta b[i] = x[i-1] * y[i-1]. Assim:

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

Então:

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

Ou seja, dados xey iniciais, podemos usar a seguinte etapa de atualização:

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

Ou, ainda mais sofisticado, podemos definir h = 0.5 * y. Esta é a inicialização:

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

E esta é a etapa de atualização:

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

Este é o algoritmo de Goldschmidt, e tem uma grande vantagem se você o estiver implementando em hardware: o "loop interno" é três multiplicação-adições e nada mais, e dois deles são independentes e podem ser pipeline.

Em 1999, as FPUs já precisavam de um circuito de adição / substrato em pipeline e um circuito de multiplicação em pipeline, caso contrário, o SSE não seria muito "streaming". Apenas um de cada circuito foi necessário em 1999 para implementar este loop interno de uma forma totalmente pipeline, sem desperdiçar muito hardware apenas na raiz quadrada.

Hoje, é claro, fundimos multiplicação-adição exposta ao programador. Novamente, o loop interno são três FMAs em pipeline, que são (novamente) geralmente úteis mesmo se você não estiver computando raízes quadradas.

É mais rápido porque essas instruções ignoram os modos de arredondamento e não lidam com exceções de ponto flutuante ou números desnormalizados. Por essas razões, é muito mais fácil pipeline, especular e executar outra instrução fp fora de ordem.