intel decimal library reference tests + speed-ups

d64/decParts.go

@@ -187,27 +187,52 @@ func (dp *decParts) unpackV2(d Decimal) {
 }
 
 // https://en.wikipedia.org/wiki/Decimal64_floating-point_format#Binary_integer_significand_field
-var flavMap = [...]flavor{
-	/* 0000xx */ flNormal53, flNormal53, flNormal53, flNormal53,
-	/* 0001xx */ flNormal53, flNormal53, flNormal53, flNormal53,
-	/* 0010xx */ flNormal53, flNormal53, flNormal53, flNormal53,
-	/* 0011xx */ flNormal53, flNormal53, flNormal53, flNormal53,
-	/* 0100xx */ flNormal53, flNormal53, flNormal53, flNormal53,
-	/* 0101xx */ flNormal53, flNormal53, flNormal53, flNormal53,
-	/* 0110xx */ flNormal53, flNormal53, flNormal53, flNormal53,
-	/* 0111xx */ flNormal53, flNormal53, flNormal53, flNormal53,
-	/* 1000xx */ flNormal53, flNormal53, flNormal53, flNormal53,
-	/* 1001xx */ flNormal53, flNormal53, flNormal53, flNormal53,
-	/* 1010xx */ flNormal53, flNormal53, flNormal53, flNormal53,
-	/* 1011xx */ flNormal53, flNormal53, flNormal53, flNormal53,
-	/* 1100xx */ flNormal51, flNormal51, flNormal51, flNormal51,
-	/* 1101xx */ flNormal51, flNormal51, flNormal51, flNormal51,
-	/* 1110xx */ flNormal51, flNormal51, flNormal51, flNormal51,
-	/* 11110x */ flInf, flInf,
-	/* 111110 */ flQNaN,
-	/* 111111 */ flSNaN,
-}
+var flavMap = func() [128]flavor {
+	fm := [...]flavor{
+		/* 0000xx */ flNormal53, flNormal53, flNormal53, flNormal53,
+		/* 0001xx */ flNormal53, flNormal53, flNormal53, flNormal53,
+		/* 0010xx */ flNormal53, flNormal53, flNormal53, flNormal53,
+		/* 0011xx */ flNormal53, flNormal53, flNormal53, flNormal53,
+		/* 0100xx */ flNormal53, flNormal53, flNormal53, flNormal53,
+		/* 0101xx */ flNormal53, flNormal53, flNormal53, flNormal53,
+		/* 0110xx */ flNormal53, flNormal53, flNormal53, flNormal53,
+		/* 0111xx */ flNormal53, flNormal53, flNormal53, flNormal53,
+		/* 1000xx */ flNormal53, flNormal53, flNormal53, flNormal53,
+		/* 1001xx */ flNormal53, flNormal53, flNormal53, flNormal53,
+		/* 1010xx */ flNormal53, flNormal53, flNormal53, flNormal53,
+		/* 1011xx */ flNormal53, flNormal53, flNormal53, flNormal53,
+		/* 1100xx */ flNormal51, flNormal51, flNormal51, flNormal51,
+		/* 1101xx */ flNormal51, flNormal51, flNormal51, flNormal51,
+		/* 1110xx */ flNormal51, flNormal51, flNormal51, flNormal51,
+		/* 11110x */ flInf, flInf,
+		/* 111110 */ flQNaN,
+		/* 111111 */ flSNaN,
+		/* 0000xx */ flNormal53, flNormal53, flNormal53, flNormal53,
+		/* 0001xx */ flNormal53, flNormal53, flNormal53, flNormal53,
+		/* 0010xx */ flNormal53, flNormal53, flNormal53, flNormal53,
+		/* 0011xx */ flNormal53, flNormal53, flNormal53, flNormal53,
+		/* 0100xx */ flNormal53, flNormal53, flNormal53, flNormal53,
+		/* 0101xx */ flNormal53, flNormal53, flNormal53, flNormal53,
+		/* 0110xx */ flNormal53, flNormal53, flNormal53, flNormal53,
+		/* 0111xx */ flNormal53, flNormal53, flNormal53, flNormal53,
+		/* 1000xx */ flNormal53, flNormal53, flNormal53, flNormal53,
+		/* 1001xx */ flNormal53, flNormal53, flNormal53, flNormal53,
+		/* 1010xx */ flNormal53, flNormal53, flNormal53, flNormal53,
+		/* 1011xx */ flNormal53, flNormal53, flNormal53, flNormal53,
+		/* 1100xx */ flNormal51, flNormal51, flNormal51, flNormal51,
+		/* 1101xx */ flNormal51, flNormal51, flNormal51, flNormal51,
+		/* 1110xx */ flNormal51, flNormal51, flNormal51, flNormal51,
+		/* 11110x */ flInf, flInf,
+		/* 111110 */ flQNaN,
+		/* 111111 */ flSNaN,
+	}
+	return fm
+}()
 
 func (d Decimal) flavor() flavor {
-	return flavMap[int(d.bits>>(64-7))%len(flavMap)]
+	return flavMap[int(d.bits>>(64-7))]
+}
+
+func (d Decimal) isFinite() bool {
+	return d.bits>>(64-5)&0b1111 < 0b1111
 }

d64/decimal.go

@@ -450,7 +450,7 @@ func (d Decimal) Signbit() bool {
 
 func (d Decimal) ScaleB(e Decimal) Decimal {
 	var dp, ep decParts
-	if nan, is := checkNan(d, e, &dp, &ep); is {
+	if nan, is := checkNan2(d, e, &dp, &ep); is {
 		return nan
 	}
 
@@ -534,23 +534,39 @@ func (d Decimal) Class() string {
 	return "+Normal-Normal"[7*dp.sign : 7*(dp.sign+1)]
 }
 
-func checkNan(d, e Decimal, dp, ep *decParts) (Decimal, bool) {
+func checkFinite2(d, e Decimal, dp, ep *decParts) bool {
+	if d.isFinite() && e.isFinite() {
+		dp.fl = d.flavor()
+		ep.fl = e.flavor()
+		dp.unpackV2(d)
+		ep.unpackV2(e)
+		return true
+	}
+	return false
+}
+
+func checkNan2(d, e Decimal, dp, ep *decParts) (Decimal, bool) {
 	dp.fl = d.flavor()
 	ep.fl = e.flavor()
 	switch {
 	case dp.fl == flSNaN:
-		return d, true
 	case ep.fl == flSNaN:
-		return e, true
+		d = e
 	case dp.fl == flQNaN:
-		return d, true
 	case ep.fl == flQNaN:
-		return e, true
+		d = e
 	default:
 		dp.unpackV2(d)
 		ep.unpackV2(e)
 		return Decimal{}, false
 	}
+	return d.qNan(), true
+}
+
+func (d Decimal) qNan() Decimal {
+	// Remove the exponent bits from the NaN.
+	// This is a hack to make sure that the NaN is not interpreted as a normal number.
+	return newDec(d.bits &^ (0xff << 50))
 }
 
 // checkNan3 returns the decimal NaN that is to be propogated and true else first decimal and false
@@ -560,21 +576,20 @@ func checkNan3(d, e, f Decimal, dp, ep, fp *decParts) (Decimal, bool) {
 	fp.fl = f.flavor()
 	switch {
 	case dp.fl == flSNaN:
-		return d, true
 	case ep.fl == flSNaN:
-		return e, true
+		d = e
 	case fp.fl == flSNaN:
-		return f, true
+		d = f
 	case dp.fl == flQNaN:
-		return d, true
 	case ep.fl == flQNaN:
-		return e, true
+		d = e
 	case fp.fl == flQNaN:
-		return f, true
+		d = f
 	default:
 		dp.unpackV2(d)
 		ep.unpackV2(e)
 		fp.unpackV2(f)
 		return Decimal{}, false
 	}
+	return d.qNan(), true
 }

d64/decimal_debug.go

@@ -9,12 +9,12 @@ import "fmt"
 // It uses the binary representation method.
 // Decimal is intentionally a struct to ensure users don't accidentally cast it to uint64.
 type Decimal struct {
+	bits        uint64
 	s           string
 	fl          flavor
 	sign        int8
 	exp         int16
 	significand uint64
-	bits        uint64
 }
 
 func newDec(bits uint64) Decimal {

d64/divconst_test.go

@@ -0,0 +1,86 @@
+package d64
+
+import (
+	"fmt"
+	"math/bits"
+	"math/rand"
+	"testing"
+)
+
+const (
+	base    = 10_000_000_000_000_000
+	tenBase = 10 * base
+	limit   = tenBase * tenBase
+	limitHi = limit / (1 << 64)
+	limitLo = limit % (1 << 64)
+)
+
+func TestU128_divrem_10_15(t *testing.T) {
+	t.Parallel()
+	testU128_divrem(t, u128_div_10_15, 1_000_000_000_000_000)
+}
+
+func TestU128_divrem_10_16(t *testing.T) {
+	t.Parallel()
+	testU128_divrem(t, u128_div_10_16, 10_000_000_000_000_000)
+}
+
+func testU128_divrem(t *testing.T, div func(hi, lo uint64) uint64, d uint64) {
+	t.Helper()
+
+	test := func(hi, lo uint64) func(t *testing.T) {
+		return func(t *testing.T) {
+			t.Helper()
+			q := div(hi, lo)
+			Q, _ := bits.Div64(hi, lo, d)
+			if q != Q {
+				t.Fatalf("U128_div_10_16(%016x_%016x) = %d (%016[3]x), want %d (%016[4]x)",
+					hi, lo, q, Q)
+			}
+		}
+	}
+	t.Run("1", test(0, 0))
+	t.Run("2", test(0, 1))
+	t.Run("3", test(0, 10_000_000_000_000_000))
+	t.Run("4", test(0, 20_000_000_000_000_000))
+	t.Run("5", test(
+		10_000_000_000_000_000_000_000/(1<<64),
+		10_000_000_000_000_000_000_000%(1<<64),
+	))
+	t.Run("6", test(1, 20_000_000_000_000_000))
+
+	r := rand.New(rand.NewSource(0))
+	n := 1_000_000
+	if testing.Short() {
+		n = 100_000
+	}
+	for i := 0; i < n; i++ {
+		// Ensure hi:lo < 10^32.
+		hi := r.Uint64() % limitHi
+		lo := r.Uint64() % limitLo
+		if hi == limitHi {
+			lo %= limitLo
+		}
+		t.Run(fmt.Sprintf("rand[%d]", i), test(hi, lo))
+	}
+}
+
+var globalUint64 uint64
+
+func BenchmarkDiv64_10_16(b *testing.B) {
+	b.ReportAllocs()
+	for i := 0; i < b.N; i++ {
+		hi := uint64(i) * uint64(5421010862428)
+		lo := uint64(i)
+		globalUint64 = u128_div_10_16(hi, lo)
+	}
+}
+
+func BenchmarkBitsDiv64_10_16(b *testing.B) {
+	b.ReportAllocs()
+	for i := 0; i < b.N; i++ {
+		hi := uint64(i)
+		lo := uint64(i)
+		globalUint64, _ = bits.Div64(hi, lo, 10_000_000_000_000_000)
+	}
+}