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 <[email protected]>
Run-TryBot: Robert Griesemer <[email protected]>
TryBot-Result: Gobot Gobot <[email protected]>

Author: mpx

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++ {

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.

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)
+}