Skip to content

Commit

Permalink
Add implementation of natural exponent function (Taylor, Hull-Abraham) (
Browse files Browse the repository at this point in the history 
#229)

* Add implementation of natural exponent function (Taylor, Hull-Abraham)
* Update docs
* Improve docs, handle multithreading race condition in exp method
  • Loading branch information
@mwoss
mwoss committed Oct 11, 2021
1 parent 9ffd602 commit 4654447
Show file tree
Hide file tree
Showing 3 changed files with 363 additions and 1 deletion.
180 changes: 180 additions & 0 deletions decimal.go
View file
Expand Up @@ -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)
Expand All @@ -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 {
Expand Down Expand Up @@ -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
Expand Down
26 changes: 25 additions & 1 deletion decimal_bench_test.go
View file
Expand Up @@ -223,4 +223,28 @@ 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)
}
}
158 changes: 158 additions & 0 deletions decimal_test.go
View file
Expand Up @@ -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)
Expand Down

0 comments on commit 4654447

Please sign in to comment.