public inbox for glibc-bugs@sourceware.org
help / color / mirror / Atom feed
* [Bug math/14032] New: sqrt is not correctly rounded on x86
@ 2012-04-29  1:21 bugdal at aerifal dot cx
  2012-04-29 17:00 ` [Bug math/14032] " ppluzhnikov at google dot com
                   ` (4 more replies)
  0 siblings, 5 replies; 6+ messages in thread
From: bugdal at aerifal dot cx @ 2012-04-29  1:21 UTC (permalink / raw)
  To: glibc-bugs

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

             Bug #: 14032
           Summary: sqrt is not correctly rounded on x86
           Product: glibc
           Version: unspecified
            Status: NEW
          Severity: normal
          Priority: P2
         Component: math
        AssignedTo: unassigned@sourceware.org
        ReportedBy: bugdal@aerifal.cx
    Classification: Unclassified


Created attachment 6381
  --> http://sourceware.org/bugzilla/attachment.cgi?id=6381
test case and example workaround code

glibc's 387 asm for sqrt() simply uses the fsqrt instruction, yielding results
that may not be correctly rounded, which violates the IEEE requirements for the
sqrt operation. The cause is double-rounding due to the fact that the fsqrt
operation is performed in extended 80-bit precision (and gives
correctly-rounded results at this precision), then rounded again when it's
converted to 64-bit double precision. Naturally this is very problematic for
any algorithm that depends on the result of sqrt being correctly rounded.

Until recently, the bug was almost always masked by the corresponding bug in
gcc (which generates fsqrt for __builtin_sqrt), gcc bug #52593:

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

But since -fexcess-precision=standard was added (and made default for -std=c99
and -std=c11) somewhere in the gcc 4.5 series, the library function will be
called when correct standards-conformance is desired, only to repeat the exact
same mistake in the library code.

Using the C sqrt code would fix the problem, but I think that's highly
undesirable. I have an algorithm we're using in musl that uses fsqrt and only
fixes up the result in cases where the double-rounding could have an effect on
the result; it uses the floating point status word to determine the direction
of the rounding that takes place and biases the 80-bit value so that the result
will be correctly rounded to double precision after store and load:

http://git.etalabs.net/cgi-bin/gitweb.cgi?p=musl;a=blob;f=src/math/i386/sqrt.s

Attached is test code by Szabolcs Nagy (nsz) which demonstrates incorrectly
rounded results and includes a C function d_sqrt() which patches up the results
of an incorrectly-rounded underlying sqrt.

Finally, note that sqrtf() is not affected. It can be shown (left as an
exercise for the reader) that the correctly-rounded 80-bit result of fsqrt
subsequently rounded to 32-bit single precision is also the correctly-rounded
single-precision square root.

-- 
Configure bugmail: http://sourceware.org/bugzilla/userprefs.cgi?tab=email
------- You are receiving this mail because: -------
You are on the CC list for the bug.


^ permalink raw reply	[flat|nested] 6+ messages in thread
end of thread, other threads:[~2014-06-13 10:49 UTC | newest]

Thread overview: 6+ messages (download: mbox.gz / follow: Atom feed)
-- links below jump to the message on this page --
2012-04-29  1:21 [Bug math/14032] New: sqrt is not correctly rounded on x86 bugdal at aerifal dot cx
2012-04-29 17:00 ` [Bug math/14032] " ppluzhnikov at google dot com
2013-10-13  6:30 ` bugdal at aerifal dot cx
2013-11-29 16:31 ` cvs-commit at gcc dot gnu.org
2013-11-29 16:34 ` jsm28 at gcc dot gnu.org
2014-06-13 10:49 ` fweimer at redhat dot com

This is a public inbox, see mirroring instructions
for how to clone and mirror all data and code used for this inbox;
as well as URLs for read-only IMAP folder(s) and NNTP newsgroup(s).