dec: rename ntz to trailingZeroDigits to use nat-style names

dec.go

@@ -46,13 +46,16 @@ func (x dec) digits() uint {
 	return 0
 }
 
-func (x dec) ntz() uint {
-	for i, w := range x {
-		if w != 0 {
-			return uint(i)*_DW + decTrailingZeros(uint(w))
-		}
+func (x dec) trailingZeroDigits() uint {
+	if len(x) == 0 {
+		return 0
 	}
-	return 0
+	var i uint
+	for x[i] == 0 {
+		i++
+	}
+	// x[i] != 0
+	return i*_DW + trailingZeroDigits(uint(x[i]))
 }
 
 func (x dec) digit(i uint) uint {
@@ -108,6 +111,50 @@ func (z dec) setUint64(x uint64) (dec, int32) {
 	return z.norm(), dig
 }
 
+func (x dec) msd64() (uint64, uint) {
+	// l := len(x)
+	// if l == 0 {
+	// 	return 0, 0
+	// }
+	// var sticky uint
+	// if _W == 64 {
+	// 	var hi, lo Word
+	// 	l -= 2
+	// 	lo = x[l+1]
+	// 	if l > 0 {
+	// 		hi, lo = mulAddWWW_g(lo, _DB, x[l])
+	// 	} else {
+	// 		hi, lo = lo, 0
+	// 	}
+	// 	// normalize result
+	// 	sh := bits.LeadingZeros64(uint64(hi))
+	// 	if lo<<sh != 0 {
+	// 		sticky = 1
+	// 	} else if l > 0 {
+	// 		sticky = x.sticky(uint(l) * _DW)
+	// 	}
+	// 	return uint64(hi<<sh | lo>>(64-sh)), sticky
+	// } else {
+	// 	panic("not implemented")
+	// }
+	t := decToNat(nil, x)
+	i := len(t)
+	for i > 0 && t[i-1] == 0 {
+		i--
+	}
+	t = t[:i]
+	sh := bits.LeadingZeros(uint(t[len(t)-1]))
+	if sh != 0 {
+		var msb big.Word
+		for i := 0; i < len(t); i++ {
+			w := t[i]
+			t[i] = w<<sh | msb
+			msb = w >> (64 - sh)
+		}
+	}
+	return uint64(t[len(t)-1]), 0
+}
+
 // toUint64 returns the low 64 bits of z or MaxUint64 and true if z <= MaxUint64.
 func (x dec) toUint64() (uint64, bool) {
 	// using a decToNat style loop would modify x

dec_arith.go

@@ -67,7 +67,7 @@ func nlz10(x Word) uint {
 	return _DW - decDigits64(uint64(x))
 }
 
-func decTrailingZeros(n uint) uint {
+func trailingZeroDigits(n uint) uint {
 	var d uint
 	if bits.UintSize > 32 {
 		if uint64(n)%10000000000000000 == 0 {

dec_test.go

@@ -566,20 +566,6 @@ func BenchmarkDecExpNN(b *testing.B) {
 	}
 }
 
-func BenchmarkDecExp3Power(b *testing.B) {
-	const x = 3
-	for _, y := range []Word{
-		0x10, 0x40, 0x100, 0x400, 0x1000, 0x4000, 0x10000, 0x40000, 0x100000, 0x400000,
-	} {
-		b.Run(fmt.Sprintf("%#x", y), func(b *testing.B) {
-			var z dec
-			for i := 0; i < b.N; i++ {
-				z.expWW(x, y)
-			}
-		})
-	}
-}
-
 func decFibo(n int) dec {
 	switch n {
 	case 0:

decimal.go

@@ -648,7 +648,7 @@ func (x *Decimal) MinPrec() uint {
 	if x.form != finite {
 		return 0
 	}
-	return uint(len(x.mant))*_DW - x.mant.ntz()
+	return uint(len(x.mant))*_DW - x.mant.trailingZeroDigits()
 }
 
 // Mode returns the rounding mode of x.
@@ -967,7 +967,7 @@ func (z *Decimal) SetInt(x *big.Int) *Decimal {
 	exp := dnorm(z.mant)
 	if z.prec == 0 {
 		// adjust precision
-		z.prec = uint32(max(len(z.mant)*_DW-int(z.mant.ntz()), DefaultDecimalPrec))
+		z.prec = uint32(max(len(z.mant)*_DW-int(z.mant.trailingZeroDigits()), DefaultDecimalPrec))
 	}
 	z.setExpAndRound(int64(len(z.mant))*_DW-exp, 0)
 	return z

decimal_test.go

@@ -1,7 +1,6 @@
 package decimal
 
 import (
-	"math"
 	"math/big"
 	"math/rand"
 	"reflect"
@@ -115,19 +114,11 @@ func BenchmarkDecimal_Sqrt(b *testing.B) {
 	}
 }
 
-func TestDecimal_String(t *testing.T) {
-	z := big.NewFloat(1.1)
-	t.Logf("z: %g!", z)
-	// x, _ := new(Decimal).SetPrec(34).SetString("9223372036854775808")
-	x := new(Decimal).SetPrec(17).SetFloat64(math.MaxFloat64)
-	// t.Log(x.Text('p', 0))
-	t.Logf("x: %g!", x)
-	x.SetPrec(17).SetFloat(z)
-	t.Logf("z = %g -> x (%d) = %g", z, x.Prec(), x)
-	x.Float(z)
-	t.Logf("x = %g -> z = %g", x, z)
-	x.Sqrt(x.SetPrec(33))
-	t.Logf("x: %g!", x)
-	t.Logf("x: %.*f!", x.Prec()-1, x)
-	t.Log("   1.41421356237309504880168872420969807856967187537694807317667974")
+func BenchmarkDecimal_Float(b *testing.B) {
+	d := new(Decimal).SetPrec(100).SetUint64(2)
+	d.Sqrt(d)
+	f := d.Float(nil)
+	for i := 0; i < b.N; i++ {
+		d.Float(f)
+	}
 }

stdlib.go

@@ -326,3 +326,67 @@ func bigEndianWord(buf []byte) Word {
 	}
 	return Word(binary.BigEndian.Uint32(buf))
 }
+
+// These powers of 5 fit into a uint64.
+//
+//	for p, q := uint64(0), uint64(1); p < q; p, q = q, q*5 {
+//		fmt.Println(q)
+//	}
+//
+var pow5tab = [...]uint64{
+	1,
+	5,
+	25,
+	125,
+	625,
+	3125,
+	15625,
+	78125,
+	390625,
+	1953125,
+	9765625,
+	48828125,
+	244140625,
+	1220703125,
+	6103515625,
+	30517578125,
+	152587890625,
+	762939453125,
+	3814697265625,
+	19073486328125,
+	95367431640625,
+	476837158203125,
+	2384185791015625,
+	11920928955078125,
+	59604644775390625,
+	298023223876953125,
+	1490116119384765625,
+	7450580596923828125,
+}
+
+// pow5 sets z to 5**n and returns z.
+// n must not be negative.
+func floatPow5(z *big.Float, n uint64) *big.Float {
+	const m = uint64(len(pow5tab) - 1)
+	if n <= m {
+		return z.SetUint64(pow5tab[n])
+	}
+	// n > m
+
+	z.SetUint64(pow5tab[m])
+	n -= m
+
+	// use more bits for f than for z
+	// TODO(gri) what is the right number?
+	f := new(big.Float).SetPrec(z.Prec() + 64).SetUint64(5)
+
+	for n > 0 {
+		if n&1 != 0 {
+			z.Mul(z, f)
+		}
+		f.Mul(f, f)
+		n >>= 1
+	}
+
+	return z
+}