Line data    Source code

       1             : /* @(#)s_scalbn.c 5.1 93/09/24 */
       2             : /*
       3             :  * ====================================================
       4             :  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
       5             :  *
       6             :  * Developed at SunPro, a Sun Microsystems, Inc. business.
       7             :  * Permission to use, copy, modify, and distribute this
       8             :  * software is freely granted, provided that this notice
       9             :  * is preserved.
      10             :  * ====================================================
      11             :  */
      12             : 
      13             : #ifndef lint
      14             : //static char rcsid[] = "$FreeBSD$";
      15             : #endif
      16             : 
      17             : /*
      18             :  * scalbn (double x, int n)
      19             :  * scalbn(x,n) returns x* 2**n  computed by  exponent
      20             :  * manipulation rather than by actually performing an
      21             :  * exponentiation or a multiplication.
      22             :  */
      23             : 
      24             : //#include <sys/cdefs.h>
      25             : #include <float.h>
      26             : 
      27             : #include "math_private.h"
      28             : 
      29             : static const double
      30             : two54   =  1.80143985094819840000e+16, /* 0x43500000, 0x00000000 */
      31             : twom54  =  5.55111512312578270212e-17, /* 0x3C900000, 0x00000000 */
      32             : huge   = 1.0e+300,
      33             : tiny   = 1.0e-300;
      34             : 
      35             : double
      36           0 : scalbn (double x, int n)
      37             : {
      38             :         int32_t k,hx,lx;
      39           0 :         EXTRACT_WORDS(hx,lx,x);
      40           0 :         k = (hx&0x7ff00000)>>20;              /* extract exponent */
      41           0 :         if (k==0) {                             /* 0 or subnormal x */
      42           0 :             if ((lx|(hx&0x7fffffff))==0) return x; /* +-0 */
      43           0 :             x *= two54;
      44           0 :             GET_HIGH_WORD(hx,x);
      45           0 :             k = ((hx&0x7ff00000)>>20) - 54;
      46           0 :             if (n< -50000) return tiny*x;    /*underflow*/
      47             :             }
      48           0 :         if (k==0x7ff) return x+x;               /* NaN or Inf */
      49           0 :         k = k+n;
      50           0 :         if (k >  0x7fe) return huge*copysign(huge,x); /* overflow  */
      51           0 :         if (k > 0)                           /* normal result */
      52           0 :             {SET_HIGH_WORD(x,(hx&0x800fffff)|(k<<20)); return x;}
      53           0 :         if (k <= -54)
      54           0 :             if (n > 50000)   /* in case integer overflow in n+k */
      55           0 :                 return huge*copysign(huge,x);   /*overflow*/
      56           0 :             else return tiny*copysign(tiny,x);  /*underflow*/
      57           0 :         k += 54;                                /* subnormal result */
      58           0 :         SET_HIGH_WORD(x,(hx&0x800fffff)|(k<<20));
      59           0 :         return x*twom54;
      60             : }