From: Tatsuro MATSUOKA <tmmingw@gm...>  20080425 03:10:47

Hello Thank you for your prompt treating the matter. Regards Tatsuro 2008/4/25, Greg Chicares <gchicares@...>: > On 20080424 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 floatingpoint 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^531] > > 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 > _______________________________________________ > MinGWusers mailing list > MinGWusers@... > > You may change your MinGW Account Options or unsubscribe at: > https://lists.sourceforge.net/lists/listinfo/mingwusers > 