lib/gnu/ieee_float/ldexp.c

   1 /*
   2 libc/ieee_float/ldexp.c
   3
   4 Created:        Oct 14, 1993 by Philip Homburg <philip@cs.vu.nl>
   5
   6 Implementation of ldexp that directly manipulates the exponent bits in an
   7 ieee float
   8 */
   9
  10 #include <sys/types.h>
  11 #include <errno.h>
  12 #include <math.h>
  13
  14 #include "ieee_float.h"
  15
  16 double ldexp(value, exp)
  17 double value;
  18 int exp;
  19 {
  20         struct f64 *f64p;
  21         int oldexp, exp_bias;
  22         double factor;
  23
  24         f64p= (struct f64 *)&value;
  25         exp_bias= 0;
  26
  27         oldexp= F64_GET_EXP(f64p);
  28         if (oldexp == F64_EXP_MAX)
  29         {       /* Either infinity or Nan */
  30                 return value;
  31         }
  32         if (oldexp == 0)
  33         {
  34                 /* Either 0 or denormal */
  35                 if (F64_GET_MANT_LOW(f64p) == 0 &&
  36                         F64_GET_MANT_HIGH(f64p) == 0)
  37                 {
  38                         return value;
  39                 }
  40         }
  41
  42         /* If exp is too large (> 2*F64_EXP_MAX) or too small
  43          * (< -2*F64_EXP_MAX) return HUGE_VAL or 0. This prevents overflows
  44          * in exp if exp is really weird
  45          */
  46         if (exp >= 2*F64_EXP_MAX)
  47         {
  48                 errno= ERANGE;
  49                 return HUGE_VAL;
  50         }
  51         if (exp <= -2*F64_EXP_MAX)
  52         {
  53                 errno= ERANGE;
  54                 return 0;
  55         }
  56
  57         /* Normalize a denormal */
  58         if (oldexp == 0)
  59         {
  60                 /* Multiply by 2^64 */
  61                 factor= 65536.0;        /* 2^16 */
  62                 factor *= factor;       /* 2^32 */
  63                 factor *= factor;       /* 2^64 */
  64                 value *= factor;
  65                 exp= -64;
  66                 oldexp= F64_GET_EXP(f64p);
  67         }
  68
  69         exp= oldexp + exp;
  70         if (exp >= F64_EXP_MAX)
  71         {       /* Overflow */
  72                 errno= ERANGE;
  73                 return HUGE_VAL;
  74         }
  75         if (exp > 0)
  76         {
  77                 /* Normal */
  78                 F64_SET_EXP(f64p, exp);
  79                 return value;
  80         }
  81         /* Denormal, or underflow. */
  82         exp += 64;
  83         F64_SET_EXP(f64p, exp);
  84         /* Divide by 2^64 */
  85         factor= 65536.0;        /* 2^16 */
  86         factor *= factor;       /* 2^32 */
  87         factor *= factor;       /* 2^64 */
  88         value /= factor;
  89         if (value == 0.0)
  90         {
  91                 /* Underflow */
  92                 errno= ERANGE;
  93         }
  94         return value;
  95 }
  96
  97 /*
  98  * $PchId: ldexp.c,v 1.3 1996/02/22 21:01:39 philip Exp $
  99  */