math/mul_split.h

   1 /* Compute full X * Y for double type.
   2    Copyright (C) 2013-2017 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    <http://www.gnu.org/licenses/>.  */
  18
  19 #ifndef _MUL_SPLIT_H
  20 #define _MUL_SPLIT_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 void
  29 mul_split (double *hi, double *lo, double x, double y)
  30 {
  31 #ifdef __FP_FAST_FMA
  32   /* Fast built-in fused multiply-add.  */
  33   *hi = x * y;
  34   *lo = __builtin_fma (x, y, -*hi);
  35 #else
  36   /* Apply Dekker's algorithm.  */
  37   *hi = x * y;
  38 # define C ((1 << (DBL_MANT_DIG + 1) / 2) + 1)
  39   double x1 = x * C;
  40   double y1 = y * C;
  41 # undef C
  42   x1 = (x - x1) + x1;
  43   y1 = (y - y1) + y1;
  44   double x2 = x - x1;
  45   double y2 = y - y1;
  46   *lo = (((x1 * y1 - *hi) + x1 * y2) + x2 * y1) + x2 * y2;
  47 #endif
  48 }
  49
  50 #endif /* _MUL_SPLIT_H */