lib/libc/gen/ldexp.c

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