Menu
▾
▴
#370 lround and llround functions are incorrect
closed-fixed
2008-06-03
2008-05-11
No
mingw lround and llround functions (and their float and long double variants) give incorrect results for some inputs, including nextafter(0.5,-1.) (the largest number smaller than 0.5) and 4.503599627370497e+15. This is due to the algorithm used, which is to add +/- 0.5 and round towards zero. The fix is to add not +/- 0.5, but +/- nextafter(0.5,-1.). See http://gcc.gnu.org/ml/gcc-patches/2006-10/msg00917.html and follow-ups for a good explaination of the issues involved. I spotted this issue because it causes failures in the gfortran testsuite.
The patch attached does that change. Please consider it for the next mingw-runtime release.
PS: a short program showing both the failure and the fix working is:
#include <math.h>
#include <limits.h>
#include <errno.h>
#include <stdio.h>
#include <stdlib.h>
long long
my_llround (double x)
{
/* Add +/- 0.5 then then round towards zero. */
double tmp = trunc (x + (x >= 0.0 ? 0.5 : -0.5));
if (!isfinite (tmp)
|| tmp > (double)LONG_LONG_MAX
|| tmp < (double)LONG_LONG_MIN)
{
errno = ERANGE;
/* Undefined behaviour, so we could return anything. */
/* return tmp > 0.0 ? LONG_LONG_MAX : LONG_LONG_MIN; */
}
return (long long)tmp;
}
long long
my_llround_fixed (double x)
{
double c;
c = nextafter (0.5, -1.);
/* Add +/- 0.5 then then round towards zero. */
double tmp = trunc (x + (x >= 0.0 ? c : -c));
if (!isfinite (tmp)
|| tmp > (double)LONG_LONG_MAX
|| tmp < (double)LONG_LONG_MIN)
{
errno = ERANGE;
/* Undefined behaviour, so we could return anything. */
/* return tmp > 0.0 ? LONG_LONG_MAX : LONG_LONG_MIN; */
}
return (long long)tmp;
}
int main (void)
{
double a;
long long int i1, i2;
a = 4.99999999999999944488848768742172978818416595459e-1L;
printf ("%lli %lli %lli\n", llround(a), my_llround(a), my_llround_fixed(a));
a = 4.503599627370497e+15L;
printf ("%lli %lli %lli\n", llround(a), my_llround(a), my_llround_fixed(a));
a = -4.503599627370497e+15L;
printf ("%lli %lli %lli\n", llround(a), my_llround(a), my_llround_fixed(a));
return 0;
}
my_llround() is the current mingw implementation, my_llround_fixed() is the one proposed in my patch. If you run it on a system with correct native llround(), you get:
0 1 0
4503599627370497 4503599627370498 4503599627370497
-4503599627370497 -4503599627370498 -4503599627370497
where the second column (current mingw) differs from the first (correct llround) and third (patched mingw).
Patch fixing the issue.
Logged In: YES
user_id=865300
Originator: YES
I should have added there's the same issue in round(), roundf() and roundl(). The same solution also applies there.
Logged In: YES
user_id=11494
Originator: NO
Your proposed patch works but ....it doesn't extend itself to fix problems with round{fl} in corner cases where we are near [FLOAT|DOUBLE|LONG_DOUBLE]_MAX. I think something on this track that uses ceil (maybe floor is more efficient, I haven't tested):
#include <math.h>
double __round (double x)
{
double res;
if (x >= 0.0)
{
res = ceil (x);
if (res - x > 0.5)
res -= 1.0;
}
else
{
res = ceil (-x);
if (res + x > 0.5)
res -= 1.0;
res = -res;;
}
return res;
}
#define __lround(x) (long)__round(x)
#define __llround(x) (long long)__round(x)
int main (void)
{
double a;
a = 4.99999999999999944488848768742172978818416595459e-1L;
printf ("%I64i %I64i\n", llround(a), __llround(a));
a = 4.503599627370497e+15L;
printf ("%I64i %I64i\n", llround(a), __llround(a));
a = -4.503599627370497e+15L;
printf ("%I64i %I64i\n", llround(a), __llround(a));
return 0;
}
Logged In: YES
user_id=865300
Originator: YES
I don't understand what you mean by "Your proposed patch works but ....it doesn't extend itself to fix problems with round{fl} in corner cases where we are near [FLOAT|DOUBLE|LONG_DOUBLE]_MAX." Are you referring to 4.503599627370497e+15L? Because the patch I propose does fix it AFAICT. Can you post an example case where it fails?
Logged In: YES
user_id=11494
Originator: NO
coudert wrote:
"I don't understand what you mean by 'Your proposed patch works but ...
doesn't extend itself to fix problems with round{fl} in corner cases...'
It would ease maintenance if round[fl] use the same basic algorithm as l[l]round[fl] so that a fix for you would also avoid the overflow problems with round
http://lists-archives.org/mingw-users/10476-about-round-function.html.
That is, I am referring to cases of roundl (float_value_near_LONG_DOUBLE_MAX_that_overflows_to_INF_when_you_add_4.99999999999999944488848768742172978818416595459e-1L)
That patch would look something like attached round.diff (testing in progress). Oops, sorry I am not allowed to attach anyhing, so here is the patch inline:
Index: mingwex/math/llround.c
RCS file: /cvs/src/src/winsup/mingw/mingwex/math/llround.c,v
retrieving revision 1.3
diff -c -3 -p -r1.3 llround.c
*** mingwex/math/llround.c 22 Apr 2004 04:38:24 -0000 1.3
--- mingwex/math/llround.c 12 May 2008 09:58:53 -0000
***************
*** 5,19 ****
long long
llround (double x)
{
! /* Add +/- 0.5, then round towards zero. */
! double tmp = trunc (x + (x >= 0.0 ? 0.5 : -0.5));
! if (!isfinite (tmp)
! || tmp > (double)LONG_LONG_MAX
! || tmp < (double)LONG_LONG_MIN)
{
errno = ERANGE;
/* Undefined behaviour, so we could return anything. */
! /* return tmp > 0.0 ? LONG_LONG_MAX : LONG_LONG_MIN; */
}
! return (long long)tmp;
}
--- 5,31 ----
long long
llround (double x)
{
! double res;
! if (x >= 0.0)
! {
! res = ceil (x);
! if (res - x > 0.5)
! res -= 1.0;
! }
! else
! {
! res = ceil (-x);
! if (res + x > 0.5)
! res -= 1.0;
! res = -res;;
! }
! if (!isfinite (res)
! || res > (double) LONG_LONG_MAX
! || res < (double) LONG_LONG_MIN)
{
errno = ERANGE;
/* Undefined behaviour, so we could return anything. */
! /* return res > 0.0 ? LONG_LONG_MAX : LONG_LONG_MIN; */
}
! return (long long) res;
}
Index: mingwex/math/llroundf.c
===================================================================
RCS file: /cvs/src/src/winsup/mingw/mingwex/math/llroundf.c,v
retrieving revision 1.2
diff -c -3 -p -r1.2 llroundf.c
*** mingwex/math/llroundf.c 22 Apr 2004 04:38:24 -0000 1.2
--- mingwex/math/llroundf.c 12 May 2008 09:58:53 -0000
***************
*** 5,19 ****
long long
llroundf (float x)
{
! /* Add +/- 0.5, then round towards zero. */
! float tmp = truncf (x + (x >= 0.0F ? 0.5F : -0.5F));
! if (!isfinite (tmp)
! || tmp > (float)LONG_LONG_MAX
! || tmp < (float)LONG_LONG_MIN)
{
errno = ERANGE;
/* Undefined behaviour, so we could return anything. */
! /* return tmp > 0.0F ? LONG_LONG_MAX : LONG_LONG_MIN; */
}
! return (long long)tmp;
! }
--- 5,31 ----
long long
llroundf (float x)
{
! float res;
! if (x >= 0.0F)
! {
! res = ceilf (x);
! if (res - x > 0.5F)
! res -= 1.0F;
! }
! else
! {
! res = ceilf (-x);
! if (res + x > 0.5F)
! res -= 1.0F;
! res = -res;
! }
! if (!isfinite (res)
! || res > (float) LONG_LONG_MAX
! || res < (float) LONG_LONG_MIN)
{
errno = ERANGE;
/* Undefined behaviour, so we could return anything. */
! /* return res > 0.0F ? LONG_LONG_MAX : LONG_LONG_MIN; */
}
! return (long long) res;
! }
Index: mingwex/math/llroundl.c
===================================================================
RCS file: /cvs/src/src/winsup/mingw/mingwex/math/llroundl.c,v
retrieving revision 1.2
diff -c -3 -p -r1.2 llroundl.c
*** mingwex/math/llroundl.c 22 Apr 2004 04:38:24 -0000 1.2
--- mingwex/math/llroundl.c 12 May 2008 09:58:53 -0000
***************
*** 5,19 ****
long long
llroundl (long double x)
{
! /* Add +/- 0.5, then round towards zero. */
! long double tmp = truncl (x + (x >= 0.0L ? 0.5L : -0.5L));
! if (!isfinite (tmp)
! || tmp > (long double)LONG_LONG_MAX
! || tmp < (long double)LONG_LONG_MIN)
{
errno = ERANGE;
/* Undefined behaviour, so we could return anything. */
! /* return tmp > 0.0L ? LONG_LONG_MAX : LONG_LONG_MIN; */
}
! return (long long)tmp;
}
--- 5,31 ----
long long
llroundl (long double x)
{
! long double res;
! if (x >= 0.0L)
! {
! res = ceill (x);
! if (res - x > 0.5L)
! res -= 1.0L;
! }
! else
! {
! res = ceill (-x);
! if (res + x > 0.5L)
! res -= 1.0L;
! res = -res;
! }
! if (!isfinite (res)
! || res > (double) LONG_LONG_MAX
! || res < (double) LONG_LONG_MIN)
{
errno = ERANGE;
/* Undefined behaviour, so we could return anything. */
! /* return res > 0.0 ? LONG_LONG_MAX : LONG_LONG_MIN; */
}
! return (long long) res;
}
Index: mingwex/math/lround.c
===================================================================
RCS file: /cvs/src/src/winsup/mingw/mingwex/math/lround.c,v
retrieving revision 1.2
diff -c -3 -p -r1.2 lround.c
*** mingwex/math/lround.c 22 Apr 2004 04:38:24 -0000 1.2
--- mingwex/math/lround.c 12 May 2008 09:58:53 -0000
***************
*** 1,19 ****
#include <math.h>
#include <limits.h>
#include <errno.h>
long
! lround (double x)
{
! /* Add +/- 0.5 then then round towards zero. */
! double tmp = trunc (x + (x >= 0.0 ? 0.5 : -0.5));
! if (!isfinite (tmp)
! || tmp > (double)LONG_MAX
! || tmp < (double)LONG_MIN)
{
errno = ERANGE;
/* Undefined behaviour, so we could return anything. */
! /* return tmp > 0.0 ? LONG_MAX : LONG_MIN; */
}
! return (long)tmp;
}
--- 1,33 ----
#include <math.h>
#include <limits.h>
#include <errno.h>
+ #include <math.h>
long
! lround (double x)
{
! double res;
! if (x >= 0.0)
! {
! res = ceil (x);
! if (res - x > 0.5)
! res -= 1.0;
! }
! else
! {
! res = ceil (-x);
! if (res + x > 0.5)
! res -= 1.0;
! res = -res;
! }
!
! if (!isfinite (res)
! || res > (double) LONG_MAX
! || res < (double) LONG_MIN)
{
errno = ERANGE;
/* Undefined behaviour, so we could return anything. */
! /* return res > 0.0 ? LONG_MAX : LONG_MIN; */
}
! return (long) res;
}
Index: mingwex/math/lroundf.c
===================================================================
RCS file: /cvs/src/src/winsup/mingw/mingwex/math/lroundf.c,v
retrieving revision 1.2
diff -c -3 -p -r1.2 lroundf.c
*** mingwex/math/lroundf.c 22 Apr 2004 04:38:24 -0000 1.2
--- mingwex/math/lroundf.c 12 May 2008 09:58:53 -0000
***************
*** 5,19 ****
long
lroundf (float x)
{
! /* Add +/- 0.5, then round towards zero. */
! float tmp = truncf (x + (x >= 0.0F ? 0.5F : -0.5F));
! if (!isfinite (tmp)
! || tmp > (float)LONG_MAX
! || tmp < (float)LONG_MIN)
{
errno = ERANGE;
/* Undefined behaviour, so we could return anything. */
! /* return tmp > 0.0F ? LONG_MAX : LONG_MIN; */
}
! return (long)tmp;
! }
--- 5,31 ----
long
lroundf (float x)
{
! float res;
! if (x >= 0.0F)
! {
! res = ceilf (x);
! if (res - x > 0.5F)
! res -= 1.0F;
! }
! else
! {
! res = ceilf (-x);
! if (res + x > 0.5F)
! res -= 1.0F;
! res = -res;
! }
! if (!isfinite (res)
! || res > (float) LONG_MAX
! || res < (float) LONG_MIN)
{
errno = ERANGE;
/* Undefined behaviour, so we could return anything. */
! /* return res > 0.0F ? LONG_MAX : LONG_MIN; */
}
! return (long) res;
! }
Index: mingwex/math/lroundl.c
===================================================================
RCS file: /cvs/src/src/winsup/mingw/mingwex/math/lroundl.c,v
retrieving revision 1.2
diff -c -3 -p -r1.2 lroundl.c
*** mingwex/math/lroundl.c 22 Apr 2004 04:38:24 -0000 1.2
--- mingwex/math/lroundl.c 12 May 2008 09:58:53 -0000
***************
*** 3,19 ****
#include <errno.h>
long
! lroundl (long double x)
{
! /* Add +/- 0.5, then round towards zero. */
! long double tmp = truncl (x + (x >= 0.0L ? 0.5L : -0.5L));
! if (!isfinite (tmp)
! || tmp > (long double)LONG_MAX
! || tmp < (long double)LONG_MIN)
{
errno = ERANGE;
/* Undefined behaviour, so we could return anything. */
! /* return tmp > 0.0L ? LONG_MAX : LONG_MIN; */
}
! return (long)tmp;
}
--- 3,32 ----
#include <errno.h>
long
! lroundl (long double x)
{
! long double res;
! if (x >= 0.0L)
! {
! res = ceill (x);
! if (res - x > 0.5L)
! res -= 1.0;
! }
! else
! {
! res = ceill (-x);
! if (res + x > 0.5L)
! res -= 1.0L;
! res = -res;;
! }
! if (!isfinite (res)
! || res > (long double)LONG_MAX
! || res < (long double)LONG_MIN)
{
errno = ERANGE;
/* Undefined behaviour, so we could return anything. */
! /* return res > 0.0L ? LONG_MAX : LONG_MIN; */
}
! return (long)res;
}
+
Index: mingwex/math/round.c
===================================================================
RCS file: /cvs/src/src/winsup/mingw/mingwex/math/round.c,v
retrieving revision 1.2
diff -c -3 -p -r1.2 round.c
*** mingwex/math/round.c 29 Mar 2004 08:22:19 -0000 1.2
--- mingwex/math/round.c 12 May 2008 09:58:53 -0000
***************
*** 1,8 ****
#include <math.h>
double
! round (double x)
{
! /* Add +/- 0.5 then then round towards zero. */
! return trunc ( x + (x >= 0.0 ? 0.5 : -0.5));
}
--- 1,22 ----
#include <math.h>
double
! round (double x)
{
! double res;
! if (x >= 0.0)
! {
! res = ceil (x);
! if (res - x > 0.5)
! res -= 1.0;
! }
! else
! {
! res = ceil (-x);
! if (res + x > 0.5)
! res -= 1.0;
! res = -res;
! }
! return res;
}
+
Index: mingwex/math/roundf.c
===================================================================
RCS file: /cvs/src/src/winsup/mingw/mingwex/math/roundf.c,v
retrieving revision 1.2
diff -c -3 -p -r1.2 roundf.c
*** mingwex/math/roundf.c 29 Mar 2004 08:22:20 -0000 1.2
--- mingwex/math/roundf.c 12 May 2008 09:58:53 -0000
***************
*** 3,8 ****
float
roundf (float x)
{
! /* Add +/- 0.5 then then round towards zero. */
! return truncf ( x + (x >= 0.0F ? 0.5F : -0.5F));
}
--- 3,22 ----
float
roundf (float x)
{
! float res;
! if (x >= 0.0f)
! {
! res = ceilf (x);
! if (res - x > 0.5f)
! res -= 1.0f;
! }
! else
! {
! res = ceilf (-x);
! if (res + x > 0.5f)
! res -= 1.0f;
! res = -res;
! }
! return res;
}
+
Index: mingwex/math/roundl.c
===================================================================
RCS file: /cvs/src/src/winsup/mingw/mingwex/math/roundl.c,v
retrieving revision 1.2
diff -c -3 -p -r1.2 roundl.c
*** mingwex/math/roundl.c 29 Mar 2004 08:22:20 -0000 1.2
--- mingwex/math/roundl.c 12 May 2008 09:58:53 -0000
***************
*** 3,8 ****
long double
roundl (long double x)
{
! /* Add +/- 0.5 then then round towards zero. */
! return truncl ( x + (x >= 0.0L ? 0.5L : -0.5L));
}
--- 3,22 ----
long double
roundl (long double x)
{
! long double res;
! if (x >= 0.0L)
! {
! res = ceill (x);
! if (res - x > 0.5L)
! res -= 1.0L;
! }
! else
! {
! res = ceill (-x);
! if (res + x > 0.5L)
! res -= 1.0L;
! res = -res;;
! }
! return res;
}
+
Logged In: YES
user_id=823908
Originator: NO
I have a number of issues with this patch; not all are related to the subject matter, but I'll list them anyway, noting those which are not immediately pertinent as "side issues".
1) The supplied test case is not valid; to fix it, I had to `#include <inttypes.h>' and replace all of the `%lli' printf format specs with `%"PRIi64"', since `%lli' is not supported by MSVCRT printf. I also had to replace the non-standard use of `LONG_LONG_MAX' and `LONG_LONG_MIN' with their standards conforming `LLONG_MAX' and `LLONG_MIN' counterparts; contrary to the comment in MinGW's limits.h, GCC on my Ubuntu 6.06 LTS box does *not* support this non-standard extension. (The issue of non-standard LONG_LONG_MAX vs. standard LLONG_MAX is, of course, a side issue in the context of this patch).
2) After fixing your test case, as described in (1), I have compiled and run it, *both* with MinGW on a Win2K box, *and* in GNU/Linux on the Ubuntu box; in *neither* case, can I reproduce your claimed results:--
2.1) MinGW GCC-3.4.5 on Win2K...
$ gcc testcase.c
$ ./a
1 1 0
4503599627370498 4503599627370498 4503599627370498
-4503599627370498 -4503599627370498 -4503599627370498
2.2) GCC-4.0.3 on Ubuntu 6.06 LTS
$ gcc -std=c99 testcase.c -lm
$ ./a.out
0 1 0
4503599627370497 4503599627370498 4503599627370498
-4503599627370497 -4503599627370498 -4503599627370498
Note that, while your patch does produce the correct result for rounding the value of `nextafter( 0.5, 0.0 )', it does *not* appear to reproducibly round correctly, (to the value you claim to be correct), in the other two cases you cite; the only function which does that, for my running of your test case, is the system supplied `llround()' on the Ubuntu box.
3) On the Ubuntu box, I have to add the `-std=c99' compiler option, else the test case will not compile; this is as it should be. However, if I add that switch for the MinGW case, then the use of LONG_LONG_MAX and LONG_LONG_MIN become invalid, and the code will not compile. Again, this is a side issue, and it goes away when LLONG_MAX and LLONG_MIN are used instead; however it still leaves the anomaly that MinGW does not require `-std=c99', where perhaps it really should.
4) Another (minor) side issue: when the rounded value lies beyond the limits of a long long int, the MinGW implementation sets `errno = ERANGE'. In this circumstance, POSIX requires that it sets `errno = EDOM'; there is no provision for any of the l[l]round*() functions to ever set `errno = ERANGE'.
I've not explored Danny's proposed alternative implementation, and while my own evaluation would seem to refute his statement that "your patch works...", I concur with his opinion that this issue needs to be resolved in a manner which simultaneously addresses https://sourceforge.net/tracker/index.php?func=detail&aid=1962656&group_id=2435&atid=302435
Logged In: YES
user_id=865300
Originator: YES
I'll reply only to the main issue, mainly that the patch I proposed doesn't behave as expected. It indeed does not work on linux while the very same testcase works fine on Darwin, where I initially tested. I think it's an issue of extra precision, and this can be fixed by performing calculations in the widest available FP mode (here, long double):
long long
my_llround_fixed (double x)
{
long double c = nextafter (0.5L, -1.L);
/* Add +/- 0.5 then then round towards zero. */
double tmp = (double) truncl((long double) x + (x >= 0.0 ? c : -c));
if (!isfinite (tmp)
|| tmp > (double)LLONG_MAX
|| tmp < (double)LLONG_MIN)
{
errno = ERANGE;
/* Undefined behaviour, so we could return anything. */
/* return tmp > 0.0 ? LLONG_MAX : LLONG_MIN; */
}
return (long long)tmp;
}
I mention this only for completeness, as it seems Danny has proposed a better patch.
Logged In: YES
user_id=823908
Originator: NO
I know you are trying to deal with imprecision in the floating point representation, but to me your proposed solution seems mathematically suspect; the more coercions between types of differing precisions you introduce, the more susceptible the implementation becomes to loss of precision in the result.
I read the thread you pointed to, and all I see is a glib statement, with no supporting mathematical basis or proof, that this `nextafter' technique is the right way to do it; without such a proof, I find it very difficult to accept that hypothesis.
Danny's approach definitely appears to be more robust. However, it also seems rather verbose. For a more compact solution, which appears to yield the correct results for each of your test inputs, and also in the limiting case near the overflow limits, please see the reference implementations of `round()' and `llround()', which I've attached to the other ticket I referred to previously.
Logged In: YES
user_id=865300
Originator: YES
A slightly more detailed explanation is there: http://gcc.gnu.org/ml/fortran/2005-04/msg00139.html
The algorithm is used in glibc for some targets; see e.g. glibc/sysdeps/alpha/fpu/s_lround.c:
long int
__lround (double x)
{
double adj;
adj = 0x1.fffffffffffffp-2; /* nextafter (0.5, 0.0) */
adj = copysign (adj, x);
return x + adj;
}
But anyway, as maintainers, it's your choice to what version you feel the most at ease with, so I'll leave you at it. Thanks!
Logged In: YES
user_id=823908
Originator: NO
I've attached a patch, (http://sourceforge.net/tracker/download.php?group_id=2435&atid=302435&file_id=278995&aid=1962656), to http://sourceforge.net/tracker/index.php?func=detail&aid=1962656&group_id=2435&atid=302435, which should fix this. As I've noted there, it isn't clear to me, how I should choose corresponding input values, to create equivalent test cases for the float and long double parameter function variants.
Logged In: YES
user_id=823908
Originator: NO
I believe this is now fixed, by the replacement implementation described here:
http://sourceforge.net/tracker/index.php?func=detail&aid=1962656&group_id=2435&atid=302435