From nobody Fri Sep 11 17:53:42 2026 X-Original-To: dev-commits-src-main@mlmmj.nyi.freebsd.org Received: from mx1.freebsd.org (mx1.freebsd.org [IPv6:2610:1c1:1:606c::19:1]) by mlmmj.nyi.freebsd.org (Postfix) with ESMTP id 4hhMbq0FTBz6sjYB for ; Fri, 11 Sep 2026 17:53:43 +0000 (UTC) (envelope-from git@FreeBSD.org) Received: from mxrelay.nyi.freebsd.org (mxrelay.nyi.freebsd.org [IPv6:2610:1c1:1:606c::19:3]) (using TLSv1.3 with cipher TLS_AES_256_GCM_SHA384 (256/256 bits) key-exchange x25519 server-signature RSA-PSS (4096 bits) server-digest SHA256 client-signature RSA-PSS (4096 bits) client-digest SHA256) (Client CN "mxrelay.nyi.freebsd.org", Issuer "YR2" (not verified)) by mx1.freebsd.org (Postfix) with ESMTPS id 4hhMbp5p5mz4tjt for ; Fri, 11 Sep 2026 17:53:42 +0000 (UTC) (envelope-from git@FreeBSD.org) DKIM-Signature: v=1; a=rsa-sha256; c=relaxed/relaxed; d=freebsd.org; s=dkim; t=1789149222; h=from:from:reply-to:subject:subject:date:date:message-id:message-id: to:to:cc:cc:mime-version:mime-version:content-type:content-type: content-transfer-encoding:content-transfer-encoding; bh=0lF6RZRb9tNLOQKFWjSWE5tUohH1JC2FA4if+76dkjA=; b=YShtsVjWw0KINHK1OA7e5ujeNg4KQfqEeQGKInNfXr9j3mXxFkkBmwjNMcd62NvS4AvAfi 4UuS9QxTqPKLpV8wgQMZ0Mvf8OFv6bVVugKbw4kcjdYkhP9y7iT2rjAw3AqyLTagnsVT9Y 3Ry4+WF4lZQAWgrYoNQjymUHs4t1+ngZhmGpekBFEEyJbs/uy1Q4kztyKL1tn1LhDoknI3 vUcY5YFO1LlbF5RhvImkbJaR0G5UEi8bllPp7aqEATn5+twvi1Uqf+s//MXge7bL41Sh9B 6ylBbGF9QTrXm7+kIUZfDPINM+X7/67gt02WxP92VBO62U9j6mv3mBQ8W4FKZg== ARC-Seal: i=1; a=rsa-sha256; d=freebsd.org; s=dkim; cv=none; t=1789149222; b=slDKGp4t+Akwg4VJXysKEe28sgAkF1w9IRWXBWxDul/Bu1qJRT6BJ0xuVLdgtuCVbpCCCo hsgqxiXlgwHlpnoHpzPW4cSRm4hNpp2UPdahv+zN8cjZriTrH0hcdsMof1bmUxmKzF9SK3 Ljgyba2xMYjlWJ9MsxUzP3v74SQ3RarZe4/34+CKt7Zg+ItaDOchYfdeCv+ZgX0OpQzLu/ rwl1KqcY9jn/xzgmbhpIq+/5dXILT6sRr1qmhD35N2ZmznWEIsaU/oHUphos9etgp393xN aTVGCfojR5Efd/nKprNumYGXv8q4s34uQmgECTMsrq7of4rVNIVJC2EwL4L7qQ== ARC-Message-Signature: i=1; a=rsa-sha256; c=relaxed/relaxed; d=freebsd.org; s=dkim; t=1789149222; h=from:from:reply-to:subject:subject:date:date:message-id:message-id: to:to:cc:cc:mime-version:mime-version:content-type:content-type: content-transfer-encoding:content-transfer-encoding; bh=0lF6RZRb9tNLOQKFWjSWE5tUohH1JC2FA4if+76dkjA=; b=GfrNANc04xnvQdJlWqkhlnPOrHZt9iink6EbNN0fgdVn1Q/zwlps1h9ZiBq+VCN4wQj+tq qLq4b98XImPbmAGzrmFhF72DM/ainUQgDAu3JWMRaY9k+j5i7cif6sNB4OY2qfhi5KXqpo rvcImhTBCEBhv3PmTcluhR+50y9C1g1VkUd2PZcnwOYbrCaT3z5mWMrddXqPG43tHu/OPL PjO0k+se8itjnDVWKWpITyIPw+TUrHmSAh4H3QcT08NcC0/1QRnEWNmVymUF9eLwOcJd4u MPixGe2pRQNI7Oc+lBWoc5/p2wQACBRJZRLPVaaOfhEc/e0Dg+1J46SKJ+E6NA== ARC-Authentication-Results: i=1; mx1.freebsd.org; none Received: from gitrepo.freebsd.org (gitrepo.freebsd.org [IPv6:2610:1c1:1:6068::e6a:5]) by mxrelay.nyi.freebsd.org (Postfix) with ESMTP id 4hhMbp4tDJzNPN for ; Fri, 11 Sep 2026 17:53:42 +0000 (UTC) (envelope-from git@FreeBSD.org) Received: from git (uid 1279) (envelope-from git@FreeBSD.org) id 459af by gitrepo.freebsd.org (DragonFly Mail Agent v0.13+ on gitrepo.freebsd.org); Fri, 11 Sep 2026 17:53:42 +0000 To: src-committers@FreeBSD.org, dev-commits-src-all@FreeBSD.org, dev-commits-src-main@FreeBSD.org Cc: Steve Kargl From: Konstantin Belousov Subject: git: b08e6a3882b4 - main - lib/msun: use the same algorithm for sw fmaf(3) as for fma(3) List-Id: Commit messages for the main branch of the src repository List-Archive: https://lists.freebsd.org/archives/dev-commits-src-main List-Help: List-Post: List-Subscribe: List-Unsubscribe: X-BeenThere: dev-commits-src-main@freebsd.org Sender: owner-dev-commits-src-main@FreeBSD.org List-Id: List-Post: List-Help: List-Subscribe: List-Unsubscribe: List-Owner: Precedence: list MIME-Version: 1.0 Content-Type: text/plain; charset=utf-8 Content-Transfer-Encoding: 8bit X-Git-Committer: kib X-Git-Repository: src X-Git-Refname: refs/heads/main X-Git-Reftype: branch X-Git-Commit: b08e6a3882b44a09fbda84dd315a431c7a585219 Auto-Submitted: auto-generated Date: Fri, 11 Sep 2026 17:53:42 +0000 Message-Id: <6aa44026.459af.718e42a6@gitrepo.freebsd.org> The branch main has been updated by kib: URL: https://cgit.FreeBSD.org/src/commit/?id=b08e6a3882b44a09fbda84dd315a431c7a585219 commit b08e6a3882b44a09fbda84dd315a431c7a585219 Author: Steve Kargl AuthorDate: 2026-09-10 23:18:58 +0000 Commit: Konstantin Belousov CommitDate: 2026-09-11 17:53:08 +0000 lib/msun: use the same algorithm for sw fmaf(3) as for fma(3) This fixes rounding at the last bit for subnormals. PR: 298260 Reviewed by: kib MFC after: 1 week Differential revision: https://reviews.freebsd.org/D59579 --- lib/msun/src/s_fmaf.c | 254 +++++++++++++++++++++++++++++++++++++++++++++----- 1 file changed, 230 insertions(+), 24 deletions(-) diff --git a/lib/msun/src/s_fmaf.c b/lib/msun/src/s_fmaf.c index 666d2bb26d91..27fbe462a3d8 100644 --- a/lib/msun/src/s_fmaf.c +++ b/lib/msun/src/s_fmaf.c @@ -27,17 +27,132 @@ */ #include +#include /* Needed for FLT_MANT_DIG. */ #include "math.h" #include "math_private.h" +#pragma STDC FENV_ACCESS ON + #ifdef USE_BUILTIN_FMAF float -fmaf(float x, float y, float z) +#ifdef _RENAME_FMAF +fmaf_sw +#else +fmaf +#endif +(float x, float y, float z) { return (__builtin_fmaf(x, y, z)); } #else + +#define _CC (0x1p12F + 1) /* Needed by _SLOW2SUM() and _MUL() below */ + +/* + * A struct dd represents a floating-point number with twice the precision + * of a float. We maintain the invariant that "hi" stores the 24 high-order + * bits of the result. + */ +struct dd { + float hi; + float lo; +}; + +/* + * Compute a+b exactly, returning the exact result in a struct dd. We assume + * that both a and b are finite, but make no assumptions about their relative + * magnitudes. + */ +static inline struct dd +dd_add(float a, float b) +{ + struct dd ret; + _SLOW2SUM(a, b, ret.hi, ret.lo); + return (ret); +} + +/* + * Compute a+b, with a small tweak: The least significant bit of the + * result is adjusted into a sticky bit summarizing all the bits that + * were lost to rounding. This adjustment negates the effects of double + * rounding when the result is added to another number with a higher + * exponent. For an explanation of round and sticky bits, see any reference + * on FPU design, e.g., + * + * J. Coonen. An Implementation Guide to a Proposed Standard for + * Floating-Point Arithmetic. Computer, vol. 13, no. 1, Jan 1980. + */ +static inline float +add_adjusted(float a, float b) +{ + struct dd sum; + uint32_t hibits, lobits; + + sum = dd_add(a, b); + if (sum.lo != 0) { + GET_FLOAT_WORD(hibits, sum.hi); + if ((hibits & 1) == 0) { + /* hibits += (int)copysignf(1.0f, sum.hi * sum.lo); */ + GET_FLOAT_WORD(lobits, sum.lo); + hibits += 1 - ((hibits ^ lobits) >> 30); + SET_FLOAT_WORD(sum.hi, hibits); + } + } + return (sum.hi); +} + +/* + * Compute ldexp(a+b, scale) with a single rounding error. It is assumed + * that the result will be subnormal, and care is taken to ensure that + * double rounding does not occur. + */ + +static inline float +add_and_denormalize(float a, float b, int scale) +{ + struct dd sum; + uint32_t hibits, lobits; + int bits_lost; + + sum = dd_add(a, b); + + /* + * If we are losing at least two bits of accuracy to denormalization, + * then the first lost bit becomes a round bit, and we adjust the + * lowest bit of sum.hi to make it a sticky bit summarizing all the + * bits in sum.lo. With the sticky bit adjusted, the hardware will + * break any ties in the correct direction. + * + * If we are losing only one bit to denormalization, however, we must + * break the ties manually. + */ + if (sum.lo != 0) { + GET_FLOAT_WORD(hibits, sum.hi); + bits_lost = -((int)(hibits >> 24) & 0x7f) - scale + 1; + if ((bits_lost != 1) ^ (int)(hibits & 1)) { + /* hibits += (int)copysign(1.0, sum.hi * sum.lo) */ + GET_FLOAT_WORD(lobits, sum.lo); + hibits += 1 - (((hibits ^ lobits) >> 30) & 2); + SET_FLOAT_WORD(sum.hi, hibits); + } + } + return (ldexpf(sum.hi, scale)); +} + +/* + * Compute a*b exactly, returning the exact result in a struct dd. We assume + * that both a and b are normalized, so no underflow or overflow will occur. + * The current rounding mode must be round-to-nearest. + */ +static inline struct dd +dd_mul(float a, float b) +{ + struct dd ret; + _MUL(a, b, ret.hi, ret.lo); + return (ret); +} + #ifdef _RENAME_FMAF float fmaf_sw(float, float, float); #endif @@ -45,9 +160,19 @@ float fmaf_sw(float, float, float); /* * Fused multiply-add: Compute x * y + z with a single rounding error. * - * A double has more than twice as much precision than a float, so - * direct double-precision arithmetic suffices, except where double - * rounding occurs. + * We use scaling to avoid overflow/underflow, along with the + * canonical precision-doubling technique adapted from: + * + * Dekker, T. A Floating-Point Technique for Extending the + * Available Precision. Numer. Math. 18, 224-242 (1971). + * + * This algorithm is sensitive to the rounding precision. FPUs such + * as the i387 must be set in double-precision mode if variables are + * to be stored in FP registers in order to avoid incorrect results. + * This is the default on FreeBSD, but not on many other systems. + * + * Hardware instructions should be used on architectures that support it, + * since this implementation will likely be several times slower. */ float #ifdef _RENAME_FMAF @@ -57,29 +182,110 @@ fmaf #endif (float x, float y, float z) { - double xy, result; - uint32_t hr, lr; - - xy = (double)x * y; - result = xy + z; - EXTRACT_WORDS(hr, lr, result); - /* Common case: The double precision result is fine. */ - if ((lr & 0x1fffffff) != 0x10000000 || /* not a halfway case */ - (hr & 0x7ff00000) == 0x7ff00000 || /* NaN */ - result - xy == z || /* exact */ - fegetround() != FE_TONEAREST) /* not round-to-nearest */ - return (result); + float xs, ys, zs, adj; + struct dd xy, r; + int oround; + int ex, ey, ez; + int spread; /* - * If result is inexact, and exactly halfway between two float values, - * we need to adjust the low-order bit in the direction of the error. + * Handle special cases. The order of operations and the particular + * return values here are crucial in handling special cases involving + * infinities, NaNs, overflows, and signed zeroes correctly. */ - fesetround(FE_TOWARDZERO); - volatile double vxy = xy; /* XXX work around gcc CSE bug */ - double adjusted_result = vxy + z; + if (x == 0 || y == 0) + return (x * y + z); + if (z == 0) + return (x * y); + if (!isfinite(x) || !isfinite(y)) + return (x * y + z); + if (!isfinite(z)) + return (z); + + xs = frexpf(x, &ex); + ys = frexpf(y, &ey); + zs = frexpf(z, &ez); + oround = fegetround(); + spread = ex + ey - ez; + + /* + * If x * y and z are many orders of magnitude apart, the scaling + * will overflow, so we handle these cases specially. Rounding + * modes other than FE_TONEAREST are painful. + */ + if (spread < -FLT_MANT_DIG) { + feraiseexcept(FE_INEXACT); + if (!isnormal(z)) + feraiseexcept(FE_UNDERFLOW); + switch (oround) { + case FE_TONEAREST: + return (z); + case FE_TOWARDZERO: + if ((x > 0) ^ (y < 0) ^ (z < 0)) + return (z); + else + return (nextafterf(z, 0)); + case FE_DOWNWARD: + if ((x > 0) ^ (y < 0)) + return (z); + else + return (nextafterf(z, -INFINITY)); + default: /* FE_UPWARD */ + if ((x > 0) ^ (y < 0)) + return (nextafterf(z, INFINITY)); + else + return (z); + } + } + if (spread <= FLT_MANT_DIG * 2) + zs = ldexpf(zs, -spread); + else + zs = copysignf(FLT_MIN, zs); + fesetround(FE_TONEAREST); - if (result == adjusted_result) - SET_LOW_WORD(adjusted_result, lr + 1); - return (adjusted_result); + /* work around clang issue #8472 */ + volatile float vxs = xs; + + /* + * Basic approach for round-to-nearest: + * + * (xy.hi, xy.lo) = x * y (exact) + * (r.hi, r.lo) = xy.hi + z (exact) + * adj = xy.lo + r.lo (inexact; low bit is sticky) + * result = r.hi + adj (correctly rounded) + */ + xy = dd_mul(vxs, ys); + r = dd_add(xy.hi, zs); + + spread = ex + ey; + + if (r.hi == 0 && xy.lo == 0) { + /* + * When the addends cancel to 0, ensure that the result has + * the correct sign. + */ + fesetround(oround); + volatile float vzs = zs; /* XXX gcc CSE bug workaround */ + return (xy.hi + vzs); + } + + if (oround != FE_TONEAREST) { + /* + * There is no need to worry about double rounding in directed + * rounding modes. + */ + fesetround(oround); + /* work around clang issue #8472 */ + volatile float vrlo = r.lo; + adj = vrlo + xy.lo; + return (ldexpf(r.hi + adj, spread)); + } + + adj = add_adjusted(r.lo, xy.lo); + + if (spread + ilogbf(r.hi) > -127) + return (ldexpf(r.hi + adj, spread)); + else + return (add_and_denormalize(r.hi, adj, spread)); } #endif /* !USE_BUILTIN_FMAF */