Hello Ben, > I'm appending a program that can be used to test the correctness > of each of the three implementations, by evaluating them for > every possible 32-bit value of float and checking the results > against the system roundf's result.
Well done! This test program can be made faster, by varying only the highest and the lowest bits; the bits in the middle can be set to 00....00 or 11...11 without weakening the check. This allows to generalize the test to 'double' (64-bit). When applying this test program to 'floor', 'ceil', 'trunc', it uncovered bugs there. I've fixed them now: http://git.sv.gnu.org/gitweb/?p=gnulib.git;a=commitdiff;h=3db9b0d802e34ec28fc924076ae4c196bd7d23a9 http://git.sv.gnu.org/gitweb/?p=gnulib.git;a=commitdiff;h=33d6a0e5ca06c6d0a8b0fddba8c7f41bbd6794ae http://git.sv.gnu.org/gitweb/?p=gnulib.git;a=commitdiff;h=c7f33ef75bc531bc875e34019324811962b001d1 > I took the liberty of > correcting a bug that had crept into round_bruno. Thanks, this was a mistake. > (By the way, is it intentional that TWO_MANT_DIG is always of type > "double"? Should it be type DOUBLE instead?) A mistake as well. > round_bruno does not: it fails in many cases, such as > 0.5-2*epsilon: > round 0.500000(0x1.fffffep-2) = 0.000000(0x0p+0) or 1.000000(0x1p+0)? Here's a better version of both. ============================================================================== DOUBLE round_bruno1 (DOUBLE x) { /* The use of 'volatile' guarantees that excess precision bits are dropped at each addition step and before the following comparison at the caller's site. It is necessary on x86 systems where double-floats are not IEEE compliant by default, to avoid that the results become platform and compiler option dependent. 'volatile' is a portable alternative to gcc's -ffloat-store option. */ volatile DOUBLE y = x; volatile DOUBLE z = y; if (z > L_(0.0)) { /* Avoid rounding error for x = 0.5 - 2^(-MANT_DIG-1). */ if (z < L_(0.5)) z = L_(0.0); /* Avoid rounding errors for values near 2^k, where k >= MANT_DIG-1. */ else if (z < TWO_MANT_DIG) { /* Add 0.5 to the absolute value. */ y = z += L_(0.5); /* Round to the next integer (nearest or up or down, doesn't matter). */ z += TWO_MANT_DIG; z -= TWO_MANT_DIG; /* Enforce rounding down. */ if (z > y) z -= L_(1.0); } } else if (z < L_(0.0)) { /* Avoid rounding error for x = -(0.5 - 2^(-MANT_DIG-1)). */ if (z > - L_(0.5)) z = L_(0.0); /* Avoid rounding errors for values near -2^k, where k >= MANT_DIG-1. */ else if (z > -TWO_MANT_DIG) { /* Add 0.5 to the absolute value. */ y = z -= L_(0.5); /* Round to the next integer (nearest or up or down, doesn't matter). */ z -= TWO_MANT_DIG; z += TWO_MANT_DIG; /* Enforce rounding up. */ if (z < y) z += L_(1.0); } } return z; } DOUBLE round_bruno2 (DOUBLE x) { if (x >= L_(0.0)) { /* Avoid rounding error for x = 0.5 - 2^(-MANT_DIG-1). */ if (x < L_(0.5)) return L_(0.0); if (x < TWO_MANT_DIG) return FLOOR (x + L_(0.5)); } else { /* Avoid rounding error for x = -(0.5 - 2^(-MANT_DIG-1)). */ if (x > - L_(0.5)) return L_(0.0); if (x > - TWO_MANT_DIG) return CEIL (x - L_(0.5)); } return x; } ============================================================================== What to do with these two functions? I would suggest to 1) test them both as reference implementations, like done in test-truncf2.c and test-trunc2.c. 2) use round_bruno1 as the implementation in lib/round.c if the corresponding functions floor, ceil (or floorf, ceilf, or floorl, ceill) are not declared. If they are declared, one can assume that they are efficiently implemented, and then round_blp looks like quite efficient. Bruno
