Annotation of src/lib/libc/gen/ldexp_ieee754.c, Revision 1.2
1.2 ! kleink 1: /* $NetBSD: ldexp_ieee754.c,v 1.1 2003/05/12 15:15:16 kleink 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: * 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)
1.2 ! kleink 41: __RCSID("$NetBSD: ldexp_ieee754.c,v 1.1 2003/05/12 15:15:16 kleink Exp $");
1.1 kleink 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;
1.2 ! kleink 56: union ieee_double_u u, mul;
1.1 kleink 57:
1.2 ! kleink 58: u.dlbu_d = val;
! 59: oldexp = u.dblu_dbl.dbl_exp;
1.1 kleink 60:
61: /*
62: * If input is zero, Inf or NaN, just return it.
63: */
1.2 ! kleink 64: if (u.dblu_d == 0.0 || oldexp == DBL_EXP_INFNAN)
1.1 kleink 65: return (val);
66:
67: if (oldexp == 0) {
68: /*
69: * u.v is denormal. We must adjust it so that the exponent
70: * arithmetic below will work.
71: */
72: if (expon <= DBL_EXP_BIAS) {
73: /*
74: * Optimization: if the scaling can be done in a single
75: * multiply, or underflows, just do it now.
76: */
77: if (expon <= -DBL_FRACBITS) {
78: errno = ERANGE;
79: return (0.0);
80: }
1.2 ! kleink 81: mul.dblu_d = 0.0;
! 82: mul.dblu_dbl.dbl_exp = expon + DBL_EXP_BIAS;
! 83: u.dblu_d *= mul.dblu_d;
! 84: if (u.dblu_d == 0.0) {
1.1 kleink 85: errno = ERANGE;
86: return (0.0);
87: }
1.2 ! kleink 88: return (u.dblu_d);
1.1 kleink 89: } else {
90: /*
91: * We know that expon is very large, and therefore the
92: * result cannot be denormal (though it may be Inf).
93: * Shift u.v by just enough to make it normal.
94: */
1.2 ! kleink 95: mul.dblu_d = 0.0;
! 96: mul.dblu_dbl.dbl_exp = DBL_FRACBITS + DBL_EXP_BIAS;
! 97: u.dblu_d *= mul.dblu_d;
1.1 kleink 98: expon -= DBL_FRACBITS;
1.2 ! kleink 99: oldexp = u.dblu_dbl.dbl_exp;
1.1 kleink 100: }
101: }
102:
103: /*
104: * u.v is now normalized and oldexp has been adjusted if necessary.
105: * Calculate the new exponent and check for underflow and overflow.
106: */
107: newexp = oldexp + expon;
108:
109: if (newexp <= 0) {
110: /*
111: * The output number is either denormal or underflows (see
112: * comments in machine/ieee.h).
113: */
114: if (newexp <= -DBL_FRACBITS) {
115: errno = ERANGE;
116: return (0.0);
117: }
118: /*
119: * Denormalize the result. We do this with a multiply. If
120: * expon is very large, it won't fit in a double, so we have
121: * to adjust the exponent first. This is safe because we know
122: * that u.v is normal at this point.
123: */
124: if (expon <= -DBL_EXP_BIAS) {
1.2 ! kleink 125: u.dblu_dbl.dbl_exp = 1;
1.1 kleink 126: expon += oldexp - 1;
127: }
1.2 ! kleink 128: mul.dblu_d = 0.0;
! 129: mul.dblu_dbl.dbl_exp = expon + DBL_EXP_BIAS;
! 130: u.dblu_d *= mul.dblu_d;
! 131: return (u.dblu_d);
1.1 kleink 132: } else if (newexp >= DBL_EXP_INFNAN) {
133: /*
134: * The result overflowed; return +/-Inf.
135: */
1.2 ! kleink 136: u.dblu_dbl.dbl_exp = DBL_EXP_INFNAN;
! 137: u.dblu_dbl.dbl_frach = 0;
! 138: u.dblu_dbl.dbl_fracl = 0;
1.1 kleink 139: errno = ERANGE;
1.2 ! kleink 140: return (u.dblu_d);
1.1 kleink 141: } else {
142: /*
143: * The result is normal; just replace the old exponent with the
144: * new one.
145: */
1.2 ! kleink 146: u.dblu_dbl.dbl_exp = newexp;
! 147: return (u.dblu_d);
1.1 kleink 148: }
149: }
CVSweb <webmaster@jp.NetBSD.org>