Fixed roundf(). The following cases never worked in FreeBSD:

- in round-towards-minus-infinity mode, on all machines, roundf(x) never
  worked for 0 < |x| < 0.5 (2*0x3effffff cases in all, or almost half of
  float space).  It was -0 for 0 < x < 0.5 and 0 for -0.5 < x < 0, but
  should be 0 and -0, respectively.  This is because t = ceilf(|x|) = 1
  for these args, and when we adjust t from 1 to 0 by subtracting 1, we
  get -0 in this rounding mode, but we want and expected to get 0.
- in round-towards-minus-infinity, round towards zero and round-to-nearest
  modes, on machines that evaluate float expressions in float precision
  (most machines except i386's), roundf(x) never worked for |x| =
  <float value immediately below 0.5> (2 cases in all).  It was +-1 but
  should have been +-0.  This is because t = ceilf(|x|) = 1 for these
  args, and when we try to classify |x| by subtracting it from 1 we
  get an unexpected rounding error -- the result is 0.5 after rounding
  to float in all 3 rounding modes, so we we have forgotten the
  difference between |x| and 0.5 and end up returning the same value
  as for +-0.5.

The fix is to use floorf() instead of ceilf() and to add 1 instead of
-1 in the adjustment.  With floorf() all the expressions used are
always evaluated exactly so there are no rounding problems, and with
adjustments of +1 we don't go near -0 when adjusting.

Attempted to fix round() and roundl() by cloning the fix for roundf().
This has only been tested for round(), only on args representable as
floats.  Double expressions are evaluated in double precision even on
i386's, so round(0.5-epsilon) was broken even on i386's.  roundl()
must be completely broken on i386's since long double precision is not
really supported.  There seem to be no other dependencies on the
precision.

3 files changed, 18 insertions, 18 deletions

diff --git a/lib/msun/src/s_round.c b/lib/msun/src/s_round.c
index 34b26f0..65de31b 100644
--- a/lib/msun/src/s_round.c
+++ b/lib/msun/src/s_round.c
@@ -38,14 +38,14 @@ round(double x)
 		return (x);
 
 	if (x >= 0.0) {
-		t = ceil(x);
-		if (t - x > 0.5)
-			t -= 1.0;
+		t = floor(x);
+		if (t - x <= -0.5)
+			t += 1.0;
 		return (t);
 	} else {
-		t = ceil(-x);
-		if (t + x > 0.5)
-			t -= 1.0;
+		t = floor(-x);
+		if (t + x <= -0.5)
+			t += 1.0;
 		return (-t);
 	}
 }
diff --git a/lib/msun/src/s_roundf.c b/lib/msun/src/s_roundf.c
index a7a3da9..952e8e7 100644
--- a/lib/msun/src/s_roundf.c
+++ b/lib/msun/src/s_roundf.c
@@ -38,14 +38,14 @@ roundf(float x)
 		return (x);
 
 	if (x >= 0.0) {
-		t = ceilf(x);
-		if (t - x > 0.5)
-			t -= 1.0;
+		t = floorf(x);
+		if (t - x <= -0.5)
+			t += 1.0;
 		return (t);
 	} else {
-		t = ceilf(-x);
-		if (t + x > 0.5)
-			t -= 1.0;
+		t = floorf(-x);
+		if (t + x <= -0.5)
+			t += 1.0;
 		return (-t);
 	}
 }
diff --git a/lib/msun/src/s_roundl.c b/lib/msun/src/s_roundl.c
index 0cb688b..a70b617 100644
--- a/lib/msun/src/s_roundl.c
+++ b/lib/msun/src/s_roundl.c
@@ -38,14 +38,14 @@ roundl(long double x)
 		return (x);
 
 	if (x >= 0.0) {
-		t = ceill(x);
-		if (t - x > 0.5)
-			t -= 1.0;
+		t = floorl(x);
+		if (t - x <= -0.5)
+			t += 1.0;
 		return (t);
 	} else {
-		t = ceill(-x);
-		if (t + x > 0.5)
-			t -= 1.0;
+		t = floorl(-x);
+		if (t + x <= -0.5)
+			t += 1.0;
 		return (-t);
 	}
 }