diff options
Diffstat (limited to 'libm')
| -rw-r--r-- | libm/double/floor.c | 78 | ||||
| -rw-r--r-- | libm/double/rint.c | 52 | ||||
| -rw-r--r-- | libm/double/round.c | 69 |
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; |