public inbox for libc-hacker@sourceware.org
 help / color / mirror / Atom feed
* [PATCH] pow{,f,l} NaN handling fixes
@ 2007-03-01 15:33 Jakub Jelinek
  0 siblings, 0 replies; only message in thread
From: Jakub Jelinek @ 2007-03-01 15:33 UTC (permalink / raw)
  To: Ulrich Drepper; +Cc: Glibc hackers

Hi!

http://www.opengroup.org/onlinepubs/009695399/functions/pow.html
says that pow{,f,l} should return NaN when either x or y is NaN,
except for the few exceptions (pow (1, NAN) and pow (NAN, +-0) should
return 1.0).  But dbl-64 generic pow was checking for NAN too late
and missed a couple of cases.  While adding tests I noticed the i386/x86_64
assembly pow* raise invalid exception when they should not, so that's
fixed in the patch as well.

2007-03-01  Jakub Jelinek  <jakub@redhat.com>

	[BZ #4096]
	* sysdeps/ieee754/dbl-64/e_pow.c (__ieee754_pow): Check for NaN earlier.
	* math/libm-test.inc (pow_test): Add more tests involving NaNs.

	* sysdeps/i386/fpu/e_powf.S (__ieee754_powf): Avoid invalid exception
	for x qNaN and y either +-inf or non-integer value.
	* sysdeps/i386/fpu/e_pow.S (__ieee754_pow): Likewise.
	* sysdeps/i386/fpu/e_powl.S (__ieee754_powl): Likewise.
	* sysdeps/x86_64/fpu/e_powl.S (__ieee754_powl): Likewise.

--- libc/math/libm-test.inc.jj	2006-09-20 21:36:26.000000000 +0200
+++ libc/math/libm-test.inc	2007-03-01 15:06:05.000000000 +0100
@@ -4614,6 +4614,17 @@ pow_test (void)
 
   /* pow (x, NaN) == NaN.  */
   TEST_ff_f (pow, 3.0, nan_value, nan_value);
+  TEST_ff_f (pow, minus_zero, nan_value, nan_value);
+  TEST_ff_f (pow, plus_infty, nan_value, nan_value);
+  TEST_ff_f (pow, -3.0, nan_value, nan_value);
+  TEST_ff_f (pow, minus_infty, nan_value, nan_value);
+
+  TEST_ff_f (pow, nan_value, 3.0, nan_value);
+  TEST_ff_f (pow, nan_value, -3.0, nan_value);
+  TEST_ff_f (pow, nan_value, plus_infty, nan_value);
+  TEST_ff_f (pow, nan_value, minus_infty, nan_value);
+  TEST_ff_f (pow, nan_value, 2.5, nan_value);
+  TEST_ff_f (pow, nan_value, -2.5, nan_value);
 
   TEST_ff_f (pow, 1, plus_infty, 1);
   TEST_ff_f (pow, -1, plus_infty, 1);
--- libc/sysdeps/i386/fpu/e_pow.S.jj	2005-05-04 19:45:15.000000000 +0200
+++ libc/sysdeps/i386/fpu/e_pow.S	2007-03-01 16:10:28.000000000 +0100
@@ -1,5 +1,5 @@
 /* ix87 specific implementation of pow function.
-   Copyright (C) 1996, 1997, 1998, 1999, 2001, 2004, 2005
+   Copyright (C) 1996, 1997, 1998, 1999, 2001, 2004, 2005, 2007
    Free Software Foundation, Inc.
    This file is part of the GNU C Library.
    Contributed by Ulrich Drepper <drepper@cygnus.com>, 1996.
@@ -161,10 +161,11 @@ ENTRY(__ieee754_pow)
 2:	/* y is a real number.  */
 	fxch			// x : y
 	fldl	MO(one)		// 1.0 : x : y
-	fld	%st(1)		// x : 1.0 : x : y
-	fsub	%st(1)		// x-1 : 1.0 : x : y
-	fabs			// |x-1| : 1.0 : x : y
-	fcompl	MO(limit)	// 1.0 : x : y
+	fldl	MO(limit)	// 0.29 : 1.0 : x : y
+	fld	%st(2)		// x : 0.29 : 1.0 : x : y
+	fsub	%st(2)		// x-1 : 0.29 : 1.0 : x : y
+	fabs			// |x-1| : 0.29 : 1.0 : x : y
+	fucompp			// 1.0 : x : y
 	fnstsw
 	fxch			// x : 1.0 : y
 	sahf
@@ -197,9 +198,10 @@ ENTRY(__ieee754_pow)
 	// y == ±inf
 	.align ALIGNARG(4)
 12:	fstp	%st(0)		// pop y
-	fldl	4(%esp)		// x
-	fabs
-	fcompl	MO(one)		// < 1, == 1, or > 1
+	fldl	MO(one)		// 1
+	fldl	4(%esp)		// x : 1
+	fabs			// abs(x) : 1
+	fucompp			// < 1, == 1, or > 1
 	fnstsw
 	andb	$0x45, %ah
 	cmpb	$0x45, %ah
--- libc/sysdeps/i386/fpu/e_powl.S.jj	2005-05-04 19:45:15.000000000 +0200
+++ libc/sysdeps/i386/fpu/e_powl.S	2007-03-01 16:02:21.000000000 +0100
@@ -1,5 +1,5 @@
 /* ix87 specific implementation of pow function.
-   Copyright (C) 1996, 1997, 1998, 1999, 2001, 2004, 2005
+   Copyright (C) 1996, 1997, 1998, 1999, 2001, 2004, 2005, 2007
    Free Software Foundation, Inc.
    This file is part of the GNU C Library.
    Contributed by Ulrich Drepper <drepper@cygnus.com>, 1996.
@@ -161,10 +161,11 @@ ENTRY(__ieee754_powl)
 2:	/* y is a real number.  */
 	fxch			// x : y
 	fldl	MO(one)		// 1.0 : x : y
-	fld	%st(1)		// x : 1.0 : x : y
-	fsub	%st(1)		// x-1 : 1.0 : x : y
-	fabs			// |x-1| : 1.0 : x : y
-	fcompl	MO(limit)	// 1.0 : x : y
+	fldl	MO(limit)	// 0.29 : 1.0 : x : y
+	fld	%st(2)		// x : 0.29 : 1.0 : x : y
+	fsub	%st(2)		// x-1 : 0.29 : 1.0 : x : y
+	fabs			// |x-1| : 0.29 : 1.0 : x : y
+	fucompp			// 1.0 : x : y
 	fnstsw
 	fxch			// x : 1.0 : y
 	sahf
@@ -210,9 +211,10 @@ ENTRY(__ieee754_powl)
 	// y == ±inf
 	.align ALIGNARG(4)
 12:	fstp	%st(0)		// pop y
-	fldt	4(%esp)		// x
-	fabs
-	fcompl	MO(one)		// < 1, == 1, or > 1
+	fldl	MO(one)		// 1
+	fldt	4(%esp)		// x : 1
+	fabs			// abs(x) : 1
+	fucompp			// < 1, == 1, or > 1
 	fnstsw
 	andb	$0x45, %ah
 	cmpb	$0x45, %ah
--- libc/sysdeps/i386/fpu/e_powf.S.jj	2005-05-04 19:45:15.000000000 +0200
+++ libc/sysdeps/i386/fpu/e_powf.S	2007-03-01 16:11:59.000000000 +0100
@@ -1,5 +1,5 @@
 /* ix87 specific implementation of pow function.
-   Copyright (C) 1996, 1997, 1999, 2001, 2004, 2005
+   Copyright (C) 1996, 1997, 1999, 2001, 2004, 2005, 2007
    Free Software Foundation, Inc.
    This file is part of the GNU C Library.
    Contributed by Ulrich Drepper <drepper@cygnus.com>, 1996.
@@ -155,10 +155,11 @@ ENTRY(__ieee754_powf)
 2:	/* y is a real number.  */
 	fxch			// x : y
 	fldl	MO(one)		// 1.0 : x : y
-	fld	%st(1)		// x : 1.0 : x : y
-	fsub	%st(1)		// x-1 : 1.0 : x : y
-	fabs			// |x-1| : 1.0 : x : y
-	fcompl	MO(limit)	// 1.0 : x : y
+	fldl	MO(limit)	// 0.29 : 1.0 : x : y
+	fld	%st(2)		// x : 0.29 : 1.0 : x : y
+	fsub	%st(2)		// x-1 : 0.29 : 1.0 : x : y
+	fabs			// |x-1| : 0.29 : 1.0 : x : y
+	fucompp			// 1.0 : x : y
 	fnstsw
 	fxch			// x : 1.0 : y
 	sahf
@@ -191,9 +192,10 @@ ENTRY(__ieee754_powf)
 	// y == ±inf
 	.align ALIGNARG(4)
 12:	fstp	%st(0)		// pop y
-	flds	4(%esp)		// x
-	fabs
-	fcompl	MO(one)		// < 1, == 1, or > 1
+	fldl	MO(one)		// 1
+	flds	4(%esp)		// x : 1
+	fabs			// abs(x) : 1
+	fucompp			// < 1, == 1, or > 1
 	fnstsw
 	andb	$0x45, %ah
 	cmpb	$0x45, %ah
--- libc/sysdeps/ieee754/dbl-64/e_pow.c.jj	2004-01-23 14:16:24.000000000 +0100
+++ libc/sysdeps/ieee754/dbl-64/e_pow.c	2007-03-01 14:23:31.000000000 +0100
@@ -106,20 +106,28 @@ double __ieee754_pow(double x, double y)
     else
       return y < 0 ? 1.0/ABS(x) : 0.0;                               /* return 0 */
   }
+
+  qx = u.i[HIGH_HALF]&0x7fffffff;  /*   no sign   */
+  qy = v.i[HIGH_HALF]&0x7fffffff;  /*   no sign   */
+
+  if (qx >= 0x7ff00000 && (qx > 0x7ff00000 || u.i[LOW_HALF] != 0)) return NaNQ.x;
+  if (qy >= 0x7ff00000 && (qy > 0x7ff00000 || v.i[LOW_HALF] != 0))
+    return x == 1.0 ? 1.0 : NaNQ.x;
+
   /* if x<0 */
   if (u.i[HIGH_HALF] < 0) {
     k = checkint(y);
     if (k==0) {
-      if ((v.i[HIGH_HALF] & 0x7fffffff) == 0x7ff00000 && v.i[LOW_HALF] == 0) {
+      if (qy == 0x7ff00000) {
 	if (x == -1.0) return 1.0;
 	else if (x > -1.0) return v.i[HIGH_HALF] < 0 ? INF.x : 0.0;
 	else return v.i[HIGH_HALF] < 0 ? 0.0 : INF.x;
       }
-      else if (u.i[HIGH_HALF] == 0xfff00000 && u.i[LOW_HALF] == 0)
+      else if (qx == 0x7ff00000)
 	return y < 0 ? 0.0 : INF.x;
       return NaNQ.x;                              /* y not integer and x<0 */
     }
-    else if (u.i[HIGH_HALF] == 0xfff00000 && u.i[LOW_HALF] == 0)
+    else if (qx == 0x7ff00000)
       {
 	if (k < 0)
 	  return y < 0 ? nZERO.x : nINF.x;
@@ -129,14 +137,6 @@ double __ieee754_pow(double x, double y)
     return (k==1)?__ieee754_pow(-x,y):-__ieee754_pow(-x,y); /* if y even or odd */
   }
   /* x>0 */
-  qx = u.i[HIGH_HALF]&0x7fffffff;  /*   no sign   */
-  qy = v.i[HIGH_HALF]&0x7fffffff;  /*   no sign   */
-
-  if (qx > 0x7ff00000 || (qx == 0x7ff00000 && u.i[LOW_HALF] != 0)) return NaNQ.x;
-                                                                 /*  if 0<x<2^-0x7fe */
-  if (qy > 0x7ff00000 || (qy == 0x7ff00000 && v.i[LOW_HALF] != 0))
-    return x == 1.0 ? 1.0 : NaNQ.x;
-                                                                 /*  if y<2^-0x7fe   */
 
   if (qx == 0x7ff00000)                              /* x= 2^-0x3ff */
     {if (y == 0) return NaNQ.x;
--- libc/sysdeps/x86_64/fpu/e_powl.S.jj	2004-07-20 09:05:58.000000000 +0200
+++ libc/sysdeps/x86_64/fpu/e_powl.S	2007-03-01 15:30:42.000000000 +0100
@@ -1,5 +1,6 @@
 /* ix87 specific implementation of pow function.
-   Copyright (C) 1996, 1997, 1998, 1999, 2001, 2004 Free Software Foundation, Inc.
+   Copyright (C) 1996, 1997, 1998, 1999, 2001, 2004, 2007
+   Free Software Foundation, Inc.
    This file is part of the GNU C Library.
    Contributed by Ulrich Drepper <drepper@cygnus.com>, 1996.
 
@@ -146,10 +147,11 @@ ENTRY(__ieee754_powl)
 2:	/* y is a real number.  */
 	fxch			// x : y
 	fldl	MO(one)		// 1.0 : x : y
-	fld	%st(1)		// x : 1.0 : x : y
-	fsub	%st(1)		// x-1 : 1.0 : x : y
-	fabs			// |x-1| : 1.0 : x : y
-	fcompl	MO(limit)	// 1.0 : x : y
+	fldl	MO(limit)	// 0.29 : 1.0 : x : y
+	fld	%st(2)		// x : 0.29 : 1.0 : x : y
+	fsub	%st(2)		// x-1 : 0.29 : 1.0 : x : y
+	fabs			// |x-1| : 0.29 : 1.0 : x : y
+	fucompp			// 1.0 : x : y
 	fnstsw
 	fxch			// x : 1.0 : y
 	test	$4500,%eax
@@ -190,9 +192,10 @@ ENTRY(__ieee754_powl)
 	// y == ±inf
 	.align ALIGNARG(4)
 12:	fstp	%st(0)		// pop y
-	fldt	8(%rsp)		// x
-	fabs
-	fcompl	MO(one)		// < 1, == 1, or > 1
+	fldl	MO(one)		// 1
+	fldt	8(%rsp)		// x : 1
+	fabs			// abs(x) : 1
+	fucompp			// < 1, == 1, or > 1
 	fnstsw
 	andb	$0x45, %ah
 	cmpb	$0x45, %ah

	Jakub

^ permalink raw reply	[flat|nested] only message in thread

only message in thread, other threads:[~2007-03-01 15:33 UTC | newest]

Thread overview: (only message) (download: mbox.gz / follow: Atom feed)
-- links below jump to the message on this page --
2007-03-01 15:33 [PATCH] pow{,f,l} NaN handling fixes Jakub Jelinek

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).