Fix x86 sqrt rounding (bug 14032)

["Joseph S. Myers" <joseph at codesourcery dot com>]

This patch, relative to a tree with
<https://sourceware.org/ml/libc-alpha/2013-11/msg00809.html> (as fixed
by <https://sourceware.org/ml/libc-alpha/2013-11/msg00830.html>) and
<https://sourceware.org/ml/libc-alpha/2013-11/msg00829.html> applied,
fixes bug 14032, x86 sqrt not producing results correctly rounded to
double in some cases where the result rounded to long double is half
way between two double values, using the approach suggested in that
bug of adjusting the long double result using the C1 status bit (which
is set by fsqrt when the result was rounded up, cleared otherwise) in
the half-way cases.

This patch also ensures sqrt actually returns a result rounded to
double rather than one with excess precision.  Whether library
functions may return results with excess precision is a bit unclear
(C11 Annex F specifies that the return statement removes excess
precision, but there's no requirement for library functions to behave
as if written in C - hence for example the requirement for a sequence
point before a library function returns, even though library functions
written in C would have such a sequence point anyway); in my view they
shouldn't, and this is the first project listed under
<https://sourceware.org/glibc/wiki/Development_Todo/Master#libm_projects>,
but in any case it seems wrong to return a result with excess
precision, for a function that is meant to be correctly rounded, where
the result is not correctly rounded for the excess precision (such as
the adjusted value computed by sqrt internally to get the rounding to
double right).

w_sqrt.c used a version of __ieee754_sqrt from mathinline.h that's
conditional on __LIBC_INTERNAL_MATH_INLINES.  That's fine for most
uses of __ieee754_sqrt in libm (e.g. in complex functions) where
correct rounding doesn't matter, but inlining in w_sqrt.c would defeat
the adjustments required for correct rounding.  So w_sqrt.c is itself
wrapped by an x86-specific version that disables the problematic
__ieee754_sqrt to ensure the out-of-line version gets used there,
while functions not needing correct rounding of sqrt can continue to
get the faster inline version.

Tested x86 (and the new tests pass on x86_64).

[Large auto-libm-test-out diffs omitted below.]

2013-11-27  Joseph Myers  <joseph@codesourcery.com>

	[BZ #14032]
	* sysdeps/i386/fpu/e_sqrt.S (__ieee754_sqrt): Adjust half-way
	results and ensure return value rounded to double.
	* 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..21fe50e 100644
--- a/sysdeps/i386/fpu/e_sqrt.S
+++ b/sysdeps/i386/fpu/e_sqrt.S
@@ -8,6 +8,24 @@
 ENTRY(__ieee754_sqrt)
 	fldl	4(%esp)
 	fsqrt
+	fstsw	%ax
+	subl	$12, %esp
+	cfi_adjust_cfa_offset (12)
+	fld	%st
+	fstpt	(%esp)
+	movl	(%esp), %edx
+	andl	$0x7ff, %edx
+	cmpl	$0x400, %edx
+	jne	1f
+	andl	$0x200, %eax
+	subl	$0x100, %eax
+	subl	%eax, (%esp)
+	fstp	%st
+	fldt	(%esp)
+1:	fstpl	(%esp)
+	fldl	(%esp)
+	addl	$12,%esp
+	cfi_adjust_cfa_offset (-12)
 	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