Remove Erik's broken implementation of rint(). Replace it by one "less broken".

Also correct rounding beharior of round() and add trunc().  Note that round()
and rint() currently don't check for infs and nans.  I decided to wait on
that until the big cleanup.

3 files changed, 117 insertions, 82 deletions

diff --git a/libm/double/floor.c b/libm/double/floor.c
index dcc1a10f1..affc7753e 100644
--- a/libm/double/floor.c
+++ b/libm/double/floor.c
@@ -451,3 +451,81 @@ else
 	return(u.y);
 	}
 }
+
+/**********************************************************************/
+/*
+ * trunc is just a slightly modified version of floor above.
+ */
+
+double trunc(double x)
+{
+	union {
+		double y;
+		unsigned short sh[4];
+	} u;
+	unsigned short *p;
+	int e;
+
+#ifdef UNK
+	mtherr( "trunc", DOMAIN );
+	return(0.0);
+#endif
+#ifdef NANS
+	if( isnan(x) )
+		return( x );
+#endif
+#ifdef INFINITIES
+	if(!isfinite(x))
+		return(x);
+#endif
+#ifdef MINUSZERO
+	if(x == 0.0L)
+		return(x);
+#endif
+	u.y = x;
+	/* find the exponent (power of 2) */
+#ifdef DEC
+	p = (unsigned short *)&u.sh[0];
+	e = (( *p  >> 7) & 0377) - 0201;
+	p += 3;
+#endif
+
+#ifdef IBMPC
+	p = (unsigned short *)&u.sh[3];
+	e = (( *p >> 4) & 0x7ff) - 0x3ff;
+	p -= 3;
+#endif
+
+#ifdef MIEEE
+	p = (unsigned short *)&u.sh[0];
+	e = (( *p >> 4) & 0x7ff) - 0x3ff;
+	p += 3;
+#endif
+
+	if( e < 0 )
+		return( 0.0 );
+
+	e = (NBITS -1) - e;
+	/* clean out 16 bits at a time */
+	while( e >= 16 )
+		{
+#ifdef IBMPC
+			*p++ = 0;
+#endif
+
+#ifdef DEC
+			*p-- = 0;
+#endif
+
+#ifdef MIEEE
+			*p-- = 0;
+#endif
+			e -= 16;
+		}
+
+	/* clear the remaining bits */
+	if( e > 0 )
+		*p &= bmask[e];
+
+	return(u.y);
+}
diff --git a/libm/double/rint.c b/libm/double/rint.c
deleted file mode 100644
index 35cf5f503..000000000
--- a/libm/double/rint.c
+++ /dev/null
@@ -1,52 +0,0 @@
-/* vi: set sw=4 ts=4: */
-/*
- * rint for uClibc
- *
- * Copyright (C) 2001 by Lineo, inc.  
- * Written by Erik Andersen <andersen@lineo.com>, <andersee@debian.org>
- *
- * This program is free software; you can redistribute it and/or modify it
- * under the terms of the GNU Library General Public License as published by
- * the Free Software Foundation; either version 2 of the License, or (at your
- * option) any later version.
- *
- * This program is distributed in the hope that it will be useful, but WITHOUT
- * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
- * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public License
- * for more details.
- *
- * You should have received a copy of the GNU Library General Public License
- * along with this program; if not, write to the Free Software Foundation,
- * Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
- *
- */
-
-#include <math.h>
-
-/* From the Linux man page:
- *
- * NAME
- *	rint - round to closest integer
- *
- * SYNOPSIS
- *	#include <math.h>
- *	double rint(double x);
- *
- * DESCRIPTION
- *	The rint() function rounds x to an integer value according
- *	to the prevalent rounding mode.  The default rounding mode
- *	is to round to the nearest integer.  
- *
- * RETURN VALUE
- *	The rint() function returns the integer value as a float­ 
- *	ing-point number.
- */
-
-double rint (double x) {
-	double low = floor(x);
-	if (fmod(x,low) >= (double)0.5)
-		return(ceil(x));
-	else
-		return(low);
-}
-
diff --git a/libm/double/round.c b/libm/double/round.c
index df4564a0f..c80ae66a0 100644
--- a/libm/double/round.c
+++ b/libm/double/round.c
@@ -1,45 +1,54 @@
-/*							round.c
- *
- *	Round double to nearest or even integer valued double
- *
- *
- *
- * SYNOPSIS:
- *
- * double x, y, round();
- *
- * y = round(x);
- *
+/*
+ * June 19, 2001       Manuel Novoa III
  *
+ * Replaced cephes round (which was actually round to nearest or even)
+ * with a (really lame actually) version that always rounds away from 0
+ * in conformance with ANSI/ISO.
  *
- * DESCRIPTION:
+ * This doesn't check for inf or nan (hence the lame part) but the
+ * cephes function it replaces didn't either.  I plan to deal with
+ * those issues when I rework things w.r.t. common code.
  *
+ * Also, for now rename the original cephes round routine to rint since
+ * it behaves the same for the default rounding mode (round to nearest).
+ * This will have to be changed off course when floating point env
+ * control functions are added.
+ */
+
+#include <math.h>
+
+double round(x)
+double x;
+{
+	double ax, fax;
+
+	ax = fabs(x);
+	fax = floor(ax);
+	if (ax - fax >= 0.5) {
+		fax += 1.0;
+	}
+	if (x < 0) {
+		x = -fax;
+	} else {
+		x = fax;
+	}
+	return x;
+}
+
+/***********************************************************************/
+/*
  * Returns the nearest integer to x as a double precision
  * floating point result.  If x ends in 0.5 exactly, the
  * nearest even integer is chosen.
- * 
- *
- *
- * ACCURACY:
- *
- * If x is greater than 1/(2*MACHEP), its closest machine
- * representation is already an integer, so rounding does
- * not change it.
- */
-
+ * /
 /*
+Originally round from
 Cephes Math Library Release 2.1:  January, 1989
 Copyright 1984, 1987, 1989 by Stephen L. Moshier
 Direct inquiries to 30 Frost Street, Cambridge, MA 02140
 */
-#include <math.h>
-#ifdef ANSIPROT
-double floor ( double );
-#else
-double floor();
-#endif
 
-double round(x)
+double rint(x)
 double x;
 {
 double y, r;