From nobody Mon Sep 14 00:45:10 2026 X-Original-To: dev-commits-src-branches@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 4hjmdg4Dw2z6rKkL for ; Mon, 14 Sep 2026 00:45:11 +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 4hjmdg0tjRz4pnQ for ; Mon, 14 Sep 2026 00:45:11 +0000 (UTC) (envelope-from git@FreeBSD.org) DKIM-Signature: v=1; a=rsa-sha256; c=relaxed/relaxed; d=freebsd.org; s=dkim; t=1789346711; 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=lVKK4Iuta8vYiBrPPwGYlVgWpPs7dy1SpWnkbQovns0=; b=EL/N4aeJ1vPk8sjxI752HB+buA0a+Cs5ahR580MQHAlQ/lfyFVyn+VBrU+cWCfaNQUZYMQ cC/WrQmrxtVa5mxGLE3MHcPcIZxct2aSTFLBbnjkptSbha9zphi4y71DG+DFYg6i/zR6GW f9dIeq3Qh0WSTM+1J1ntyWcnsx4onDDCSH6AJifmT8KmJQqsj7YM/RTDs0m4xPWITJcpPz 7Ut/1s15Fy7RHu91+GlcLIQVcBYFlobWDZRtsZN2cmYZJdO+QLE3MCNVTlBM04WsgfuzOd Aas+EBHD1i1vBgB5AWd5Yc8yzAEejWFPYO817lYSCXdoV88PxedVqY3ZpieKIA== ARC-Seal: i=1; a=rsa-sha256; d=freebsd.org; s=dkim; cv=none; t=1789346711; b=aK7loUyzmIzd4u+/PMXNQRT4+zYyP8tptdcCLtevL8nyL+4dOwzgErFFhMnE2DgQabO/Cx zYyqLHa3CPof7rwk/ZTV6X0qcwUfNfkfgoVzKxCYzpCIKCGxHPINy9Cdeffm1i3SG+mpLI j1/P8Z0+VGpJ8TdRYlcQm48LeeMdrqJLx9lR4pJFCg6WHk0H8PaTZC4R4FS+LpLEcIjjCC VqZQEZexzxmbLm79YaPTdpR2IJ26M1Gr4yvZIDPNKDpaBlBw9lRRLWMnmfiNzOtwB1+FlL 7uscHBLtsWfn+FHQprs9UI1/2sl5lY+SJMHTKeSvZlitwlscZxKULWInWVc9fQ== ARC-Message-Signature: i=1; a=rsa-sha256; c=relaxed/relaxed; d=freebsd.org; s=dkim; t=1789346711; 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=lVKK4Iuta8vYiBrPPwGYlVgWpPs7dy1SpWnkbQovns0=; b=YYNhMcWtot5C8dCfkz0Vj4ig+lQtQiwH1vmjHahnZlF78QujhjgLuzY6NBt//qkXlSx+gS SKFrj+egOAVghupinheil+SqDLOnGyQgDnzp7y8W46568olmUTXAA6kKTk5aNL0hjV6AZH jWJIM9We+jgyrT5h8Vt6zijVZa9VeCK4RIqRtcAe+WzahIHUIGKd1aehs2Q/q7IYXSHuGo 5fJvXPERcFQ3be35n3fWzA+vnjqbZ46rKqJpbUwIjTTtWzLLHpqnxU8hZAU5Uj5vHwxpvw 4dqAwqOZ8VJYp50EjunKBd021iVN+MBGFXXyB5r7gRPC60vkraSq+kwPbOpZbQ== 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 4hjmdf7388z1631 for ; Mon, 14 Sep 2026 00:45:10 +0000 (UTC) (envelope-from git@FreeBSD.org) Received: from git (uid 1279) (envelope-from git@FreeBSD.org) id 1a364 by gitrepo.freebsd.org (DragonFly Mail Agent v0.13+ on gitrepo.freebsd.org); Mon, 14 Sep 2026 00:45:10 +0000 To: src-committers@FreeBSD.org, dev-commits-src-all@FreeBSD.org, dev-commits-src-branches@FreeBSD.org Cc: Steve Kargl From: Konstantin Belousov Subject: git: 4564e1f715d8 - stable/15 - lib/msun: use the same algorithm for sw fmaf(3) as for fma(3) List-Id: Commits to the stable branches of the FreeBSD src repository List-Archive: https://lists.freebsd.org/archives/dev-commits-src-branches List-Help: List-Post: List-Subscribe: List-Unsubscribe: X-BeenThere: dev-commits-src-branches@freebsd.org Sender: owner-dev-commits-src-branches@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/stable/15 X-Git-Reftype: branch X-Git-Commit: 4564e1f715d808ebfa9199879b4611ed6f4856eb Auto-Submitted: auto-generated Date: Mon, 14 Sep 2026 00:45:10 +0000 Message-Id: <6aa74396.1a364.263c08b7@gitrepo.freebsd.org> The branch stable/15 has been updated by kib: URL: https://cgit.FreeBSD.org/src/commit/?id=4564e1f715d808ebfa9199879b4611ed6f4856eb commit 4564e1f715d808ebfa9199879b4611ed6f4856eb Author: Steve Kargl AuthorDate: 2026-09-10 23:18:58 +0000 Commit: Konstantin Belousov CommitDate: 2026-09-14 00:44:42 +0000 lib/msun: use the same algorithm for sw fmaf(3) as for fma(3) PR: 298260 (cherry picked from commit b08e6a3882b44a09fbda84dd315a431c7a585219) --- 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 */