[BZ4997]

	* sysdeps/powerpc/powerpc32/fpu/s_lround.S (__lround): Fixed erroneous
	result when x is +/-nextafter(+/-0.5,-/+1) i.e. all 1's in the
	mantissa.
	* sysdeps/powerpc/powerpc32/power4/fpu/s_llround.S (__llround):
	Likewise.  Also account for when x is an odd number between 2^52
	and 2^53-1.
	* sysdeps/powerpc/powerpc64/fpu/s_llround.S (__llround): Likewise.
	* sysdeps/powerpc/powerpc64/fpu/s_llroundf.S (__llroundf): Likewise.
	* math/libm-test.inc (lround_test, llround_test): Added test cases to
	detect aforementioned erroneous conditions.

6 files changed, 203 insertions, 87 deletions

diff --git a/ChangeLog b/ChangeLog
index dc385f2..d60b35a 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,3 +1,17 @@
+2007-11-20  Ryan S. Arnold  <rsa@us.ibm.com>
+
+	[BZ4997]
+	* sysdeps/powerpc/powerpc32/fpu/s_lround.S (__lround): Fixed erroneous
+	result when x is +/-nextafter(+/-0.5,-/+1) i.e. all 1's in the
+	mantissa.
+	* sysdeps/powerpc/powerpc32/power4/fpu/s_llround.S (__llround):
+	Likewise.  Also account for when x is an odd number between 2^52
+	and 2^53-1.
+	* sysdeps/powerpc/powerpc64/fpu/s_llround.S (__llround): Likewise.
+	* sysdeps/powerpc/powerpc64/fpu/s_llroundf.S (__llroundf): Likewise.
+	* math/libm-test.inc (lround_test, llround_test): Added test cases to
+	detect aforementioned erroneous conditions.
+
 2008-04-11  Jakub Jelinek  <jakub@redhat.com>
 
 	* configure.in: Check for -fno-section-anchors in addition to
diff --git a/math/libm-test.inc b/math/libm-test.inc
index a33a182..b8a73ae 100644
--- a/math/libm-test.inc
+++ b/math/libm-test.inc
@@ -4300,6 +4300,17 @@ lround_test (void)
 # endif
   TEST_f_l (lround, 2097152.5, 2097153);
   TEST_f_l (lround, -2097152.5, -2097153);
+  /* nextafter(0.5,-1)  */
+  TEST_f_l (lround, 0x1.fffffffffffffp-2, 0);
+  /* nextafter(-0.5,1)  */
+  TEST_f_l (lround, -0x1.fffffffffffffp-2, 0);
+#else
+  /* nextafter(0.5,-1)  */
+  TEST_f_l (lround, 0x1.fffffp-2, 0);
+  /* nextafter(-0.5,1)  */
+  TEST_f_l (lround, -0x1.fffffp-2, 0);
+  TEST_f_l (lround, 0x1.fffffep+23, 16777215);
+  TEST_f_l (lround, -0x1.fffffep+23, -16777215);
 #endif
   END (lround);
 }
@@ -4359,8 +4370,40 @@ llround_test (void)
   TEST_f_L (llround, 4294967295.5, 4294967296LL);
   /* 0x200000000 */
   TEST_f_L (llround, 8589934591.5, 8589934592LL);
+
+  /* nextafter(0.5,-1)  */
+  TEST_f_L (llround, 0x1.fffffffffffffp-2, 0);
+  /* nextafter(-0.5,1)  */
+  TEST_f_L (llround, -0x1.fffffffffffffp-2, 0);
+  /* On PowerPC an exponent of '52' is the largest incrementally
+   * representable sequence of whole-numbers in the 'double' range.  We test
+   * lround to make sure that a guard bit set during the lround operation
+   * hasn't forced an erroneous shift giving us an incorrect result.  The odd
+   * numbers between +-(2^52+1 and 2^53-1) are affected since they have the
+   * rightmost bit set.  */
+  /* +-(2^52+1)  */
+  TEST_f_L (llround, 0x1.0000000000001p+52,4503599627370497LL);
+  TEST_f_L (llround, -0x1.0000000000001p+52,-4503599627370497LL);
+  /* +-(2^53-1): Input is the last (positive and negative) incrementally
+   * representable whole-number in the 'double' range that might round
+   * erroneously.  */
+  TEST_f_L (llround, 0x1.fffffffffffffp+52, 9007199254740991LL);
+  TEST_f_L (llround, -0x1.fffffffffffffp+52, -9007199254740991LL);
+#else
+  /* nextafter(0.5,-1)  */
+  TEST_f_L (llround, 0x1.fffffep-2, 0);
+  /* nextafter(-0.5,1)  */
+  TEST_f_L (llround, -0x1.fffffep-2, 0);
+  /* As above, on PowerPC an exponent of '23' is the largest incrementally
+   * representable sequence of whole-numbers in the 'float' range.
+   * Likewise, numbers between +-(2^23+1 and 2^24-1) are affected.  */
+  TEST_f_L (llround, 0x1.000002p+23,8388609);
+  TEST_f_L (llround, -0x1.000002p+23,-8388609);
+  TEST_f_L (llround, 0x1.fffffep+23, 16777215);
+  TEST_f_L (llround, -0x1.fffffep+23, -16777215);
 #endif
 
+
 #ifdef TEST_LDOUBLE
   /* The input can only be represented in long double.  */
   TEST_f_L (llround, 4503599627370495.5L, 4503599627370496LL);
diff --git a/sysdeps/powerpc/powerpc32/fpu/s_lround.S b/sysdeps/powerpc/powerpc32/fpu/s_lround.S
index 9c534ec..ebacccc 100644
--- a/sysdeps/powerpc/powerpc32/fpu/s_lround.S
+++ b/sysdeps/powerpc/powerpc32/fpu/s_lround.S
@@ -1,5 +1,5 @@
 /* lround function.  PowerPC32 version.
-   Copyright (C) 2004, 2006 Free Software Foundation, Inc.
+   Copyright (C) 2004, 2006, 2007 Free Software Foundation, Inc.
    This file is part of the GNU C Library.
 
    The GNU C Library is free software; you can redistribute it and/or
@@ -20,13 +20,10 @@
 #include <sysdep.h>
 #include <math_ldbl_opt.h>
 
-	.section	.rodata.cst8,"aM",@progbits,8
+	.section	.rodata.cst4,"aM",@progbits,4
 	.align	2
-.LC0:	/* 0.0 */
-	.long 0x00000000
-.LC1:	/* 0.5 */
+.LC0:	/* 0.5 */
 	.long 0x3f000000
-
 	.section	".text"
 	
 /* long [r3] lround (float x [fp1])
@@ -37,7 +34,10 @@
    tie, choose the one that is even (least significant bit o).". 
    So we can't use the PowerPC "round to Nearest" mode. Instead we set
    "round toward Zero" mode and round by adding +-0.5 before rounding
-   to the integer value.  */
+   to the integer value.  It is necessary to detect when x is
+   (+-)0x1.fffffffffffffp-2 because adding +-0.5 in this case will
+   cause an erroneous shift, carry and round.  We simply return 0 if
+   0.5 > x > -0.5.  */
 
 ENTRY (__lround)
 	stwu    r1,-16(r1)
@@ -49,40 +49,40 @@ ENTRY (__lround)
 	bcl	20,31,1f
 1:	mflr	r9
 	addis	r9,r9,.LC0-1b@ha
-	addi	r9,r9,.LC0-1b@l
+	lfs	fp10,.LC0-1b@l(r9)
 # else
 	bl	_GLOBAL_OFFSET_TABLE_@local-4
 	mflr	r10
 	lwz	r9,.LC0@got(10)
+	lfs	fp10,0(r9)
 # endif
 	mtlr	r11
 	cfi_same_value (lr)
-	lfs	fp12,0(r9)
 #else
 	lis	r9,.LC0@ha
-	lfs	fp12,.LC0@l(r9)
-#endif
-#ifdef SHARED
-	lfs	fp10,.LC1-.LC0(r9)
-#else
-	lis	r9,.LC1@ha
-	lfs	fp10,.LC1@l(r9)
+	lfs	fp10,.LC0@l(r9)
 #endif
-	fcmpu	cr6,fp1,fp12	/* if (x > 0.0)  */
-	ble-	cr6,.L4
-	fadd	fp1,fp1,fp10	/* x+= 0.5;  */
-.L9:
-	fctiwz	fp2,fp1		/* Convert To Integer DW lround toward 0.  */
-	stfd	fp2,8(r1)
+	fabs	fp2, fp1	/* Get the absolute value of x.  */
+	fsub	fp12,fp10,fp10	/* Compute 0.0.  */
+	fcmpu	cr6, fp2, fp10	/* if |x| < 0.5  */
+	fcmpu	cr3, fp1, fp12	/* x is negative? x < 0.0  */
+	blt-	cr6,.Lretzero
+	fadd	fp3,fp2,fp10	/* |x|+=0.5 bias to prepare to round.  */
+	bge	cr3,.Lconvert	/* x is positive so don't negate x.  */
+	fnabs	fp3,fp3		/* -(|x|+=0.5)  */ 
+.Lconvert:
+	fctiwz	fp4,fp3		/* Convert to Integer word lround toward 0.  */
+	stfd	fp4,8(r1)
 	nop	/* Ensure the following load is in a different dispatch  */
 	nop	/* group to avoid pipe stall on POWER4&5.  */
 	nop
-	lwz	r3,12(r1)
+	lwz	r3,12(r1)	/* Load return as integer.  */
+.Lout:
 	addi	r1,r1,16
 	blr
-.L4:
-	fsub	fp1,fp1,fp10	/* x-= 0.5;  */
-	b	.L9
+.Lretzero:			/* when 0.5 > x > -0.5  */
+	li	r3,0		/* return 0.  */
+	b	.Lout
 	END (__lround)
 
 weak_alias (__lround, lround)
diff --git a/sysdeps/powerpc/powerpc32/power4/fpu/s_llround.S b/sysdeps/powerpc/powerpc32/power4/fpu/s_llround.S
index 952d2aa..4b1691e 100644
--- a/sysdeps/powerpc/powerpc32/power4/fpu/s_llround.S
+++ b/sysdeps/powerpc/powerpc32/power4/fpu/s_llround.S
@@ -1,5 +1,5 @@
 /* llround function.  PowerPC32 on PowerPC64 version.
-   Copyright (C) 2004, 2006 Free Software Foundation, Inc.
+   Copyright (C) 2004, 2006, 2007 Free Software Foundation, Inc.
    This file is part of the GNU C Library.
 
    The GNU C Library is free software; you can redistribute it and/or
@@ -20,15 +20,15 @@
 #include <sysdep.h>
 #include <math_ldbl_opt.h>
 
-	.section	.rodata.cst8,"aM",@progbits,8
-	.align	2
-.LC0:	/* 0.0 */
+ .section .rodata.cst12,"aM",@progbits,12
+ .align 3
+ .LC0:   /* 0x1.0000000000000p+52 == 2^52 */
+	.long 0x43300000
 	.long 0x00000000
-.LC1:	/* 0.5 */
-	.long 0x3f000000
+	.long 0x3f000000 /* Use this for 0.5  */
 
 	.section	".text"
-	
+
 /* long [r3] lround (float x [fp1])
    IEEE 1003.1 lround function.  IEEE specifies "round to the nearest 
    integer value, rounding halfway cases away from zero, regardless of
@@ -37,7 +37,15 @@
    tie, choose the one that is even (least significant bit o).". 
    So we can't use the PowerPC "round to Nearest" mode. Instead we set
    "round toward Zero" mode and round by adding +-0.5 before rounding
-   to the integer value.  */
+   to the integer value.
+
+   It is necessary to detect when x is (+-)0x1.fffffffffffffp-2
+   because adding +-0.5 in this case will cause an erroneous shift,
+   carry and round.  We simply return 0 if 0.5 > x > -0.5.  Likewise
+   if x is and odd number between +-(2^52 and 2^53-1) a shift and
+   carry will erroneously round if biased with +-0.5.  Therefore if x
+   is greater/less than +-2^52 we don't need to bias the number with
+   +-0.5.  */
 
 ENTRY (__llround)
 	stwu    r1,-16(r1)
@@ -57,30 +65,41 @@ ENTRY (__llround)
 # endif
 	mtlr	r11
 	cfi_same_value (lr)
-	lfs	fp12,0(r9)
-	lfs	fp10,.LC1-.LC0(r9)
+	lfd	fp9,0(r9)
+	lfs	fp10,8(r9)
 #else
-	lis	r9,.LC0@ha
-	lis	r10,.LC1@ha
-	lfs	fp12,.LC0@l(r9)
-	lfs	fp10,.LC1@l(r10)
+	lis r9,.LC0@ha
+	lfd fp9,.LC0@l(r9)	/* Load 2^52 into fpr9.  */
+	lfs fp10,.LC0@l+8(r9)	/* Load 0.5 into fpr10.  */
 #endif
-	fcmpu	cr6,fp1,fp12	/* if (x > 0.0)  */
-	ble-	cr6,.L4
-	fadd	fp1,fp1,fp10	/* x+= 0.5;  */
-.L9:
-	fctidz	fp2,fp1		/* Convert To Integer DW round toward 0.  */
-	stfd	fp2,8(r1)
-	nop	/* Ensure the following load is in a different dispatch  */
-	nop	/* group to avoid pipe stall on POWER4&5.  */
+	fabs	fp2,fp1		/* Get the absolute value of x.  */
+	fsub	fp12,fp10,fp10	/* Compute 0.0 into fpr12.  */
+	fcmpu	cr6,fp2,fp10	/* if |x| < 0.5  */
+	fcmpu	cr4,fp2,fp9	/* if |x| >= 2^52  */
+	fcmpu	cr3,fp1,fp12	/* x is negative? x < 0.0  */
+	blt-	cr6,.Lretzero	/* 0.5 > x < -0.5 so just return 0.  */
+	bge-	cr4,.Lnobias	/* 2^52 > x < -2^52 just convert with no bias.  */
+	fadd	fp3,fp2,fp10	/* |x|+=0.5 bias to prepare to round.  */
+	bge	cr3,.Lconvert	/* x is positive so don't negate x.  */
+	fnabs	fp3,fp3		/* -(|x|+=0.5)  */
+.Lconvert:
+	fctidz	fp4,fp3		/* Convert to Integer double word round toward 0.  */
+	stfd	fp4,8(r1)
+	nop
+	nop
 	nop
-	lwz	r4,12(r1)
+	lwz	r4,12(r1)	/* Load return as integer.  */
 	lwz	r3,8(r1)
+.Lout:
 	addi	r1,r1,16
 	blr
-.L4:
-	fsub	fp1,fp1,fp10	/* x-= 0.5;  */
-	b	.L9
+.Lretzero:			/* 0.5 > x > -0.5  */
+	li	r3,0		/* return 0.  */
+	li	r4,0
+	b	.Lout
+.Lnobias:
+	fmr	fp3,fp1
+	b	.Lconvert
 	END (__llround)
 
 weak_alias (__llround, llround)
diff --git a/sysdeps/powerpc/powerpc64/fpu/s_llround.S b/sysdeps/powerpc/powerpc64/fpu/s_llround.S
index d023b8f..4134847 100644
--- a/sysdeps/powerpc/powerpc64/fpu/s_llround.S
+++ b/sysdeps/powerpc/powerpc64/fpu/s_llround.S
@@ -1,5 +1,5 @@
 /* llround function.  PowerPC64 version.
-   Copyright (C) 2004, 2006 Free Software Foundation, Inc.
+   Copyright (C) 2004, 2006, 2007 Free Software Foundation, Inc.
    This file is part of the GNU C Library.
 
    The GNU C Library is free software; you can redistribute it and/or
@@ -21,13 +21,13 @@
 #include <math_ldbl_opt.h>
 
 	.section	".toc","aw"
-.LC0:	/* -0.0 */
-	.tc FD_00000000_0[TC],0x0000000000000000
+.LC0:	/* 2^52 */
+	.tc FD_43300000_0[TC],0x4330000000000000
 .LC1:	/* 0.5 */
 	.tc FD_3fe00000_0[TC],0x3fe0000000000000
 	.section	".text"
 	
-/* long long [r3] llround (float x [fp1])
+/* long long [r3] llround (double x [fp1])
    IEEE 1003.1 llround function.  IEEE specifies "round to the nearest 
    integer value, rounding halfway cases away from zero, regardless of
    the current rounding mode."  However PowerPC Architecture defines
@@ -35,26 +35,45 @@
    tie, choose the one that is even (least significant bit o).". 
    So we can't use the PowerPC "round to Nearest" mode. Instead we set
    "round toward Zero" mode and round by adding +-0.5 before rounding
-   to the integer value.  */
+   to the integer value.
+
+   It is necessary to detect when x is (+-)0x1.fffffffffffffp-2
+   because adding +-0.5 in this case will cause an erroneous shift,
+   carry and round.  We simply return 0 if 0.5 > x > -0.5.  Likewise
+   if x is and odd number between +-(2^52 and 2^53-1) a shift and
+   carry will erroneously round if biased with +-0.5.  Therefore if x
+   is greater/less than +-2^52 we don't need to bias the number with
+   +-0.5.  */
 
 ENTRY (__llround)
 	CALL_MCOUNT 0
-	lfd	fp12,.LC0@toc(2)
-	lfd	fp10,.LC1@toc(2)
-	fcmpu	cr6,fp1,fp12	/* if (x > 0.0)  */
-	ble-	cr6,.L4
-	fadd	fp1,fp1,fp10	/* x+= 0.5;  */
-.L9:
-	fctidz	fp2,fp1		/* Convert To Integer DW llround toward 0.  */
-	stfd	fp2,-16(r1)
-	nop	/* Insure the following load is in a different dispatch group */
-	nop	/* to avoid pipe stall on POWER4&5.  */
+	lfd	fp9,.LC0@toc(2)	/* Load 2^52 into fpr9.  */
+	lfd	fp10,.LC1@toc(2)/* Load 0.5 into fpr10.  */
+	fabs	fp2,fp1		/* Get the absolute value of x.  */
+	fsub	fp12,fp10,fp10	/* Compute 0.0 into fp12.  */
+	fcmpu	cr6,fp2,fp10	/* if |x| < 0.5  */
+	fcmpu	cr4,fp2,fp9	/* if |x| >= 2^52  */
+	fcmpu	cr3,fp1,fp12	/* x is negative? x < 0.0  */
+	blt-	cr6,.Lretzero	/* 0.5 > x < -0.5 so just return 0.  */
+	bge-	cr4,.Lnobias	/* 2^52 > x < -2^52 just convert with no bias.  */
+	fadd	fp3,fp2,fp10	/* |x|+=0.5 bias to prepare to round.  */
+	bge	cr3,.Lconvert	/* x is positive so don't negate x.  */
+	fnabs	fp3,fp3		/* -(|x|+=0.5)  */
+.Lconvert:
+	fctidz	fp4,fp3		/* Convert to Integer double word round toward 0.  */
+	stfd	fp4,-16(r1)
+	nop
+	nop
 	nop
-	ld	r3,-16(r1)
+	ld	r3,-16(r1)	/* Load return as integer.  */
+.Lout:
 	blr
-.L4:
-	fsub	fp1,fp1,fp10	/* x-= 0.5;  */
-	b	.L9
+.Lretzero:			/* 0.5 > x > -0.5  */
+	li	r3,0		/* return 0.  */
+	b	.Lout
+.Lnobias:
+	fmr	fp3,fp1
+	b	.Lconvert
 	END (__llround)
 
 strong_alias (__llround, __lround)
diff --git a/sysdeps/powerpc/powerpc64/fpu/s_llroundf.S b/sysdeps/powerpc/powerpc64/fpu/s_llroundf.S
index bbbd054..a211879 100644
--- a/sysdeps/powerpc/powerpc64/fpu/s_llroundf.S
+++ b/sysdeps/powerpc/powerpc64/fpu/s_llroundf.S
@@ -1,5 +1,5 @@
 /* llroundf function.  PowerPC64 version.
-   Copyright (C) 2004, 2006 Free Software Foundation, Inc.
+   Copyright (C) 2004, 2006, 2007 Free Software Foundation, Inc.
    This file is part of the GNU C Library.
 
    The GNU C Library is free software; you can redistribute it and/or
@@ -20,8 +20,8 @@
 #include <sysdep.h>
 
 	.section	".toc","aw"
-.LC0:	/* -0.0 */
-	.tc FD_00000000_0[TC],0x0000000000000000
+.LC0:	/* 2^23 */
+	.tc FD_41600000_0[TC],0x4160000000000000
 .LC1:	/* 0.5 */
 	.tc FD_3fe00000_0[TC],0x3fe0000000000000
 	.section	".text"
@@ -34,24 +34,45 @@
    tie, choose the one that is even (least significant bit o).". 
    So we can't use the PowerPC "round to Nearest" mode. Instead we set
    "round toward Zero" mode and round by adding +-0.5 before rounding
-   to the integer value.  */
+   to the integer value.
+
+   It is necessary to detect when x is (+-)0x1.fffffffffffffp-2
+   because adding +-0.5 in this case will cause an erroneous shift,
+   carry and round.  We simply return 0 if 0.5 > x > -0.5.  Likewise
+   if x is and odd number between +-(2^23 and 2^24-1) a shift and
+   carry will erroneously round if biased with +-0.5.  Therefore if x
+   is greater/less than +-2^23 we don't need to bias the number with
+   +-0.5.  */
 
 ENTRY (__llroundf)
 	CALL_MCOUNT 0
-	lfd	fp12,.LC0@toc(2)
-	lfd	fp10,.LC1@toc(2)
-	fcmpu	cr6,fp1,fp12	/* if (x < 0.0)  */
-	fsubs	fp3,fp1,fp10	/* x-= 0.5;  */
-	ble-	cr6,.L9
-	fadds	fp3,fp1,fp10	/* x+= 0.5;  */
-.L9:
-	fctidz	fp2,fp3		/* Convert To Integer DW round toward 0.  */
-	stfd	fp2,-16(r1)
-	nop	/* Insure the following load is in a different dispatch group */
-	nop	/* to avoid pipe stall on POWER4&5.  */
+	lfd	fp9,.LC0@toc(2)	/* Load 2^23 into fpr9.  */
+	lfd	fp10,.LC1@toc(2)/* Load 0.5 into fpr10.  */
+	fabs	fp2,fp1		/* Get the absolute value of x.  */
+	fsub	fp12,fp10,fp10	/* Compute 0.0 into fp12.  */
+	fcmpu	cr6,fp2,fp10	/* if |x| < 0.5  */
+	fcmpu	cr4,fp2,fp9	/* if |x| >= 2^23  */
+	fcmpu	cr3,fp1,fp12	/* x is negative? x < 0.0  */
+	blt-	cr6,.Lretzero	/* 0.5 > x < -0.5 so just return 0.  */
+	bge-	cr4,.Lnobias	/* 2^23 > x < -2^23 just convert with no bias.  */
+	fadd	fp3,fp2,fp10	/* |x|+=0.5 bias to prepare to round.  */
+	bge	cr3,.Lconvert	/* x is positive so don't negate x.  */
+	fnabs	fp3,fp3		/* -(|x|+=0.5)  */
+.Lconvert:
+	fctidz	fp4,fp3		/* Convert to Integer double word round toward 0.  */
+	stfd	fp4,-16(r1)
+	nop
+	nop
 	nop
-	ld	r3,-16(r1)
+	ld	r3,-16(r1)	/* Load return as integer.  */
+.Lout:
 	blr
+.Lretzero:			/* 0.5 > x > -0.5  */
+	li	r3,0		/* return 0.  */
+	b	.Lout
+.Lnobias:
+	fmr	fp3,fp1
+	b	.Lconvert
 	END (__llroundf)
 
 strong_alias (__llroundf, __lroundf)