optimize reciprocal square root with fast-math (x86)

[rotateright]

Bugzilla Link	20900
Version	trunk
OS	All
Blocks	#31672
CC	@hfinkel,@RKSimon,@tycho

Extended Description

$ ./clang -v
clang version 3.6.0 (217530)
Target: x86_64-apple-darwin13.3.0
Thread model: posix

$ cat rsqrt.c
#include <math.h>
float reciprocal_square_root(float x) {
return 1.0f / sqrtf(x);
}

$ ./clang -O2 -ffast-math -S -o - rsqrt.c
...
sqrtss %xmm0, %xmm1
movss LCPI0_0(%rip), %xmm0
divss %xmm1, %xmm0

---

This should be optimized to use 'rsqrtss'.

ICC 14 does this at -O2:

    rsqrtss   %xmm0, %xmm2
    mulss     %xmm2, %xmm0
    mulss     %xmm2, %xmm0
    movss     L_2il0floatpacket.2(%rip), %xmm1
    mulss     %xmm1, %xmm2
    subss     L_2il0floatpacket.1(%rip), %xmm0
    mulss     %xmm2, %xmm0
    ret       

L_2il0floatpacket.1:
.long 0x40400000
L_2il0floatpacket.2:
.long 0xbf000000

[rotateright]

assigned to @rotateright

[hfinkel]

In case it helps, the PowerPC backend has code in PPCISelLowering that does this reasonably.

[rotateright]

Very cool. Thanks, Hal!

[rotateright]

gcc 4.9 with fast-math produces:
rsqrtss %xmm0, %xmm1
mulss %xmm1, %xmm0
mulss %xmm1, %xmm0
addss LC0(%rip), %xmm0
mulss LC1(%rip), %xmm1
mulss %xmm1, %xmm0
ret
LC0:
.long 3225419776 ## 0xC0400000
LC1:
.long 3204448256 ## 0xBF000000

[rotateright]

Just for reference - and relates to the discussion here:
http://reviews.llvm.org/D5222

I was temporarily thwarted in generating the optimized PPC code because I wasn't explicitly passing "-enable-unsafe-fp-math" to llc as a parameter. I thought the 'unsafe' function attribute on this IR would be picked up, but that's not enough to trigger the PPC optimization:

; Function Attrs: nounwind readnone ssp uwtable
define float @​reciprocal_square_root(float %x) #​0 {
entry:
%0 = tail call float @​llvm.sqrt.f32(float %x)
%div = fdiv fast float 1.000000e+00, %0
ret float %div
}

; Function Attrs: nounwind readnone
declare float @​llvm.sqrt.f32(float) #​1

attributes #​0 = { nounwind readnone ssp uwtable "less-precise-fpmad"="false" "no-frame-pointer-elim"="true" "no-frame-pointer-elim-non-leaf" "no-infs-fp-math"="true" "no-nans-fp-math"="true" "stack-protector-buffer-size"="8" "unsafe-fp-math"="true" "use-soft-float"="false" }
attributes #​1 = { nounwind readnone }

[tycho]

This should also apply to the ARM backend. I wrote a patch in February which adapts the PowerPC reciprocal square root lowering to the ARM backend, but nothing happened with it:

http://lists.cs.uiuc.edu/pipermail/llvm-commits/Week-of-Mon-20140224/206860.html

[tycho]

This should also apply to the ARM backend. I wrote a patch in February which
adapts the PowerPC reciprocal square root lowering to the ARM backend, but
nothing happened with it:

http://lists.cs.uiuc.edu/pipermail/llvm-commits/Week-of-Mon-20140224/206860.
html

Er, sorry, the above was for X86, and could apply for this bug. The ARM one was not what I submitted, though I have a patch somewhere locally to implement it...

[rotateright]

Thanks, Steven! And I agree with what Hal suggested then - we should be able to refactor some of this logic into a common area. I'll start working on it this week.

[RKSimon]

http://reviews.llvm.org/D5370

I recently posted this patch to improve X86 scheduling for rsqrt instructions - SSE rsqrt instructions are current classed the same as sqrt, meaning that the scheduler assumes latencies of >20 cycles. The patch splits them off into their own class to drastically improves scheduled latency acccuracy (often <5 cycles).

The equivalent SSE rcp instruction has already been split off from div in a similar manner.

[rotateright]

Preliminary patch (just refactoring) for review here:
http://reviews.llvm.org/D5425

[rotateright]

Even more refactoring here:
http://reviews.llvm.org/D5484

My hope is that this will allow any backend to control generation of the estimate algorithm using just 2 knobs:

1. Existence of an estimate instruction for a given data type
2. Number of N-R steps needed to refine the estimate to something useful

[rotateright]

We need to do some benchmarking. Hal's rsqrt implementation for PPC trades a constant FP load for an additional FP subtract operation. This is different than what gcc and icc generate and also different from what Steven proposed in his X86 patch.

The difference arises from implementing the Newton-Raphson refinement as:
(est * est * X - 3.0f) * est * -0.5 ; use 2 constants
((est * est * (X - (1.5 * X)) + 1.5) * est ; use 1 constant

As AVX asm, the difference looks like this:

$ cat rsqrt_gcc.s
...
vrsqrtss %xmm0, %xmm1, %xmm1
vmulss %xmm0, %xmm1, %xmm0
vmulss %xmm1, %xmm0, %xmm0
vaddss LC0(%rip), %xmm0, %xmm0
vmulss LC1(%rip), %xmm1, %xmm1
vmulss %xmm1, %xmm0, %xmm0
ret

$ cat rsqrt_llvm.s
...
vmovss LCPI0_0(%rip), %xmm1
vmulss %xmm1, %xmm0, %xmm2
vsubss %xmm2, %xmm0, %xmm2
vrsqrtss %xmm0, %xmm0, %xmm0
vmulss %xmm0, %xmm0, %xmm3
vmulss %xmm3, %xmm2, %xmm2
vaddss %xmm1, %xmm2, %xmm1
vmulss %xmm1, %xmm0, %xmm0
retq

Here's a rough test harness for checking the speed of each:

$ cat test_rsqrt.c
#include <stdio.h>
#define ITERATIONS (2 * 1100 * 1000 * 1000ULL)

extern float foo(float x);

int main() {
float x = 1.0f;

unsigned long long i;
for (i=0; i<ITERATIONS; i++)
x = foo(x);

printf("x = %f\n", x);

return 0;
}

$ cat rsqrt.c
#include <math.h>
float foo(float x) {
return 1.0f/sqrtf(x);
}

After enabling the rsqrt codegen for x86, I tested on a 2.2GHz SandyBridge and got:

$ gcc -O1 -ffast-math -fomit-frame-pointer -march=corei7-avx rsqrt.c -S -o rsqrt_gcc.s

$ clang -O1 -ffast-math -fomit-frame-pointer -march=corei7-avx rsqrt.c -S -o rsqrt_llvm.s

$ clang -O1 test_rsqrt.c rsqrt_llvm.s -march=corei7-avx -o rsqrt_llvm

$ clang -O1 test_rsqrt.c rsqrt_gcc.s -march=corei7-avx -o rsqrt_gcc

$ time ./rsqrt_gcc
x = 1.000000

real 0m15.657s
user 0m15.645s
sys 0m0.009s

$ time ./rsqrt_llvm
x = 1.000000

real 0m16.239s
user 0m16.229s
sys 0m0.007s

[rotateright]

The 2 constant implementation is also faster on AMD Jaguar 1.5GHz:

$ time ./rsqrt_gcc
x = 1.000000

real 0m16.212s
user 0m16.175s
sys 0m0.028s

$ time ./rsqrt_llvm
x = 1.000000

real 0m17.687s
user 0m17.647s
sys 0m0.030s

[rotateright]

For reference, the current -ffast-math clang codegen that uses full-precision but slow sqrtss and divss running on the same 1.5GHz AMD Jaguar system:

~$ time ./rsqrt_llvm_nottoofast
x = 1.000000

real 1m22.929s
user 1m22.615s
sys 0m0.215s

So using the estimate is about 5x faster.

[rotateright]

Hal, is there a particular PPC system that you benchmarked that was faster using the single constant version of the N-R refinement? Or did you use a different test loop? My test is attempting to measure the average latency of a single execution of the rsqrt path by making the result dependent on the previous iteration.

Note: I've applied for access to https://gcc.gnu.org/wiki/CompileFarm (thx, Bill Schmidt!) for more testability, but have not received a reply.