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>

diff --git a/src/math/all_test.go b/src/math/all_test.go
index d0630aef4464750da14feffa5c4f1459869354e3..7598d88570d93be4e784164ec002f99556266c28 100644 (file)
--- 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 d85ee9cb137775ad39aaafc9c2e72b851449abdd..77bcdbe1ce649a662049c005e34c06dfbecd6e86 100644 (file)
--- 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 d03c4bcdadbc4b3dc3760ac70b1afba4ad018326..18e89ef89f3d6bfa4cfa72a4352fcc9f6eb9f196 100644 (file)
--- 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)
+}