The builtin sqrt() on x86 (i387) should be disabled except with -ffast-math because it is not correctly rounded. For example, sqrt(0x1.fffffffffffffp-1) yields 1 instead of 0x1.fffffffffffffp-1. Using -fno-builtin-sqrt will give the correct value assuming your C library/libm is correctly rounded.
Unfortunately bugs like this seem endemic in gcc. I would really like to see all dubious builtins and other dubious floating point optimizations disabled except with -ffast-math until somebody takes the time to rigorously test them and prove their correctness.
Which GCC version did you test? Please provide a compilable testcase that can
be executed and shows the error.
Tested with gcc 4.6.2.
volatile double x = 0x1.fffffffffffffp-1;
volatile double y = sqrt(x);
Compile with -O2 -ffloat-store, and this program gives an output of 0x1p+0, rather than the correct output of 0x1.fffffffffffffp-1. Unfortunately it's impossible to get the correct output with glibc's -lm version of sqrt either since it has the exact same bug. But you can look at the generated assembly with and without -fno-builtin-sqrt and see that the version with builtin sqrt is wrong (using the fsqrt opcode directly).
On x86_64-apple-darwin10 (default '-mfpmath=sse'), I get '0x1.fffffffffffffp-1' for all the revisions I have tested (from 4.4 to 4.8) unless I compile the test with '-mfpmath=387'.
Of course. This bug is 387-math-specific.
Confirmed. I think sqrt is special (compared to sin, cos, etc.) because it's
one of the core IEEE arithmetic functions. I suppose correct rounding
is only ensured for 80bit long double.
It will of course be an unexpected performance drop for most people with
no additional benefit as the libm implementation is wrong as well :/
The 387 FPU ensures correct rounding for the currently selected precision mode, which per the ABI is always extended precision.
As for the usefulness of fixing this, I found the bug while working on my correct sqrt implementation in musl libc, because despite the existence of a correct version of the function, I was still getting wrong results. It turned out gcc was replacing it with a buggy builtin. I don't think "glibc gets it wrong anyway" is a reason not to fix the problem, especially now that glibc seems to be under new maintainership and actually fixing longstanding WONTFIX bugs.
Folks who just care about speed and want to throw correctness away should already be using -ffast-math and similar.
Actually since this bug is rounding-related, perhaps it would suffice to make -frounding-math turn off the builtin sqrt when using 387 math...?
This bug seems to have been fixed with the addition of the -fexcess-precision=standard feature, which is now set by default with -std=c99 or c11, and which disables the builtin sqrt based on 387 fsqrt. So apparently it had already been fixed at the time I reported this, but I was unaware of the right options to enable the fix and did not even think to try just using -std=c99.
Note that for buggy libm (including glibc's), the fact that gcc has fixed the issue will not fix the incorrect results, since the code in libm makes exactly the same mistake gcc was making. But at least it's possible to fix it there.
If you have a bug in glibc's libm, please make sure there is an open bug
report for it in glibc Bugzilla, component "math"; I don't see anything
there about sqrt.
Reported to glibc bug tracker as bug #14032:
fixed in glibc