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