This is the mail archive of the
libc-alpha@sourceware.org
mailing list for the glibc project.
Fix rintf rounding (bug 2548)
- From: "Joseph S. Myers" <joseph at codesourcery dot com>
- To: libc-alpha at sourceware dot org
- Date: Sun, 19 Feb 2012 19:52:24 +0000 (UTC)
- Subject: Fix rintf rounding (bug 2548)
Bug 2548 reports a case where rintf incorrectly rounds some values.
This is a case of the implementation trying to be overly clever, and
getting it wrong.
Adding +-2**23 to a finite float value of the same sign and smaller
magnitude, then subtracting it, always results in the correct value
rounded to an integer according to the current rounding direction,
with the correct exceptions raised, except that the sign of a zero
result may be wrong. So the only bit manipulation needed in this case
is to fix the sign of zero. Accordingly I propose this patch which
removes the unnecessary code (which the test in the bug shows to be
incorrect) and adds testcases for both the case in the bug and a range
of similar cases (including powers of 2 plus 0.75, from 2**46 to
2**50, to make sure there are no similar bugs for "double").
Tested on x86_64 --disable-multi-arch (without --disable-multi-arch,
the problem won't reproduce on a system with SSE4.1 because that uses
an appropriate hardware instruction directly).
Carlos, I think this will be appropriate for 2.15 branch (and older
branches) as a bug fix - once we've allowed time to see if the new
tests show up bugs in any other architectures' implementations of
these functions.
2012-02-19 Joseph Myers <joseph@codesourcery.com>
[BZ #2548]
* sysdeps/ieee754/flt-32/s_rintf.c (__rintf): Do not manipulate
bits before adding and subtracting TWO23[sx].
* math/libm-test.inc (rint_test): Add more tests.
(rint_test_tonearest): Likewise.
(rint_test_towardzero): Likewise.
(rint_test_downward): Likewise.
(rint_test_upward: Likewise.
diff --git a/math/libm-test.inc b/math/libm-test.inc
index 6243e1e..c8186c8 100644
--- a/math/libm-test.inc
+++ b/math/libm-test.inc
@@ -5037,6 +5037,22 @@ rint_test (void)
TEST_f_f (rint, 262142.75, 262143.0);
TEST_f_f (rint, 524286.75, 524287.0);
TEST_f_f (rint, 524288.75, 524289.0);
+ TEST_f_f (rint, 1048576.75, 1048577.0);
+ TEST_f_f (rint, 2097152.75, 2097153.0);
+ TEST_f_f (rint, -1048576.75, -1048577.0);
+ TEST_f_f (rint, -2097152.75, -2097153.0);
+#ifndef TEST_FLOAT
+ TEST_f_f (rint, 70368744177664.75, 70368744177665.0);
+ TEST_f_f (rint, 140737488355328.75, 140737488355329.0);
+ TEST_f_f (rint, 281474976710656.75, 281474976710657.0);
+ TEST_f_f (rint, 562949953421312.75, 562949953421313.0);
+ TEST_f_f (rint, 1125899906842624.75, 1125899906842625.0);
+ TEST_f_f (rint, -70368744177664.75, -70368744177665.0);
+ TEST_f_f (rint, -140737488355328.75, -140737488355329.0);
+ TEST_f_f (rint, -281474976710656.75, -281474976710657.0);
+ TEST_f_f (rint, -562949953421312.75, -562949953421313.0);
+ TEST_f_f (rint, -1125899906842624.75, -1125899906842625.0);
+#endif
#ifdef TEST_LDOUBLE
/* The result can only be represented in long double. */
TEST_f_f (rint, 4503599627370495.5L, 4503599627370496.0L);
@@ -5137,6 +5153,22 @@ rint_test_tonearest (void)
TEST_f_f (rint, -0.1, -0.0);
TEST_f_f (rint, -0.25, -0.0);
TEST_f_f (rint, -0.625, -1.0);
+ TEST_f_f (rint, 1048576.75, 1048577.0);
+ TEST_f_f (rint, 2097152.75, 2097153.0);
+ TEST_f_f (rint, -1048576.75, -1048577.0);
+ TEST_f_f (rint, -2097152.75, -2097153.0);
+#ifndef TEST_FLOAT
+ TEST_f_f (rint, 70368744177664.75, 70368744177665.0);
+ TEST_f_f (rint, 140737488355328.75, 140737488355329.0);
+ TEST_f_f (rint, 281474976710656.75, 281474976710657.0);
+ TEST_f_f (rint, 562949953421312.75, 562949953421313.0);
+ TEST_f_f (rint, 1125899906842624.75, 1125899906842625.0);
+ TEST_f_f (rint, -70368744177664.75, -70368744177665.0);
+ TEST_f_f (rint, -140737488355328.75, -140737488355329.0);
+ TEST_f_f (rint, -281474976710656.75, -281474976710657.0);
+ TEST_f_f (rint, -562949953421312.75, -562949953421313.0);
+ TEST_f_f (rint, -1125899906842624.75, -1125899906842625.0);
+#endif
#ifdef TEST_LDOUBLE
/* The result can only be represented in long double. */
TEST_f_f (rint, 4503599627370495.5L, 4503599627370496.0L);
@@ -5207,6 +5239,22 @@ rint_test_towardzero (void)
TEST_f_f (rint, -0.1, -0.0);
TEST_f_f (rint, -0.25, -0.0);
TEST_f_f (rint, -0.625, -0.0);
+ TEST_f_f (rint, 1048576.75, 1048576.0);
+ TEST_f_f (rint, 2097152.75, 2097152.0);
+ TEST_f_f (rint, -1048576.75, -1048576.0);
+ TEST_f_f (rint, -2097152.75, -2097152.0);
+#ifndef TEST_FLOAT
+ TEST_f_f (rint, 70368744177664.75, 70368744177664.0);
+ TEST_f_f (rint, 140737488355328.75, 140737488355328.0);
+ TEST_f_f (rint, 281474976710656.75, 281474976710656.0);
+ TEST_f_f (rint, 562949953421312.75, 562949953421312.0);
+ TEST_f_f (rint, 1125899906842624.75, 1125899906842624.0);
+ TEST_f_f (rint, -70368744177664.75, -70368744177664.0);
+ TEST_f_f (rint, -140737488355328.75, -140737488355328.0);
+ TEST_f_f (rint, -281474976710656.75, -281474976710656.0);
+ TEST_f_f (rint, -562949953421312.75, -562949953421312.0);
+ TEST_f_f (rint, -1125899906842624.75, -1125899906842624.0);
+#endif
#ifdef TEST_LDOUBLE
/* The result can only be represented in long double. */
TEST_f_f (rint, 4503599627370495.5L, 4503599627370495.0L);
@@ -5277,6 +5325,22 @@ rint_test_downward (void)
TEST_f_f (rint, -0.1, -1.0);
TEST_f_f (rint, -0.25, -1.0);
TEST_f_f (rint, -0.625, -1.0);
+ TEST_f_f (rint, 1048576.75, 1048576.0);
+ TEST_f_f (rint, 2097152.75, 2097152.0);
+ TEST_f_f (rint, -1048576.75, -1048577.0);
+ TEST_f_f (rint, -2097152.75, -2097153.0);
+#ifndef TEST_FLOAT
+ TEST_f_f (rint, 70368744177664.75, 70368744177664.0);
+ TEST_f_f (rint, 140737488355328.75, 140737488355328.0);
+ TEST_f_f (rint, 281474976710656.75, 281474976710656.0);
+ TEST_f_f (rint, 562949953421312.75, 562949953421312.0);
+ TEST_f_f (rint, 1125899906842624.75, 1125899906842624.0);
+ TEST_f_f (rint, -70368744177664.75, -70368744177665.0);
+ TEST_f_f (rint, -140737488355328.75, -140737488355329.0);
+ TEST_f_f (rint, -281474976710656.75, -281474976710657.0);
+ TEST_f_f (rint, -562949953421312.75, -562949953421313.0);
+ TEST_f_f (rint, -1125899906842624.75, -1125899906842625.0);
+#endif
#ifdef TEST_LDOUBLE
/* The result can only be represented in long double. */
TEST_f_f (rint, 4503599627370495.5L, 4503599627370495.0L);
@@ -5347,6 +5411,22 @@ rint_test_upward (void)
TEST_f_f (rint, -0.1, -0.0);
TEST_f_f (rint, -0.25, -0.0);
TEST_f_f (rint, -0.625, -0.0);
+ TEST_f_f (rint, 1048576.75, 1048577.0);
+ TEST_f_f (rint, 2097152.75, 2097153.0);
+ TEST_f_f (rint, -1048576.75, -1048576.0);
+ TEST_f_f (rint, -2097152.75, -2097152.0);
+#ifndef TEST_FLOAT
+ TEST_f_f (rint, 70368744177664.75, 70368744177665.0);
+ TEST_f_f (rint, 140737488355328.75, 140737488355329.0);
+ TEST_f_f (rint, 281474976710656.75, 281474976710657.0);
+ TEST_f_f (rint, 562949953421312.75, 562949953421313.0);
+ TEST_f_f (rint, 1125899906842624.75, 1125899906842625.0);
+ TEST_f_f (rint, -70368744177664.75, -70368744177664.0);
+ TEST_f_f (rint, -140737488355328.75, -140737488355328.0);
+ TEST_f_f (rint, -281474976710656.75, -281474976710656.0);
+ TEST_f_f (rint, -562949953421312.75, -562949953421312.0);
+ TEST_f_f (rint, -1125899906842624.75, -1125899906842624.0);
+#endif
#ifdef TEST_LDOUBLE
/* The result can only be represented in long double. */
TEST_f_f (rint, 4503599627370495.5L, 4503599627370496.0L);
diff --git a/sysdeps/ieee754/flt-32/s_rintf.c b/sysdeps/ieee754/flt-32/s_rintf.c
index 9ea9b6f..9ba6b57 100644
--- a/sysdeps/ieee754/flt-32/s_rintf.c
+++ b/sysdeps/ieee754/flt-32/s_rintf.c
@@ -26,34 +26,22 @@ float
__rintf(float x)
{
int32_t i0,j0,sx;
- u_int32_t i,i1;
float w,t;
GET_FLOAT_WORD(i0,x);
sx = (i0>>31)&1;
j0 = ((i0>>23)&0xff)-0x7f;
if(j0<23) {
if(j0<0) {
- if((i0&0x7fffffff)==0) return x;
- i1 = (i0&0x07fffff);
- i0 &= 0xfff00000;
- i0 |= ((i1|-i1)>>9)&0x400000;
- SET_FLOAT_WORD(x,i0);
w = TWO23[sx]+x;
t = w-TWO23[sx];
GET_FLOAT_WORD(i0,t);
SET_FLOAT_WORD(t,(i0&0x7fffffff)|(sx<<31));
return t;
- } else {
- i = (0x007fffff)>>j0;
- if((i0&i)==0) return x; /* x is integral */
- i>>=1;
- if((i0&i)!=0) i0 = (i0&(~i))|((0x100000)>>j0);
}
} else {
if(j0==0x80) return x+x; /* inf or NaN */
else return x; /* x is integral */
}
- SET_FLOAT_WORD(x,i0);
w = TWO23[sx]+x;
return w-TWO23[sx];
}
--
Joseph S. Myers
joseph@codesourcery.com