diff options
Diffstat (limited to 'libm/powerpc')
| -rw-r--r-- | libm/powerpc/classic/s_ceil.c | 113 | ||||
| -rw-r--r-- | libm/powerpc/classic/s_copysign.c | 59 | ||||
| -rw-r--r-- | libm/powerpc/classic/s_floor.c | 113 | ||||
| -rw-r--r-- | libm/powerpc/classic/s_frexp.c | 65 | ||||
| -rw-r--r-- | libm/powerpc/classic/s_ldexp.c | 49 | ||||
| -rw-r--r-- | libm/powerpc/classic/s_logb.c | 107 | ||||
| -rw-r--r-- | libm/powerpc/classic/s_modf.c | 344 | ||||
| -rw-r--r-- | libm/powerpc/classic/s_nearbyint.c | 38 | ||||
| -rw-r--r-- | libm/powerpc/classic/s_rint.c | 159 | ||||
| -rw-r--r-- | libm/powerpc/classic/s_round.c | 115 | ||||
| -rw-r--r-- | libm/powerpc/classic/s_trunc.c | 89 | ||||
| -rw-r--r-- | libm/powerpc/classic/w_scalb.c | 94 |
12 files changed, 0 insertions, 1345 deletions
diff --git a/libm/powerpc/classic/s_ceil.c b/libm/powerpc/classic/s_ceil.c deleted file mode 100644 index ee4ceb5fc..000000000 --- a/libm/powerpc/classic/s_ceil.c +++ /dev/null @@ -1,113 +0,0 @@ -/******************************************************************************* -* * -* File ceilfloor.c, * -* Function ceil(x) and floor(x), * -* Implementation of ceil and floor for the PowerPC. * -* * -* Copyright © 1991 Apple Computer, Inc. All rights reserved. * -* * -* Written by Ali Sazegari, started on November 1991, * -* * -* based on math.h, library code for Macintoshes with a 68881/68882 * -* by Jim Thomas. * -* * -* W A R N I N G: This routine expects a 64 bit double model. * -* * -* December 03 1992: first rs6000 port. * -* July 14 1993: comment changes and addition of #pragma fenv_access. * -* May 06 1997: port of the ibm/taligent ceil and floor routines. * -* April 11 2001: first port to os x using gcc. * -* June 13 2001: replaced __setflm with in-line assembly * -* * -*******************************************************************************/ - -#include <math.h> -#include <endian.h> - -static const double twoTo52 = 4503599627370496.0; -static const unsigned long signMask = 0x80000000ul; - -typedef union - { - struct { -#if (__BYTE_ORDER == __BIG_ENDIAN) - unsigned long int hi; - unsigned long int lo; -#else - unsigned long int lo; - unsigned long int hi; -#endif - } words; - double dbl; - } DblInHex; - -/******************************************************************************* -* Functions needed for the computation. * -*******************************************************************************/ - -/******************************************************************************* -* Ceil(x) returns the smallest integer not less than x. * -*******************************************************************************/ - -libm_hidden_proto(ceil) -double ceil ( double x ) - { - DblInHex xInHex,OldEnvironment; - register double y; - register unsigned long int xhi; - register int target; - - xInHex.dbl = x; - xhi = xInHex.words.hi & 0x7fffffffUL; // xhi is the high half of |x| - target = ( xInHex.words.hi < signMask ); - - if ( xhi < 0x43300000ul ) -/******************************************************************************* -* Is |x| < 2.0^52? * -*******************************************************************************/ - { - if ( xhi < 0x3ff00000ul ) -/******************************************************************************* -* Is |x| < 1.0? * -*******************************************************************************/ - { - if ( ( xhi | xInHex.words.lo ) == 0ul ) // zero x is exact case - return ( x ); - else - { // inexact case - __asm__ ("mffs %0" : "=f" (OldEnvironment.dbl)); - OldEnvironment.words.lo |= 0x02000000ul; - __asm__ ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl )); - if ( target ) - return ( 1.0 ); - else - return ( -0.0 ); - } - } -/******************************************************************************* -* Is 1.0 < |x| < 2.0^52? * -*******************************************************************************/ - if ( target ) - { - y = ( x + twoTo52 ) - twoTo52; // round at binary pt. - if ( y < x ) - return ( y + 1.0 ); - else - return ( y ); - } - - else - { - y = ( x - twoTo52 ) + twoTo52; // round at binary pt. - if ( y < x ) - return ( y + 1.0 ); - else - return ( y ); - } - } -/******************************************************************************* -* |x| >= 2.0^52 or x is a NaN. * -*******************************************************************************/ - return ( x ); - } -libm_hidden_def(ceil) diff --git a/libm/powerpc/classic/s_copysign.c b/libm/powerpc/classic/s_copysign.c deleted file mode 100644 index c6f1307a3..000000000 --- a/libm/powerpc/classic/s_copysign.c +++ /dev/null @@ -1,59 +0,0 @@ -/******************************************************************************* -* * -* File sign.c, * -* Functions copysign and __signbitd. * -* For PowerPC based machines. * -* * -* Copyright © 1991, 2001 Apple Computer, Inc. All rights reserved. * -* * -* Written by Ali Sazegari, started on June 1991. * -* * -* August 26 1991: no CFront Version 1.1d17 warnings. * -* September 06 1991: passes the test suite with invalid raised on * -* signaling nans. sane rom code behaves the same. * -* September 24 1992: took the ̉#include support.hÓ out. * -* Dcember 02 1992: PowerPC port. * -* July 20 1994: __fabs added * -* July 21 1994: deleted unnecessary functions: neg, COPYSIGNnew, * -* and SIGNNUMnew. * -* April 11 2001: first port to os x using gcc. * -* removed fabs and deffered to gcc for direct * -* instruction generation. * -* * -*******************************************************************************/ - -#include <math.h> -#include "../../fp_private.h" - -/******************************************************************************* -* * -* Function copysign. * -* Implementation of copysign for the PowerPC. * -* * -******************************************************************************** -* Note: The order of the operands in this function is reversed from that * -* suggested in the IEEE standard 754. * -*******************************************************************************/ - -libm_hidden_proto(copysign) -double copysign ( double arg2, double arg1 ) - { - union - { - dHexParts hex; - double dbl; - } x, y; - -/******************************************************************************* -* No need to flush NaNs out. * -*******************************************************************************/ - - x.dbl = arg1; - y.dbl = arg2; - - y.hex.high = y.hex.high & 0x7FFFFFFF; - y.hex.high = ( y.hex.high | ( x.hex.high & dSgnMask ) ); - - return y.dbl; - } -libm_hidden_def(copysign) diff --git a/libm/powerpc/classic/s_floor.c b/libm/powerpc/classic/s_floor.c deleted file mode 100644 index 2cd720eff..000000000 --- a/libm/powerpc/classic/s_floor.c +++ /dev/null @@ -1,113 +0,0 @@ -/******************************************************************************* -* * -* File ceilfloor.c, * -* Function ceil(x) and floor(x), * -* Implementation of ceil and floor for the PowerPC. * -* * -* Copyright © 1991 Apple Computer, Inc. All rights reserved. * -* * -* Written by Ali Sazegari, started on November 1991, * -* * -* based on math.h, library code for Macintoshes with a 68881/68882 * -* by Jim Thomas. * -* * -* W A R N I N G: This routine expects a 64 bit double model. * -* * -* December 03 1992: first rs6000 port. * -* July 14 1993: comment changes and addition of #pragma fenv_access. * -* May 06 1997: port of the ibm/taligent ceil and floor routines. * -* April 11 2001: first port to os x using gcc. * -* June 13 2001: replaced __setflm with in-line assembly * -* * -*******************************************************************************/ - -#include <math.h> -#include <endian.h> - -static const double twoTo52 = 4503599627370496.0; -static const unsigned long signMask = 0x80000000ul; - -typedef union - { - struct { -#if (__BYTE_ORDER == __BIG_ENDIAN) - unsigned long int hi; - unsigned long int lo; -#else - unsigned long int lo; - unsigned long int hi; -#endif - } words; - double dbl; - } DblInHex; - -/******************************************************************************* -* Functions needed for the computation. * -*******************************************************************************/ - -/******************************************************************************* -* Floor(x) returns the largest integer not greater than x. * -*******************************************************************************/ - -libm_hidden_proto(floor) -double floor ( double x ) - { - DblInHex xInHex,OldEnvironment; - register double y; - register unsigned long int xhi; - register long int target; - - xInHex.dbl = x; - xhi = xInHex.words.hi & 0x7fffffffUL; // xhi is the high half of |x| - target = ( xInHex.words.hi < signMask ); - - if ( xhi < 0x43300000ul ) -/******************************************************************************* -* Is |x| < 2.0^52? * -*******************************************************************************/ - { - if ( xhi < 0x3ff00000ul ) -/******************************************************************************* -* Is |x| < 1.0? * -*******************************************************************************/ - { - if ( ( xhi | xInHex.words.lo ) == 0ul ) // zero x is exact case - return ( x ); - else - { // inexact case - __asm__ ("mffs %0" : "=f" (OldEnvironment.dbl)); - OldEnvironment.words.lo |= 0x02000000ul; - __asm__ ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl )); - if ( target ) - return ( 0.0 ); - else - return ( -1.0 ); - } - } -/******************************************************************************* -* Is 1.0 < |x| < 2.0^52? * -*******************************************************************************/ - if ( target ) - { - y = ( x + twoTo52 ) - twoTo52; // round at binary pt. - if ( y > x ) - return ( y - 1.0 ); - else - return ( y ); - } - - else - { - y = ( x - twoTo52 ) + twoTo52; // round at binary pt. - if ( y > x ) - return ( y - 1.0 ); - else - return ( y ); - } - } -/******************************************************************************* -* |x| >= 2.0^52 or x is a NaN. * -*******************************************************************************/ - return ( x ); - } -libm_hidden_def(floor) diff --git a/libm/powerpc/classic/s_frexp.c b/libm/powerpc/classic/s_frexp.c deleted file mode 100644 index 001aaf708..000000000 --- a/libm/powerpc/classic/s_frexp.c +++ /dev/null @@ -1,65 +0,0 @@ -/******************************************************************************* -* * -* File frexpldexp.c, * -* Functions frexp(x) and ldexp(x), * -* Implementation of frexp and ldexp functions for the PowerPC. * -* * -* Copyright © 1991 Apple Computer, Inc. All rights reserved. * -* * -* Written by Ali Sazegari, started on January 1991, * -* * -* W A R N I N G: This routine expects a 64 bit double model. * -* * -* December03 1992: first rs6000 implementation. * -* October 05 1993: added special cases for NaN and ° in frexp. * -* May 27 1997: improved the performance of frexp by eliminating the * -* switch statement. * -* June 13 2001: (ram) rewrote frexp to eliminate calls to scalb and * -* logb. * -* * -*******************************************************************************/ - -#include <limits.h> -#include <math.h> -#include <endian.h> - -static const double two54 = 1.80143985094819840000e+16; /* 0x43500000, 0x00000000 */ - -typedef union - { - struct { -#if (__BYTE_ORDER == __BIG_ENDIAN) - unsigned long int hi; - unsigned long int lo; -#else - unsigned long int lo; - unsigned long int hi; -#endif - } words; - double dbl; - } DblInHex; - -libm_hidden_proto(frexp) -double frexp ( double value, int *eptr ) - { - DblInHex argument; - unsigned long int valueHead; - - argument.dbl = value; - valueHead = argument.words.hi & 0x7fffffffUL; // valueHead <- |x| - - *eptr = 0; - if ( valueHead >= 0x7ff00000 || ( valueHead | argument.words.lo ) == 0 ) - return value; // 0, inf, or NaN - - if ( valueHead < 0x00100000 ) - { // denorm - argument.dbl = two54 * value; - valueHead = argument.words.hi &0x7fffffff; - *eptr = -54; - } - *eptr += ( valueHead >> 20 ) - 1022; - argument.words.hi = ( argument.words.hi & 0x800fffff ) | 0x3fe00000; - return argument.dbl; - } -libm_hidden_def(frexp) diff --git a/libm/powerpc/classic/s_ldexp.c b/libm/powerpc/classic/s_ldexp.c deleted file mode 100644 index 10100d7c2..000000000 --- a/libm/powerpc/classic/s_ldexp.c +++ /dev/null @@ -1,49 +0,0 @@ -/******************************************************************************* -* * -* File frexpldexp.c, * -* Functions frexp(x) and ldexp(x), * -* Implementation of frexp and ldexp functions for the PowerPC. * -* * -* Copyright © 1991 Apple Computer, Inc. All rights reserved. * -* * -* Written by Ali Sazegari, started on January 1991, * -* * -* W A R N I N G: This routine expects a 64 bit double model. * -* * -* December03 1992: first rs6000 implementation. * -* October 05 1993: added special cases for NaN and ° in frexp. * -* May 27 1997: improved the performance of frexp by eliminating the * -* switch statement. * -* June 13 2001: (ram) rewrote frexp to eliminate calls to scalb and * -* logb. * -* * -*******************************************************************************/ - -#include <limits.h> -#include <math.h> -#include <endian.h> - -typedef union - { - struct { -#if (__BYTE_ORDER == __BIG_ENDIAN) - unsigned long int hi; - unsigned long int lo; -#else - unsigned long int lo; - unsigned long int hi; -#endif - } words; - double dbl; - } DblInHex; - -libm_hidden_proto(ldexp) -double ldexp ( double value, int exp ) - { - if ( exp > SHRT_MAX ) - exp = SHRT_MAX; - else if ( exp < -SHRT_MAX ) - exp = -SHRT_MAX; - return scalb ( value, exp ); - } -libm_hidden_def(ldexp) diff --git a/libm/powerpc/classic/s_logb.c b/libm/powerpc/classic/s_logb.c deleted file mode 100644 index 81daa412e..000000000 --- a/libm/powerpc/classic/s_logb.c +++ /dev/null @@ -1,107 +0,0 @@ -/******************************************************************************* -* * -* File logb.c, * -* Functions logb. * -* Implementation of logb for the PowerPC. * -* * -* Copyright © 1991 Apple Computer, Inc. All rights reserved. * -* * -* Written by Ali Sazegari, started on June 1991, * -* * -* August 26 1991: removed CFront Version 1.1d17 warnings. * -* August 27 1991: no errors reported by the test suite. * -* November 11 1991: changed CLASSEXTENDED to the macro CLASSIFY and * -* + or - infinity to constants. * -* November 18 1991: changed the macro CLASSIFY to CLASSEXTENDEDint to * -* improve performance. * -* February 07 1992: changed bit operations to macros ( object size is * -* unchanged ). * -* September24 1992: took the "#include support.h" out. * -* December 03 1992: first rs/6000 port. * -* August 30 1992: set the divide by zero for the zero argument case. * -* October 05 1993: corrected the environment. * -* October 17 1994: replaced all environmental functions with __setflm. * -* May 28 1997: made speed improvements. * -* April 30 2001: forst mac os x port using gcc. * -* * -******************************************************************************** -* The C math library offers a similar function called "frexp". It is * -* different in details from logb, but similar in spirit. This current * -* implementation of logb follows the recommendation in IEEE Standard 854 * -* which is different in its handling of denormalized numbers from the IEEE * -* Standard 754. * -*******************************************************************************/ - -#include <math.h> -#include <endian.h> - -typedef union - { - struct { -#if (__BYTE_ORDER == __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 double twoTo52 = 4.50359962737049600e15; // 0x1p52 -static const double klTod = 4503601774854144.0; // 0x1.000008p52 -static const unsigned long int signMask = 0x80000000ul; -static const DblInHex minusInf = {{ 0xFFF00000, 0x00000000 }}; - - -/******************************************************************************* -******************************************************************************** -* L O G B * -******************************************************************************** -*******************************************************************************/ - -libm_hidden_proto(logb) -double logb ( double x ) - { - DblInHex xInHex; - long int shiftedExp; - - xInHex.dbl = x; - shiftedExp = ( xInHex.words.hi & 0x7ff00000UL ) >> 20; - - if ( shiftedExp == 2047 ) - { // NaN or INF - if ( ( ( xInHex.words.hi & signMask ) == 0 ) || ( x != x ) ) - return x; // NaN or +INF return x - else - return -x; // -INF returns +INF - } - - if ( shiftedExp != 0 ) // normal number - shiftedExp -= 1023; // unbias exponent - - else if ( x == 0.0 ) - { // zero - xInHex.words.hi = 0x0UL; // return -infinity - return ( minusInf.dbl ); - } - - else - { // subnormal number - xInHex.dbl *= twoTo52; // scale up - shiftedExp = ( xInHex.words.hi & 0x7ff00000UL ) >> 20; - shiftedExp -= 1075; // unbias exponent - } - - if ( shiftedExp == 0 ) // zero result - return ( 0.0 ); - - else - { // nonzero result - xInHex.dbl = klTod; - xInHex.words.lo += shiftedExp; - return ( xInHex.dbl - klTod ); - } - } -libm_hidden_def(logb) diff --git a/libm/powerpc/classic/s_modf.c b/libm/powerpc/classic/s_modf.c deleted file mode 100644 index c2221bc0b..000000000 --- a/libm/powerpc/classic/s_modf.c +++ /dev/null @@ -1,344 +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 rinttol, roundtol, -** modf and modfl. This file targets PowrPC 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> -#include <endian.h> - -#define SET_INVALID 0x01000000UL - -typedef union - { - struct { -#if (__BYTE_ORDER == __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 TOWARDZERO = {{ 0x00000000, 0x00000001 }}; - - -/******************************************************************************* -* * -* 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 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; - const DblInHex kTZ = {{ 0x0, 0x1 }}; - const DblInHex kUP = {{ 0x0, 0x2 }}; - - 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 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. * -*******************************************************************************/ - -libm_hidden_proto(modf) -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 ); - } - } -libm_hidden_def(modf) diff --git a/libm/powerpc/classic/s_nearbyint.c b/libm/powerpc/classic/s_nearbyint.c deleted file mode 100644 index d08430dc6..000000000 --- a/libm/powerpc/classic/s_nearbyint.c +++ /dev/null @@ -1,38 +0,0 @@ -#include <limits.h> -#include <math.h> - -/******************************************************************************* -* * -* 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. * -* * -*******************************************************************************/ - -static const double twoTo52 = 4503599627370496.0; - -libm_hidden_proto(nearbyint) -double nearbyint ( double x ) - { - double y; - double OldEnvironment; - - 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 ); - } -libm_hidden_def(nearbyint) diff --git a/libm/powerpc/classic/s_rint.c b/libm/powerpc/classic/s_rint.c deleted file mode 100644 index dd2ae585e..000000000 --- a/libm/powerpc/classic/s_rint.c +++ /dev/null @@ -1,159 +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> -#include <endian.h> - -#define SET_INVALID 0x01000000UL - -typedef union - { - struct { -#if (__BYTE_ORDER == __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 }}; - -libm_hidden_proto(rint) -/******************************************************************************* -* * -* 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 ); - } -libm_hidden_def(rint) diff --git a/libm/powerpc/classic/s_round.c b/libm/powerpc/classic/s_round.c deleted file mode 100644 index 9fd73801d..000000000 --- a/libm/powerpc/classic/s_round.c +++ /dev/null @@ -1,115 +0,0 @@ -#include <limits.h> -#include <math.h> -#include <endian.h> - -typedef union - { - struct { -#if (__BYTE_ORDER == __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; - -/******************************************************************************* -* * -* 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. * -* * -*******************************************************************************/ - -libm_hidden_proto(round) -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 ); - } -libm_hidden_def(round) diff --git a/libm/powerpc/classic/s_trunc.c b/libm/powerpc/classic/s_trunc.c deleted file mode 100644 index e9c668127..000000000 --- a/libm/powerpc/classic/s_trunc.c +++ /dev/null @@ -1,89 +0,0 @@ -#include <limits.h> -#include <math.h> -#include <endian.h> - -typedef union - { - struct { -#if (__BYTE_ORDER == __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; - -/******************************************************************************* -* * -* 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. * -* * -*******************************************************************************/ - -libm_hidden_proto(trunc) -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 ); - } -libm_hidden_def(trunc) diff --git a/libm/powerpc/classic/w_scalb.c b/libm/powerpc/classic/w_scalb.c deleted file mode 100644 index 408136001..000000000 --- a/libm/powerpc/classic/w_scalb.c +++ /dev/null @@ -1,94 +0,0 @@ -/*********************************************************************** -** File: scalb.c -** -** Contains: C source code for implementations of floating-point -** scalb functions defined in header <fp.h>. In -** particular, this file contains implementations of -** functions scalb and scalbl for double and long double -** formats on PowerPC platforms. -** -** Written by: Jon Okada, SANEitation Engineer, ext. 4-4838 -** -** Copyright: © 1992 by Apple Computer, Inc., all rights reserved -** -** Change History ( most recent first ): -** -** 28 May 97 ali made an speed improvement for large n, -** removed scalbl. -** 12 Dec 92 JPO First created. -** -***********************************************************************/ - -#include <math.h> -#include <endian.h> - -typedef union - { - struct { -#if (__BYTE_ORDER == __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 double twoTo1023 = 8.988465674311579539e307; // 0x1p1023 -static const double twoToM1022 = 2.225073858507201383e-308; // 0x1p-1022 - - -/*********************************************************************** - double scalb( double x, long int n ) returns its argument x scaled - by the factor 2^m. NaNs, signed zeros, and infinities are propagated - by this function regardless of the value of n. - - Exceptions: OVERFLOW/INEXACT or UNDERFLOW inexact may occur; - INVALID for signaling NaN inputs ( quiet NaN returned ). - - Calls: none. -***********************************************************************/ - -libm_hidden_proto(scalb) -#ifdef _SCALB_INT -double scalb ( double x, int n ) -#else -double scalb ( double x, double n ) -#endif - { - DblInHex xInHex; - - xInHex.words.lo = 0UL; // init. low half of xInHex - - if ( n > 1023 ) - { // large positive scaling - if ( n > 2097 ) // huge scaling - return ( ( x * twoTo1023 ) * twoTo1023 ) * twoTo1023; - while ( n > 1023 ) - { // scale reduction loop - x *= twoTo1023; // scale x by 2^1023 - n -= 1023; // reduce n by 1023 - } - } - - else if ( n < -1022 ) - { // large negative scaling - if ( n < -2098 ) // huge negative scaling - return ( ( x * twoToM1022 ) * twoToM1022 ) * twoToM1022; - while ( n < -1022 ) - { // scale reduction loop - x *= twoToM1022; // scale x by 2^( -1022 ) - n += 1022; // incr n by 1022 - } - } - -/******************************************************************************* -* -1022 <= n <= 1023; convert n to double scale factor. * -*******************************************************************************/ - - xInHex.words.hi = ( ( unsigned long ) ( n + 1023 ) ) << 20; - return ( x * xInHex.dbl ); - } -libm_hidden_def(scalb) |