runtime: fix softfloat64 add/sub for normal operands cancelling to subnormal

When fadd64/fsub64 cancel two near-equal, opposite-sign normal operands
into a subnormal result, fpack64 is handed a mantissa far below
1<<mantbits64 while exp is still a normal-range exponent. fpack64 saved
mant0/exp0/trunc0 before its normalization loop, so the denormal branch
restored an un-normalized mantissa and right-shifted it to align to the
subnormal exponent (the wrong direction here), returning the
un-normalized cancellation mantissa at the wrong scale instead of the
correctly rounded subnormal.

Save mant0/exp0/trunc0 after the normalization loop so the denormal path
restores a normalized mantissa and aligns correctly. This is a no-op for
callers that already pass a normalized mantissa (fmul64, fdiv64,
conversions). fpack32 is updated identically for parity, though its
denormal branch is not reachable with an un-normalized mantissa through
any current caller.

This only manifests on softfloat targets (e.g. GOMIPS=softfloat), which
is why it has gone unnoticed on hardware-float platforms. A randomized
differential check against hardware found the previous code wrong on
>50% of normal pairs that cancel into a subnormal; with this change
those cases match hardware.

Fixes #79964

Change-Id: I29f7ec79905ddb98af4f997a69b0a2f0574c3c9d
GitHub-Last-Rev: 42a73cd890f798e432c1aa5a279c758cb465e323
GitHub-Pull-Request: golang/go#79965
Reviewed-on: https://go-review.googlesource.com/c/go/+/789861
Reviewed-by: Sayer Turner <sayerturner65@gmail.com>
Reviewed-by: Keith Randall <khr@google.com>
Reviewed-by: Dmitri Shuralyov <dmitshur@google.com>
Reviewed-by: Keith Randall <khr@golang.org>
Auto-Submit: Keith Randall <khr@golang.org>
LUCI-TryBot-Result: golang-scoped@luci-project-accounts.iam.gserviceaccount.com <golang-scoped@luci-project-accounts.iam.gserviceaccount.com>
diff --git a/src/runtime/softfloat64.go b/src/runtime/softfloat64.go
index 42ef009..c230c64 100644
--- a/src/runtime/softfloat64.go
+++ b/src/runtime/softfloat64.go
@@ -91,7 +91,6 @@
 }
 
 func fpack64(sign, mant uint64, exp int, trunc uint64) uint64 {
-	mant0, exp0, trunc0 := mant, exp, trunc
 	if mant == 0 {
 		return sign
 	}
@@ -99,6 +98,11 @@
 		mant <<= 1
 		exp--
 	}
+	// Save the normalized mantissa; the denormal path below restores it and
+	// re-aligns to the subnormal exponent. Saving before this loop (as the
+	// code originally did) left a heavily-cancelled add/sub mantissa
+	// un-normalized, which that path then shifted the wrong way. See #79964.
+	mant0, exp0, trunc0 := mant, exp, trunc
 	for mant >= 4<<mantbits64 {
 		trunc |= mant & 1
 		mant >>= 1
@@ -142,7 +146,6 @@
 }
 
 func fpack32(sign, mant uint32, exp int, trunc uint32) uint32 {
-	mant0, exp0, trunc0 := mant, exp, trunc
 	if mant == 0 {
 		return sign
 	}
@@ -150,6 +153,8 @@
 		mant <<= 1
 		exp--
 	}
+	// See fpack64: save the normalized mantissa for the denormal path below.
+	mant0, exp0, trunc0 := mant, exp, trunc
 	for mant >= 4<<mantbits32 {
 		trunc |= mant & 1
 		mant >>= 1
diff --git a/src/runtime/softfloat64_test.go b/src/runtime/softfloat64_test.go
index 3f53e8b..1a0a502 100644
--- a/src/runtime/softfloat64_test.go
+++ b/src/runtime/softfloat64_test.go
@@ -36,10 +36,12 @@
 		math.Inf(-1),
 		0.1,
 		1.5,
-		1.9999999999999998,     // all 1s mantissa
-		1.3333333333333333,     // 1.010101010101...
-		1.1428571428571428,     // 1.001001001001...
-		1.112536929253601e-308, // first normal
+		1.9999999999999998,      // all 1s mantissa
+		1.3333333333333333,      // 1.010101010101...
+		1.1428571428571428,      // 1.001001001001...
+		1.112536929253601e-308,  // first normal
+		-2.662107816930723e-301, // #79964: two near-equal opposite-sign
+		2.662107858822336e-301,  // normals cancelling to a subnormal
 		2,
 		4,
 		8,