From 88e814e997c9bc4fd633583d4d6d4779c5b18cfb Mon Sep 17 00:00:00 2001 From: Eric Andersen Date: Sat, 24 Nov 2001 04:17:07 +0000 Subject: Fix naming so things will actually work -Erik --- libm/powerpc/Makefile | 2 +- libm/powerpc/ceilfloor.c | 174 ------------- libm/powerpc/frexpldexp.c | 72 ------ libm/powerpc/logb.c | 103 -------- libm/powerpc/rndint.c | 620 ---------------------------------------------- libm/powerpc/s_ceil.c | 109 ++++++++ libm/powerpc/s_copysign.c | 56 +++++ libm/powerpc/s_floor.c | 109 ++++++++ libm/powerpc/s_frexp.c | 63 +++++ libm/powerpc/s_ldexp.c | 47 ++++ libm/powerpc/s_logb.c | 103 ++++++++ libm/powerpc/s_modf.c | 535 +++++++++++++++++++++++++++++++++++++++ libm/powerpc/s_rint.c | 157 ++++++++++++ libm/powerpc/scalb.c | 85 ------- libm/powerpc/sign.c | 56 ----- libm/powerpc/w_scalb.c | 85 +++++++ 16 files changed, 1265 insertions(+), 1111 deletions(-) delete mode 100644 libm/powerpc/ceilfloor.c delete mode 100644 libm/powerpc/frexpldexp.c delete mode 100644 libm/powerpc/logb.c delete mode 100644 libm/powerpc/rndint.c create mode 100644 libm/powerpc/s_ceil.c create mode 100644 libm/powerpc/s_copysign.c create mode 100644 libm/powerpc/s_floor.c create mode 100644 libm/powerpc/s_frexp.c create mode 100644 libm/powerpc/s_ldexp.c create mode 100644 libm/powerpc/s_logb.c create mode 100644 libm/powerpc/s_modf.c create mode 100644 libm/powerpc/s_rint.c delete mode 100644 libm/powerpc/scalb.c delete mode 100644 libm/powerpc/sign.c create mode 100644 libm/powerpc/w_scalb.c (limited to 'libm/powerpc') diff --git a/libm/powerpc/Makefile b/libm/powerpc/Makefile index 9872d6077..dec4c49f3 100644 --- a/libm/powerpc/Makefile +++ b/libm/powerpc/Makefile @@ -27,7 +27,7 @@ TARGET_CC= $(TOPDIR)extra/gcc-uClibc/$(TARGET_ARCH)-uclibc-gcc TARGET_CFLAGS+=-D_IEEE_LIBM -D_ISOC99_SOURCE -D_SVID_SOURCE ifeq ($(strip $(DO_C99_MATH)),true) -CSRC = ceilfloor.c frexpldexp.c logb.c rndint.c scalb.c sign.c +CSRC = s_ceil.c s_floor.c s_ldexp.c s_frexp.c s_logb.c s_modf.c w_scalb.c s_copysign.c s_rint.c else CSRC = endif diff --git a/libm/powerpc/ceilfloor.c b/libm/powerpc/ceilfloor.c deleted file mode 100644 index d5947ab50..000000000 --- a/libm/powerpc/ceilfloor.c +++ /dev/null @@ -1,174 +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 * -* * -*******************************************************************************/ - -static const double twoTo52 = 4503599627370496.0; -static const unsigned long signMask = 0x80000000ul; - -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; - -/******************************************************************************* -* Functions needed for the computation. * -*******************************************************************************/ - -/******************************************************************************* -* Ceil(x) returns the smallest integer not less than x. * -*******************************************************************************/ - -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 ); - } - -/******************************************************************************* -* Floor(x) returns the largest integer not greater than x. * -*******************************************************************************/ - -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 ); - } - diff --git a/libm/powerpc/frexpldexp.c b/libm/powerpc/frexpldexp.c deleted file mode 100644 index 6b759f91b..000000000 --- a/libm/powerpc/frexpldexp.c +++ /dev/null @@ -1,72 +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 -#include - -static const double two54 = 1.80143985094819840000e+16; /* 0x43500000, 0x00000000 */ - -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; - -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 ); - } - -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; - } - diff --git a/libm/powerpc/logb.c b/libm/powerpc/logb.c deleted file mode 100644 index 3a410ba5c..000000000 --- a/libm/powerpc/logb.c +++ /dev/null @@ -1,103 +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. * -*******************************************************************************/ - -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 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 * -******************************************************************************** -*******************************************************************************/ - -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 ); - } - } - diff --git a/libm/powerpc/rndint.c b/libm/powerpc/rndint.c deleted file mode 100644 index 1b7a9ec81..000000000 --- a/libm/powerpc/rndint.c +++ /dev/null @@ -1,620 +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 . 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 -#include - -#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 ); - } - -/******************************************************************************* -* * -* 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; - 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 ); - } - -/******************************************************************************* -* * -* 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; - 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 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. * -*******************************************************************************/ - -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 ); - } - } diff --git a/libm/powerpc/s_ceil.c b/libm/powerpc/s_ceil.c new file mode 100644 index 000000000..790dd7ba0 --- /dev/null +++ b/libm/powerpc/s_ceil.c @@ -0,0 +1,109 @@ +/******************************************************************************* +* * +* 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 * +* * +*******************************************************************************/ + +static const double twoTo52 = 4503599627370496.0; +static const unsigned long signMask = 0x80000000ul; + +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; + +/******************************************************************************* +* Functions needed for the computation. * +*******************************************************************************/ + +/******************************************************************************* +* Ceil(x) returns the smallest integer not less than x. * +*******************************************************************************/ + +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 ); + } + diff --git a/libm/powerpc/s_copysign.c b/libm/powerpc/s_copysign.c new file mode 100644 index 000000000..18e91998a --- /dev/null +++ b/libm/powerpc/s_copysign.c @@ -0,0 +1,56 @@ +/******************************************************************************* +* * +* 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 "../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. * +*******************************************************************************/ + +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; + } diff --git a/libm/powerpc/s_floor.c b/libm/powerpc/s_floor.c new file mode 100644 index 000000000..f7b4534bb --- /dev/null +++ b/libm/powerpc/s_floor.c @@ -0,0 +1,109 @@ +/******************************************************************************* +* * +* 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 * +* * +*******************************************************************************/ + +static const double twoTo52 = 4503599627370496.0; +static const unsigned long signMask = 0x80000000ul; + +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; + +/******************************************************************************* +* Functions needed for the computation. * +*******************************************************************************/ + +/******************************************************************************* +* Floor(x) returns the largest integer not greater than x. * +*******************************************************************************/ + +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 ); + } + diff --git a/libm/powerpc/s_frexp.c b/libm/powerpc/s_frexp.c new file mode 100644 index 000000000..0ca7c1514 --- /dev/null +++ b/libm/powerpc/s_frexp.c @@ -0,0 +1,63 @@ +/******************************************************************************* +* * +* 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 +#include + +static const double two54 = 1.80143985094819840000e+16; /* 0x43500000, 0x00000000 */ + +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; + +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; + } + diff --git a/libm/powerpc/s_ldexp.c b/libm/powerpc/s_ldexp.c new file mode 100644 index 000000000..0e9a66611 --- /dev/null +++ b/libm/powerpc/s_ldexp.c @@ -0,0 +1,47 @@ +/******************************************************************************* +* * +*