Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

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

Merged
merged 5 commits into from
Oct 11, 2021
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
180 changes: 180 additions & 0 deletions decimal.go
Original file line number Diff line number Diff line change
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 @@ -630,6 +636,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
Original file line number Diff line number Diff line change
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
Original file line number Diff line number Diff line change
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