123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157 |
- /*******************************************************************************
- ** 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 );
- }
|