1 /*- 2 * Copyright (c) 2005-2011 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 "cdefs-compat.h" 28 //__FBSDID("$FreeBSD: src/lib/msun/src/s_fma.c,v 1.8 2011/10/21 06:30:43 das Exp $"); 29 30 #include <float.h> 31 #include <openlibm_fenv.h> 32 #include <openlibm_math.h> 33 34 #include "math_private.h" 35 36 /* 37 * A struct dd represents a floating-point number with twice the precision 38 * of a double. We maintain the invariant that "hi" stores the 53 high-order 39 * bits of the result. 40 */ 41 struct dd { 42 double hi; 43 double lo; 44 }; 45 46 /* 47 * Compute a+b exactly, returning the exact result in a struct dd. We assume 48 * that both a and b are finite, but make no assumptions about their relative 49 * magnitudes. 50 */ 51 static inline struct dd 52 dd_add(double a, double b) 53 { 54 struct dd ret; 55 double s; 56 57 ret.hi = a + b; 58 s = ret.hi - a; 59 ret.lo = (a - (ret.hi - s)) + (b - s); 60 return (ret); 61 } 62 63 /* 64 * Compute a+b, with a small tweak: The least significant bit of the 65 * result is adjusted into a sticky bit summarizing all the bits that 66 * were lost to rounding. This adjustment negates the effects of double 67 * rounding when the result is added to another number with a higher 68 * exponent. For an explanation of round and sticky bits, see any reference 69 * on FPU design, e.g., 70 * 71 * J. Coonen. An Implementation Guide to a Proposed Standard for 72 * Floating-Point Arithmetic. Computer, vol. 13, no. 1, Jan 1980. 73 */ 74 static inline double 75 add_adjusted(double a, double b) 76 { 77 struct dd sum; 78 u_int64_t hibits, lobits; 79 80 sum = dd_add(a, b); 81 if (sum.lo != 0) { 82 EXTRACT_WORD64(hibits, sum.hi); 83 if ((hibits & 1) == 0) { 84 /* hibits += (int)copysign(1.0, sum.hi * sum.lo) */ 85 EXTRACT_WORD64(lobits, sum.lo); 86 hibits += 1 - ((hibits ^ lobits) >> 62); 87 INSERT_WORD64(sum.hi, hibits); 88 } 89 } 90 return (sum.hi); 91 } 92 93 /* 94 * Compute ldexp(a+b, scale) with a single rounding error. It is assumed 95 * that the result will be subnormal, and care is taken to ensure that 96 * double rounding does not occur. 97 */ 98 static inline double 99 add_and_denormalize(double a, double b, int scale) 100 { 101 struct dd sum; 102 u_int64_t hibits, lobits; 103 int bits_lost; 104 105 sum = dd_add(a, b); 106 107 /* 108 * If we are losing at least two bits of accuracy to denormalization, 109 * then the first lost bit becomes a round bit, and we adjust the 110 * lowest bit of sum.hi to make it a sticky bit summarizing all the 111 * bits in sum.lo. With the sticky bit adjusted, the hardware will 112 * break any ties in the correct direction. 113 * 114 * If we are losing only one bit to denormalization, however, we must 115 * break the ties manually. 116 */ 117 if (sum.lo != 0) { 118 EXTRACT_WORD64(hibits, sum.hi); 119 bits_lost = -((int)(hibits >> 52) & 0x7ff) - scale + 1; 120 if ((bits_lost != 1) ^ (int)(hibits & 1)) { 121 /* hibits += (int)copysign(1.0, sum.hi * sum.lo) */ 122 EXTRACT_WORD64(lobits, sum.lo); 123 hibits += 1 - (((hibits ^ lobits) >> 62) & 2); 124 INSERT_WORD64(sum.hi, hibits); 125 } 126 } 127 return (ldexp(sum.hi, scale)); 128 } 129 130 /* 131 * Compute a*b exactly, returning the exact result in a struct dd. We assume 132 * that both a and b are normalized, so no underflow or overflow will occur. 133 * The current rounding mode must be round-to-nearest. 134 */ 135 static inline struct dd 136 dd_mul(double a, double b) 137 { 138 static const double split = 0x1p27 + 1.0; 139 struct dd ret; 140 double ha, hb, la, lb, p, q; 141 142 p = a * split; 143 ha = a - p; 144 ha += p; 145 la = a - ha; 146 147 p = b * split; 148 hb = b - p; 149 hb += p; 150 lb = b - hb; 151 152 p = ha * hb; 153 q = ha * lb + la * hb; 154 155 ret.hi = p + q; 156 ret.lo = p - ret.hi + q + la * lb; 157 return (ret); 158 } 159 160 /* 161 * Fused multiply-add: Compute x * y + z with a single rounding error. 162 * 163 * We use scaling to avoid overflow/underflow, along with the 164 * canonical precision-doubling technique adapted from: 165 * 166 * Dekker, T. A Floating-Point Technique for Extending the 167 * Available Precision. Numer. Math. 18, 224-242 (1971). 168 * 169 * This algorithm is sensitive to the rounding precision. FPUs such 170 * as the i387 must be set in double-precision mode if variables are 171 * to be stored in FP registers in order to avoid incorrect results. 172 * This is the default on FreeBSD, but not on many other systems. 173 * 174 * Hardware instructions should be used on architectures that support it, 175 * since this implementation will likely be several times slower. 176 */ 177 DLLEXPORT double 178 fma(double x, double y, double z) 179 { 180 double xs, ys, zs, adj; 181 struct dd xy, r; 182 int oround; 183 int ex, ey, ez; 184 int spread; 185 186 /* 187 * Handle special cases. The order of operations and the particular 188 * return values here are crucial in handling special cases involving 189 * infinities, NaNs, overflows, and signed zeroes correctly. 190 */ 191 if (x == 0.0 || y == 0.0) 192 return (x * y + z); 193 if (z == 0.0) 194 return (x * y); 195 if (!isfinite(x) || !isfinite(y)) 196 return (x * y + z); 197 if (!isfinite(z)) 198 return (z); 199 200 xs = frexp(x, &ex); 201 ys = frexp(y, &ey); 202 zs = frexp(z, &ez); 203 oround = fegetround(); 204 spread = ex + ey - ez; 205 206 /* 207 * If x * y and z are many orders of magnitude apart, the scaling 208 * will overflow, so we handle these cases specially. Rounding 209 * modes other than FE_TONEAREST are painful. 210 */ 211 if (spread < -DBL_MANT_DIG) { 212 feraiseexcept(FE_INEXACT); 213 if (!isnormal(z)) 214 feraiseexcept(FE_UNDERFLOW); 215 switch (oround) { 216 case FE_TONEAREST: 217 return (z); 218 case FE_TOWARDZERO: 219 if ((x > 0.0) ^ (y < 0.0) ^ (z < 0.0)) 220 return (z); 221 else 222 return (nextafter(z, 0)); 223 case FE_DOWNWARD: 224 if ((x > 0.0) ^ (y < 0.0)) 225 return (z); 226 else 227 return (nextafter(z, -INFINITY)); 228 default: /* FE_UPWARD */ 229 if ((x > 0.0) ^ (y < 0.0)) 230 return (nextafter(z, INFINITY)); 231 else 232 return (z); 233 } 234 } 235 if (spread <= DBL_MANT_DIG * 2) 236 zs = ldexp(zs, -spread); 237 else 238 zs = copysign(DBL_MIN, zs); 239 240 fesetround(FE_TONEAREST); 241 242 /* 243 * Basic approach for round-to-nearest: 244 * 245 * (xy.hi, xy.lo) = x * y (exact) 246 * (r.hi, r.lo) = xy.hi + z (exact) 247 * adj = xy.lo + r.lo (inexact; low bit is sticky) 248 * result = r.hi + adj (correctly rounded) 249 */ 250 xy = dd_mul(xs, ys); 251 r = dd_add(xy.hi, zs); 252 253 spread = ex + ey; 254 255 if (r.hi == 0.0) { 256 /* 257 * When the addends cancel to 0, ensure that the result has 258 * the correct sign. 259 */ 260 fesetround(oround); 261 volatile double vzs = zs; /* XXX gcc CSE bug workaround */ 262 return (xy.hi + vzs + ldexp(xy.lo, spread)); 263 } 264 265 if (oround != FE_TONEAREST) { 266 /* 267 * There is no need to worry about double rounding in directed 268 * rounding modes. 269 */ 270 fesetround(oround); 271 adj = r.lo + xy.lo; 272 return (ldexp(r.hi + adj, spread)); 273 } 274 275 adj = add_adjusted(r.lo, xy.lo); 276 if (spread + ilogb(r.hi) > -1023) 277 return (ldexp(r.hi + adj, spread)); 278 else 279 return (add_and_denormalize(r.hi, adj, spread)); 280 } 281 282 #if (LDBL_MANT_DIG == 53) 283 __weak_reference(fma, fmal); 284 #endif 285