cmd/compile: unstable arithmetic

[changkun]

What version of Go are you using (go version)?

$ go version
go version go1.18 darwin/arm64

Does this issue reproduce with the latest release?

Yes

What operating system and processor architecture are you using (go env)?

go env Output

$ go env
GO111MODULE="auto"
GOARCH="arm64"
GOBIN=""
GOCACHE="/Users/changkun/Library/Caches/go-build"
GOENV="/Users/changkun/Library/Application Support/go/env"
GOEXE=""
GOEXPERIMENT=""
GOFLAGS=""
GOHOSTARCH="arm64"
GOHOSTOS="darwin"
GOINSECURE=""
GOMODCACHE="/Users/changkun/go/pkg/mod"
GONOPROXY=""
GONOSUMDB=""
GOOS="darwin"
GOPATH="/Users/changkun/go"
GOPRIVATE=""
GOPROXY="https://proxy.golang.org,direct"
GOROOT="/Users/changkun/goes/go"
GOSUMDB="sum.golang.org"
GOTMPDIR=""
GOTOOLDIR="/Users/changkun/goes/go/pkg/tool/darwin_arm64"
GOVCS=""
GOVERSION="go1.18"
GCCGO="gccgo"
AR="ar"
CC="clang"
CXX="clang++"
CGO_ENABLED="1"
GOMOD="/Users/changkun/dev/poly.red/polyred/go.mod"
GOWORK=""
CGO_CFLAGS="-g -O2"
CGO_CPPFLAGS=""
CGO_CXXFLAGS="-g -O2"
CGO_FFLAGS="-g -O2"
CGO_LDFLAGS="-g -O2"
PKG_CONFIG="pkg-config"
GOGCCFLAGS="-fPIC -arch arm64 -pthread -fno-caret-diagnostics -Qunused-arguments -fmessage-length=0 -fdebug-prefix-map=/var/folders/g2/6fmr1qzx0ns3shq74zrp6bd40000gn/T/go-build1155872957=/tmp/go-build -gno-record-gcc-switches -fno-common"

What did you do?

https://go.dev/play/p/5zmfDpA4ylc

package main

import (
	"math"
)

type Vec3 struct {
	X, Y, Z float32
}

func (v Vec3) Unit() Vec3 {
	n := 1.0 / float32(math.Sqrt(float64(v.X*v.X+v.Y*v.Y+v.Z*v.Z)))
	return Vec3{v.X * n, v.Y * n, v.Z * n}
}

func (v Vec3) Cross(u Vec3) Vec3 {
	x := v.Y*u.Z - v.Z*u.Y
	y := v.Z*u.X - v.X*u.Z
	z := v.X*u.Y - v.Y*u.X
	return Vec3{x, y, z}
}

func compute1(a, b, c Vec3) float32 {
	l := Vec3{b.X - a.X, b.Y - a.Y, b.Z - a.Z}.Unit()
	x := l.Cross(c).Unit()
	return -x.X*a.X - x.Y*a.Y - x.Z*a.Z // v1
}

func compute2(a, b, c Vec3) float32 {
	l := Vec3{b.X - a.X, b.Y - a.Y, b.Z - a.Z}.Unit()
	x := l.Cross(c).Unit()
	return -(x.X*a.X + x.Y*a.Y + x.Z*a.Z) // v2
}

func compute3(x, a Vec3) float32 {
	return -x.X*a.X - x.Y*a.Y - x.Z*a.Z // v3
}

func compute4(x, a Vec3) float32 {
	return -(x.X*a.X + x.Y*a.Y + x.Z*a.Z) // v4
}

func main() {
	a, b, c := Vec3{-550, 194, 734}, Vec3{-1000, 0, 0}, Vec3{0, 1, 1}
	v1 := compute1(a, b, c)
	v2 := compute2(a, b, c)

	l := Vec3{b.X - a.X, b.Y - a.Y, b.Z - a.Z}.Unit()
	x := l.Cross(c).Unit()
	v3 := compute3(x, a)
	v4 := compute4(x, a)

	println(v1 == v2)
	println(v1 == v3)
	println(v1 == v4)
	println(v2 == v3)
	println(v2 == v4)
	println(v3 == v4)
}

What did you expect to see?

true
true
true
true
true
true

This result is what go.dev/play produces

What did you see instead?

false
false
false
true
true
true

This is the result produced on a M1 mac.

[cherrymui]

This is expected, at least according to the current spec. See https://go.dev/ref/spec#Arithmetic_operators , look for "fused multiply and add". We do use FMA on ARM64 (and some other architectures).

[changkun]

This is expected, at least according to the current spec. See https://go.dev/ref/spec#Arithmetic_operators , look for "fused multiply and add". We do use FMA on ARM64 (and some other architectures).

Thanks! That's indeed good to know. But if I understand the behavior correctly, and it is all about FMA. Shouldn't we at least have consistent results for each of the return statements on a specific architecture? Meaning for all compute functions:

func compute1(a, b, c Vec3) float32 {
	...
	return -x.X*a.X - x.Y*a.Y - x.Z*a.Z // v1
}

func compute2(a, b, c Vec3) float32 {
	...
	return -(x.X*a.X + x.Y*a.Y + x.Z*a.Z) // v2
}

func compute3(x, a Vec3) float32 {
	return -x.X*a.X - x.Y*a.Y - x.Z*a.Z // v3
}

func compute4(x, a Vec3) float32 {
	return -(x.X*a.X + x.Y*a.Y + x.Z*a.Z) // v4
}

they at least all produce the same results (does not necessarily require linux/amd64 == darwin/arm64)?

[cherrymui]

The + and - are the difference. They are not the same order of operation, e.g. compute1 has multiply, negate, multiply, subtract, ..., compute2 has multiply, add, multiply, ..., negate. This would also cause different operations being fused in different cases.

If it really matters, you can try adding explicit rounding, e.g. -float32(x.X*a.X) - ...

[changkun]

Isn't compute1 and compute3 exactly the same in its arthmetic order? They are still result in a false comparison.

I can't certainly answer regarding the influences yet without enough in-depth practice. But as far as I could imagine, it might influence the numeric stability of different solvers and convey hidden issues, which eventually cause unexpected results on different platforms (e.g., Laplacian and Poisson solvers that compute differently optimized geometric structures where one is a good geometry and the other is a faulty geometry).

[ianlancetaylor]

The compiler is permitted but not required to use FMA for any given expression.

People who need precise control over floating-point operations can and should use conversions to float32 as @cherrymui suggested and/or call math.FMA.

[changkun]

Thanks for the suggestions! I will close this for now and need further develop some expertise here and see how to tacle this precisely.