Branch data Line data Source code
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 "cdefs-compat.h"
28 : : //__FBSDID("$FreeBSD: src/lib/msun/src/s_fmal.c,v 1.7 2011/10/21 06:30:43 das Exp $");
29 : :
30 : : #include <float.h>
31 : : #include <openlibm_fenv.h>
32 : : #include <openlibm_math.h>
33 : :
34 : : #include "fpmath.h"
35 : : #include "math_private.h"
36 : :
37 : : /*
38 : : * A struct dd represents a floating-point number with twice the precision
39 : : * of a long double. We maintain the invariant that "hi" stores the high-order
40 : : * bits of the result.
41 : : */
42 : : struct dd {
43 : : long double hi;
44 : : long double lo;
45 : : };
46 : :
47 : : /*
48 : : * Compute a+b exactly, returning the exact result in a struct dd. We assume
49 : : * that both a and b are finite, but make no assumptions about their relative
50 : : * magnitudes.
51 : : */
52 : : static inline struct dd
53 : 0 : dd_add(long double a, long double b)
54 : : {
55 : : struct dd ret;
56 : : long double s;
57 : :
58 : 0 : ret.hi = a + b;
59 : 0 : s = ret.hi - a;
60 : 0 : ret.lo = (a - (ret.hi - s)) + (b - s);
61 : 0 : return (ret);
62 : : }
63 : :
64 : : /*
65 : : * Compute a+b, with a small tweak: The least significant bit of the
66 : : * result is adjusted into a sticky bit summarizing all the bits that
67 : : * were lost to rounding. This adjustment negates the effects of double
68 : : * rounding when the result is added to another number with a higher
69 : : * exponent. For an explanation of round and sticky bits, see any reference
70 : : * on FPU design, e.g.,
71 : : *
72 : : * J. Coonen. An Implementation Guide to a Proposed Standard for
73 : : * Floating-Point Arithmetic. Computer, vol. 13, no. 1, Jan 1980.
74 : : */
75 : : static inline long double
76 : 0 : add_adjusted(long double a, long double b)
77 : : {
78 : : struct dd sum;
79 : : union IEEEl2bits u;
80 : :
81 : 0 : sum = dd_add(a, b);
82 [ # # ]: 0 : if (sum.lo != 0) {
83 : 0 : u.e = sum.hi;
84 [ # # ]: 0 : if ((u.bits.manl & 1) == 0)
85 : 0 : sum.hi = nextafterl(sum.hi, INFINITY * sum.lo);
86 : : }
87 : 0 : return (sum.hi);
88 : : }
89 : :
90 : : /*
91 : : * Compute ldexp(a+b, scale) with a single rounding error. It is assumed
92 : : * that the result will be subnormal, and care is taken to ensure that
93 : : * double rounding does not occur.
94 : : */
95 : : static inline long double
96 : 0 : add_and_denormalize(long double a, long double b, int scale)
97 : : {
98 : : struct dd sum;
99 : : int bits_lost;
100 : : union IEEEl2bits u;
101 : :
102 : 0 : sum = dd_add(a, b);
103 : :
104 : : /*
105 : : * If we are losing at least two bits of accuracy to denormalization,
106 : : * then the first lost bit becomes a round bit, and we adjust the
107 : : * lowest bit of sum.hi to make it a sticky bit summarizing all the
108 : : * bits in sum.lo. With the sticky bit adjusted, the hardware will
109 : : * break any ties in the correct direction.
110 : : *
111 : : * If we are losing only one bit to denormalization, however, we must
112 : : * break the ties manually.
113 : : */
114 [ # # ]: 0 : if (sum.lo != 0) {
115 : 0 : u.e = sum.hi;
116 : 0 : bits_lost = -u.bits.exp - scale + 1;
117 [ # # ]: 0 : if ((bits_lost != 1) ^ (int)(u.bits.manl & 1))
118 : 0 : sum.hi = nextafterl(sum.hi, INFINITY * sum.lo);
119 : : }
120 : 0 : return (ldexp(sum.hi, scale));
121 : : }
122 : :
123 : : /*
124 : : * Compute a*b exactly, returning the exact result in a struct dd. We assume
125 : : * that both a and b are normalized, so no underflow or overflow will occur.
126 : : * The current rounding mode must be round-to-nearest.
127 : : */
128 : : static inline struct dd
129 : 0 : dd_mul(long double a, long double b)
130 : : {
131 : : #if LDBL_MANT_DIG == 64
132 : : static const long double split = 0x1p32L + 1.0;
133 : : #elif LDBL_MANT_DIG == 113
134 : : static const long double split = 0x1p57L + 1.0;
135 : : #endif
136 : : struct dd ret;
137 : : long double ha, hb, la, lb, p, q;
138 : :
139 : 0 : p = a * split;
140 : 0 : ha = a - p;
141 : 0 : ha += p;
142 : 0 : la = a - ha;
143 : :
144 : 0 : p = b * split;
145 : 0 : hb = b - p;
146 : 0 : hb += p;
147 : 0 : lb = b - hb;
148 : :
149 : 0 : p = ha * hb;
150 : 0 : q = ha * lb + la * hb;
151 : :
152 : 0 : ret.hi = p + q;
153 : 0 : ret.lo = p - ret.hi + q + la * lb;
154 : 0 : return (ret);
155 : : }
156 : :
157 : : /*
158 : : * Fused multiply-add: Compute x * y + z with a single rounding error.
159 : : *
160 : : * We use scaling to avoid overflow/underflow, along with the
161 : : * canonical precision-doubling technique adapted from:
162 : : *
163 : : * Dekker, T. A Floating-Point Technique for Extending the
164 : : * Available Precision. Numer. Math. 18, 224-242 (1971).
165 : : */
166 : : OLM_DLLEXPORT long double
167 : 0 : fmal(long double x, long double y, long double z)
168 : : {
169 : : long double xs, ys, zs, adj;
170 : : struct dd xy, r;
171 : : int oround;
172 : : int ex, ey, ez;
173 : : int spread;
174 : :
175 : : /*
176 : : * Handle special cases. The order of operations and the particular
177 : : * return values here are crucial in handling special cases involving
178 : : * infinities, NaNs, overflows, and signed zeroes correctly.
179 : : */
180 [ # # # # ]: 0 : if (x == 0.0 || y == 0.0)
181 : 0 : return (x * y + z);
182 [ # # ]: 0 : if (z == 0.0)
183 : 0 : return (x * y);
184 [ # # # # ]: 0 : if (!isfinite(x) || !isfinite(y))
185 : 0 : return (x * y + z);
186 [ # # ]: 0 : if (!isfinite(z))
187 : 0 : return (z);
188 : :
189 : 0 : xs = frexpl(x, &ex);
190 : 0 : ys = frexpl(y, &ey);
191 : 0 : zs = frexpl(z, &ez);
192 : 0 : oround = fegetround();
193 : 0 : spread = ex + ey - ez;
194 : :
195 : : /*
196 : : * If x * y and z are many orders of magnitude apart, the scaling
197 : : * will overflow, so we handle these cases specially. Rounding
198 : : * modes other than FE_TONEAREST are painful.
199 : : */
200 [ # # ]: 0 : if (spread < -LDBL_MANT_DIG) {
201 : 0 : feraiseexcept(FE_INEXACT);
202 [ # # ]: 0 : if (!isnormal(z))
203 : 0 : feraiseexcept(FE_UNDERFLOW);
204 [ # # # # ]: 0 : switch (oround) {
205 : 0 : case FE_TONEAREST:
206 : 0 : return (z);
207 : 0 : case FE_TOWARDZERO:
208 [ # # ]: 0 : if ((x > 0.0) ^ (y < 0.0) ^ (z < 0.0))
209 : 0 : return (z);
210 : : else
211 : 0 : return (nextafterl(z, 0));
212 : 0 : case FE_DOWNWARD:
213 [ # # ]: 0 : if ((x > 0.0) ^ (y < 0.0))
214 : 0 : return (z);
215 : : else
216 : 0 : return (nextafterl(z, -INFINITY));
217 : 0 : default: /* FE_UPWARD */
218 [ # # ]: 0 : if ((x > 0.0) ^ (y < 0.0))
219 : 0 : return (nextafterl(z, INFINITY));
220 : : else
221 : 0 : return (z);
222 : : }
223 : : }
224 [ # # ]: 0 : if (spread <= LDBL_MANT_DIG * 2)
225 : 0 : zs = ldexpl(zs, -spread);
226 : : else
227 : 0 : zs = copysignl(LDBL_MIN, zs);
228 : :
229 : 0 : fesetround(FE_TONEAREST);
230 : :
231 : : /*
232 : : * Basic approach for round-to-nearest:
233 : : *
234 : : * (xy.hi, xy.lo) = x * y (exact)
235 : : * (r.hi, r.lo) = xy.hi + z (exact)
236 : : * adj = xy.lo + r.lo (inexact; low bit is sticky)
237 : : * result = r.hi + adj (correctly rounded)
238 : : */
239 : 0 : xy = dd_mul(xs, ys);
240 : 0 : r = dd_add(xy.hi, zs);
241 : :
242 : 0 : spread = ex + ey;
243 : :
244 [ # # ]: 0 : if (r.hi == 0.0) {
245 : : /*
246 : : * When the addends cancel to 0, ensure that the result has
247 : : * the correct sign.
248 : : */
249 : 0 : fesetround(oround);
250 : 0 : volatile long double vzs = zs; /* XXX gcc CSE bug workaround */
251 : 0 : return (xy.hi + vzs + ldexpl(xy.lo, spread));
252 : : }
253 : :
254 [ # # ]: 0 : if (oround != FE_TONEAREST) {
255 : : /*
256 : : * There is no need to worry about double rounding in directed
257 : : * rounding modes.
258 : : */
259 : 0 : fesetround(oround);
260 : 0 : adj = r.lo + xy.lo;
261 : 0 : return (ldexpl(r.hi + adj, spread));
262 : : }
263 : :
264 : 0 : adj = add_adjusted(r.lo, xy.lo);
265 [ # # ]: 0 : if (spread + ilogbl(r.hi) > -16383)
266 : 0 : return (ldexpl(r.hi + adj, spread));
267 : : else
268 : 0 : return (add_and_denormalize(r.hi, adj, spread));
269 : : }
|
