math: fix portable FMA when x*y < 0 and x*y == -z

When x*y == -z the portable implementation of FMA copied the sign
bit from x*y into the result. This meant that when x*y == -z and
x*y < 0 the result was -0 which is incorrect.

Fixes #61130.

Change-Id: Ib93a568b7bdb9031e2aedfa1bdfa9bddde90851d
Reviewed-on: https://go-review.googlesource.com/c/go/+/507376
Reviewed-by: Dmitri Shuralyov <dmitshur@google.com>
Reviewed-by: Keith Randall <khr@golang.org>
Run-TryBot: Michael Munday <mike.munday@lowrisc.org>
TryBot-Result: Gopher Robot <gobot@golang.org>
Reviewed-by: Joedian Reid <joedian@golang.org>

diff --git a/src/math/all_test.go b/src/math/all_test.go
index 886267bc1756c108fc87f1551eeda06b60574b08..96a398e9c6f3c214496ef48651d49c5d476a151f 100644 (file)
--- a/src/math/all_test.go
+++ b/src/math/all_test.go
@@ -2059,6 +2059,9 @@ var fmaC = []struct{ x, y, z, want float64 }{
 
        // Special
        {0, 0, 0, 0},
+       {Copysign(0, -1), 0, 0, 0},
+       {0, 0, Copysign(0, -1), 0},
+       {Copysign(0, -1), 0, Copysign(0, -1), Copysign(0, -1)},
        {-1.1754226043408471e-38, NaN(), Inf(0), NaN()},
        {0, 0, 2.22507385643494e-308, 2.22507385643494e-308},
        {-8.65697792e+09, NaN(), -7.516192799999999e+09, NaN()},
@@ -2077,6 +2080,10 @@ var fmaC = []struct{ x, y, z, want float64 }{
        {4.612811918325842e+18, 1.4901161193847641e-08, 2.6077032311277997e-08, 6.873625395187494e+10},
        {-9.094947033611148e-13, 4.450691014249257e-308, 2.086006742350485e-308, 2.086006742346437e-308},
        {-7.751454006381804e-05, 5.588653777189071e-308, -2.2207280111272877e-308, -2.2211612130544025e-308},
+
+       // Issue #61130
+       {-1, 1, 1, 0},
+       {1, 1, -1, 0},
 }
 
 var sqrt32 = []float32{

diff --git a/src/math/fma.go b/src/math/fma.go
index ca0bf99f215603579286bf250cf6ae31c51cb193..ba03fbe8a93b27a2834b3248cc1c1fbb0804177a 100644 (file)
--- a/src/math/fma.go
+++ b/src/math/fma.go
@@ -132,6 +132,11 @@ func FMA(x, y, z float64) float64 {
                ps, pe, pm1, pm2, zs, ze, zm1, zm2 = zs, ze, zm1, zm2, ps, pe, pm1, pm2
        }
 
+       // Special case: if p == -z the result is always +0 since neither operand is zero.
+       if ps != zs && pe == ze && pm1 == zm1 && pm2 == zm2 {
+               return 0
+       }
+
        // Align significands
        zm1, zm2 = shrcompress(zm1, zm2, uint(pe-ze))