diff options
Diffstat (limited to 'libm/powerpc/s_rint.c')
| -rw-r--r-- | libm/powerpc/s_rint.c | 157 |
1 files changed, 157 insertions, 0 deletions
diff --git a/libm/powerpc/s_rint.c b/libm/powerpc/s_rint.c new file mode 100644 index 000000000..47a79ae5e --- /dev/null +++ b/libm/powerpc/s_rint.c @@ -0,0 +1,157 @@ +/******************************************************************************* +** File: rndint.c +** +** Contains: C source code for implementations of floating-point +** functions which round to integral value or format, as +** defined in header <fp.h>. In particular, this file +** contains implementations of functions rint, nearbyint, +** rinttol, round, roundtol, trunc, modf and modfl. This file +** targets PowerPC or Power platforms. +** +** Written by: A. Sazegari, Apple AltiVec Group +** Created originally by Jon Okada, Apple Numerics Group +** +** Copyright: © 1992-2001 by Apple Computer, Inc., all rights reserved +** +** Change History (most recent first): +** +** 13 Jul 01 ram replaced --setflm calls with inline assembly +** 03 Mar 01 ali first port to os x using gcc, added the crucial __setflm +** definition. +** 1. removed double_t, put in double for now. +** 2. removed iclass from nearbyint. +** 3. removed wrong comments intrunc. +** 4. +** 13 May 97 ali made performance improvements in rint, rinttol, roundtol +** and trunc by folding some of the taligent ideas into this +** implementation. nearbyint is faster than the one in taligent, +** rint is more elegant, but slower by %30 than the taligent one. +** 09 Apr 97 ali deleted modfl and deferred to AuxiliaryDD.c +** 15 Sep 94 ali Major overhaul and performance improvements of all functions. +** 20 Jul 94 PAF New faster version +** 16 Jul 93 ali Added the modfl function. +** 18 Feb 93 ali Changed the return value of fenv functions +** feclearexcept and feraiseexcept to their new +** NCEG X3J11.1/93-001 definitions. +** 16 Dec 92 JPO Removed __itrunc implementation to a +** separate file. +** 15 Dec 92 JPO Added __itrunc implementation and modified +** rinttol to include conversion from double +** to long int format. Modified roundtol to +** call __itrunc. +** 10 Dec 92 JPO Added modf (double) implementation. +** 04 Dec 92 JPO First created. +** +*******************************************************************************/ + +#include <limits.h> +#include <math.h> + +#define SET_INVALID 0x01000000UL + +typedef union + { + struct { +#if defined(__BIG_ENDIAN__) + unsigned long int hi; + unsigned long int lo; +#else + unsigned long int lo; + unsigned long int hi; +#endif + } words; + double dbl; + } DblInHex; + +static const unsigned long int signMask = 0x80000000ul; +static const double twoTo52 = 4503599627370496.0; +static const double doubleToLong = 4503603922337792.0; // 2^52 +static const DblInHex Huge = {{ 0x7FF00000, 0x00000000 }}; +static const DblInHex TOWARDZERO = {{ 0x00000000, 0x00000001 }}; + +/******************************************************************************* +* * +* The function rint rounds its double argument to integral value * +* according to the current rounding direction and returns the result in * +* double format. This function signals inexact if an ordered return * +* value is not equal to the operand. * +* * +******************************************************************************** +* * +* This function calls: fabs. * +* * +*******************************************************************************/ + +/******************************************************************************* +* First, an elegant implementation. * +******************************************************************************** +* +*double rint ( double x ) +* { +* double y; +* +* y = twoTo52.fval; +* +* if ( fabs ( x ) >= y ) // huge case is exact +* return x; +* if ( x < 0 ) y = -y; // negative case +* y = ( x + y ) - y; // force rounding +* if ( y == 0.0 ) // zero results mirror sign of x +* y = copysign ( y, x ); +* return ( y ); +* } +******************************************************************************** +* Now a bit twidling version that is about %30 faster. * +*******************************************************************************/ + +double rint ( double x ) + { + DblInHex argument; + register double y; + unsigned long int xHead; + register long int target; + + argument.dbl = x; + xHead = argument.words.hi & 0x7fffffffUL; // xHead <- high half of |x| + target = ( argument.words.hi < signMask ); // flags positive sign + + if ( xHead < 0x43300000ul ) +/******************************************************************************* +* Is |x| < 2.0^52? * +*******************************************************************************/ + { + if ( xHead < 0x3ff00000ul ) +/******************************************************************************* +* Is |x| < 1.0? * +*******************************************************************************/ + { + if ( target ) + y = ( x + twoTo52 ) - twoTo52; // round at binary point + else + y = ( x - twoTo52 ) + twoTo52; // round at binary point + if ( y == 0.0 ) + { // fix sign of zero result + if ( target ) + return ( 0.0 ); + else + return ( -0.0 ); + } + return y; + } + +/******************************************************************************* +* Is 1.0 < |x| < 2.0^52? * +*******************************************************************************/ + + if ( target ) + return ( ( x + twoTo52 ) - twoTo52 ); // round at binary pt. + else + return ( ( x - twoTo52 ) + twoTo52 ); + } + +/******************************************************************************* +* |x| >= 2.0^52 or x is a NaN. * +*******************************************************************************/ + return ( x ); + } + |