diff mbox

Fix sysdeps/ieee754 pow handling of sNaN arguments (bug 20916) [committed]

Message ID alpine.DEB.2.20.1612022321370.1078@digraph.polyomino.org.uk
State New
Headers show

Commit Message

Joseph Myers Dec. 2, 2016, 11:21 p.m. UTC
Various pow function implementations mishandle sNaN arguments in
various ways.  This includes returning sNaN instead of qNaN for sNaN
arguments.  For arguments (1, sNaN) and (sNaN, 0), TS 18661-1
semantics are also that the result should be qNaN, whereas with a qNaN
argument there the result should be 1, but for the dbl-64
implementation of pow there are issues with sNaN arguments beyond not
implementing the TS 18661-1 semantics in those special cases.

This patch makes the implementations in sysdeps/ieee754 follow the TS
18661-1 semantics consistently.  Because x86 / x86_64 implementations
still need fixing, testcases are not included with this patch; they
will be included with the fix for the x86 / x86_64 versions.

Tested for x86_64, x86, mips64 and powerpc (with such testcases, which
pass in the mips64 and powerpc cases).  Committed.

2016-12-02  Joseph Myers  <joseph@codesourcery.com>

	[BZ #20916]
	* sysdeps/ieee754/dbl-64/e_pow.c (__ieee754_pow): Do not return 1
	for arguments (sNaN, 0) or (1, sNaN).  Do arithmetic on NaN
	arguments to compute result.
	* sysdeps/ieee754/flt-32/e_powf.c (__ieee754_powf): Do not return
	1 for arguments (sNaN, 0) or (1, sNaN).
	* sysdeps/ieee754/ldbl-128/e_powl.c (__ieee754_powl): Likewise.
	* sysdeps/ieee754/ldbl-128ibm/e_powl.c (__ieee754_powl): Likewise.

Comments

Stefan Liebler Dec. 12, 2016, 4:43 p.m. UTC | #1
On 12/03/2016 12:21 AM, Joseph Myers wrote:
> Various pow function implementations mishandle sNaN arguments in
> various ways.  This includes returning sNaN instead of qNaN for sNaN
> arguments.  For arguments (1, sNaN) and (sNaN, 0), TS 18661-1
> semantics are also that the result should be qNaN, whereas with a qNaN
> argument there the result should be 1, but for the dbl-64
> implementation of pow there are issues with sNaN arguments beyond not
> implementing the TS 18661-1 semantics in those special cases.
>
> This patch makes the implementations in sysdeps/ieee754 follow the TS
> 18661-1 semantics consistently.  Because x86 / x86_64 implementations
> still need fixing, testcases are not included with this patch; they
> will be included with the fix for the x86 / x86_64 versions.
>
> Tested for x86_64, x86, mips64 and powerpc (with such testcases, which
> pass in the mips64 and powerpc cases).  Committed.

Hi Joseph,
on s390, I get test fails in math/test-double / math/test-idouble:
testing double (without inline functions)
Failure: Test: pow (1, sNaN)
Result:
  is:          1.0000000000000000e+00   0x1.0000000000000p+0
  should be:  qNaN
Failure: Test: pow (1, -sNaN)
Result:
  is:          1.0000000000000000e+00   0x1.0000000000000p+0
  should be:  qNaN
Failure: Test: pow_downward (1, sNaN)
Result:
  is:          1.0000000000000000e+00   0x1.0000000000000p+0
  should be:  qNaN
Failure: Test: pow_downward (1, -sNaN)
Result:
  is:          1.0000000000000000e+00   0x1.0000000000000p+0
  should be:  qNaN
Failure: Test: pow_towardzero (1, sNaN)
Result:
  is:          1.0000000000000000e+00   0x1.0000000000000p+0
  should be:  qNaN
Failure: Test: pow_towardzero (1, -sNaN)
Result:
  is:          1.0000000000000000e+00   0x1.0000000000000p+0
  should be:  qNaN
Failure: Test: pow_upward (1, sNaN)
Result:
  is:          1.0000000000000000e+00   0x1.0000000000000p+0
  should be:  qNaN
Failure: Test: pow_upward (1, -sNaN)
Result:
  is:          1.0000000000000000e+00   0x1.0000000000000p+0
  should be:  qNaN

Test suite completed:
   97258 test cases plus 86314 tests for exception flags and
     86314 tests for errno executed.
   8 errors occurred.


I've debugged in sysdeps/ieee754/dbl-64/e_pow.c and recognized that sNaN 
is converted to qNaN in line 85:
       if (y == 0)
	return 1.0;
This comparison is done with a load-and-test instruction from and to the 
same register, which results in a qNaN.
This value is passed to issignaling (y) in line 148:
   if (qy >= 0x7ff00000 && (qy > 0x7ff00000 || v.i[LOW_HALF] != 0))	/* 
NaN */
     return x == 1.0 && !issignaling (y) ? 1.0 : y + y;


 From ieee 754-2008 "6.2 Operations with NaNs":
"Under default exception handling, any operation signaling an invalid 
operation exception and for which a floating-point result is to be 
delivered shall deliver a quiet NaN."

As the used load-and-test instruction delivers a result, I think qNaN is 
correct. But is the compiler allowed to use this instruction for a 
comparision against zero?

How can I get rid of this issue?
Shall I use "if ((v.i[HIGH_HALF] & 0x7fffffff) == 0x0)" on s390?

Bye
Stefan
Joseph Myers Dec. 13, 2016, 9:42 p.m. UTC | #2
On Mon, 12 Dec 2016, Stefan Liebler wrote:

> I've debugged in sysdeps/ieee754/dbl-64/e_pow.c and recognized that sNaN is
> converted to qNaN in line 85:
>       if (y == 0)
> 	return 1.0;
> This comparison is done with a load-and-test instruction from and to the same
> register, which results in a qNaN.
> This value is passed to issignaling (y) in line 148:
>   if (qy >= 0x7ff00000 && (qy > 0x7ff00000 || v.i[LOW_HALF] != 0))	/* NaN
> */
>     return x == 1.0 && !issignaling (y) ? 1.0 : y + y;
> 
> 
> From ieee 754-2008 "6.2 Operations with NaNs":
> "Under default exception handling, any operation signaling an invalid
> operation exception and for which a floating-point result is to be delivered
> shall deliver a quiet NaN."
> 
> As the used load-and-test instruction delivers a result, I think qNaN is
> correct. But is the compiler allowed to use this instruction for a comparision
> against zero?

TS 18661-1 allows for floating-point assignments and argument passing 
(etc.) to act as convertFormat operations, so converting signaling NaNs to 
quiet.

Obviously when implemented that way, issignaling cannot work reliably.  As 
a quality-of-implementation issue it's probably best to avoid such 
instructions for simple loads of stored values, at least when 
-fsignaling-nans is in use.

So if -fsignaling-nans avoids the issue, compiling pow with 
-fsignaling-nans on s390 would make sense (most of libm isn't built with 
-fsignaling-nans and should still work fine).  Otherwise the tests of 
sNaNs can be disabled in a math-tests.h file for your architecture (like 
sysdeps/i386/fpu/math-tests.h), which should have a comment pointing to a 
bug report in GCC Bugzilla about the issue.
Andreas Krebbel Dec. 14, 2016, 9:53 a.m. UTC | #3
On Tue, Dec 13, 2016 at 09:42:14PM +0000, Joseph Myers wrote:
> On Mon, 12 Dec 2016, Stefan Liebler wrote:
> 
> > I've debugged in sysdeps/ieee754/dbl-64/e_pow.c and recognized that sNaN is
> > converted to qNaN in line 85:
> >       if (y == 0)
> > 	return 1.0;
> > This comparison is done with a load-and-test instruction from and to the same
> > register, which results in a qNaN.
> > This value is passed to issignaling (y) in line 148:
> >   if (qy >= 0x7ff00000 && (qy > 0x7ff00000 || v.i[LOW_HALF] != 0))	/* NaN
> > */
> >     return x == 1.0 && !issignaling (y) ? 1.0 : y + y;
> > 
> > 
> > From ieee 754-2008 "6.2 Operations with NaNs":
> > "Under default exception handling, any operation signaling an invalid
> > operation exception and for which a floating-point result is to be delivered
> > shall deliver a quiet NaN."
> > 
> > As the used load-and-test instruction delivers a result, I think qNaN is
> > correct. But is the compiler allowed to use this instruction for a comparision
> > against zero?
> 
> TS 18661-1 allows for floating-point assignments and argument passing 
> (etc.) to act as convertFormat operations, so converting signaling NaNs to 
> quiet.
> 
> Obviously when implemented that way, issignaling cannot work reliably.  As 
> a quality-of-implementation issue it's probably best to avoid such 
> instructions for simple loads of stored values, at least when 
> -fsignaling-nans is in use.
> 
> So if -fsignaling-nans avoids the issue, compiling pow with 
> -fsignaling-nans on s390 would make sense (most of libm isn't built with 
> -fsignaling-nans and should still work fine).  Otherwise the tests of 
> sNaNs can be disabled in a math-tests.h file for your architecture (like 
> sysdeps/i386/fpu/math-tests.h), which should have a comment pointing to a 
> bug report in GCC Bugzilla about the issue.

The problem with the load and test instruction is that it is not
consistent to doing a load and a compare separately.  Our floating
point loads do not quiet a NaN.

While the IEEE standard leaves it open to quiet NaNs with loads it is
probably supposed to be the same for all loads.

The plan is to disable the use of the FPR result of a load-and-test
instruction in the compiler.  We would leave the Glibc tests failing
since this appears to be a real finding to me.

Does this sound reasonable?

-Andreas-
Joseph Myers Dec. 14, 2016, 5:21 p.m. UTC | #4
On Wed, 14 Dec 2016, Andreas Krebbel wrote:

> The problem with the load and test instruction is that it is not
> consistent to doing a load and a compare separately.  Our floating
> point loads do not quiet a NaN.
> 
> While the IEEE standard leaves it open to quiet NaNs with loads it is
> probably supposed to be the same for all loads.
> 
> The plan is to disable the use of the FPR result of a load-and-test
> instruction in the compiler.  We would leave the Glibc tests failing
> since this appears to be a real finding to me.
> 
> Does this sound reasonable?

It sounds reasonable, given the quality-of-implementation goal to allow 
signaling NaNs to be copied without quieting.  If using that FPR result is 
useful for optimization, and harmless in the absence of signaling NaNs, of 
course you could restrict the disabling to the -fsignaling-nans case (and 
arrange for pow to be built with that option in glibc for s390).
diff mbox

Patch

diff --git a/sysdeps/ieee754/dbl-64/e_pow.c b/sysdeps/ieee754/dbl-64/e_pow.c
index db6ecf7..8f9b1c0 100644
--- a/sysdeps/ieee754/dbl-64/e_pow.c
+++ b/sysdeps/ieee754/dbl-64/e_pow.c
@@ -74,8 +74,8 @@  __ieee754_pow (double x, double y)
       qx = u.i[HIGH_HALF] & 0x7fffffff;
       /* Is x a NaN?  */
       if ((((qx == 0x7ff00000) && (u.i[LOW_HALF] != 0)) || (qx > 0x7ff00000))
-	  && y != 0)
-	return x;
+	  && (y != 0 || issignaling (x)))
+	return x + x;
       if (y == 1.0)
 	return x;
       if (y == 2.0)
@@ -129,7 +129,7 @@  __ieee754_pow (double x, double y)
     {
       if (((v.i[HIGH_HALF] & 0x7fffffff) == 0x7ff00000 && v.i[LOW_HALF] != 0)
 	  || (v.i[HIGH_HALF] & 0x7fffffff) > 0x7ff00000)	/* NaN */
-	return y;
+	return y + y;
       if (fabs (y) > 1.0e20)
 	return (y > 0) ? 0 : 1.0 / 0.0;
       k = checkint (y);
@@ -143,9 +143,9 @@  __ieee754_pow (double x, double y)
   qy = v.i[HIGH_HALF] & 0x7fffffff;	/*   no sign   */
 
   if (qx >= 0x7ff00000 && (qx > 0x7ff00000 || u.i[LOW_HALF] != 0))	/* NaN */
-    return x;
+    return x + y;
   if (qy >= 0x7ff00000 && (qy > 0x7ff00000 || v.i[LOW_HALF] != 0))	/* NaN */
-    return x == 1.0 ? 1.0 : y;
+    return x == 1.0 && !issignaling (y) ? 1.0 : y + y;
 
   /* if x<0 */
   if (u.i[HIGH_HALF] < 0)
diff --git a/sysdeps/ieee754/flt-32/e_powf.c b/sysdeps/ieee754/flt-32/e_powf.c
index c72fe37..d9470f1 100644
--- a/sysdeps/ieee754/flt-32/e_powf.c
+++ b/sysdeps/ieee754/flt-32/e_powf.c
@@ -62,10 +62,10 @@  __ieee754_powf(float x, float y)
 	ix = hx&0x7fffffff;  iy = hy&0x7fffffff;
 
     /* y==zero: x**0 = 1 */
-	if(iy==0) return one;
+	if(iy==0 && !issignaling (x)) return one;
 
     /* x==+-1 */
-	if(x == 1.0) return one;
+	if(x == 1.0 && !issignaling (y)) return one;
 	if(x == -1.0 && isinf(y)) return one;
 
     /* +-NaN return x+y */
diff --git a/sysdeps/ieee754/ldbl-128/e_powl.c b/sysdeps/ieee754/ldbl-128/e_powl.c
index e6cd975..a344840 100644
--- a/sysdeps/ieee754/ldbl-128/e_powl.c
+++ b/sysdeps/ieee754/ldbl-128/e_powl.c
@@ -165,11 +165,12 @@  __ieee754_powl (_Float128 x, _Float128 y)
 
 
   /* y==zero: x**0 = 1 */
-  if ((iy | q.parts32.w1 | q.parts32.w2 | q.parts32.w3) == 0)
+  if ((iy | q.parts32.w1 | q.parts32.w2 | q.parts32.w3) == 0
+      && !issignaling (x))
     return one;
 
   /* 1.0**y = 1; -1.0**+-Inf = 1 */
-  if (x == one)
+  if (x == one && !issignaling (y))
     return one;
   if (x == -1 && iy == 0x7fff0000
       && (q.parts32.w1 | q.parts32.w2 | q.parts32.w3) == 0)
diff --git a/sysdeps/ieee754/ldbl-128ibm/e_powl.c b/sysdeps/ieee754/ldbl-128ibm/e_powl.c
index 861b44a..d6fbef6 100644
--- a/sysdeps/ieee754/ldbl-128ibm/e_powl.c
+++ b/sysdeps/ieee754/ldbl-128ibm/e_powl.c
@@ -165,11 +165,11 @@  __ieee754_powl (long double x, long double y)
   iy = hy & 0x7fffffff;
 
   /* y==zero: x**0 = 1 */
-  if ((iy | ly) == 0)
+  if ((iy | ly) == 0 && !issignaling (x))
     return one;
 
   /* 1.0**y = 1; -1.0**+-Inf = 1 */
-  if (x == one)
+  if (x == one && !issignaling (y))
     return one;
   if (x == -1.0L && ((iy - 0x7ff00000) | ly) == 0)
     return one;