lib/msun/src/s_fma.c

   1 /*-
   2  * Copyright (c) 2005 David Schultz <das@FreeBSD.ORG>
   3  * All rights reserved.
   4  *
   5  * Redistribution and use in source and binary forms, with or without
   6  * modification, are permitted provided that the following conditions
   7  * are met:
   8  * 1. Redistributions of source code must retain the above copyright
   9  *    notice, this list of conditions and the following disclaimer.
  10  * 2. Redistributions in binary form must reproduce the above copyright
  11  *    notice, this list of conditions and the following disclaimer in the
  12  *    documentation and/or other materials provided with the distribution.
  13  *
  14  * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
  15  * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
  16  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
  17  * ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
  18  * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
  19  * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
  20  * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
  21  * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
  22  * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
  23  * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
  24  * SUCH DAMAGE.
  25  */
  26
  27 #include <sys/cdefs.h>
  28 __FBSDID("$FreeBSD$");
  29
  30 #include <fenv.h>
  31 #include <float.h>
  32 #include <math.h>
  33
  34 /*
  35  * Fused multiply-add: Compute x * y + z with a single rounding error.
  36  *
  37  * We use scaling to avoid overflow/underflow, along with the
  38  * canonical precision-doubling technique adapted from:
  39  *
  40  *      Dekker, T.  A Floating-Point Technique for Extending the
  41  *      Available Precision.  Numer. Math. 18, 224-242 (1971).
  42  *
  43  * This algorithm is sensitive to the rounding precision.  FPUs such
  44  * as the i387 must be set in double-precision mode if variables are
  45  * to be stored in FP registers in order to avoid incorrect results.
  46  * This is the default on FreeBSD, but not on many other systems.
  47  *
  48  * Hardware instructions should be used on architectures that support it,
  49  * since this implementation will likely be several times slower.
  50  */
  51 #if LDBL_MANT_DIG != 113
  52 double
  53 fma(double x, double y, double z)
  54 {
  55         static const double split = 0x1p27 + 1.0;
  56         double xs, ys, zs;
  57         double c, cc, hx, hy, p, q, tx, ty;
  58         double r, rr, s;
  59         int oround;
  60         int ex, ey, ez;
  61         int spread;
  62
  63         /*
  64          * Handle special cases. The order of operations and the particular
  65          * return values here are crucial in handling special cases involving
  66          * infinities, NaNs, overflows, and signed zeroes correctly.
  67          */
  68         if (x == 0.0 || y == 0.0)
  69                 return (x * y + z);
  70         if (z == 0.0)
  71                 return (x * y);
  72         if (!isfinite(x) || !isfinite(y))
  73                 return (x * y + z);
  74         if (!isfinite(z))
  75                 return (z);
  76
  77         xs = frexp(x, &ex);
  78         ys = frexp(y, &ey);
  79         zs = frexp(z, &ez);
  80         oround = fegetround();
  81         spread = ex + ey - ez;
  82
  83         /*
  84          * If x * y and z are many orders of magnitude apart, the scaling
  85          * will overflow, so we handle these cases specially.  Rounding
  86          * modes other than FE_TONEAREST are painful.
  87          */
  88         if (spread > DBL_MANT_DIG * 2) {
  89                 fenv_t env;
  90                 feraiseexcept(FE_INEXACT);
  91                 switch(oround) {
  92                 case FE_TONEAREST:
  93                         return (x * y);
  94                 case FE_TOWARDZERO:
  95                         if (x > 0.0 ^ y < 0.0 ^ z < 0.0)
  96                                 return (x * y);
  97                         feholdexcept(&env);
  98                         r = x * y;
  99                         if (!fetestexcept(FE_INEXACT))
 100                                 r = nextafter(r, 0);
 101                         feupdateenv(&env);
 102                         return (r);
 103                 case FE_DOWNWARD:
 104                         if (z > 0.0)
 105                                 return (x * y);
 106                         feholdexcept(&env);
 107                         r = x * y;
 108                         if (!fetestexcept(FE_INEXACT))
 109                                 r = nextafter(r, -INFINITY);
 110                         feupdateenv(&env);
 111                         return (r);
 112                 default:        /* FE_UPWARD */
 113                         if (z < 0.0)
 114                                 return (x * y);
 115                         feholdexcept(&env);
 116                         r = x * y;
 117                         if (!fetestexcept(FE_INEXACT))
 118                                 r = nextafter(r, INFINITY);
 119                         feupdateenv(&env);
 120                         return (r);
 121                 }
 122         }
 123         if (spread < -DBL_MANT_DIG) {
 124                 feraiseexcept(FE_INEXACT);
 125                 if (!isnormal(z))
 126                         feraiseexcept(FE_UNDERFLOW);
 127                 switch (oround) {
 128                 case FE_TONEAREST:
 129                         return (z);
 130                 case FE_TOWARDZERO:
 131                         if (x > 0.0 ^ y < 0.0 ^ z < 0.0)
 132                                 return (z);
 133                         else
 134                                 return (nextafter(z, 0));
 135                 case FE_DOWNWARD:
 136                         if (x > 0.0 ^ y < 0.0)
 137                                 return (z);
 138                         else
 139                                 return (nextafter(z, -INFINITY));
 140                 default:        /* FE_UPWARD */
 141                         if (x > 0.0 ^ y < 0.0)
 142                                 return (nextafter(z, INFINITY));
 143                         else
 144                                 return (z);
 145                 }
 146         }
 147
 148         /*
 149          * Use Dekker's algorithm to perform the multiplication and
 150          * subsequent addition in twice the machine precision.
 151          * Arrange so that x * y = c + cc, and x * y + z = r + rr.
 152          */
 153         fesetround(FE_TONEAREST);
 154
 155         p = xs * split;
 156         hx = xs - p;
 157         hx += p;
 158         tx = xs - hx;
 159
 160         p = ys * split;
 161         hy = ys - p;
 162         hy += p;
 163         ty = ys - hy;
 164
 165         p = hx * hy;
 166         q = hx * ty + tx * hy;
 167         c = p + q;
 168         cc = p - c + q + tx * ty;
 169
 170         zs = ldexp(zs, -spread);
 171         r = c + zs;
 172         s = r - c;
 173         rr = (c - (r - s)) + (zs - s) + cc;
 174
 175         spread = ex + ey;
 176         if (spread + ilogb(r) > -1023) {
 177                 fesetround(oround);
 178                 r = r + rr;
 179         } else {
 180                 /*
 181                  * The result is subnormal, so we round before scaling to
 182                  * avoid double rounding.
 183                  */
 184                 p = ldexp(copysign(0x1p-1022, r), -spread);
 185                 c = r + p;
 186                 s = c - r;
 187                 cc = (r - (c - s)) + (p - s) + rr;
 188                 fesetround(oround);
 189                 r = (c + cc) - p;
 190         }
 191         return (ldexp(r, spread));
 192 }
 193 #else   /* LDBL_MANT_DIG == 113 */
 194 /*
 195  * 113 bits of precision is more than twice the precision of a double,
 196  * so it is enough to represent the intermediate product exactly.
 197  */
 198 double
 199 fma(double x, double y, double z)
 200 {
 201         return ((long double)x * y + z);
 202 }
 203 #endif  /* LDBL_MANT_DIG != 113 */
 204
 205 #if (LDBL_MANT_DIG == 53)
 206 __weak_reference(fma, fmal);
 207 #endif