1 | /* Compute x * y + z as ternary operation. |
2 | Copyright (C) 2010-2018 Free Software Foundation, Inc. |
3 | This file is part of the GNU C Library. |
4 | Contributed by Jakub Jelinek <jakub@redhat.com>, 2010. |
5 | |
6 | The GNU C Library is free software; you can redistribute it and/or |
7 | modify it under the terms of the GNU Lesser General Public |
8 | License as published by the Free Software Foundation; either |
9 | version 2.1 of the License, or (at your option) any later version. |
10 | |
11 | The GNU C Library is distributed in the hope that it will be useful, |
12 | but WITHOUT ANY WARRANTY; without even the implied warranty of |
13 | MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU |
14 | Lesser General Public License for more details. |
15 | |
16 | You should have received a copy of the GNU Lesser General Public |
17 | License along with the GNU C Library; if not, see |
18 | <http://www.gnu.org/licenses/>. */ |
19 | |
20 | #include <float.h> |
21 | #include <math.h> |
22 | #include <fenv.h> |
23 | #include <ieee754.h> |
24 | #include <math_private.h> |
25 | #include <libm-alias-ldouble.h> |
26 | #include <tininess.h> |
27 | |
28 | /* This implementation uses rounding to odd to avoid problems with |
29 | double rounding. See a paper by Boldo and Melquiond: |
30 | http://www.lri.fr/~melquion/doc/08-tc.pdf */ |
31 | |
32 | _Float128 |
33 | __fmal (_Float128 x, _Float128 y, _Float128 z) |
34 | { |
35 | union ieee854_long_double u, v, w; |
36 | int adjust = 0; |
37 | u.d = x; |
38 | v.d = y; |
39 | w.d = z; |
40 | if (__builtin_expect (u.ieee.exponent + v.ieee.exponent |
41 | >= 0x7fff + IEEE854_LONG_DOUBLE_BIAS |
42 | - LDBL_MANT_DIG, 0) |
43 | || __builtin_expect (u.ieee.exponent >= 0x7fff - LDBL_MANT_DIG, 0) |
44 | || __builtin_expect (v.ieee.exponent >= 0x7fff - LDBL_MANT_DIG, 0) |
45 | || __builtin_expect (w.ieee.exponent >= 0x7fff - LDBL_MANT_DIG, 0) |
46 | || __builtin_expect (u.ieee.exponent + v.ieee.exponent |
47 | <= IEEE854_LONG_DOUBLE_BIAS + LDBL_MANT_DIG, 0)) |
48 | { |
49 | /* If z is Inf, but x and y are finite, the result should be |
50 | z rather than NaN. */ |
51 | if (w.ieee.exponent == 0x7fff |
52 | && u.ieee.exponent != 0x7fff |
53 | && v.ieee.exponent != 0x7fff) |
54 | return (z + x) + y; |
55 | /* If z is zero and x are y are nonzero, compute the result |
56 | as x * y to avoid the wrong sign of a zero result if x * y |
57 | underflows to 0. */ |
58 | if (z == 0 && x != 0 && y != 0) |
59 | return x * y; |
60 | /* If x or y or z is Inf/NaN, or if x * y is zero, compute as |
61 | x * y + z. */ |
62 | if (u.ieee.exponent == 0x7fff |
63 | || v.ieee.exponent == 0x7fff |
64 | || w.ieee.exponent == 0x7fff |
65 | || x == 0 |
66 | || y == 0) |
67 | return x * y + z; |
68 | /* If fma will certainly overflow, compute as x * y. */ |
69 | if (u.ieee.exponent + v.ieee.exponent |
70 | > 0x7fff + IEEE854_LONG_DOUBLE_BIAS) |
71 | return x * y; |
72 | /* If x * y is less than 1/4 of LDBL_TRUE_MIN, neither the |
73 | result nor whether there is underflow depends on its exact |
74 | value, only on its sign. */ |
75 | if (u.ieee.exponent + v.ieee.exponent |
76 | < IEEE854_LONG_DOUBLE_BIAS - LDBL_MANT_DIG - 2) |
77 | { |
78 | int neg = u.ieee.negative ^ v.ieee.negative; |
79 | _Float128 tiny = neg ? L(-0x1p-16494) : L(0x1p-16494); |
80 | if (w.ieee.exponent >= 3) |
81 | return tiny + z; |
82 | /* Scaling up, adding TINY and scaling down produces the |
83 | correct result, because in round-to-nearest mode adding |
84 | TINY has no effect and in other modes double rounding is |
85 | harmless. But it may not produce required underflow |
86 | exceptions. */ |
87 | v.d = z * L(0x1p114) + tiny; |
88 | if (TININESS_AFTER_ROUNDING |
89 | ? v.ieee.exponent < 115 |
90 | : (w.ieee.exponent == 0 |
91 | || (w.ieee.exponent == 1 |
92 | && w.ieee.negative != neg |
93 | && w.ieee.mantissa3 == 0 |
94 | && w.ieee.mantissa2 == 0 |
95 | && w.ieee.mantissa1 == 0 |
96 | && w.ieee.mantissa0 == 0))) |
97 | { |
98 | _Float128 force_underflow = x * y; |
99 | math_force_eval (force_underflow); |
100 | } |
101 | return v.d * L(0x1p-114); |
102 | } |
103 | if (u.ieee.exponent + v.ieee.exponent |
104 | >= 0x7fff + IEEE854_LONG_DOUBLE_BIAS - LDBL_MANT_DIG) |
105 | { |
106 | /* Compute 1p-113 times smaller result and multiply |
107 | at the end. */ |
108 | if (u.ieee.exponent > v.ieee.exponent) |
109 | u.ieee.exponent -= LDBL_MANT_DIG; |
110 | else |
111 | v.ieee.exponent -= LDBL_MANT_DIG; |
112 | /* If x + y exponent is very large and z exponent is very small, |
113 | it doesn't matter if we don't adjust it. */ |
114 | if (w.ieee.exponent > LDBL_MANT_DIG) |
115 | w.ieee.exponent -= LDBL_MANT_DIG; |
116 | adjust = 1; |
117 | } |
118 | else if (w.ieee.exponent >= 0x7fff - LDBL_MANT_DIG) |
119 | { |
120 | /* Similarly. |
121 | If z exponent is very large and x and y exponents are |
122 | very small, adjust them up to avoid spurious underflows, |
123 | rather than down. */ |
124 | if (u.ieee.exponent + v.ieee.exponent |
125 | <= IEEE854_LONG_DOUBLE_BIAS + 2 * LDBL_MANT_DIG) |
126 | { |
127 | if (u.ieee.exponent > v.ieee.exponent) |
128 | u.ieee.exponent += 2 * LDBL_MANT_DIG + 2; |
129 | else |
130 | v.ieee.exponent += 2 * LDBL_MANT_DIG + 2; |
131 | } |
132 | else if (u.ieee.exponent > v.ieee.exponent) |
133 | { |
134 | if (u.ieee.exponent > LDBL_MANT_DIG) |
135 | u.ieee.exponent -= LDBL_MANT_DIG; |
136 | } |
137 | else if (v.ieee.exponent > LDBL_MANT_DIG) |
138 | v.ieee.exponent -= LDBL_MANT_DIG; |
139 | w.ieee.exponent -= LDBL_MANT_DIG; |
140 | adjust = 1; |
141 | } |
142 | else if (u.ieee.exponent >= 0x7fff - LDBL_MANT_DIG) |
143 | { |
144 | u.ieee.exponent -= LDBL_MANT_DIG; |
145 | if (v.ieee.exponent) |
146 | v.ieee.exponent += LDBL_MANT_DIG; |
147 | else |
148 | v.d *= L(0x1p113); |
149 | } |
150 | else if (v.ieee.exponent >= 0x7fff - LDBL_MANT_DIG) |
151 | { |
152 | v.ieee.exponent -= LDBL_MANT_DIG; |
153 | if (u.ieee.exponent) |
154 | u.ieee.exponent += LDBL_MANT_DIG; |
155 | else |
156 | u.d *= L(0x1p113); |
157 | } |
158 | else /* if (u.ieee.exponent + v.ieee.exponent |
159 | <= IEEE854_LONG_DOUBLE_BIAS + LDBL_MANT_DIG) */ |
160 | { |
161 | if (u.ieee.exponent > v.ieee.exponent) |
162 | u.ieee.exponent += 2 * LDBL_MANT_DIG + 2; |
163 | else |
164 | v.ieee.exponent += 2 * LDBL_MANT_DIG + 2; |
165 | if (w.ieee.exponent <= 4 * LDBL_MANT_DIG + 6) |
166 | { |
167 | if (w.ieee.exponent) |
168 | w.ieee.exponent += 2 * LDBL_MANT_DIG + 2; |
169 | else |
170 | w.d *= L(0x1p228); |
171 | adjust = -1; |
172 | } |
173 | /* Otherwise x * y should just affect inexact |
174 | and nothing else. */ |
175 | } |
176 | x = u.d; |
177 | y = v.d; |
178 | z = w.d; |
179 | } |
180 | |
181 | /* Ensure correct sign of exact 0 + 0. */ |
182 | if (__glibc_unlikely ((x == 0 || y == 0) && z == 0)) |
183 | { |
184 | x = math_opt_barrier (x); |
185 | return x * y + z; |
186 | } |
187 | |
188 | fenv_t env; |
189 | feholdexcept (&env); |
190 | fesetround (FE_TONEAREST); |
191 | |
192 | /* Multiplication m1 + m2 = x * y using Dekker's algorithm. */ |
193 | #define C ((1LL << (LDBL_MANT_DIG + 1) / 2) + 1) |
194 | _Float128 x1 = x * C; |
195 | _Float128 y1 = y * C; |
196 | _Float128 m1 = x * y; |
197 | x1 = (x - x1) + x1; |
198 | y1 = (y - y1) + y1; |
199 | _Float128 x2 = x - x1; |
200 | _Float128 y2 = y - y1; |
201 | _Float128 m2 = (((x1 * y1 - m1) + x1 * y2) + x2 * y1) + x2 * y2; |
202 | |
203 | /* Addition a1 + a2 = z + m1 using Knuth's algorithm. */ |
204 | _Float128 a1 = z + m1; |
205 | _Float128 t1 = a1 - z; |
206 | _Float128 t2 = a1 - t1; |
207 | t1 = m1 - t1; |
208 | t2 = z - t2; |
209 | _Float128 a2 = t1 + t2; |
210 | /* Ensure the arithmetic is not scheduled after feclearexcept call. */ |
211 | math_force_eval (m2); |
212 | math_force_eval (a2); |
213 | feclearexcept (FE_INEXACT); |
214 | |
215 | /* If the result is an exact zero, ensure it has the correct sign. */ |
216 | if (a1 == 0 && m2 == 0) |
217 | { |
218 | feupdateenv (&env); |
219 | /* Ensure that round-to-nearest value of z + m1 is not reused. */ |
220 | z = math_opt_barrier (z); |
221 | return z + m1; |
222 | } |
223 | |
224 | fesetround (FE_TOWARDZERO); |
225 | /* Perform m2 + a2 addition with round to odd. */ |
226 | u.d = a2 + m2; |
227 | |
228 | if (__glibc_likely (adjust == 0)) |
229 | { |
230 | if ((u.ieee.mantissa3 & 1) == 0 && u.ieee.exponent != 0x7fff) |
231 | u.ieee.mantissa3 |= fetestexcept (FE_INEXACT) != 0; |
232 | feupdateenv (&env); |
233 | /* Result is a1 + u.d. */ |
234 | return a1 + u.d; |
235 | } |
236 | else if (__glibc_likely (adjust > 0)) |
237 | { |
238 | if ((u.ieee.mantissa3 & 1) == 0 && u.ieee.exponent != 0x7fff) |
239 | u.ieee.mantissa3 |= fetestexcept (FE_INEXACT) != 0; |
240 | feupdateenv (&env); |
241 | /* Result is a1 + u.d, scaled up. */ |
242 | return (a1 + u.d) * L(0x1p113); |
243 | } |
244 | else |
245 | { |
246 | if ((u.ieee.mantissa3 & 1) == 0) |
247 | u.ieee.mantissa3 |= fetestexcept (FE_INEXACT) != 0; |
248 | v.d = a1 + u.d; |
249 | /* Ensure the addition is not scheduled after fetestexcept call. */ |
250 | math_force_eval (v.d); |
251 | int j = fetestexcept (FE_INEXACT) != 0; |
252 | feupdateenv (&env); |
253 | /* Ensure the following computations are performed in default rounding |
254 | mode instead of just reusing the round to zero computation. */ |
255 | asm volatile ("" : "=m" (u) : "m" (u)); |
256 | /* If a1 + u.d is exact, the only rounding happens during |
257 | scaling down. */ |
258 | if (j == 0) |
259 | return v.d * L(0x1p-228); |
260 | /* If result rounded to zero is not subnormal, no double |
261 | rounding will occur. */ |
262 | if (v.ieee.exponent > 228) |
263 | return (a1 + u.d) * L(0x1p-228); |
264 | /* If v.d * 0x1p-228L with round to zero is a subnormal above |
265 | or equal to LDBL_MIN / 2, then v.d * 0x1p-228L shifts mantissa |
266 | down just by 1 bit, which means v.ieee.mantissa3 |= j would |
267 | change the round bit, not sticky or guard bit. |
268 | v.d * 0x1p-228L never normalizes by shifting up, |
269 | so round bit plus sticky bit should be already enough |
270 | for proper rounding. */ |
271 | if (v.ieee.exponent == 228) |
272 | { |
273 | /* If the exponent would be in the normal range when |
274 | rounding to normal precision with unbounded exponent |
275 | range, the exact result is known and spurious underflows |
276 | must be avoided on systems detecting tininess after |
277 | rounding. */ |
278 | if (TININESS_AFTER_ROUNDING) |
279 | { |
280 | w.d = a1 + u.d; |
281 | if (w.ieee.exponent == 229) |
282 | return w.d * L(0x1p-228); |
283 | } |
284 | /* v.ieee.mantissa3 & 2 is LSB bit of the result before rounding, |
285 | v.ieee.mantissa3 & 1 is the round bit and j is our sticky |
286 | bit. */ |
287 | w.d = 0; |
288 | w.ieee.mantissa3 = ((v.ieee.mantissa3 & 3) << 1) | j; |
289 | w.ieee.negative = v.ieee.negative; |
290 | v.ieee.mantissa3 &= ~3U; |
291 | v.d *= L(0x1p-228); |
292 | w.d *= L(0x1p-2); |
293 | return v.d + w.d; |
294 | } |
295 | v.ieee.mantissa3 |= j; |
296 | return v.d * L(0x1p-228); |
297 | } |
298 | } |
299 | libm_alias_ldouble (__fma, fma) |
300 | |