[<prev] [next>] [thread-next>] [day] [month] [year] [list]
Message-ID: <CABaUXi6TkZBYLYeF=mD5ew8+1eLYUUuGF8f_T64t9C_p5sgFcA@mail.gmail.com>
Date: Mon, 10 Aug 2026 16:11:19 +0100
From: Sergey Davidoff <shnatsel@...il.com>
To: musl@...ts.openwall.com
Subject: Bug in fmaf: subnormals are rounded incorrectly
Hi all,
The algorithm in src/math/fmaf.c (software FMA implementation for systems
without hardware support) rounds subnormal results incorrectly.
For example, these bit patterns reinterpreted (not cast) into floats:
x = 0x97000800
y = 0x1cfff001
z = 0x00010002
evaluate to 0x00010002 in musl fmaf, but hardware FMA instructions evaluate
to 0x00010001.
This was tested on x86_64 so I'm sure it's not running into x87 float
precision issues. The algorithm is simply incorrect. The root cause is this
line:
if ((u.i & 0x1fffffff) != 0x10000000 /* not a halfway case */
This only checks halfway cases for normal floats, but fails to consider
subnormals.
A correct algorithm that rounds subnormals accurately can be found in the 2008
paper "Emulation of FMA and correctly-rounded sums: proved algorithms using
rounding to odd" by Sylvie Boldo and Guillaume Melquiond. The paper is freely
available online but not linked here to avoid aggravating spam filters. It
includes a proof of correctness in Coq for their algorithm.
Cheers,
--
Shnatsel
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.
