gonzui

    1: /*      $OpenBSD: ldexp.c,v 1.5 2011/07/26 11:43:01 martynas Exp $   */
    2: /* @(#)s_scalbn.c 5.1 93/09/24 */
    3: /* @(#)fdlibm.h 5.1 93/09/24 */
    4: /*
    5:  * ====================================================
    6:  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
    7:  *
    8:  * Developed at SunPro, a Sun Microsystems, Inc. business.
    9:  * Permission to use, copy, modify, and distribute this
   10:  * software is freely granted, provided that this notice
   11:  * is preserved.
   12:  * ====================================================
   13:  */
   14: 
   15: /* LINTLIBRARY */
   16: 
   17: #include <sys/types.h>
   18: #include <sys/cdefs.h>
   19: #include <machine/endian.h>
   20: #include <float.h>
   21: #include <math.h>
   22: 
   23: /* Bit fiddling routines copied from msun/src/math_private.h,v 1.15 */
   24: 
   25: #if (BYTE_ORDER == BIG_ENDIAN) || (defined(__arm__) && !defined(__VFP_FP__))
   26: 
   27: typedef union
   28: {
   29:   double value;
   30:   struct
   31:   {
   32:     u_int32_t msw;
   33:     u_int32_t lsw;
   34:   } parts;
   35: } ieee_double_shape_type;
   36: 
   37: #endif
   38: 
   39: #if (BYTE_ORDER == LITTLE_ENDIAN) && !(defined(__arm__) && !defined(__VFP_FP__))
   40: 
   41: typedef union
   42: {
   43:   double value;
   44:   struct
   45:   {
   46:     u_int32_t lsw;
   47:     u_int32_t msw;
   48:   } parts;
   49: } ieee_double_shape_type;
   50: 
   51: #endif
   52: 
   53: /* Get two 32 bit ints from a double.  */
   54: 
   55: #define EXTRACT_WORDS(ix0,ix1,d)                                \
   56: do {                                                            \
   57:   ieee_double_shape_type ew_u;                                  \
   58:   ew_u.value = (d);                                             \
   59:   (ix0) = ew_u.parts.msw;                                       \
   60:   (ix1) = ew_u.parts.lsw;                                       \
   61: } while (0)
   62: 
   63: /* Get the more significant 32 bit int from a double.  */
   64: 
   65: #define GET_HIGH_WORD(i,d)                                      \
   66: do {                                                            \
   67:   ieee_double_shape_type gh_u;                                  \
   68:   gh_u.value = (d);                                             \
   69:   (i) = gh_u.parts.msw;                                         \
   70: } while (0)
   71: 
   72: /* Set the more significant 32 bits of a double from an int.  */
   73: 
   74: #define SET_HIGH_WORD(d,v)                                      \
   75: do {                                                            \
   76:   ieee_double_shape_type sh_u;                                  \
   77:   sh_u.value = (d);                                             \
   78:   sh_u.parts.msw = (v);                                         \
   79:   (d) = sh_u.value;                                             \
   80: } while (0)
   81: 
   82: 
   83: static const double
   84: two54   =  1.80143985094819840000e+16, /* 0x43500000, 0x00000000 */
   85: twom54  =  5.55111512312578270212e-17, /* 0x3C900000, 0x00000000 */
   86: huge   = 1.0e+300,
   87: tiny   = 1.0e-300;
   88: 
   89: static double
   90: _copysign(double x, double y)
   91: {
   92:         u_int32_t hx,hy;
   93:         GET_HIGH_WORD(hx,x);
   94:         GET_HIGH_WORD(hy,y);
   95:         SET_HIGH_WORD(x,(hx&0x7fffffff)|(hy&0x80000000));
   96:         return x;
   97: }
   98: 
   99: double
  100: ldexp(double x, int n)
  101: {
  102:         int32_t k,hx,lx;
  103:         EXTRACT_WORDS(hx,lx,x);
  104:         k = (hx&0x7ff00000)>>20;                /* extract exponent */
  105:         if (k==0) {                             /* 0 or subnormal x */
  106:             if ((lx|(hx&0x7fffffff))==0) return x; /* +-0 */
  107:             x *= two54;
  108:             GET_HIGH_WORD(hx,x);
  109:             k = ((hx&0x7ff00000)>>20) - 54;
  110:             if (n< -50000) return tiny*x;       /*underflow*/
  111:             }
  112:         if (k==0x7ff) return x+x;               /* NaN or Inf */
  113:         k = k+n;
  114:         if (k >  0x7fe) return huge*_copysign(huge,x); /* overflow  */
  115:         if (k > 0)                              /* normal result */
  116:             {SET_HIGH_WORD(x,(hx&0x800fffff)|(k<<20)); return x;}
  117:         if (k <= -54) {
  118:             if (n > 50000)      /* in case integer overflow in n+k */
  119:                 return huge*_copysign(huge,x);        /*overflow*/
  120:             else return tiny*_copysign(tiny,x);        /*underflow*/
  121:         }
  122:         k += 54;                                /* subnormal result */
  123:         SET_HIGH_WORD(x,(hx&0x800fffff)|(k<<20));
  124:         return x*twom54;
  125: }
  126: 
  127: #if     LDBL_MANT_DIG == 53
  128: #ifdef  lint
  129: /* PROTOLIB1 */
  130: long double ldexpl(long double, int);
  131: #else   /* lint */
  132: __weak_alias(ldexpl, ldexp);
  133: #endif  /* lint */
  134: #endif  /* LDBL_MANT_DIG == 53 */