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

106

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 10 (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.)

Crashworks
fonte
13
x / sqrt (x) = sqrt (x). Ou, dito de outra forma: x ^ 1 * x ^ (- 1/2) = x ^ (1 - 1/2) = x ^ (1/2) = sqrt (x)
Crashworks
6
claro inline float SSESqrt( float restrict fIn ) { float fOut; _mm_store_ss( &fOut, _mm_sqrt_ss( _mm_load_ss( &fIn ) ) ); return fOut; },. Mas esta é uma má ideia porque pode facilmente induzir um bloqueio de carga-acerto-armazenamento se a CPU grava os flutuadores na pilha e os lê de volta imediatamente - fazendo malabarismo do registrador vetorial para um registrador flutuante para o valor de retorno em particular são más notícias. Além disso, os opcodes subjacentes da máquina que os intrínsecos SSE representam levam operandos de endereço de qualquer maneira.
Crashworks de
4
O quanto LHS importa depende da geração particular e da etapa de um determinado x86: minha experiência é que em qualquer coisa até i7, mover dados entre conjuntos de registros (por exemplo, FPU para SSE para eax) é muito ruim, enquanto uma viagem de ida e volta entre xmm0 e pilha e de volta não, por causa do encaminhamento da Intel para as lojas. Você mesmo pode cronometrar para ver com certeza. Geralmente, a maneira mais fácil de ver o potencial LHS é olhar o conjunto emitido e ver onde os dados são manipulados entre os conjuntos de registros; seu compilador pode fazer a coisa certa ou não. Quanto aos vetores normalizando, escrevi os meus resultados aqui: bit.ly/9W5zoU
Crashworks
2
Para o PowerPC, sim: a IBM tem um simulador de CPU que pode prever LHS e muitas outras bolhas de pipeline por meio de análise estática. Alguns PPCs também têm um contador de hardware para LHS que você pode consultar. É mais difícil para o x86; boas ferramentas de criação de perfil são mais escassas (o VTune está um tanto quebrado atualmente) e os pipelines reordenados são menos determinísticos. Você pode tentar medi-lo empiricamente medindo instruções por ciclo, o que pode ser feito precisamente com os contadores de desempenho de hardware. Os registros de "instruções retiradas" e "ciclos totais" podem ser lidos com, por exemplo, PAPI ou PerfSuite ( bit.ly/an6cMt ).
Crashworks de
2
Você também pode simplesmente escrever algumas permutações em uma função e cronometrá-las para ver se alguma sofre particularmente de paralisações. A Intel não publica muitos detalhes sobre a forma como seus pipelines funcionam (que eles LHS em tudo é uma espécie de segredo sujo), então muito do que aprendi foi olhando para um cenário que causa uma paralisação em outros archs (por exemplo, PPC ) e, em seguida, construir um experimento controlado para ver se o x86 também tem.
Crashworks de

Respostas:

216

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.

Stephen Canon
fonte
3
A etapa n / r oferece 22 bits de precisão (ela a duplica); 23 bits seriam exatamente a precisão total.
Jasper Bekkers
7
@Jasper Bekkers: Não, não faria. Primeiro, float tem 24 bits de precisão. Em segundo lugar, sqrtssestá corretamente arredondado , o que requer ~ 50 bits antes do arredondamento e não pode ser obtido usando uma iteração N / R simples com precisão única.
Stephen Canon
1
Esse é definitivamente o motivo. Para estender este resultado: o projeto Embree da Intel ( software.intel.com/en-us/articles/… ), usa a vetorização para sua matemática. Você pode baixar a fonte nesse link e ver como eles fazem seus vetores 3/4 D. A normalização do vetor deles usa rsqrt seguido por uma iteração de newton-raphson, que é muito preciso e ainda mais rápido do que 1 / ssqrt!
Brandon Pelfrey
7
Uma pequena advertência: x rsqrt (x) resulta em NaN se x for zero ou infinito. 0 * rsqrt (0) = 0 * INF = NaN. INF rsqrt (INF) = INF * 0 = NaN. Por esse motivo, o CUDA nas GPUs NVIDIA calcula raízes quadradas de precisão única aproximadas como receitas (rsqrt (x)), com o hardware fornecendo uma aproximação rápida da raiz quadrada recíproca e recíproca. Obviamente, verificações explícitas que tratam dos dois casos especiais também são possíveis (mas seriam mais lentas na GPU).
njuffa
@BrandonPelfrey Em qual arquivo você encontrou a etapa do Newton Rhapson?
fredoverflow
7

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.

Spat
fonte
1
Broadwell e posteriores têm melhor desempenho de divisão de FP, então compiladores como o clang optam por não usar recíproco + Newton para escalar em CPUs recentes, porque geralmente não é mais rápido. Na maioria dos loops, divnão é a única operação, então a taxa de transferência total do uop costuma ser o gargalo, mesmo quando há um divpsou divss. Consulte Divisão de ponto flutuante vs multiplicação de ponto flutuante , onde minha resposta contém uma seção sobre por que rcppsnão é mais uma vitória de rendimento. (Ou um ganho de latência) e números na divisão de taxa de transferência / latência.
Peter Cordes,
Se seus requisitos de precisão são tão baixos que você pode pular uma iteração de Newton, então sim a * rcpss(b)pode ser mais rápido, mas ainda é mais uops do que a/b!
Peter Cordes,
5

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 :).

Marcin Deptuła
fonte
4

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.

skal
fonte
2

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.

Pseudônimo
fonte
1
Relacionado: Como sqrt () do GCC funciona depois de compilado? Qual método de root é usado? Newton-Raphson? tem alguns links para designs de unidade de execução de div / sqrt de hardware. Rsqrt vetorizado rápido e recíproco com SSE / AVX dependendo da precisão - uma iteração de Newton no software, com ou sem FMA, para uso com _mm256_rsqrt_ps, com análise de desempenho Haswell. Normalmente, apenas uma boa ideia se você não tiver outro trabalho no loop e isso prejudicaria fortemente a taxa de transferência do divisor. HW sqrt é único uop, então está ok misturado com outro trabalho.
Peter Cordes
-2

É 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.

Witek
fonte
Obviamente errado. FMA depende do modo de arredondamento atual, mas tem uma taxa de transferência de dois por relógio no Haswell e posterior. Com duas unidades FMA totalmente tubuladas, Haswell pode ter até 10 FMAs em vôo ao mesmo tempo. A resposta certa é rsqrt's muito menor precisão, o que significa muito menos trabalho a fazer (ou nenhum?) Depois de uma tabela de lookup para obter uma estimativa inicial.
Peter Cordes