1 /* Compute x * y + z as ternary operation.
2 Copyright (C) 2010-2022 Free Software Foundation, Inc.
3 This file is part of the GNU C Library.
4
5 The GNU C Library is free software; you can redistribute it and/or
6 modify it under the terms of the GNU Lesser General Public
7 License as published by the Free Software Foundation; either
8 version 2.1 of the License, or (at your option) any later version.
9
10 The GNU C Library is distributed in the hope that it will be useful,
11 but WITHOUT ANY WARRANTY; without even the implied warranty of
12 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
13 Lesser General Public License for more details.
14
15 You should have received a copy of the GNU Lesser General Public
16 License along with the GNU C Library; if not, see
17 <https://www.gnu.org/licenses/>. */
18
19 #define NO_MATH_REDIRECT
20 #define dfmal __hide_dfmal
21 #define f32xfmaf64 __hide_f32xfmaf64
22 #include <math.h>
23 #undef dfmal
24 #undef f32xfmaf64
25 #include <fenv.h>
26 #include <ieee754.h>
27 #include <libm-alias-double.h>
28 #include <math-narrow-alias.h>
29 #include <math-use-builtins.h>
30
31 /* This implementation relies on long double being more than twice as
32 precise as double and uses rounding to odd in order to avoid problems
33 with double rounding.
34 See a paper by Boldo and Melquiond:
35 http://www.lri.fr/~melquion/doc/08-tc.pdf */
36
37 double
__fma(double x,double y,double z)38 __fma (double x, double y, double z)
39 {
40 #if USE_FMA_BUILTIN
41 return __builtin_fma (x, y, z);
42 #else
43 fenv_t env;
44 /* Multiplication is always exact. */
45 long double temp = (long double) x * (long double) y;
46
47 /* Ensure correct sign of an exact zero result by performing the
48 addition in the original rounding mode in that case. */
49 if (temp == -z)
50 return (double) temp + z;
51
52 union ieee854_long_double u;
53 feholdexcept (&env);
54 fesetround (FE_TOWARDZERO);
55 /* Perform addition with round to odd. */
56 u.d = temp + (long double) z;
57 if ((u.ieee.mantissa3 & 1) == 0 && u.ieee.exponent != 0x7fff)
58 u.ieee.mantissa3 |= fetestexcept (FE_INEXACT) != 0;
59 feupdateenv (&env);
60 /* And finally truncation with round to nearest. */
61 return (double) u.d;
62 #endif /* ! USE_FMA_BUILTIN */
63 }
64 #ifndef __fma
65 libm_alias_double (__fma, fma)
66 libm_alias_double_narrow (__fma, fma)
67 #endif
68