]>
git.ipfire.org Git - thirdparty/glibc.git/blob - sysdeps/ieee754/dbl-64/s_fma.c
14c6503de27597f02bd14426af27a1da6a68c11b
1 /* Compute x * y + z as ternary operation.
2 Copyright (C) 2010, 2011 Free Software Foundation, Inc.
3 This file is part of the GNU C Library.
4 Contributed by Jakub Jelinek <jakub@redhat.com>, 2010.
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.
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.
16 You should have received a copy of the GNU Lesser General Public
17 License along with the GNU C Library; if not, write to the Free
18 Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
25 #include <math_private.h>
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 */
32 __fma (double x
, double y
, double z
)
34 union ieee754_double u
, v
, w
;
39 if (__builtin_expect (u
.ieee
.exponent
+ v
.ieee
.exponent
40 >= 0x7ff + IEEE754_DOUBLE_BIAS
- DBL_MANT_DIG
, 0)
41 || __builtin_expect (u
.ieee
.exponent
>= 0x7ff - DBL_MANT_DIG
, 0)
42 || __builtin_expect (v
.ieee
.exponent
>= 0x7ff - DBL_MANT_DIG
, 0)
43 || __builtin_expect (w
.ieee
.exponent
>= 0x7ff - DBL_MANT_DIG
, 0)
44 || __builtin_expect (u
.ieee
.exponent
+ v
.ieee
.exponent
45 <= IEEE754_DOUBLE_BIAS
+ DBL_MANT_DIG
, 0))
47 /* If z is Inf, but x and y are finite, the result should be
49 if (w
.ieee
.exponent
== 0x7ff
50 && u
.ieee
.exponent
!= 0x7ff
51 && v
.ieee
.exponent
!= 0x7ff)
53 /* If x or y or z is Inf/NaN, or if fma will certainly overflow,
54 or if x * y is less than half of DBL_DENORM_MIN,
55 compute as x * y + z. */
56 if (u
.ieee
.exponent
== 0x7ff
57 || v
.ieee
.exponent
== 0x7ff
58 || w
.ieee
.exponent
== 0x7ff
59 || u
.ieee
.exponent
+ v
.ieee
.exponent
60 > 0x7ff + IEEE754_DOUBLE_BIAS
61 || u
.ieee
.exponent
+ v
.ieee
.exponent
62 < IEEE754_DOUBLE_BIAS
- DBL_MANT_DIG
- 2)
64 if (u
.ieee
.exponent
+ v
.ieee
.exponent
65 >= 0x7ff + IEEE754_DOUBLE_BIAS
- DBL_MANT_DIG
)
67 /* Compute 1p-53 times smaller result and multiply
69 if (u
.ieee
.exponent
> v
.ieee
.exponent
)
70 u
.ieee
.exponent
-= DBL_MANT_DIG
;
72 v
.ieee
.exponent
-= DBL_MANT_DIG
;
73 /* If x + y exponent is very large and z exponent is very small,
74 it doesn't matter if we don't adjust it. */
75 if (w
.ieee
.exponent
> DBL_MANT_DIG
)
76 w
.ieee
.exponent
-= DBL_MANT_DIG
;
79 else if (w
.ieee
.exponent
>= 0x7ff - DBL_MANT_DIG
)
82 If z exponent is very large and x and y exponents are
83 very small, it doesn't matter if we don't adjust it. */
84 if (u
.ieee
.exponent
> v
.ieee
.exponent
)
86 if (u
.ieee
.exponent
> DBL_MANT_DIG
)
87 u
.ieee
.exponent
-= DBL_MANT_DIG
;
89 else if (v
.ieee
.exponent
> DBL_MANT_DIG
)
90 v
.ieee
.exponent
-= DBL_MANT_DIG
;
91 w
.ieee
.exponent
-= DBL_MANT_DIG
;
94 else if (u
.ieee
.exponent
>= 0x7ff - DBL_MANT_DIG
)
96 u
.ieee
.exponent
-= DBL_MANT_DIG
;
98 v
.ieee
.exponent
+= DBL_MANT_DIG
;
102 else if (v
.ieee
.exponent
>= 0x7ff - DBL_MANT_DIG
)
104 v
.ieee
.exponent
-= DBL_MANT_DIG
;
106 u
.ieee
.exponent
+= DBL_MANT_DIG
;
110 else /* if (u.ieee.exponent + v.ieee.exponent
111 <= IEEE754_DOUBLE_BIAS + DBL_MANT_DIG) */
113 if (u
.ieee
.exponent
> v
.ieee
.exponent
)
114 u
.ieee
.exponent
+= 2 * DBL_MANT_DIG
;
116 v
.ieee
.exponent
+= 2 * DBL_MANT_DIG
;
117 if (w
.ieee
.exponent
<= 4 * DBL_MANT_DIG
+ 4)
120 w
.ieee
.exponent
+= 2 * DBL_MANT_DIG
;
125 /* Otherwise x * y should just affect inexact
132 /* Multiplication m1 + m2 = x * y using Dekker's algorithm. */
133 #define C ((1 << (DBL_MANT_DIG + 1) / 2) + 1)
141 double m2
= (((x1
* y1
- m1
) + x1
* y2
) + x2
* y1
) + x2
* y2
;
143 /* Addition a1 + a2 = z + m1 using Knuth's algorithm. */
152 libc_feholdexcept_setround (&env
, FE_TOWARDZERO
);
153 /* Perform m2 + a2 addition with round to odd. */
156 if (__builtin_expect (adjust
== 0, 1))
158 if ((u
.ieee
.mantissa1
& 1) == 0 && u
.ieee
.exponent
!= 0x7ff)
159 u
.ieee
.mantissa1
|= libc_fetestexcept (FE_INEXACT
) != 0;
160 libc_feupdateenv (&env
);
161 /* Result is a1 + u.d. */
164 else if (__builtin_expect (adjust
> 0, 1))
166 if ((u
.ieee
.mantissa1
& 1) == 0 && u
.ieee
.exponent
!= 0x7ff)
167 u
.ieee
.mantissa1
|= libc_fetestexcept (FE_INEXACT
) != 0;
168 libc_feupdateenv (&env
);
169 /* Result is a1 + u.d, scaled up. */
170 return (a1
+ u
.d
) * 0x1p
53;
174 if ((u
.ieee
.mantissa1
& 1) == 0)
175 u
.ieee
.mantissa1
|= libc_fetestexcept (FE_INEXACT
) != 0;
177 int j
= libc_fetestexcept (FE_INEXACT
) != 0;
178 libc_feupdateenv (&env
);
179 /* Ensure the following computations are performed in default rounding
180 mode instead of just reusing the round to zero computation. */
181 asm volatile ("" : "=m" (u
) : "m" (u
));
182 /* If a1 + u.d is exact, the only rounding happens during
185 return v
.d
* 0x1p
-106;
186 /* If result rounded to zero is not subnormal, no double
187 rounding will occur. */
188 if (v
.ieee
.exponent
> 106)
189 return (a1
+ u
.d
) * 0x1p
-106;
190 /* If v.d * 0x1p-106 with round to zero is a subnormal above
191 or equal to DBL_MIN / 2, then v.d * 0x1p-106 shifts mantissa
192 down just by 1 bit, which means v.ieee.mantissa1 |= j would
193 change the round bit, not sticky or guard bit.
194 v.d * 0x1p-106 never normalizes by shifting up,
195 so round bit plus sticky bit should be already enough
196 for proper rounding. */
197 if (v
.ieee
.exponent
== 106)
199 /* v.ieee.mantissa1 & 2 is LSB bit of the result before rounding,
200 v.ieee.mantissa1 & 1 is the round bit and j is our sticky
201 bit. In round-to-nearest 001 rounds down like 00,
202 011 rounds up, even though 01 rounds down (thus we need
203 to adjust), 101 rounds down like 10 and 111 rounds up
205 if ((v
.ieee
.mantissa1
& 3) == 1)
209 return v
.d
- 0x1p
-1074 /* __DBL_DENORM_MIN__ */;
211 return v
.d
+ 0x1p
-1074 /* __DBL_DENORM_MIN__ */;
214 return v
.d
* 0x1p
-106;
216 v
.ieee
.mantissa1
|= j
;
217 return v
.d
* 0x1p
-106;
221 weak_alias (__fma
, fma
)
224 #ifdef NO_LONG_DOUBLE
225 strong_alias (__fma
, __fmal
)
226 weak_alias (__fmal
, fmal
)