1/* $NetBSD: compat_ldexp_ieee754.c,v 1.4 2008/09/28 18:54:30 christos Exp $ */
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)
34__RCSID("$NetBSD: compat_ldexp_ieee754.c,v 1.4 2008/09/28 18:54:30 christos Exp $");
35#endif /* LIBC_SCCS and not lint */
36
37#include <sys/types.h>
38#include <machine/ieee.h>
39#include <errno.h>
40
41double ldexp(double, int);
42
43/*
44 * Multiply the given value by 2^expon.
45 */
46double
47ldexp(double val, int expon)
48{
49	int oldexp, newexp;
50	union ieee_double_u u, mul;
51
52	u.dblu_d = val;
53	oldexp = u.dblu_dbl.dbl_exp;
54
55	/*
56	 * If input is zero, Inf or NaN, just return it.
57	 */
58	if (u.dblu_d == 0.0 || oldexp == DBL_EXP_INFNAN)
59		return (val);
60
61	if (oldexp == 0) {
62		/*
63		 * u.v is denormal.  We must adjust it so that the exponent
64		 * arithmetic below will work.
65		 */
66		if (expon <= DBL_EXP_BIAS) {
67			/*
68			 * Optimization: if the scaling can be done in a single
69			 * multiply, or underflows, just do it now.
70			 */
71			if (expon <= -DBL_FRACBITS) {
72				errno = ERANGE;
73				return (val < 0.0 ? -0.0 : 0.0);
74			}
75			mul.dblu_d = 0.0;
76			mul.dblu_dbl.dbl_exp = expon + DBL_EXP_BIAS;
77			u.dblu_d *= mul.dblu_d;
78			if (u.dblu_d == 0.0) {
79				errno = ERANGE;
80				return (val < 0.0 ? -0.0 : 0.0);
81			}
82			return (u.dblu_d);
83		} else {
84			/*
85			 * We know that expon is very large, and therefore the
86			 * result cannot be denormal (though it may be Inf).
87			 * Shift u.v by just enough to make it normal.
88			 */
89			mul.dblu_d = 0.0;
90			mul.dblu_dbl.dbl_exp = DBL_FRACBITS + DBL_EXP_BIAS;
91			u.dblu_d *= mul.dblu_d;
92			expon -= DBL_FRACBITS;
93			oldexp = u.dblu_dbl.dbl_exp;
94		}
95	}
96
97	/*
98	 * u.v is now normalized and oldexp has been adjusted if necessary.
99	 * Calculate the new exponent and check for underflow and overflow.
100	 */
101	newexp = oldexp + expon;
102
103	if (newexp <= 0) {
104		/*
105		 * The output number is either denormal or underflows (see
106		 * comments in machine/ieee.h).
107		 */
108		if (newexp <= -DBL_FRACBITS) {
109			errno = ERANGE;
110			return (val < 0.0 ? -0.0 : 0.0);
111		}
112		/*
113		 * Denormalize the result.  We do this with a multiply.  If
114		 * expon is very large, it won't fit in a double, so we have
115		 * to adjust the exponent first.  This is safe because we know
116		 * that u.v is normal at this point.
117		 */
118		if (expon <= -DBL_EXP_BIAS) {
119			u.dblu_dbl.dbl_exp = 1;
120			expon += oldexp - 1;
121		}
122		mul.dblu_d = 0.0;
123		mul.dblu_dbl.dbl_exp = expon + DBL_EXP_BIAS;
124		u.dblu_d *= mul.dblu_d;
125		return (u.dblu_d);
126	} else if (newexp >= DBL_EXP_INFNAN) {
127		/*
128		 * The result overflowed; return +/-Inf.
129		 */
130		u.dblu_dbl.dbl_exp = DBL_EXP_INFNAN;
131		u.dblu_dbl.dbl_frach = 0;
132		u.dblu_dbl.dbl_fracl = 0;
133		errno = ERANGE;
134		return (u.dblu_d);
135	} else {
136		/*
137		 * The result is normal; just replace the old exponent with the
138		 * new one.
139		 */
140		u.dblu_dbl.dbl_exp = newexp;
141		return (u.dblu_d);
142	}
143}
144