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