|
|
Message-ID: <20261029215739.166761-1-bernie@codewiz.org>
Date: Fri, 30 Oct 2026 06:57:38 +0900
From: Bernie Innocenti <bernie@...ewiz.org>
To: musl@...ts.openwall.com
Cc: Bernie Innocenti <bernie@...ewiz.org>
Subject: [PATCH] math: fix fma(x,y,0) double rounding with excess precision
with FLT_EVAL_METHOD==2 (i386 x87, m68k with 68881) x*y is rounded
to long double first and then to double, e.g.
fma(0x1.27365a6456ba5p+0, 0x1.63fdaf7759a4fp+0, 0.0)
returned 0x1.9a84ea6842dccp+0 instead of 0x1.9a84ea6842dcdp+0.
add the zero z exactly on the integer path instead: then x*y is
rounded once in the current rounding mode with its own sign, which
keeps the fix of commit 9683bd62414604d3bd56cf6bd7be8f54aa31e7d3.
found by arczi84 (GitHub) reviewing a port of this code to newlib.
---
src/math/fma.c | 10 +++++++---
1 file changed, 7 insertions(+), 3 deletions(-)
diff --git a/src/math/fma.c b/src/math/fma.c
index adfadca8..072e8e5d 100644
--- a/src/math/fma.c
+++ b/src/math/fma.c
@@ -52,9 +52,13 @@ double fma(double x, double y, double z)
if (nx.e >= ZEROINFNAN || ny.e >= ZEROINFNAN)
return x*y + z;
if (nz.e >= ZEROINFNAN) {
- if (nz.e > ZEROINFNAN) /* z==0 */
- return x*y;
- return z;
+ if (nz.e == ZEROINFNAN)
+ return z;
+ /* z==0: x*y would round twice with FLT_EVAL_METHOD==2,
+ so add z exactly, x*y then rounds once with its own sign */
+ nz.m = 0;
+ nz.e = nx.e + ny.e;
+ nz.sign = nx.sign ^ ny.sign;
}
/* mul: r = x*y */
--
2.56.0
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.