This is the mail archive of the
libc-alpha@sourceware.org
mailing list for the glibc project.
[patch] Fix lround() loses precision.
- From: ppluzhnikov at google dot com (Paul Pluzhnikov)
- To: libc-alpha at sourceware dot org
- Cc: ppluzhnikov at google dot com, iant at google dot com
- Date: Thu, 18 Aug 2011 12:19:26 -0700 (PDT)
- Subject: [patch] Fix lround() loses precision.
Greetings,
When using glibc-2.11, we've noticed that lround() loses precision in
64-bit mode:
//--- cut ---
#include <math.h>
#include <stdio.h>
void
a ()
{
double d = -3.65309740835E17;
printf ("%lld\n", llround (d));
}
void
b ()
{
printf ("%lld\n", llround (-3.65309740835E17));
}
int
main (int argc, char **argv)
{
a ();
b ();
return 0;
}
//--- cut ---
./a.out
-365309715065196224
-365309740835000000
(In b(), GCC optimizes lround call away and computes correct output itself.)
Analysis by Ian Lance Taylor:
The bug is this line in sysdeps/ieee754/dbl-64/s_lround.c in glibc:
result = ((long int) i0 << (j0 - 20)) | (i1 << (j0 - 52));
I believe it should be
result = ((long int) i0 << (j0 - 20)) | ((long int) i1 << (j0 - 52));
The bug will occur for doubles with an exponent >= 52 and < 63.
This doesn't affect current glibc on x86_64, as it uses
sysdeps/ieee754/dbl-64/wordsize-64/s_lround.c instead, which doesn't suffer
from that bug.
However sysdeps/ieee754/dbl-64/s_lround.c is still broken and likely is
used on some architectures.
Attached patch adds a test case and a fix.
Thanks,
--
Paul Pluzhnikov
2011-08-18 Paul Pluzhnikov <ppluzhnikov@google.com>
Ian Lance Taylor <iant@google.com>
* math/libm-test.inc (lround_test): New testcase.
* sysdeps/ieee754/dbl-64/s_lround.c (__lround): Don't lose precision.
diff --git a/math/libm-test.inc b/math/libm-test.inc
index c6ed7a3..301d4a8 100644
--- a/math/libm-test.inc
+++ b/math/libm-test.inc
@@ -4386,6 +4386,7 @@ lround_test (void)
TEST_f_l (lround, 1073741824.01, 1073741824);
# if LONG_MAX > 281474976710656
TEST_f_l (lround, 281474976710656.025, 281474976710656);
+ TEST_f_l (llround, -3.65309740835E17, -365309740835000000);
# endif
TEST_f_l (lround, 2097152.5, 2097153);
TEST_f_l (lround, -2097152.5, -2097153);
diff --git a/sysdeps/ieee754/dbl-64/s_lround.c b/sysdeps/ieee754/dbl-64/s_lround.c
index 4e1302a..a849997 100644
--- a/sysdeps/ieee754/dbl-64/s_lround.c
+++ b/sysdeps/ieee754/dbl-64/s_lround.c
@@ -51,7 +51,7 @@ __lround (double x)
else if (j0 < (int32_t) (8 * sizeof (long int)) - 1)
{
if (j0 >= 52)
- result = ((long int) i0 << (j0 - 20)) | (i1 << (j0 - 52));
+ result = ((long int) i0 << (j0 - 20)) | ((long int) i1 << (j0 - 52));
else
{
u_int32_t j = i1 + (0x80000000 >> (j0 - 20));