This is the mail archive of the libc-alpha@sourceware.org mailing list for the glibc project.


Index Nav: [Date Index] [Subject Index] [Author Index] [Thread Index]
Message Nav: [Date Prev] [Date Next] [Thread Prev] [Thread Next]
Other format: [Raw text]

Re: Fix x86 sqrt rounding (bug 14032)


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


Index Nav: [Date Index] [Subject Index] [Author Index] [Thread Index]
Message Nav: [Date Prev] [Date Next] [Thread Prev] [Thread Next]