math: add RoundToEven function

Rounding ties to even is statistically useful for some applications.
This implementation completes IEEE float64 rounding mode support (in
addition to Round, Ceil, Floor, Trunc).

This function avoids subtle faults found in ad-hoc implementations, and
is simple enough to be inlined by the compiler.

Fixes #21748

Change-Id: I09415df2e42435f9e7dabe3bdc0148e9b9ebd609
Reviewed-on: https://go-review.googlesource.com/61211
Reviewed-by: Robert Griesemer <gri@golang.org>
Run-TryBot: Robert Griesemer <gri@golang.org>
TryBot-Result: Gobot Gobot <gobot@golang.org>

3 files changed, 81 insertions, 8 deletions

diff --git a/src/math/all_test.go b/src/math/all_test.go
index d0630aef44..7598d88570 100644
--- a/src/math/all_test.go
+++ b/src/math/all_test.go
@@ -1774,9 +1774,26 @@ var vfroundSC = [][2]float64{
 	{0.5, 1},
 	{0.5000000000000001, 1}, // 0.5+epsilon
 	{-1.5, -2},
+	{-2.5, -3},
 	{NaN(), NaN()},
 	{Inf(1), Inf(1)},
 	{2251799813685249.5, 2251799813685250}, // 1 bit fraction
+	{2251799813685250.5, 2251799813685251},
+	{4503599627370495.5, 4503599627370496}, // 1 bit fraction, rounding to 0 bit fraction
+	{4503599627370497, 4503599627370497},   // large integer
+}
+var vfroundEvenSC = [][2]float64{
+	{0, 0},
+	{1.390671161567e-309, 0}, // denormal
+	{0.49999999999999994, 0}, // 0.5-epsilon
+	{0.5, 0},
+	{0.5000000000000001, 1}, // 0.5+epsilon
+	{-1.5, -2},
+	{-2.5, -2},
+	{NaN(), NaN()},
+	{Inf(1), Inf(1)},
+	{2251799813685249.5, 2251799813685250}, // 1 bit fraction
+	{2251799813685250.5, 2251799813685250},
 	{4503599627370495.5, 4503599627370496}, // 1 bit fraction, rounding to 0 bit fraction
 	{4503599627370497, 4503599627370497},   // large integer
 }
@@ -2752,6 +2769,19 @@ func TestRound(t *testing.T) {
 	}
 }
 
+func TestRoundToEven(t *testing.T) {
+	for i := 0; i < len(vf); i++ {
+		if f := RoundToEven(vf[i]); !alike(round[i], f) {
+			t.Errorf("RoundToEven(%g) = %g, want %g", vf[i], f, round[i])
+		}
+	}
+	for i := 0; i < len(vfroundEvenSC); i++ {
+		if f := RoundToEven(vfroundEvenSC[i][0]); !alike(vfroundEvenSC[i][1], f) {
+			t.Errorf("RoundToEven(%g) = %g, want %g", vfroundEvenSC[i][0], f, vfroundEvenSC[i][1])
+		}
+	}
+}
+
 func TestSignbit(t *testing.T) {
 	for i := 0; i < len(vf); i++ {
 		if f := Signbit(vf[i]); signbit[i] != f {
@@ -3413,6 +3443,14 @@ func BenchmarkRound(b *testing.B) {
 	GlobalF = x
 }
 
+func BenchmarkRoundToEven(b *testing.B) {
+	x := 0.0
+	for i := 0; i < b.N; i++ {
+		x = RoundToEven(roundNeg)
+	}
+	GlobalF = x
+}
+
 func BenchmarkRemainder(b *testing.B) {
 	x := 0.0
 	for i := 0; i < b.N; i++ {
diff --git a/src/math/bits.go b/src/math/bits.go
index d85ee9cb13..77bcdbe1ce 100644
--- a/src/math/bits.go
+++ b/src/math/bits.go
@@ -8,9 +8,12 @@ const (
 	uvnan    = 0x7FF8000000000001
 	uvinf    = 0x7FF0000000000000
 	uvneginf = 0xFFF0000000000000
+	uvone    = 0x3FF0000000000000
 	mask     = 0x7FF
 	shift    = 64 - 11 - 1
 	bias     = 1023
+	signMask = 1 << 63
+	fracMask = 1<<shift - 1
 )
 
 // Inf returns positive infinity if sign >= 0, negative infinity if sign < 0.
diff --git a/src/math/floor.go b/src/math/floor.go
index d03c4bcdad..18e89ef89f 100644
--- a/src/math/floor.go
+++ b/src/math/floor.go
@@ -71,29 +71,61 @@ func Round(x float64) float64 {
 	//   }
 	//   return t
 	// }
-	const (
-		signMask = 1 << 63
-		fracMask = 1<<shift - 1
-		half     = 1 << (shift - 1)
-		one      = bias << shift
-	)
-
 	bits := Float64bits(x)
 	e := uint(bits>>shift) & mask
 	if e < bias {
 		// Round abs(x) < 1 including denormals.
 		bits &= signMask // +-0
 		if e == bias-1 {
-			bits |= one // +-1
+			bits |= uvone // +-1
 		}
 	} else if e < bias+shift {
 		// Round any abs(x) >= 1 containing a fractional component [0,1).
 		//
 		// Numbers with larger exponents are returned unchanged since they
 		// must be either an integer, infinity, or NaN.
+		const half = 1 << (shift - 1)
 		e -= bias
 		bits += half >> e
 		bits &^= fracMask >> e
 	}
 	return Float64frombits(bits)
 }
+
+// RoundToEven returns the nearest integer, rounding ties to even.
+//
+// Special cases are:
+//	RoundToEven(±0) = ±0
+//	RoundToEven(±Inf) = ±Inf
+//	RoundToEven(NaN) = NaN
+func RoundToEven(x float64) float64 {
+	// RoundToEven is a faster implementation of:
+	//
+	// func RoundToEven(x float64) float64 {
+	//   t := math.Trunc(x)
+	//   odd := math.Remainder(t, 2) != 0
+	//   if d := math.Abs(x - t); d > 0.5 || (d == 0.5 && odd) {
+	//     return t + math.Copysign(1, x)
+	//   }
+	//   return t
+	// }
+	bits := Float64bits(x)
+	e := uint(bits>>shift) & mask
+	if e >= bias {
+		// Round abs(x) >= 1.
+		// - Large numbers without fractional components, infinity, and NaN are unchanged.
+		// - Add 0.499.. or 0.5 before truncating depending on whether the truncated
+		//   number is even or odd (respectively).
+		const halfMinusULP = (1 << (shift - 1)) - 1
+		e -= bias
+		bits += (halfMinusULP + (bits>>(shift-e))&1) >> e
+		bits &^= fracMask >> e
+	} else if e == bias-1 && bits&fracMask != 0 {
+		// Round 0.5 < abs(x) < 1.
+		bits = bits&signMask | uvone // +-1
+	} else {
+		// Round abs(x) <= 0.5 including denormals.
+		bits &= signMask // +-0
+	}
+	return Float64frombits(bits)
+}