1 /* 2 * Copyright (c) 1992, 1993 3 * The Regents of the University of California. All rights reserved. 4 * 5 * This software was developed by the Computer Systems Engineering group 6 * at Lawrence Berkeley Laboratory under DARPA contract BG 91-66 and 7 * contributed to Berkeley. 8 * 9 * %sccs.include.redist.c% 10 * 11 * from: $Header: ldexp.c,v 1.1 91/07/07 04:28:19 torek Exp $ 12 */ 13 14 #if defined(LIBC_SCCS) && !defined(lint) 15 static const char sccsid[] = "@(#)ldexp.c 8.1 (Berkeley) 06/04/93"; 16 #endif /* LIBC_SCCS and not lint */ 17 18 #include <sys/types.h> 19 #include <machine/ieee.h> 20 #include <errno.h> 21 22 /* 23 * double ldexp(double val, int exp) 24 * returns: val * (2**exp) 25 */ 26 double 27 ldexp(val, exp) 28 double val; 29 int exp; 30 { 31 register int oldexp, newexp, mulexp; 32 union doub { 33 double v; 34 struct ieee_double s; 35 } u, mul; 36 37 /* 38 * If input is zero, or no change, just return input. 39 * Likewise, if input is Inf or NaN, just return it. 40 */ 41 u.v = val; 42 oldexp = u.s.dbl_exp; 43 if (val == 0 || exp == 0 || oldexp == DBL_EXP_INFNAN) 44 return (val); 45 46 /* 47 * Compute new exponent and check for over/under flow. 48 * Underflow, unfortunately, could mean switching to denormal. 49 * If result out of range, set ERANGE and return 0 if too small 50 * or Inf if too big, with the same sign as the input value. 51 */ 52 newexp = oldexp + exp; 53 if (newexp >= DBL_EXP_INFNAN) { 54 /* u.s.dbl_sign = val < 0; -- already set */ 55 u.s.dbl_exp = DBL_EXP_INFNAN; 56 u.s.dbl_frach = u.s.dbl_fracl = 0; 57 errno = ERANGE; 58 return (u.v); /* Inf */ 59 } 60 if (newexp <= 0) { 61 /* 62 * The output number is either a denormal or underflows 63 * (see comments in machine/ieee.h). 64 */ 65 if (newexp <= -DBL_FRACBITS) { 66 /* u.s.dbl_sign = val < 0; -- already set */ 67 u.s.dbl_exp = 0; 68 u.s.dbl_frach = u.s.dbl_fracl = 0; 69 errno = ERANGE; 70 return (u.v); /* zero */ 71 } 72 /* 73 * We are going to produce a denorm. Our `exp' argument 74 * might be as small as -2097, and we cannot compute 75 * 2^-2097, so we may have to do this as many as three 76 * steps (not just two, as for positive `exp's below). 77 */ 78 mul.v = 0; 79 while (exp <= -DBL_EXP_BIAS) { 80 mul.s.dbl_exp = 1; 81 val *= mul.v; 82 exp += DBL_EXP_BIAS - 1; 83 } 84 mul.s.dbl_exp = exp + DBL_EXP_BIAS; 85 val *= mul.v; 86 return (val); 87 } 88 89 /* 90 * Newexp is positive. 91 * 92 * If oldexp is zero, we are starting with a denorm, and simply 93 * adjusting the exponent will produce bogus answers. We need 94 * to fix that first. 95 */ 96 if (oldexp == 0) { 97 /* 98 * Multiply by 2^mulexp to make the number normalizable. 99 * We cannot multiply by more than 2^1023, but `exp' 100 * argument might be as large as 2046. A single 101 * adjustment, however, will normalize the number even 102 * for huge `exp's, and then we can use exponent 103 * arithmetic just as for normal `double's. 104 */ 105 mulexp = exp <= DBL_EXP_BIAS ? exp : DBL_EXP_BIAS; 106 mul.v = 0; 107 mul.s.dbl_exp = mulexp + DBL_EXP_BIAS; 108 val *= mul.v; 109 if (mulexp == exp) 110 return (val); 111 u.v = val; 112 newexp -= mulexp; 113 } 114 115 /* 116 * Both oldexp and newexp are positive; just replace the 117 * old exponent with the new one. 118 */ 119 u.s.dbl_exp = newexp; 120 return (u.v); 121 } 122