From: Szabolcs Nagy <nsz@port70.net>
To: Alba Pompeo <albapompeo@gmail.com>
Cc: musl@lists.openwall.com
Subject: Re: musl libmath (implied algorithm issue)
Date: Sun, 1 Apr 2018 22:42:27 +0200 [thread overview]
Message-ID: <20180401204227.GG4418@port70.net> (raw)
In-Reply-To: <CAJDAfTDD+x3zS0kbbNU=5XzGqmNTene1HmZiG44feDXpkb=s2A@mail.gmail.com>
[-- Attachment #1: Type: text/plain, Size: 497 bytes --]
* Alba Pompeo <albapompeo@gmail.com> [2018-04-01 10:54:54 -0300]:
> I found a commit elsewhere that implies musl's fma algorithm is wrong.
> You may be interested in investigating.
> https://git.eigenstate.org/npnth/mc.git/commit/?h=libmath&id=debd13efe12eaca072229e9859aaa6df8708a09f
thanks, it is a wrong result bug in fmaf, copied from freebsd.
near half way cases may get incorrectly rounded under some
circumstances, i didnt try to understand their fix, but
used a minimal change for now.
[-- Attachment #2: 0001-fix-fmaf-wrong-result.patch --]
[-- Type: text/x-diff, Size: 1418 bytes --]
From 2a38f5d89f7a69d03861a47a7f981136ee60630f 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)
this example was reported in
https://git.eigenstate.org/npnth/mc.git/commit/?h=libmath&id=debd13efe12eaca072229e9859aaa6df8708a09f
---
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.16.3
parent reply other threads:[~2018-04-01 20:42 UTC|newest]
Thread overview: expand[flat|nested] mbox.gz Atom feed
[parent not found: <CAJDAfTDD+x3zS0kbbNU=5XzGqmNTene1HmZiG44feDXpkb=s2A@mail.gmail.com>]
Reply instructions:
You may reply publicly to this message via plain-text email
using any one of the following methods:
* Save the following mbox file, import it into your mail client,
and reply-to-all from there: mbox
Avoid top-posting and favor interleaved quoting:
https://en.wikipedia.org/wiki/Posting_style#Interleaved_style
* Reply using the --to, --cc, and --in-reply-to
switches of git-send-email(1):
git send-email \
--in-reply-to=20180401204227.GG4418@port70.net \
--to=nsz@port70.net \
--cc=albapompeo@gmail.com \
--cc=musl@lists.openwall.com \
/path/to/YOUR_REPLY
https://kernel.org/pub/software/scm/git/docs/git-send-email.html
* If your mail client supports setting the In-Reply-To header
via mailto: links, try the mailto: link
Be sure your reply has a Subject: header at the top and a blank line
before the message body.
Code repositories for project(s) associated with this public inbox
https://git.vuxu.org/mirror/musl/
This is a public inbox, see mirroring instructions
for how to clone and mirror all data and code used for this inbox;
as well as URLs for NNTP newsgroup(s).