1/* Compute x * y + z as ternary operation.
2 Copyright (C) 2010-2024 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#include <float.h>
21#define dfmal __hide_dfmal
22#define f32xfmaf64 __hide_f32xfmaf64
23#include <math.h>
24#undef dfmal
25#undef f32xfmaf64
26#include <fenv.h>
27#include <ieee754.h>
28#include <math-barriers.h>
29#include <libm-alias-double.h>
30#include <math-narrow-alias.h>
31
32/* This implementation uses rounding to odd to avoid problems with
33 double rounding. See a paper by Boldo and Melquiond:
34 http://www.lri.fr/~melquion/doc/08-tc.pdf */
35
36double
37__fma (double x, double y, double z)
38{
39 if (__glibc_unlikely (!isfinite (x) || !isfinite (y)))
40 return x * y + z;
41 else if (__glibc_unlikely (!isfinite (z)))
42 /* If z is Inf, but x and y are finite, the result should be z
43 rather than NaN. */
44 return (z + x) + y;
45
46 /* Ensure correct sign of exact 0 + 0. */
47 if (__glibc_unlikely ((x == 0 || y == 0) && z == 0))
48 {
49 x = math_opt_barrier (x);
50 return x * y + z;
51 }
52
53 fenv_t env;
54 feholdexcept (envp: &env);
55 fesetround (FE_TONEAREST);
56
57 /* Multiplication m1 + m2 = x * y using Dekker's algorithm. */
58#define C ((1ULL << (LDBL_MANT_DIG + 1) / 2) + 1)
59 long double x1 = (long double) x * C;
60 long double y1 = (long double) y * C;
61 long double m1 = (long double) x * y;
62 x1 = (x - x1) + x1;
63 y1 = (y - y1) + y1;
64 long double x2 = x - x1;
65 long double y2 = y - y1;
66 long double m2 = (((x1 * y1 - m1) + x1 * y2) + x2 * y1) + x2 * y2;
67
68 /* Addition a1 + a2 = z + m1 using Knuth's algorithm. */
69 long double a1 = z + m1;
70 long double t1 = a1 - z;
71 long double t2 = a1 - t1;
72 t1 = m1 - t1;
73 t2 = z - t2;
74 long double a2 = t1 + t2;
75 /* Ensure the arithmetic is not scheduled after feclearexcept call. */
76 math_force_eval (m2);
77 math_force_eval (a2);
78 feclearexcept (FE_INEXACT);
79
80 /* If the result is an exact zero, ensure it has the correct sign. */
81 if (a1 == 0 && m2 == 0)
82 {
83 feupdateenv (envp: &env);
84 /* Ensure that round-to-nearest value of z + m1 is not reused. */
85 z = math_opt_barrier (z);
86 return z + m1;
87 }
88
89 fesetround (FE_TOWARDZERO);
90 /* Perform m2 + a2 addition with round to odd. */
91 a2 = a2 + m2;
92
93 /* Add that to a1 again using rounding to odd. */
94 union ieee854_long_double u;
95 u.d = a1 + a2;
96 if ((u.ieee.mantissa1 & 1) == 0 && u.ieee.exponent != 0x7fff)
97 u.ieee.mantissa1 |= fetestexcept (FE_INEXACT) != 0;
98 feupdateenv (envp: &env);
99
100 /* Add finally round to double precision. */
101 return u.d;
102}
103#ifndef __fma
104libm_alias_double (__fma, fma)
105libm_alias_double_narrow (__fma, fma)
106#endif
107

source code of glibc/sysdeps/ieee754/ldbl-96/s_fma.c