>From 0fde44871320b7b957d9f0d3c2e269d0fdad4651 Mon Sep 17 00:00:00 2001 From: Szabolcs Nagy Date: Wed, 12 Aug 2026 17:53:01 +0000 Subject: [PATCH 2/2] math: fmaf optimization and cleanup rewrite the algorithm after accumulated fixups, simplify the logic and optimize the generated code: - layout the code so the hot path has minimal clutter and the special cases are branched away efficiently (on both sf and hf). - handle underflow with simple fp ops instead of fenv. - use fp ops instead of fegetround so fmaf is a leaf function on hf targets and drop fenv.h. - reuse the residual t for inexact checks. subnormal and halfway handling is substantially different from the orig freebsd code so dropped references to it. fmaf.lo code size: x86_64: 723 -> 318 armhf: 486 -> 280 (v7 thumb, no vfma op) arm: 536 -> 544 (no fenv support, so uflow handling was ifdefed) --- src/math/fmaf.c | 117 ++++++++++++------------------------------------ 1 file changed, 29 insertions(+), 88 deletions(-) diff --git a/src/math/fmaf.c b/src/math/fmaf.c index e14dd194..e8b3a590 100644 --- a/src/math/fmaf.c +++ b/src/math/fmaf.c @@ -1,102 +1,43 @@ -/* origin: FreeBSD /usr/src/lib/msun/src/s_fmaf.c */ -/*- - * Copyright (c) 2005-2011 David Schultz - * All rights reserved. - * - * Redistribution and use in source and binary forms, with or without - * modification, are permitted provided that the following conditions - * are met: - * 1. Redistributions of source code must retain the above copyright - * notice, this list of conditions and the following disclaimer. - * 2. Redistributions in binary form must reproduce the above copyright - * notice, this list of conditions and the following disclaimer in the - * documentation and/or other materials provided with the distribution. - * - * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND - * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE - * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE - * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE - * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL - * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS - * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) - * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT - * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY - * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF - * SUCH DAMAGE. - */ - -#include #include #include -/* - * 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. - */ float fmaf(float x, float y, float z) { - #pragma STDC FENV_ACCESS ON - double xy, result; + double xy, r; union {double f; uint64_t i;} u; - int e, halfway; + int e, halfway, subnormal; xy = (double)x * y; - result = xy + z; - u.f = result; + u.f = r = xy + z; e = u.i>>52 & 0x7ff; halfway = (u.i & 0x1fffffff) == 0x10000000; - - /* subnormal range */ - if (e < 0x3ff-126 && e >= 0x3ff-149) { - /* fix halfway for subnormals */ - uint64_t m = 1; - m <<= 52 + 0x3ff-149 - e; - halfway = (u.i & m-1) == m/2; - } - - /* Common case: The double precision result is fine. */ - if (!halfway || - e == 0x7ff || /* NaN */ - (result - xy == z && result - z == xy) || /* exact */ - fegetround() != FE_TONEAREST) /* not round-to-nearest */ - { - /* - underflow may not be raised correctly, example: - fmaf(0x1p-120f, 0x1p-120f, 0x1p-149f) - */ -#if defined(FE_INEXACT) && defined(FE_UNDERFLOW) - if (e < 0x3ff-126 && e >= 0x3ff-149 && fetestexcept(FE_INEXACT)) { - feclearexcept(FE_INEXACT); - /* TODO: gcc and clang bug workaround */ - volatile float vz = z; - result = xy + vz; - if (fetestexcept(FE_INEXACT)) - feraiseexcept(FE_UNDERFLOW); + subnormal = e < 0x3ff-126 && e >= 0x3ff-149; + if (halfway || subnormal) { + /* subnormal or halfway normal or nan */ + int neg = u.i >> 63; + double t; + if (e == 0x7ff) return (float)r; /* nan */ + /* r + t == exact x*y+z */ + t = neg == (z > xy) ? xy - r + z : z - r + xy; + if (subnormal) { + uint64_t m = 1; + volatile float v; + /* ensure underflow is raised for inexact subnormal */ + v = (float)t; + /* check for halfway subnormal */ + m <<= 52 + 0x3ff-149 - e; + halfway = (u.i & m-1) == m/2; + if (!halfway) return (float)r; + } + /* big-t == big+t detects nearest rounding without calls */ + if (t && 0x1p999-t == 0x1p999+t) { + /* nearest rounding, inexact halfway case */ + if (neg == (t < 0)) + u.i++; else - feraiseexcept(FE_INEXACT); + u.i--; + r = u.f; } -#endif - z = result; - return z; } - - /* - * 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. - */ - double err; - int neg = u.i >> 63; - if (neg == (z > xy)) - err = xy - result + z; - else - err = z - result + xy; - if (neg == (err < 0)) - u.i++; - else - u.i--; - z = u.f; - return z; + return (float)r; } -- 2.52.0