[Mingw-notify] [ mingw-Patches-1961860 ] lround and llround functions are incorrect

[SourceForge.net <no...@so...>]

Patches item #1961860, was opened at 2008-05-12 00:30
Message generated for change (Comment added) made by dannysmith
You can respond by visiting: 
https://sourceforge.net/tracker/?func=detail&atid=302435&aid=1961860&group_id=2435

Please note that this message will contain a full copy of the comment thread,
including the initial issue submission, for this request,
not just the latest update.
Category: runtime
Group: None
Status: Open
Resolution: None
Priority: 5
Private: No
Submitted By: François-Xavier Coudert (coudert)
Assigned to: Nobody/Anonymous (nobody)
Summary: lround and llround functions are incorrect

Initial Comment:
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).

----------------------------------------------------------------------

Comment By: Danny Smith (dannysmith)
Date: 2008-05-12 22:07

Message:
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;
  }
+ 


----------------------------------------------------------------------

Comment By: François-Xavier Coudert (coudert)
Date: 2008-05-12 21:14

Message:
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?

----------------------------------------------------------------------

Comment By: Danny Smith (dannysmith)
Date: 2008-05-12 20:06

Message:
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;
}

----------------------------------------------------------------------

Comment By: François-Xavier Coudert (coudert)
Date: 2008-05-12 00:36

Message:
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.

----------------------------------------------------------------------

You can respond by visiting: 
https://sourceforge.net/tracker/?func=detail&atid=302435&aid=1961860&group_id=2435