summaryrefslogtreecommitdiff
path: root/libm/rndint.c
diff options
context:
space:
mode:
authorEric Andersen <andersen@codepoet.org>2001-11-24 03:46:25 +0000
committerEric Andersen <andersen@codepoet.org>2001-11-24 03:46:25 +0000
commit26d7ea91124a9405dff7755a78cfb6232dd15d80 (patch)
treec14436f6013bd76a2a339deb1424e02e73419957 /libm/rndint.c
parent683c13fcc85276e9a030d6a98d50366bef03a6b6 (diff)
Move powerpc specific optimizations (courtesy of apple) to powerpc
subdir. Put together a theoretical framework for adding arch specific optimizations. Havn't tried this yet but it looks correct... -Erik
Diffstat (limited to 'libm/rndint.c')
-rw-r--r--libm/rndint.c632
1 files changed, 0 insertions, 632 deletions
diff --git a/libm/rndint.c b/libm/rndint.c
deleted file mode 100644
index 7f8c183d4..000000000
--- a/libm/rndint.c
+++ /dev/null
@@ -1,632 +0,0 @@
-/*******************************************************************************
-** 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>
-
-#if !defined(__ppc__)
-#define asm(x)
-#endif
-
-#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. *
-*******************************************************************************/
-
-#if defined(__ppc__)
-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 );
- }
-#endif /* __ppc__ */
-
-/*******************************************************************************
-* *
-* The function nearbyint rounds its double argument to integral value *
-* according to the current rounding direction and returns the result in *
-* double format. This function does not signal inexact. *
-* *
-********************************************************************************
-* *
-* This function calls fabs and copysign. *
-* *
-*******************************************************************************/
-
-double nearbyint ( double x )
- {
- double y;
-#if defined(__ppc__)
- double OldEnvironment;
-#endif /* __ppc__ */
-
- y = twoTo52;
-
- asm ("mffs %0" : "=f" (OldEnvironment)); /* get the environement */
-
- 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 );
-// restore old flags
- asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment ));
- return ( y );
- }
-
-/*******************************************************************************
-* *
-* The function rinttol converts its double argument to integral value *
-* according to the current rounding direction and returns the result in *
-* long int format. This conversion signals invalid if the argument is a *
-* NaN or the rounded intermediate result is out of range of the *
-* destination long int format, and it delivers an unspecified result in *
-* this case. This function signals inexact if the rounded result is *
-* within range of the long int format but unequal to the operand. *
-* *
-*******************************************************************************/
-
-long int rinttol ( double x )
- {
- register double y;
- DblInHex argument, OldEnvironment;
- unsigned long int xHead;
- register long int target;
-
- argument.dbl = x;
- target = ( argument.words.hi < signMask ); // flag positive sign
- xHead = argument.words.hi & 0x7ffffffful; // high 32 bits of x
-
- if ( target )
-/*******************************************************************************
-* Sign of x is positive. *
-*******************************************************************************/
- {
- if ( xHead < 0x41dffffful )
- { // x is safely in long range
- y = ( x + twoTo52 ) - twoTo52; // round at binary point
- argument.dbl = y + doubleToLong; // force result into argument.words.lo
- return ( ( long ) argument.words.lo );
- }
-
- asm ("mffs %0" : "=f" (OldEnvironment.dbl)); // get environment
-
- if ( xHead > 0x41dffffful )
- { // x is safely out of long range
- OldEnvironment.words.lo |= SET_INVALID;
- asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl ));
- return ( LONG_MAX );
- }
-
-/*******************************************************************************
-* x > 0.0 and may or may not be out of range of long. *
-*******************************************************************************/
-
- y = ( x + twoTo52 ) - twoTo52; // do rounding
- if ( y > ( double ) LONG_MAX )
- { // out of range of long
- OldEnvironment.words.lo |= SET_INVALID;
- asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl ));
- return ( LONG_MAX );
- }
- argument.dbl = y + doubleToLong; // in range
- return ( ( long ) argument.words.lo ); // return result & flags
- }
-
-/*******************************************************************************
-* Sign of x is negative. *
-*******************************************************************************/
- if ( xHead < 0x41e00000ul )
- { // x is safely in long range
- y = ( x - twoTo52 ) + twoTo52;
- argument.dbl = y + doubleToLong;
- return ( ( long ) argument.words.lo );
- }
-
- asm ("mffs %0" : "=f" (OldEnvironment.dbl)); // get environment
-
- if ( xHead > 0x41e00000ul )
- { // x is safely out of long range
- OldEnvironment.words.lo |= SET_INVALID;
- asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl ));
- return ( LONG_MIN );
- }
-
-/*******************************************************************************
-* x < 0.0 and may or may not be out of range of long. *
-*******************************************************************************/
-
- y = ( x - twoTo52 ) + twoTo52; // do rounding
- if ( y < ( double ) LONG_MIN )
- { // out of range of long
- OldEnvironment.words.lo |= SET_INVALID;
- asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl ));
- return ( LONG_MIN );
- }
- argument.dbl = y + doubleToLong; // in range
- return ( ( long ) argument.words.lo ); // return result & flags
- }
-
-/*******************************************************************************
-* *
-* The function round rounds its double argument to integral value *
-* according to the "add half to the magnitude and truncate" rounding of *
-* Pascal's Round function and FORTRAN's ANINT function and returns the *
-* result in double format. This function signals inexact if an ordered *
-* return value is not equal to the operand. *
-* *
-*******************************************************************************/
-
-double round ( double x )
- {
- DblInHex argument, OldEnvironment;
- register double y, z;
- register 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 ); // flag positive sign
-
- if ( xHead < 0x43300000ul )
-/*******************************************************************************
-* Is |x| < 2.0^52? *
-*******************************************************************************/
- {
- if ( xHead < 0x3ff00000ul )
-/*******************************************************************************
-* Is |x| < 1.0? *
-*******************************************************************************/
- {
- asm ("mffs %0" : "=f" (OldEnvironment.dbl)); // get environment
- if ( xHead < 0x3fe00000ul )
-/*******************************************************************************
-* Is |x| < 0.5? *
-*******************************************************************************/
- {
- if ( ( xHead | argument.words.lo ) != 0ul )
- OldEnvironment.words.lo |= 0x02000000ul;
- asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl ));
- if ( target )
- return ( 0.0 );
- else
- return ( -0.0 );
- }
-/*******************************************************************************
-* Is 0.5 ² |x| < 1.0? *
-*******************************************************************************/
- OldEnvironment.words.lo |= 0x02000000ul;
- asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl ));
- if ( target )
- return ( 1.0 );
- else
- return ( -1.0 );
- }
-/*******************************************************************************
-* Is 1.0 < |x| < 2.0^52? *
-*******************************************************************************/
- if ( target )
- { // positive x
- y = ( x + twoTo52 ) - twoTo52; // round at binary point
- if ( y == x ) // exact case
- return ( x );
- z = x + 0.5; // inexact case
- y = ( z + twoTo52 ) - twoTo52; // round at binary point
- if ( y > z )
- return ( y - 1.0 );
- else
- return ( y );
- }
-
-/*******************************************************************************
-* Is x < 0? *
-*******************************************************************************/
- else
- {
- y = ( x - twoTo52 ) + twoTo52; // round at binary point
- if ( y == x )
- return ( x );
- z = x - 0.5;
- y = ( z - twoTo52 ) + twoTo52; // round at binary point
- if ( y < z )
- return ( y + 1.0 );
- else
- return ( y );
- }
- }
-/*******************************************************************************
-* |x| >= 2.0^52 or x is a NaN. *
-*******************************************************************************/
- return ( x );
- }
-
-/*******************************************************************************
-* *
-* The function roundtol converts its double argument to integral format *
-* according to the "add half to the magnitude and chop" rounding mode of *
-* Pascal's Round function and FORTRAN's NINT function. This conversion *
-* signals invalid if the argument is a NaN or the rounded intermediate *
-* result is out of range of the destination long int format, and it *
-* delivers an unspecified result in this case. This function signals *
-* inexact if the rounded result is within range of the long int format but *
-* unequal to the operand. *
-* *
-*******************************************************************************/
-
-long int roundtol ( double x )
- {
- register double y, z;
- DblInHex argument, OldEnvironment;
- register unsigned long int xhi;
- register long int target;
-#if defined(__ppc__)
- const DblInHex kTZ = {{ 0x0, 0x1 }};
- const DblInHex kUP = {{ 0x0, 0x2 }};
-#endif /* __ppc__ */
-
- argument.dbl = x;
- xhi = argument.words.hi & 0x7ffffffful; // high 32 bits of x
- target = ( argument.words.hi < signMask ); // flag positive sign
-
- if ( xhi > 0x41e00000ul )
-/*******************************************************************************
-* Is x is out of long range or NaN? *
-*******************************************************************************/
- {
- asm ("mffs %0" : "=f" (OldEnvironment.dbl)); // get environment
- OldEnvironment.words.lo |= SET_INVALID;
- asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl ));
- if ( target ) // pin result
- return ( LONG_MAX );
- else
- return ( LONG_MIN );
- }
-
- if ( target )
-/*******************************************************************************
-* Is sign of x is "+"? *
-*******************************************************************************/
- {
- if ( x < 2147483647.5 )
-/*******************************************************************************
-* x is in the range of a long. *
-*******************************************************************************/
- {
- y = ( x + doubleToLong ) - doubleToLong; // round at binary point
- if ( y != x )
- { // inexact case
- asm ("mffs %0" : "=f" (OldEnvironment.dbl)); // save environment
- asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( kTZ.dbl )); // truncate rounding
- z = x + 0.5; // truncate x + 0.5
- argument.dbl = z + doubleToLong;
- asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl ));
- return ( ( long ) argument.words.lo );
- }
-
- argument.dbl = y + doubleToLong; // force result into argument.words.lo
- return ( ( long ) argument.words.lo ); // return long result
- }
-/*******************************************************************************
-* Rounded positive x is out of the range of a long. *
-*******************************************************************************/
- asm ("mffs %0" : "=f" (OldEnvironment.dbl));
- OldEnvironment.words.lo |= SET_INVALID;
- asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl ));
- return ( LONG_MAX ); // return pinned result
- }
-/*******************************************************************************
-* x < 0.0 and may or may not be out of the range of a long. *
-*******************************************************************************/
- if ( x > -2147483648.5 )
-/*******************************************************************************
-* x is in the range of a long. *
-*******************************************************************************/
- {
- y = ( x + doubleToLong ) - doubleToLong; // round at binary point
- if ( y != x )
- { // inexact case
- asm ("mffs %0" : "=f" (OldEnvironment.dbl)); // save environment
- asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( kUP.dbl )); // round up
- z = x - 0.5; // truncate x - 0.5
- argument.dbl = z + doubleToLong;
- asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl ));
- return ( ( long ) argument.words.lo );
- }
-
- argument.dbl = y + doubleToLong;
- return ( ( long ) argument.words.lo ); // return long result
- }
-/*******************************************************************************
-* Rounded negative x is out of the range of a long. *
-*******************************************************************************/
- asm ("mffs %0" : "=f" (OldEnvironment.dbl));
- OldEnvironment.words.lo |= SET_INVALID;
- asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl ));
- return ( LONG_MIN ); // return pinned result
- }
-
-/*******************************************************************************
-* *
-* The function trunc truncates its double argument to integral value *
-* and returns the result in double format. This function signals *
-* inexact if an ordered return value is not equal to the operand. *
-* *
-*******************************************************************************/
-
-double trunc ( double x )
- {
- DblInHex argument,OldEnvironment;
- register double y;
- register unsigned long int xhi;
- register long int target;
-
- argument.dbl = x;
- xhi = argument.words.hi & 0x7fffffffUL; // xhi <- high half of |x|
- target = ( argument.words.hi < signMask ); // flag positive sign
-
- if ( xhi < 0x43300000ul )
-/*******************************************************************************
-* Is |x| < 2.0^53? *
-*******************************************************************************/
- {
- if ( xhi < 0x3ff00000ul )
-/*******************************************************************************
-* Is |x| < 1.0? *
-*******************************************************************************/
- {
- if ( ( xhi | argument.words.lo ) != 0ul )
- { // raise deserved INEXACT
- asm ("mffs %0" : "=f" (OldEnvironment.dbl));
- OldEnvironment.words.lo |= 0x02000000ul;
- asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl ));
- }
- if ( target ) // return properly signed zero
- return ( 0.0 );
- else
- return ( -0.0 );
- }
-/*******************************************************************************
-* Is 1.0 < |x| < 2.0^52? *
-*******************************************************************************/
- if ( target )
- {
- y = ( x + twoTo52 ) - twoTo52; // round at binary point
- if ( y > x )
- return ( y - 1.0 );
- else
- return ( y );
- }
-
- else
- {
- y = ( x - twoTo52 ) + twoTo52; // round at binary point.
- if ( y < x )
- return ( y + 1.0 );
- else
- return ( y );
- }
- }
-/*******************************************************************************
-* Is |x| >= 2.0^52 or x is a NaN. *
-*******************************************************************************/
- return ( x );
- }
-
-/*******************************************************************************
-* The modf family of functions separate a floating-point number into its *
-* fractional and integral parts, returning the fractional part and writing *
-* the integral part in floating-point format to the object pointed to by a *
-* pointer argument. If the input argument is integral or infinite in *
-* value, the return value is a zero with the sign of the input argument. *
-* The modf family of functions raises no floating-point exceptions. older *
-* implemenation set the INVALID flag due to signaling NaN input. *
-* *
-*******************************************************************************/
-
-/*******************************************************************************
-* modf is the double implementation. *
-*******************************************************************************/
-
-#if defined(__ppc__)
-double modf ( double x, double *iptr )
- {
- register double OldEnvironment, xtrunc;
- register unsigned long int xHead, signBit;
- DblInHex argument;
-
- argument.dbl = x;
- xHead = argument.words.hi & 0x7ffffffful; // |x| high bit pattern
- signBit = ( argument.words.hi & 0x80000000ul ); // isolate sign bit
- if (xHead == 0x7ff81fe0)
- signBit = signBit | 0;
-
- if ( xHead < 0x43300000ul )
-/*******************************************************************************
-* Is |x| < 2.0^53? *
-*******************************************************************************/
- {
- if ( xHead < 0x3ff00000ul )
-/*******************************************************************************
-* Is |x| < 1.0? *
-*******************************************************************************/
- {
- argument.words.hi = signBit; // truncate to zero
- argument.words.lo = 0ul;
- *iptr = argument.dbl;
- return ( x );
- }
-/*******************************************************************************
-* Is 1.0 < |x| < 2.0^52? *
-*******************************************************************************/
- asm ("mffs %0" : "=f" (OldEnvironment)); // save environment
- // round toward zero
- asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( TOWARDZERO.dbl ));
- if ( signBit == 0ul ) // truncate to integer
- xtrunc = ( x + twoTo52 ) - twoTo52;
- else
- xtrunc = ( x - twoTo52 ) + twoTo52;
- // restore caller's env
- asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment ));
- *iptr = xtrunc; // store integral part
- if ( x != xtrunc ) // nonzero fraction
- return ( x - xtrunc );
- else
- { // zero with x's sign
- argument.words.hi = signBit;
- argument.words.lo = 0ul;
- return ( argument.dbl );
- }
- }
-
- *iptr = x; // x is integral or NaN
- if ( x != x ) // NaN is returned
- return x;
- else
- { // zero with x's sign
- argument.words.hi = signBit;
- argument.words.lo = 0ul;
- return ( argument.dbl );
- }
- }
-#endif /* __ppc__ */