From a24228919f12b9c68add647de03d4a8c478653bf Mon Sep 17 00:00:00 2001 From: mwoss Date: Sat, 1 May 2021 01:35:19 +0200 Subject: [PATCH 1/5] Add implementation of exponential function --- decimal.go | 125 ++++++++++++++++++++++++++++++++++++++++++++++++ decimal_test.go | 103 +++++++++++++++++++++++++++++++++++++++ 2 files changed, 228 insertions(+) diff --git a/decimal.go b/decimal.go index 4876f8c9..3f612a61 100644 --- a/decimal.go +++ b/decimal.go @@ -52,6 +52,8 @@ var DivisionPrecision = 16 // silently lose precision. var MarshalJSONWithoutQuotes = false +var ExpMaxIterations = 1000 + // Zero constant, to make computations faster. // Zero should never be compared with == or != directly, please use decimal.Equal or decimal.Cmp instead. var Zero = New(0, 1) @@ -64,6 +66,8 @@ var fiveInt = big.NewInt(5) var tenInt = big.NewInt(10) var twentyInt = big.NewInt(20) +var factorials = []Decimal{New(1, 0)} + // Decimal represents a fixed-point decimal. It is immutable. // number = value * 10 ^ exp type Decimal struct { @@ -630,6 +634,127 @@ func (d Decimal) Pow(d2 Decimal) Decimal { return temp.Mul(temp).Div(d) } +// Exp calculates the natural exponent of decimal (e to the power of d) +func (d Decimal) Exp(overallPrecision uint32) (Decimal, error) { + // Algorithm based on Variable precision exponential function. + // ACM Transactions on Mathematical Software by Hull, T. E., & Abrham, A. + if d.IsZero() { + return Decimal{oneInt, 0}, nil + } + + // special cases preventing under and overflow + overflowThreshold := New(23*int64(overallPrecision), 0) + + // fail if abs(d) beyond an over/underflow threshold + if d.Abs().Cmp(overflowThreshold) > 0 { + // TODO: Increase the precision + return Decimal{}, fmt.Errorf("over/underflow threshold") + } + + overflowThreshold2 := New(9, -int32(overallPrecision)-1) + // return 1 if abs(d) small enough; this also avoids later over/underflow + if d.Abs().Cmp(overflowThreshold2) <= 0 { + return Decimal{oneInt, d.exp}, nil + } + + // t is the smallest integer >= 0 such that the corresponding abs(d/k) < 1 + t := d.exp + int32(d.NumDigits()) // Add d.NumDigits because the paper assumes that d.value [0.1, 1) + + if t < 0 { + t = 0 + } + + k := New(1, t) // reduction factor + r := Decimal{new(big.Int).Set(d.value), d.exp - t} // reduced argument + p := int32(overallPrecision) + t + 2 // precision for calculating the sum + + // determine n, the number of therms for calculating sum + // use first Newton step (1.435p - 1.182) / log10(p/abs(r)) + // for solving appropriate equation, along with directed + // roundings and simple rational bound for log10(p/abs(r)) + rf := r.Abs().InexactFloat64() + pf := float64(p) + nf := math.Ceil((1.453*pf - 1.182) / math.Log10(pf/rf)) + if nf > float64(ExpMaxIterations) || math.IsNaN(nf) { + return Decimal{}, fmt.Errorf("exact value cannot be calculated in <=ExpMaxIterations iterations") + } + n := int64(nf) + + tmp := New(0, 0) + sum := New(1, 0) + one := New(1, 0) + for i := n - 1; i > 0; i-- { + tmp.value.SetInt64(i) + sum = sum.Mul(r.Div(tmp)) + sum = sum.Add(one) + } + + ki := k.IntPart() + + res := New(1, 0) + for i := ki; i > 0; i-- { + res = res.Mul(sum) + } + res = res.Round(-d.exp + 1) + + return res, nil +} + +// ExpTaylor calculates the natural exponent of decimal (e to the power of d) using Taylor series expansion +func (d Decimal) ExpTaylor(precision int32) (Decimal, error) { + // Note(mwoss): Implementation can be optimized by exclusively using big.Int API only + if d.IsZero() { + return Decimal{oneInt, 0}, nil + } + + epsilon := New(1, precision-1) + result := New(1, 0) + decAbs := d.Abs() + pow := d.Abs() + factorial := New(1, 0) + + for i := int64(1); ; { + step := pow.Div(factorial) + result = result.Add(step) + + // Stop Taylor series when current step is smaller than epsilon + if step.Cmp(epsilon) < 0 { + break + } + + pow = pow.Mul(decAbs) + + i++ + + // Calculate next factorial number or retrieve cached value + if len(factorials) >= int(i) { + factorial = factorials[i-1] + } else { + factorial = factorials[i-2].Mul(New(i, 0)) + factorials = append(factorials, factorial) + } + + //fmt.Println(i, result) + } + + //fmt.Println(factorials) + if d.Sign() < 0 { + result = New(1, 0).Div(result) + } + + result = result.Round(-precision) + + //fmt.Println(factorials) + + return result, nil +} + +// NumDigits return the number of digits coefficient (d.Value) +func (d Decimal) NumDigits() int { + // Note(mwoss): It can be optimized, unnecessary cast of bit.Int to string + return len(d.value.String()) +} + // IsInteger returns true when decimal can be represented as an integer value, otherwise, it returns false. func (d Decimal) IsInteger() bool { // The most typical case, all decimal with exponent higher or equal 0 can be represented as integer diff --git a/decimal_test.go b/decimal_test.go index 36677548..55ec5d98 100644 --- a/decimal_test.go +++ b/decimal_test.go @@ -2596,6 +2596,109 @@ func TestDecimal_IsInteger(t *testing.T) { } } +func TestDecimal_Exp(t *testing.T) { + for _, testCase := range []struct { + Dec string + OverallPrecision uint32 + ExpectedDec string + }{ + {"0", 1, "1"}, + {"0.00", 2, "1"}, + {"1.0", 5, ""}, + {"3.0", 2, ""}, + {"5.26", 2, ""}, + {"5.2663117716", 2, ""}, + {"26.1", 2, ""}, + {"529.1591", 2, ""}, + {"-1.0", 2, ""}, + {"-3.0", 2, ""}, + {"-5.26", 2, ""}, + {"-5.2663117716", 2, ""}, + {"-26.1", 2, ""}, + } { + d, _ := NewFromString(testCase.Dec) + expected, _ := NewFromString(testCase.ExpectedDec) + + exp, err := d.Exp(testCase.OverallPrecision) + if err != nil { + t.Fatal(err) + } + if exp.Cmp(expected) != 0 { + t.Errorf("expected %s, got %s", testCase.ExpectedDec, exp.String()) + } + + } +} + +func TestDecimal_ExpTaylor(t *testing.T) { + for _, testCase := range []struct { + Dec string + Precision int32 + ExpectedDec string + }{ + {"0", 1, "1"}, + {"0", -5, "1"}, + {"0.00", 5, "1"}, + {"1.0", -5, "2.71828"}, + {"1.0", -1, "2.7"}, + {"1.0", 5, "0"}, + {"1.0", 1, "0"}, + {"1.0", 0, "3"}, + {"3.0", -1, "20.1"}, + {"3.0", -2, "20.09"}, + //{"5.26", 2, ""}, + //{"5.2663117716", 2, ""}, + {"26.1", -2, "216314672147.06"}, + //{"529.1591", 2, ""}, + //{"-1.0", 2, ""}, + //{"-3.0", 2, ""}, + //{"-5.26", 2, ""}, + //{"-5.2663117716", 2, ""}, + //{"-26.1", 2, ""}, + } { + d, _ := NewFromString(testCase.Dec) + expected, _ := NewFromString(testCase.ExpectedDec) + + exp, err := d.ExpTaylor(testCase.Precision) + if err != nil { + t.Fatal(err) + } + + if exp.Cmp(expected) != 0 { + t.Errorf("expected %s, got %s", testCase.ExpectedDec, exp.String()) + } + + } +} + +func TestDecimal_NumDigits(t *testing.T) { + for _, testCase := range []struct { + Dec string + ExpectedNumDigits int32 + }{ + {"0", 1}, + {"0.00", 2}, + {"1.0", 5}, + {"3.0", 2}, + {"5.26", 2}, + {"5.2663117716", 2}, + {"26.1", 2}, + {"529.1591", 2}, + {"-1.0", 2}, + {"-3.0", 2}, + {"-5.26", 2}, + {"-5.2663117716", 2}, + {"-26.1", 2}, + } { + d, _ := NewFromString(testCase.Dec) + + nums := d.NumDigits() + if nums != 0 { + t.Errorf("expected %d, got %d", testCase.ExpectedNumDigits, nums) + } + } +} + func TestDecimal_Sign(t *testing.T) { if Zero.Sign() != 0 { t.Errorf("%q should have sign 0", Zero) From 87f8d50bef64c4a65a0f155016ef5e4dcf7fbe0b Mon Sep 17 00:00:00 2001 From: mwoss Date: Sun, 30 May 2021 02:14:45 +0200 Subject: [PATCH 2/5] Add implementation of natural exponent function (Taylor, Hull-Abraham) --- decimal.go | 96 +++++++++++++++++++++++---------- decimal_bench_test.go | 26 ++++++++- decimal_test.go | 122 +++++++++++++++++++++++++++++++----------- 3 files changed, 184 insertions(+), 60 deletions(-) diff --git a/decimal.go b/decimal.go index 3f612a61..19ebd1f0 100644 --- a/decimal.go +++ b/decimal.go @@ -634,25 +634,41 @@ func (d Decimal) Pow(d2 Decimal) Decimal { return temp.Mul(temp).Div(d) } -// Exp calculates the natural exponent of decimal (e to the power of d) -func (d Decimal) Exp(overallPrecision uint32) (Decimal, error) { +// ExpHullAbraham calculates the natural exponent of decimal (e to the power of d) using Hull-Abraham algorithm. +// OverallPrecision argument specifies the overall precision of the result (integer part + decimal part). +// +// ExpHullAbraham is faster than ExpTaylor for small precision values, but it is much slower for large precision values. +// +// Example: +// +// NewFromFloat(26.1).ExpHullAbraham(2).String() // output: "220000000000" +// NewFromFloat(26.1).ExpHullAbraham(20).String() // output: "216314672147.05767284" +// +func (d Decimal) ExpHullAbraham(overallPrecision uint32) (Decimal, error) { // Algorithm based on Variable precision exponential function. - // ACM Transactions on Mathematical Software by Hull, T. E., & Abrham, A. + // ACM Transactions on Mathematical Software by T. E. Hull & A. Abrham. if d.IsZero() { return Decimal{oneInt, 0}, nil } - // special cases preventing under and overflow - overflowThreshold := New(23*int64(overallPrecision), 0) + currentPrecision := overallPrecision + + // Algorithm does not work if currentPrecision * 23 < |x|. + // Precision is automatically increased in such cases, so the value can be calculated precisely. + // If newly calculated precision is higher than ExpMaxIterations the currentPrecision will not be changed. + f := d.Abs().InexactFloat64() + if ncp := f / 23; ncp > float64(currentPrecision) && ncp < float64(ExpMaxIterations) { + currentPrecision = uint32(math.Ceil(ncp)) + } // fail if abs(d) beyond an over/underflow threshold + overflowThreshold := New(23*int64(currentPrecision), 0) if d.Abs().Cmp(overflowThreshold) > 0 { - // TODO: Increase the precision - return Decimal{}, fmt.Errorf("over/underflow threshold") + return Decimal{}, fmt.Errorf("over/underflow threshold, exp(x) cannot be calculated precisely") } - overflowThreshold2 := New(9, -int32(overallPrecision)-1) - // return 1 if abs(d) small enough; this also avoids later over/underflow + // Return 1 if abs(d) small enough; this also avoids later over/underflow + overflowThreshold2 := New(9, -int32(currentPrecision)-1) if d.Abs().Cmp(overflowThreshold2) <= 0 { return Decimal{oneInt, d.exp}, nil } @@ -666,9 +682,9 @@ func (d Decimal) Exp(overallPrecision uint32) (Decimal, error) { k := New(1, t) // reduction factor r := Decimal{new(big.Int).Set(d.value), d.exp - t} // reduced argument - p := int32(overallPrecision) + t + 2 // precision for calculating the sum + p := int32(currentPrecision) + t + 2 // precision for calculating the sum - // determine n, the number of therms for calculating sum + // Determine n, the number of therms for calculating sum // use first Newton step (1.435p - 1.182) / log10(p/abs(r)) // for solving appropriate equation, along with directed // roundings and simple rational bound for log10(p/abs(r)) @@ -685,36 +701,66 @@ func (d Decimal) Exp(overallPrecision uint32) (Decimal, error) { one := New(1, 0) for i := n - 1; i > 0; i-- { tmp.value.SetInt64(i) - sum = sum.Mul(r.Div(tmp)) + sum = sum.Mul(r.DivRound(tmp, p)) sum = sum.Add(one) } ki := k.IntPart() - res := New(1, 0) for i := ki; i > 0; i-- { res = res.Mul(sum) } - res = res.Round(-d.exp + 1) + + resNumDigits := int32(res.NumDigits()) + + var roundDigits int32 + if resNumDigits > abs(res.exp) { + roundDigits = int32(currentPrecision) - resNumDigits - res.exp + } else { + roundDigits = int32(currentPrecision) + } + + res = res.Round(roundDigits) return res, nil } -// ExpTaylor calculates the natural exponent of decimal (e to the power of d) using Taylor series expansion +// ExpTaylor calculates the natural exponent of decimal (e to the power of d) using Taylor series expansion. +// Precision argument specifies how precise the result must be (number of digits after decimal point). +// Negative precision is allowed. +// +// ExpTaylor is much faster for large precision values than ExpHullAbraham. +// +// Example: +// +// NewFromFloat(26.1).ExpHullAbraham(2).String() // output: "216314672147.06" +// NewFromFloat(26.1).ExpHullAbraham(20).String() // output: "216314672147.05767284062928674083" +// NewFromFloat(26.1).ExpHullAbraham(-10).String() // output: "220000000000" +// func (d Decimal) ExpTaylor(precision int32) (Decimal, error) { // Note(mwoss): Implementation can be optimized by exclusively using big.Int API only if d.IsZero() { - return Decimal{oneInt, 0}, nil + return Decimal{oneInt, 0}.Round(precision), nil + } + + var epsilon Decimal + var divPrecision int32 + if precision < 0 { + epsilon = New(1, -1) + divPrecision = 8 + } else { + epsilon = New(1, -precision-1) + divPrecision = precision + 1 } - epsilon := New(1, precision-1) - result := New(1, 0) decAbs := d.Abs() pow := d.Abs() factorial := New(1, 0) + result := New(1, 0) + for i := int64(1); ; { - step := pow.Div(factorial) + step := pow.DivRound(factorial, divPrecision) result = result.Add(step) // Stop Taylor series when current step is smaller than epsilon @@ -733,23 +779,17 @@ func (d Decimal) ExpTaylor(precision int32) (Decimal, error) { factorial = factorials[i-2].Mul(New(i, 0)) factorials = append(factorials, factorial) } - - //fmt.Println(i, result) } - //fmt.Println(factorials) if d.Sign() < 0 { - result = New(1, 0).Div(result) + result = New(1, 0).DivRound(result, precision+1) } - result = result.Round(-precision) - - //fmt.Println(factorials) - + result = result.Round(precision) return result, nil } -// NumDigits return the number of digits coefficient (d.Value) +// NumDigits return the number of digits of decimal coefficient (d.Value) func (d Decimal) NumDigits() int { // Note(mwoss): It can be optimized, unnecessary cast of bit.Int to string return len(d.value.String()) diff --git a/decimal_bench_test.go b/decimal_bench_test.go index fd3ddb6f..18e5043a 100644 --- a/decimal_bench_test.go +++ b/decimal_bench_test.go @@ -223,4 +223,28 @@ func BenchmarkDecimal_NewFromString_large_number(b *testing.B) { } } } -} \ No newline at end of file +} + +func BenchmarkDecimal_ExpHullAbraham(b *testing.B) { + b.ResetTimer() + + d := RequireFromString("30.412346346346") + + b.ReportAllocs() + b.ResetTimer() + for i := 0; i < b.N; i++ { + _, _ = d.ExpHullAbraham(10) + } +} + +func BenchmarkDecimal_ExpTaylor(b *testing.B) { + b.ResetTimer() + + d := RequireFromString("30.412346346346") + + b.ReportAllocs() + b.ResetTimer() + for i := 0; i < b.N; i++ { + _, _ = d.ExpTaylor(10) + } +} diff --git a/decimal_test.go b/decimal_test.go index 55ec5d98..d337a01d 100644 --- a/decimal_test.go +++ b/decimal_test.go @@ -2596,33 +2596,60 @@ func TestDecimal_IsInteger(t *testing.T) { } } -func TestDecimal_Exp(t *testing.T) { +func TestDecimal_ExpHullAbraham(t *testing.T) { for _, testCase := range []struct { Dec string OverallPrecision uint32 ExpectedDec string }{ {"0", 1, "1"}, - {"0.00", 2, "1"}, - {"1.0", 5, ""}, - {"3.0", 2, ""}, - {"5.26", 2, ""}, - {"5.2663117716", 2, ""}, - {"26.1", 2, ""}, - {"529.1591", 2, ""}, - {"-1.0", 2, ""}, - {"-3.0", 2, ""}, - {"-5.26", 2, ""}, - {"-5.2663117716", 2, ""}, - {"-26.1", 2, ""}, + {"0.00", 5, "1"}, + {"0.5", 5, "1.6487"}, + {"0.569297", 10, "1.767024397"}, + {"0.569297", 16, "1.76702439654095"}, + {"0.569297", 20, "1.7670243965409496521"}, + {"1.0", 0, "3"}, + {"1.0", 1, "3"}, + {"1.0", 5, "2.7183"}, + {"1.0", 10, "2.718281828"}, + {"3.0", 0, "20"}, + {"3.0", 2, "20"}, + {"5.26", 0, "200"}, + {"5.26", 2, "190"}, + {"5.26", 10, "192.4814913"}, + {"5.2663117716", 2, "190"}, + {"5.2663117716", 10, "193.7002327"}, + {"26.1", 2, "220000000000"}, + {"26.1", 10, "216314672100"}, + {"26.1", 20, "216314672147.05767284"}, + {"50.1591", 10, "6078834887000000000000"}, + {"50.1591", 30, "6078834886623434464595.07937141"}, + {"-0.5", 5, "0.60653"}, + {"-0.569297", 10, "0.5659231429"}, + {"-0.569297", 16, "0.565923142859576"}, + {"-0.569297", 20, "0.56592314285957604443"}, + {"-1.0", 1, "0.4"}, + {"-1.0", 5, "0.36788"}, + {"-3.0", 1, "0"}, + {"-3.0", 2, "0.05"}, + {"-3.0", 10, "0.0497870684"}, + {"-5.26", 2, "0.01"}, + {"-5.26", 10, "0.0051953047"}, + {"-5.2663117716", 2, "0.01"}, + {"-5.2663117716", 10, "0.0051626164"}, + {"-26.1", 2, "0"}, + {"-26.1", 15, "0.000000000004623"}, + {"-50.1591", 10, "0"}, + {"-50.1591", 30, "0.000000000000000000000164505208"}, } { d, _ := NewFromString(testCase.Dec) expected, _ := NewFromString(testCase.ExpectedDec) - exp, err := d.Exp(testCase.OverallPrecision) + exp, err := d.ExpHullAbraham(testCase.OverallPrecision) if err != nil { t.Fatal(err) } + if exp.Cmp(expected) != 0 { t.Errorf("expected %s, got %s", testCase.ExpectedDec, exp.String()) } @@ -2637,24 +2664,52 @@ func TestDecimal_ExpTaylor(t *testing.T) { ExpectedDec string }{ {"0", 1, "1"}, - {"0", -5, "1"}, {"0.00", 5, "1"}, - {"1.0", -5, "2.71828"}, - {"1.0", -1, "2.7"}, - {"1.0", 5, "0"}, - {"1.0", 1, "0"}, + {"0", -1, "0"}, + {"0.5", 5, "1.64872"}, + {"0.569297", 10, "1.7670243965"}, + {"0.569297", 16, "1.7670243965409497"}, + {"0.569297", 20, "1.76702439654094965215"}, {"1.0", 0, "3"}, - {"3.0", -1, "20.1"}, - {"3.0", -2, "20.09"}, - //{"5.26", 2, ""}, - //{"5.2663117716", 2, ""}, - {"26.1", -2, "216314672147.06"}, - //{"529.1591", 2, ""}, - //{"-1.0", 2, ""}, - //{"-3.0", 2, ""}, - //{"-5.26", 2, ""}, - //{"-5.2663117716", 2, ""}, - //{"-26.1", 2, ""}, + {"1.0", 1, "2.7"}, + {"1.0", 5, "2.71828"}, + {"1.0", -1, "0"}, + {"1.0", -5, "0"}, + {"3.0", 1, "20.1"}, + {"3.0", 2, "20.09"}, + {"5.26", 0, "192"}, + {"5.26", 2, "192.48"}, + {"5.26", 10, "192.4814912972"}, + {"5.26", -2, "200"}, + {"5.2663117716", 2, "193.70"}, + {"5.2663117716", 10, "193.7002326701"}, + {"26.1", 2, "216314672147.06"}, + {"26.1", 20, "216314672147.05767284062928674083"}, + {"26.1", -2, "216314672100"}, + {"26.1", -10, "220000000000"}, + {"50.1591", 10, "6078834886623434464595.0793714061"}, + {"-0.5", 5, "0.60653"}, + {"-0.569297", 10, "0.5659231429"}, + {"-0.569297", 16, "0.565923142859576"}, + {"-0.569297", 20, "0.56592314285957604443"}, + {"-1.0", 1, "0.4"}, + {"-1.0", 5, "0.36788"}, + {"-1.0", -1, "0"}, + {"-1.0", -5, "0"}, + {"-3.0", 1, "0.1"}, + {"-3.0", 2, "0.05"}, + {"-3.0", 10, "0.0497870684"}, + {"-5.26", 2, "0.01"}, + {"-5.26", 10, "0.0051953047"}, + {"-5.26", -2, "0"}, + {"-5.2663117716", 2, "0.01"}, + {"-5.2663117716", 10, "0.0051626164"}, + {"-26.1", 2, "0"}, + {"-26.1", 15, "0.000000000004623"}, + {"-26.1", -2, "0"}, + {"-26.1", -10, "0"}, + {"-50.1591", 10, "0"}, + {"-50.1591", 30, "0.000000000000000000000164505208"}, } { d, _ := NewFromString(testCase.Dec) expected, _ := NewFromString(testCase.ExpectedDec) @@ -2667,7 +2722,6 @@ func TestDecimal_ExpTaylor(t *testing.T) { if exp.Cmp(expected) != 0 { t.Errorf("expected %s, got %s", testCase.ExpectedDec, exp.String()) } - } } @@ -3185,3 +3239,9 @@ func ExampleNewFromFloat() { //0.123123123123123 //-10000000000000 } + +func Test2(t *testing.T) { + d1 := RequireFromString("3104646.0") + d2 := RequireFromString("3.0") + t.Log(d1.DivRound(d2, -6)) +} From 346274570b39abceae7cfa80b614c02ebd4c54be Mon Sep 17 00:00:00 2001 From: mwoss Date: Sun, 30 May 2021 02:19:45 +0200 Subject: [PATCH 3/5] Update docs --- decimal.go | 2 ++ 1 file changed, 2 insertions(+) diff --git a/decimal.go b/decimal.go index 19ebd1f0..028ffae5 100644 --- a/decimal.go +++ b/decimal.go @@ -52,6 +52,8 @@ var DivisionPrecision = 16 // silently lose precision. var MarshalJSONWithoutQuotes = false +// ExpMaxIterations specifies the maximum number of iterations needed to calculate +// precise natural exponent value using ExpHullAbraham method. var ExpMaxIterations = 1000 // Zero constant, to make computations faster. From 65bffbbce95163b0305aab7370ab4a8b5de64cf0 Mon Sep 17 00:00:00 2001 From: mwoss Date: Sun, 30 May 2021 02:36:07 +0200 Subject: [PATCH 4/5] Fix failing NumDigits tests --- decimal.go | 5 ++++- decimal_test.go | 23 ++++++++++++----------- 2 files changed, 16 insertions(+), 12 deletions(-) diff --git a/decimal.go b/decimal.go index 028ffae5..1c68119b 100644 --- a/decimal.go +++ b/decimal.go @@ -793,7 +793,10 @@ func (d Decimal) ExpTaylor(precision int32) (Decimal, error) { // NumDigits return the number of digits of decimal coefficient (d.Value) func (d Decimal) NumDigits() int { - // Note(mwoss): It can be optimized, unnecessary cast of bit.Int to string + // Note(mwoss): It can be optimized, unnecessary cast of big.Int to string + if d.IsNegative() { + return len(d.value.String()) - 1 + } return len(d.value.String()) } diff --git a/decimal_test.go b/decimal_test.go index d337a01d..7ec3924f 100644 --- a/decimal_test.go +++ b/decimal_test.go @@ -2728,26 +2728,27 @@ func TestDecimal_ExpTaylor(t *testing.T) { func TestDecimal_NumDigits(t *testing.T) { for _, testCase := range []struct { Dec string - ExpectedNumDigits int32 + ExpectedNumDigits int }{ {"0", 1}, - {"0.00", 2}, - {"1.0", 5}, + {"0.00", 1}, + {"1.0", 2}, {"3.0", 2}, - {"5.26", 2}, - {"5.2663117716", 2}, - {"26.1", 2}, - {"529.1591", 2}, + {"5.26", 3}, + {"5.2663117716", 11}, + {"3164836416948884.2162426426426267863", 35}, + {"26.1", 3}, + {"529.1591", 7}, {"-1.0", 2}, {"-3.0", 2}, - {"-5.26", 2}, - {"-5.2663117716", 2}, - {"-26.1", 2}, + {"-5.26", 3}, + {"-5.2663117716", 11}, + {"-26.1", 3}, } { d, _ := NewFromString(testCase.Dec) nums := d.NumDigits() - if nums != 0 { + if nums != testCase.ExpectedNumDigits { t.Errorf("expected %d, got %d", testCase.ExpectedNumDigits, nums) } } From a82be57b52c19c8f2d1caabd5a0b8f703f9ccbbd Mon Sep 17 00:00:00 2001 From: mwoss Date: Mon, 11 Oct 2021 18:12:55 +0200 Subject: [PATCH 5/5] Improve docs, handle multithreading race condition in exp method --- decimal.go | 36 +++++++++++++++++++++++------------- decimal_bench_test.go | 2 +- decimal_test.go | 14 ++++---------- 3 files changed, 28 insertions(+), 24 deletions(-) diff --git a/decimal.go b/decimal.go index 1c68119b..62a8c769 100644 --- a/decimal.go +++ b/decimal.go @@ -53,7 +53,7 @@ var DivisionPrecision = 16 var MarshalJSONWithoutQuotes = false // ExpMaxIterations specifies the maximum number of iterations needed to calculate -// precise natural exponent value using ExpHullAbraham method. +// precise natural exponent value using ExpHullAbrham method. var ExpMaxIterations = 1000 // Zero constant, to make computations faster. @@ -636,17 +636,17 @@ func (d Decimal) Pow(d2 Decimal) Decimal { return temp.Mul(temp).Div(d) } -// ExpHullAbraham calculates the natural exponent of decimal (e to the power of d) using Hull-Abraham algorithm. +// ExpHullAbrham calculates the natural exponent of decimal (e to the power of d) using Hull-Abraham algorithm. // OverallPrecision argument specifies the overall precision of the result (integer part + decimal part). // -// ExpHullAbraham is faster than ExpTaylor for small precision values, but it is much slower for large precision values. +// ExpHullAbrham is faster than ExpTaylor for small precision values, but it is much slower for large precision values. // // Example: // -// NewFromFloat(26.1).ExpHullAbraham(2).String() // output: "220000000000" -// NewFromFloat(26.1).ExpHullAbraham(20).String() // output: "216314672147.05767284" +// NewFromFloat(26.1).ExpHullAbrham(2).String() // output: "220000000000" +// NewFromFloat(26.1).ExpHullAbrham(20).String() // output: "216314672147.05767284" // -func (d Decimal) ExpHullAbraham(overallPrecision uint32) (Decimal, error) { +func (d Decimal) ExpHullAbrham(overallPrecision uint32) (Decimal, error) { // Algorithm based on Variable precision exponential function. // ACM Transactions on Mathematical Software by T. E. Hull & A. Abrham. if d.IsZero() { @@ -731,13 +731,18 @@ func (d Decimal) ExpHullAbraham(overallPrecision uint32) (Decimal, error) { // Precision argument specifies how precise the result must be (number of digits after decimal point). // Negative precision is allowed. // -// ExpTaylor is much faster for large precision values than ExpHullAbraham. +// ExpTaylor is much faster for large precision values than ExpHullAbrham. // // Example: // -// NewFromFloat(26.1).ExpHullAbraham(2).String() // output: "216314672147.06" -// NewFromFloat(26.1).ExpHullAbraham(20).String() // output: "216314672147.05767284062928674083" -// NewFromFloat(26.1).ExpHullAbraham(-10).String() // output: "220000000000" +// d, err := NewFromFloat(26.1).ExpTaylor(2).String() +// d.String() // output: "216314672147.06" +// +// NewFromFloat(26.1).ExpTaylor(20).String() +// d.String() // output: "216314672147.05767284062928674083" +// +// NewFromFloat(26.1).ExpTaylor(-10).String() +// d.String() // output: "220000000000" // func (d Decimal) ExpTaylor(precision int32) (Decimal, error) { // Note(mwoss): Implementation can be optimized by exclusively using big.Int API only @@ -775,11 +780,15 @@ func (d Decimal) ExpTaylor(precision int32) (Decimal, error) { i++ // Calculate next factorial number or retrieve cached value - if len(factorials) >= int(i) { + if len(factorials) >= int(i) && !factorials[i-1].IsZero() { factorial = factorials[i-1] } else { + // To avoid any race conditions, firstly the zero value is appended to a slice to create + // a spot for newly calculated factorial. After that, the zero value is replaced by calculated + // factorial using the index notation. factorial = factorials[i-2].Mul(New(i, 0)) - factorials = append(factorials, factorial) + factorials = append(factorials, Zero) + factorials[i-1] = factorial } } @@ -791,7 +800,8 @@ func (d Decimal) ExpTaylor(precision int32) (Decimal, error) { return result, nil } -// NumDigits return the number of digits of decimal coefficient (d.Value) +// NumDigits returns the number of digits of the decimal coefficient (d.Value) +// Note: Current implementation is extremely slow for large decimals and/or decimals with large fractional part func (d Decimal) NumDigits() int { // Note(mwoss): It can be optimized, unnecessary cast of big.Int to string if d.IsNegative() { diff --git a/decimal_bench_test.go b/decimal_bench_test.go index 18e5043a..269a9f61 100644 --- a/decimal_bench_test.go +++ b/decimal_bench_test.go @@ -233,7 +233,7 @@ func BenchmarkDecimal_ExpHullAbraham(b *testing.B) { b.ReportAllocs() b.ResetTimer() for i := 0; i < b.N; i++ { - _, _ = d.ExpHullAbraham(10) + _, _ = d.ExpHullAbrham(10) } } diff --git a/decimal_test.go b/decimal_test.go index 7ec3924f..fbd188fd 100644 --- a/decimal_test.go +++ b/decimal_test.go @@ -2596,7 +2596,7 @@ func TestDecimal_IsInteger(t *testing.T) { } } -func TestDecimal_ExpHullAbraham(t *testing.T) { +func TestDecimal_ExpHullAbrham(t *testing.T) { for _, testCase := range []struct { Dec string OverallPrecision uint32 @@ -2645,13 +2645,13 @@ func TestDecimal_ExpHullAbraham(t *testing.T) { d, _ := NewFromString(testCase.Dec) expected, _ := NewFromString(testCase.ExpectedDec) - exp, err := d.ExpHullAbraham(testCase.OverallPrecision) + exp, err := d.ExpHullAbrham(testCase.OverallPrecision) if err != nil { t.Fatal(err) } if exp.Cmp(expected) != 0 { - t.Errorf("expected %s, got %s", testCase.ExpectedDec, exp.String()) + t.Errorf("expected %s, got %s, for decimal %s", testCase.ExpectedDec, exp.String(), testCase.Dec) } } @@ -2749,7 +2749,7 @@ func TestDecimal_NumDigits(t *testing.T) { nums := d.NumDigits() if nums != testCase.ExpectedNumDigits { - t.Errorf("expected %d, got %d", testCase.ExpectedNumDigits, nums) + t.Errorf("expected %d digits for decimal %s, got %d", testCase.ExpectedNumDigits, testCase.Dec, nums) } } } @@ -3240,9 +3240,3 @@ func ExampleNewFromFloat() { //0.123123123123123 //-10000000000000 } - -func Test2(t *testing.T) { - d1 := RequireFromString("3104646.0") - d2 := RequireFromString("3.0") - t.Log(d1.DivRound(d2, -6)) -}