libquadmath/math/fmaq.c

   1 /* Compute x * y + z as ternary operation.
   2    Copyright (C) 2010-2012 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               volatile __float128 force_underflow = x * y;
 101               (void) 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;
 165           else
 166             v.ieee.exponent += 2 * FLT128_MANT_DIG;
 167           if (w.ieee.exponent <= 4 * FLT128_MANT_DIG + 4)
 168             {
 169               if (w.ieee.exponent)
 170                 w.ieee.exponent += 2 * FLT128_MANT_DIG;
 171               else
 172                 w.value *= 0x1p226Q;
 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     return x * y + z;
 186
 187 #ifdef USE_FENV_H
 188   fenv_t env;
 189   feholdexcept (&env);
 190   fesetround (FE_TONEAREST);
 191 #endif
 192
 193   /* Multiplication m1 + m2 = x * y using Dekker's algorithm.  */
 194 #define C ((1LL << (FLT128_MANT_DIG + 1) / 2) + 1)
 195   __float128 x1 = x * C;
 196   __float128 y1 = y * C;
 197   __float128 m1 = x * y;
 198   x1 = (x - x1) + x1;
 199   y1 = (y - y1) + y1;
 200   __float128 x2 = x - x1;
 201   __float128 y2 = y - y1;
 202   __float128 m2 = (((x1 * y1 - m1) + x1 * y2) + x2 * y1) + x2 * y2;
 203
 204   /* Addition a1 + a2 = z + m1 using Knuth's algorithm.  */
 205   __float128 a1 = z + m1;
 206   __float128 t1 = a1 - z;
 207   __float128 t2 = a1 - t1;
 208   t1 = m1 - t1;
 209   t2 = z - t2;
 210   __float128 a2 = t1 + t2;
 211 #ifdef USE_FENV_H
 212   feclearexcept (FE_INEXACT);
 213 #endif
 214
 215   /* If the result is an exact zero, ensure it has the correct
 216      sign.  */
 217   if (a1 == 0 && m2 == 0)
 218     {
 219 #ifdef USE_FENV_H
 220       feupdateenv (&env);
 221 #endif
 222       /* Ensure that round-to-nearest value of z + m1 is not
 223          reused.  */
 224       asm volatile ("" : "=m" (z) : "m" (z));
 225       return z + m1;
 226     }
 227
 228
 229 #ifdef USE_FENV_H
 230   fesetround (FE_TOWARDZERO);
 231 #endif
 232   /* Perform m2 + a2 addition with round to odd.  */
 233   u.value = a2 + m2;
 234
 235   if (__builtin_expect (adjust == 0, 1))
 236     {
 237 #ifdef USE_FENV_H
 238       if ((u.ieee.mant_low & 1) == 0 && u.ieee.exponent != 0x7fff)
 239         u.ieee.mant_low |= fetestexcept (FE_INEXACT) != 0;
 240       feupdateenv (&env);
 241 #endif
 242       /* Result is a1 + u.value.  */
 243       return a1 + u.value;
 244     }
 245   else if (__builtin_expect (adjust > 0, 1))
 246     {
 247 #ifdef USE_FENV_H
 248       if ((u.ieee.mant_low & 1) == 0 && u.ieee.exponent != 0x7fff)
 249         u.ieee.mant_low |= fetestexcept (FE_INEXACT) != 0;
 250       feupdateenv (&env);
 251 #endif
 252       /* Result is a1 + u.value, scaled up.  */
 253       return (a1 + u.value) * 0x1p113Q;
 254     }
 255   else
 256     {
 257 #ifdef USE_FENV_H
 258       if ((u.ieee.mant_low & 1) == 0)
 259         u.ieee.mant_low |= fetestexcept (FE_INEXACT) != 0;
 260 #endif
 261       v.value = a1 + u.value;
 262       /* Ensure the addition is not scheduled after fetestexcept call.  */
 263       asm volatile ("" : : "m" (v.value));
 264 #ifdef USE_FENV_H
 265       int j = fetestexcept (FE_INEXACT) != 0;
 266       feupdateenv (&env);
 267 #else
 268       int j = 0;
 269 #endif
 270       /* Ensure the following computations are performed in default rounding
 271          mode instead of just reusing the round to zero computation.  */
 272       asm volatile ("" : "=m" (u) : "m" (u));
 273       /* If a1 + u.value is exact, the only rounding happens during
 274          scaling down.  */
 275       if (j == 0)
 276         return v.value * 0x1p-226Q;
 277       /* If result rounded to zero is not subnormal, no double
 278          rounding will occur.  */
 279       if (v.ieee.exponent > 226)
 280         return (a1 + u.value) * 0x1p-226Q;
 281       /* If v.value * 0x1p-226Q with round to zero is a subnormal above
 282          or equal to FLT128_MIN / 2, then v.value * 0x1p-226Q shifts mantissa
 283          down just by 1 bit, which means v.ieee.mant_low |= j would
 284          change the round bit, not sticky or guard bit.
 285          v.value * 0x1p-226Q never normalizes by shifting up,
 286          so round bit plus sticky bit should be already enough
 287          for proper rounding.  */
 288       if (v.ieee.exponent == 226)
 289         {
 290           /* If the exponent would be in the normal range when
 291              rounding to normal precision with unbounded exponent
 292              range, the exact result is known and spurious underflows
 293              must be avoided on systems detecting tininess after
 294              rounding.  */
 295           if (TININESS_AFTER_ROUNDING)
 296             {
 297               w.value = a1 + u.value;
 298               if (w.ieee.exponent == 227)
 299                 return w.value * 0x1p-226Q;
 300             }
 301           /* v.ieee.mant_low & 2 is LSB bit of the result before rounding,
 302              v.ieee.mant_low & 1 is the round bit and j is our sticky
 303              bit. */
 304           w.value = 0.0Q;
 305           w.ieee.mant_low = ((v.ieee.mant_low & 3) << 1) | j;
 306           w.ieee.negative = v.ieee.negative;
 307           v.ieee.mant_low &= ~3U;
 308           v.value *= 0x1p-226Q;
 309           w.value *= 0x1p-2Q;
 310           return v.value + w.value;
 311         }
 312       v.ieee.mant_low |= j;
 313       return v.value * 0x1p-226Q;
 314     }
 315 }