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