Re: [Mingw-users] about 'round' function From: Greg Chicares - 2008-04-25 01:01:15 ```On 2008-04-25 00:31Z, Tatsuro MATSUOKA wrote: > > octave:6> help bitmax > -- Built-in Function: bitmax () > Return the largest integer that can be represented as a floating > point value. On IEEE-754 compatiable systems, `bitmax' is `2^53 - > 1'. [...] > Is this information OK? Yes, thank you. ```
 [Mingw-users] about 'round' function From: Tatsuro MATSUOKA - 2008-04-24 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/help-octave/2008-April/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 floating-point 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 io-mappers.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,Amax-1),bitshift(1,Amax-2))); > > !!!!! 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 re-configration. 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 ```
 Re: [Mingw-users] about 'round' function From: Greg Chicares - 2008-04-24 22:34:33 ```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. 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 floating-point number, but without seeing its hexadecimal representation, I can't be sure. ```
 Re: [Mingw-users] about 'round' function From: Tatsuro MATSUOKA - 2008-04-24 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 : > 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. > > > 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 floating-point 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 > _______________________________________________ > MinGW-users mailing list > MinGW-users@... > > You may change your MinGW Account Options or unsubscribe at: > https://lists.sourceforge.net/lists/listinfo/mingw-users > ```
 Re: [Mingw-users] about 'round' function From: Tatsuro MATSUOKA - 2008-04-25 00:31:30 ```Hello I've searched by myself octave:6> help bitmax -- Built-in Function: bitmax () Return the largest integer that can be represented as a floating point value. On IEEE-754 compatiable systems, `bitmax' is `2^53 - 1'. bitmax is a built-in function Additional help for built-in functions and operators is available in the on-line version of the manual. Use the command `doc ' 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 : > 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 : > > > 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. > > > > > > 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 floating-point 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 > > _______________________________________________ > > MinGW-users mailing list > > MinGW-users@... > > > > You may change your MinGW Account Options or unsubscribe at: > > https://lists.sourceforge.net/lists/listinfo/mingw-users > > > ```
 Re: [Mingw-users] about 'round' function From: Greg Chicares - 2008-04-25 01:01:15 ```On 2008-04-25 00:31Z, Tatsuro MATSUOKA wrote: > > octave:6> help bitmax > -- Built-in Function: bitmax () > Return the largest integer that can be represented as a floating > point value. On IEEE-754 compatiable systems, `bitmax' is `2^53 - > 1'. [...] > Is this information OK? Yes, thank you. ```
 Re: [Mingw-users] about 'round' function From: Greg Chicares - 2008-04-25 02:55:02 ```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 < #include /* 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 #include 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 ```
 Re: [Mingw-users] about 'round' function From: Tatsuro MATSUOKA - 2008-04-25 03:10:47 ```Hello Thank you for your prompt treating the matter. Regards Tatsuro 2008/4/25, Greg Chicares : > 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 < > #include > #include > > /* 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 > #include > > 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@... > > You may change your MinGW Account Options or unsubscribe at: > https://lists.sourceforge.net/lists/listinfo/mingw-users > ```