Add implementation of natural exponent function (Taylor, Hull-Abraham) (#229)

* Add implementation of natural exponent function (Taylor, Hull-Abraham)
* Update docs
* Improve docs, handle multithreading race condition in exp method

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

decimal_bench_test.go

@@ -224,3 +224,27 @@ func BenchmarkDecimal_NewFromString_large_number(b *testing.B) {
 		}
 	}
 }
+
+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)
+	}
+}

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)