NetBSD Problem Report #58742

From www@netbsd.org  Fri Oct 11 20:21:40 2024
Return-Path: <www@netbsd.org>
Received: from mail.netbsd.org (mail.netbsd.org [199.233.217.200])
	(using TLSv1.3 with cipher TLS_AES_256_GCM_SHA384 (256/256 bits)
	 key-exchange X25519 server-signature RSA-PSS (2048 bits)
	 client-signature RSA-PSS (2048 bits))
	(Client CN "mail.NetBSD.org", Issuer "mail.NetBSD.org CA" (not verified))
	by mollari.NetBSD.org (Postfix) with ESMTPS id B1B541A923B
	for <gnats-bugs@gnats.NetBSD.org>; Fri, 11 Oct 2024 20:21:40 +0000 (UTC)
Message-Id: <20241011202139.815DF1A923D@mollari.NetBSD.org>
Date: Fri, 11 Oct 2024 20:21:39 +0000 (UTC)
From: furkanonder@protonmail.com
Reply-To: furkanonder@protonmail.com
To: gnats-bugs@NetBSD.org
Subject: Sign of Zero Issue with libc fma(): Positive Zero Returned Instead of Negative
X-Send-Pr-Version: www-1.0

>Number:         58742
>Category:       lib
>Synopsis:       Sign of Zero Issue with libc fma(): Positive Zero Returned Instead of Negative
>Confidential:   no
>Severity:       serious
>Priority:       medium
>Responsible:    riastradh
>State:          open
>Class:          sw-bug
>Submitter-Id:   net
>Arrival-Date:   Fri Oct 11 20:25:00 +0000 2024
>Last-Modified:  Sat Oct 12 19:50:01 +0000 2024
>Originator:     Furkan Onder
>Release:        NetBSD 10.0
>Organization:
>Environment:
NetBSD home.localhost 10.0 NetBSD 10.0 (GENERIC) #0: Thu Mar 28 08:33:33 UTC 2024  mkrepro@mkrepro.NetBSD.org:/usr/src/sys/arch/amd64/compile/GENERIC amd64
>Description:
There is an inconsistency in the behavior of the fma() function in the libc implementation, where it returns positive zero instead of the expected negative zero in certain floating-point arithmetic scenarios. The problem is that the Python test suite fails on NetBSD x86_64 on tests on the test_fma_zero_result().

https://github.com/python/cpython/blob/92760bd85b8f48b88df5b81100a757048979de83/Lib/test/test_math.py#L2724-L2776

>How-To-Repeat:
Example to Reproduce:
```
#include <math.h>
#include <stdio.h>

int main() {
  volatile double x = 0x1p-500, y = 0x1p-550, z = 0x1p-1000;
  double a, b, c, r;

  a = x - y;
  b = x + y;
  c = -z;
  r = fma(a, b, c);

  printf("fma(%+a, %+a, %+a) = %+g\n", a, b, c, r);
  return 0;
}
```

On Linux x86_64:
```
$ gcc fma-test.c -lm -o fma-test
$ ./fma-test
fma(+0x1.ffffffffffff8p-501, +0x1.0000000000004p-500, -0x1p-1000) = -0
```


On NetBSD x86_64:
```
$ gcc fma-test.c -lm -o fma-test
$ ./fma-test
fma(+0x1.ffffffffffff8p-501, +0x1.0000000000004p-500, -0x1p-1000) = +0
$
```

Program should print -0 instead of +0.
>Fix:

>Release-Note:

>Audit-Trail:

Responsible-Changed-From-To: lib-bug-people->riastradh
Responsible-Changed-By: riastradh@NetBSD.org
Responsible-Changed-When: Fri, 11 Oct 2024 20:28:29 +0000
Responsible-Changed-Why:
I'll take a look, thanks!


From: Taylor R Campbell <riastradh@NetBSD.org>
To: furkanonder@protonmail.com
Cc: gnats-bugs@NetBSD.org, netbsd-bugs@NetBSD.org
Subject: Re: lib/58742: Sign of Zero Issue with libc fma(): Positive Zero Returned Instead of Negative
Date: Sat, 12 Oct 2024 19:46:34 +0000

 I started adding FreeBSD's fma(3) tests, and adapted this one to use
 slightly less magic constants (so we can do the same for float and
 long double):

 	volatile double a = ldexp(1 - DBL_EPSILON, (DBL_MIN_EXP - 1)/2);
 	volatile double b = ldexp(1 + DBL_EPSILON, (DBL_MIN_EXP - 1)/2);
 	volatile double c = -DBL_MIN;

 That is,

 	a = 0x1.ffffffffffffep-512
 	b = 0x1.0000000000001p-511
 	c = -0x1p-1022

 for which a*b + c is a very small number below zero, -0x1p-1126, much
 closer to zero than to the smallest negative subnormal, -0x1p-1074.

 Stepping through fma(3) with gdb reveals that for these inputs, we end
 up at:

     268 	if (r.hi == 0.0) {
     269 		/*
     270 		 * When the addends cancel to 0, ensure that the result has
     271 		 * the correct sign.
     272 		 */
     273 		fesetround(oround);
     274 		{
     275 		volatile double vzs = zs; /* XXX gcc CSE bug workaround */
     276 		return (xy.hi + vzs + ldexp(xy.lo, spread));
     277 		}
     278 	}

 https://nxr.netbsd.org/xref/src/lib/libm/src/s_fma.c?r=1.7#268

 The local variables at this point are [hexadecimal notation added]:

 (gdb) info locals
 xs = 0.99999999999999978 [0x1.ffffffffffffep-1]
 ys = 0.50000000000000011 [0x1.0000000000001p-1]
 zs = -0.5
 adj = <optimized out>
 xy = {hi = 0.5, lo = -2.4651903288156619e-32 [-0x1p-105]}
 r = {hi = 0, lo = 0}
 oround = 0
 ex = -511
 ey = -510
 ez = -1021
 spread = -1021

 Summing xy.hi + zs yields 0 because xy.hi = zs, and scaling
 ldexp(xy.lo, spread) is rounded to zero, so it's zero plus zero plus
 zero.  The trouble is that the true sum (without rounding in ldexp)
 isn't actually zero, so we're not rounding correctly to plus or minus
 zero -- or, in non-default rounding modes like FE_DOWNWARD, away from
 zero.

 It's tempting, in this case, to try something like copysign(0, xy.lo)
 or, to handle rounding modes other than FE_TONEAREST, something like
 copysign(DBL_DENORM_MIN, xy.lo)/2, but that breaks some other cases.

>Unformatted:

Home	PR Database Search