sysdeps/ieee754/ldbl-96/s_fmal.c

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