You can subscribe to this list here.
2000 
_{Jan}

_{Feb}

_{Mar}

_{Apr}

_{May}

_{Jun}

_{Jul}

_{Aug}

_{Sep}

_{Oct}

_{Nov}
(4) 
_{Dec}
(47) 

2001 
_{Jan}
(49) 
_{Feb}
(45) 
_{Mar}
(35) 
_{Apr}
(75) 
_{May}
(30) 
_{Jun}
(67) 
_{Jul}
(53) 
_{Aug}
(70) 
_{Sep}
(33) 
_{Oct}
(30) 
_{Nov}
(21) 
_{Dec}
(29) 
2002 
_{Jan}
(43) 
_{Feb}
(28) 
_{Mar}
(43) 
_{Apr}
(23) 
_{May}
(86) 
_{Jun}
(67) 
_{Jul}
(55) 
_{Aug}
(116) 
_{Sep}
(87) 
_{Oct}
(27) 
_{Nov}
(48) 
_{Dec}
(93) 
2003 
_{Jan}
(122) 
_{Feb}
(139) 
_{Mar}
(170) 
_{Apr}
(46) 
_{May}
(84) 
_{Jun}
(60) 
_{Jul}
(60) 
_{Aug}
(86) 
_{Sep}
(106) 
_{Oct}
(42) 
_{Nov}
(24) 
_{Dec}
(43) 
2004 
_{Jan}
(63) 
_{Feb}
(134) 
_{Mar}
(95) 
_{Apr}
(98) 
_{May}
(85) 
_{Jun}
(44) 
_{Jul}
(142) 
_{Aug}
(71) 
_{Sep}
(45) 
_{Oct}
(88) 
_{Nov}
(46) 
_{Dec}
(50) 
2005 
_{Jan}
(100) 
_{Feb}
(72) 
_{Mar}
(71) 
_{Apr}
(55) 
_{May}
(76) 
_{Jun}
(86) 
_{Jul}
(158) 
_{Aug}
(142) 
_{Sep}
(51) 
_{Oct}
(37) 
_{Nov}
(59) 
_{Dec}
(79) 
2006 
_{Jan}
(61) 
_{Feb}
(34) 
_{Mar}
(95) 
_{Apr}
(170) 
_{May}
(66) 
_{Jun}
(37) 
_{Jul}
(29) 
_{Aug}
(28) 
_{Sep}
(59) 
_{Oct}
(48) 
_{Nov}
(72) 
_{Dec}
(50) 
2007 
_{Jan}
(68) 
_{Feb}
(49) 
_{Mar}
(38) 
_{Apr}
(79) 
_{May}
(63) 
_{Jun}
(29) 
_{Jul}
(64) 
_{Aug}
(47) 
_{Sep}
(67) 
_{Oct}
(101) 
_{Nov}
(42) 
_{Dec}
(29) 
2008 
_{Jan}
(37) 
_{Feb}
(44) 
_{Mar}
(64) 
_{Apr}
(87) 
_{May}
(132) 
_{Jun}
(92) 
_{Jul}
(135) 
_{Aug}
(70) 
_{Sep}
(72) 
_{Oct}
(30) 
_{Nov}
(21) 
_{Dec}
(32) 
2009 
_{Jan}
(101) 
_{Feb}
(65) 
_{Mar}
(82) 
_{Apr}
(38) 
_{May}
(29) 
_{Jun}
(75) 
_{Jul}
(70) 
_{Aug}
(69) 
_{Sep}
(82) 
_{Oct}
(28) 
_{Nov}
(51) 
_{Dec}
(19) 
2010 
_{Jan}
(46) 
_{Feb}
(67) 
_{Mar}
(66) 
_{Apr}
(54) 
_{May}
(55) 
_{Jun}
(50) 
_{Jul}
(84) 
_{Aug}
(86) 
_{Sep}
(43) 
_{Oct}
(63) 
_{Nov}
(33) 
_{Dec}
(27) 
2011 
_{Jan}
(70) 
_{Feb}
(29) 
_{Mar}
(54) 
_{Apr}
(50) 
_{May}
(105) 
_{Jun}
(45) 
_{Jul}
(30) 
_{Aug}
(83) 
_{Sep}
(38) 
_{Oct}
(71) 
_{Nov}
(124) 
_{Dec}
(61) 
2012 
_{Jan}
(33) 
_{Feb}
(37) 
_{Mar}
(60) 
_{Apr}
(60) 
_{May}
(51) 
_{Jun}
(137) 
_{Jul}
(80) 
_{Aug}
(156) 
_{Sep}
(32) 
_{Oct}
(168) 
_{Nov}
(56) 
_{Dec}
(50) 
2013 
_{Jan}
(54) 
_{Feb}

_{Mar}

_{Apr}

_{May}

_{Jun}

_{Jul}

_{Aug}

_{Sep}

_{Oct}

_{Nov}

_{Dec}

S  M  T  W  T  F  S 





1
(6) 
2
(9) 
3
(3) 
4
(6) 
5
(4) 
6
(4) 
7
(7) 
8
(15) 
9
(6) 
10

11
(9) 
12
(5) 
13
(3) 
14
(6) 
15
(5) 
16
(8) 
17
(2) 
18

19

20

21
(2) 
22
(9) 
23

24

25
(5) 
26

27
(1) 
28
(4) 
29
(5) 
30
(8) 
31

From: SourceForge.net <noreply@so...>  20080525 18:11:07

Patches item #1962656, was opened at 20080512 21:36 Message generated for change (Comment added) made by chicares 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 integralas 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 nonintegral arguments, and a couple percent faster for integral arguments. The F.9.6.6 implementation would run at about onethird 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=freebsdportsbugs&a=200504&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: Greg Chicares (chicares) Date: 20080525 18:11 Message: Logged In: YES user_id=231935 Originator: YES > I wanted to give > Danny, Greg, FX and others an opportunity to comment first. I'll try to make time for testing soon, but it'll probably take until Tuesday before I can do that. Monday's a national holiday here, and right now I'm wrestling with a production problem in my day job. > For {,l,ll}roundf() the overflow > case equivalent to 2^53 would be 2^24, (I believe), and I think it > should be 2^64 for {,l,ll}roundl() I think that's all correct: $echo '#include <float.h>'  cpp dM  grep MANT grep '[19]' #define __FLT_MANT_DIG__ 24 #define __LDBL_MANT_DIG__ 64 #define __DBL_MANT_DIG__ 53 > corner case values suggested by FX, I don't know his basis for > those particular choices I guess 4.99999999999999944488848768742172978818416595459e1L is the nextlower neighbor of onehalf; I think he got the other one 4.503599627370497e+15L from this message: http://gcc.gnu.org/ml/gccpatches/200610/msg00917.html which gives no exact rationale for it, but it appears to be 1+2^52.  Comment By: Keith Marshall (keithmarshall) Date: 20080525 13:53 Message: Logged In: YES user_id=823908 Originator: NO Chris, I was, and still am, willing to commit this myself, but I wanted to give Danny, Greg, FX and others an opportunity to comment first. I'm hoping no further patch will be required. I haven't tested with any other input values than x = +/0.5 +/ nextafter(0.5, 0.0), (for which all cases yield correct results), Greg's 2^53 value, and the two corner case values for which FX reported a problem. For {,l,ll}roundf() the overflow case equivalent to 2^53 would be 2^24, (I believe), and I think it should be 2^64 for {,l,ll}roundl(), (although I haven't found any definitive confirmation that 80bit long doubles do use a 64bit mantissa, but http://tinyurl.com/3jn5hl hints that it may be so). For the two corner case values suggested by FX, I don't know his basis for those particular choices, so it isn't clear to me how to choose equivalents for the float and long double cases.  Comment By: Chris Sutcliffe (ir0nh34d) Date: 20080525 11: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: 20080525 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: altround.patch.gz  Comment By: Keith Marshall (keithmarshall) Date: 20080522 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: 20080522 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: 20080522 06:44 Message: Logged In: YES user_id=823908 Originator: NO Danny wrote: > Keith's approach in altround.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 midway 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/dede/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: 20080522 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 altround.tar.gz looks good then. Is someone willing extend that to float and long double versions? Danny  Comment By: Greg Chicares (chicares) Date: 20080516 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: 20080516 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.99999999999999940000e001 round(x): 0.00000000000000000000e+000 llround(x): 0 Input value: 4.99999999999999940000e001 round(x): 0.00000000000000000000e+000 llround(x): 0 Input value: 5.00000000000000000000e001 round(x): 0.00000000000000000000e+000 llround(x): 0 Input value: 5.00000000000000000000e001 round(x): 0.00000000000000000000e+000 llround(x): 0 Input value: 5.00000000000000110000e001 round(x): 1.00000000000000000000e+000 llround(x): 1 Input value: 5.00000000000000110000e001 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.00000000000000000000e001) == 0.0, using the roundtoeven rule for 0.5 Danny  Comment By: Greg Chicares (chicares) Date: 20080516 00:00 Message: Logged In: YES user_id=231935 Originator: YES I have no objection to anything in your 20080514 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: 20080514 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: altround.tar.gz  Comment By: Keith Marshall (keithmarshall) Date: 20080513 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 
From: SourceForge.net <noreply@so...>  20080525 13:59:03

Patches item #1961860, was opened at 20080511 12:30 Message generated for change (Comment added) made by keithmarshall You can respond by visiting: https://sourceforge.net/tracker/?func=detail&atid=302435&aid=1961860&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: runtime Group: None Status: Open Resolution: None Priority: 5 Private: No Submitted By: FrançoisXavier Coudert (coudert) Assigned to: Nobody/Anonymous (nobody) Summary: lround and llround functions are incorrect Initial Comment: mingw lround and llround functions (and their float and long double variants) give incorrect results for some inputs, including nextafter(0.5,1.) (the largest number smaller than 0.5) and 4.503599627370497e+15. This is due to the algorithm used, which is to add +/ 0.5 and round towards zero. The fix is to add not +/ 0.5, but +/ nextafter(0.5,1.). See http://gcc.gnu.org/ml/gccpatches/200610/msg00917.html and followups for a good explaination of the issues involved. I spotted this issue because it causes failures in the gfortran testsuite. The patch attached does that change. Please consider it for the next mingwruntime release. PS: a short program showing both the failure and the fix working is: #include <math.h> #include <limits.h> #include <errno.h> #include <stdio.h> #include <stdlib.h> long long my_llround (double x) { /* Add +/ 0.5 then then round towards zero. */ double tmp = trunc (x + (x >= 0.0 ? 0.5 : 0.5)); if (!isfinite (tmp)  tmp > (double)LONG_LONG_MAX  tmp < (double)LONG_LONG_MIN) { errno = ERANGE; /* Undefined behaviour, so we could return anything. */ /* return tmp > 0.0 ? LONG_LONG_MAX : LONG_LONG_MIN; */ } return (long long)tmp; } long long my_llround_fixed (double x) { double c; c = nextafter (0.5, 1.); /* Add +/ 0.5 then then round towards zero. */ double tmp = trunc (x + (x >= 0.0 ? c : c)); if (!isfinite (tmp)  tmp > (double)LONG_LONG_MAX  tmp < (double)LONG_LONG_MIN) { errno = ERANGE; /* Undefined behaviour, so we could return anything. */ /* return tmp > 0.0 ? LONG_LONG_MAX : LONG_LONG_MIN; */ } return (long long)tmp; } int main (void) { double a; long long int i1, i2; a = 4.99999999999999944488848768742172978818416595459e1L; printf ("%lli %lli %lli\n", llround(a), my_llround(a), my_llround_fixed(a)); a = 4.503599627370497e+15L; printf ("%lli %lli %lli\n", llround(a), my_llround(a), my_llround_fixed(a)); a = 4.503599627370497e+15L; printf ("%lli %lli %lli\n", llround(a), my_llround(a), my_llround_fixed(a)); return 0; } my_llround() is the current mingw implementation, my_llround_fixed() is the one proposed in my patch. If you run it on a system with correct native llround(), you get: 0 1 0 4503599627370497 4503599627370498 4503599627370497 4503599627370497 4503599627370498 4503599627370497 where the second column (current mingw) differs from the first (correct llround) and third (patched mingw).  >Comment By: Keith Marshall (keithmarshall) Date: 20080525 13:59 Message: Logged In: YES user_id=823908 Originator: NO I've attached a patch, (http://sourceforge.net/tracker/download.php?group_id=2435&atid=302435&file_id=278995&aid=1962656), to http://sourceforge.net/tracker/index.php?func=detail&aid=1962656&group_id=2435&atid=302435, which should fix this. As I've noted there, it isn't clear to me, how I should choose corresponding input values, to create equivalent test cases for the float and long double parameter function variants.  Comment By: FrançoisXavier Coudert (coudert) Date: 20080514 13:48 Message: Logged In: YES user_id=865300 Originator: YES A slightly more detailed explanation is there: http://gcc.gnu.org/ml/fortran/200504/msg00139.html The algorithm is used in glibc for some targets; see e.g. glibc/sysdeps/alpha/fpu/s_lround.c: long int __lround (double x) { double adj; adj = 0x1.fffffffffffffp2; /* nextafter (0.5, 0.0) */ adj = copysign (adj, x); return x + adj; } But anyway, as maintainers, it's your choice to what version you feel the most at ease with, so I'll leave you at it. Thanks!  Comment By: Keith Marshall (keithmarshall) Date: 20080514 13:30 Message: Logged In: YES user_id=823908 Originator: NO I know you are trying to deal with imprecision in the floating point representation, but to me your proposed solution seems mathematically suspect; the more coercions between types of differing precisions you introduce, the more susceptible the implementation becomes to loss of precision in the result. I read the thread you pointed to, and all I see is a glib statement, with no supporting mathematical basis or proof, that this `nextafter' technique is the right way to do it; without such a proof, I find it very difficult to accept that hypothesis. Danny's approach definitely appears to be more robust. However, it also seems rather verbose. For a more compact solution, which appears to yield the correct results for each of your test inputs, and also in the limiting case near the overflow limits, please see the reference implementations of `round()' and `llround()', which I've attached to the other ticket I referred to previously.  Comment By: FrançoisXavier Coudert (coudert) Date: 20080514 09:25 Message: Logged In: YES user_id=865300 Originator: YES I'll reply only to the main issue, mainly that the patch I proposed doesn't behave as expected. It indeed does not work on linux while the very same testcase works fine on Darwin, where I initially tested. I think it's an issue of extra precision, and this can be fixed by performing calculations in the widest available FP mode (here, long double): long long my_llround_fixed (double x) { long double c = nextafter (0.5L, 1.L); /* Add +/ 0.5 then then round towards zero. */ double tmp = (double) truncl((long double) x + (x >= 0.0 ? c : c)); if (!isfinite (tmp)  tmp > (double)LLONG_MAX  tmp < (double)LLONG_MIN) { errno = ERANGE; /* Undefined behaviour, so we could return anything. */ /* return tmp > 0.0 ? LLONG_MAX : LLONG_MIN; */ } return (long long)tmp; } I mention this only for completeness, as it seems Danny has proposed a better patch.  Comment By: Keith Marshall (keithmarshall) Date: 20080514 09:07 Message: Logged In: YES user_id=823908 Originator: NO I have a number of issues with this patch; not all are related to the subject matter, but I'll list them anyway, noting those which are not immediately pertinent as "side issues". 1) The supplied test case is not valid; to fix it, I had to `#include <inttypes.h>' and replace all of the `%lli' printf format specs with `%"PRIi64"', since `%lli' is not supported by MSVCRT printf. I also had to replace the nonstandard use of `LONG_LONG_MAX' and `LONG_LONG_MIN' with their standards conforming `LLONG_MAX' and `LLONG_MIN' counterparts; contrary to the comment in MinGW's limits.h, GCC on my Ubuntu 6.06 LTS box does *not* support this nonstandard extension. (The issue of nonstandard LONG_LONG_MAX vs. standard LLONG_MAX is, of course, a side issue in the context of this patch). 2) After fixing your test case, as described in (1), I have compiled and run it, *both* with MinGW on a Win2K box, *and* in GNU/Linux on the Ubuntu box; in *neither* case, can I reproduce your claimed results: 2.1) MinGW GCC3.4.5 on Win2K... $ gcc testcase.c $ ./a 1 1 0 4503599627370498 4503599627370498 4503599627370498 4503599627370498 4503599627370498 4503599627370498 2.2) GCC4.0.3 on Ubuntu 6.06 LTS $ gcc std=c99 testcase.c lm $ ./a.out 0 1 0 4503599627370497 4503599627370498 4503599627370498 4503599627370497 4503599627370498 4503599627370498 Note that, while your patch does produce the correct result for rounding the value of `nextafter( 0.5, 0.0 )', it does *not* appear to reproducibly round correctly, (to the value you claim to be correct), in the other two cases you cite; the only function which does that, for my running of your test case, is the system supplied `llround()' on the Ubuntu box. 3) On the Ubuntu box, I have to add the `std=c99' compiler option, else the test case will not compile; this is as it should be. However, if I add that switch for the MinGW case, then the use of LONG_LONG_MAX and LONG_LONG_MIN become invalid, and the code will not compile. Again, this is a side issue, and it goes away when LLONG_MAX and LLONG_MIN are used instead; however it still leaves the anomaly that MinGW does not require `std=c99', where perhaps it really should. 4) Another (minor) side issue: when the rounded value lies beyond the limits of a long long int, the MinGW implementation sets `errno = ERANGE'. In this circumstance, POSIX requires that it sets `errno = EDOM'; there is no provision for any of the l[l]round*() functions to ever set `errno = ERANGE'. I've not explored Danny's proposed alternative implementation, and while my own evaluation would seem to refute his statement that "your patch works...", I concur with his opinion that this issue needs to be resolved in a manner which simultaneously addresses https://sourceforge.net/tracker/index.php?func=detail&aid=1962656&group_id=2435&atid=302435  Comment By: Danny Smith (dannysmith) Date: 20080512 10:07 Message: Logged In: YES user_id=11494 Originator: NO coudert wrote: "I don't understand what you mean by 'Your proposed patch works but ... doesn't extend itself to fix problems with round{fl} in corner cases...' It would ease maintenance if round[fl] use the same basic algorithm as l[l]round[fl] so that a fix for you would also avoid the overflow problems with round http://listsarchives.org/mingwusers/10476aboutroundfunction.html. That is, I am referring to cases of roundl (float_value_near_LONG_DOUBLE_MAX_that_overflows_to_INF_when_you_add_4.99999999999999944488848768742172978818416595459e1L) That patch would look something like attached round.diff (testing in progress). Oops, sorry I am not allowed to attach anyhing, so here is the patch inline: Index: mingwex/math/llround.c =================================================================== RCS file: /cvs/src/src/winsup/mingw/mingwex/math/llround.c,v retrieving revision 1.3 diff c 3 p r1.3 llround.c *** mingwex/math/llround.c 22 Apr 2004 04:38:24 0000 1.3  mingwex/math/llround.c 12 May 2008 09:58:53 0000 *************** *** 5,19 **** long long llround (double x) { ! /* Add +/ 0.5, then round towards zero. */ ! double tmp = trunc (x + (x >= 0.0 ? 0.5 : 0.5)); ! if (!isfinite (tmp) !  tmp > (double)LONG_LONG_MAX !  tmp < (double)LONG_LONG_MIN) { errno = ERANGE; /* Undefined behaviour, so we could return anything. */ ! /* return tmp > 0.0 ? LONG_LONG_MAX : LONG_LONG_MIN; */ } ! return (long long)tmp; }  5,31  long long llround (double x) { ! double res; ! if (x >= 0.0) ! { ! res = ceil (x); ! if (res  x > 0.5) ! res = 1.0; ! } ! else ! { ! res = ceil (x); ! if (res + x > 0.5) ! res = 1.0; ! res = res;; ! } ! 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; } Index: mingwex/math/llroundf.c =================================================================== RCS file: /cvs/src/src/winsup/mingw/mingwex/math/llroundf.c,v retrieving revision 1.2 diff c 3 p r1.2 llroundf.c *** mingwex/math/llroundf.c 22 Apr 2004 04:38:24 0000 1.2  mingwex/math/llroundf.c 12 May 2008 09:58:53 0000 *************** *** 5,19 **** long long llroundf (float x) { ! /* Add +/ 0.5, then round towards zero. */ ! float tmp = truncf (x + (x >= 0.0F ? 0.5F : 0.5F)); ! if (!isfinite (tmp) !  tmp > (float)LONG_LONG_MAX !  tmp < (float)LONG_LONG_MIN) { errno = ERANGE; /* Undefined behaviour, so we could return anything. */ ! /* return tmp > 0.0F ? LONG_LONG_MAX : LONG_LONG_MIN; */ } ! return (long long)tmp; ! }  5,31  long long llroundf (float x) { ! float res; ! if (x >= 0.0F) ! { ! res = ceilf (x); ! if (res  x > 0.5F) ! res = 1.0F; ! } ! else ! { ! res = ceilf (x); ! if (res + x > 0.5F) ! res = 1.0F; ! res = res; ! } ! if (!isfinite (res) !  res > (float) LONG_LONG_MAX !  res < (float) LONG_LONG_MIN) { errno = ERANGE; /* Undefined behaviour, so we could return anything. */ ! /* return res > 0.0F ? LONG_LONG_MAX : LONG_LONG_MIN; */ } ! return (long long) res; ! } Index: mingwex/math/llroundl.c =================================================================== RCS file: /cvs/src/src/winsup/mingw/mingwex/math/llroundl.c,v retrieving revision 1.2 diff c 3 p r1.2 llroundl.c *** mingwex/math/llroundl.c 22 Apr 2004 04:38:24 0000 1.2  mingwex/math/llroundl.c 12 May 2008 09:58:53 0000 *************** *** 5,19 **** long long llroundl (long double x) { ! /* Add +/ 0.5, then round towards zero. */ ! long double tmp = truncl (x + (x >= 0.0L ? 0.5L : 0.5L)); ! if (!isfinite (tmp) !  tmp > (long double)LONG_LONG_MAX !  tmp < (long double)LONG_LONG_MIN) { errno = ERANGE; /* Undefined behaviour, so we could return anything. */ ! /* return tmp > 0.0L ? LONG_LONG_MAX : LONG_LONG_MIN; */ } ! return (long long)tmp; }  5,31  long long llroundl (long double x) { ! long double res; ! if (x >= 0.0L) ! { ! res = ceill (x); ! if (res  x > 0.5L) ! res = 1.0L; ! } ! else ! { ! res = ceill (x); ! if (res + x > 0.5L) ! res = 1.0L; ! res = res; ! } ! 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; } Index: mingwex/math/lround.c =================================================================== RCS file: /cvs/src/src/winsup/mingw/mingwex/math/lround.c,v retrieving revision 1.2 diff c 3 p r1.2 lround.c *** mingwex/math/lround.c 22 Apr 2004 04:38:24 0000 1.2  mingwex/math/lround.c 12 May 2008 09:58:53 0000 *************** *** 1,19 **** #include <math.h> #include <limits.h> #include <errno.h> long ! lround (double x) { ! /* Add +/ 0.5 then then round towards zero. */ ! double tmp = trunc (x + (x >= 0.0 ? 0.5 : 0.5)); ! if (!isfinite (tmp) !  tmp > (double)LONG_MAX !  tmp < (double)LONG_MIN) { errno = ERANGE; /* Undefined behaviour, so we could return anything. */ ! /* return tmp > 0.0 ? LONG_MAX : LONG_MIN; */ } ! return (long)tmp; }  1,33  #include <math.h> #include <limits.h> #include <errno.h> + #include <math.h> long ! lround (double x) { ! double res; ! if (x >= 0.0) ! { ! res = ceil (x); ! if (res  x > 0.5) ! res = 1.0; ! } ! else ! { ! res = ceil (x); ! if (res + x > 0.5) ! res = 1.0; ! res = res; ! } ! ! 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; } Index: mingwex/math/lroundf.c =================================================================== RCS file: /cvs/src/src/winsup/mingw/mingwex/math/lroundf.c,v retrieving revision 1.2 diff c 3 p r1.2 lroundf.c *** mingwex/math/lroundf.c 22 Apr 2004 04:38:24 0000 1.2  mingwex/math/lroundf.c 12 May 2008 09:58:53 0000 *************** *** 5,19 **** long lroundf (float x) { ! /* Add +/ 0.5, then round towards zero. */ ! float tmp = truncf (x + (x >= 0.0F ? 0.5F : 0.5F)); ! if (!isfinite (tmp) !  tmp > (float)LONG_MAX !  tmp < (float)LONG_MIN) { errno = ERANGE; /* Undefined behaviour, so we could return anything. */ ! /* return tmp > 0.0F ? LONG_MAX : LONG_MIN; */ } ! return (long)tmp; ! }  5,31  long lroundf (float x) { ! float res; ! if (x >= 0.0F) ! { ! res = ceilf (x); ! if (res  x > 0.5F) ! res = 1.0F; ! } ! else ! { ! res = ceilf (x); ! if (res + x > 0.5F) ! res = 1.0F; ! res = res; ! } ! if (!isfinite (res) !  res > (float) LONG_MAX !  res < (float) LONG_MIN) { errno = ERANGE; /* Undefined behaviour, so we could return anything. */ ! /* return res > 0.0F ? LONG_MAX : LONG_MIN; */ } ! return (long) res; ! } Index: mingwex/math/lroundl.c =================================================================== RCS file: /cvs/src/src/winsup/mingw/mingwex/math/lroundl.c,v retrieving revision 1.2 diff c 3 p r1.2 lroundl.c *** mingwex/math/lroundl.c 22 Apr 2004 04:38:24 0000 1.2  mingwex/math/lroundl.c 12 May 2008 09:58:53 0000 *************** *** 3,19 **** #include <errno.h> long ! lroundl (long double x) { ! /* Add +/ 0.5, then round towards zero. */ ! long double tmp = truncl (x + (x >= 0.0L ? 0.5L : 0.5L)); ! if (!isfinite (tmp) !  tmp > (long double)LONG_MAX !  tmp < (long double)LONG_MIN) { errno = ERANGE; /* Undefined behaviour, so we could return anything. */ ! /* return tmp > 0.0L ? LONG_MAX : LONG_MIN; */ } ! return (long)tmp; }  3,32  #include <errno.h> long ! lroundl (long double x) { ! long double res; ! if (x >= 0.0L) ! { ! res = ceill (x); ! if (res  x > 0.5L) ! res = 1.0; ! } ! else ! { ! res = ceill (x); ! if (res + x > 0.5L) ! res = 1.0L; ! res = res;; ! } ! if (!isfinite (res) !  res > (long double)LONG_MAX !  res < (long double)LONG_MIN) { errno = ERANGE; /* Undefined behaviour, so we could return anything. */ ! /* return res > 0.0L ? LONG_MAX : LONG_MIN; */ } ! return (long)res; } + Index: mingwex/math/round.c =================================================================== RCS file: /cvs/src/src/winsup/mingw/mingwex/math/round.c,v retrieving revision 1.2 diff c 3 p r1.2 round.c *** mingwex/math/round.c 29 Mar 2004 08:22:19 0000 1.2  mingwex/math/round.c 12 May 2008 09:58:53 0000 *************** *** 1,8 **** #include <math.h> double ! round (double x) { ! /* Add +/ 0.5 then then round towards zero. */ ! return trunc ( x + (x >= 0.0 ? 0.5 : 0.5)); }  1,22  #include <math.h> double ! round (double x) { ! double res; ! if (x >= 0.0) ! { ! res = ceil (x); ! if (res  x > 0.5) ! res = 1.0; ! } ! else ! { ! res = ceil (x); ! if (res + x > 0.5) ! res = 1.0; ! res = res; ! } ! return res; } + Index: mingwex/math/roundf.c =================================================================== RCS file: /cvs/src/src/winsup/mingw/mingwex/math/roundf.c,v retrieving revision 1.2 diff c 3 p r1.2 roundf.c *** mingwex/math/roundf.c 29 Mar 2004 08:22:20 0000 1.2  mingwex/math/roundf.c 12 May 2008 09:58:53 0000 *************** *** 3,8 **** float roundf (float x) { ! /* Add +/ 0.5 then then round towards zero. */ ! return truncf ( x + (x >= 0.0F ? 0.5F : 0.5F)); }  3,22  float roundf (float x) { ! float res; ! if (x >= 0.0f) ! { ! res = ceilf (x); ! if (res  x > 0.5f) ! res = 1.0f; ! } ! else ! { ! res = ceilf (x); ! if (res + x > 0.5f) ! res = 1.0f; ! res = res; ! } ! return res; } + Index: mingwex/math/roundl.c =================================================================== RCS file: /cvs/src/src/winsup/mingw/mingwex/math/roundl.c,v retrieving revision 1.2 diff c 3 p r1.2 roundl.c *** mingwex/math/roundl.c 29 Mar 2004 08:22:20 0000 1.2  mingwex/math/roundl.c 12 May 2008 09:58:53 0000 *************** *** 3,8 **** long double roundl (long double x) { ! /* Add +/ 0.5 then then round towards zero. */ ! return truncl ( x + (x >= 0.0L ? 0.5L : 0.5L)); }  3,22  long double roundl (long double x) { ! long double res; ! if (x >= 0.0L) ! { ! res = ceill (x); ! if (res  x > 0.5L) ! res = 1.0L; ! } ! else ! { ! res = ceill (x); ! if (res + x > 0.5L) ! res = 1.0L; ! res = res;; ! } ! return res; } +  Comment By: FrançoisXavier Coudert (coudert) Date: 20080512 09:14 Message: Logged In: YES user_id=865300 Originator: YES I don't understand what you mean by "Your proposed patch works but ....it doesn't extend itself to fix problems with round{fl} in corner cases where we are near [FLOATDOUBLELONG_DOUBLE]_MAX." Are you referring to 4.503599627370497e+15L? Because the patch I propose does fix it AFAICT. Can you post an example case where it fails?  Comment By: Danny Smith (dannysmith) Date: 20080512 08:06 Message: Logged In: YES user_id=11494 Originator: NO Your proposed patch works but ....it doesn't extend itself to fix problems with round{fl} in corner cases where we are near [FLOATDOUBLELONG_DOUBLE]_MAX. I think something on this track that uses ceil (maybe floor is more efficient, I haven't tested): #include <math.h> double __round (double x) { double res; if (x >= 0.0) { res = ceil (x); if (res  x > 0.5) res = 1.0; } else { res = ceil (x); if (res + x > 0.5) res = 1.0; res = res;; } return res; } #define __lround(x) (long)__round(x) #define __llround(x) (long long)__round(x) int main (void) { double a; a = 4.99999999999999944488848768742172978818416595459e1L; printf ("%I64i %I64i\n", llround(a), __llround(a)); a = 4.503599627370497e+15L; printf ("%I64i %I64i\n", llround(a), __llround(a)); a = 4.503599627370497e+15L; printf ("%I64i %I64i\n", llround(a), __llround(a)); return 0; }  Comment By: FrançoisXavier Coudert (coudert) Date: 20080511 12:36 Message: Logged In: YES user_id=865300 Originator: YES I should have added there's the same issue in round(), roundf() and roundl(). The same solution also applies there.  You can respond by visiting: https://sourceforge.net/tracker/?func=detail&atid=302435&aid=1961860&group_id=2435 
From: SourceForge.net <noreply@so...>  20080525 13:53:20

Patches item #1962656, was opened at 20080512 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 integralas 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 nonintegral arguments, and a couple percent faster for integral arguments. The F.9.6.6 implementation would run at about onethird 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=freebsdportsbugs&a=200504&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: 20080525 13:53 Message: Logged In: YES user_id=823908 Originator: NO Chris, I was, and still am, willing to commit this myself, but I wanted to give Danny, Greg, FX and others an opportunity to comment first. I'm hoping no further patch will be required. I haven't tested with any other input values than x = +/0.5 +/ nextafter(0.5, 0.0), (for which all cases yield correct results), Greg's 2^53 value, and the two corner case values for which FX reported a problem. For {,l,ll}roundf() the overflow case equivalent to 2^53 would be 2^24, (I believe), and I think it should be 2^64 for {,l,ll}roundl(), (although I haven't found any definitive confirmation that 80bit long doubles do use a 64bit mantissa, but http://tinyurl.com/3jn5hl hints that it may be so). For the two corner case values suggested by FX, I don't know his basis for those particular choices, so it isn't clear to me how to choose equivalents for the float and long double cases.  Comment By: Chris Sutcliffe (ir0nh34d) Date: 20080525 11: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: 20080525 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: altround.patch.gz  Comment By: Keith Marshall (keithmarshall) Date: 20080522 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: 20080522 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: 20080522 06:44 Message: Logged In: YES user_id=823908 Originator: NO Danny wrote: > Keith's approach in altround.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 midway 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/dede/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: 20080522 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 altround.tar.gz looks good then. Is someone willing extend that to float and long double versions? Danny  Comment By: Greg Chicares (chicares) Date: 20080516 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: 20080516 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.99999999999999940000e001 round(x): 0.00000000000000000000e+000 llround(x): 0 Input value: 4.99999999999999940000e001 round(x): 0.00000000000000000000e+000 llround(x): 0 Input value: 5.00000000000000000000e001 round(x): 0.00000000000000000000e+000 llround(x): 0 Input value: 5.00000000000000000000e001 round(x): 0.00000000000000000000e+000 llround(x): 0 Input value: 5.00000000000000110000e001 round(x): 1.00000000000000000000e+000 llround(x): 1 Input value: 5.00000000000000110000e001 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.00000000000000000000e001) == 0.0, using the roundtoeven rule for 0.5 Danny  Comment By: Greg Chicares (chicares) Date: 20080516 00:00 Message: Logged In: YES user_id=231935 Originator: YES I have no objection to anything in your 20080514 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: 20080514 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: altround.tar.gz  Comment By: Keith Marshall (keithmarshall) Date: 20080513 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 
From: SourceForge.net <noreply@so...>  20080525 11:55:31

Patches item #1962656, was opened at 20080512 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 integralas 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 nonintegral arguments, and a couple percent faster for integral arguments. The F.9.6.6 implementation would run at about onethird 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=freebsdportsbugs&a=200504&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: 20080525 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: 20080525 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: altround.patch.gz  Comment By: Keith Marshall (keithmarshall) Date: 20080522 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: 20080522 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: 20080522 02:44 Message: Logged In: YES user_id=823908 Originator: NO Danny wrote: > Keith's approach in altround.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 midway 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/dede/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: 20080521 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 altround.tar.gz looks good then. Is someone willing extend that to float and long double versions? Danny  Comment By: Greg Chicares (chicares) Date: 20080516 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: 20080515 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.99999999999999940000e001 round(x): 0.00000000000000000000e+000 llround(x): 0 Input value: 4.99999999999999940000e001 round(x): 0.00000000000000000000e+000 llround(x): 0 Input value: 5.00000000000000000000e001 round(x): 0.00000000000000000000e+000 llround(x): 0 Input value: 5.00000000000000000000e001 round(x): 0.00000000000000000000e+000 llround(x): 0 Input value: 5.00000000000000110000e001 round(x): 1.00000000000000000000e+000 llround(x): 1 Input value: 5.00000000000000110000e001 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.00000000000000000000e001) == 0.0, using the roundtoeven rule for 0.5 Danny  Comment By: Greg Chicares (chicares) Date: 20080515 20:00 Message: Logged In: YES user_id=231935 Originator: YES I have no objection to anything in your 20080514 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: 20080514 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: altround.tar.gz  Comment By: Keith Marshall (keithmarshall) Date: 20080513 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 
From: SourceForge.net <noreply@so...>  20080525 10:01:58

Patches item #1962656, was opened at 20080512 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 integralas 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 nonintegral arguments, and a couple percent faster for integral arguments. The F.9.6.6 implementation would run at about onethird 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=freebsdportsbugs&a=200504&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: 20080525 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: altround.patch.gz  Comment By: Keith Marshall (keithmarshall) Date: 20080522 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: 20080522 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: 20080522 06:44 Message: Logged In: YES user_id=823908 Originator: NO Danny wrote: > Keith's approach in altround.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 midway 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/dede/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: 20080522 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 altround.tar.gz looks good then. Is someone willing extend that to float and long double versions? Danny  Comment By: Greg Chicares (chicares) Date: 20080516 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: 20080516 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.99999999999999940000e001 round(x): 0.00000000000000000000e+000 llround(x): 0 Input value: 4.99999999999999940000e001 round(x): 0.00000000000000000000e+000 llround(x): 0 Input value: 5.00000000000000000000e001 round(x): 0.00000000000000000000e+000 llround(x): 0 Input value: 5.00000000000000000000e001 round(x): 0.00000000000000000000e+000 llround(x): 0 Input value: 5.00000000000000110000e001 round(x): 1.00000000000000000000e+000 llround(x): 1 Input value: 5.00000000000000110000e001 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.00000000000000000000e001) == 0.0, using the roundtoeven rule for 0.5 Danny  Comment By: Greg Chicares (chicares) Date: 20080516 00:00 Message: Logged In: YES user_id=231935 Originator: YES I have no objection to anything in your 20080514 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: 20080514 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: altround.tar.gz  Comment By: Keith Marshall (keithmarshall) Date: 20080513 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 