From 3934d7607b6db7498e2e36f213cc60d6214afefa Mon Sep 17 00:00:00 2001 From: Denis Bernard Date: Thu, 21 May 2020 16:12:41 +0200 Subject: [PATCH] dec: rename ntz to trailingZeroDigits to use nat-style names --- dec.go | 59 ++++++++++++++++++++++++++++++++++++++++----- dec_arith.go | 2 +- dec_test.go | 14 ----------- decimal.go | 4 ++-- decimal_test.go | 23 ++++++------------ stdlib.go | 64 +++++++++++++++++++++++++++++++++++++++++++++++++ 6 files changed, 127 insertions(+), 39 deletions(-) diff --git a/dec.go b/dec.go index b44d5cb..6c358e5 100644 --- a/dec.go +++ b/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< 0 { + // sticky = x.sticky(uint(l) * _DW) + // } + // return uint64(hi<>(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<> (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 diff --git a/dec_arith.go b/dec_arith.go index daf5471..39a87c2 100644 --- a/dec_arith.go +++ b/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 { diff --git a/dec_test.go b/dec_test.go index e64797f..4f8502b 100644 --- a/dec_test.go +++ b/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: diff --git a/decimal.go b/decimal.go index dd9996b..bd8f681 100644 --- a/decimal.go +++ b/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 diff --git a/decimal_test.go b/decimal_test.go index 7d2786e..7847ec3 100644 --- a/decimal_test.go +++ b/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) + } } diff --git a/stdlib.go b/stdlib.go index 702c28a..84c68ff 100644 --- a/stdlib.go +++ b/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 +}