libquadmath/math/fmaq.c

   1 /* Compute x * y + z as ternary operation.
   2    Copyright (C) 2010-2017 Free Software Foundation, Inc.
   3    This file is part of the GNU C Library.
   4    Contributed by Jakub Jelinek <jakub@redhat.com>, 2010.
   5
   6    The GNU C Library is free software; you can redistribute it and/or
   7    modify it under the terms of the GNU Lesser General Public
   8    License as published by the Free Software Foundation; either
   9    version 2.1 of the License, or (at your option) any later version.
  10
  11    The GNU C Library is distributed in the hope that it will be useful,
  12    but WITHOUT ANY WARRANTY; without even the implied warranty of
  13    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  14    Lesser General Public License for more details.
  15
  16    You should have received a copy of the GNU Lesser General Public
  17    License along with the GNU C Library; if not, see
  18    <http://www.gnu.org/licenses/>.  */
  19
  20 #include "quadmath-imp.h"
  21 #include <math.h>
  22 #include <float.h>
  23 #ifdef HAVE_FENV_H
  24 # include <fenv.h>
  25 # if defined HAVE_FEHOLDEXCEPT && defined HAVE_FESETROUND \
  26      && defined HAVE_FEUPDATEENV && defined HAVE_FETESTEXCEPT \
  27      && defined FE_TOWARDZERO && defined FE_INEXACT
  28 #  define USE_FENV_H
  29 # endif
  30 #endif
  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 __float128
  37 fmaq (__float128 x, __float128 y, __float128 z)
  38 {
  39   ieee854_float128 u, v, w;
  40   int adjust = 0;
  41   u.value = x;
  42   v.value = y;
  43   w.value = z;
  44   if (__builtin_expect (u.ieee.exponent + v.ieee.exponent
  45                         >= 0x7fff + IEEE854_FLOAT128_BIAS
  46                            - FLT128_MANT_DIG, 0)
  47       || __builtin_expect (u.ieee.exponent >= 0x7fff - FLT128_MANT_DIG, 0)
  48       || __builtin_expect (v.ieee.exponent >= 0x7fff - FLT128_MANT_DIG, 0)
  49       || __builtin_expect (w.ieee.exponent >= 0x7fff - FLT128_MANT_DIG, 0)
  50       || __builtin_expect (u.ieee.exponent + v.ieee.exponent
  51                            <= IEEE854_FLOAT128_BIAS + FLT128_MANT_DIG, 0))
  52     {
  53       /* If z is Inf, but x and y are finite, the result should be
  54          z rather than NaN.  */
  55       if (w.ieee.exponent == 0x7fff
  56           && u.ieee.exponent != 0x7fff
  57           && v.ieee.exponent != 0x7fff)
  58         return (z + x) + y;
  59       /* If z is zero and x are y are nonzero, compute the result
  60          as x * y to avoid the wrong sign of a zero result if x * y
  61          underflows to 0.  */
  62       if (z == 0 && x != 0 && y != 0)
  63         return x * y;
  64       /* If x or y or z is Inf/NaN, or if x * y is zero, compute as
  65          x * y + z.  */
  66       if (u.ieee.exponent == 0x7fff
  67           || v.ieee.exponent == 0x7fff
  68           || w.ieee.exponent == 0x7fff
  69           || x == 0
  70           || y == 0)
  71         return x * y + z;
  72       /* If fma will certainly overflow, compute as x * y.  */
  73       if (u.ieee.exponent + v.ieee.exponent
  74           > 0x7fff + IEEE854_FLOAT128_BIAS)
  75         return x * y;
  76       /* If x * y is less than 1/4 of FLT128_DENORM_MIN, neither the
  77          result nor whether there is underflow depends on its exact
  78          value, only on its sign.  */
  79       if (u.ieee.exponent + v.ieee.exponent
  80           < IEEE854_FLOAT128_BIAS - FLT128_MANT_DIG - 2)
  81         {
  82           int neg = u.ieee.negative ^ v.ieee.negative;
  83           __float128 tiny = neg ? -0x1p-16494Q : 0x1p-16494Q;
  84           if (w.ieee.exponent >= 3)
  85             return tiny + z;
  86           /* Scaling up, adding TINY and scaling down produces the
  87              correct result, because in round-to-nearest mode adding
  88              TINY has no effect and in other modes double rounding is
  89              harmless.  But it may not produce required underflow
  90              exceptions.  */
  91           v.value = z * 0x1p114Q + tiny;
  92           if (TININESS_AFTER_ROUNDING
  93               ? v.ieee.exponent < 115
  94               : (w.ieee.exponent == 0
  95                  || (w.ieee.exponent == 1
  96                      && w.ieee.negative != neg
  97                      && w.ieee.mant_low == 0
  98                      && w.ieee.mant_high == 0)))
  99             {
 100               __float128 force_underflow = x * y;
 101               math_force_eval (force_underflow);
 102             }
 103           return v.value * 0x1p-114Q;
 104         }
 105       if (u.ieee.exponent + v.ieee.exponent
 106           >= 0x7fff + IEEE854_FLOAT128_BIAS - FLT128_MANT_DIG)
 107         {
 108           /* Compute 1p-113 times smaller result and multiply
 109              at the end.  */
 110           if (u.ieee.exponent > v.ieee.exponent)
 111             u.ieee.exponent -= FLT128_MANT_DIG;
 112           else
 113             v.ieee.exponent -= FLT128_MANT_DIG;
 114           /* If x + y exponent is very large and z exponent is very small,
 115              it doesn't matter if we don't adjust it.  */
 116           if (w.ieee.exponent > FLT128_MANT_DIG)
 117             w.ieee.exponent -= FLT128_MANT_DIG;
 118           adjust = 1;
 119         }
 120       else if (w.ieee.exponent >= 0x7fff - FLT128_MANT_DIG)
 121         {
 122           /* Similarly.
 123              If z exponent is very large and x and y exponents are
 124              very small, adjust them up to avoid spurious underflows,
 125              rather than down.  */
 126           if (u.ieee.exponent + v.ieee.exponent
 127               <= IEEE854_FLOAT128_BIAS + FLT128_MANT_DIG)
 128             {
 129               if (u.ieee.exponent > v.ieee.exponent)
 130                 u.ieee.exponent += 2 * FLT128_MANT_DIG + 2;
 131               else
 132                 v.ieee.exponent += 2 * FLT128_MANT_DIG + 2;
 133             }
 134           else if (u.ieee.exponent > v.ieee.exponent)
 135             {
 136               if (u.ieee.exponent > FLT128_MANT_DIG)
 137                 u.ieee.exponent -= FLT128_MANT_DIG;
 138             }
 139           else if (v.ieee.exponent > FLT128_MANT_DIG)
 140             v.ieee.exponent -= FLT128_MANT_DIG;
 141           w.ieee.exponent -= FLT128_MANT_DIG;
 142           adjust = 1;
 143         }
 144       else if (u.ieee.exponent >= 0x7fff - FLT128_MANT_DIG)
 145         {
 146           u.ieee.exponent -= FLT128_MANT_DIG;
 147           if (v.ieee.exponent)
 148             v.ieee.exponent += FLT128_MANT_DIG;
 149           else
 150             v.value *= 0x1p113Q;
 151         }
 152       else if (v.ieee.exponent >= 0x7fff - FLT128_MANT_DIG)
 153         {
 154           v.ieee.exponent -= FLT128_MANT_DIG;
 155           if (u.ieee.exponent)
 156             u.ieee.exponent += FLT128_MANT_DIG;
 157           else
 158             u.value *= 0x1p113Q;
 159         }
 160       else /* if (u.ieee.exponent + v.ieee.exponent
 161                   <= IEEE854_FLOAT128_BIAS + FLT128_MANT_DIG) */
 162         {
 163           if (u.ieee.exponent > v.ieee.exponent)
 164             u.ieee.exponent += 2 * FLT128_MANT_DIG + 2;
 165           else
 166             v.ieee.exponent += 2 * FLT128_MANT_DIG + 2;
 167           if (w.ieee.exponent <= 4 * FLT128_MANT_DIG + 6)
 168             {
 169               if (w.ieee.exponent)
 170                 w.ieee.exponent += 2 * FLT128_MANT_DIG + 2;
 171               else
 172                 w.value *= 0x1p228Q;
 173               adjust = -1;
 174             }
 175           /* Otherwise x * y should just affect inexact
 176              and nothing else.  */
 177         }
 178       x = u.value;
 179       y = v.value;
 180       z = w.value;
 181     }
 182
 183   /* Ensure correct sign of exact 0 + 0.  */
 184   if (__builtin_expect ((x == 0 || y == 0) && z == 0, 0))
 185     {
 186       x = math_opt_barrier (x);
 187       return x * y + z;
 188     }
 189
 190 #ifdef USE_FENV_H
 191   fenv_t env;
 192   feholdexcept (&env);
 193   fesetround (FE_TONEAREST);
 194 #endif
 195
 196   /* Multiplication m1 + m2 = x * y using Dekker's algorithm.  */
 197 #define C ((1LL << (FLT128_MANT_DIG + 1) / 2) + 1)
 198   __float128 x1 = x * C;
 199   __float128 y1 = y * C;
 200   __float128 m1 = x * y;
 201   x1 = (x - x1) + x1;
 202   y1 = (y - y1) + y1;
 203   __float128 x2 = x - x1;
 204   __float128 y2 = y - y1;
 205   __float128 m2 = (((x1 * y1 - m1) + x1 * y2) + x2 * y1) + x2 * y2;
 206
 207   /* Addition a1 + a2 = z + m1 using Knuth's algorithm.  */
 208   __float128 a1 = z + m1;
 209   __float128 t1 = a1 - z;
 210   __float128 t2 = a1 - t1;
 211   t1 = m1 - t1;
 212   t2 = z - t2;
 213   __float128 a2 = t1 + t2;
 214   /* Ensure the arithmetic is not scheduled after feclearexcept call.  */
 215   math_force_eval (m2);
 216   math_force_eval (a2);
 217 #ifdef USE_FENV_H
 218   feclearexcept (FE_INEXACT);
 219 #endif
 220
 221   /* If the result is an exact zero, ensure it has the correct sign.  */
 222   if (a1 == 0 && m2 == 0)
 223     {
 224 #ifdef USE_FENV_H
 225       feupdateenv (&env);
 226 #endif
 227       /* Ensure that round-to-nearest value of z + m1 is not reused.  */
 228       z = math_opt_barrier (z);
 229       return z + m1;
 230     }
 231
 232 #ifdef USE_FENV_H
 233   fesetround (FE_TOWARDZERO);
 234 #endif
 235   /* Perform m2 + a2 addition with round to odd.  */
 236   u.value = a2 + m2;
 237
 238   if (__builtin_expect (adjust == 0, 1))
 239     {
 240 #ifdef USE_FENV_H
 241       if ((u.ieee.mant_low & 1) == 0 && u.ieee.exponent != 0x7fff)
 242         u.ieee.mant_low |= fetestexcept (FE_INEXACT) != 0;
 243       feupdateenv (&env);
 244 #endif
 245       /* Result is a1 + u.value.  */
 246       return a1 + u.value;
 247     }
 248   else if (__builtin_expect (adjust > 0, 1))
 249     {
 250 #ifdef USE_FENV_H
 251       if ((u.ieee.mant_low & 1) == 0 && u.ieee.exponent != 0x7fff)
 252         u.ieee.mant_low |= fetestexcept (FE_INEXACT) != 0;
 253       feupdateenv (&env);
 254 #endif
 255       /* Result is a1 + u.value, scaled up.  */
 256       return (a1 + u.value) * 0x1p113Q;
 257     }
 258   else
 259     {
 260 #ifdef USE_FENV_H
 261       if ((u.ieee.mant_low & 1) == 0)
 262         u.ieee.mant_low |= fetestexcept (FE_INEXACT) != 0;
 263 #endif
 264       v.value = a1 + u.value;
 265       /* Ensure the addition is not scheduled after fetestexcept call.  */
 266       asm volatile ("" : : "m" (v.value));
 267 #ifdef USE_FENV_H
 268       int j = fetestexcept (FE_INEXACT) != 0;
 269       feupdateenv (&env);
 270 #else
 271       int j = 0;
 272 #endif
 273       /* Ensure the following computations are performed in default rounding
 274          mode instead of just reusing the round to zero computation.  */
 275       asm volatile ("" : "=m" (u) : "m" (u));
 276       /* If a1 + u.value is exact, the only rounding happens during
 277          scaling down.  */
 278       if (j == 0)
 279         return v.value * 0x1p-228Q;
 280       /* If result rounded to zero is not subnormal, no double
 281          rounding will occur.  */
 282       if (v.ieee.exponent > 228)
 283         return (a1 + u.value) * 0x1p-228Q;
 284       /* If v.value * 0x1p-228Q with round to zero is a subnormal above
 285          or equal to FLT128_MIN / 2, then v.value * 0x1p-228Q shifts mantissa
 286          down just by 1 bit, which means v.ieee.mant_low |= j would
 287          change the round bit, not sticky or guard bit.
 288          v.value * 0x1p-228Q never normalizes by shifting up,
 289          so round bit plus sticky bit should be already enough
 290          for proper rounding.  */
 291       if (v.ieee.exponent == 228)
 292         {
 293           /* If the exponent would be in the normal range when
 294              rounding to normal precision with unbounded exponent
 295              range, the exact result is known and spurious underflows
 296              must be avoided on systems detecting tininess after
 297              rounding.  */
 298           if (TININESS_AFTER_ROUNDING)
 299             {
 300               w.value = a1 + u.value;
 301               if (w.ieee.exponent == 229)
 302                 return w.value * 0x1p-228Q;
 303             }
 304           /* v.ieee.mant_low & 2 is LSB bit of the result before rounding,
 305              v.ieee.mant_low & 1 is the round bit and j is our sticky
 306              bit. */
 307           w.value = 0.0Q;
 308           w.ieee.mant_low = ((v.ieee.mant_low & 3) << 1) | j;
 309           w.ieee.negative = v.ieee.negative;
 310           v.ieee.mant_low &= ~3U;
 311           v.value *= 0x1p-228Q;
 312           w.value *= 0x1p-2Q;
 313           return v.value + w.value;
 314         }
 315       v.ieee.mant_low |= j;
 316       return v.value * 0x1p-228Q;
 317     }
 318 }