summaryrefslogtreecommitdiff
path: root/libm
diff options
context:
space:
mode:
Diffstat (limited to 'libm')
-rw-r--r--libm/double/floor.c78
-rw-r--r--libm/double/rint.c52
-rw-r--r--libm/double/round.c69
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;