math/big: cleaner handling of exponent under/overflow

Fixed several corner-case bugs and added corresponding tests.

Change-Id: I23096b9caeeff0956f65ab59fa91e168d0e47bb8
Reviewed-on: https://go-review.googlesource.com/7001
Reviewed-by: Alan Donovan <adonovan@google.com>

1 files changed, 100 insertions, 102 deletions

diff --git a/src/math/big/float.go b/src/math/big/float.go
index d716c8ca59..a86471e2a5 100644
--- a/src/math/big/float.go
+++ b/src/math/big/float.go
@@ -154,7 +154,7 @@ func (z *Float) SetPrec(prec uint) *Float {
 		z.prec = 0
 		if z.form == finite {
 			// truncate z to 0
-			z.acc = z.cmpZero()
+			z.acc = makeAcc(z.neg)
 			z.form = zero
 		}
 		return z
@@ -172,8 +172,8 @@ func (z *Float) SetPrec(prec uint) *Float {
 	return z
 }
 
-func (x *Float) cmpZero() Accuracy {
-	if x.neg {
+func makeAcc(above bool) Accuracy {
+	if above {
 		return Above
 	}
 	return Below
@@ -265,22 +265,24 @@ func (x *Float) MantExp(mant *Float) (exp int) {
 	return
 }
 
-// setExp sets the exponent for z.
-// If e < MinExp, z becomes ±0; if e > MaxExp, z becomes ±Inf.
-func (z *Float) setExp(e int64) {
-	if debugFloat && z.form != finite {
-		panic("setExp called for non-finite Float")
+func (z *Float) setExpAndRound(exp int64, sbit uint) {
+	if exp < MinExp {
+		// underflow
+		z.acc = makeAcc(z.neg)
+		z.form = zero
+		return
 	}
-	switch {
-	case e < MinExp:
-		// TODO(gri) check that accuracy is adjusted if necessary
-		z.form = zero // underflow
-	default:
-		z.exp = int32(e)
-	case e > MaxExp:
-		// TODO(gri) check that accuracy is adjusted if necessary
-		z.form = inf // overflow
+
+	if exp > MaxExp {
+		// overflow
+		z.acc = makeAcc(!z.neg)
+		z.form = inf
+		return
 	}
+
+	z.form = finite
+	z.exp = int32(exp)
+	z.round(sbit)
 }
 
 // SetMantExp sets z to mant × 2**exp and and returns z.
@@ -308,7 +310,7 @@ func (z *Float) SetMantExp(mant *Float, exp int) *Float {
 	if z.form != finite {
 		return z
 	}
-	z.setExp(int64(z.exp) + int64(exp))
+	z.setExpAndRound(int64(z.exp)+int64(exp), 0)
 	return z
 }
 
@@ -368,14 +370,14 @@ func (x *Float) validate() {
 	}
 	m := len(x.mant)
 	if m == 0 {
-		panic("nonzero finite x with empty mantissa")
+		panic("nonzero finite number with empty mantissa")
 	}
 	const msb = 1 << (_W - 1)
 	if x.mant[m-1]&msb == 0 {
 		panic(fmt.Sprintf("msb not set in last word %#x of %s", x.mant[m-1], x.Format('p', 0)))
 	}
-	if x.prec <= 0 {
-		panic(fmt.Sprintf("invalid precision %d", x.prec))
+	if x.prec == 0 {
+		panic("zero precision finite number")
 	}
 }
 
@@ -507,7 +509,14 @@ func (z *Float) round(sbit uint) {
 			shrVU(z.mant, z.mant, 1)
 			z.mant[n-1] |= 1 << (_W - 1)
 			// adjust exponent
-			z.exp++
+			if z.exp < MaxExp {
+				z.exp++
+			} else {
+				// exponent overflow
+				z.acc = makeAcc(!z.neg)
+				z.form = inf
+				return
+			}
 		}
 		z.acc = Above
 	}
@@ -515,8 +524,6 @@ func (z *Float) round(sbit uint) {
 	// zero out trailing bits in least-significant word
 	z.mant[0] &^= lsb - 1
 
-	// TODO(gri) can z.mant be all 0s at this point?
-
 	// update accuracy
 	if z.acc != Exact && z.neg {
 		z.acc ^= Below | Above
@@ -655,13 +662,9 @@ func (z *Float) SetInt(x *Int) *Float {
 		return z
 	}
 	// x != 0
-	z.form = finite
 	z.mant = z.mant.set(x.abs)
 	fnorm(z.mant)
-	z.setExp(int64(bits))
-	if z.prec < bits {
-		z.round(0)
-	}
+	z.setExpAndRound(int64(bits), 0)
 	return z
 }
 
@@ -692,7 +695,7 @@ func (z *Float) SetInf(sign int) *Float {
 }
 
 // SetNaN sets z to a NaN value, and returns z.
-// The precision of z is unchanged and the result is always Undef.
+// The precision of z is unchanged and the result accuracy is always Undef.
 func (z *Float) SetNaN() *Float {
 	z.acc = Undef
 	z.form = nan
@@ -711,14 +714,15 @@ func (z *Float) Set(x *Float) *Float {
 	}
 	z.acc = Exact
 	if z != x {
-		if z.prec == 0 {
-			z.prec = x.prec
-		}
 		z.form = x.form
 		z.neg = x.neg
-		z.exp = x.exp
-		z.mant = z.mant.set(x.mant)
-		if z.prec < x.prec {
+		if x.form == finite {
+			z.exp = x.exp
+			z.mant = z.mant.set(x.mant)
+		}
+		if z.prec == 0 {
+			z.prec = x.prec
+		} else if z.prec < x.prec {
 			z.round(0)
 		}
 	}
@@ -738,8 +742,10 @@ func (z *Float) Copy(x *Float) *Float {
 		z.acc = x.acc
 		z.form = x.form
 		z.neg = x.neg
-		z.mant = z.mant.set(x.mant)
-		z.exp = x.exp
+		if z.form == finite {
+			z.mant = z.mant.set(x.mant)
+			z.exp = x.exp
+		}
 	}
 	return z
 }
@@ -821,7 +827,7 @@ func (x *Float) Int64() (int64, Accuracy) {
 	switch x.form {
 	case finite:
 		// 0 < |x| < +Inf
-		acc := x.cmpZero()
+		acc := makeAcc(x.neg)
 		if x.exp <= 0 {
 			// 0 < |x| < 1
 			return 0, acc
@@ -927,7 +933,7 @@ func (x *Float) Int(z *Int) (*Int, Accuracy) {
 	switch x.form {
 	case finite:
 		// 0 < |x| < +Inf
-		acc := x.cmpZero()
+		acc := makeAcc(x.neg)
 		if x.exp <= 0 {
 			// 0 < |x| < 1
 			return z.SetInt64(0), acc
@@ -960,7 +966,7 @@ func (x *Float) Int(z *Int) (*Int, Accuracy) {
 		return z.SetInt64(0), Exact
 
 	case inf:
-		return nil, x.cmpZero()
+		return nil, makeAcc(x.neg)
 
 	case nan:
 		return nil, Undef
@@ -1010,7 +1016,7 @@ func (x *Float) Rat(z *Rat) (*Rat, Accuracy) {
 		return z.SetInt64(0), Exact
 
 	case inf:
-		return nil, x.cmpZero()
+		return nil, makeAcc(x.neg)
 
 	case nan:
 		return nil, Undef
@@ -1035,8 +1041,22 @@ func (z *Float) Neg(x *Float) *Float {
 	return z
 }
 
-// z = x + y, ignoring signs of x and y.
-// x.form and y.form must be finite.
+func validateBinaryOperands(x, y *Float) {
+	if !debugFloat {
+		// avoid performance bugs
+		panic("validateBinaryOperands called but debugFloat is not set")
+	}
+	if len(x.mant) == 0 {
+		panic("empty mantissa for x")
+	}
+	if len(y.mant) == 0 {
+		panic("empty mantissa for y")
+	}
+}
+
+// z = x + y, ignoring signs of x and y for the addition
+// but using the sign of z for rounding the result.
+// x and y must have a non-empty mantissa and valid exponent.
 func (z *Float) uadd(x, y *Float) {
 	// Note: This implementation requires 2 shifts most of the
 	// time. It is also inefficient if exponents or precisions
@@ -1048,8 +1068,8 @@ func (z *Float) uadd(x, y *Float) {
 	// Point Addition With Exact Rounding (as in the MPFR Library)"
 	// http://www.vinc17.net/research/papers/rnc6.pdf
 
-	if debugFloat && (len(x.mant) == 0 || len(y.mant) == 0) {
-		panic("uadd called with empty mantissa")
+	if debugFloat {
+		validateBinaryOperands(x, y)
 	}
 
 	// compute exponents ex, ey for mantissa with "binary point"
@@ -1075,20 +1095,20 @@ func (z *Float) uadd(x, y *Float) {
 	}
 	// len(z.mant) > 0
 
-	z.setExp(ex + int64(len(z.mant))*_W - fnorm(z.mant))
-	z.round(0)
+	z.setExpAndRound(ex+int64(len(z.mant))*_W-fnorm(z.mant), 0)
 }
 
-// z = x - y for x >= y, ignoring signs of x and y.
-// x.form and y.form must be finite.
+// z = x - y for |x| > |y|, ignoring signs of x and y for the subtraction
+// but using the sign of z for rounding the result.
+// x and y must have a non-empty mantissa and valid exponent.
 func (z *Float) usub(x, y *Float) {
 	// This code is symmetric to uadd.
 	// We have not factored the common code out because
 	// eventually uadd (and usub) should be optimized
 	// by special-casing, and the code will diverge.
 
-	if debugFloat && (len(x.mant) == 0 || len(y.mant) == 0) {
-		panic("usub called with empty mantissa")
+	if debugFloat {
+		validateBinaryOperands(x, y)
 	}
 
 	ex := int64(x.exp) - int64(len(x.mant))*_W
@@ -1113,19 +1133,20 @@ func (z *Float) usub(x, y *Float) {
 	if len(z.mant) == 0 {
 		z.acc = Exact
 		z.form = zero
+		z.neg = false
 		return
 	}
 	// len(z.mant) > 0
 
-	z.setExp(ex + int64(len(z.mant))*_W - fnorm(z.mant))
-	z.round(0)
+	z.setExpAndRound(ex+int64(len(z.mant))*_W-fnorm(z.mant), 0)
 }
 
-// z = x * y, ignoring signs of x and y.
-// x.form and y.form must be finite.
+// z = x * y, ignoring signs of x and y for the multiplication
+// but using the sign of z for rounding the result.
+// x and y must have a non-empty mantissa and valid exponent.
 func (z *Float) umul(x, y *Float) {
-	if debugFloat && (len(x.mant) == 0 || len(y.mant) == 0) {
-		panic("umul called with empty mantissa")
+	if debugFloat {
+		validateBinaryOperands(x, y)
 	}
 
 	// Note: This is doing too much work if the precision
@@ -1137,16 +1158,15 @@ func (z *Float) umul(x, y *Float) {
 	e := int64(x.exp) + int64(y.exp)
 	z.mant = z.mant.mul(x.mant, y.mant)
 
-	// normalize mantissa
-	z.setExp(e - fnorm(z.mant))
-	z.round(0)
+	z.setExpAndRound(e-fnorm(z.mant), 0)
 }
 
-// z = x / y, ignoring signs of x and y.
-// x.form and y.form must be finite.
+// z = x / y, ignoring signs of x and y for the division
+// but using the sign of z for rounding the result.
+// x and y must have a non-empty mantissa and valid exponent.
 func (z *Float) uquo(x, y *Float) {
-	if debugFloat && (len(x.mant) == 0 || len(y.mant) == 0) {
-		panic("uquo called with empty mantissa")
+	if debugFloat {
+		validateBinaryOperands(x, y)
 	}
 
 	// mantissa length in words for desired result precision + 1
@@ -1172,13 +1192,8 @@ func (z *Float) uquo(x, y *Float) {
 	// divide
 	var r nat
 	z.mant, r = z.mant.div(nil, xadj, y.mant)
-
-	// determine exponent
 	e := int64(x.exp) - int64(y.exp) - int64(d-len(z.mant))*_W
 
-	// normalize mantissa
-	z.setExp(e - fnorm(z.mant))
-
 	// The result is long enough to include (at least) the rounding bit.
 	// If there's a non-zero remainder, the corresponding fractional part
 	// (if it were computed), would have a non-zero sticky bit (if it were
@@ -1187,15 +1202,16 @@ func (z *Float) uquo(x, y *Float) {
 	if len(r) > 0 {
 		sbit = 1
 	}
-	z.round(sbit)
+
+	z.setExpAndRound(e-fnorm(z.mant), sbit)
 }
 
 // ucmp returns Below, Exact, or Above, depending
-// on whether x < y, x == y, or x > y.
-// x.form and y.form must be finite.
+// on whether |x| < |y|, |x| == |y|, or |x| > |y|.
+// x and y must have a non-empty mantissa and valid exponent.
 func (x *Float) ucmp(y *Float) Accuracy {
-	if debugFloat && (len(x.mant) == 0 || len(y.mant) == 0) {
-		panic("ucmp called with empty mantissa")
+	if debugFloat {
+		validateBinaryOperands(x, y)
 	}
 
 	switch {
@@ -1284,7 +1300,6 @@ func (z *Float) Add(x, y *Float) *Float {
 	}
 
 	// x, y != 0
-	z.form = finite
 	z.neg = x.neg
 	if x.neg == y.neg {
 		// x + y == x + y
@@ -1301,11 +1316,6 @@ func (z *Float) Add(x, y *Float) *Float {
 		}
 	}
 
-	// -0 is only possible for -0 + -0
-	if z.form == zero {
-		z.neg = false
-	}
-
 	return z
 }
 
@@ -1340,7 +1350,6 @@ func (z *Float) Sub(x, y *Float) *Float {
 	}
 
 	// x, y != 0
-	z.form = finite
 	z.neg = x.neg
 	if x.neg != y.neg {
 		// x - (-y) == x + y
@@ -1357,11 +1366,6 @@ func (z *Float) Sub(x, y *Float) *Float {
 		}
 	}
 
-	// -0 is only possible for -0 - 0
-	if z.form == zero {
-		z.neg = false
-	}
-
 	return z
 }
 
@@ -1392,15 +1396,9 @@ func (z *Float) Mul(x, y *Float) *Float {
 		return z
 	}
 
-	if x.form == zero || y.form == zero {
-		z.acc = Exact
-		z.form = zero
-		return z
-	}
-
 	// x, y != 0
-	z.form = finite
 	z.umul(x, y)
+
 	return z
 }
 
@@ -1426,6 +1424,7 @@ func (z *Float) Quo(x, y *Float) *Float {
 			// TODO(gri) handle Inf separately
 			return z.SetNaN()
 		}
+		// x == ±0 || y == ±0
 		if x.form == zero {
 			if y.form == zero {
 				return z.SetNaN()
@@ -1433,16 +1432,14 @@ func (z *Float) Quo(x, y *Float) *Float {
 			z.form = zero
 			return z
 		}
-		// x != 0
-		if y.form == zero {
-			z.form = inf
-			return z
-		}
+		// y == ±0
+		z.form = inf
+		return z
 	}
 
 	// x, y != 0
-	z.form = finite
 	z.uquo(x, y)
+
 	return z
 }
 
@@ -1505,6 +1502,7 @@ func (res cmpResult) Geq() bool     { return res.acc&Below == 0 }
 //	+1 if 0 < x < +Inf
 //	+2 if x == +Inf
 //
+// x must not be NaN.
 func (x *Float) ord() int {
 	var m int
 	switch x.form {
@@ -1514,8 +1512,8 @@ func (x *Float) ord() int {
 		return 0
 	case inf:
 		m = 2
-	case nan:
-		panic("unimplemented")
+	default:
+		panic("unreachable")
 	}
 	if x.neg {
 		m = -m