diff --git a/dec_arith.go b/dec_arith.go index 2140b60..daf5471 100644 --- a/dec_arith.go +++ b/dec_arith.go @@ -19,13 +19,13 @@ const ( _DWb = _DW*100000/30103 + 1 ) -var pow10s = [...]uint64{ +var pow10tab = [...]uint64{ 1, 10, 100, 1000, 10000, 100000, 1000000, 10000000, 100000000, 1000000000, 10000000000, 100000000000, 1000000000000, 10000000000000, 100000000000000, 1000000000000000, 10000000000000000, 100000000000000000, 1000000000000000000, 10000000000000000000, } -func pow10(n uint) Word { return Word(pow10s[n]) } +func pow10(n uint) Word { return Word(pow10tab[n]) } var maxDigits = [...]uint{ 1, 1, 1, 1, 2, 2, 2, 3, 3, 3, 4, 4, 4, 4, 5, 5, diff --git a/dec_arith_amd64.s b/dec_arith_amd64.s index d19225b..beaad84 100644 --- a/dec_arith_amd64.s +++ b/dec_arith_amd64.s @@ -528,7 +528,7 @@ TEXT ·shl10VU(SB),NOSPLIT,$0 JEQ X8c // copy if s = 0 MOVQ $_DW, CX - LEAQ ·pow10s(SB), DI + LEAQ ·pow10tab(SB), DI SUBQ BX, CX MOVQ 0(DI)(BX*8), R12 // m = pow10(s) MOVQ 0(DI)(CX*8), R11 // d = pow10(_DW-s) @@ -582,7 +582,7 @@ TEXT ·shr10VU(SB),NOSPLIT,$0 JEQ X9c // copy if s = 0 MOVQ $_DW, CX - LEAQ ·pow10s(SB), SI + LEAQ ·pow10tab(SB), SI SUBQ BX, CX MOVQ 0(SI)(BX*8), R11 // d = pow10(s) MOVQ 0(SI)(CX*8), R12 // m = pow10(_DW-s) diff --git a/decimal.go b/decimal.go index ae58293..7426775 100644 --- a/decimal.go +++ b/decimal.go @@ -33,9 +33,9 @@ type Decimal struct { neg bool } -// NewDecimal allocates and returns a new Float set to x, -// with precision DefaultDecimalPrec and rounding mode ToNearestEven. -// NewDecimal panics with ErrNaN if x is a NaN. +// NewDecimal allocates and returns a new Float set to x, with precision +// DefaultDecimalPrec and rounding mode ToNearestEven. NewDecimal panics with +// ErrNaN if x is a NaN. func NewDecimal(x float64) *Decimal { if math.IsNaN(x) { panic(ErrNaN{"NewDecimal(NaN)"}) @@ -376,35 +376,25 @@ func (z *Decimal) Copy(x *Decimal) *Decimal { return z } -// Float returns the big.Float value nearest to x with precision prec and -// RoundingMode set to that of x. The returned accuracy is the accuracy if the -// conversion from base 10 to base 2. -// -// If prec is 0, the result precision will be set to the precision of x -// (converting precision in decimal digits to bits). -// -// Note that for high enough exponents, the result might overflow and be set to -// ±Inf. In this case, accuracy will be either Above or Below, depending on the -// sign of x. -func (x *Decimal) Float(prec uint) (*big.Float, Accuracy) { - p := uint64(prec) - if p == 0 { - p = uint64(float64(x.prec)*log2_10) + 1 - if p < 64 { - p = 64 - } +// Float sets z to the (possibly rounded) value of x. If a non-nil *big.Float +// argument z is provided, Float stores the result in z instead of allocating a +// new big.Float. +// If z's precision is 0, it is changed to max(⌈x.Prec() * log2(10)⌉, 64). +func (x *Decimal) Float(z *big.Float) *big.Float { + if z == nil { + z = new(big.Float).SetMode(big.RoundingMode(x.mode)) } - if p > MaxPrec { - p = MaxPrec + p := uint64(z.Prec()) + if p == 0 { + p = uint64(max(int(math.Ceil(float64(x.prec)*log2_10)), 64)) + z.SetPrec(0).SetPrec(uint(p)) } - z := new(big.Float).SetMode(big.RoundingMode(x.mode)).SetPrec(uint(p)) - switch x.form { case zero: - return z, Exact + return z.SetPrec(0).SetPrec(uint(p)) case inf: - return z.SetInf(x.neg), Exact + return z.SetInf(x.neg) } // big.Float has no SetBits. Need to use a temp Int. @@ -412,34 +402,25 @@ func (x *Decimal) Float(prec uint) (*big.Float, Accuracy) { i.SetBits(decToNat(nil, x.mant)) exp := int64(x.exp) - int64(len(x.mant)*_DW) - if exp < MinExp { - // overflow. - return z, makeAcc(x.neg) - } - z = z.SetInt(&i) - // z = x.mant * 2**exp * 5**exp - // Set 2 exponent - z.SetMantExp(z, int(exp)) - - // now multiply/divide by 5**exp - // add a full Word of precision for exponent conversion - tp := ((p+_W-1)/_W + 1) * _W - if tp > MaxPrec { - tp = MaxPrec - } - t := new(big.Float).SetPrec(uint(tp)) - if exp < 0 { - z.Quo(z, pow5Float(t, uint64(-exp))) - } else { - z.Mul(z, pow5Float(t, uint64(exp))) - } - if z.IsInf() { - // inaccurate result - return z, makeAcc(!x.neg) + // now multiply/divide by 10**exp + if exp != 0 { + // add a full Word of precision for exponent conversion + t := new(big.Float).SetPrec(uint(p + _W)) + if exp < 0 { + if exp < MinExp { + // exponent overflow, convert mantissa first + l := uint64(len(x.mant) * _DW) + z.Quo(z, pow10Float(t, l)) + exp += int64(l) + } + z.Quo(z, pow10Float(t, uint64(-exp))) + } else { + z.Mul(z, pow10Float(t, uint64(exp))) + } } - return z, Accuracy(z.Acc()) + return z } // Float32 returns the float32 value nearest to x. If x is too small to be @@ -448,13 +429,13 @@ func (x *Decimal) Float(prec uint) (*big.Float, Accuracy) { // If x is too large to be represented by a float32 (|x| > math.MaxFloat32), // the result is (+Inf, Above) or (-Inf, Below), depending on the sign of x. func (x *Decimal) Float32() (float32, Accuracy) { - z, da := x.Float(32) - f, fa := z.Float32() + z := x.Float(new(big.Float).SetPrec(32)) + f, a := z.Float32() // If big.Float -> float64 conversion is accurate, use Decimal->Float accuracy. - if fa == big.Exact { - fa = big.Accuracy(da) + if a == big.Exact { + a = z.Acc() } - return f, Accuracy(fa) + return f, Accuracy(a) } // Float64 returns the float64 value nearest to x. If x is too small to be @@ -463,13 +444,13 @@ func (x *Decimal) Float32() (float32, Accuracy) { // If x is too large to be represented by a float64 (|x| > math.MaxFloat64), // the result is (+Inf, Above) or (-Inf, Below), depending on the sign of x. func (x *Decimal) Float64() (float64, Accuracy) { - z, da := x.Float(64) - f, fa := z.Float64() + z := x.Float(new(big.Float).SetPrec(64)) + f, a := z.Float64() // If big.Float -> float64 conversion is accurate, use Decimal->Float accuracy. - if fa == big.Exact { - fa = big.Accuracy(da) + if a == big.Exact { + a = z.Acc() } - return f, Accuracy(fa) + return f, Accuracy(a) } func (z *Decimal) GobDecode(buf []byte) error { @@ -510,7 +491,7 @@ func (x *Decimal) Int(z *big.Int) (*big.Int, Accuracy) { acc = Exact } // TODO(db47h): decToNat incurs a copy of its x parameter. - // here we do not care about trashing it. + // Here we do not care about trashing it. z.SetBits(decToNat(z.Bits(), x.intMant())) return z, acc @@ -824,17 +805,64 @@ func (z *Decimal) Set(x *Decimal) *Decimal { return z } +// SetFloat sets z to the (possibly rounded) value of x and returns z. If z's +// precision is 0, it is changed to max(⌈x.MinPrec() * Log10(2)⌉, +// DefaultDecimalPrec). Conversion is performed using x's precision. func (z *Decimal) SetFloat(x *big.Float) *Decimal { - panic("not implemented") + pp := z.prec + // convert using x's precision + z.prec = uint32(math.Ceil(float64(x.Prec()) * log10_2)) + if pp == 0 { + pp = uint32(max(int(z.prec), DefaultDecimalPrec)) + } + z.acc = Exact + z.neg = x.Signbit() + if z.IsInf() { + z.form = inf + z.prec = pp + return z + } + + // TODO(db47h): the conversion is somewhat contrieved, + // but we don't have access to x's mantissa. + f := new(big.Float).Copy(x) + // get int value of mantissa + exp2 := int64(f.MantExp(f)) + fprec := f.MinPrec() + f.SetMantExp(f, int(fprec)) + i, _ := f.Int(nil) + z.SetInt(i) + exp2 -= int64(fprec) + if exp2 != 0 { + t := new(Decimal).SetPrec(uint(z.prec)) + if exp2 < 0 { + if exp2 < MinExp { + // handle exponent overflow. Can only happen if exp2 < 0 + // convert mantissa first, then apply exponent. + exp2 += int64(fprec) + z = z.Quo(z, t.pow2(uint64(fprec))) + } + z = z.Quo(z, t.pow2(uint64(-exp2))) + } else { + z = z.Mul(z, t.pow2(uint64(exp2))) + } + } + z.prec = pp + z.round(0) + return z } -// SetFloat64 sets z to the (possibly rounded) value of x and returns z. -// If z's precision is 0, it is changed to DefaultDecimalPrec (and rounding will have -// no effect). SetFloat64 panics with ErrNaN if x is a NaN. +// SetFloat64 sets z to the (possibly rounded) value of x and returns z. If z's +// precision is 0, it is changed to DefaultDecimalPrec (and rounding will have no effect). +// SetFloat64 panics with ErrNaN if x is a NaN. Conversion is performed using +// x's precision. func (z *Decimal) SetFloat64(x float64) *Decimal { - op := z.prec + pp := z.prec // perform the conversion with 17 decimal digits z.prec = 17 + if pp == 0 { + pp = DefaultDecimalPrec + } if math.IsNaN(x) { panic(ErrNaN{"Decimal.SetFloat64(NaN)"}) } @@ -856,19 +884,16 @@ func (z *Decimal) SetFloat64(x float64) *Decimal { dnorm(z.mant) // multiply / divide by 2**exp if exp2 != 0 { - t := new(Decimal) + t := new(Decimal).SetPrec(uint(z.prec)) if exp2 < 0 { z = z.Quo(z, t.pow2(uint64(-exp2))) } else { z = z.Mul(z, t.pow2(uint64(exp2))) } } - if op == 0 { - op = DefaultDecimalPrec - } - z.SetPrec(uint(op)) + z.prec = pp + z.round(0) return z - } // SetInf sets z to the infinite Decimal -Inf if signbit is @@ -883,33 +908,37 @@ func (z *Decimal) SetInf(signbit bool) *Decimal { } const log2_10 = math.Ln10 / math.Ln2 +const log10_2 = math.Ln2 / math.Ln10 -// SetInt sets z to the (possibly rounded) value of x and returns z. -// If z's precision is 0, it is changed to the larger of x.BitLen() -// or DefaultDecimalPrec (and rounding will have no effect). +// SetInt sets z to the (possibly rounded) value of x and returns z. If z's +// precision is 0, it is changed to max(⌈x.BitLen() * Log10(2)⌉, +// DefaultDecimalPrec) (and rounding will have no effect). func (z *Decimal) SetInt(x *big.Int) *Decimal { bits := uint32(x.BitLen()) - prec := uint32(float64(bits)/log2_10) + 1 // off by 1 at most - // TODO(db47h): adjust precision if needed - if z.prec == 0 { - z.prec = umax32(prec, DefaultDecimalPrec) - } - // TODO(db47h) truncating x could be more efficient if z.prec > 0 - // but small compared to the size of x, or if there are many trailing 0's. z.acc = Exact z.neg = x.Sign() < 0 if bits == 0 { z.form = zero + z.prec = uint32(max(int(z.prec), DefaultDecimalPrec)) return z } // x != 0 - z.mant = z.mant.make((int(prec) + _DW - 1) / _DW).setNat(x.Bits()) + + prec := int(math.Ceil(float64(bits) * log10_2)) // off by 1 at most + // TODO(db47h) truncating x could be more efficient if z.prec > 0 + // but small compared to the size of x. + z.mant = z.mant.make((prec + _DW - 1) / _DW).setNat(x.Bits()) exp := dnorm(z.mant) + if z.prec == 0 { + // adjust precision + z.prec = uint32(max(len(z.mant)*_DW-int(z.mant.ntz()), DefaultDecimalPrec)) + } z.setExpAndRound(int64(len(z.mant))*_DW-exp, 0) return z } func (z *Decimal) setBits64(neg bool, x uint64) *Decimal { + prec := uint32(20) // TODO(db47h): should we adjust to 19 for int64? if z.prec == 0 { z.prec = DefaultDecimalPrec } @@ -923,7 +952,7 @@ func (z *Decimal) setBits64(neg bool, x uint64) *Decimal { z.form = finite z.mant, z.exp = z.mant.setUint64(x) dnorm(z.mant) - if z.prec < 20 { + if z.prec < prec { z.round(0) } return z diff --git a/decimal_test.go b/decimal_test.go index 26ef8d8..6fb57e6 100644 --- a/decimal_test.go +++ b/decimal_test.go @@ -1,6 +1,7 @@ package decimal import ( + "math" "math/big" "math/rand" "reflect" @@ -106,3 +107,16 @@ func BenchmarkDecimal_dnorm(b *testing.B) { benchU = uint(dnorm(d)) } } + +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(128).SetFloat(z) + t.Logf("z = %g -> x = %g", z, x) + x.Float(z) + t.Logf("x = %g -> z = %g", x, z) +} diff --git a/stdlib.go b/stdlib.go index 8a4ffce..3943c42 100644 --- a/stdlib.go +++ b/stdlib.go @@ -280,56 +280,19 @@ func greaterThan(x1, x2, y1, y2 Word) bool { return x1 > y1 || x1 == y1 && x2 > y2 } -// 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. +// pow10 sets z to 10**n and returns z. // n must not be negative. -func pow5Float(z *big.Float, n uint64) *big.Float { - const m = uint64(len(pow5tab) - 1) +func pow10Float(z *big.Float, n uint64) *big.Float { + const m = uint64(len(pow10tab) - 1) if n <= m { - return z.SetUint64(pow5tab[n]) + return z.SetUint64(pow10tab[n]) } // n > m - z.SetUint64(pow5tab[m]) + z.SetUint64(pow10tab[m]) n -= m - f := new(big.Float).SetPrec(z.Prec()).SetUint64(5) + f := new(big.Float).SetPrec(z.Prec() + _W).SetUint64(10) for n > 0 { if n&1 != 0 {