sysdeps/ieee754/ldbl-128/s_fmal.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 <float.h>
  21 #include <math.h>
  22 #include <fenv.h>
  23 #include <ieee754.h>
  24 #include <math_private.h>
  25 #include <libm-alias-ldouble.h>
  26 #include <tininess.h>
  27
  28 /* This implementation uses rounding to odd to avoid problems with
  29    double rounding.  See a paper by Boldo and Melquiond:
  30    http://www.lri.fr/~melquion/doc/08-tc.pdf  */
  31
  32 _Float128
  33 __fmal (_Float128 x, _Float128 y, _Float128 z)
  34 {
  35   union ieee854_long_double u, v, w;
  36   int adjust = 0;
  37   u.d = x;
  38   v.d = y;
  39   w.d = z;
  40   if (__builtin_expect (u.ieee.exponent + v.ieee.exponent
  41                         >= 0x7fff + IEEE854_LONG_DOUBLE_BIAS
  42                            - LDBL_MANT_DIG, 0)
  43       || __builtin_expect (u.ieee.exponent >= 0x7fff - LDBL_MANT_DIG, 0)
  44       || __builtin_expect (v.ieee.exponent >= 0x7fff - LDBL_MANT_DIG, 0)
  45       || __builtin_expect (w.ieee.exponent >= 0x7fff - LDBL_MANT_DIG, 0)
  46       || __builtin_expect (u.ieee.exponent + v.ieee.exponent
  47                            <= IEEE854_LONG_DOUBLE_BIAS + LDBL_MANT_DIG, 0))
  48     {
  49       /* If z is Inf, but x and y are finite, the result should be
  50          z rather than NaN.  */
  51       if (w.ieee.exponent == 0x7fff
  52           && u.ieee.exponent != 0x7fff
  53           && v.ieee.exponent != 0x7fff)
  54         return (z + x) + y;
  55       /* If z is zero and x are y are nonzero, compute the result
  56          as x * y to avoid the wrong sign of a zero result if x * y
  57          underflows to 0.  */
  58       if (z == 0 && x != 0 && y != 0)
  59         return x * y;
  60       /* If x or y or z is Inf/NaN, or if x * y is zero, compute as
  61          x * y + z.  */
  62       if (u.ieee.exponent == 0x7fff
  63           || v.ieee.exponent == 0x7fff
  64           || w.ieee.exponent == 0x7fff
  65           || x == 0
  66           || y == 0)
  67         return x * y + z;
  68       /* If fma will certainly overflow, compute as x * y.  */
  69       if (u.ieee.exponent + v.ieee.exponent
  70           > 0x7fff + IEEE854_LONG_DOUBLE_BIAS)
  71         return x * y;
  72       /* If x * y is less than 1/4 of LDBL_TRUE_MIN, neither the
  73          result nor whether there is underflow depends on its exact
  74          value, only on its sign.  */
  75       if (u.ieee.exponent + v.ieee.exponent
  76           < IEEE854_LONG_DOUBLE_BIAS - LDBL_MANT_DIG - 2)
  77         {
  78           int neg = u.ieee.negative ^ v.ieee.negative;
  79           _Float128 tiny = neg ? L(-0x1p-16494) : L(0x1p-16494);
  80           if (w.ieee.exponent >= 3)
  81             return tiny + z;
  82           /* Scaling up, adding TINY and scaling down produces the
  83              correct result, because in round-to-nearest mode adding
  84              TINY has no effect and in other modes double rounding is
  85              harmless.  But it may not produce required underflow
  86              exceptions.  */
  87           v.d = z * L(0x1p114) + tiny;
  88           if (TININESS_AFTER_ROUNDING
  89               ? v.ieee.exponent < 115
  90               : (w.ieee.exponent == 0
  91                  || (w.ieee.exponent == 1
  92                      && w.ieee.negative != neg
  93                      && w.ieee.mantissa3 == 0
  94                      && w.ieee.mantissa2 == 0
  95                      && w.ieee.mantissa1 == 0
  96                      && w.ieee.mantissa0 == 0)))
  97             {
  98               _Float128 force_underflow = x * y;
  99               math_force_eval (force_underflow);
 100             }
 101           return v.d * L(0x1p-114);
 102         }
 103       if (u.ieee.exponent + v.ieee.exponent
 104           >= 0x7fff + IEEE854_LONG_DOUBLE_BIAS - LDBL_MANT_DIG)
 105         {
 106           /* Compute 1p-113 times smaller result and multiply
 107              at the end.  */
 108           if (u.ieee.exponent > v.ieee.exponent)
 109             u.ieee.exponent -= LDBL_MANT_DIG;
 110           else
 111             v.ieee.exponent -= LDBL_MANT_DIG;
 112           /* If x + y exponent is very large and z exponent is very small,
 113              it doesn't matter if we don't adjust it.  */
 114           if (w.ieee.exponent > LDBL_MANT_DIG)
 115             w.ieee.exponent -= LDBL_MANT_DIG;
 116           adjust = 1;
 117         }
 118       else if (w.ieee.exponent >= 0x7fff - LDBL_MANT_DIG)
 119         {
 120           /* Similarly.
 121              If z exponent is very large and x and y exponents are
 122              very small, adjust them up to avoid spurious underflows,
 123              rather than down.  */
 124           if (u.ieee.exponent + v.ieee.exponent
 125               <= IEEE854_LONG_DOUBLE_BIAS + 2 * LDBL_MANT_DIG)
 126             {
 127               if (u.ieee.exponent > v.ieee.exponent)
 128                 u.ieee.exponent += 2 * LDBL_MANT_DIG + 2;
 129               else
 130                 v.ieee.exponent += 2 * LDBL_MANT_DIG + 2;
 131             }
 132           else if (u.ieee.exponent > v.ieee.exponent)
 133             {
 134               if (u.ieee.exponent > LDBL_MANT_DIG)
 135                 u.ieee.exponent -= LDBL_MANT_DIG;
 136             }
 137           else if (v.ieee.exponent > LDBL_MANT_DIG)
 138             v.ieee.exponent -= LDBL_MANT_DIG;
 139           w.ieee.exponent -= LDBL_MANT_DIG;
 140           adjust = 1;
 141         }
 142       else if (u.ieee.exponent >= 0x7fff - LDBL_MANT_DIG)
 143         {
 144           u.ieee.exponent -= LDBL_MANT_DIG;
 145           if (v.ieee.exponent)
 146             v.ieee.exponent += LDBL_MANT_DIG;
 147           else
 148             v.d *= L(0x1p113);
 149         }
 150       else if (v.ieee.exponent >= 0x7fff - LDBL_MANT_DIG)
 151         {
 152           v.ieee.exponent -= LDBL_MANT_DIG;
 153           if (u.ieee.exponent)
 154             u.ieee.exponent += LDBL_MANT_DIG;
 155           else
 156             u.d *= L(0x1p113);
 157         }
 158       else /* if (u.ieee.exponent + v.ieee.exponent
 159                   <= IEEE854_LONG_DOUBLE_BIAS + LDBL_MANT_DIG) */
 160         {
 161           if (u.ieee.exponent > v.ieee.exponent)
 162             u.ieee.exponent += 2 * LDBL_MANT_DIG + 2;
 163           else
 164             v.ieee.exponent += 2 * LDBL_MANT_DIG + 2;
 165           if (w.ieee.exponent <= 4 * LDBL_MANT_DIG + 6)
 166             {
 167               if (w.ieee.exponent)
 168                 w.ieee.exponent += 2 * LDBL_MANT_DIG + 2;
 169               else
 170                 w.d *= L(0x1p228);
 171               adjust = -1;
 172             }
 173           /* Otherwise x * y should just affect inexact
 174              and nothing else.  */
 175         }
 176       x = u.d;
 177       y = v.d;
 178       z = w.d;
 179     }
 180
 181   /* Ensure correct sign of exact 0 + 0.  */
 182   if (__glibc_unlikely ((x == 0 || y == 0) && z == 0))
 183     {
 184       x = math_opt_barrier (x);
 185       return x * y + z;
 186     }
 187
 188   fenv_t env;
 189   feholdexcept (&env);
 190   fesetround (FE_TONEAREST);
 191
 192   /* Multiplication m1 + m2 = x * y using Dekker's algorithm.  */
 193 #define C ((1LL << (LDBL_MANT_DIG + 1) / 2) + 1)
 194   _Float128 x1 = x * C;
 195   _Float128 y1 = y * C;
 196   _Float128 m1 = x * y;
 197   x1 = (x - x1) + x1;
 198   y1 = (y - y1) + y1;
 199   _Float128 x2 = x - x1;
 200   _Float128 y2 = y - y1;
 201   _Float128 m2 = (((x1 * y1 - m1) + x1 * y2) + x2 * y1) + x2 * y2;
 202
 203   /* Addition a1 + a2 = z + m1 using Knuth's algorithm.  */
 204   _Float128 a1 = z + m1;
 205   _Float128 t1 = a1 - z;
 206   _Float128 t2 = a1 - t1;
 207   t1 = m1 - t1;
 208   t2 = z - t2;
 209   _Float128 a2 = t1 + t2;
 210   /* Ensure the arithmetic is not scheduled after feclearexcept call.  */
 211   math_force_eval (m2);
 212   math_force_eval (a2);
 213   feclearexcept (FE_INEXACT);
 214
 215   /* If the result is an exact zero, ensure it has the correct sign.  */
 216   if (a1 == 0 && m2 == 0)
 217     {
 218       feupdateenv (&env);
 219       /* Ensure that round-to-nearest value of z + m1 is not reused.  */
 220       z = math_opt_barrier (z);
 221       return z + m1;
 222     }
 223
 224   fesetround (FE_TOWARDZERO);
 225   /* Perform m2 + a2 addition with round to odd.  */
 226   u.d = a2 + m2;
 227
 228   if (__glibc_likely (adjust == 0))
 229     {
 230       if ((u.ieee.mantissa3 & 1) == 0 && u.ieee.exponent != 0x7fff)
 231         u.ieee.mantissa3 |= fetestexcept (FE_INEXACT) != 0;
 232       feupdateenv (&env);
 233       /* Result is a1 + u.d.  */
 234       return a1 + u.d;
 235     }
 236   else if (__glibc_likely (adjust > 0))
 237     {
 238       if ((u.ieee.mantissa3 & 1) == 0 && u.ieee.exponent != 0x7fff)
 239         u.ieee.mantissa3 |= fetestexcept (FE_INEXACT) != 0;
 240       feupdateenv (&env);
 241       /* Result is a1 + u.d, scaled up.  */
 242       return (a1 + u.d) * L(0x1p113);
 243     }
 244   else
 245     {
 246       if ((u.ieee.mantissa3 & 1) == 0)
 247         u.ieee.mantissa3 |= fetestexcept (FE_INEXACT) != 0;
 248       v.d = a1 + u.d;
 249       /* Ensure the addition is not scheduled after fetestexcept call.  */
 250       math_force_eval (v.d);
 251       int j = fetestexcept (FE_INEXACT) != 0;
 252       feupdateenv (&env);
 253       /* Ensure the following computations are performed in default rounding
 254          mode instead of just reusing the round to zero computation.  */
 255       asm volatile ("" : "=m" (u) : "m" (u));
 256       /* If a1 + u.d is exact, the only rounding happens during
 257          scaling down.  */
 258       if (j == 0)
 259         return v.d * L(0x1p-228);
 260       /* If result rounded to zero is not subnormal, no double
 261          rounding will occur.  */
 262       if (v.ieee.exponent > 228)
 263         return (a1 + u.d) * L(0x1p-228);
 264       /* If v.d * 0x1p-228L with round to zero is a subnormal above
 265          or equal to LDBL_MIN / 2, then v.d * 0x1p-228L shifts mantissa
 266          down just by 1 bit, which means v.ieee.mantissa3 |= j would
 267          change the round bit, not sticky or guard bit.
 268          v.d * 0x1p-228L never normalizes by shifting up,
 269          so round bit plus sticky bit should be already enough
 270          for proper rounding.  */
 271       if (v.ieee.exponent == 228)
 272         {
 273           /* If the exponent would be in the normal range when
 274              rounding to normal precision with unbounded exponent
 275              range, the exact result is known and spurious underflows
 276              must be avoided on systems detecting tininess after
 277              rounding.  */
 278           if (TININESS_AFTER_ROUNDING)
 279             {
 280               w.d = a1 + u.d;
 281               if (w.ieee.exponent == 229)
 282                 return w.d * L(0x1p-228);
 283             }
 284           /* v.ieee.mantissa3 & 2 is LSB bit of the result before rounding,
 285              v.ieee.mantissa3 & 1 is the round bit and j is our sticky
 286              bit.  */
 287           w.d = 0;
 288           w.ieee.mantissa3 = ((v.ieee.mantissa3 & 3) << 1) | j;
 289           w.ieee.negative = v.ieee.negative;
 290           v.ieee.mantissa3 &= ~3U;
 291           v.d *= L(0x1p-228);
 292           w.d *= L(0x1p-2);
 293           return v.d + w.d;
 294         }
 295       v.ieee.mantissa3 |= j;
 296       return v.d * L(0x1p-228);
 297     }
 298 }
 299 libm_alias_ldouble (__fma, fma)