Re: [Mingw-users] about 'round' function

[Danny Smith <dannysmith@xxxxxxxxxxxx>]

> -----Original Message-----
> From: mingw-users-bounces@xxxxxxxxxxxxxxxxxxxxx 
> [mailto:mingw-users-bounces@xxxxxxxxxxxxxxxxxxxxx] On Behalf 
> Of Greg Chicares
> Sent: Friday, 25 April 2008 2:55 p.m.
> To: MinGW Users List
> Cc: matsuoka@xxxxxxxxxxxxxxxxxxx
> Subject: Re: [Mingw-users] about 'round' function
> 
> 
> On 2008-04-24 20:59Z, Tatsuro MATSUOKA wrote:
> > 
> >> | The code used internally by MinGW for "round" is the following
> >> | 
> (http://www.koders.com/c/fidB63EAFA2C117AAFC1CF9FE9691F8DDBE4D
> D01A22.aspx)
> >> |
> >> | double
> >> | round (double x)
> >> | {
> >> |   /* Add +/- 0.5 then then round towards zero.  */
> >> |   return trunc ( x + (x >= 0.0 ?  0.5 : -0.5));
> >> | }
> >> |
> >> | The problem with this implementation is when x is the 
> vicinity of bitmax: it
> >> | leads to floating-point overflow (in the mantissa) and 
> produces rounding error,
> >> | For instance if x = bitmax, round(x) => bitmax+1, while 
> it should return
> >> | bitmax.
> 
> [You later confirmed that bitmax is 2^53-1]
> 
> Here's a replacement. I'll have to do some more testing, but it
> should be okay: it's just mingwex's 'trunc.c' with a couple of
> lines changed.
> 
> If no one else wants to claim this, then I'll try to find time
> next week to put together a patch (including 'round[fl].c' and
> 'ChangeLog'). It looks like the lround() family already tests
> limits correctly.
> 
> cat >round_m53.c <<EOF
> 
> #include <math.h>
> #include <stdio.h>
> 
> /* Present libmingwex implementation */
> 
> double
> round0 (double x)
> {
>   /* Add +/- 0.5 then then round towards zero.  */
>   return trunc ( x + (x >= 0.0 ?  0.5 : -0.5));
> }
> 
> /* Proposed replacement */
> 
> /* --- 'round.c' begins --- */
> #include <fenv.h>
> #include <math.h>
> 
> double
> round1 (double _x){
>   double retval;
>   unsigned short saved_cw;
>   unsigned short tmp_cw;
>   __asm__ ("fnstcw %0;" : "=m" (saved_cw)); /* save FPU 
> control word */
>   tmp_cw = (saved_cw & ~(FE_TONEAREST | FE_DOWNWARD | 
> FE_UPWARD | FE_TOWARDZERO))
> 	    | (_x < 0.0 ? FE_DOWNWARD : FE_UPWARD);
>   __asm__ ("fldcw %0;" : : "m" (tmp_cw));
>   __asm__ ("frndint;" : "=t" (retval)  : "0" (_x)); /* round 
> towards zero */
>   __asm__ ("fldcw %0;" : : "m" (saved_cw) ); /* restore saved 
> control word */
>   return retval;
> }
> /* --- 'round.c' ends --- */
> 
> int main()
> {
>     double const M53 = 6361.0 * 69431.0 * 20394401.0;
> 
>     printf("%20.f\n", round (M53));
>     printf("%20.f\n", round0(M53));
>     printf("%20.f\n", round1(M53));
>     printf("%20.f\n", round1(9007199254740991.0));
> 
>     return 0;
> }
> EOF
> 
> /MinGW_/bin/gcc -W -Wall -pedantic -std=c99 round_m53.c
> ./a
>     9007199254740992
>     9007199254740992
>     9007199254740991
>     9007199254740991
> 


Greg, how is this testing going.  Your patch fixes this testcase in gfortran.
! Test that NINT gives right results even in corner cases
!
! PR 31202
! http://gcc.gnu.org/ml/fortran/2005-04/msg00139.html
!
! { dg-do run { xfail powerpc-ibm-aix* powerpc*-*-linux* } }
  real(kind=8) :: a
  integer(kind=8) :: i1, i2
  real :: b
  integer :: j1, j2

  a = nearest(0.5_8,-1.0_8)
  i2 = nint(nearest(0.5_8,-1.0_8))
  i1 = nint(a)
  if (i1 /= 0 .or. i2 /= 0) call abort

  a = 0.5_8
  i2 = nint(0.5_8)
  i1 = nint(a)
  if (i1 /= 1 .or. i2 /= 1) call abort

  a = nearest(0.5_8,1.0_8)
  i2 = nint(nearest(0.5_8,1.0_8))
  i1 = nint(a)
  if (i1 /= 1 .or. i2 /= 1) call abort

  b = nearest(0.5,-1.0)
  j2 = nint(nearest(0.5,-1.0))
  j1 = nint(b)
  if (j1 /= 0 .or. j2 /= 0) call abort

  b = 0.5
  j2 = nint(0.5)
  j1 = nint(b)
  if (j1 /= 1 .or. j2 /= 1) call abort

  b = nearest(0.5,1.0)
  j2 = nint(nearest(0.5,1.0))
  j1 = nint(b)
  if (j1 /= 1 .or. j2 /= 1) call abort

  a = 4503599627370497.0_8
  i1 = nint(a,kind=8)
  i2 = nint(4503599627370497.0_8,kind=8)
  if (i1 /= i2 .or. i1 /= 4503599627370497_8) call abort

  a = -4503599627370497.0_8
  i1 = nint(a,kind=8)
  i2 = nint(-4503599627370497.0_8,kind=8)
  if (i1 /= i2 .or. i1 /= -4503599627370497_8) call abort
  end


-------------------------------------------------------------------------
This SF.net email is sponsored by the 2008 JavaOne(SM) Conference 
Don't miss this year's exciting event. There's still time to save $100. 
Use priority code J8TL2D2. 
http://ad.doubleclick.net/clk;198757673;13503038;p?http://java.sun.com/javaone
_______________________________________________
MinGW-users mailing list
MinGW-users@xxxxxxxxxxxxxxxxxxxxx

You may change your MinGW Account Options or unsubscribe at:
https://lists.sourceforge.net/lists/listinfo/mingw-users