1140609Sdas/*- 2226245Sdas * Copyright (c) 2005-2011 David Schultz <das@FreeBSD.ORG> 3140609Sdas * All rights reserved. 4140609Sdas * 5140609Sdas * Redistribution and use in source and binary forms, with or without 6140609Sdas * modification, are permitted provided that the following conditions 7140609Sdas * are met: 8140609Sdas * 1. Redistributions of source code must retain the above copyright 9140609Sdas * notice, this list of conditions and the following disclaimer. 10140609Sdas * 2. Redistributions in binary form must reproduce the above copyright 11140609Sdas * notice, this list of conditions and the following disclaimer in the 12140609Sdas * documentation and/or other materials provided with the distribution. 13140609Sdas * 14140609Sdas * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND 15140609Sdas * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 16140609Sdas * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE 17140609Sdas * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE 18140609Sdas * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL 19140609Sdas * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS 20140609Sdas * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) 21140609Sdas * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT 22140609Sdas * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY 23140609Sdas * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF 24140609Sdas * SUCH DAMAGE. 25140609Sdas */ 26140609Sdas 27140609Sdas#include <sys/cdefs.h> 28140609Sdas__FBSDID("$FreeBSD: releng/10.2/lib/msun/src/s_fma.c 252170 2013-06-24 19:12:17Z eadler $"); 29140609Sdas 30140609Sdas#include <fenv.h> 31140609Sdas#include <float.h> 32140609Sdas#include <math.h> 33140609Sdas 34226371Sdas#include "math_private.h" 35226371Sdas 36140609Sdas/* 37226245Sdas * A struct dd represents a floating-point number with twice the precision 38226245Sdas * of a double. We maintain the invariant that "hi" stores the 53 high-order 39226245Sdas * bits of the result. 40226245Sdas */ 41226245Sdasstruct dd { 42226245Sdas double hi; 43226245Sdas double lo; 44226245Sdas}; 45226245Sdas 46226245Sdas/* 47226245Sdas * Compute a+b exactly, returning the exact result in a struct dd. We assume 48226245Sdas * that both a and b are finite, but make no assumptions about their relative 49226245Sdas * magnitudes. 50226245Sdas */ 51226245Sdasstatic inline struct dd 52226245Sdasdd_add(double a, double b) 53226245Sdas{ 54226245Sdas struct dd ret; 55226245Sdas double s; 56226245Sdas 57226245Sdas ret.hi = a + b; 58226245Sdas s = ret.hi - a; 59226245Sdas ret.lo = (a - (ret.hi - s)) + (b - s); 60226245Sdas return (ret); 61226245Sdas} 62226245Sdas 63226245Sdas/* 64226371Sdas * Compute a+b, with a small tweak: The least significant bit of the 65226371Sdas * result is adjusted into a sticky bit summarizing all the bits that 66226371Sdas * were lost to rounding. This adjustment negates the effects of double 67226371Sdas * rounding when the result is added to another number with a higher 68226371Sdas * exponent. For an explanation of round and sticky bits, see any reference 69226371Sdas * on FPU design, e.g., 70226371Sdas * 71226371Sdas * J. Coonen. An Implementation Guide to a Proposed Standard for 72226371Sdas * Floating-Point Arithmetic. Computer, vol. 13, no. 1, Jan 1980. 73226371Sdas */ 74226371Sdasstatic inline double 75226371Sdasadd_adjusted(double a, double b) 76226371Sdas{ 77226371Sdas struct dd sum; 78226371Sdas uint64_t hibits, lobits; 79226371Sdas 80226371Sdas sum = dd_add(a, b); 81226371Sdas if (sum.lo != 0) { 82226371Sdas EXTRACT_WORD64(hibits, sum.hi); 83226371Sdas if ((hibits & 1) == 0) { 84226371Sdas /* hibits += (int)copysign(1.0, sum.hi * sum.lo) */ 85226371Sdas EXTRACT_WORD64(lobits, sum.lo); 86226371Sdas hibits += 1 - ((hibits ^ lobits) >> 62); 87226371Sdas INSERT_WORD64(sum.hi, hibits); 88226371Sdas } 89226371Sdas } 90226371Sdas return (sum.hi); 91226371Sdas} 92226371Sdas 93226371Sdas/* 94226371Sdas * Compute ldexp(a+b, scale) with a single rounding error. It is assumed 95226371Sdas * that the result will be subnormal, and care is taken to ensure that 96226371Sdas * double rounding does not occur. 97226371Sdas */ 98226371Sdasstatic inline double 99226371Sdasadd_and_denormalize(double a, double b, int scale) 100226371Sdas{ 101226371Sdas struct dd sum; 102226371Sdas uint64_t hibits, lobits; 103226371Sdas int bits_lost; 104226371Sdas 105226371Sdas sum = dd_add(a, b); 106226371Sdas 107226371Sdas /* 108226371Sdas * If we are losing at least two bits of accuracy to denormalization, 109226371Sdas * then the first lost bit becomes a round bit, and we adjust the 110226371Sdas * lowest bit of sum.hi to make it a sticky bit summarizing all the 111226371Sdas * bits in sum.lo. With the sticky bit adjusted, the hardware will 112226371Sdas * break any ties in the correct direction. 113226371Sdas * 114226371Sdas * If we are losing only one bit to denormalization, however, we must 115226371Sdas * break the ties manually. 116226371Sdas */ 117226371Sdas if (sum.lo != 0) { 118226371Sdas EXTRACT_WORD64(hibits, sum.hi); 119226371Sdas bits_lost = -((int)(hibits >> 52) & 0x7ff) - scale + 1; 120252170Seadler if ((bits_lost != 1) ^ (int)(hibits & 1)) { 121226371Sdas /* hibits += (int)copysign(1.0, sum.hi * sum.lo) */ 122226371Sdas EXTRACT_WORD64(lobits, sum.lo); 123226371Sdas hibits += 1 - (((hibits ^ lobits) >> 62) & 2); 124226371Sdas INSERT_WORD64(sum.hi, hibits); 125226371Sdas } 126226371Sdas } 127226371Sdas return (ldexp(sum.hi, scale)); 128226371Sdas} 129226371Sdas 130226371Sdas/* 131226245Sdas * Compute a*b exactly, returning the exact result in a struct dd. We assume 132226245Sdas * that both a and b are normalized, so no underflow or overflow will occur. 133226245Sdas * The current rounding mode must be round-to-nearest. 134226245Sdas */ 135226245Sdasstatic inline struct dd 136226245Sdasdd_mul(double a, double b) 137226245Sdas{ 138226245Sdas static const double split = 0x1p27 + 1.0; 139226245Sdas struct dd ret; 140226245Sdas double ha, hb, la, lb, p, q; 141226245Sdas 142226245Sdas p = a * split; 143226245Sdas ha = a - p; 144226245Sdas ha += p; 145226245Sdas la = a - ha; 146226245Sdas 147226245Sdas p = b * split; 148226245Sdas hb = b - p; 149226245Sdas hb += p; 150226245Sdas lb = b - hb; 151226245Sdas 152226245Sdas p = ha * hb; 153226245Sdas q = ha * lb + la * hb; 154226245Sdas 155226245Sdas ret.hi = p + q; 156226245Sdas ret.lo = p - ret.hi + q + la * lb; 157226245Sdas return (ret); 158226245Sdas} 159226245Sdas 160226245Sdas/* 161140609Sdas * Fused multiply-add: Compute x * y + z with a single rounding error. 162140609Sdas * 163140609Sdas * We use scaling to avoid overflow/underflow, along with the 164140609Sdas * canonical precision-doubling technique adapted from: 165140609Sdas * 166140609Sdas * Dekker, T. A Floating-Point Technique for Extending the 167140609Sdas * Available Precision. Numer. Math. 18, 224-242 (1971). 168140609Sdas * 169140609Sdas * This algorithm is sensitive to the rounding precision. FPUs such 170140609Sdas * as the i387 must be set in double-precision mode if variables are 171140609Sdas * to be stored in FP registers in order to avoid incorrect results. 172140609Sdas * This is the default on FreeBSD, but not on many other systems. 173140609Sdas * 174143780Sdas * Hardware instructions should be used on architectures that support it, 175143780Sdas * since this implementation will likely be several times slower. 176140609Sdas */ 177140609Sdasdouble 178140609Sdasfma(double x, double y, double z) 179140609Sdas{ 180226371Sdas double xs, ys, zs, adj; 181226371Sdas struct dd xy, r; 182140609Sdas int oround; 183140609Sdas int ex, ey, ez; 184140609Sdas int spread; 185140609Sdas 186177875Sdas /* 187177875Sdas * Handle special cases. The order of operations and the particular 188177875Sdas * return values here are crucial in handling special cases involving 189177875Sdas * infinities, NaNs, overflows, and signed zeroes correctly. 190177875Sdas */ 191177875Sdas if (x == 0.0 || y == 0.0) 192177875Sdas return (x * y + z); 193140609Sdas if (z == 0.0) 194140609Sdas return (x * y); 195177875Sdas if (!isfinite(x) || !isfinite(y)) 196143230Sdas return (x * y + z); 197177875Sdas if (!isfinite(z)) 198177875Sdas return (z); 199140609Sdas 200140609Sdas xs = frexp(x, &ex); 201140609Sdas ys = frexp(y, &ey); 202140609Sdas zs = frexp(z, &ez); 203140609Sdas oround = fegetround(); 204140609Sdas spread = ex + ey - ez; 205140609Sdas 206140609Sdas /* 207140609Sdas * If x * y and z are many orders of magnitude apart, the scaling 208140609Sdas * will overflow, so we handle these cases specially. Rounding 209140609Sdas * modes other than FE_TONEAREST are painful. 210140609Sdas */ 211140609Sdas if (spread < -DBL_MANT_DIG) { 212140609Sdas feraiseexcept(FE_INEXACT); 213140609Sdas if (!isnormal(z)) 214140609Sdas feraiseexcept(FE_UNDERFLOW); 215140609Sdas switch (oround) { 216140609Sdas case FE_TONEAREST: 217140609Sdas return (z); 218140609Sdas case FE_TOWARDZERO: 219140609Sdas if (x > 0.0 ^ y < 0.0 ^ z < 0.0) 220140609Sdas return (z); 221140609Sdas else 222140609Sdas return (nextafter(z, 0)); 223140609Sdas case FE_DOWNWARD: 224140609Sdas if (x > 0.0 ^ y < 0.0) 225140609Sdas return (z); 226140609Sdas else 227140609Sdas return (nextafter(z, -INFINITY)); 228140609Sdas default: /* FE_UPWARD */ 229140609Sdas if (x > 0.0 ^ y < 0.0) 230140609Sdas return (nextafter(z, INFINITY)); 231140609Sdas else 232140609Sdas return (z); 233140609Sdas } 234140609Sdas } 235226371Sdas if (spread <= DBL_MANT_DIG * 2) 236226371Sdas zs = ldexp(zs, -spread); 237226371Sdas else 238226371Sdas zs = copysign(DBL_MIN, zs); 239140609Sdas 240140609Sdas fesetround(FE_TONEAREST); 241251024Sdas /* work around clang bug 8100 */ 242251024Sdas volatile double vxs = xs; 243140609Sdas 244226371Sdas /* 245226371Sdas * Basic approach for round-to-nearest: 246226371Sdas * 247226371Sdas * (xy.hi, xy.lo) = x * y (exact) 248226371Sdas * (r.hi, r.lo) = xy.hi + z (exact) 249226371Sdas * adj = xy.lo + r.lo (inexact; low bit is sticky) 250226371Sdas * result = r.hi + adj (correctly rounded) 251226371Sdas */ 252251024Sdas xy = dd_mul(vxs, ys); 253226245Sdas r = dd_add(xy.hi, zs); 254140609Sdas 255226601Sdas spread = ex + ey; 256226601Sdas 257226371Sdas if (r.hi == 0.0) { 258226371Sdas /* 259226371Sdas * When the addends cancel to 0, ensure that the result has 260226371Sdas * the correct sign. 261226371Sdas */ 262226371Sdas fesetround(oround); 263226371Sdas volatile double vzs = zs; /* XXX gcc CSE bug workaround */ 264226601Sdas return (xy.hi + vzs + ldexp(xy.lo, spread)); 265226371Sdas } 266226371Sdas 267226371Sdas if (oround != FE_TONEAREST) { 268143780Sdas /* 269226371Sdas * There is no need to worry about double rounding in directed 270226371Sdas * rounding modes. 271143780Sdas */ 272143780Sdas fesetround(oround); 273251024Sdas /* work around clang bug 8100 */ 274251024Sdas volatile double vrlo = r.lo; 275251024Sdas adj = vrlo + xy.lo; 276226371Sdas return (ldexp(r.hi + adj, spread)); 277143780Sdas } 278226371Sdas 279226371Sdas adj = add_adjusted(r.lo, xy.lo); 280226371Sdas if (spread + ilogb(r.hi) > -1023) 281226371Sdas return (ldexp(r.hi + adj, spread)); 282226371Sdas else 283226371Sdas return (add_and_denormalize(r.hi, adj, spread)); 284140609Sdas} 285143230Sdas 286143230Sdas#if (LDBL_MANT_DIG == 53) 287143264Sdas__weak_reference(fma, fmal); 288143230Sdas#endif 289