1 /* $NetBSD: s_fma.c,v 1.7 2017/05/06 18:02:52 christos Exp $ */ 2 3 /*- 4 * Copyright (c) 2005-2011 David Schultz <das@FreeBSD.ORG> 5 * All rights reserved. 6 * 7 * Redistribution and use in source and binary forms, with or without 8 * modification, are permitted provided that the following conditions 9 * are met: 10 * 1. Redistributions of source code must retain the above copyright 11 * notice, this list of conditions and the following disclaimer. 12 * 2. Redistributions in binary form must reproduce the above copyright 13 * notice, this list of conditions and the following disclaimer in the 14 * documentation and/or other materials provided with the distribution. 15 * 16 * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND 17 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 18 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE 19 * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE 20 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL 21 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS 22 * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) 23 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT 24 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY 25 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF 26 * SUCH DAMAGE. 27 */ 28 29 #include <sys/cdefs.h> 30 #if 0 31 __FBSDID("$FreeBSD: src/lib/msun/src/s_fma.c,v 1.8 2011/10/21 06:30:43 das Exp $"); 32 #else 33 __RCSID("$NetBSD: s_fma.c,v 1.7 2017/05/06 18:02:52 christos Exp $"); 34 #endif 35 36 #include "namespace.h" 37 38 #include <machine/ieee.h> 39 #include <fenv.h> 40 #include <float.h> 41 #include <math.h> 42 43 #include "math_private.h" 44 45 #ifndef __HAVE_LONG_DOUBLE 46 __strong_alias(fmal, fma) 47 #endif 48 49 /* 50 * A struct dd represents a floating-point number with twice the precision 51 * of a double. We maintain the invariant that "hi" stores the 53 high-order 52 * bits of the result. 53 */ 54 struct dd { 55 double hi; 56 double lo; 57 }; 58 59 /* 60 * Compute a+b exactly, returning the exact result in a struct dd. We assume 61 * that both a and b are finite, but make no assumptions about their relative 62 * magnitudes. 63 */ 64 static inline struct dd 65 dd_add(double a, double b) 66 { 67 struct dd ret; 68 double s; 69 70 ret.hi = a + b; 71 s = ret.hi - a; 72 ret.lo = (a - (ret.hi - s)) + (b - s); 73 return (ret); 74 } 75 76 /* 77 * Compute a+b, with a small tweak: The least significant bit of the 78 * result is adjusted into a sticky bit summarizing all the bits that 79 * were lost to rounding. This adjustment negates the effects of double 80 * rounding when the result is added to another number with a higher 81 * exponent. For an explanation of round and sticky bits, see any reference 82 * on FPU design, e.g., 83 * 84 * J. Coonen. An Implementation Guide to a Proposed Standard for 85 * Floating-Point Arithmetic. Computer, vol. 13, no. 1, Jan 1980. 86 */ 87 static inline double 88 add_adjusted(double a, double b) 89 { 90 struct dd sum; 91 uint64_t hibits, lobits; 92 93 sum = dd_add(a, b); 94 if (sum.lo != 0) { 95 EXTRACT_WORD64(hibits, sum.hi); 96 if ((hibits & 1) == 0) { 97 /* hibits += (int)copysign(1.0, sum.hi * sum.lo) */ 98 EXTRACT_WORD64(lobits, sum.lo); 99 hibits += 1 - ((hibits ^ lobits) >> 62); 100 INSERT_WORD64(sum.hi, hibits); 101 } 102 } 103 return (sum.hi); 104 } 105 106 /* 107 * Compute ldexp(a+b, scale) with a single rounding error. It is assumed 108 * that the result will be subnormal, and care is taken to ensure that 109 * double rounding does not occur. 110 */ 111 static inline double 112 add_and_denormalize(double a, double b, int scale) 113 { 114 struct dd sum; 115 uint64_t hibits, lobits; 116 int bits_lost; 117 118 sum = dd_add(a, b); 119 120 /* 121 * If we are losing at least two bits of accuracy to denormalization, 122 * then the first lost bit becomes a round bit, and we adjust the 123 * lowest bit of sum.hi to make it a sticky bit summarizing all the 124 * bits in sum.lo. With the sticky bit adjusted, the hardware will 125 * break any ties in the correct direction. 126 * 127 * If we are losing only one bit to denormalization, however, we must 128 * break the ties manually. 129 */ 130 if (sum.lo != 0) { 131 EXTRACT_WORD64(hibits, sum.hi); 132 bits_lost = -((int)(hibits >> 52) & 0x7ff) - scale + 1; 133 if ((bits_lost != 1) ^ (int)(hibits & 1)) { 134 /* hibits += (int)copysign(1.0, sum.hi * sum.lo) */ 135 EXTRACT_WORD64(lobits, sum.lo); 136 hibits += 1 - (((hibits ^ lobits) >> 62) & 2); 137 INSERT_WORD64(sum.hi, hibits); 138 } 139 } 140 return (ldexp(sum.hi, scale)); 141 } 142 143 /* 144 * Compute a*b exactly, returning the exact result in a struct dd. We assume 145 * that both a and b are normalized, so no underflow or overflow will occur. 146 * The current rounding mode must be round-to-nearest. 147 */ 148 static inline struct dd 149 dd_mul(double a, double b) 150 { 151 static const double split = 0x1p27 + 1.0; 152 struct dd ret; 153 double ha, hb, la, lb, p, q; 154 155 p = a * split; 156 ha = a - p; 157 ha += p; 158 la = a - ha; 159 160 p = b * split; 161 hb = b - p; 162 hb += p; 163 lb = b - hb; 164 165 p = ha * hb; 166 q = ha * lb + la * hb; 167 168 ret.hi = p + q; 169 ret.lo = p - ret.hi + q + la * lb; 170 return (ret); 171 } 172 173 /* 174 * Fused multiply-add: Compute x * y + z with a single rounding error. 175 * 176 * We use scaling to avoid overflow/underflow, along with the 177 * canonical precision-doubling technique adapted from: 178 * 179 * Dekker, T. A Floating-Point Technique for Extending the 180 * Available Precision. Numer. Math. 18, 224-242 (1971). 181 * 182 * This algorithm is sensitive to the rounding precision. FPUs such 183 * as the i387 must be set in double-precision mode if variables are 184 * to be stored in FP registers in order to avoid incorrect results. 185 * This is the default on FreeBSD, but not on many other systems. 186 * 187 * Hardware instructions should be used on architectures that support it, 188 * since this implementation will likely be several times slower. 189 */ 190 double 191 fma(double x, double y, double z) 192 { 193 double xs, ys, zs, adj; 194 struct dd xy, r; 195 int oround; 196 int ex, ey, ez; 197 int spread; 198 199 /* 200 * Handle special cases. The order of operations and the particular 201 * return values here are crucial in handling special cases involving 202 * infinities, NaNs, overflows, and signed zeroes correctly. 203 */ 204 if (x == 0.0 || y == 0.0) 205 return (x * y + z); 206 if (z == 0.0) 207 return (x * y); 208 if (!isfinite(x) || !isfinite(y)) 209 return (x * y + z); 210 if (!isfinite(z)) 211 return (z); 212 213 xs = frexp(x, &ex); 214 ys = frexp(y, &ey); 215 zs = frexp(z, &ez); 216 oround = fegetround(); 217 spread = ex + ey - ez; 218 219 /* 220 * If x * y and z are many orders of magnitude apart, the scaling 221 * will overflow, so we handle these cases specially. Rounding 222 * modes other than FE_TONEAREST are painful. 223 */ 224 if (spread < -DBL_MANT_DIG) { 225 feraiseexcept(FE_INEXACT); 226 if (!isnormal(z)) 227 feraiseexcept(FE_UNDERFLOW); 228 switch (oround) { 229 case FE_TONEAREST: 230 return (z); 231 case FE_TOWARDZERO: 232 if ((x > 0.0) ^ (y < 0.0) ^ (z < 0.0)) 233 return (z); 234 else 235 return (nextafter(z, 0)); 236 case FE_DOWNWARD: 237 if ((x > 0.0) ^ (y < 0.0)) 238 return (z); 239 else 240 return (nextafter(z, -INFINITY)); 241 default: /* FE_UPWARD */ 242 if ((x > 0.0) ^ (y < 0.0)) 243 return (nextafter(z, INFINITY)); 244 else 245 return (z); 246 } 247 } 248 if (spread <= DBL_MANT_DIG * 2) 249 zs = ldexp(zs, -spread); 250 else 251 zs = copysign(DBL_MIN, zs); 252 253 fesetround(FE_TONEAREST); 254 255 /* 256 * Basic approach for round-to-nearest: 257 * 258 * (xy.hi, xy.lo) = x * y (exact) 259 * (r.hi, r.lo) = xy.hi + z (exact) 260 * adj = xy.lo + r.lo (inexact; low bit is sticky) 261 * result = r.hi + adj (correctly rounded) 262 */ 263 xy = dd_mul(xs, ys); 264 r = dd_add(xy.hi, zs); 265 266 spread = ex + ey; 267 268 if (r.hi == 0.0) { 269 /* 270 * When the addends cancel to 0, ensure that the result has 271 * the correct sign. 272 */ 273 fesetround(oround); 274 { 275 volatile double vzs = zs; /* XXX gcc CSE bug workaround */ 276 return (xy.hi + vzs + ldexp(xy.lo, spread)); 277 } 278 } 279 280 if (oround != FE_TONEAREST) { 281 /* 282 * There is no need to worry about double rounding in directed 283 * rounding modes. 284 */ 285 fesetround(oround); 286 adj = r.lo + xy.lo; 287 return (ldexp(r.hi + adj, spread)); 288 } 289 290 adj = add_adjusted(r.lo, xy.lo); 291 if (spread + ilogb(r.hi) > -1023) 292 return (ldexp(r.hi + adj, spread)); 293 else 294 return (add_and_denormalize(r.hi, adj, spread)); 295 } 296