This is the mail archive of the
libc-alpha@sourceware.org
mailing list for the glibc project.
Re: Fix x86 sqrt rounding (bug 14032)
- From: "Joseph S. Myers" <joseph at codesourcery dot com>
- To: Richard Henderson <rth at twiddle dot net>
- Cc: <libc-alpha at sourceware dot org>
- Date: Wed, 27 Nov 2013 23:13:25 +0000
- Subject: Re: Fix x86 sqrt rounding (bug 14032)
- Authentication-results: sourceware.org; auth=none
- References: <Pine dot LNX dot 4 dot 64 dot 1311271803540 dot 7837 at digraph dot polyomino dot org dot uk> <52966555 dot 20603 at twiddle dot net>
On Thu, 28 Nov 2013, Richard Henderson wrote:
> Surely this is much more complex than simply setting the rounding precision to
> double before using fsqrt.
Here's a patch using that approach. Tested x86.
2013-11-27 Joseph Myers <joseph@codesourcery.com>
[BZ #14032]
* sysdeps/i386/fpu/e_sqrt.S (__ieee754_sqrt): Do fsqrt with
precision control set to double precision.
* sysdeps/i386/fpu/w_sqrt.c: New file.
* math/auto-libm-test-in: Add more tests.
* math/auto-libm-test-out: Update.
diff --git a/math/auto-libm-test-in b/math/auto-libm-test-in
index 5dc790c..12cb27d 100644
--- a/math/auto-libm-test-in
+++ b/math/auto-libm-test-in
@@ -25,3 +25,60 @@ sqrt 0.25
sqrt 6642.25
sqrt 15190.5625
sqrt 0.75
+sqrt 0x1.fffffffffffffp+1023
+sqrt 0x1.ffffffffffffbp+1023
+sqrt 0x1.ffffffffffff7p+1023
+sqrt 0x1.ffffffffffff3p+1023
+sqrt 0x1.fffffffffffefp+1023
+sqrt 0x1.fffffffffffebp+1023
+sqrt 0x1.fffffffffffe7p+1023
+sqrt 0x1.fffffffffffe3p+1023
+sqrt 0x1.fffffffffffdfp+1023
+sqrt 0x1.fffffffffffdbp+1023
+sqrt 0x1.fffffffffffd7p+1023
+sqrt 0x1.0000000000003p-1022
+sqrt 0x1.0000000000007p-1022
+sqrt 0x1.000000000000bp-1022
+sqrt 0x1.000000000000fp-1022
+sqrt 0x1.0000000000013p-1022
+sqrt 0x1.0000000000017p-1022
+sqrt 0x1.000000000001bp-1022
+sqrt 0x1.000000000001fp-1022
+sqrt 0x1.0000000000023p-1022
+sqrt 0x1.0000000000027p-1022
+sqrt 0x1.000000000002bp-1022
+sqrt 0x1.000000000002fp-1022
+sqrt 0x1.0000000000033p-1022
+sqrt 0x1.0000000000037p-1022
+sqrt 0x1.7167bc36eaa3bp+6
+sqrt 0x1.7570994273ad7p+6
+sqrt 0x1.7dae969442fe6p+6
+sqrt 0x1.7f8444fcf67e5p+6
+sqrt 0x1.8364650e63a54p+6
+sqrt 0x1.85bedd274edd8p+6
+sqrt 0x1.8609cf496ab77p+6
+sqrt 0x1.873849c70a375p+6
+sqrt 0x1.8919c962cbaaep+6
+sqrt 0x1.8de4493e22dc6p+6
+sqrt 0x1.924829a17a288p+6
+sqrt 0x1.92702cd992f12p+6
+sqrt 0x1.92b763a8311fdp+6
+sqrt 0x1.947da013c7293p+6
+sqrt 0x1.9536091c494d2p+6
+sqrt 0x1.61b04c6p-1019
+sqrt 0x1.93789f1p-1018
+sqrt 0x1.a1989b4p-1018
+sqrt 0x1.f93bc9p-1018
+sqrt 0x1.2f675e3p-1017
+sqrt 0x1.a158508p-1017
+sqrt 0x1.cd31f078p-1017
+sqrt 0x1.33b43b08p-1016
+sqrt 0x1.6e66a858p-1016
+sqrt 0x1.8661cbf8p-1016
+sqrt 0x1.bbb221b4p-1016
+sqrt 0x1.c4942f3cp-1016
+sqrt 0x1.dbb258c8p-1016
+sqrt 0x1.57103ea4p-1015
+sqrt 0x1.9b294f88p-1015
+sqrt 0x1.0000000000001p+0
+sqrt 0x1.fffffffffffffp-1
diff --git a/sysdeps/i386/fpu/e_sqrt.S b/sysdeps/i386/fpu/e_sqrt.S
index 1054ba4..fba5833 100644
--- a/sysdeps/i386/fpu/e_sqrt.S
+++ b/sysdeps/i386/fpu/e_sqrt.S
@@ -7,7 +7,17 @@
ENTRY(__ieee754_sqrt)
fldl 4(%esp)
+ subl $8, %esp
+ cfi_adjust_cfa_offset (8)
+ fstcw 4(%esp)
+ movl $0xfeff, %edx
+ andl 4(%esp), %edx
+ movl %edx, (%esp)
+ fldcw (%esp)
fsqrt
+ fldcw 4(%esp)
+ addl $8, %esp
+ cfi_adjust_cfa_offset (-8)
ret
END (__ieee754_sqrt)
strong_alias (__ieee754_sqrt, __sqrt_finite)
diff --git a/sysdeps/i386/fpu/w_sqrt.c b/sysdeps/i386/fpu/w_sqrt.c
new file mode 100644
index 0000000..19b5074
--- /dev/null
+++ b/sysdeps/i386/fpu/w_sqrt.c
@@ -0,0 +1,8 @@
+/* The inline __ieee754_sqrt is not correctly rounding; it's OK for
+ most internal uses in glibc, but not for sqrt itself. */
+#define __ieee754_sqrt __avoid_ieee754_sqrt
+#include <math.h>
+#include <math_private.h>
+#undef __ieee754_sqrt
+extern double __ieee754_sqrt (double);
+#include <math/w_sqrt.c>
--
Joseph S. Myers
joseph@codesourcery.com