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