diff --git a/decimal.go b/decimal.go index 8067efd3..acc8a9ba 100644 --- a/decimal.go +++ b/decimal.go @@ -52,6 +52,10 @@ var DivisionPrecision = 16 // silently lose precision. var MarshalJSONWithoutQuotes = false +// ExpMaxIterations specifies the maximum number of iterations needed to calculate +// precise natural exponent value using ExpHullAbrham method. +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 +68,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 { @@ -633,6 +639,180 @@ func (d Decimal) Pow(d2 Decimal) Decimal { return temp.Mul(temp).Div(d) } +// 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). +// +// ExpHullAbrham is faster than ExpTaylor for small precision values, but it is much slower for large precision values. +// +// Example: +// +// NewFromFloat(26.1).ExpHullAbrham(2).String() // output: "220000000000" +// NewFromFloat(26.1).ExpHullAbrham(20).String() // output: "216314672147.05767284" +// +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() { + return Decimal{oneInt, 0}, nil + } + + 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 { + return Decimal{}, fmt.Errorf("over/underflow threshold, exp(x) cannot be calculated precisely") + } + + // 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 + } + + // 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(currentPrecision) + 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.DivRound(tmp, p)) + sum = sum.Add(one) + } + + ki := k.IntPart() + res := New(1, 0) + for i := ki; i > 0; i-- { + res = res.Mul(sum) + } + + 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. +// 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 ExpHullAbrham. +// +// Example: +// +// 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 + if d.IsZero() { + 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 + } + + decAbs := d.Abs() + pow := d.Abs() + factorial := New(1, 0) + + result := New(1, 0) + + for i := int64(1); ; { + step := pow.DivRound(factorial, divPrecision) + 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) && !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, Zero) + factorials[i-1] = factorial + } + } + + if d.Sign() < 0 { + result = New(1, 0).DivRound(result, precision+1) + } + + result = result.Round(precision) + return result, nil +} + +// 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() { + return len(d.value.String()) - 1 + } + 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_bench_test.go b/decimal_bench_test.go index fd3ddb6f..269a9f61 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.ExpHullAbrham(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 b0e5ccba..25b664d4 100644 --- a/decimal_test.go +++ b/decimal_test.go @@ -2596,6 +2596,164 @@ func TestDecimal_IsInteger(t *testing.T) { } } +func TestDecimal_ExpHullAbrham(t *testing.T) { + for _, testCase := range []struct { + Dec string + OverallPrecision uint32 + ExpectedDec string + }{ + {"0", 1, "1"}, + {"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.ExpHullAbrham(testCase.OverallPrecision) + if err != nil { + t.Fatal(err) + } + + if exp.Cmp(expected) != 0 { + t.Errorf("expected %s, got %s, for decimal %s", testCase.ExpectedDec, exp.String(), testCase.Dec) + } + + } +} + +func TestDecimal_ExpTaylor(t *testing.T) { + for _, testCase := range []struct { + Dec string + Precision int32 + ExpectedDec string + }{ + {"0", 1, "1"}, + {"0.00", 5, "1"}, + {"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"}, + {"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) + + 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 int + }{ + {"0", 1}, + {"0.00", 1}, + {"1.0", 2}, + {"3.0", 2}, + {"5.26", 3}, + {"5.2663117716", 11}, + {"3164836416948884.2162426426426267863", 35}, + {"26.1", 3}, + {"529.1591", 7}, + {"-1.0", 2}, + {"-3.0", 2}, + {"-5.26", 3}, + {"-5.2663117716", 11}, + {"-26.1", 3}, + } { + d, _ := NewFromString(testCase.Dec) + + nums := d.NumDigits() + if nums != testCase.ExpectedNumDigits { + t.Errorf("expected %d digits for decimal %s, got %d", testCase.ExpectedNumDigits, testCase.Dec, nums) + } + } +} + func TestDecimal_Sign(t *testing.T) { if Zero.Sign() != 0 { t.Errorf("%q should have sign 0", Zero)