summaryrefslogtreecommitdiff
path: root/libm/logb.c
diff options
context:
space:
mode:
Diffstat (limited to 'libm/logb.c')
-rw-r--r--libm/logb.c104
1 files changed, 104 insertions, 0 deletions
diff --git a/libm/logb.c b/libm/logb.c
new file mode 100644
index 000000000..da2a27d72
--- /dev/null
+++ b/libm/logb.c
@@ -0,0 +1,104 @@
+#if defined(__ppc__)
+/*******************************************************************************
+* *
+* 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 );
+ }
+ }
+#endif /* __ppc__ */