[Bug 277783] libc fma() doesn't not return the correct zero sign

https://bugs.freebsd.org/bugzilla/show_bug.cgi?id=277783

            Bug ID: 277783
           Summary: libc fma() doesn't not return the correct zero sign
           Product: Base System
           Version: 14.0-RELEASE
          Hardware: Any
                OS: Any
            Status: New
          Severity: Affects Only Me
          Priority: ---
         Component: kern
          Assignee: bugs@FreeBSD.org
          Reporter: victor.stinner@gmail.com
 Attachment #249266 text/plain
         mime type:

Created attachment 249266
  --> https://bugs.freebsd.org/bugzilla/attachment.cgi?id=249266&action=edit
Example reproducing the fma() zero sign issue

Hi,

I exposed the libc fma() function in Python:

*
https://github.com/python/cpython/commit/8e3c953b3acd7c656f66696806c9ae917e816492
* https://github.com/python/cpython/pull/116667
* https://github.com/python/cpython/issues/73468

The problem is that the Python test suite fails on FreeBSD x86-64 on tests on
the zero sign: test_fma_zero_result()
https://github.com/python/cpython/blob/cd2ed917801b93fb46d1dcf19dd480e5146932d8/Lib/test/test_math.py#L2694-L2748

Example:
---
#include <math.h>
#include <string.h>
#include <stdio.h>

void set_double(double *x, const char *bytes)
{
    memcpy(x, bytes, sizeof(*x));
}

int main(int argc, char **argv)
{
    double x = 0; //float.fromhex('0x1p-500')
    double y = 0; //float.fromhex('0x1p-550')
    double z = 0; //float.fromhex('0x1p-1000')
#if 1
    // condition always true, but trick the compiler, otherwise it will not
    // call fma() at runtime, but just compute the result at build time.
    if (argc) {
        set_double(&x, "\x00\x00\x00\x00\x00\x00\xb0\x20");
    }
    else {
        x = 123;
    }
#else
    set_double(&x, "\x00\x00\x00\x00\x00\x00\xb0\x20");
#endif
    set_double(&y, "\x00\x00\x00\x00\x00\x00\x90\x1d");
    set_double(&z, "\x00\x00\x00\x00\x00\x00p\x01");

    printf("x %g\n", x);
    printf("y %g\n", y);
    printf("z %g\n", z);

    double a, b, c, r;
    a = x-y;
    b = x+y;
    c = -z;
    r = fma(a, b, c);
    printf("fma(x-y, x+y, -z):\n");
    printf("fma(%+g, %+g, %+g) = %+g\n", a, b, c, r);
    return 0;
}
---

Output on x86-64:
---
$ clang x.c -Og -o x -lm && ./x
x 3.05494e-151
y 2.71333e-166
z 9.33264e-302
fma(x-y, x+y, -z):
fma(+3.05494e-151, +3.05494e-151, -9.33264e-302) = +0
---

fma() result is "+0" which is wrong. For example, if you replace "#if 1" with
"#if 0", you get:
---
(...)
fma(+3.05494e-151, +3.05494e-151, -9.33264e-302) = -0
---

fma() result is now "-0" which is correct. This time, the result was computed
at build time.


Extract of the Python tests which fail on FreeBSD:
---
        # Corner case where rounding the multiplication would
        # give the wrong result.
        x = float.fromhex('0x1p-500')
        y = float.fromhex('0x1p-550')
        z = float.fromhex('0x1p-1000')
        self.assertIsNegativeZero(math.fma(x-y, x+y, -z))
        self.assertIsPositiveZero(math.fma(y-x, x+y, z))
        self.assertIsNegativeZero(math.fma(y-x, -(x+y), -z))
        self.assertIsPositiveZero(math.fma(x-y, -(x+y), z))
---

For example, on Linux x86-64, the example works as expected with GCC and clang:
---
vstinner@mona$ gcc x.c -o x -O2 -lm && ./x
x 3.05494e-151
y 2.71333e-166
z 9.33264e-302
fma(x-y, x+y, -z):
fma(+3.05494e-151, +3.05494e-151, -9.33264e-302) = -0


vstinner@mona$ clang x.c -o x -O2 -lm && ./x
x 3.05494e-151
y 2.71333e-166
z 9.33264e-302
fma(x-y, x+y, -z):
fma(+3.05494e-151, +3.05494e-151, -9.33264e-302) = -0
---

fma() returns -0 which is correct.


On the bug, I selected "kern" component. I expect that the bug is more in the
libc, but I cannot see "libc" choice.

-- 
You are receiving this mail because:
You are the assignee for the bug.