From: SourceForge.net <no...@so...> - 2008-05-25 11:55:31
|
Patches item #1962656, was opened at 2008-05-12 17:36 Message generated for change (Comment added) made by ir0nh34d 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: Chris Sutcliffe (ir0nh34d) Date: 2008-05-25 07:55 Message: Logged In: YES user_id=570619 Originator: NO Keith, should I hold off applying the patch and wait for the patch that addesses the adjustment of extremes, or is this good to go? ---------------------------------------------------------------------- Comment By: Keith Marshall (keithmarshall) Date: 2008-05-25 06: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 08: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 06: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 02: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-21 22: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 00: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-15 23: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-15 20: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 08: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 09: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 |