Fix pow of negative numbers to integer exponents (bugs 369, 2678, 3866).

5 files changed, 158 insertions, 17 deletions

diff --git a/sysdeps/i386/fpu/e_pow.S b/sysdeps/i386/fpu/e_pow.S
index 1abedf6..b61a946 100644
--- a/sysdeps/i386/fpu/e_pow.S
+++ b/sysdeps/i386/fpu/e_pow.S
@@ -114,7 +114,7 @@ ENTRY(__ieee754_pow)
 	fucomp	%st(1)		// y : x
 	fnstsw
 	sahf
-	jne	2f
+	jne	3f
 
 	/* OK, we have an integer value for y.  */
 	popl	%eax
@@ -157,7 +157,12 @@ ENTRY(__ieee754_pow)
 
 	cfi_adjust_cfa_offset (8)
 	.align ALIGNARG(4)
-2:	/* y is a real number.  */
+2:	/* y is a large integer (so even).  */
+	fxch			// x : y
+	fabs			// |x| : y
+	fxch			// y : x
+	.align ALIGNARG(4)
+3:	/* y is a real number.  */
 	fxch			// x : y
 	fldl	MO(one)		// 1.0 : x : y
 	fldl	MO(limit)	// 0.29 : 1.0 : x : y
diff --git a/sysdeps/i386/fpu/e_powf.S b/sysdeps/i386/fpu/e_powf.S
index aa58ed2..529a96f 100644
--- a/sysdeps/i386/fpu/e_powf.S
+++ b/sysdeps/i386/fpu/e_powf.S
@@ -114,7 +114,7 @@ ENTRY(__ieee754_powf)
 	fucomp	%st(1)		// y : x
 	fnstsw
 	sahf
-	jne	2f
+	jne	3f
 
 	/* OK, we have an integer value for y.  */
 	popl	%edx
@@ -151,7 +151,12 @@ ENTRY(__ieee754_powf)
 
 	cfi_adjust_cfa_offset (4)
 	.align ALIGNARG(4)
-2:	/* y is a real number.  */
+2:	/* y is a large integer (so even).  */
+	fxch			// x : y
+	fabs			// |x| : y
+	fxch			// y : x
+	.align ALIGNARG(4)
+3:	/* y is a real number.  */
 	fxch			// x : y
 	fldl	MO(one)		// 1.0 : x : y
 	fldl	MO(limit)	// 0.29 : 1.0 : x : y
diff --git a/sysdeps/i386/fpu/e_powl.S b/sysdeps/i386/fpu/e_powl.S
index c0aa194c..0e7c05b 100644
--- a/sysdeps/i386/fpu/e_powl.S
+++ b/sysdeps/i386/fpu/e_powl.S
@@ -117,7 +117,7 @@ ENTRY(__ieee754_powl)
 	fucomp	%st(1)		// y : x
 	fnstsw
 	sahf
-	jne	2f
+	jne	3f
 
 	/* OK, we have an integer value for y.  */
 	popl	%eax
@@ -160,7 +160,14 @@ ENTRY(__ieee754_powl)
 
 	cfi_adjust_cfa_offset (8)
 	.align ALIGNARG(4)
-2:	/* y is a real number.  */
+2:	// y is a large integer (absolute value at least 1L<<63), but
+	// may be odd unless at least 1L<<64.  So it may be necessary
+	// to adjust the sign of a negative result afterwards.
+	fxch			// x : y
+	fabs			// |x| : y
+	fxch			// y : |x|
+	.align ALIGNARG(4)
+3:	/* y is a real number.  */
 	fxch			// x : y
 	fldl	MO(one)		// 1.0 : x : y
 	fldl	MO(limit)	// 0.29 : 1.0 : x : y
@@ -190,18 +197,51 @@ ENTRY(__ieee754_powl)
 	f2xm1			// 2^fract(y*log2(x))-1 : int(y*log2(x))
 	faddl	MO(one)		// 2^fract(y*log2(x)) : int(y*log2(x))
 	fscale			// 2^fract(y*log2(x))*2^int(y*log2(x)) : int(y*log2(x))
-	addl	$8, %esp
-	cfi_adjust_cfa_offset (-8)
 	fstp	%st(1)		// 2^fract(y*log2(x))*2^int(y*log2(x))
-	ret
+	jmp	29f
 
-	cfi_adjust_cfa_offset (8)
 28:	fstp	%st(1)		// y*log2(x)
 	fldl	MO(one)		// 1 : y*log2(x)
 	fscale			// 2^(y*log2(x)) : y*log2(x)
-	addl	$8, %esp
-	cfi_adjust_cfa_offset (-8)
 	fstp	%st(1)		// 2^(y*log2(x))
+29:	testb	$2, %dh
+	jz	292f
+	// x is negative.  If y is an odd integer, negate the result.
+	fldt	24(%esp)	// y : abs(result)
+	fld	%st		// y : y : abs(result)
+	fabs			// |y| : y : abs(result)
+	fcompl	MO(p64)		// y : abs(result)
+	fnstsw
+	sahf
+	jnc	291f
+	fldl	MO(p63)		// p63 : y : abs(result)
+	fxch			// y : p63 : abs(result)
+	fprem			// y%p63 : p63 : abs(result)
+	fstp	%st(1)		// y%p63 : abs(result)
+
+	// We must find out whether y is an odd integer.
+	fld	%st		// y : y : abs(result)
+	fistpll	(%esp)		// y : abs(result)
+	fildll	(%esp)		// int(y) : y : abs(result)
+	fucompp			// abs(result)
+	fnstsw
+	sahf
+	jne	292f
+
+	// OK, the value is an integer, but is it odd?
+	popl	%eax
+	cfi_adjust_cfa_offset (-4)
+	popl	%edx
+	cfi_adjust_cfa_offset (-4)
+	andb	$1, %al
+	jz	290f		// jump if not odd
+	// It's an odd integer.
+	fchs
+290:	ret
+	cfi_adjust_cfa_offset (8)
+291:	fstp	%st(0)		// abs(result)
+292:	addl	$8, %esp
+	cfi_adjust_cfa_offset (-8)
 	ret
 
 	// pow(x,±0) = 1
diff --git a/sysdeps/ieee754/k_standard.c b/sysdeps/ieee754/k_standard.c
index c3326d9..4e65bb1 100644
--- a/sysdeps/ieee754/k_standard.c
+++ b/sysdeps/ieee754/k_standard.c
@@ -508,6 +508,9 @@ __kernel_standard(double x, double y, int type)
 		exc.type = UNDERFLOW;
 		exc.name = type < 100 ? "pow" : (type < 200 ? "powf" : "powl");
 		exc.retval =  zero;
+		y *= 0.5;
+		if (x < zero && __rint (y) != y)
+		  exc.retval = -zero;
 		if (_LIB_VERSION == _POSIX_)
 		  __set_errno (ERANGE);
 		else if (!matherr(&exc)) {
@@ -1004,6 +1007,8 @@ long double
 __kernel_standard_l (long double x, long double y, int type)
 {
   double dx, dy;
+  struct exception exc;
+
   if (isfinite (x))
     {
       long double ax = fabsl (x);
@@ -1028,5 +1033,52 @@ __kernel_standard_l (long double x, long double y, int type)
     }
   else
     dy = y;
-  return __kernel_standard (dx, dy, type);
+
+  switch (type)
+    {
+    case 221:
+      /* powl (x, y) overflow.  */
+      exc.arg1 = dx;
+      exc.arg2 = dy;
+      exc.type = OVERFLOW;
+      exc.name = "powl";
+      if (_LIB_VERSION == _SVID_)
+	{
+	  exc.retval = HUGE;
+	  y *= 0.5;
+	  if (x < zero && __rintl (y) != y)
+	    exc.retval = -HUGE;
+	}
+      else
+	{
+	  exc.retval = HUGE_VAL;
+	  y *= 0.5;
+	  if (x < zero && __rintl (y) != y)
+	    exc.retval = -HUGE_VAL;
+	}
+      if (_LIB_VERSION == _POSIX_)
+	__set_errno (ERANGE);
+      else if (!matherr (&exc))
+	__set_errno (ERANGE);
+      return exc.retval;
+
+    case 222:
+      /* powl (x, y) underflow.  */
+      exc.arg1 = dx;
+      exc.arg2 = dy;
+      exc.type = UNDERFLOW;
+      exc.name = "powl";
+      exc.retval = zero;
+      y *= 0.5;
+      if (x < zero && __rintl (y) != y)
+	exc.retval = -zero;
+      if (_LIB_VERSION == _POSIX_)
+	__set_errno (ERANGE);
+      else if (!matherr (&exc))
+	__set_errno (ERANGE);
+      return exc.retval;
+
+    default:
+      return __kernel_standard (dx, dy, type);
+    }
 }
diff --git a/sysdeps/x86_64/fpu/e_powl.S b/sysdeps/x86_64/fpu/e_powl.S
index 562791d..0626ce4 100644
--- a/sysdeps/x86_64/fpu/e_powl.S
+++ b/sysdeps/x86_64/fpu/e_powl.S
@@ -107,7 +107,7 @@ ENTRY(__ieee754_powl)
 	fistpll	-8(%rsp)	// y : x
 	fildll	-8(%rsp)	// int(y) : y : x
 	fucomip	%st(1),%st	// y : x
-	jne	2f
+	jne	3f
 
 	/* OK, we have an integer value for y.  */
 	mov	-8(%rsp),%eax
@@ -145,7 +145,14 @@ ENTRY(__ieee754_powl)
 	ret
 
 	.align ALIGNARG(4)
-2:	/* y is a real number.  */
+2:	// y is a large integer (absolute value at least 1L<<63), but
+	// may be odd unless at least 1L<<64.  So it may be necessary
+	// to adjust the sign of a negative result afterwards.
+	fxch			// x : y
+	fabs			// |x| : y
+	fxch			// y : |x|
+	.align ALIGNARG(4)
+3:	/* y is a real number.  */
 	fxch			// x : y
 	fldl	MO(one)		// 1.0 : x : y
 	fldl	MO(limit)	// 0.29 : 1.0 : x : y
@@ -176,13 +183,45 @@ ENTRY(__ieee754_powl)
 	faddl	MO(one)		// 2^fract(y*log2(x)) : int(y*log2(x))
 	fscale			// 2^fract(y*log2(x))*2^int(y*log2(x)) : int(y*log2(x))
 	fstp	%st(1)		// 2^fract(y*log2(x))*2^int(y*log2(x))
-	ret
+	jmp	29f
 
 28:	fstp	%st(1)		// y*log2(x)
 	fldl	MO(one)		// 1 : y*log2(x)
 	fscale			// 2^(y*log2(x)) : y*log2(x)
 	fstp	%st(1)		// 2^(y*log2(x))
-	ret
+29:	testb	$2, %dh
+	jz	292f
+	// x is negative.  If y is an odd integer, negate the result.
+	fldt	24(%rsp)	// y : abs(result)
+	fldl	MO(p64)		// 1L<<64 : y : abs(result)
+	fld	%st(1)		// y : 1L<<64 : y : abs(result)
+	fabs			// |y| : 1L<<64 : y : abs(result)
+	fcomip	%st(1), %st	// 1L<<64 : y : abs(result)
+	fstp	%st(0)		// y : abs(result)
+	jnc	291f
+	fldl	MO(p63)		// p63 : y : abs(result)
+	fxch			// y : p63 : abs(result)
+	fprem			// y%p63 : p63 : abs(result)
+	fstp	%st(1)		// y%p63 : abs(result)
+
+	// We must find out whether y is an odd integer.
+	fld	%st		// y : y : abs(result)
+	fistpll	-8(%rsp)	// y : abs(result)
+	fildll	-8(%rsp)	// int(y) : y : abs(result)
+	fucomip	%st(1),%st	// y : abs(result)
+	ffreep	%st		// abs(result)
+	jne	292f
+
+	// OK, the value is an integer, but is it odd?
+	mov	-8(%rsp), %eax
+	mov	-4(%rsp), %edx
+	andb	$1, %al
+	jz	290f		// jump if not odd
+	// It's an odd integer.
+	fchs
+290:	ret
+291:	fstp	%st(0)		// abs(result)
+292:	ret
 
 	// pow(x,±0) = 1
 	.align ALIGNARG(4)