This is the mail archive of the
libc-alpha@sourceware.org
mailing list for the glibc project.
Fix lrint, llrint (bug 2592)
- From: "Joseph S. Myers" <joseph at codesourcery dot com>
- To: libc-alpha at sourceware dot org
- Date: Sun, 18 Jun 2006 06:57:17 +0000 (UTC)
- Subject: Fix lrint, llrint (bug 2592)
The C implementations of lrint* and llrint* always return 0 for inputs
close to 0, when depending on the rounding mode they should sometimes
return 1 or -1 instead. In addition, the constants used for 0x1p23
and 0x1p52 were sometimes double instead of float or long double
instead of double, leading to inefficient arithmetic in the wider type
and in some cases to incorrect results through double rounding. The
incorrect results for inputs close to 0 were reported for SPARC as bug
2592, and also appear on MIPS64 and presumably any target using the C
implementations. This patch fixes the implementations and adds
rounding mode tests for lrint like those already present for llrint.
2006-06-17 Joseph S. Myers <joseph@codesourcery.com>
[BZ #2592]
* math/libm-test.inc (lrint_test_tonearest): New.
(lrint_test_towardzero): New.
(lrint_test_downward): New.
(lrint_test_upward): New.
(main): Run these new tests.
* sysdeps/ieee754/dbl-64/s_llrint.c (__llrint): Correct rounding
of values near to 0.
(two52): Use double not long double.
* sysdeps/ieee754/dbl-64/s_lrint.c (__lrint): Likewise.
* sysdeps/ieee754/flt-32/s_llrintf.c (__llrintf): Likewise.
(two23): Use float not double.
* sysdeps/ieee754/flt-32/s_lrintf.c (__lrintf): Likewise.
(two23): Use float not double.
* sysdeps/ieee754/ldbl-128/s_llrintl.c (__llrintl): Likewise.
* sysdeps/ieee754/ldbl-128/s_lrintl.c (__lrintl): Likewise.
* sysdeps/ieee754/ldbl-96/s_llrintl.c (__llrintl): Likewise.
* sysdeps/ieee754/ldbl-96/s_lrintl.c (__lrintl): Likewise.
Index: math/libm-test.inc
===================================================================
RCS file: /cvs/glibc/libc/math/libm-test.inc,v
retrieving revision 1.71
diff -u -p -r1.71 libm-test.inc
--- math/libm-test.inc 16 Mar 2006 23:16:56 -0000 1.71
+++ math/libm-test.inc 17 Jun 2006 16:47:18 -0000
@@ -3273,6 +3273,166 @@ lrint_test (void)
static void
+lrint_test_tonearest (void)
+{
+ int save_round_mode;
+ START (lrint_tonearest);
+
+ save_round_mode = fegetround ();
+
+ if (!fesetround (FE_TONEAREST))
+ {
+ TEST_f_l (lrint, 0.0, 0);
+ TEST_f_l (lrint, minus_zero, 0);
+ TEST_f_l (lrint, 0.2L, 0);
+ TEST_f_l (lrint, -0.2L, 0);
+ TEST_f_l (lrint, 0.5L, 0);
+ TEST_f_l (lrint, -0.5L, 0);
+ TEST_f_l (lrint, 0.8L, 1);
+ TEST_f_l (lrint, -0.8L, -1);
+
+ TEST_f_l (lrint, 1.4L, 1);
+ TEST_f_l (lrint, -1.4L, -1);
+
+ TEST_f_l (lrint, 8388600.3L, 8388600);
+ TEST_f_l (lrint, -8388600.3L, -8388600);
+
+ TEST_f_l (lrint, 1071930.0008, 1071930);
+#ifndef TEST_FLOAT
+ TEST_f_l (lrint, 1073741824.01, 1073741824);
+# if LONG_MAX > 281474976710656
+ TEST_f_l (lrint, 281474976710656.025, 281474976710656);
+# endif
+#endif
+ }
+
+ fesetround (save_round_mode);
+
+ END (lrint_tonearest);
+}
+
+
+static void
+lrint_test_towardzero (void)
+{
+ int save_round_mode;
+ START (lrint_towardzero);
+
+ save_round_mode = fegetround ();
+
+ if (!fesetround (FE_TOWARDZERO))
+ {
+ TEST_f_l (lrint, 0.0, 0);
+ TEST_f_l (lrint, minus_zero, 0);
+ TEST_f_l (lrint, 0.2L, 0);
+ TEST_f_l (lrint, -0.2L, 0);
+ TEST_f_l (lrint, 0.5L, 0);
+ TEST_f_l (lrint, -0.5L, 0);
+ TEST_f_l (lrint, 0.8L, 0);
+ TEST_f_l (lrint, -0.8L, 0);
+
+ TEST_f_l (lrint, 1.4L, 1);
+ TEST_f_l (lrint, -1.4L, -1);
+
+ TEST_f_l (lrint, 8388600.3L, 8388600);
+ TEST_f_l (lrint, -8388600.3L, -8388600);
+
+ TEST_f_l (lrint, 1071930.0008, 1071930);
+#ifndef TEST_FLOAT
+ TEST_f_l (lrint, 1073741824.01, 1073741824);
+# if LONG_MAX > 281474976710656
+ TEST_f_l (lrint, 281474976710656.025, 281474976710656);
+# endif
+#endif
+ }
+
+ fesetround (save_round_mode);
+
+ END (lrint_towardzero);
+}
+
+
+static void
+lrint_test_downward (void)
+{
+ int save_round_mode;
+ START (lrint_downward);
+
+ save_round_mode = fegetround ();
+
+ if (!fesetround (FE_DOWNWARD))
+ {
+ TEST_f_l (lrint, 0.0, 0);
+ TEST_f_l (lrint, minus_zero, 0);
+ TEST_f_l (lrint, 0.2L, 0);
+ TEST_f_l (lrint, -0.2L, -1);
+ TEST_f_l (lrint, 0.5L, 0);
+ TEST_f_l (lrint, -0.5L, -1);
+ TEST_f_l (lrint, 0.8L, 0);
+ TEST_f_l (lrint, -0.8L, -1);
+
+ TEST_f_l (lrint, 1.4L, 1);
+ TEST_f_l (lrint, -1.4L, -2);
+
+ TEST_f_l (lrint, 8388600.3L, 8388600);
+ TEST_f_l (lrint, -8388600.3L, -8388601);
+
+ TEST_f_l (lrint, 1071930.0008, 1071930);
+#ifndef TEST_FLOAT
+ TEST_f_l (lrint, 1073741824.01, 1073741824);
+# if LONG_MAX > 281474976710656
+ TEST_f_l (lrint, 281474976710656.025, 281474976710656);
+# endif
+#endif
+ }
+
+ fesetround (save_round_mode);
+
+ END (lrint_downward);
+}
+
+
+static void
+lrint_test_upward (void)
+{
+ int save_round_mode;
+ START (lrint_upward);
+
+ save_round_mode = fegetround ();
+
+ if (!fesetround (FE_UPWARD))
+ {
+ TEST_f_l (lrint, 0.0, 0);
+ TEST_f_l (lrint, minus_zero, 0);
+ TEST_f_l (lrint, 0.2L, 1);
+ TEST_f_l (lrint, -0.2L, 0);
+ TEST_f_l (lrint, 0.5L, 1);
+ TEST_f_l (lrint, -0.5L, 0);
+ TEST_f_l (lrint, 0.8L, 1);
+ TEST_f_l (lrint, -0.8L, 0);
+
+ TEST_f_l (lrint, 1.4L, 2);
+ TEST_f_l (lrint, -1.4L, -1);
+
+ TEST_f_l (lrint, 8388600.3L, 8388601);
+ TEST_f_l (lrint, -8388600.3L, -8388600);
+
+#ifndef TEST_FLOAT
+ TEST_f_l (lrint, 1071930.0008, 1071931);
+ TEST_f_l (lrint, 1073741824.01, 1073741825);
+# if LONG_MAX > 281474976710656 && defined (TEST_LDOUBLE)
+ TEST_f_l (lrint, 281474976710656.025, 28147497671065);
+# endif
+#endif
+ }
+
+ fesetround (save_round_mode);
+
+ END (lrint_upward);
+}
+
+
+static void
llrint_test (void)
{
/* XXX this test is incomplete. We need to have a way to specifiy
@@ -5937,6 +6097,10 @@ main (int argc, char **argv)
rint_test_downward ();
rint_test_upward ();
lrint_test ();
+ lrint_test_tonearest ();
+ lrint_test_towardzero ();
+ lrint_test_downward ();
+ lrint_test_upward ();
llrint_test ();
llrint_test_tonearest ();
llrint_test_towardzero ();
Index: sysdeps/ieee754/dbl-64/s_llrint.c
===================================================================
RCS file: /cvs/glibc/libc/sysdeps/ieee754/dbl-64/s_llrint.c,v
retrieving revision 1.3
diff -u -p -r1.3 s_llrint.c
--- sysdeps/ieee754/dbl-64/s_llrint.c 1 Feb 2004 19:20:22 -0000 1.3
+++ sysdeps/ieee754/dbl-64/s_llrint.c 17 Jun 2006 16:47:19 -0000
@@ -1,6 +1,6 @@
/* Round argument to nearest integral value according to current rounding
direction.
- Copyright (C) 1997, 2004 Free Software Foundation, Inc.
+ Copyright (C) 1997, 2004, 2006 Free Software Foundation, Inc.
This file is part of the GNU C Library.
Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997.
@@ -23,7 +23,7 @@
#include "math_private.h"
-static const long double two52[2] =
+static const double two52[2] =
{
4.50359962737049600000e+15, /* 0x43300000, 0x00000000 */
-4.50359962737049600000e+15, /* 0xC3300000, 0x00000000 */
@@ -48,19 +48,14 @@ __llrint (double x)
if (j0 < 20)
{
- if (j0 < -1)
- return 0;
- else
- {
- w = two52[sx] + x;
- t = w - two52[sx];
- EXTRACT_WORDS (i0, i1, t);
- j0 = ((i0 >> 20) & 0x7ff) - 0x3ff;
- i0 &= 0xfffff;
- i0 |= 0x100000;
+ w = two52[sx] + x;
+ t = w - two52[sx];
+ EXTRACT_WORDS (i0, i1, t);
+ j0 = ((i0 >> 20) & 0x7ff) - 0x3ff;
+ i0 &= 0xfffff;
+ i0 |= 0x100000;
- result = i0 >> (20 - j0);
- }
+ result = (j0 < 0 ? 0 : i0 >> (20 - j0));
}
else if (j0 < (int32_t) (8 * sizeof (long long int)) - 1)
{
Index: sysdeps/ieee754/dbl-64/s_lrint.c
===================================================================
RCS file: /cvs/glibc/libc/sysdeps/ieee754/dbl-64/s_lrint.c,v
retrieving revision 1.3
diff -u -p -r1.3 s_lrint.c
--- sysdeps/ieee754/dbl-64/s_lrint.c 1 Feb 2004 19:01:03 -0000 1.3
+++ sysdeps/ieee754/dbl-64/s_lrint.c 17 Jun 2006 16:47:19 -0000
@@ -1,6 +1,6 @@
/* Round argument to nearest integral value according to current rounding
direction.
- Copyright (C) 1997, 2004 Free Software Foundation, Inc.
+ Copyright (C) 1997, 2004, 2006 Free Software Foundation, Inc.
This file is part of the GNU C Library.
Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997.
@@ -48,19 +48,14 @@ __lrint (double x)
if (j0 < 20)
{
- if (j0 < -1)
- return 0;
- else
- {
- w = two52[sx] + x;
- t = w - two52[sx];
- EXTRACT_WORDS (i0, i1, t);
- j0 = ((i0 >> 20) & 0x7ff) - 0x3ff;
- i0 &= 0xfffff;
- i0 |= 0x100000;
+ w = two52[sx] + x;
+ t = w - two52[sx];
+ EXTRACT_WORDS (i0, i1, t);
+ j0 = ((i0 >> 20) & 0x7ff) - 0x3ff;
+ i0 &= 0xfffff;
+ i0 |= 0x100000;
- result = i0 >> (20 - j0);
- }
+ result = (j0 < 0 ? 0 : i0 >> (20 - j0));
}
else if (j0 < (int32_t) (8 * sizeof (long int)) - 1)
{
Index: sysdeps/ieee754/flt-32/s_llrintf.c
===================================================================
RCS file: /cvs/glibc/libc/sysdeps/ieee754/flt-32/s_llrintf.c,v
retrieving revision 1.2
diff -u -p -r1.2 s_llrintf.c
--- sysdeps/ieee754/flt-32/s_llrintf.c 6 Jul 2001 04:55:55 -0000 1.2
+++ sysdeps/ieee754/flt-32/s_llrintf.c 17 Jun 2006 16:47:19 -0000
@@ -1,6 +1,6 @@
/* Round argument to nearest integral value according to current rounding
direction.
- Copyright (C) 1997 Free Software Foundation, Inc.
+ Copyright (C) 1997, 2006 Free Software Foundation, Inc.
This file is part of the GNU C Library.
Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997.
@@ -23,7 +23,7 @@
#include "math_private.h"
-static const double two23[2] =
+static const float two23[2] =
{
8.3886080000e+06, /* 0x4B000000 */
-8.3886080000e+06, /* 0xCB000000 */
@@ -49,9 +49,7 @@ __llrintf (float x)
if (j0 < (int32_t) (sizeof (long long int) * 8) - 1)
{
- if (j0 < -1)
- return 0;
- else if (j0 >= 23)
+ if (j0 >= 23)
result = (long long int) i0 << (j0 - 23);
else
{
@@ -62,7 +60,7 @@ __llrintf (float x)
i0 &= 0x7fffff;
i0 |= 0x800000;
- result = i0 >> (23 - j0);
+ result = (j0 < 0 ? 0 : i0 >> (23 - j0));
}
}
else
Index: sysdeps/ieee754/flt-32/s_lrintf.c
===================================================================
RCS file: /cvs/glibc/libc/sysdeps/ieee754/flt-32/s_lrintf.c,v
retrieving revision 1.2
diff -u -p -r1.2 s_lrintf.c
--- sysdeps/ieee754/flt-32/s_lrintf.c 6 Jul 2001 04:55:55 -0000 1.2
+++ sysdeps/ieee754/flt-32/s_lrintf.c 17 Jun 2006 16:47:19 -0000
@@ -1,6 +1,6 @@
/* Round argument to nearest integral value according to current rounding
direction.
- Copyright (C) 1997 Free Software Foundation, Inc.
+ Copyright (C) 1997, 2006 Free Software Foundation, Inc.
This file is part of the GNU C Library.
Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997.
@@ -23,7 +23,7 @@
#include "math_private.h"
-static const double two23[2] =
+static const float two23[2] =
{
8.3886080000e+06, /* 0x4B000000 */
-8.3886080000e+06, /* 0xCB000000 */
@@ -49,9 +49,7 @@ __lrintf (float x)
if (j0 < (int32_t) (sizeof (long int) * 8) - 1)
{
- if (j0 < -1)
- return 0;
- else if (j0 >= 23)
+ if (j0 >= 23)
result = (long int) i0 << (j0 - 23);
else
{
@@ -62,7 +60,7 @@ __lrintf (float x)
i0 &= 0x7fffff;
i0 |= 0x800000;
- result = i0 >> (23 - j0);
+ result = (j0 < 0 ? 0 : i0 >> (23 - j0));
}
}
else
Index: sysdeps/ieee754/ldbl-128/s_llrintl.c
===================================================================
RCS file: /cvs/glibc/libc/sysdeps/ieee754/ldbl-128/s_llrintl.c,v
retrieving revision 1.3
diff -u -p -r1.3 s_llrintl.c
--- sysdeps/ieee754/ldbl-128/s_llrintl.c 1 Feb 2006 19:42:43 -0000 1.3
+++ sysdeps/ieee754/ldbl-128/s_llrintl.c 17 Jun 2006 16:47:19 -0000
@@ -48,8 +48,6 @@ __llrintl (long double x)
if (j0 < (int32_t) (8 * sizeof (long long int)) - 1)
{
- if (j0 < -1)
- return 0;
w = two112[sx] + x;
t = w - two112[sx];
GET_LDOUBLE_WORDS64 (i0, i1, t);
@@ -57,7 +55,9 @@ __llrintl (long double x)
i0 &= 0x0000ffffffffffffLL;
i0 |= 0x0001000000000000LL;
- if (j0 <= 48)
+ if (j0 < 0)
+ result = 0;
+ else if (j0 <= 48)
result = i0 >> (48 - j0);
else
result = ((long long int) i0 << (j0 - 48)) | (i1 >> (112 - j0));
Index: sysdeps/ieee754/ldbl-128/s_lrintl.c
===================================================================
RCS file: /cvs/glibc/libc/sysdeps/ieee754/ldbl-128/s_lrintl.c,v
retrieving revision 1.3
diff -u -p -r1.3 s_lrintl.c
--- sysdeps/ieee754/ldbl-128/s_lrintl.c 12 Feb 2004 20:57:40 -0000 1.3
+++ sysdeps/ieee754/ldbl-128/s_lrintl.c 17 Jun 2006 16:47:19 -0000
@@ -1,6 +1,6 @@
/* Round argument to nearest integral value according to current rounding
direction.
- Copyright (C) 1997, 1999, 2004 Free Software Foundation, Inc.
+ Copyright (C) 1997, 1999, 2004, 2006 Free Software Foundation, Inc.
This file is part of the GNU C Library.
Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997 and
Jakub Jelinek <jj@ultra.linux.cz>, 1999.
@@ -48,19 +48,14 @@ __lrintl (long double x)
if (j0 < 48)
{
- if (j0 < -1)
- return 0;
- else
- {
- w = two112[sx] + x;
- t = w - two112[sx];
- GET_LDOUBLE_WORDS64 (i0, i1, x);
- j0 = ((i0 >> 48) & 0x7fff) - 0x3fff;
- i0 &= 0x0000ffffffffffffLL;
- i0 |= 0x0001000000000000LL;
+ w = two112[sx] + x;
+ t = w - two112[sx];
+ GET_LDOUBLE_WORDS64 (i0, i1, x);
+ j0 = ((i0 >> 48) & 0x7fff) - 0x3fff;
+ i0 &= 0x0000ffffffffffffLL;
+ i0 |= 0x0001000000000000LL;
- result = i0 >> (48 - j0);
- }
+ result = (j0 < 0 ? 0 : i0 >> (48 - j0));
}
else if (j0 < (int32_t) (8 * sizeof (long int)) - 1)
{
Index: sysdeps/ieee754/ldbl-96/s_llrintl.c
===================================================================
RCS file: /cvs/glibc/libc/sysdeps/ieee754/ldbl-96/s_llrintl.c,v
retrieving revision 1.3
diff -u -p -r1.3 s_llrintl.c
--- sysdeps/ieee754/ldbl-96/s_llrintl.c 12 Feb 2004 20:57:52 -0000 1.3
+++ sysdeps/ieee754/ldbl-96/s_llrintl.c 17 Jun 2006 16:47:19 -0000
@@ -1,6 +1,6 @@
/* Round argument to nearest integral value according to current rounding
direction.
- Copyright (C) 1997, 2004 Free Software Foundation, Inc.
+ Copyright (C) 1997, 2004, 2006 Free Software Foundation, Inc.
This file is part of the GNU C Library.
Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997.
@@ -47,9 +47,7 @@ __llrintl (long double x)
if (j0 < (int32_t) (8 * sizeof (long long int)) - 1)
{
- if (j0 < -1)
- return 0;
- else if (j0 >= 63)
+ if (j0 >= 63)
result = (((long long int) i0 << 32) | i1) << (j0 - 63);
else
{
@@ -58,7 +56,9 @@ __llrintl (long double x)
GET_LDOUBLE_WORDS (se, i0, i1, t);
j0 = (se & 0x7fff) - 0x3fff;
- if (j0 <= 31)
+ if (j0 < 0)
+ result = 0;
+ else if (j0 <= 31)
result = i0 >> (31 - j0);
else
result = ((long long int) i0 << (j0 - 31)) | (i1 >> (63 - j0));
Index: sysdeps/ieee754/ldbl-96/s_lrintl.c
===================================================================
RCS file: /cvs/glibc/libc/sysdeps/ieee754/ldbl-96/s_lrintl.c,v
retrieving revision 1.3
diff -u -p -r1.3 s_lrintl.c
--- sysdeps/ieee754/ldbl-96/s_lrintl.c 1 Feb 2004 19:11:52 -0000 1.3
+++ sysdeps/ieee754/ldbl-96/s_lrintl.c 17 Jun 2006 16:47:19 -0000
@@ -1,6 +1,6 @@
/* Round argument to nearest integral value according to current rounding
direction.
- Copyright (C) 1997, 2004 Free Software Foundation, Inc.
+ Copyright (C) 1997, 2004, 2006 Free Software Foundation, Inc.
This file is part of the GNU C Library.
Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997.
@@ -47,17 +47,12 @@ __lrintl (long double x)
if (j0 < 31)
{
- if (j0 < -1)
- return 0;
- else
- {
- w = two63[sx] + x;
- t = w - two63[sx];
- GET_LDOUBLE_WORDS (se, i0, i1, t);
- j0 = (se & 0x7fff) - 0x3fff;
+ w = two63[sx] + x;
+ t = w - two63[sx];
+ GET_LDOUBLE_WORDS (se, i0, i1, t);
+ j0 = (se & 0x7fff) - 0x3fff;
- result = i0 >> (31 - j0);
- }
+ result = (j0 < 0 ? 0 : i0 >> (31 - j0));
}
else if (j0 < (int32_t) (8 * sizeof (long int)) - 1)
{
--
Joseph S. Myers
joseph@codesourcery.com