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
return z
}
-func (x *Float) cmpZero() Accuracy {
- if x.neg {
+func makeAcc(above bool) Accuracy {
+ if above {
return Above
}
return Below
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.
if z.form != finite {
return z
}
- z.setExp(int64(z.exp) + int64(exp))
+ z.setExpAndRound(int64(z.exp)+int64(exp), 0)
return z
}
}
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")
}
}
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
}
// 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
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
}
}
// 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
}
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)
}
}
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
}
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
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
return z.SetInt64(0), Exact
case inf:
- return nil, x.cmpZero()
+ return nil, makeAcc(x.neg)
case nan:
return nil, Undef
return z.SetInt64(0), Exact
case inf:
- return nil, x.cmpZero()
+ return nil, makeAcc(x.neg)
case nan:
return nil, Undef
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
// 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"
}
// 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
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
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
// 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
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 {
}
// x, y != 0
- z.form = finite
z.neg = x.neg
if x.neg == y.neg {
// x + y == x + y
}
}
- // -0 is only possible for -0 + -0
- if z.form == zero {
- z.neg = false
- }
-
return z
}
}
// x, y != 0
- z.form = finite
z.neg = x.neg
if x.neg != y.neg {
// x - (-y) == x + y
}
}
- // -0 is only possible for -0 - 0
- if z.form == zero {
- z.neg = false
- }
-
return z
}
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
}
// TODO(gri) handle Inf separately
return z.SetNaN()
}
+ // x == ±0 || y == ±0
if x.form == zero {
if y.form == zero {
return z.SetNaN()
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
}
// +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 {
return 0
case inf:
m = 2
- case nan:
- panic("unimplemented")
+ default:
+ panic("unreachable")
}
if x.neg {
m = -m
}
}
+func TestFloatArithmeticOverflow(t *testing.T) {
+ for _, test := range []struct {
+ prec uint
+ mode RoundingMode
+ op byte
+ x, y, want string
+ acc Accuracy
+ }{
+ {4, ToNearestEven, '+', "0", "0", "0", Exact}, // smoke test
+ {4, ToNearestEven, '+', "0x.8p0", "0x.8p0", "0x.8p1", Exact}, // smoke test
+
+ {4, ToNearestEven, '+', "0", "0x.8p2147483647", "0x.8p2147483647", Exact},
+ {4, ToNearestEven, '+', "0x.8p2147483500", "0x.8p2147483647", "0x.8p2147483647", Below}, // rounded to zero
+ {4, ToNearestEven, '+', "0x.8p2147483647", "0x.8p2147483647", "+Inf", Above}, // exponent overflow in +
+ {4, ToNearestEven, '+', "-0x.8p2147483647", "-0x.8p2147483647", "-Inf", Below}, // exponent overflow in +
+ {4, ToNearestEven, '-', "-0x.8p2147483647", "0x.8p2147483647", "-Inf", Below}, // exponent overflow in -
+
+ {4, ToZero, '+', "0x.fp2147483647", "0x.8p2147483643", "0x.fp2147483647", Below}, // rounded to zero
+ {4, ToNearestEven, '+', "0x.fp2147483647", "0x.8p2147483643", "+Inf", Above}, // exponent overflow in rounding
+ {4, AwayFromZero, '+', "0x.fp2147483647", "0x.8p2147483643", "+Inf", Above}, // exponent overflow in rounding
+
+ {4, AwayFromZero, '-', "-0x.fp2147483647", "0x.8p2147483644", "-Inf", Below}, // exponent overflow in rounding
+ {4, ToNearestEven, '-', "-0x.fp2147483647", "0x.8p2147483643", "-Inf", Below}, // exponent overflow in rounding
+ {4, ToZero, '-', "-0x.fp2147483647", "0x.8p2147483643", "-0x.fp2147483647", Above}, // rounded to zero
+
+ {4, ToNearestEven, '+', "0", "0x.8p-2147483648", "0x.8p-2147483648", Exact},
+ {4, ToNearestEven, '+', "0x.8p-2147483648", "0x.8p-2147483648", "0x.8p-2147483647", Exact},
+
+ {4, ToNearestEven, '*', "1", "0x.8p2147483647", "0x.8p2147483647", Exact},
+ {4, ToNearestEven, '*', "2", "0x.8p2147483647", "+Inf", Above}, // exponent overflow in *
+ {4, ToNearestEven, '*', "-2", "0x.8p2147483647", "-Inf", Below}, // exponent overflow in *
+
+ {4, ToNearestEven, '/', "0.5", "0x.8p2147483647", "0x.8p-2147483646", Exact},
+ {4, ToNearestEven, '/', "0x.8p0", "0x.8p2147483647", "0x.8p-2147483646", Exact},
+ {4, ToNearestEven, '/', "0x.8p-1", "0x.8p2147483647", "0x.8p-2147483647", Exact},
+ {4, ToNearestEven, '/', "0x.8p-2", "0x.8p2147483647", "0x.8p-2147483648", Exact},
+ {4, ToNearestEven, '/', "0x.8p-3", "0x.8p2147483647", "0", Below}, // exponent underflow in /
+ } {
+ x := makeFloat(test.x)
+ y := makeFloat(test.y)
+ z := new(Float).SetPrec(test.prec).SetMode(test.mode)
+ switch test.op {
+ case '+':
+ z.Add(x, y)
+ case '-':
+ z.Sub(x, y)
+ case '*':
+ z.Mul(x, y)
+ case '/':
+ z.Quo(x, y)
+ default:
+ panic("unreachable")
+ }
+ if got := z.Format('p', 0); got != test.want || z.Acc() != test.acc {
+ t.Errorf(
+ "prec = %d (%s): %s %c %s = %s (%s); want %s (%s)",
+ test.prec, test.mode, x.Format('p', 0), test.op, y.Format('p', 0), got, z.Acc(), test.want, test.acc,
+ )
+ }
+ }
+}
+
// TODO(gri) Add tests that check correctness in the presence of aliasing.
// For rounding modes ToNegativeInf and ToPositiveInf, rounding is affected