[Bug target/52593] New: Builtin sqrt on x86 is not correctly rounded

[Bug target/52593] New: Builtin sqrt on x86 is not correctly rounded

[bugdal at aerifal dot cx]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=52593

             Bug #: 52593
           Summary: Builtin sqrt on x86 is not correctly rounded
    Classification: Unclassified
           Product: gcc
           Version: unknown
            Status: UNCONFIRMED
          Severity: normal
          Priority: P3
         Component: target
        AssignedTo: unassigned@gcc.gnu.org
        ReportedBy: bugdal@aerifal.cx


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.

[Bug target/52593] Builtin sqrt on x86 is not correctly rounded

[rguenth at gcc dot gnu.org]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=52593

Richard Guenther <rguenth at gcc dot gnu.org> changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
             Status|UNCONFIRMED                 |WAITING
   Last reconfirmed|                            |2012-03-15
     Ever Confirmed|0                           |1

--- Comment #1 from Richard Guenther <rguenth at gcc dot gnu.org> 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.

[Bug target/52593] Builtin sqrt on x86 is not correctly rounded

[bugdal at aerifal dot cx]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=52593

--- Comment #2 from Rich Felker <bugdal at aerifal dot cx> 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).

[Bug target/52593] Builtin sqrt on x86 is not correctly rounded

[dominiq at lps dot ens.fr]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=52593

--- Comment #3 from Dominique d'Humieres <dominiq at lps dot ens.fr> 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'.

[Bug target/52593] Builtin sqrt on x86 is not correctly rounded

[bugdal at aerifal dot cx]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=52593

--- Comment #4 from Rich Felker <bugdal at aerifal dot cx> 2012-03-15 23:53:51 UTC ---
Of course. This bug is 387-math-specific.

[Bug target/52593] Builtin sqrt on x86 is not correctly rounded

[rguenth at gcc dot gnu.org]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=52593

Richard Guenther <rguenth at gcc dot gnu.org> changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
           Keywords|                            |wrong-code
             Target|                            |i?86-*-*
             Status|WAITING                     |NEW
            Version|unknown                     |4.6.4

--- Comment #5 from Richard Guenther <rguenth at gcc dot gnu.org> 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 :/

[Bug target/52593] Builtin sqrt on x86 is not correctly rounded

[bugdal at aerifal dot cx]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=52593

--- Comment #6 from Rich Felker <bugdal at aerifal dot cx> 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...?

[Bug target/52593] Builtin sqrt on x86 is not correctly rounded

[bugdal at aerifal dot cx]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=52593

--- Comment #7 from Rich Felker <bugdal at aerifal dot cx> 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.

[Bug target/52593] Builtin sqrt on x86 is not correctly rounded

[joseph at codesourcery dot com]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=52593

--- Comment #8 from joseph at codesourcery dot com <joseph at codesourcery dot 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.

[Bug target/52593] Builtin sqrt on x86 is not correctly rounded

[bugdal at aerifal dot cx]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=52593

--- Comment #9 from Rich Felker <bugdal at aerifal dot cx> 2012-04-29 01:21:59 UTC ---
Reported to glibc bug tracker as bug #14032:

http://sourceware.org/bugzilla/show_bug.cgi?id=14032

[Bug target/52593] Builtin sqrt on x86 is not correctly rounded

[david.heidelberger at ixit dot cz]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=52593

David Heidelberger (okias) <david.heidelberger at ixit dot cz> changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
                 CC|                            |david.heidelberger at ixit dot cz

--- Comment #10 from David Heidelberger (okias) <david.heidelberger at ixit dot cz> ---
fixed in glibc