>From ba9080aa2932becc68f030e833bd7189449e59bd Mon Sep 17 00:00:00 2001 From: Szabolcs Nagy Date: Tue, 11 Aug 2026 19:00:25 +0000 Subject: [PATCH v4 2/2] math: fmaf rewrite the new code uses that for all real |x| in [0x1p-999,0x1p999] y = (float)x is the same as r = (double)x t = x - r if (t!=0 && r.bits%2==0) r.bits += (r<0)==(t<0) ? 1 : -1 y = (float)r in all rounding modes, with the same fenv effects. this can be interpreted as a round to odd adjustment[1]. in fmaf t is computed with a fast2sum variant. - optimized common case. - implicit uflow handling instead of fenv calls. - no fegetround check. - no api calls on hf targets. - fixed missed uflow on targets that signal it before rounding: fmaf(-0x1p-100f, 0x1p-100f, 0x1p-126f) - removed freebsd code references and comments. - fmaf.o code size vs before the halfway subnormal fix: x86_64: 514 -> 210 armhf: 304 -> 156 (v7 thumb, no vfma op) arm: 400 -> 348 (soft float, no fenv) round to odd paper (suggested by Sergey Davidoff): [1] S. Boldo et al., Emulation of a FMA and correctly-rounded sums: proved algorithms using rounding to odd, 2008 --- src/math/fmaf.c | 116 +++++++++--------------------------------------- 1 file changed, 22 insertions(+), 94 deletions(-) diff --git a/src/math/fmaf.c b/src/math/fmaf.c index e14dd194..0bc35ec8 100644 --- a/src/math/fmaf.c +++ b/src/math/fmaf.c @@ -1,102 +1,30 @@ -/* 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; - union {double f; uint64_t i;} u; - int e, halfway; - - xy = (double)x * y; - result = xy + z; - u.f = result; - 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); - else - feraiseexcept(FE_INEXACT); + double xy = (double)x * y; + union {double r; uint64_t i;} u = {xy + z}; + int e = u.i>>52 & 0x7ff; + /* covers |r| > 0x1p-126 halfway cases (may round incorrectly) */ + int halfway = (u.i & 0x1fffffff) == 0x10000000; + /* covers tiny inexact (may miss uflow) and tiny halfway cases */ + int tiny = e <= 0x3ff-126 && e >= 0x3ff-149; + if (!halfway && !tiny) + /* common case, optimization only */ + return (float)u.r; + if (e != 0x7ff) { + /* r+t == x*y+z exactly in nearest rounding, otherwise + * rounding does not affect t!=0, t<0 and inexact flag. */ + int s = u.i >> 63; + double t = s == (xy < z) ? xy - u.r + z : z - u.r + xy; + if (t) { + /* adjust r toward r+t if r%2==0 (round to odd). + * may be needed if r is halfway or exact float, + * wrong if r becomes halfway or exact float. */ + u.i -= s ^ (t<0); + u.i |= 1; } -#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)u.r; } -- 2.52.0