1  /* Compute x * y + z as ternary operation.
       2     Copyright (C) 2010-2023 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
      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