Bug 52593 - Builtin sqrt on x86 is not correctly rounded
Summary: Builtin sqrt on x86 is not correctly rounded
Status: NEW
Alias: None
Product: gcc
Classification: Unclassified
Component: target (show other bugs)
Version: 4.6.4
: P3 normal
Target Milestone: ---
Assignee: Not yet assigned to anyone
URL:
Keywords: wrong-code
Depends on:
Blocks:
 
Reported: 2012-03-15 05:40 UTC by Rich Felker
Modified: 2014-03-10 22:36 UTC (History)
1 user (show)

See Also:
Host:
Target: i?86-*-*
Build:
Known to work:
Known to fail:
Last reconfirmed: 2012-03-15 00:00:00


Attachments

Note You need to log in before you can comment on or make changes to this bug.
Description Rich Felker 2012-03-15 05:40:34 UTC
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.
Comment 1 Richard Biener 2012-03-15 09:51:25 UTC
Which GCC version did you test?  Please provide a compilable testcase that can
be executed and shows the error.
Comment 2 Rich Felker 2012-03-15 16:37:25 UTC
Tested with gcc 4.6.2.

#include <stdio.h>
#include <math.h>
int main()
{
    volatile double x = 0x1.fffffffffffffp-1;
    volatile double y = sqrt(x);
    printf("%a\n", y);
}

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).
Comment 3 Dominique d'Humieres 2012-03-15 19:24:53 UTC
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'.
Comment 4 Rich Felker 2012-03-15 23:53:51 UTC
Of course. This bug is 387-math-specific.
Comment 5 Richard Biener 2012-03-16 10:53:12 UTC
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 :/
Comment 6 Rich Felker 2012-03-16 14:23:09 UTC
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...?
Comment 7 Rich Felker 2012-04-28 23:21:51 UTC
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.
Comment 8 joseph@codesourcery.com 2012-04-29 00:16:38 UTC
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.
Comment 9 Rich Felker 2012-04-29 01:21:59 UTC
Reported to glibc bug tracker as bug #14032:

http://sourceware.org/bugzilla/show_bug.cgi?id=14032
Comment 10 David Heidelberger (okias) 2014-03-10 22:36:20 UTC
fixed in glibc