s_fmal.c revision 226245 
1/*-
2 * Copyright (c) 2005-2011 David Schultz <das@FreeBSD.ORG>
3 * All rights reserved.
4 *
5 * Redistribution and use in source and binary forms, with or without
6 * modification, are permitted provided that the following conditions
7 * are met:
8 * 1. Redistributions of source code must retain the above copyright
9 *    notice, this list of conditions and the following disclaimer.
10 * 2. Redistributions in binary form must reproduce the above copyright
11 *    notice, this list of conditions and the following disclaimer in the
12 *    documentation and/or other materials provided with the distribution.
13 *
14 * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
15 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
16 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
17 * ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
18 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
19 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
20 * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
21 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
22 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
23 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
24 * SUCH DAMAGE.
25 */
26
27#include <sys/cdefs.h>
28__FBSDID("$FreeBSD: head/lib/msun/src/s_fmal.c 226245 2011-10-11 05:17:45Z das $");
29
30#include <fenv.h>
31#include <float.h>
32#include <math.h>
33
34/*
35 * A struct dd represents a floating-point number with twice the precision
36 * of a long double.  We maintain the invariant that "hi" stores the high-order
37 * bits of the result.
38 */
39struct dd {
40	long double hi;
41	long double lo;
42};
43
44/*
45 * Compute a+b exactly, returning the exact result in a struct dd.  We assume
46 * that both a and b are finite, but make no assumptions about their relative
47 * magnitudes.
48 */
49static inline struct dd
50dd_add(long double a, long double b)
51{
52	struct dd ret;
53	long double s;
54
55	ret.hi = a + b;
56	s = ret.hi - a;
57	ret.lo = (a - (ret.hi - s)) + (b - s);
58	return (ret);
59}
60
61/*
62 * Compute a*b exactly, returning the exact result in a struct dd.  We assume
63 * that both a and b are normalized, so no underflow or overflow will occur.
64 * The current rounding mode must be round-to-nearest.
65 */
66static inline struct dd
67dd_mul(long double a, long double b)
68{
69#if LDBL_MANT_DIG == 64
70	static const long double split = 0x1p32L + 1.0;
71#elif LDBL_MANT_DIG == 113
72	static const long double split = 0x1p57L + 1.0;
73#endif
74	struct dd ret;
75	long double ha, hb, la, lb, p, q;
76
77	p = a * split;
78	ha = a - p;
79	ha += p;
80	la = a - ha;
81
82	p = b * split;
83	hb = b - p;
84	hb += p;
85	lb = b - hb;
86
87	p = ha * hb;
88	q = ha * lb + la * hb;
89
90	ret.hi = p + q;
91	ret.lo = p - ret.hi + q + la * lb;
92	return (ret);
93}
94
95/*
96 * Fused multiply-add: Compute x * y + z with a single rounding error.
97 *
98 * We use scaling to avoid overflow/underflow, along with the
99 * canonical precision-doubling technique adapted from:
100 *
101 *	Dekker, T.  A Floating-Point Technique for Extending the
102 *	Available Precision.  Numer. Math. 18, 224-242 (1971).
103 */
104long double
105fmal(long double x, long double y, long double z)
106{
107	long double xs, ys, zs;
108	struct dd xy, r, r2;
109	long double p;
110	long double s;
111	int oround;
112	int ex, ey, ez;
113	int spread;
114
115	/*
116	 * Handle special cases. The order of operations and the particular
117	 * return values here are crucial in handling special cases involving
118	 * infinities, NaNs, overflows, and signed zeroes correctly.
119	 */
120	if (x == 0.0 || y == 0.0)
121		return (x * y + z);
122	if (z == 0.0)
123		return (x * y);
124	if (!isfinite(x) || !isfinite(y))
125		return (x * y + z);
126	if (!isfinite(z))
127		return (z);
128
129	xs = frexpl(x, &ex);
130	ys = frexpl(y, &ey);
131	zs = frexpl(z, &ez);
132	oround = fegetround();
133	spread = ex + ey - ez;
134
135	/*
136	 * If x * y and z are many orders of magnitude apart, the scaling
137	 * will overflow, so we handle these cases specially.  Rounding
138	 * modes other than FE_TONEAREST are painful.
139	 */
140	if (spread > LDBL_MANT_DIG * 2) {
141		fenv_t env;
142		feraiseexcept(FE_INEXACT);
143		switch(oround) {
144		case FE_TONEAREST:
145			return (x * y);
146		case FE_TOWARDZERO:
147			if (x > 0.0 ^ y < 0.0 ^ z < 0.0)
148				return (x * y);
149			feholdexcept(&env);
150			s = x * y;
151			if (!fetestexcept(FE_INEXACT))
152				s = nextafterl(s, 0);
153			feupdateenv(&env);
154			return (s);
155		case FE_DOWNWARD:
156			if (z > 0.0)
157				return (x * y);
158			feholdexcept(&env);
159			s = x * y;
160			if (!fetestexcept(FE_INEXACT))
161				s = nextafterl(s, -INFINITY);
162			feupdateenv(&env);
163			return (s);
164		default:	/* FE_UPWARD */
165			if (z < 0.0)
166				return (x * y);
167			feholdexcept(&env);
168			s = x * y;
169			if (!fetestexcept(FE_INEXACT))
170				s = nextafterl(s, INFINITY);
171			feupdateenv(&env);
172			return (s);
173		}
174	}
175	if (spread < -LDBL_MANT_DIG) {
176		feraiseexcept(FE_INEXACT);
177		if (!isnormal(z))
178			feraiseexcept(FE_UNDERFLOW);
179		switch (oround) {
180		case FE_TONEAREST:
181			return (z);
182		case FE_TOWARDZERO:
183			if (x > 0.0 ^ y < 0.0 ^ z < 0.0)
184				return (z);
185			else
186				return (nextafterl(z, 0));
187		case FE_DOWNWARD:
188			if (x > 0.0 ^ y < 0.0)
189				return (z);
190			else
191				return (nextafterl(z, -INFINITY));
192		default:	/* FE_UPWARD */
193			if (x > 0.0 ^ y < 0.0)
194				return (nextafterl(z, INFINITY));
195			else
196				return (z);
197		}
198	}
199
200	fesetround(FE_TONEAREST);
201
202	xy = dd_mul(xs, ys);
203	zs = ldexpl(zs, -spread);
204	r = dd_add(xy.hi, zs);
205	r.lo += xy.lo;
206
207	spread = ex + ey;
208	if (spread + ilogbl(r.hi) > -16383) {
209		fesetround(oround);
210		r.hi = r.hi + r.lo;
211	} else {
212		/*
213		 * The result is subnormal, so we round before scaling to
214		 * avoid double rounding.
215		 */
216		p = ldexpl(copysignl(0x1p-16382L, r.hi), -spread);
217		r2 = dd_add(r.hi, p);
218		r2.lo += r.lo;
219		fesetround(oround);
220		r.hi = (r2.hi + r2.lo) - p;
221	}
222	return (ldexpl(r.hi, spread));
223}
224