libm/logb.c

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