From: SourceForge.net <no...@so...> - 2008-05-16 03:03:54
|
Patches item #1962656, was opened at 2008-05-13 09:36 Message generated for change (Comment added) made by dannysmith 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: Danny Smith (dannysmith) Date: 2008-05-16 15: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 12: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-15 00: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-14 01: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 |