Annotation of src/lib/libc/gen/ldexp_ieee754.c, Revision 1.1
1.1 ! kleink 1: /* $NetBSD$ */
! 2:
! 3: /*-
! 4: * Copyright (c) 1999 The NetBSD Foundation, Inc.
! 5: * All rights reserved.
! 6: *
! 7: * This code is derived from software contributed to The NetBSD Foundation
! 8: * by Charles M. Hannum.
! 9: *
! 10: * Redistribution and use in source and binary forms, with or without
! 11: * modification, are permitted provided that the following conditions
! 12: * are met:
! 13: * 1. Redistributions of source code must retain the above copyright
! 14: * notice, this list of conditions and the following disclaimer.
! 15: * 2. Redistributions in binary form must reproduce the above copyright
! 16: * notice, this list of conditions and the following disclaimer in the
! 17: * documentation and/or other materials provided with the distribution.
! 18: * 3. All advertising materials mentioning features or use of this software
! 19: * must display the following acknowledgement:
! 20: * This product includes software developed by the NetBSD
! 21: * Foundation, Inc. and its contributors.
! 22: * 4. Neither the name of The NetBSD Foundation nor the names of its
! 23: * contributors may be used to endorse or promote products derived
! 24: * from this software without specific prior written permission.
! 25: *
! 26: * THIS SOFTWARE IS PROVIDED BY THE NETBSD FOUNDATION, INC. AND CONTRIBUTORS
! 27: * ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
! 28: * TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
! 29: * PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE FOUNDATION OR CONTRIBUTORS
! 30: * BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
! 31: * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
! 32: * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
! 33: * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
! 34: * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
! 35: * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
! 36: * POSSIBILITY OF SUCH DAMAGE.
! 37: */
! 38:
! 39: #include <sys/cdefs.h>
! 40: #if defined(LIBC_SCCS) && !defined(lint)
! 41: __RCSID("$NetBSD$");
! 42: #endif /* LIBC_SCCS and not lint */
! 43:
! 44: #include <sys/types.h>
! 45: #include <machine/ieee.h>
! 46: #include <errno.h>
! 47: #include <math.h>
! 48:
! 49: /*
! 50: * Multiply the given value by 2^expon.
! 51: */
! 52: double
! 53: ldexp(double val, int expon)
! 54: {
! 55: int oldexp, newexp;
! 56: union {
! 57: double v;
! 58: struct ieee_double s;
! 59: } u, mul;
! 60:
! 61: u.v = val;
! 62: oldexp = u.s.dbl_exp;
! 63:
! 64: /*
! 65: * If input is zero, Inf or NaN, just return it.
! 66: */
! 67: if (u.v == 0.0 || oldexp == DBL_EXP_INFNAN)
! 68: return (val);
! 69:
! 70: if (oldexp == 0) {
! 71: /*
! 72: * u.v is denormal. We must adjust it so that the exponent
! 73: * arithmetic below will work.
! 74: */
! 75: if (expon <= DBL_EXP_BIAS) {
! 76: /*
! 77: * Optimization: if the scaling can be done in a single
! 78: * multiply, or underflows, just do it now.
! 79: */
! 80: if (expon <= -DBL_FRACBITS) {
! 81: errno = ERANGE;
! 82: return (0.0);
! 83: }
! 84: mul.v = 0.0;
! 85: mul.s.dbl_exp = expon + DBL_EXP_BIAS;
! 86: u.v *= mul.v;
! 87: if (u.v == 0.0) {
! 88: errno = ERANGE;
! 89: return (0.0);
! 90: }
! 91: return (u.v);
! 92: } else {
! 93: /*
! 94: * We know that expon is very large, and therefore the
! 95: * result cannot be denormal (though it may be Inf).
! 96: * Shift u.v by just enough to make it normal.
! 97: */
! 98: mul.v = 0.0;
! 99: mul.s.dbl_exp = DBL_FRACBITS + DBL_EXP_BIAS;
! 100: u.v *= mul.v;
! 101: expon -= DBL_FRACBITS;
! 102: oldexp = u.s.dbl_exp;
! 103: }
! 104: }
! 105:
! 106: /*
! 107: * u.v is now normalized and oldexp has been adjusted if necessary.
! 108: * Calculate the new exponent and check for underflow and overflow.
! 109: */
! 110: newexp = oldexp + expon;
! 111:
! 112: if (newexp <= 0) {
! 113: /*
! 114: * The output number is either denormal or underflows (see
! 115: * comments in machine/ieee.h).
! 116: */
! 117: if (newexp <= -DBL_FRACBITS) {
! 118: errno = ERANGE;
! 119: return (0.0);
! 120: }
! 121: /*
! 122: * Denormalize the result. We do this with a multiply. If
! 123: * expon is very large, it won't fit in a double, so we have
! 124: * to adjust the exponent first. This is safe because we know
! 125: * that u.v is normal at this point.
! 126: */
! 127: if (expon <= -DBL_EXP_BIAS) {
! 128: u.s.dbl_exp = 1;
! 129: expon += oldexp - 1;
! 130: }
! 131: mul.v = 0.0;
! 132: mul.s.dbl_exp = expon + DBL_EXP_BIAS;
! 133: u.v *= mul.v;
! 134: return (u.v);
! 135: } else if (newexp >= DBL_EXP_INFNAN) {
! 136: /*
! 137: * The result overflowed; return +/-Inf.
! 138: */
! 139: u.s.dbl_exp = DBL_EXP_INFNAN;
! 140: u.s.dbl_frach = 0;
! 141: u.s.dbl_fracl = 0;
! 142: errno = ERANGE;
! 143: return (u.v);
! 144: } else {
! 145: /*
! 146: * The result is normal; just replace the old exponent with the
! 147: * new one.
! 148: */
! 149: u.s.dbl_exp = newexp;
! 150: return (u.v);
! 151: }
! 152: }
CVSweb <webmaster@jp.NetBSD.org>