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

["Tatsuro MATSUOKA" <tmmingw@xxxxxxxxx>]

Hello
Thank you for your prompt treating the matter.

Regards

Tatsuro

2008/4/25, Greg Chicares <gchicares@xxxxxxxxxxxxx>:
> 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/fidB63EAFA2C117AAFC1CF9FE9691F8DDBE4DD01A22.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
>
>
>
>  -------------------------------------------------------------------------
>  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
>

-------------------------------------------------------------------------
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