Merge pull request #12 from raverkamp/master

replaced divsion code

decimal.go

@@ -281,18 +281,81 @@ func (d Decimal) Mul(d2 Decimal) Decimal {
 // Div returns d / d2. If it doesn't divide exactly, the result will have
 // DivisionPrecision digits after the decimal point.
 func (d Decimal) Div(d2 Decimal) Decimal {
-	// NOTE(vadim): division is hard, use Rat to do it
+	return d.DivRound(d2, int32(DivisionPrecision))
-	ratNum := d.Rat()
+}
-	ratDenom := d2.Rat()
 
-	quoRat := big.NewRat(0, 1).Quo(ratNum, ratDenom)
+// QuoRem does divsion with remainder
-
+// d.QuoRem(d2,precision) returns quotient q and remainder r such that
-	// HACK(vadim): converting from Rat to Decimal inefficiently for now
+//   d = d2 * q + r, q an integer multiple of 10^(-precision)
-	ret, err := NewFromString(quoRat.FloatString(DivisionPrecision))
+//   0 <= r < abs(d2) * 10 ^(-precision) if d>=0
-	if err != nil {
+//   0 >= r > -abs(d2) * 10 ^(-precision) if d<0
-		panic(err) // this should never happen
+// Note that precision<0 is allowed as input.
+func (d Decimal) QuoRem(d2 Decimal, precision int32) (Decimal, Decimal) {
+	d.ensureInitialized()
+	d2.ensureInitialized()
+	if d2.value.Sign() == 0 {
+		panic("decimal division by 0")
+	}
+	scale := -precision
+	e := int64(d.exp - d2.exp - scale)
+	if e > math.MaxInt32 || e < math.MinInt32 {
+		panic("overflow in decimal QuoRem")
+	}
+	var aa, bb, expo big.Int
+	var scalerest int32
+	// d = a 10^ea
+	// d2 = b 10^eb
+	if e < 0 {
+		aa = *d.value
+		expo.SetInt64(-e)
+		bb.Exp(tenInt, &expo, nil)
+		bb.Mul(d2.value, &bb)
+		scalerest = d.exp
+		// now aa = a
+		//     bb = b 10^(scale + eb - ea)
+	} else {
+		expo.SetInt64(e)
+		aa.Exp(tenInt, &expo, nil)
+		aa.Mul(d.value, &aa)
+		bb = *d2.value
+		scalerest = scale + d2.exp
+		// now aa = a ^ (ea - eb - scale)
+		//     bb = b
+	}
+	var q, r big.Int
+	q.QuoRem(&aa, &bb, &r)
+	dq := Decimal{value: &q, exp: scale}
+	dr := Decimal{value: &r, exp: scalerest}
+	return dq, dr
+}
+
+// DivRound divides and rounds to a given precision
+// i.e. to an integer multiple of 10^(-precision)
+//   for a positive quotient digit 5 is rounded up, away from 0
+//   if the quotient is negative then digit 5 is rounded down, away from 0
+// Note that precision<0 is allowed as input.
+func (d Decimal) DivRound(d2 Decimal, precision int32) Decimal {
+	// QuoRem already checks initialization
+	q, r := d.QuoRem(d2, precision)
+	// the actual rounding decision is based on comparing r*10^precision and d2/2
+	// instead compare 2 r 10 ^precision and d2
+	var rv2 big.Int
+	rv2.Abs(r.value)
+	rv2.Lsh(&rv2, 1)
+	// now rv2 = abs(r.value) * 2
+	r2 := Decimal{value: &rv2, exp: r.exp + precision}
+	// r2 is now 2 * r * 10 ^ precision
+	var c = r2.Cmp(d2.Abs())
+
+	if c < 0 {
+		return q
+	} else {
+		if d.value.Sign()*d2.value.Sign() < 0 {
+			return q.Sub(New(1, -precision))
+		} else {
+			return q.Add(New(1, -precision))
+		}
 	}
-	return ret
 }
 
 // Mod returns d % d2.

decimal_test.go

@@ -4,6 +4,7 @@ import (
 	"encoding/json"
 	"encoding/xml"
 	"math"
+        "math/big"
 	"sort"
 	"strconv"
 	"strings"
@@ -643,7 +644,7 @@ func TestDecimal_Div(t *testing.T) {
 		Inp{"-4612301402398.4753343454", "23.5"}: "-196268144782.9138440146978723",
 	}
 
-	for inp, expected := range inputs {
+	for inp, expectedStr := range inputs {
 		num, err := NewFromString(inp.a)
 		if err != nil {
 			t.FailNow()
@@ -653,10 +654,15 @@ func TestDecimal_Div(t *testing.T) {
 			t.FailNow()
 		}
 		got := num.Div(denom)
-		if got.String() != expected {
+		expected, _ := NewFromString(expectedStr)
-			t.Errorf("expected %s when dividing %v by %v, got %v",
+		if !got.Equals(expected) {
+			t.Errorf("expected %v when dividing %v by %v, got %v",
 				expected, num, denom, got)
 		}
+		got2 := num.DivRound(denom, int32(DivisionPrecision))
+		if !got2.Equals(expected) {
+			t.Errorf("expected %v on DivRound (%v,%v), got %v", expected, num, denom, got2)
+		}
 	}
 
 	type Inp2 struct {
@@ -696,6 +702,258 @@ func TestDecimal_Div(t *testing.T) {
 	}
 }
 
+func TestDecimal_QuoRem(t *testing.T) {
+	type Inp4 struct {
+		d   string
+		d2  string
+		exp int32
+		q   string
+		r   string
+	}
+	cases := []Inp4{
+		Inp4{"10", "1", 0, "10", "0"},
+		Inp4{"1", "10", 0, "0", "1"},
+		Inp4{"1", "4", 2, "0.25", "0"},
+		Inp4{"1", "8", 2, "0.12", "0.04"},
+		Inp4{"10", "3", 1, "3.3", "0.1"},
+		Inp4{"100", "3", 1, "33.3", "0.1"},
+		Inp4{"1000", "10", -3, "0", "1000"},
+		Inp4{"1e-3", "2e-5", 0, "50", "0"},
+		Inp4{"1e-3", "2e-3", 1, "0.5", "0"},
+		Inp4{"4e-3", "0.8", 4, "5e-3", "0"},
+		Inp4{"4.1e-3", "0.8", 3, "5e-3", "1e-4"},
+		Inp4{"-4", "-3", 0, "1", "-1"},
+		Inp4{"-4", "3", 0, "-1", "-1"},
+	}
+
+	for _, inp4 := range cases {
+		d, _ := NewFromString(inp4.d)
+		d2, _ := NewFromString(inp4.d2)
+		prec := inp4.exp
+		q, r := d.QuoRem(d2, prec)
+		expectedQ, _ := NewFromString(inp4.q)
+		expectedR, _ := NewFromString(inp4.r)
+		if !q.Equals(expectedQ) || !r.Equals(expectedR) {
+			t.Errorf("bad QuoRem division %s , %s , %d got %v, %v expected %s , %s",
+				inp4.d, inp4.d2, prec, q, r, inp4.q, inp4.r)
+		}
+		if !d.Equals(d2.Mul(q).Add(r)) {
+			t.Errorf("not fitting: d=%v, d2= %v, prec=%d, q=%v, r=%v",
+				d, d2, prec, q, r)
+		}
+		if !q.Equals(q.Truncate(prec)) {
+			t.Errorf("quotient wrong precision: d=%v, d2= %v, prec=%d, q=%v, r=%v",
+				d, d2, prec, q, r)
+		}
+		if r.Abs().Cmp(d2.Abs().Mul(New(1, -prec))) >= 0 {
+			t.Errorf("remainder too large: d=%v, d2= %v, prec=%d, q=%v, r=%v",
+				d, d2, prec, q, r)
+		}
+		if r.value.Sign()*d.value.Sign() < 0 {
+			t.Errorf("signum of divisor and rest do not match: d=%v, d2= %v, prec=%d, q=%v, r=%v",
+				d, d2, prec, q, r)
+		}
+	}
+}
+
+type DivTestCase struct {
+	d    Decimal
+	d2   Decimal
+	prec int32
+}
+
+func createDivTestCases() []DivTestCase {
+	res := make([]DivTestCase, 0)
+	var n int32 = 5
+	a := []int{1, 2, 3, 6, 7, 10, 100, 14, 5, 400, 0, 1000000, 1000000 + 1, 1000000 - 1}
+	for s := -1; s < 2; s = s + 2 { // 2
+		for s2 := -1; s2 < 2; s2 = s2 + 2 { // 2
+			for e1 := -n; e1 <= n; e1++ { // 2n+1
+				for e2 := -n; e2 <= n; e2++ { // 2n+1
+					var prec int32
+					for prec = -n; prec <= n; prec++ { // 2n+1
+						for _, v1 := range a { // 11
+							for _, v2 := range a { // 11, even if 0 is skipped
+								sign1 := New(int64(s), 0)
+								sign2 := New(int64(s2), 0)
+								d := sign1.Mul(New(int64(v1), int32(e1)))
+								d2 := sign2.Mul(New(int64(v2), int32(e2)))
+								res = append(res, DivTestCase{d, d2, prec})
+							}
+						}
+					}
+				}
+			}
+		}
+	}
+	return res
+}
+
+func TestDecimal_QuoRem2(t *testing.T) {
+	for _, tc := range createDivTestCases() {
+		d := tc.d
+		if sign(tc.d2) == 0 {
+			continue
+		}
+		d2 := tc.d2
+		prec := tc.prec
+		q, r := d.QuoRem(d2, prec)
+		// rule 1: d = d2*q +r
+		if !d.Equals(d2.Mul(q).Add(r)) {
+			t.Errorf("not fitting, d=%v, d2=%v, prec=%d, q=%v, r=%v",
+				d, d2, prec, q, r)
+		}
+		// rule 2: q is integral multiple of 10^(-prec)
+		if !q.Equals(q.Truncate(prec)) {
+			t.Errorf("quotient wrong precision, d=%v, d2=%v, prec=%d, q=%v, r=%v",
+				d, d2, prec, q, r)
+		}
+		// rule 3: abs(r)<abs(d) * 10^(-prec)
+		if r.Abs().Cmp(d2.Abs().Mul(New(1, -prec))) >= 0 {
+			t.Errorf("remainder too large, d=%v, d2=%v, prec=%d, q=%v, r=%v",
+				d, d2, prec, q, r)
+		}
+		// rule 4: r and d have the same sign
+		if r.value.Sign()*d.value.Sign() < 0 {
+			t.Errorf("signum of divisor and rest do not match, "+
+				"d=%v, d2=%v, prec=%d, q=%v, r=%v",
+				d, d2, prec, q, r)
+		}
+	}
+}
+
+// this is the old Div method from decimal
+// Div returns d / d2. If it doesn't divide exactly, the result will have
+// DivisionPrecision digits after the decimal point.
+func (d Decimal) DivOld(d2 Decimal, prec int) Decimal {
+	// NOTE(vadim): division is hard, use Rat to do it
+	ratNum := d.Rat()
+	ratDenom := d2.Rat()
+
+	quoRat := big.NewRat(0, 1).Quo(ratNum, ratDenom)
+
+	// HACK(vadim): converting from Rat to Decimal inefficiently for now
+	ret, err := NewFromString(quoRat.FloatString(prec))
+	if err != nil {
+		panic(err) // this should never happen
+	}
+	return ret
+}
+
+func Benchmark_DivideOriginal(b *testing.B) {
+	tcs := createDivTestCases()
+	b.ResetTimer()
+	for i := 0; i < b.N; i++ {
+		for _, tc := range tcs {
+			d := tc.d
+			if sign(tc.d2) == 0 {
+				continue
+			}
+			d2 := tc.d2
+			prec := tc.prec
+			a := d.DivOld(d2, int(prec))
+			if sign(a) > 2 {
+				panic("dummy panic")
+			}
+		}
+	}
+}
+
+func Benchmark_DivideNew(b *testing.B) {
+	tcs := createDivTestCases()
+	b.ResetTimer()
+	for i := 0; i < b.N; i++ {
+		for _, tc := range tcs {
+			d := tc.d
+			if sign(tc.d2) == 0 {
+				continue
+			}
+			d2 := tc.d2
+			prec := tc.prec
+			a := d.DivRound(d2, prec)
+			if sign(a) > 2 {
+				panic("dummy panic")
+			}
+		}
+	}
+}
+
+func sign(d Decimal) int {
+	return d.value.Sign()
+}
+
+// rules for rounded divide, rounded to integer
+// rounded_divide(d,d2) = q
+// sign q * sign (d/d2) >= 0
+// for d and d2 >0 :
+// q is already rounded
+// q = d/d2 + r , with r > -0.5 and r <= 0.5
+// thus q-d/d2 = r, with r > -0.5 and r <= 0.5
+// and d2 q -d = r d2 with r d2 > -d2/2 and r d2 <= d2/2
+// and 2 (d2 q -d) = x with x > -d2 and x <= d2
+// if we factor in precision then x > -d2 * 10^(-precision) and x <= d2 * 10(-precision)
+
+func TestDecimal_DivRound(t *testing.T) {
+	cases := []struct {
+		d      string
+		d2     string
+		prec   int32
+		result string
+	}{
+		{"2", "2", 0, "1"},
+		{"1", "2", 0, "1"},
+		{"-1", "2", 0, "-1"},
+		{"-1", "-2", 0, "1"},
+		{"1", "-2", 0, "-1"},
+		{"1", "-20", 1, "-0.1"},
+		{"1", "-20", 2, "-0.05"},
+		{"1", "20.0000000000000000001", 1, "0"},
+		{"1", "19.9999999999999999999", 1, "0.1"},
+	}
+	for _, s := range cases {
+		d, _ := NewFromString(s.d)
+		d2, _ := NewFromString(s.d2)
+		result, _ := NewFromString(s.result)
+		prec := s.prec
+		q := d.DivRound(d2, prec)
+		if sign(q)*sign(d)*sign(d2) < 0 {
+			t.Errorf("sign of quotient wrong, got: %v/%v is about %v", d, d2, q)
+		}
+		x := q.Mul(d2).Abs().Sub(d.Abs()).Mul(New(2, 0))
+		if x.Cmp(d2.Abs().Mul(New(1, -prec))) > 0 {
+			t.Errorf("wrong rounding, got: %v/%v prec=%d is about %v", d, d2, prec, q)
+		}
+		if x.Cmp(d2.Abs().Mul(New(-1, -prec))) <= 0 {
+			t.Errorf("wrong rounding, got: %v/%v prec=%d is about %v", d, d2, prec, q)
+		}
+		if !q.Equals(result) {
+			t.Errorf("rounded division wrong %s / %s scale %d = %s, got %v", s.d, s.d2, prec, s.result, q)
+		}
+	}
+}
+
+func TestDecimal_DivRound2(t *testing.T) {
+	for _, tc := range createDivTestCases() {
+		d := tc.d
+		if sign(tc.d2) == 0 {
+			continue
+		}
+		d2 := tc.d2
+		prec := tc.prec
+		q := d.DivRound(d2, prec)
+		if sign(q)*sign(d)*sign(d2) < 0 {
+			t.Errorf("sign of quotient wrong, got: %v/%v is about %v", d, d2, q)
+		}
+		x := q.Mul(d2).Abs().Sub(d.Abs()).Mul(New(2, 0))
+		if x.Cmp(d2.Abs().Mul(New(1, -prec))) > 0 {
+			t.Errorf("wrong rounding, got: %v/%v prec=%d is about %v", d, d2, prec, q)
+		}
+		if x.Cmp(d2.Abs().Mul(New(-1, -prec))) <= 0 {
+			t.Errorf("wrong rounding, got: %v/%v prec=%d is about %v", d, d2, prec, q)
+		}
+	}
+}
+
 func TestDecimal_Mod(t *testing.T) {
 	type Inp struct {
 		a string