Work at SourceForge, help us to make it a better place! We have an immediate need for a Support Technician in our San Francisco or Denver office.
Close
From: Greg Chicares <gchicares@sb...>  20080424 22:34:33

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. Please say very precisely what "bitmax" means, so that I can write a minimal standalone testcase for your exact problem without trying to build your package. I'd guess you mean the largest normalized floatingpoint number, but without seeing its hexadecimal representation, I can't be sure. 
From: Tatsuro MATSUOKA <tmmingw@gm...>  20080424 20:59:19

Hello I'm a member of the octave team and am dealing with the octave building on mingw. In the octave test , there found the question about the implementation of function of 'round'. To tell the truth, I do not fully understand the problem. The Octave thread is the following http://www.cae.wisc.edu/pipermail/helpoctave/2008April/008958.html The digest is >  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. >  >  Octave used to use the above implementation for MSVC, but because of >  this problem, I switched it to use the (more robust) gnulib implementation. In the octave iomappers.cc the following is exist #if defined (HAVE_ROUND) return round (x); #else if (x >= 0) { double y = floor (x); if ((x  y) >= 0.5) y += 1.0; return y; } else { double y = ceil (x); if ((y  x) >= 0.5) y = 1.0; return y; } #endif } ********************************************************************** In octave building on mingw, it has fonction 'round' and HAVE_ROUND is set and round is used. In this cases, octave bitmap test failed. > > ***** assert(bitcmp(A,Amax),bitor(bitshift(1,Amax1),bitshift(1,Amax2))); > > !!!!! test failed > > error: assert (bitcmp (A, Amax),bitor (bitshift (1, Amax  1), bitshift (1, Amax  2))) > expected > > 6.7554e+015 > > but got > > 1.1259e+016 > > values do not match > > shared variables { > > Amax = 53 > > Bmax = 9.0072e+015 > > A = 2.2518e+015 > > } > > 6.7554e+015 So I modified config.h and undfined HAVE_ROUND and make again without reconfigration. The test was succcessful because the following code was usesd #else if (x >= 0) { double y = floor (x); if ((x  y) >= 0.5) y += 1.0; return y; } else { double y = ceil (x); if ((y  x) >= 0.5) y = 1.0; return y; } #endif ====================== Regards Tatsuro 
From: Greg Chicares <gchicares@sb...>  20080424 22:34:33

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. Please say very precisely what "bitmax" means, so that I can write a minimal standalone testcase for your exact problem without trying to build your package. I'd guess you mean the largest normalized floatingpoint number, but without seeing its hexadecimal representation, I can't be sure. 
From: Tatsuro MATSUOKA <tmmingw@gm...>  20080424 22:58:12

Hello Greg Thank you for your reply. Sorry. I'm an maintainer without enough knowldge for codes. I'll post your message to Octave ML and come back here. 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. > > > Please say very precisely what "bitmax" means, so that I can > write a minimal standalone testcase for your exact problem > without trying to build your package. I'd guess you mean the > largest normalized floatingpoint number, but without seeing > its hexadecimal representation, I can't be sure. > >  > 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 > 
From: Tatsuro MATSUOKA <tmmingw@gm...>  20080425 00:31:30

Hello I've searched by myself octave:6> help bitmax  Builtin Function: bitmax () Return the largest integer that can be represented as a floating point value. On IEEE754 compatiable systems, `bitmax' is `2^53  1'. bitmax is a builtin function Additional help for builtin functions and operators is available in the online version of the manual. Use the command `doc <topic>' to search the manual index. Help and information about Octave is also available on the WWW at http://www.octave.org and via the help@... mailing list. octave:7> bitmax ans = 9.00719925474099e+15 Is this information OK? Regards Tatsuro 2008/4/25, Tatsuro MATSUOKA <tmmingw@...>: > Hello Greg > Thank you for your reply. > Sorry. I'm an maintainer without enough knowldge for codes. > I'll post your message to Octave ML and come back here. > > 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. > > > > > > Please say very precisely what "bitmax" means, so that I can > > write a minimal standalone testcase for your exact problem > > without trying to build your package. I'd guess you mean the > > largest normalized floatingpoint number, but without seeing > > its hexadecimal representation, I can't be sure. > > > >  > > 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 > > > 
From: Greg Chicares <gchicares@sb...>  20080425 01:01:15

On 20080425 00:31Z, Tatsuro MATSUOKA wrote: > > octave:6> help bitmax >  Builtin Function: bitmax () > Return the largest integer that can be represented as a floating > point value. On IEEE754 compatiable systems, `bitmax' is `2^53  > 1'. [...] > Is this information OK? Yes, thank you. 
From: Greg Chicares <gchicares@sb...>  20080425 02:55:02

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