| 1 |
1325 |
phoenix |
/*******************************************************************************
|
| 2 |
|
|
* *
|
| 3 |
|
|
* File logb.c, *
|
| 4 |
|
|
* Functions logb. *
|
| 5 |
|
|
* Implementation of logb for the PowerPC. *
|
| 6 |
|
|
* *
|
| 7 |
|
|
* Copyright © 1991 Apple Computer, Inc. All rights reserved. *
|
| 8 |
|
|
* *
|
| 9 |
|
|
* Written by Ali Sazegari, started on June 1991, *
|
| 10 |
|
|
* *
|
| 11 |
|
|
* August 26 1991: removed CFront Version 1.1d17 warnings. *
|
| 12 |
|
|
* August 27 1991: no errors reported by the test suite. *
|
| 13 |
|
|
* November 11 1991: changed CLASSEXTENDED to the macro CLASSIFY and *
|
| 14 |
|
|
* + or - infinity to constants. *
|
| 15 |
|
|
* November 18 1991: changed the macro CLASSIFY to CLASSEXTENDEDint to *
|
| 16 |
|
|
* improve performance. *
|
| 17 |
|
|
* February 07 1992: changed bit operations to macros ( object size is *
|
| 18 |
|
|
* unchanged ). *
|
| 19 |
|
|
* September24 1992: took the "#include support.h" out. *
|
| 20 |
|
|
* December 03 1992: first rs/6000 port. *
|
| 21 |
|
|
* August 30 1992: set the divide by zero for the zero argument case. *
|
| 22 |
|
|
* October 05 1993: corrected the environment. *
|
| 23 |
|
|
* October 17 1994: replaced all environmental functions with __setflm. *
|
| 24 |
|
|
* May 28 1997: made speed improvements. *
|
| 25 |
|
|
* April 30 2001: forst mac os x port using gcc. *
|
| 26 |
|
|
* *
|
| 27 |
|
|
********************************************************************************
|
| 28 |
|
|
* The C math library offers a similar function called "frexp". It is *
|
| 29 |
|
|
* different in details from logb, but similar in spirit. This current *
|
| 30 |
|
|
* implementation of logb follows the recommendation in IEEE Standard 854 *
|
| 31 |
|
|
* which is different in its handling of denormalized numbers from the IEEE *
|
| 32 |
|
|
* Standard 754. *
|
| 33 |
|
|
*******************************************************************************/
|
| 34 |
|
|
|
| 35 |
|
|
typedef union
|
| 36 |
|
|
{
|
| 37 |
|
|
struct {
|
| 38 |
|
|
#if defined(__BIG_ENDIAN__)
|
| 39 |
|
|
unsigned long int hi;
|
| 40 |
|
|
unsigned long int lo;
|
| 41 |
|
|
#else
|
| 42 |
|
|
unsigned long int lo;
|
| 43 |
|
|
unsigned long int hi;
|
| 44 |
|
|
#endif
|
| 45 |
|
|
} words;
|
| 46 |
|
|
double dbl;
|
| 47 |
|
|
} DblInHex;
|
| 48 |
|
|
|
| 49 |
|
|
static const double twoTo52 = 4.50359962737049600e15; // 0x1p52
|
| 50 |
|
|
static const double klTod = 4503601774854144.0; // 0x1.000008p52
|
| 51 |
|
|
static const unsigned long int signMask = 0x80000000ul;
|
| 52 |
|
|
static const DblInHex minusInf = {{ 0xFFF00000, 0x00000000 }};
|
| 53 |
|
|
|
| 54 |
|
|
|
| 55 |
|
|
/*******************************************************************************
|
| 56 |
|
|
********************************************************************************
|
| 57 |
|
|
* L O G B *
|
| 58 |
|
|
********************************************************************************
|
| 59 |
|
|
*******************************************************************************/
|
| 60 |
|
|
|
| 61 |
|
|
double logb ( double x )
|
| 62 |
|
|
{
|
| 63 |
|
|
DblInHex xInHex;
|
| 64 |
|
|
long int shiftedExp;
|
| 65 |
|
|
|
| 66 |
|
|
xInHex.dbl = x;
|
| 67 |
|
|
shiftedExp = ( xInHex.words.hi & 0x7ff00000UL ) >> 20;
|
| 68 |
|
|
|
| 69 |
|
|
if ( shiftedExp == 2047 )
|
| 70 |
|
|
{ // NaN or INF
|
| 71 |
|
|
if ( ( ( xInHex.words.hi & signMask ) == 0 ) || ( x != x ) )
|
| 72 |
|
|
return x; // NaN or +INF return x
|
| 73 |
|
|
else
|
| 74 |
|
|
return -x; // -INF returns +INF
|
| 75 |
|
|
}
|
| 76 |
|
|
|
| 77 |
|
|
if ( shiftedExp != 0 ) // normal number
|
| 78 |
|
|
shiftedExp -= 1023; // unbias exponent
|
| 79 |
|
|
|
| 80 |
|
|
else if ( x == 0.0 )
|
| 81 |
|
|
{ // zero
|
| 82 |
|
|
xInHex.words.hi = 0x0UL; // return -infinity
|
| 83 |
|
|
return ( minusInf.dbl );
|
| 84 |
|
|
}
|
| 85 |
|
|
|
| 86 |
|
|
else
|
| 87 |
|
|
{ // subnormal number
|
| 88 |
|
|
xInHex.dbl *= twoTo52; // scale up
|
| 89 |
|
|
shiftedExp = ( xInHex.words.hi & 0x7ff00000UL ) >> 20;
|
| 90 |
|
|
shiftedExp -= 1075; // unbias exponent
|
| 91 |
|
|
}
|
| 92 |
|
|
|
| 93 |
|
|
if ( shiftedExp == 0 ) // zero result
|
| 94 |
|
|
return ( 0.0 );
|
| 95 |
|
|
|
| 96 |
|
|
else
|
| 97 |
|
|
{ // nonzero result
|
| 98 |
|
|
xInHex.dbl = klTod;
|
| 99 |
|
|
xInHex.words.lo += shiftedExp;
|
| 100 |
|
|
return ( xInHex.dbl - klTod );
|
| 101 |
|
|
}
|
| 102 |
|
|
}
|
| 103 |
|
|
|
