lib/libc/alpha/gen/ldexp.c

   1 /*      $NetBSD: ldexp.c,v 1.1 1995/02/10 17:50:24 cgd Exp $    */
   2 /* $FreeBSD: src/lib/libc/alpha/gen/ldexp.c,v 1.1.1.1.6.1 2000/08/21 21:09:29 jhb Exp $ */
   3
   4 /*
   5  * Copyright (c) 1994, 1995 Carnegie-Mellon University.
   6  * All rights reserved.
   7  *
   8  * Author: Chris G. Demetriou
   9  *
  10  * Permission to use, copy, modify and distribute this software and
  11  * its documentation is hereby granted, provided that both the copyright
  12  * notice and this permission notice appear in all copies of the
  13  * software, derivative works or modified versions, and any portions
  14  * thereof, and that both notices appear in supporting documentation.
  15  *
  16  * CARNEGIE MELLON ALLOWS FREE USE OF THIS SOFTWARE IN ITS "AS IS"
  17  * CONDITION.  CARNEGIE MELLON DISCLAIMS ANY LIABILITY OF ANY KIND
  18  * FOR ANY DAMAGES WHATSOEVER RESULTING FROM THE USE OF THIS SOFTWARE.
  19  *
  20  * Carnegie Mellon requests users of this software to return to
  21  *
  22  *  Software Distribution Coordinator  or  Software.Distribution@CS.CMU.EDU
  23  *  School of Computer Science
  24  *  Carnegie Mellon University
  25  *  Pittsburgh PA 15213-3890
  26  *
  27  * any improvements or extensions that they make and grant Carnegie the
  28  * rights to redistribute these changes.
  29  */
  30
  31 #include <sys/types.h>
  32 #include <machine/ieee.h>
  33 #include <errno.h>
  34 #include <math.h>
  35
  36 /*
  37  * double ldexp(double val, int exp)
  38  * returns: val * (2**exp)
  39  */
  40 double
  41 ldexp(val, exp)
  42         double val;
  43         int exp;
  44 {
  45         register int oldexp, newexp, mulexp;
  46         union doub {
  47                 double v;
  48                 struct ieee_double s;
  49         } u, mul;
  50
  51         /*
  52          * If input is zero, or no change, just return input.
  53          * Likewise, if input is Inf or NaN, just return it.
  54          */
  55         u.v = val;
  56         oldexp = u.s.dbl_exp;
  57         if (val == 0 || exp == 0 || oldexp == DBL_EXP_INFNAN)
  58                 return (val);
  59
  60         /*
  61          * Compute new exponent and check for over/under flow.
  62          * Underflow, unfortunately, could mean switching to denormal.
  63          * If result out of range, set ERANGE and return 0 if too small
  64          * or Inf if too big, with the same sign as the input value.
  65          */
  66         newexp = oldexp + exp;
  67         if (newexp >= DBL_EXP_INFNAN) {
  68                 /* u.s.dbl_sign = val < 0; -- already set */
  69                 u.s.dbl_exp = DBL_EXP_INFNAN;
  70                 u.s.dbl_frach = u.s.dbl_fracl = 0;
  71                 errno = ERANGE;
  72                 return (u.v);           /* Inf */
  73         }
  74         if (newexp <= 0) {
  75                 /*
  76                  * The output number is either a denormal or underflows
  77                  * (see comments in machine/ieee.h).
  78                  */
  79                 if (newexp <= -DBL_FRACBITS) {
  80                         /* u.s.dbl_sign = val < 0; -- already set */
  81                         u.s.dbl_exp = 0;
  82                         u.s.dbl_frach = u.s.dbl_fracl = 0;
  83                         errno = ERANGE;
  84                         return (u.v);           /* zero */
  85                 }
  86                 /*
  87                  * We are going to produce a denorm.  Our `exp' argument
  88                  * might be as small as -2097, and we cannot compute
  89                  * 2^-2097, so we may have to do this as many as three
  90                  * steps (not just two, as for positive `exp's below).
  91                  */
  92                 mul.v = 0;
  93                 while (exp <= -DBL_EXP_BIAS) {
  94                         mul.s.dbl_exp = 1;
  95                         val *= mul.v;
  96                         exp += DBL_EXP_BIAS - 1;
  97                 }
  98                 mul.s.dbl_exp = exp + DBL_EXP_BIAS;
  99                 val *= mul.v;
 100                 return (val);
 101         }
 102
 103         /*
 104          * Newexp is positive.
 105          *
 106          * If oldexp is zero, we are starting with a denorm, and simply
 107          * adjusting the exponent will produce bogus answers.  We need
 108          * to fix that first.
 109          */
 110         if (oldexp == 0) {
 111                 /*
 112                  * Multiply by 2^mulexp to make the number normalizable.
 113                  * We cannot multiply by more than 2^1023, but `exp'
 114                  * argument might be as large as 2046.  A single
 115                  * adjustment, however, will normalize the number even
 116                  * for huge `exp's, and then we can use exponent
 117                  * arithmetic just as for normal `double's.
 118                  */
 119                 mulexp = exp <= DBL_EXP_BIAS ? exp : DBL_EXP_BIAS;
 120                 mul.v = 0;
 121                 mul.s.dbl_exp = mulexp + DBL_EXP_BIAS;
 122                 val *= mul.v;
 123                 if (mulexp == exp)
 124                         return (val);
 125                 u.v = val;
 126                 newexp -= mulexp;
 127         }
 128
 129         /*
 130          * Both oldexp and newexp are positive; just replace the
 131          * old exponent with the new one.
 132          */
 133         u.s.dbl_exp = newexp;
 134         return (u.v);
 135 }