1  /* Compute full X * Y for long double type.
       2     Copyright (C) 2013-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  #ifndef _MUL_SPLITL_H
      20  #define _MUL_SPLITL_H
      21  
      22  #include <float.h>
      23  
      24  /* Calculate X * Y exactly and store the result in *HI + *LO.  It is
      25     given that the values are small enough that no overflow occurs and
      26     large enough (or zero) that no underflow occurs.  */
      27  
      28  static inline void
      29  mul_splitl (long double *hi, long double *lo, long double x, long double y)
      30  {
      31  #ifdef __FP_FAST_FMAL
      32    /* Fast built-in fused multiply-add.  */
      33    *hi = x * y;
      34    *lo = __builtin_fmal (x, y, -*hi);
      35  #else
      36    /* Apply Dekker's algorithm.  */
      37    *hi = x * y;
      38  # define C ((1LL << (LDBL_MANT_DIG + 1) / 2) + 1)
      39    long double x1 = x * C;
      40    long double y1 = y * C;
      41  # undef C
      42    x1 = (x - x1) + x1;
      43    y1 = (y - y1) + y1;
      44    long double x2 = x - x1;
      45    long double y2 = y - y1;
      46    *lo = (((x1 * y1 - *hi) + x1 * y2) + x2 * y1) + x2 * y2;
      47  #endif
      48  }
      49  
      50  #endif /* _MUL_SPLITL_H */