[musl] Bug in fma(): should return negative zero

[musl] Bug in fma(): should return negative zero

[Victor Stinner]

Hi,

test_math.test_fma_zero_result() of Python 3.13 fails when Python is
built with musl:
https://github.com/python/cpython/issues/131032

The problem is that fma() returns +0.0 instead of -0.0 on some tests.

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

int main(int argc, char **argv)
{
    double tiny = 1e-300;
    double zero = 0.0;
    double result = fma(tiny, -tiny, zero);
    printf("fma(%+g, %+g, %+g) = %+g\n", tiny, -tiny, zero, result);
    return 0;
}
---

Output on Alpine with musl-1.2.5-r9:
---
$ gcc x.c -Og -o x -lm && ./x
fma(+1e-300, -1e-300, +0) = +0
---

Expected output:
---
fma(+1e-300, -1e-300, +0) = -0
---

Last year, I reported a similar bug to FreeBSD which was fixed:
https://bugs.freebsd.org/bugzilla/show_bug.cgi?id=277783

Link to Python test_fma_zero_result():
https://github.com/python/cpython/blob/425e0af74fb882b95f81923123bd4afad0cda157/Lib/test/test_math.py#L2762-L2816

Tell me if you need more details.

Victor
-- 
Night gathers, and now my watch begins. It shall not end until my death.

Re: [musl] Bug in fma(): should return negative zero

[Rich Felker]

On Tue, Mar 11, 2025 at 12:30:58PM +0100, Victor Stinner wrote:
> Hi,
> 
> test_math.test_fma_zero_result() of Python 3.13 fails when Python is
> built with musl:
> https://github.com/python/cpython/issues/131032
> 
> The problem is that fma() returns +0.0 instead of -0.0 on some tests.
> 
> Example:
> ---
> #include <math.h>
> #include <stdio.h>
> 
> int main(int argc, char **argv)
> {
>     double tiny = 1e-300;
>     double zero = 0.0;
>     double result = fma(tiny, -tiny, zero);
>     printf("fma(%+g, %+g, %+g) = %+g\n", tiny, -tiny, zero, result);
>     return 0;
> }
> ---
> 
> Output on Alpine with musl-1.2.5-r9:
> ---
> $ gcc x.c -Og -o x -lm && ./x
> fma(+1e-300, -1e-300, +0) = +0
> ---
> 
> Expected output:
> ---
> fma(+1e-300, -1e-300, +0) = -0
> ---
> 
> Last year, I reported a similar bug to FreeBSD which was fixed:
> https://bugs.freebsd.org/bugzilla/show_bug.cgi?id=277783
> 
> Link to Python test_fma_zero_result():
> https://github.com/python/cpython/blob/425e0af74fb882b95f81923123bd4afad0cda157/Lib/test/test_math.py#L2762-L2816
> 
> Tell me if you need more details.

What archs have you tested it on? It doesn't happen for me and the
general C code looks like it handles this case correctly, so I think
there must be broken asm for one or more archs which doesn't honor
this part of the fma contract right (due to broken cpu insn that
doesn't match IEEE semantics).

Rich

Re: [musl] Bug in fma(): should return negative zero

[Victor Stinner]

Hi,

On Tue, Mar 11, 2025 at 2:25 PM Rich Felker <dalias@libc.org> wrote:
> What archs have you tested it on?

I tested on x86-64.

> It doesn't happen for me and the general C code looks like it handles
> this case correctly,

Oh. I can reproduce the issue with my C example on:

* Fedora 41: musl 1.2.5-1, musl-gcc version "gcc (GCC) 14.2.1"
* Alpine Linux v3.21: musl-1.2.5-r9, gcc version 14.2.0.

Example on Fedora:
---
$ musl-gcc x.c -o x && ./x
fma(+1e-300, -1e-300, +0) = +0
---

The result is wrong (+0), it must be "-0": fma(x, y, z) = (+1e-300) *
(-1e-300) + (+0) = -0.

Using glibc, I get the correct result (-0):
---
$ gcc x.c -o x -lm && ./x
fma(+1e-300, -1e-300, +0) = -0
---

Note: when I wrote a reproducer for FreeBSD, I add to trick the clang
compiler to prevent it to call fma() at build time. But I didn't need
to do it for this musl reproducer (using gcc). I checked using gdb
that main() actually calls fma() at runtime:
--------
$ gdb ./x
(gdb) disassemble main
...
   0x0000000000401192 <+76>:    call   0x401030 <fma@plt>
...
--------

Victor
-- 
Night gathers, and now my watch begins. It shall not end until my death.

Re: [musl] Bug in fma(): should return negative zero

[Rich Felker]

On Tue, Mar 11, 2025 at 05:19:26PM +0100, Victor Stinner wrote:
> Hi,
> 
> On Tue, Mar 11, 2025 at 2:25 PM Rich Felker <dalias@libc.org> wrote:
> > What archs have you tested it on?
> 
> I tested on x86-64.
> 
> > It doesn't happen for me and the general C code looks like it handles
> > this case correctly,
> 
> Oh. I can reproduce the issue with my C example on:
> 
> * Fedora 41: musl 1.2.5-1, musl-gcc version "gcc (GCC) 14.2.1"
> * Alpine Linux v3.21: musl-1.2.5-r9, gcc version 14.2.0.
> 
> Example on Fedora:
> ---
> $ musl-gcc x.c -o x && ./x
> fma(+1e-300, -1e-300, +0) = +0
> ---
> 
> The result is wrong (+0), it must be "-0": fma(x, y, z) = (+1e-300) *
> (-1e-300) + (+0) = -0.

Indeed, the correct result is -0 and +0 is wrong. But there is no bug
in musl. The bug is in the version of gcc used to compile musl. Using
gcc 13.2.1 on Alpine 3.20 and building musl myself, fma correctly
returns -0. But the musl binary shipped with Alpine 3.20 has been
miscompiled so that the "return x*y;" on line 56 incorrectly just
branches to line 53's "return x*y+z;".

I'm not sure what the buggy version of gcc is, but I would guess it's
a new regression in gcc 14. So I guess we need to track this down and
report it. My guess is that someone wrongly thought they could infer
from the branches on the bit representation of the float that it's
zero, and thereby "of course adding zero is a no-op" (no! no! NO!) so
they folded the two code paths to the same thing...

Rich

Re: [musl] Bug in fma(): should return negative zero

[Rich Felker]

On Tue, Mar 11, 2025 at 01:09:23PM -0400, Rich Felker wrote:
> On Tue, Mar 11, 2025 at 05:19:26PM +0100, Victor Stinner wrote:
> > Hi,
> > 
> > On Tue, Mar 11, 2025 at 2:25 PM Rich Felker <dalias@libc.org> wrote:
> > > What archs have you tested it on?
> > 
> > I tested on x86-64.
> > 
> > > It doesn't happen for me and the general C code looks like it handles
> > > this case correctly,
> > 
> > Oh. I can reproduce the issue with my C example on:
> > 
> > * Fedora 41: musl 1.2.5-1, musl-gcc version "gcc (GCC) 14.2.1"
> > * Alpine Linux v3.21: musl-1.2.5-r9, gcc version 14.2.0.
> > 
> > Example on Fedora:
> > ---
> > $ musl-gcc x.c -o x && ./x
> > fma(+1e-300, -1e-300, +0) = +0
> > ---
> > 
> > The result is wrong (+0), it must be "-0": fma(x, y, z) = (+1e-300) *
> > (-1e-300) + (+0) = -0.
> 
> Indeed, the correct result is -0 and +0 is wrong. But there is no bug
> in musl. The bug is in the version of gcc used to compile musl. Using
> gcc 13.2.1 on Alpine 3.20 and building musl myself, fma correctly
> returns -0. But the musl binary shipped with Alpine 3.20 has been
> miscompiled so that the "return x*y;" on line 56 incorrectly just
> branches to line 53's "return x*y+z;".
> 
> I'm not sure what the buggy version of gcc is, but I would guess it's
> a new regression in gcc 14. So I guess we need to track this down and
> report it. My guess is that someone wrongly thought they could infer
> from the branches on the bit representation of the float that it's
> zero, and thereby "of course adding zero is a no-op" (no! no! NO!) so
> they folded the two code paths to the same thing...

Ah, no, I should have checked. This is a bug in musl that was fixed by
nsz last year in git, but the fix hasn't made it into any release.

Rich

Re: [musl] Bug in fma(): should return negative zero

[Khem Raj]

On Tue, Mar 11, 2025 at 10:35 AM Rich Felker <dalias@libc.org> wrote:
>
> On Tue, Mar 11, 2025 at 01:09:23PM -0400, Rich Felker wrote:
> > On Tue, Mar 11, 2025 at 05:19:26PM +0100, Victor Stinner wrote:
> > > Hi,
> > >
> > > On Tue, Mar 11, 2025 at 2:25 PM Rich Felker <dalias@libc.org> wrote:
> > > > What archs have you tested it on?
> > >
> > > I tested on x86-64.
> > >
> > > > It doesn't happen for me and the general C code looks like it handles
> > > > this case correctly,
> > >
> > > Oh. I can reproduce the issue with my C example on:
> > >
> > > * Fedora 41: musl 1.2.5-1, musl-gcc version "gcc (GCC) 14.2.1"
> > > * Alpine Linux v3.21: musl-1.2.5-r9, gcc version 14.2.0.
> > >
> > > Example on Fedora:
> > > ---
> > > $ musl-gcc x.c -o x && ./x
> > > fma(+1e-300, -1e-300, +0) = +0
> > > ---
> > >
> > > The result is wrong (+0), it must be "-0": fma(x, y, z) = (+1e-300) *
> > > (-1e-300) + (+0) = -0.
> >
> > Indeed, the correct result is -0 and +0 is wrong. But there is no bug
> > in musl. The bug is in the version of gcc used to compile musl. Using
> > gcc 13.2.1 on Alpine 3.20 and building musl myself, fma correctly
> > returns -0. But the musl binary shipped with Alpine 3.20 has been
> > miscompiled so that the "return x*y;" on line 56 incorrectly just
> > branches to line 53's "return x*y+z;".
> >
> > I'm not sure what the buggy version of gcc is, but I would guess it's
> > a new regression in gcc 14. So I guess we need to track this down and
> > report it. My guess is that someone wrongly thought they could infer
> > from the branches on the bit representation of the float that it's
> > zero, and thereby "of course adding zero is a no-op" (no! no! NO!) so
> > they folded the two code paths to the same thing...
>
> Ah, no, I should have checked. This is a bug in musl that was fixed by
> nsz last year in git, but the fix hasn't made it into any release.
>

In yocto we use the tip of the trunk and I was not able to reproduce as well.
This is the commit which is helping I think.

https://git.musl-libc.org/cgit/musl/commit/?id=9683bd62414604d3bd56cf6bd7be8f54aa31e7d3

> Rich

Re: [musl] Bug in fma(): should return negative zero

[Rich Felker]

On Tue, Mar 11, 2025 at 10:59:52AM -0700, Khem Raj wrote:
> On Tue, Mar 11, 2025 at 10:35 AM Rich Felker <dalias@libc.org> wrote:
> >
> > On Tue, Mar 11, 2025 at 01:09:23PM -0400, Rich Felker wrote:
> > > On Tue, Mar 11, 2025 at 05:19:26PM +0100, Victor Stinner wrote:
> > > > Hi,
> > > >
> > > > On Tue, Mar 11, 2025 at 2:25 PM Rich Felker <dalias@libc.org> wrote:
> > > > > What archs have you tested it on?
> > > >
> > > > I tested on x86-64.
> > > >
> > > > > It doesn't happen for me and the general C code looks like it handles
> > > > > this case correctly,
> > > >
> > > > Oh. I can reproduce the issue with my C example on:
> > > >
> > > > * Fedora 41: musl 1.2.5-1, musl-gcc version "gcc (GCC) 14.2.1"
> > > > * Alpine Linux v3.21: musl-1.2.5-r9, gcc version 14.2.0.
> > > >
> > > > Example on Fedora:
> > > > ---
> > > > $ musl-gcc x.c -o x && ./x
> > > > fma(+1e-300, -1e-300, +0) = +0
> > > > ---
> > > >
> > > > The result is wrong (+0), it must be "-0": fma(x, y, z) = (+1e-300) *
> > > > (-1e-300) + (+0) = -0.
> > >
> > > Indeed, the correct result is -0 and +0 is wrong. But there is no bug
> > > in musl. The bug is in the version of gcc used to compile musl. Using
> > > gcc 13.2.1 on Alpine 3.20 and building musl myself, fma correctly
> > > returns -0. But the musl binary shipped with Alpine 3.20 has been
> > > miscompiled so that the "return x*y;" on line 56 incorrectly just
> > > branches to line 53's "return x*y+z;".
> > >
> > > I'm not sure what the buggy version of gcc is, but I would guess it's
> > > a new regression in gcc 14. So I guess we need to track this down and
> > > report it. My guess is that someone wrongly thought they could infer
> > > from the branches on the bit representation of the float that it's
> > > zero, and thereby "of course adding zero is a no-op" (no! no! NO!) so
> > > they folded the two code paths to the same thing...
> >
> > Ah, no, I should have checked. This is a bug in musl that was fixed by
> > nsz last year in git, but the fix hasn't made it into any release.
> >
> 
> In yocto we use the tip of the trunk and I was not able to reproduce as well.
> This is the commit which is helping I think.
> 
> https://git.musl-libc.org/cgit/musl/commit/?id=9683bd62414604d3bd56cf6bd7be8f54aa31e7d3

Yep, that's the one and it fixes exactly this bug.

Rich