[Mingw-notify] [ mingw-Patches-1962656 ] Fix round(2^53)

SourceForge Headquarters 225 Broadway Suite 1600 San Diego, CA 92101 +1 (858) 454-5900

Patches item #1962656, was opened at 2008-05-12 21:36
Message generated for change (Comment added) made by keithmarshall
You can respond by visiting: 
https://sourceforge.net/tracker/?func=detail&atid=302435&aid=1962656&group_id=2435

Please note that this message will contain a full copy of the comment thread,
including the initial issue submission, for this request,
not just the latest update.
Category: None
Group: None
Status: Open
Resolution: None
Priority: 5
Private: No
Submitted By: Greg Chicares (chicares)
Assigned to: Nobody/Anonymous (nobody)
Summary: Fix round(2^53)

Initial Comment:
C99 round() failed with an argument of 2^53, because
adding 0.5 to that number caused overflow. With this
patch, it returns its argument if the argument is
already integral--as the C99 reference implementation
(F.9.6.6 paragraph 2) does. Both roundf() and roundl()
are patched similarly.

Speed is almost unaffected: a few percent slower for
non-integral arguments, and a couple percent faster
for integral arguments. The F.9.6.6 implementation
would run at about one-third this speed because it
saves and restores the whole FP environment.

With this patch, I get the same results as with this
NetBSD code:

http://archive.netbsd.se/?ml=freebsd-ports-bugs&a=2005-04&m=870458

(search for "round(double x)"), but that code runs
about ten percent slower; and a minimal change to
the 2004 mingwex sources seemed preferable anyway.

----------------------------------------------------------------------

>Comment By: Keith Marshall (keithmarshall)
Date: 2008-05-25 10:01

Message:
Logged In: YES 
user_id=823908
Originator: NO

I've attached an updated patch, against current CVS, which addresses the
rounding issue, as previously discussed, for all of {,l,ll}round{,f,l}()
functions.  This reproduces the correct results from my prior test case,
for each of the {,l,ll}round() functions, and also appears to DTRT for each
the other six; (testing on GNU/Linux, rounding is correct for _all_ cases,
in the vicinity of 0.5; however, I have not (yet) adjusted the extremes to
allow for differences in precision of the input data types).

File Added: alt-round.patch.gz

----------------------------------------------------------------------

Comment By: Keith Marshall (keithmarshall)
Date: 2008-05-22 12:20

Message:
Logged In: YES 
user_id=823908
Originator: NO

Greg wrote:
> IIRC, msvcrt's printf() assumes 'long double' is a synonym
> for 'double';

Yuk! While in mingw32, the two are distinct...

> but snprintf() should work, because we've implemented it
> via __mingw_snprintf().

All the more reason, IMO, why we should extend the __mingw_snprintf()
formatting engine for use in replacements for the entire gamut of printf
functions, (as I've suggested in
https://sourceforge.net/support/tracker.php?aid=1961893); native MSVCRT
printf doesn't appear to adequately support mingw32 data types.

But, that's a subject for another patch...

----------------------------------------------------------------------

Comment By: Greg Chicares (chicares)
Date: 2008-05-22 10:53

Message:
Logged In: YES 
user_id=231935
Originator: YES

> I don't seem to be able to printf a long double value on Win2K

IIRC, msvcrt's printf() assumes 'long double' is a synonym for
'double'; but snprintf() should work, because we've implemented
it via __mingw_snprintf().

----------------------------------------------------------------------

Comment By: Keith Marshall (keithmarshall)
Date: 2008-05-22 06:44

Message:
Logged In: YES 
user_id=823908
Originator: NO

Danny wrote:
> Keith's  approach in alt-round.tar.gz looks good then.  Is someone
> willing extend that to float and long double versions?  

Thanks Danny.  Yes, I'll follow up with that.  I already have a generic
implementation working for all nine of the affected functions; (it uses
your FE_TONEAREST technique, with a modification of Greg's proposed mid-way
adjustment).  I'll tidy it up, WRT comments, ChangeLog, etc., and post it
here for review this weekend; it reproduces the correct results from my
preceding test case, (with one caveat; I don't seem to be able to printf a
long double value on Win2K -- neither the `l' nor the `L' modifiers
documented at http://msdn.microsoft.com/de-de/library/tcxf1dw6.aspx, for
use with the `f' format spec, (but conspicuously not apparently with `e' or
`g'), seem to work as documented, for me).

----------------------------------------------------------------------

Comment By: Danny Smith (dannysmith)
Date: 2008-05-22 02:57

Message:
Logged In: YES 
user_id=11494
Originator: NO

Thanks Greg, I had forgotten about C99 vs i386 FE_TONEAREST difference.  
Keith's  approach in alt-round.tar.gz looks good then.  Is someone willing
extend that to float and long double versions?  

Danny

----------------------------------------------------------------------

Comment By: Greg Chicares (chicares)
Date: 2008-05-16 04:20

Message:
Logged In: YES 
user_id=231935
Originator: YES

Cloning trunc(), but setting FE_TONEAREST, is probably the fastest
and most reliable way to round to nearest, because all the work is
done by FRNDINT. But FE_TONEAREST resolves exact ties in favor of
the even integer, whereas C99's round() picks the integer farther
away from zero.

You could get C99 round() semantics by adding this at the bottom:

         if(0.0 < _x &&  0.5 == _x - res) return res + 1.0;
    else if(_x < 0.0 && -0.5 == _x - res) return res - 1.0;
    else                                  return res;

That seems to be only slightly faster than the trunc() approach.
Instead of adjusting half the real numbers, it adjusts only a
much smaller number of points; but it has to test every number,
and that overhead is about the same as with trunc().

----------------------------------------------------------------------

Comment By: Danny Smith (dannysmith)
Date: 2008-05-16 03:03

Message:
Logged In: YES 
user_id=11494
Originator: NO

If you look at the implementation of trunc , you will see that it just
uses a temporary change in FPU control word to get FE_TOWARDZERO.  So
rather than fiddle with trunc and FP addition, why not just do the obvious
for round, eg:

double
round (double _x){
  double res;
  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))
	    | FE_TONEAREST;
  __asm__ ("fldcw %0;" : : "m" (tmp_cw));
  __asm__ ("frndint;" : "=t" (res)  : "0" (_x)); /* round to nearest */
  __asm__ ("fldcw %0;" : : "m" (saved_cw) ); /* restore saved control word
*/
  return res;
}

long
lround (double _x)
{
  double res = round (_x);

  if (!isfinite (res) 
      || res > (double) LONG_MAX
      || res < (double) LONG_MIN)
    { 
      errno = ERANGE;
      /* Undefined behaviour, so we could return anything.  */
      /* return res > 0.0 ? LONG_MAX : LONG_MIN; */
    }
  return (long) res;
}

long long
llround (double _x)
{
  double res = round (_x);

  if (!isfinite (res) 
      || res > (double) LONG_LONG_MAX
      || res < (double) LONG_LONG_MIN)
    { 
      errno = ERANGE;
      /* Undefined behaviour, so we could return anything.  */
      /* return res > 0.0 ? LONG_LONG_MAX : LONG_LONG_MIN; */
    }
  return (long long) res;
}

This should be more efficient.

With that, and Keith test case I get:

Input value: 4.99999999999999940000e-001
round(x):    0.00000000000000000000e+000
llround(x):  0

Input value: -4.99999999999999940000e-001
round(x):    0.00000000000000000000e+000
llround(x):  0

Input value: 5.00000000000000000000e-001
round(x):    0.00000000000000000000e+000
llround(x):  0

Input value: -5.00000000000000000000e-001
round(x):    0.00000000000000000000e+000
llround(x):  0

Input value: 5.00000000000000110000e-001
round(x):    1.00000000000000000000e+000
llround(x):  1

Input value: -5.00000000000000110000e-001
round(x):    -1.00000000000000000000e+000
llround(x):  -1

Input value: 9.00719925474099100000e+015
round(x):    9.00719925474099100000e+015
llround(x):  9007199254740991

Input value: 4.50359962737049700000e+015
round(x):    4.50359962737049700000e+015
llround(x):  4503599627370497

Input value: -4.50359962737049700000e+015
round(x):    -4.50359962737049700000e+015
llround(x):  -4503599627370497

NB: round (5.00000000000000000000e-001) == 0.0, using the round-to-even
rule for 0.5

Danny

----------------------------------------------------------------------

Comment By: Greg Chicares (chicares)
Date: 2008-05-16 00:00

Message:
Logged In: YES 
user_id=231935
Originator: YES

I have no objection to anything in your 2008-05-14 12:59 posting.

It's sound engineering practice to factor out a core "helper" routine.
The one you've written is clear, fast, and passes all the tests I've
been developing. You seem to have the matter well in hand, so perhaps
it's best for me to focus on testing. I've long thought that we should
have unit tests for libmingwex, and the round() family is a good place
to start. I'll have little time before September to propose a general
framework worthy of consideration, but I'll make time to test anything
you write for round() and its kin.

----------------------------------------------------------------------

Comment By: Keith Marshall (keithmarshall)
Date: 2008-05-14 12:59

Message:
Logged In: YES 
user_id=823908
Originator: NO

Danny has proposed alternative implementations, which appear similar to
the NetBSD example cited.  To me, these seem unnecessarily verbose.  I've
attached an alternative reference implementation, for `round()' and
`llround()' only, together with a test case which seems to yield correct
results in each of the corner cases reported in these two related tickets.
File Added: alt-round.tar.gz

----------------------------------------------------------------------

Comment By: Keith Marshall (keithmarshall)
Date: 2008-05-13 13:26

Message:
Logged In: YES 
user_id=823908
Originator: NO

This fails for `round( nextafter( 0.5, 0.0 ));': it returns a value of
`1.0', where it should be `0.0'.

Please see also
https://sourceforge.net/tracker/index.php?func=detail&aid=1961860&group_id=2435&atid=302435

The two defects are related, and as Danny has suggested, require
consistent solutions.

----------------------------------------------------------------------

You can respond by visiting: 
https://sourceforge.net/tracker/?func=detail&atid=302435&aid=1962656&group_id=2435

[Mingw-notify] [ mingw-Patches-1962656 ] Fix round(2^53)

A native Windows port of the GNU Compiler Collection (GCC)

[Mingw-notify] [ mingw-Patches-1962656 ] Fix round(2^53)