[<prev] [next>] [<thread-prev] [thread-next>] [day] [month] [year] [list]
Message-ID: <20260816222336.GB23438@brightrain.aerifal.cx>
Date: Sun, 16 Aug 2026 18:23:36 -0400
From: Rich Felker <dalias@...c.org>
To: musl@...ts.openwall.com, Sergey Davidoff <shnatsel@...il.com>
Subject: Re: [PATCH v2] math: fmaf rewrite and bug fixes
On Sun, Aug 16, 2026 at 10:35:44PM +0200, Szabolcs Nagy wrote:
> i no longer think a separate bugfix patch is useful.
Unless you think it's wrong or incomplete in fixing the wrong-value
bug (wrong flags is far less of an issue), I'd like to go ahead and
apply the direct fix first. It's friendlier to bugfix-only backporting
which users may want to do for an old version they're using, and
documents the bug clearly and allows verifying any differences with
bisect.
> new code is simpler and similar to the mentioned
> round to odd approach (that proof does not apply,
> but the complexity is not too bad).
>
> a big chunk of the code is just optimization, i
> don't have real fmaf workload, so it is a guess
> that extra int arith & predictable branches help.
> From f71c18023a93fea8539e29380f4e123357f448a9 Mon Sep 17 00:00:00 2001
> From: Szabolcs Nagy <nsz@...t70.net>
> Date: Tue, 11 Aug 2026 19:00:25 +0000
> Subject: [PATCH v2] math: fmaf rewrite and bug fixes
>
> 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&0xfffffff)==0)
> r.bits += (r<0)==(t<0) ? 1 : -1
> y = (float)r
>
> in all rounding modes, with the same fenv effects. note: using
> r.bits&1 this can be interpreted as a round to odd adjustment,
> y = (float)ro(x), but testing more bits avoids adjustment in
> more cases. 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 a double rounding error of some subnormal halfway cases:
> fmaf(0x20201p-92f, 0x1fe01p-92f, 0x1p-130f)
> = (float)(0x1p-130 + 0x1.000000004p-150)
> was 0x1p-130 instead of 0x1.00002p-130.
> - 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:
> x86_64: 514 -> 221
> armhf: 304 -> 170 (v7 thumb, no vfma op)
> arm: 400 -> 380 (soft float, no fenv)
>
> subnormal double rounding bug, and round to odd reference [1]:
> Reported-by: Sergey Davidoff <shnatsel@...il.com>
>
> [1]: S. Boldo et al., Emulation of a FMA and correctly-rounded sums:
> proved algorithms using rounding to odd, 2008
> ---
> src/math/fmaf.c | 106 ++++++++++--------------------------------------
> 1 file changed, 21 insertions(+), 85 deletions(-)
>
> diff --git a/src/math/fmaf.c b/src/math/fmaf.c
> index 7c65acf1..2d2d4f48 100644
> --- a/src/math/fmaf.c
> +++ b/src/math/fmaf.c
> @@ -1,92 +1,28 @@
> -/* origin: FreeBSD /usr/src/lib/msun/src/s_fmaf.c */
> -/*-
> - * Copyright (c) 2005-2011 David Schultz <das@...eBSD.ORG>
> - * 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 <fenv.h>
> #include <math.h>
> #include <stdint.h>
>
> -/*
> - * 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;
> -
> - xy = (double)x * y;
> - result = xy + z;
> - u.f = result;
> - e = u.i>>52 & 0x7ff;
> - /* Common case: The double precision result is fine. */
> - if ((u.i & 0x1fffffff) != 0x10000000 || /* not a halfway case */
> - 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);
> - }
> -#endif
> - z = result;
> - return z;
> + 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 ((u.i & 0xfffffff) == 0 && e < 0x3ff+128) {
> + /* correct except for inf, nan and some (r&1)==1 cases.
> + r+t == x*y+z exactly in nearest rounding mode and in
> + other modes t!=0 and t<0 checks are not affected.
> + if t==0 then inexact is not signaled. */
> + 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 */
> + u.i += s == (t < 0) ? 1 : -1;
> }
Looks good. Small style nit: *'s on the continuation lines of comment
to be easily distinguishable from code.
Does this all sound ok?
Rich
Powered by blists - more mailing lists
Confused about mailing lists and their use?
Read about mailing lists on Wikipedia
and check out these
guidelines on proper formatting of your messages.
