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

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

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

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 */
+.LC0:	/* 0.5 */
-	.long 0x00000000
-.LC1:	/* 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)
+	lfs	fp10,.LC0@l(r9)
 #endif
-#ifdef SHARED
+	fabs	fp2, fp1	/* Get the absolute value of x.  */
-	lfs	fp10,.LC1-.LC0(r9)
+	fsub	fp12,fp10,fp10	/* Compute 0.0.  */
-#else
+	fcmpu	cr6, fp2, fp10	/* if |x| < 0.5  */
-	lis	r9,.LC1@ha
+	fcmpu	cr3, fp1, fp12	/* x is negative? x < 0.0  */
-	lfs	fp10,.LC1@l(r9)
+	blt-	cr6,.Lretzero
-#endif
+	fadd	fp3,fp2,fp10	/* |x|+=0.5 bias to prepare to round.  */
-	fcmpu	cr6,fp1,fp12	/* if (x > 0.0)  */
+	bge	cr3,.Lconvert	/* x is positive so don't negate x.  */
-	ble-	cr6,.L4
+	fnabs	fp3,fp3		/* -(|x|+=0.5)  */ 
-	fadd	fp1,fp1,fp10	/* x+= 0.5;  */
+.Lconvert:
-.L9:
+	fctiwz	fp4,fp3		/* Convert to Integer word lround toward 0.  */
-	fctiwz	fp2,fp1		/* Convert To Integer DW lround toward 0.  */
+	stfd	fp4,8(r1)
-	stfd	fp2,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:
+.Lretzero:			/* when 0.5 > x > -0.5  */
-	fsub	fp1,fp1,fp10	/* x-= 0.5;  */
+	li	r3,0		/* return 0.  */
-	b	.L9
+	b	.Lout
 	END (__lround)
 
 weak_alias (__lround, lround)

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
+ .section .rodata.cst12,"aM",@progbits,12
-	.align	2
+ .align 3
-.LC0:	/* 0.0 */
+ .LC0:   /* 0x1.0000000000000p+52 == 2^52 */
+	.long 0x43300000
 	.long 0x00000000
-.LC1:	/* 0.5 */
+	.long 0x3f000000 /* Use this for 0.5  */
-	.long 0x3f000000
 
 	.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)
+	lfd	fp9,0(r9)
-	lfs	fp10,.LC1-.LC0(r9)
+	lfs	fp10,8(r9)
 #else
-	lis	r9,.LC0@ha
+	lis r9,.LC0@ha
-	lis	r10,.LC1@ha
+	lfd fp9,.LC0@l(r9)	/* Load 2^52 into fpr9.  */
-	lfs	fp12,.LC0@l(r9)
+	lfs fp10,.LC0@l+8(r9)	/* Load 0.5 into fpr10.  */
-	lfs	fp10,.LC1@l(r10)
 #endif
-	fcmpu	cr6,fp1,fp12	/* if (x > 0.0)  */
+	fabs	fp2,fp1		/* Get the absolute value of x.  */
-	ble-	cr6,.L4
+	fsub	fp12,fp10,fp10	/* Compute 0.0 into fpr12.  */
-	fadd	fp1,fp1,fp10	/* x+= 0.5;  */
+	fcmpu	cr6,fp2,fp10	/* if |x| < 0.5  */
-.L9:
+	fcmpu	cr4,fp2,fp9	/* if |x| >= 2^52  */
-	fctidz	fp2,fp1		/* Convert To Integer DW round toward 0.  */
+	fcmpu	cr3,fp1,fp12	/* x is negative? x < 0.0  */
-	stfd	fp2,8(r1)
+	blt-	cr6,.Lretzero	/* 0.5 > x < -0.5 so just return 0.  */
-	nop	/* Ensure the following load is in a different dispatch  */
+	bge-	cr4,.Lnobias	/* 2^52 > x < -2^52 just convert with no bias.  */
-	nop	/* group to avoid pipe stall on POWER4&5.  */
+	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
-	lwz	r4,12(r1)
+	nop
+	nop
+	lwz	r4,12(r1)	/* Load return as integer.  */
 	lwz	r3,8(r1)
+.Lout:
 	addi	r1,r1,16
 	blr
-.L4:
+.Lretzero:			/* 0.5 > x > -0.5  */
-	fsub	fp1,fp1,fp10	/* x-= 0.5;  */
+	li	r3,0		/* return 0.  */
-	b	.L9
+	li	r4,0
+	b	.Lout
+.Lnobias:
+	fmr	fp3,fp1
+	b	.Lconvert
 	END (__llround)
 
 weak_alias (__llround, llround)

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 */
+.LC0:	/* 2^52 */
-	.tc FD_00000000_0[TC],0x0000000000000000
+	.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	fp9,.LC0@toc(2)	/* Load 2^52 into fpr9.  */
-	lfd	fp10,.LC1@toc(2)
+	lfd	fp10,.LC1@toc(2)/* Load 0.5 into fpr10.  */
-	fcmpu	cr6,fp1,fp12	/* if (x > 0.0)  */
+	fabs	fp2,fp1		/* Get the absolute value of x.  */
-	ble-	cr6,.L4
+	fsub	fp12,fp10,fp10	/* Compute 0.0 into fp12.  */
-	fadd	fp1,fp1,fp10	/* x+= 0.5;  */
+	fcmpu	cr6,fp2,fp10	/* if |x| < 0.5  */
-.L9:
+	fcmpu	cr4,fp2,fp9	/* if |x| >= 2^52  */
-	fctidz	fp2,fp1		/* Convert To Integer DW llround toward 0.  */
+	fcmpu	cr3,fp1,fp12	/* x is negative? x < 0.0  */
-	stfd	fp2,-16(r1)
+	blt-	cr6,.Lretzero	/* 0.5 > x < -0.5 so just return 0.  */
-	nop	/* Insure the following load is in a different dispatch group */
+	bge-	cr4,.Lnobias	/* 2^52 > x < -2^52 just convert with no bias.  */
-	nop	/* to avoid pipe stall on POWER4&5.  */
+	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
-	ld	r3,-16(r1)
+	nop
+	nop
+	ld	r3,-16(r1)	/* Load return as integer.  */
+.Lout:
 	blr
-.L4:
+.Lretzero:			/* 0.5 > x > -0.5  */
-	fsub	fp1,fp1,fp10	/* x-= 0.5;  */
+	li	r3,0		/* return 0.  */
-	b	.L9
+	b	.Lout
+.Lnobias:
+	fmr	fp3,fp1
+	b	.Lconvert
 	END (__llround)
 
 strong_alias (__llround, __lround)

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 */
+.LC0:	/* 2^23 */
-	.tc FD_00000000_0[TC],0x0000000000000000
+	.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	fp9,.LC0@toc(2)	/* Load 2^23 into fpr9.  */
-	lfd	fp10,.LC1@toc(2)
+	lfd	fp10,.LC1@toc(2)/* Load 0.5 into fpr10.  */
-	fcmpu	cr6,fp1,fp12	/* if (x < 0.0)  */
+	fabs	fp2,fp1		/* Get the absolute value of x.  */
-	fsubs	fp3,fp1,fp10	/* x-= 0.5;  */
+	fsub	fp12,fp10,fp10	/* Compute 0.0 into fp12.  */
-	ble-	cr6,.L9
+	fcmpu	cr6,fp2,fp10	/* if |x| < 0.5  */
-	fadds	fp3,fp1,fp10	/* x+= 0.5;  */
+	fcmpu	cr4,fp2,fp9	/* if |x| >= 2^23  */
-.L9:
+	fcmpu	cr3,fp1,fp12	/* x is negative? x < 0.0  */
-	fctidz	fp2,fp3		/* Convert To Integer DW round toward 0.  */
+	blt-	cr6,.Lretzero	/* 0.5 > x < -0.5 so just return 0.  */
-	stfd	fp2,-16(r1)
+	bge-	cr4,.Lnobias	/* 2^23 > x < -2^23 just convert with no bias.  */
-	nop	/* Insure the following load is in a different dispatch group */
+	fadd	fp3,fp2,fp10	/* |x|+=0.5 bias to prepare to round.  */
-	nop	/* to avoid pipe stall on POWER4&5.  */
+	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
-	ld	r3,-16(r1)
+	nop
+	nop
+	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)