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

[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