Re: about 'round' function on mingw

[Tatsuro MATSUOKA]

 Hello
 This is just a forwarded mail that Greg has posted to mingw ML.
 Thank you here again to Greg!!

 Until mingw complier will be fixed, I modify config.h after configure and use 
Michael's routine.
So nothing is to do for the octave on this matter.

Thank you agian for all related people on this matter.

Regards
Tatsuro 

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



--------------------------------------
GANBARE! NIPPON! Win your ticket to Olympic Games 2008.
http://pr.mail.yahoo.co.jp/ganbare-nippon/