]> git.ipfire.org Git - thirdparty/glibc.git/blame - sysdeps/ieee754/ldbl-96/s_fmal.c
Update copyright notices with scripts/update-copyrights
[thirdparty/glibc.git] / sysdeps / ieee754 / ldbl-96 / s_fmal.c
CommitLineData
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
31long 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-16445L : 0x1p-16445L;
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 * 0x1p65L + tiny;
87 if (TININESS_AFTER_ROUNDING
88 ? v.ieee.exponent < 66
89 : (w.ieee.exponent == 0
90 || (w.ieee.exponent == 1
91 && w.ieee.negative != neg
92 && w.ieee.mantissa1 == 0
93 && w.ieee.mantissa0 == 0x80000000)))
94 {
95 volatile long double force_underflow = x * y;
96 (void) force_underflow;
97 }
98 return v.d * 0x1p-65L;
99 }
3e692e05
JJ
100 if (u.ieee.exponent + v.ieee.exponent
101 >= 0x7fff + IEEE854_LONG_DOUBLE_BIAS - LDBL_MANT_DIG)
102 {
103 /* Compute 1p-64 times smaller result and multiply
104 at the end. */
105 if (u.ieee.exponent > v.ieee.exponent)
106 u.ieee.exponent -= LDBL_MANT_DIG;
107 else
108 v.ieee.exponent -= LDBL_MANT_DIG;
109 /* If x + y exponent is very large and z exponent is very small,
110 it doesn't matter if we don't adjust it. */
111 if (w.ieee.exponent > LDBL_MANT_DIG)
112 w.ieee.exponent -= LDBL_MANT_DIG;
113 adjust = 1;
114 }
115 else if (w.ieee.exponent >= 0x7fff - LDBL_MANT_DIG)
116 {
117 /* Similarly.
118 If z exponent is very large and x and y exponents are
82477c28
JM
119 very small, adjust them up to avoid spurious underflows,
120 rather than down. */
121 if (u.ieee.exponent + v.ieee.exponent
122 <= IEEE854_LONG_DOUBLE_BIAS + LDBL_MANT_DIG)
123 {
124 if (u.ieee.exponent > v.ieee.exponent)
125 u.ieee.exponent += 2 * LDBL_MANT_DIG + 2;
126 else
127 v.ieee.exponent += 2 * LDBL_MANT_DIG + 2;
128 }
129 else if (u.ieee.exponent > v.ieee.exponent)
3e692e05
JJ
130 {
131 if (u.ieee.exponent > LDBL_MANT_DIG)
132 u.ieee.exponent -= LDBL_MANT_DIG;
133 }
134 else if (v.ieee.exponent > LDBL_MANT_DIG)
135 v.ieee.exponent -= LDBL_MANT_DIG;
136 w.ieee.exponent -= LDBL_MANT_DIG;
137 adjust = 1;
138 }
139 else if (u.ieee.exponent >= 0x7fff - LDBL_MANT_DIG)
140 {
141 u.ieee.exponent -= LDBL_MANT_DIG;
142 if (v.ieee.exponent)
143 v.ieee.exponent += LDBL_MANT_DIG;
144 else
145 v.d *= 0x1p64L;
146 }
147 else if (v.ieee.exponent >= 0x7fff - LDBL_MANT_DIG)
148 {
149 v.ieee.exponent -= LDBL_MANT_DIG;
150 if (u.ieee.exponent)
151 u.ieee.exponent += LDBL_MANT_DIG;
152 else
153 u.d *= 0x1p64L;
154 }
155 else /* if (u.ieee.exponent + v.ieee.exponent
156 <= IEEE854_LONG_DOUBLE_BIAS + LDBL_MANT_DIG) */
157 {
158 if (u.ieee.exponent > v.ieee.exponent)
82477c28 159 u.ieee.exponent += 2 * LDBL_MANT_DIG + 2;
3e692e05 160 else
82477c28
JM
161 v.ieee.exponent += 2 * LDBL_MANT_DIG + 2;
162 if (w.ieee.exponent <= 4 * LDBL_MANT_DIG + 6)
3e692e05
JJ
163 {
164 if (w.ieee.exponent)
82477c28 165 w.ieee.exponent += 2 * LDBL_MANT_DIG + 2;
3e692e05 166 else
82477c28 167 w.d *= 0x1p130L;
3e692e05
JJ
168 adjust = -1;
169 }
170 /* Otherwise x * y should just affect inexact
171 and nothing else. */
172 }
173 x = u.d;
174 y = v.d;
175 z = w.d;
176 }
8ec5b013
JM
177
178 /* Ensure correct sign of exact 0 + 0. */
179 if (__builtin_expect ((x == 0 || y == 0) && z == 0, 0))
180 return x * y + z;
181
5b5b04d6
JM
182 fenv_t env;
183 feholdexcept (&env);
184 fesetround (FE_TONEAREST);
185
3e692e05
JJ
186 /* Multiplication m1 + m2 = x * y using Dekker's algorithm. */
187#define C ((1LL << (LDBL_MANT_DIG + 1) / 2) + 1)
188 long double x1 = x * C;
189 long double y1 = y * C;
190 long double m1 = x * y;
191 x1 = (x - x1) + x1;
192 y1 = (y - y1) + y1;
193 long double x2 = x - x1;
194 long double y2 = y - y1;
195 long double m2 = (((x1 * y1 - m1) + x1 * y2) + x2 * y1) + x2 * y2;
196
197 /* Addition a1 + a2 = z + m1 using Knuth's algorithm. */
198 long double a1 = z + m1;
199 long double t1 = a1 - z;
200 long double t2 = a1 - t1;
201 t1 = m1 - t1;
202 t2 = z - t2;
203 long double a2 = t1 + t2;
5b5b04d6
JM
204 feclearexcept (FE_INEXACT);
205
206 /* If the result is an exact zero, ensure it has the correct
207 sign. */
208 if (a1 == 0 && m2 == 0)
209 {
210 feupdateenv (&env);
211 /* Ensure that round-to-nearest value of z + m1 is not
212 reused. */
213 asm volatile ("" : "=m" (z) : "m" (z));
214 return z + m1;
215 }
3e692e05 216
3e692e05
JJ
217 fesetround (FE_TOWARDZERO);
218 /* Perform m2 + a2 addition with round to odd. */
219 u.d = a2 + m2;
220
221 if (__builtin_expect (adjust == 0, 1))
222 {
223 if ((u.ieee.mantissa1 & 1) == 0 && u.ieee.exponent != 0x7fff)
224 u.ieee.mantissa1 |= fetestexcept (FE_INEXACT) != 0;
225 feupdateenv (&env);
226 /* Result is a1 + u.d. */
227 return a1 + u.d;
228 }
229 else if (__builtin_expect (adjust > 0, 1))
230 {
231 if ((u.ieee.mantissa1 & 1) == 0 && u.ieee.exponent != 0x7fff)
232 u.ieee.mantissa1 |= fetestexcept (FE_INEXACT) != 0;
233 feupdateenv (&env);
234 /* Result is a1 + u.d, scaled up. */
235 return (a1 + u.d) * 0x1p64L;
236 }
237 else
238 {
239 if ((u.ieee.mantissa1 & 1) == 0)
240 u.ieee.mantissa1 |= fetestexcept (FE_INEXACT) != 0;
241 v.d = a1 + u.d;
4842e4fe
JM
242 /* Ensure the addition is not scheduled after fetestexcept call. */
243 math_force_eval (v.d);
3e692e05
JJ
244 int j = fetestexcept (FE_INEXACT) != 0;
245 feupdateenv (&env);
246 /* Ensure the following computations are performed in default rounding
247 mode instead of just reusing the round to zero computation. */
248 asm volatile ("" : "=m" (u) : "m" (u));
249 /* If a1 + u.d is exact, the only rounding happens during
250 scaling down. */
251 if (j == 0)
82477c28 252 return v.d * 0x1p-130L;
3e692e05
JJ
253 /* If result rounded to zero is not subnormal, no double
254 rounding will occur. */
82477c28
JM
255 if (v.ieee.exponent > 130)
256 return (a1 + u.d) * 0x1p-130L;
257 /* If v.d * 0x1p-130L with round to zero is a subnormal above
258 or equal to LDBL_MIN / 2, then v.d * 0x1p-130L shifts mantissa
3e692e05
JJ
259 down just by 1 bit, which means v.ieee.mantissa1 |= j would
260 change the round bit, not sticky or guard bit.
82477c28 261 v.d * 0x1p-130L never normalizes by shifting up,
3e692e05
JJ
262 so round bit plus sticky bit should be already enough
263 for proper rounding. */
82477c28 264 if (v.ieee.exponent == 130)
3e692e05 265 {
ef82f4da
JM
266 /* If the exponent would be in the normal range when
267 rounding to normal precision with unbounded exponent
268 range, the exact result is known and spurious underflows
269 must be avoided on systems detecting tininess after
270 rounding. */
271 if (TININESS_AFTER_ROUNDING)
272 {
273 w.d = a1 + u.d;
82477c28
JM
274 if (w.ieee.exponent == 131)
275 return w.d * 0x1p-130L;
ef82f4da 276 }
3e692e05
JJ
277 /* v.ieee.mantissa1 & 2 is LSB bit of the result before rounding,
278 v.ieee.mantissa1 & 1 is the round bit and j is our sticky
8627a232
JM
279 bit. */
280 w.d = 0.0L;
281 w.ieee.mantissa1 = ((v.ieee.mantissa1 & 3) << 1) | j;
282 w.ieee.negative = v.ieee.negative;
283 v.ieee.mantissa1 &= ~3U;
82477c28 284 v.d *= 0x1p-130L;
8627a232
JM
285 w.d *= 0x1p-2L;
286 return v.d + w.d;
3e692e05
JJ
287 }
288 v.ieee.mantissa1 |= j;
82477c28 289 return v.d * 0x1p-130L;
3e692e05
JJ
290 }
291}
292weak_alias (__fmal, fmal)