libm/powerpc/classic/w_scalb.c

   1 /***********************************************************************
   2 **      File:    scalb.c
   3 **
   4 **      Contains: C source code for implementations of floating-point
   5 **                scalb functions defined in header <fp.h>.  In
   6 **                particular, this file contains implementations of
   7 **                functions scalb and scalbl for double and long double
   8 **                formats on PowerPC platforms.
   9 **
  10 **      Written by: Jon Okada, SANEitation Engineer, ext. 4-4838
  11 **
  12 **      Copyright: © 1992 by Apple Computer, Inc., all rights reserved
  13 **
  14 **      Change History ( most recent first ):
  15 **
  16 **      28 May 97  ali   made an speed improvement for large n,
  17 **                       removed scalbl.
  18 **      12 Dec 92  JPO   First created.
  19 **
  20 ***********************************************************************/
  21
  22 #include <math.h>
  23 #include <endian.h>
  24
  25 typedef union
  26       {
  27       struct {
  28 #if (__BYTE_ORDER == __BIG_ENDIAN)
  29         unsigned long int hi;
  30         unsigned long int lo;
  31 #else
  32         unsigned long int lo;
  33         unsigned long int hi;
  34 #endif
  35       } words;
  36       double dbl;
  37       } DblInHex;
  38
  39 static const double twoTo1023  = 8.988465674311579539e307;   // 0x1p1023
  40 static const double twoToM1022 = 2.225073858507201383e-308;  // 0x1p-1022
  41
  42
  43 /***********************************************************************
  44       double  scalb( double  x, long int n ) returns its argument x scaled
  45       by the factor 2^m.  NaNs, signed zeros, and infinities are propagated
  46       by this function regardless of the value of n.
  47
  48       Exceptions:  OVERFLOW/INEXACT or UNDERFLOW inexact may occur;
  49                          INVALID for signaling NaN inputs ( quiet NaN returned ).
  50
  51       Calls:  none.
  52 ***********************************************************************/
  53
  54 libm_hidden_proto(scalb)
  55 #ifdef _SCALB_INT
  56 double scalb ( double x, int n  )
  57 #else
  58 double scalb ( double x, double n  )
  59 #endif
  60       {
  61       DblInHex xInHex;
  62
  63       xInHex.words.lo = 0UL;                     // init. low half of xInHex
  64
  65       if ( n > 1023 )
  66             {                                   // large positive scaling
  67             if ( n > 2097 )                     // huge scaling
  68                 return ( ( x * twoTo1023 ) * twoTo1023 ) * twoTo1023;
  69             while ( n > 1023 )
  70                   {                             // scale reduction loop
  71                   x *= twoTo1023;               // scale x by 2^1023
  72                   n -= 1023;                    // reduce n by 1023
  73                   }
  74             }
  75
  76       else if ( n < -1022 )
  77             {                                   // large negative scaling
  78             if ( n < -2098 )                    // huge negative scaling
  79                 return ( ( x * twoToM1022 ) * twoToM1022 ) * twoToM1022;
  80             while ( n < -1022 )
  81                   {                             // scale reduction loop
  82                   x *= twoToM1022;              // scale x by 2^( -1022 )
  83                   n += 1022;                    // incr n by 1022
  84                   }
  85             }
  86
  87 /*******************************************************************************
  88 *      -1022 <= n <= 1023; convert n to double scale factor.                   *
  89 *******************************************************************************/
  90
  91       xInHex.words.hi = ( ( unsigned long ) ( n + 1023 ) ) << 20;
  92       return ( x * xInHex.dbl );
  93       }
  94 libm_hidden_def(scalb)