math: rewrite rounding functions (ceil, floor, trunc, round, rint)

* faster, smaller, cleaner implementation than the bit hacks of fdlibm
* use arithmetics like y=(double)(x+0x1p52)-0x1p52, which is an integer
neighbor of x in all rounding modes (0<=x<0x1p52) and only use bithacks
when that's faster and smaller (for float it usually is)
* the code assumes standard excess precision handling for casts
* long double code supports both ld80 and ld128
* nearbyint is not changed (it is a wrapper around rint)

1 files changed, 15 insertions, 43 deletions

diff --git a/src/math/rintf.c b/src/math/rintf.c
index e8d44969..9cfc2a26 100644
--- a/src/math/rintf.c
+++ b/src/math/rintf.c
@@ -1,48 +1,20 @@
-/* origin: FreeBSD /usr/src/lib/msun/src/s_rintf.c */
-/*
- * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
- */
-/*
- * ====================================================
- * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
- *
- * Developed at SunPro, a Sun Microsystems, Inc. business.
- * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice
- * is preserved.
- * ====================================================
- */
-
-#include "libm.h"
-
-static const float
-TWO23[2] = {
-  8.3886080000e+06, /* 0x4b000000 */
- -8.3886080000e+06, /* 0xcb000000 */
-};
+#include <math.h>
+#include <stdint.h>
 
 float rintf(float x)
 {
-	int32_t i0,j0,sx;
-	float w,t;
+	union {float f; uint32_t i;} u = {x};
+	int e = u.i>>23 & 0xff;
+	int s = u.i>>31;
+	float_t y;
 
-	GET_FLOAT_WORD(i0, x);
-	sx = (i0>>31) & 1;
-	j0 = ((i0>>23)&0xff) - 0x7f;
-	if (j0 < 23) {
-		if (j0 < 0) {
-			if ((i0&0x7fffffff) == 0)
-				return x;
-			STRICT_ASSIGN(float, w, TWO23[sx] + x);
-			t = w - TWO23[sx];
-			GET_FLOAT_WORD(i0, t);
-			SET_FLOAT_WORD(t, (i0&0x7fffffff)|(sx<<31));
-			return t;
-		}
-		STRICT_ASSIGN(float, w, TWO23[sx] + x);
-		return w - TWO23[sx];
-	}
-	if (j0 == 0x80)
-		return x+x;  /* inf or NaN */
-	return x;            /* x is integral */
+	if (e >= 0x7f+23)
+		return x;
+	if (s)
+		y = (float)(x - 0x1p23f) + 0x1p23f;
+	else
+		y = (float)(x + 0x1p23f) - 0x1p23f;
+	if (y == 0)
+		return s ? -0.0f : 0.0f;
+	return y;
 }