Annotation of src/lib/libc/gen/ldexp_ieee754.c, Revision 1.4.30.1
1.4.30.1! yamt 1: /* $NetBSD: ldexp_ieee754.c,v 1.4 2003/12/30 23:19:07 martin Exp $ */
1.1 kleink 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: *
19: * THIS SOFTWARE IS PROVIDED BY THE NETBSD FOUNDATION, INC. AND CONTRIBUTORS
20: * ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
21: * TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
22: * PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE FOUNDATION OR CONTRIBUTORS
23: * BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
24: * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
25: * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
26: * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
27: * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
28: * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
29: * POSSIBILITY OF SUCH DAMAGE.
30: */
31:
32: #include <sys/cdefs.h>
33: #if defined(LIBC_SCCS) && !defined(lint)
1.4.30.1! yamt 34: __RCSID("$NetBSD: ldexp_ieee754.c,v 1.4 2003/12/30 23:19:07 martin Exp $");
1.1 kleink 35: #endif /* LIBC_SCCS and not lint */
36:
37: #include <sys/types.h>
38: #include <machine/ieee.h>
39: #include <errno.h>
40: #include <math.h>
41:
42: /*
43: * Multiply the given value by 2^expon.
44: */
45: double
46: ldexp(double val, int expon)
47: {
48: int oldexp, newexp;
1.2 kleink 49: union ieee_double_u u, mul;
1.1 kleink 50:
1.3 mrg 51: u.dblu_d = val;
1.2 kleink 52: oldexp = u.dblu_dbl.dbl_exp;
1.1 kleink 53:
54: /*
55: * If input is zero, Inf or NaN, just return it.
56: */
1.2 kleink 57: if (u.dblu_d == 0.0 || oldexp == DBL_EXP_INFNAN)
1.1 kleink 58: return (val);
59:
60: if (oldexp == 0) {
61: /*
62: * u.v is denormal. We must adjust it so that the exponent
63: * arithmetic below will work.
64: */
65: if (expon <= DBL_EXP_BIAS) {
66: /*
67: * Optimization: if the scaling can be done in a single
68: * multiply, or underflows, just do it now.
69: */
70: if (expon <= -DBL_FRACBITS) {
71: errno = ERANGE;
1.4 martin 72: return (val < 0.0 ? -0.0 : 0.0);
1.1 kleink 73: }
1.2 kleink 74: mul.dblu_d = 0.0;
75: mul.dblu_dbl.dbl_exp = expon + DBL_EXP_BIAS;
76: u.dblu_d *= mul.dblu_d;
77: if (u.dblu_d == 0.0) {
1.1 kleink 78: errno = ERANGE;
1.4 martin 79: return (val < 0.0 ? -0.0 : 0.0);
1.1 kleink 80: }
1.2 kleink 81: return (u.dblu_d);
1.1 kleink 82: } else {
83: /*
84: * We know that expon is very large, and therefore the
85: * result cannot be denormal (though it may be Inf).
86: * Shift u.v by just enough to make it normal.
87: */
1.2 kleink 88: mul.dblu_d = 0.0;
89: mul.dblu_dbl.dbl_exp = DBL_FRACBITS + DBL_EXP_BIAS;
90: u.dblu_d *= mul.dblu_d;
1.1 kleink 91: expon -= DBL_FRACBITS;
1.2 kleink 92: oldexp = u.dblu_dbl.dbl_exp;
1.1 kleink 93: }
94: }
95:
96: /*
97: * u.v is now normalized and oldexp has been adjusted if necessary.
98: * Calculate the new exponent and check for underflow and overflow.
99: */
100: newexp = oldexp + expon;
101:
102: if (newexp <= 0) {
103: /*
104: * The output number is either denormal or underflows (see
105: * comments in machine/ieee.h).
106: */
107: if (newexp <= -DBL_FRACBITS) {
108: errno = ERANGE;
1.4 martin 109: return (val < 0.0 ? -0.0 : 0.0);
1.1 kleink 110: }
111: /*
112: * Denormalize the result. We do this with a multiply. If
113: * expon is very large, it won't fit in a double, so we have
114: * to adjust the exponent first. This is safe because we know
115: * that u.v is normal at this point.
116: */
117: if (expon <= -DBL_EXP_BIAS) {
1.2 kleink 118: u.dblu_dbl.dbl_exp = 1;
1.1 kleink 119: expon += oldexp - 1;
120: }
1.2 kleink 121: mul.dblu_d = 0.0;
122: mul.dblu_dbl.dbl_exp = expon + DBL_EXP_BIAS;
123: u.dblu_d *= mul.dblu_d;
124: return (u.dblu_d);
1.1 kleink 125: } else if (newexp >= DBL_EXP_INFNAN) {
126: /*
127: * The result overflowed; return +/-Inf.
128: */
1.2 kleink 129: u.dblu_dbl.dbl_exp = DBL_EXP_INFNAN;
130: u.dblu_dbl.dbl_frach = 0;
131: u.dblu_dbl.dbl_fracl = 0;
1.1 kleink 132: errno = ERANGE;
1.2 kleink 133: return (u.dblu_d);
1.1 kleink 134: } else {
135: /*
136: * The result is normal; just replace the old exponent with the
137: * new one.
138: */
1.2 kleink 139: u.dblu_dbl.dbl_exp = newexp;
140: return (u.dblu_d);
1.1 kleink 141: }
142: }
CVSweb <webmaster@jp.NetBSD.org>