Jump to content
Sign in to follow this  
fredericopissarra

Raiz quadrada em ponto flutuante. Já reparam nisso?

Recommended Posts

Eis uma simples função:

float f( float x ) { return sqrtf( x ); }

Já que os processadores desde o 486 têm o co-processador matemático incorporado, é de se supor que a rotina acima "use" apenas algumas instruções, não é? Mas não é isso o que acontece, mesmo na otimização máxima, obtemos algo assim:

f:
  pxor    xmm2, xmm2
  sqrtss  xmm1, xmm0
  ucomiss xmm2, xmm0
  ja      .L8
  movaps  xmm0, xmm1
  ret
.L8:
  sub     rsp, 24
  movss   dword [rsp+12], xmm1
  call    sqrtf
  movss   xmm1, dword [rsp+12]
  add     rsp, 24
  movaps  xmm0, xmm1
  ret

Essencialmente isso ai em cima faz chama a função sqrtf() se sqrtss retornar algo menor que zero (note que a lógica está invertida, UCOMISS compara zero [em xmm2] com o valor retornado por sqrtss), ou seja se CF=0 e ZF=0.

Ainda não entendi porque o GCC prefere fazer algo assim, já que sqrtss retorna os mesmos valores que sqrtf() para valores inválidos:

#include <stdio.h>
#include <x86intrin.h>
#include <math.h>

int main ( void )
{
  float a[] = { -0.0, -1.0, NAN, -NAN, INFINITY, -INFINITY };
  int i;
  __m128 f;

  for ( i = 0; i < sizeof a / sizeof a[0]; i++ )
  {
    f = _mm_set_ss ( &a[i] );

    printf ( "sqrtss(%1$f) = %2$f, sqrtf(%1$f) = %3$f\n", 
      a[i], 
      _mm_cvtss_f32 ( _mm_sqrt_ss( f ) ), 
      sqrtf ( a[i] ) );
  }
}

Compilando e executando, temos:

$ cc -O2 -o test test.c -lm
$ ./test
sqrtss(-0.000000) = -0.000000, sqrtf(-0.000000) = -0.000000
sqrtss(-1.000000) = -nan, sqrtf(-1.000000) = -nan
sqrtss(nan) = nan, sqrtf(nan) = nan
sqrtss(-nan) = -nan, sqrtf(-nan) = -nan
sqrtss(inf) = inf, sqrtf(inf) = inf
sqrtss(-inf) = -nan, sqrtf(-inf) = -nan

Até que eu descubra porquê o GCC faz esse teste, prefiro lidar com SSE diretamente para ter apenas uma instrução sendo executada. Só que temos um problema:

Na função abaixo, _mm_set_ss() não deveria fazer coisa alguma, já que o argumento x é um float e estará em XMM0 ela está ali apenas para compatibilizar os tipos float e __m128. No entanto, MOVSS zera os bits superiores do registrador  quando o operando fonte é memória. Logo depois _mm_sqrt_ss() extrairá a raiz quadrada e a colocará de volta em XMM0. A função _mm_cvtss_f32() garantidamente não faz nada aqui:

float f( float x ) { return _mm_cvtss_f32( _mm_sqrt_ss( _mm_set_ss( x ) ) ); }

Isso fica:

f:
  movss   dword [rsp-12], xmm0   ;
  movss   xmm0, dword [rsp-12]   ; ZERA os bits superiores.
  sqrtss  xmm0, xmm0
  ret

Não estranhe a escrita na pilha sem o devido deslocamento de RSP. Provavelmente o compilador está usando a red zone e escritas até 128 bytes abaixo de RSP são permitidas sem alocação.

Em essência, as duas primeiras instruções MOVSS são desnecessárias, mas usar algum macete do tipo abaixo não ajuda também:

__m128 f;
*(float *)&f = x;  // copia x para XMM0[0-31].

Isso acaba criando os mesmos dois MOVSS. E essa "zeração" dos bits superiores não estão ai porque o tipo de retorno é um escalar float simples. Mesmo mudando o retorno para __m128, os MOVSS continuam lá. Assim, a única forma de garantir uma única instrução é usando assembly inline. Embora o "constraint" de entrada fique estranho:

inline float mysqrtf( float x )
{
  __m128 r;

  __asm__ __volatile__ (
    "sqrtss %0,%1"  
    : "=x" (r) 
    : "x" (x)   // estranho, mas válido!
  );

  return _mm_cvtss_f32( r );
}

Agora sim, quando chamarmos mysqrtf() apenas a instrução SQRTSS será usada.

Ahhh... os mesmos testes são feitas, em sqrtf(), se -mfpmath=387 for usado...

  • Curtir 2

Share this post


Link to post
Share on other sites
Posted (edited)

Aliás, para o x86-64, esse último código é o usado pela rotina __ieee_sqrtf(), na glibc, exceto que o contraint de entrada é diferente ("xm").

Isso torna o código original ainda mais estranho... testa-se se x < 0 para chamar a mesma rotina (instrução)?!

Muito estranho...

[]s
Fred

Edited by fredericopissarra

Share this post


Link to post
Share on other sites

Confirmado! Isso é um BUG do GCC, pelo menos até a versão 7.3 (do Debian/Ubuntu, pelo menos) para o x86 e ARM AArch64 e versão 6.3 para ARM AArch32... Embora a versão 7.3 do MinGW-64 (do linux) não tenha esse BUG...
Não sei se foi resolvida na versão 8.

Felizmente é um bug inofensivo, já que o resultado será o da instrução SQRTSS (ou SQRTSD, se for double) e a chamda à função da libm sqrtf() (ou sqrt()) será feita, mas ignorada... Se for lidar com raiz quadrada no x86 e usa o GCC, prefira a função inline em asm, como mostrei acima. ou modifique-a para não ignorar o returno de sqrtf() se x < 0...

  • Curtir 1

Share this post


Link to post
Share on other sites

Testei, agorinha, com o GCC 8.2 e o bug ainda está lá... Testei com o CLANG 6 e o bug não está lá. Ele gera até um código mais elegante:

f:
  xorps xmm1,xmm1
  ucomiss xmm0,xmm1
  jb .lessthan0
  sqrtss xmm0,xmm0
  ret
.lessthan0:
  jmp sqrtf

Mas, essa chamada a sqrtf() é inútil, já que sqrtss retorna os mesmos valores para x < 0 e até valores inválidos como NAN e INFINITY.

  • Curtir 1

Share this post


Link to post
Share on other sites
Posted (edited)

Um detalhe: O código do CLANG não só é mais elegante, quanto mais rápido (se ambos chamarem sqrtf()). Saltos condicionais para frente são assumidos, pelo algoritmo do branch predictor, do processador, como não tomados. Como a função builtin sqrt() geralmente será chamada com valores válidos, o processador não perde tempo corrigindo o salto, que, na maioria das vezes, não será feito.

Mas, como eu disse, isso ai é um exagero dos compiladores... Um simples "SQRTSS XMM0,XMM0" resolve todo o problema...

Anyway... por favor, notem que não estou dizendo que clang tende a gerar sempre código melhor que o GCC. Pelos meus testes ele costuma criar código próximo da performance do gerado pelo GCC, mas, neste caso, ele ganhou a parada (embora pudesse ser melhor!).

Edited by fredericopissarra
  • Curtir 1

Share this post


Link to post
Share on other sites

Um follow-up no que concerne esse bug:

Não é um bug! Recebi uma explicação convincente do Bugzilla do GCC, onde me foi dito que a segunda chamada é feita apenas para ajustar o errno (errno = EDOM, quando x < 0). Mas, de qualquer modo, o código é ineficiente (comentário meu)...

Em essência, o GCC está fazendo o que deveria mesmo fazer.

Só mais um aviso: Não é necessário o uso de assembly inline para obter uma única instrução SQRTSS, SQRTSD ou FSQRT. Basta adicionar a opção -ffast-math (que não é habilitada com -O2 ou -O3 -- mas é habilitada com a opção -Ofast)...

  • Curtir 1

Share this post


Link to post
Share on other sites

Join the conversation

You can post now and register later. If you have an account, sign in now to post with your account.

Guest
Reply to this topic...

×   Pasted as rich text.   Paste as plain text instead

  Only 75 emoji are allowed.

×   Your link has been automatically embedded.   Display as a link instead

×   Your previous content has been restored.   Clear editor

×   You cannot paste images directly. Upload or insert images from URL.

Sign in to follow this  

  • Recently Browsing   0 members

    No registered users viewing this page.

×
×
  • Create New...