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  #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  
      36  double
      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 (&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 (&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 (&env);
      99  
     100    /* Add finally round to double precision.  */
     101    return u.d;
     102  }
     103  #ifndef __fma
     104  libm_alias_double (__fma, fma)
     105  libm_alias_double_narrow (__fma, fma)
     106  #endif