diff options
Diffstat (limited to 'libm/powerpc/s_logb.c')
| -rw-r--r-- | libm/powerpc/s_logb.c | 103 |
1 files changed, 103 insertions, 0 deletions
diff --git a/libm/powerpc/s_logb.c b/libm/powerpc/s_logb.c new file mode 100644 index 000000000..3a410ba5c --- /dev/null +++ b/libm/powerpc/s_logb.c @@ -0,0 +1,103 @@ +/******************************************************************************* +* * +* 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 ); + } + } + |