math: portable FMA implementation incorrectly returns +0 when x*y ~ 0, x*y < 0 and z = 0

[shogo82148]

Go version

go version go1.24.3 linux/amd64

Output of go env in your module/workspace:

AR='ar'
CC='gcc'
CGO_CFLAGS='-O2 -g'
CGO_CPPFLAGS=''
CGO_CXXFLAGS='-O2 -g'
CGO_ENABLED='0'
CGO_FFLAGS='-O2 -g'
CGO_LDFLAGS='-O2 -g'
CXX='g++'
GCCGO='gccgo'
GO111MODULE=''
GOAMD64='v1'
GOARCH='amd64'
GOAUTH='netrc'
GOBIN=''
GOCACHE='/home/ec2-user/.cache/go-build'
GOCACHEPROG=''
GODEBUG=''
GOENV='/home/ec2-user/.config/go/env'
GOEXE=''
GOEXPERIMENT=''
GOFIPS140='off'
GOFLAGS=''
GOGCCFLAGS='-fPIC -m64 -fno-caret-diagnostics -Qunused-arguments -Wl,--no-gc-sections -fmessage-length=0 -ffile-prefix-map=/tmp/go-build2696373625=/tmp/go-build -gno-record-gcc-switches'
GOHOSTARCH='amd64'
GOHOSTOS='linux'
GOINSECURE=''
GOMOD='/home/ec2-user/fma/go.mod'
GOMODCACHE='/home/ec2-user/go/pkg/mod'
GONOPROXY=''
GONOSUMDB=''
GOOS='linux'
GOPATH='/home/ec2-user/go'
GOPRIVATE=''
GOPROXY='https://proxy.golang.org,direct'
GOROOT='/usr/local/go'
GOSUMDB='sum.golang.org'
GOTELEMETRY='local'
GOTELEMETRYDIR='/home/ec2-user/.config/go/telemetry'
GOTMPDIR=''
GOTOOLCHAIN='auto'
GOTOOLDIR='/usr/local/go/pkg/tool/linux_amd64'
GOVCS=''
GOVERSION='go1.24.3'
GOWORK=''
PKG_CONFIG='pkg-config'

What did you do?

https://go.dev/play/p/DgFhrLf1CuF

package main

import (
	"fmt"
	"math"
)

var portableFMA = math.FMA

func main() {
	fmt.Println(math.FMA(0x1p-1022, -0x1p-1022, 0))
	fmt.Println(portableFMA(0x1p-1022, -0x1p-1022, 0))
}

What did you see happen?

-0
0

What did you expect to see?

-0
-0

The exact result of the calculation is 0x1p-1022 * (-0x1p-1022) + 0 = -0x1p-2044. Since -0x1p-2044 cannot be represented in float64, it underflows and the result becomes -0.

[gabyhelp]

Related Issues

• math: portable FMA implementation incorrectly returns -0 in some situations #61130 (closed)
• cmd/compile: possibly incorrect fused-multiply-add bypassing floating-point to integer conversions #54533 (closed)
• x/arch: incorrect disassembly for x86 FMA instruction #41043
• cmd/compile: it is not possible to prevent FMA with complex values #36971
• math: FMA is slower than non-FMA calculation #36196 (closed)
• cmd/compile: miscompilation of simple arithmetic program targeting darwin/arm64 #62723 (closed)

Related Code Changes

• math: fix portable FMA when x*y < 0 and x*y == -z
• cmd/compile: fix FMA negative commutativity of riscv64
• cmd/compile: introduce generic ssa intrinsic for fused-multiply-add
• math: add guaranteed-precision FMA implementation

_{(Emoji vote if this was helpful or unhelpful; more detailed feedback welcome in this discussion.)}

[gopherbot]

Change https://go.dev/cl/673856 mentions this issue: math: fix portable FMA implementation when x*y ~ 0, x*y < 0 and z = 0

[randall77]

I think the fix you proposed is only correct if x*y is -0 because of rounding a very small result to -0. If x*y is -0 for another reason (-0 * 1, say), then the fix is wrong - we should get +0 in that case.

[shogo82148]

If x*y is -0 for another reason (-0 * 1, say), then the fix is wrong - we should get +0 in that case.

That case will be handled in these lines:

go/src/math/fma.go

Lines 98 to 101 in 6885cc0

	// Inf or NaN or zero involved. At most one rounding will occur.
	if x == 0.0 || y == 0.0 || bx&uvinf == uvinf || by&uvinf == uvinf {
	return x*y + z
	}

Therefore, x * y will result in -0 in the following code only if the result becomes very small and causes an underflow.

go/src/math/fma.go

Lines 102 to 108 in 6885cc0

	// Handle z == 0.0 separately.
	// Adding zero usually does not change the original value.
	// However, there is an exception with negative zero. (e.g. (-0) + (+0) = (+0))
	// This applies when x * y is negative and underflows.
	if z == 0.0 {
	return x * y
	}

[randall77]

Ah, ok. I guess the only way to get -0 from x*y without rounding is to have x or y be ==0.

[adonovan]

There's arguably a second bug here: the compiler should treat f := math.Foo; f(x) the same as math.Foo(x), even if calls to math.Foo are handled intrinsically. That would mean that taking the address of an intrinsic returns an eta abstraction.