Raiz quadrada em ponto flutuante. Já reparam nisso?

[fredericopissarra]

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

[fredericopissarra]

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

[fredericopissarra]

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

[fredericopissarra]

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.

[fredericopissarra]

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

[fredericopissarra]

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

[fredericopissarra]

Outro follow-up: A geração de código ineficiente foi aceita como bug a ser corrigido (até o momento) no bugzilla do GCC.