blob: 1fdaf2571af4e47c78491488b6f42925615f2a45 (
plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
|
From 282b1cd26649d69de038111f5876853df6ddc345 Mon Sep 17 00:00:00 2001
From: Szabolcs Nagy <nsz@port70.net>
Date: Sun, 1 Apr 2018 20:02:01 +0000
Subject: [PATCH] fix fmaf wrong result
if double precision r=x*y+z is not a half way case between two single
precision floats or it is an exact result then fmaf returns (float)r.
however the exactness check was wrong when |x*y| < |z| and could cause
incorrectly rounded result in nearest rounding mode when r is a half
way case.
fmaf(-0x1.26524ep-54, -0x1.cb7868p+11, 0x1.d10f5ep-29)
was incorrectly rounded up to 0x1.d117ap-29 instead of 0x1.d1179ep-29.
(exact result is 0x1.d1179efffffffecp-29, r is 0x1.d1179fp-29)
---
src/math/fmaf.c | 2 +-
1 file changed, 1 insertion(+), 1 deletion(-)
diff --git a/src/math/fmaf.c b/src/math/fmaf.c
index aa57feb6..80f5cd8a 100644
--- a/src/math/fmaf.c
+++ b/src/math/fmaf.c
@@ -50,7 +50,7 @@ float fmaf(float x, float y, float z)
/* Common case: The double precision result is fine. */
if ((u.i & 0x1fffffff) != 0x10000000 || /* not a halfway case */
e == 0x7ff || /* NaN */
- result - xy == z || /* exact */
+ (result - xy == z && result - z == xy) || /* exact */
fegetround() != FE_TONEAREST) /* not round-to-nearest */
{
/*
--
2.17.0
|
