123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103 |
- /*******************************************************************************
- * *
- * 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 );
- }
- }
|