From 8d267b9b592001da4c7a4d73996f92094fd23152 Mon Sep 17 00:00:00 2001 From: Robert Griesemer Date: Tue, 24 Mar 2015 16:36:16 -0700 Subject: [PATCH] math/big: fixed Float.Float64, implemented Float.Float32 - fix bounds checks for exponent range of denormalized numbers - use correct rounding precision for denormalized numbers - added extra tests Change-Id: I6be56399afd0d9a603300a2e44b5539e08d6f592 Reviewed-on: https://go-review.googlesource.com/8096 Reviewed-by: Alan Donovan --- src/math/big/float.go | 198 ++++++++++++++++++++++++++++++------- src/math/big/float_test.go | 107 +++++++++++++++++--- 2 files changed, 259 insertions(+), 46 deletions(-) diff --git a/src/math/big/float.go b/src/math/big/float.go index fa3751d0c7..629510a18e 100644 --- a/src/math/big/float.go +++ b/src/math/big/float.go @@ -750,6 +750,11 @@ func (z *Float) Copy(x *Float) *Float { return z } +func high32(x nat) uint32 { + // TODO(gri) This can be done more efficiently on 32bit platforms. + return uint32(high64(x) >> 32) +} + func high64(x nat) uint64 { i := len(x) if i == 0 { @@ -872,15 +877,16 @@ func (x *Float) Int64() (int64, Accuracy) { panic("unreachable") } -// Float64 returns the float64 value nearest to x by rounding ToNearestEven -// with 53 bits of precision. -// If x is too small to be represented by a float64 -// (|x| < math.SmallestNonzeroFloat64), the result is (0, Below) or -// (-0, Above), respectively, depending on the sign of x. -// If x is too large to be represented by a float64 (|x| > math.MaxFloat64), +// TODO(gri) Float32 and Float64 are very similar internally but for the +// floatxx parameters and some conversions. Should factor out shared code. + +// Float32 returns the float32 value nearest to x. If x is too small to be +// represented by a float32 (|x| < math.SmallestNonzeroFloat32), the result +// is (0, Below) or (-0, Above), respectively, depending on the sign of x. +// If x is too large to be represented by a float32 (|x| > math.MaxFloat32), // the result is (+Inf, Above) or (-Inf, Below), depending on the sign of x. // The result is (NaN, Undef) for NaNs. -func (x *Float) Float64() (float64, Accuracy) { +func (x *Float) Float32() (float32, Accuracy) { if debugFloat { x.validate() } @@ -888,61 +894,183 @@ func (x *Float) Float64() (float64, Accuracy) { switch x.form { case finite: // 0 < |x| < +Inf + + const ( + fbits = 32 // float size + mbits = 23 // mantissa size (excluding implicit msb) + ebits = fbits - mbits - 1 // 8 exponent size + bias = 1<<(ebits-1) - 1 // 127 exponent bias + dmin = 1 - bias - mbits // -149 smallest unbiased exponent (denormal) + emin = 1 - bias // -126 smallest unbiased exponent (normal) + emax = bias // 127 largest unbiased exponent (normal) + ) + + // Float mantissae m have an explicit msb and are in the range 0.5 <= m < 1.0. + // floatxx mantissae have an implicit msb and are in the range 1.0 <= m < 2.0. + // For a given mantissa m, we need to add 1 to a floatxx exponent to get the + // corresponding Float exponent. + // (see also implementation of math.Ldexp for similar code) + + if x.exp < dmin+1 { + // underflow + if x.neg { + var z float32 + return -z, Above + } + return 0.0, Below + } + // x.exp >= dmin+1 + var r Float - r.prec = 53 + r.prec = mbits + 1 // +1 for implicit msb + if x.exp < emin+1 { + // denormal number - round to fewer bits + r.prec = uint32(x.exp - dmin) + } r.Set(x) - // Rounding via Set may have caused r to overflow - // to ±Inf (rounding never causes underflows to 0). + // Rounding may have caused r to overflow to ±Inf + // (rounding never causes underflows to 0). if r.form == inf { - r.exp = 10000 // cause overflow below + r.exp = emax + 2 // cause overflow below } - // see also implementation of math.Ldexp + if r.exp > emax+1 { + // overflow + if x.neg { + return float32(math.Inf(-1)), Below + } + return float32(math.Inf(+1)), Above + } + // dmin+1 <= r.exp <= emax+1 + + var s uint32 + if r.neg { + s = 1 << (fbits - 1) + } + + m := high32(r.mant) >> ebits & (1< math.MaxFloat64), +// the result is (+Inf, Above) or (-Inf, Below), depending on the sign of x. +// The result is (NaN, Undef) for NaNs. +func (x *Float) Float64() (float64, Accuracy) { + if debugFloat { + x.validate() + } + + switch x.form { + case finite: + // 0 < |x| < +Inf + + const ( + fbits = 64 // float size + mbits = 52 // mantissa size (excluding implicit msb) + ebits = fbits - mbits - 1 // 11 exponent size + bias = 1<<(ebits-1) - 1 // 1023 exponent bias + dmin = 1 - bias - mbits // -1074 smallest unbiased exponent (denormal) + emin = 1 - bias // -1022 smallest unbiased exponent (normal) + emax = bias // 1023 largest unbiased exponent (normal) + ) + + // Float mantissae m have an explicit msb and are in the range 0.5 <= m < 1.0. + // floatxx mantissae have an implicit msb and are in the range 1.0 <= m < 2.0. + // For a given mantissa m, we need to add 1 to a floatxx exponent to get the + // corresponding Float exponent. + // (see also implementation of math.Ldexp for similar code) + + if x.exp < dmin+1 { // underflow if x.neg { - z := 0.0 + var z float64 return -z, Above } return 0.0, Below } - // e > -52 + // x.exp >= dmin+1 + + var r Float + r.prec = mbits + 1 // +1 for implicit msb + if x.exp < emin+1 { + // denormal number - round to fewer bits + r.prec = uint32(x.exp - dmin) + } + r.Set(x) + + // Rounding may have caused r to overflow to ±Inf + // (rounding never causes underflows to 0). + if r.form == inf { + r.exp = emax + 2 // cause overflow below + } - if e >= 2047 { + if r.exp > emax+1 { // overflow if x.neg { return math.Inf(-1), Below } return math.Inf(+1), Above } - // -52 < e < 2047 - - denormal := false - if e < 0 { - denormal = true - e += 52 - } - // 0 < e < 2047 + // dmin+1 <= r.exp <= emax+1 - s := uint64(0) + var s uint64 if r.neg { - s = 1 << 63 + s = 1 << (fbits - 1) } - m := high64(r.mant) >> 11 & (1<<52 - 1) // cut off msb (implicit 1 bit) - z := math.Float64frombits(s | uint64(e)<<52 | m) - if denormal { - // adjust for denormal - // TODO(gri) does this change accuracy? - z /= 1 << 52 + + m := high64(r.mant) >> ebits & (1<