nsz Git - libm/blob - src/math/scalbn.c

   1 /* origin: FreeBSD /usr/src/lib/msun/src/s_scalbn.c */
   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  * scalbn (double x, int n)
  14  * scalbn(x,n) returns x* 2**n  computed by  exponent
  15  * manipulation rather than by actually performing an
  16  * exponentiation or a multiplication.
  17  */
  18
  19 #include "libm.h"
  20
  21 static const double
  22 two54  = 1.80143985094819840000e+16, /* 0x43500000, 0x00000000 */
  23 twom54 = 5.55111512312578270212e-17, /* 0x3C900000, 0x00000000 */
  24 huge   = 1.0e+300,
  25 tiny   = 1.0e-300;
  26
  27 double scalbn(double x, int n)
  28 {
  29 // FIXME: k+n check depends on signed int overflow.. use unsigned hx
  30 // TODO: when long != int:
  31 //       scalbln(x,long n) { if(n>9999)n=9999; else if(n<-9999)n=-9999; return scalbn(x,n); }
  32 // TODO: n < -50000 ...
  33         int32_t k,hx,lx;
  34
  35         EXTRACT_WORDS(hx, lx, x);
  36         k = (hx&0x7ff00000)>>20;                /* extract exponent */
  37         if (k == 0) {                           /* 0 or subnormal x */
  38                 if ((lx|(hx&0x7fffffff)) == 0)  /* +-0 */
  39                         return x;
  40                 x *= two54;
  41                 GET_HIGH_WORD(hx, x);
  42                 k = ((hx&0x7ff00000)>>20) - 54;
  43                 if (n < -50000)
  44                         return tiny*x;          /*underflow*/
  45         }
  46         if (k == 0x7ff)                         /* NaN or Inf */
  47                 return x + x;
  48         k = k + n;
  49         if (k > 0x7fe)
  50                 return huge*copysign(huge, x);  /* overflow  */
  51         if (k > 0) {                            /* normal result */
  52                 SET_HIGH_WORD(x, (hx&0x800fffff)|(k<<20));
  53                 return x;
  54         }
  55         if (k <= -54)
  56                 if (n > 50000)      /* in case integer overflow in n+k */
  57                         return huge*copysign(huge, x);  /*overflow*/
  58                 return tiny*copysign(tiny, x);  /*underflow*/
  59         k += 54;                                /* subnormal result */
  60         SET_HIGH_WORD(x, (hx&0x800fffff)|(k<<20));
  61         return x*twom54;
  62 }