From 12e9241da2830ba4b0d5fc5959db09010620b369 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mateusz=20Wo=C5=9B?= Date: Fri, 29 Dec 2023 01:59:48 +0100 Subject: [PATCH] Add implementation of natural logarithm (#339) * Add initial implementation of natural logarithm * Add constApproxmation struct to represent mathematical constants with their approximations --- const.go | 63 ++++++++++++++++++++++++++ const_test.go | 34 ++++++++++++++ decimal.go | 117 ++++++++++++++++++++++++++++++++++++++++++++++++ decimal_test.go | 61 +++++++++++++++++++++++++ 4 files changed, 275 insertions(+) create mode 100644 const.go create mode 100644 const_test.go diff --git a/const.go b/const.go new file mode 100644 index 00000000..e5d6fa87 --- /dev/null +++ b/const.go @@ -0,0 +1,63 @@ +package decimal + +import ( + "strings" +) + +const ( + strLn10 = "2.302585092994045684017991454684364207601101488628772976033327900967572609677352480235997205089598298341967784042286248633409525465082806756666287369098781689482907208325554680843799894826233198528393505308965377732628846163366222287698219886746543667474404243274365155048934314939391479619404400222105101714174800368808401264708068556774321622835522011480466371565912137345074785694768346361679210180644507064800027750268491674655058685693567342067058113642922455440575892572420824131469568901675894025677631135691929203337658714166023010570308963457207544037084746994016826928280848118428931484852494864487192780967627127577539702766860595249671667418348570442250719796500471495105049221477656763693866297697952211071826454973477266242570942932258279850258550978526538320760672631716430950599508780752371033310119785754733154142180842754386359177811705430982748238504564801909561029929182431823752535770975053956518769751037497088869218020518933950723853920514463419726528728696511086257149219884997874887377134568620916705849807828059751193854445009978131146915934666241071846692310107598438319191292230792503747298650929009880391941702654416816335727555703151596113564846546190897042819763365836983716328982174407366009162177850541779276367731145041782137660111010731042397832521894898817597921798666394319523936855916447118246753245630912528778330963604262982153040874560927760726641354787576616262926568298704957954913954918049209069438580790032763017941503117866862092408537949861264933479354871737451675809537088281067452440105892444976479686075120275724181874989395971643105518848195288330746699317814634930000321200327765654130472621883970596794457943468343218395304414844803701305753674262153675579814770458031413637793236291560128185336498466942261465206459942072917119370602444929358037007718981097362533224548366988505528285966192805098447175198503666680874970496982273220244823343097169111136813588418696549323714996941979687803008850408979618598756579894836445212043698216415292987811742973332588607915912510967187510929248475023930572665446276200923068791518135803477701295593646298412366497023355174586195564772461857717369368404676577047874319780573853271810933883496338813069945569399346101090745616033312247949360455361849123333063704751724871276379140924398331810164737823379692265637682071706935846394531616949411701841938119405416449466111274712819705817783293841742231409930022911502362192186723337268385688273533371925103412930705632544426611429765388301822384091026198582888433587455960453004548370789052578473166283701953392231047527564998119228742789713715713228319641003422124210082180679525276689858180956119208391760721080919923461516952599099473782780648128058792731993893453415320185969711021407542282796298237068941764740642225757212455392526179373652434440560595336591539160312524480149313234572453879524389036839236450507881731359711238145323701508413491122324390927681724749607955799151363982881058285740538000653371655553014196332241918087621018204919492651483892" +) + +var ( + ln10 = newConstApproximation(strLn10) +) + +type constApproximation struct { + exact Decimal + approximations []Decimal +} + +func newConstApproximation(value string) constApproximation { + parts := strings.Split(value, ".") + coeff, fractional := parts[0], parts[1] + + coeffLen := len(coeff) + maxPrecision := len(fractional) + + var approximations []Decimal + for p := 1; p < maxPrecision; p *= 2 { + r := RequireFromString(value[:coeffLen+p]) + approximations = append(approximations, r) + } + + return constApproximation{ + RequireFromString(value), + approximations, + } +} + +// Returns the smallest approximation available that's at least as precise +// as the passed precision (places after decimal point), i.e. Floor[ log2(precision) ] + 1 +func (c constApproximation) withPrecision(precision int32) Decimal { + i := 0 + + if precision >= 1 { + i++ + } + + for precision >= 16 { + precision /= 16 + i += 4 + } + + for precision >= 2 { + precision /= 2 + i++ + } + + if i >= len(c.approximations) { + return c.exact + } + + return c.approximations[i] +} diff --git a/const_test.go b/const_test.go new file mode 100644 index 00000000..67a27a70 --- /dev/null +++ b/const_test.go @@ -0,0 +1,34 @@ +package decimal + +import "testing" + +func TestConstApproximation(t *testing.T) { + for _, testCase := range []struct { + Const string + Precision int32 + ExpectedApproximation string + }{ + {"2.3025850929940456840179914546", 0, "2"}, + {"2.3025850929940456840179914546", 1, "2.3"}, + {"2.3025850929940456840179914546", 3, "2.302"}, + {"2.3025850929940456840179914546", 5, "2.302585"}, + {"2.3025850929940456840179914546", 10, "2.302585092994045"}, + {"2.3025850929940456840179914546", 100, "2.3025850929940456840179914546"}, + {"2.3025850929940456840179914546", -1, "2"}, + {"2.3025850929940456840179914546", -5, "2"}, + {"3.14159265359", 0, "3"}, + {"3.14159265359", 1, "3.1"}, + {"3.14159265359", 2, "3.141"}, + {"3.14159265359", 4, "3.1415926"}, + {"3.14159265359", 13, "3.14159265359"}, + } { + ca := newConstApproximation(testCase.Const) + expected, _ := NewFromString(testCase.ExpectedApproximation) + + approximation := ca.withPrecision(testCase.Precision) + + if approximation.Cmp(expected) != 0 { + t.Errorf("expected approximation %s, got %s - for const with %s precision %d", testCase.ExpectedApproximation, approximation.String(), testCase.Const, testCase.Precision) + } + } +} diff --git a/decimal.go b/decimal.go index ada12998..aa0f364a 100644 --- a/decimal.go +++ b/decimal.go @@ -808,6 +808,123 @@ func (d Decimal) ExpTaylor(precision int32) (Decimal, error) { return result, nil } +// Ln calculates natural logarithm of d. +// Precision argument specifies how precise the result must be (number of digits after decimal point). +// Negative precision is allowed. +// +// Example: +// +// d1, err := NewFromFloat(13.3).Ln(2) +// d1.String() // output: "2.59" +// +// d2, err := NewFromFloat(579.161).Ln(10) +// d2.String() // output: "6.3615805046" +func (d Decimal) Ln(precision int32) (Decimal, error) { + // Algorithm based on The Use of Iteration Methods for Approximating the Natural Logarithm, + // James F. Epperson, The American Mathematical Monthly, Vol. 96, No. 9, November 1989, pp. 831-835. + if d.IsNegative() { + return Decimal{}, fmt.Errorf("cannot calculate natural logarithm for negative decimals") + } + + if d.IsZero() { + return Decimal{}, fmt.Errorf("cannot represent natural logarithm of 0, result: -infinity") + } + + calcPrecision := precision + 2 + z := d.Copy() + + var comp1, comp3, comp2, comp4, reduceAdjust Decimal + comp1 = z.Sub(Decimal{oneInt, 0}) + comp3 = Decimal{oneInt, -1} + + // for decimal in range [0.9, 1.1] where ln(d) is close to 0 + usePowerSeries := false + + if comp1.Abs().Cmp(comp3) <= 0 { + usePowerSeries = true + } else { + // reduce input decimal to range [0.1, 1) + expDelta := int32(z.NumDigits()) + z.exp + z.exp -= expDelta + + // Input decimal was reduced by factor of 10^expDelta, thus we will need to add + // ln(10^expDelta) = expDelta * ln(10) + // to the result to compensate that + ln10 := ln10.withPrecision(calcPrecision) + reduceAdjust = NewFromInt32(expDelta) + reduceAdjust = reduceAdjust.Mul(ln10) + + comp1 = z.Sub(Decimal{oneInt, 0}) + + if comp1.Abs().Cmp(comp3) <= 0 { + usePowerSeries = true + } else { + // initial estimate using floats + zFloat := z.InexactFloat64() + comp1 = NewFromFloat(math.Log(zFloat)) + } + } + + epsilon := Decimal{oneInt, -calcPrecision} + + if usePowerSeries { + // Power Series - https://en.wikipedia.org/wiki/Logarithm#Power_series + // Calculating n-th term of formula: ln(z+1) = 2 sum [ 1 / (2n+1) * (z / (z+2))^(2n+1) ] + // until the difference between current and next term is smaller than epsilon. + // Coverage quite fast for decimals close to 1.0 + + // z + 2 + comp2 = comp1.Add(Decimal{twoInt, 0}) + // z / (z + 2) + comp3 = comp1.DivRound(comp2, calcPrecision) + // 2 * (z / (z + 2)) + comp1 = comp3.Add(comp3) + comp2 = comp1.Copy() + + for n := 1; ; n++ { + // 2 * (z / (z+2))^(2n+1) + comp2 = comp2.Mul(comp3).Mul(comp3) + + // 1 / (2n+1) * 2 * (z / (z+2))^(2n+1) + comp4 = NewFromInt(int64(2*n + 1)) + comp4 = comp2.DivRound(comp4, calcPrecision) + + // comp1 = 2 sum [ 1 / (2n+1) * (z / (z+2))^(2n+1) ] + comp1 = comp1.Add(comp4) + + if comp4.Abs().Cmp(epsilon) <= 0 { + break + } + } + } else { + // Halley's Iteration. + // Calculating n-th term of formula: a_(n+1) = a_n - 2 * (exp(a_n) - z) / (exp(a_n) + z), + // until the difference between current and next term is smaller than epsilon + for { + // exp(a_n) + comp3, _ = comp1.ExpTaylor(calcPrecision) + // exp(a_n) - z + comp2 = comp3.Sub(z) + // 2 * (exp(a_n) - z) + comp2 = comp2.Add(comp2) + // exp(a_n) + z + comp4 = comp3.Add(z) + // 2 * (exp(a_n) - z) / (exp(a_n) + z) + comp3 = comp2.DivRound(comp4, calcPrecision) + // comp1 = a_(n+1) = a_n - 2 * (exp(a_n) - z) / (exp(a_n) + z) + comp1 = comp1.Sub(comp3) + + if comp3.Abs().Cmp(epsilon) <= 0 { + break + } + } + } + + comp1 = comp1.Add(reduceAdjust) + + return comp1.Round(precision), 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 { diff --git a/decimal_test.go b/decimal_test.go index c6d2b249..2df95915 100644 --- a/decimal_test.go +++ b/decimal_test.go @@ -2749,6 +2749,67 @@ func TestDecimal_ExpTaylor(t *testing.T) { } } +func TestDecimal_Ln(t *testing.T) { + for _, testCase := range []struct { + Dec string + Precision int32 + Expected string + }{ + {"0.1", 25, "-2.3025850929940456840179915"}, + {"0.01", 25, "-4.6051701859880913680359829"}, + {"0.001", 25, "-6.9077552789821370520539744"}, + {"0.00000001", 25, "-18.4206807439523654721439316"}, + {"1.0", 10, "0.0"}, + {"1.01", 25, "0.0099503308531680828482154"}, + {"1.001", 25, "0.0009995003330835331668094"}, + {"1.0001", 25, "0.0000999950003333083353332"}, + {"1.1", 25, "0.0953101798043248600439521"}, + {"1.13", 25, "0.1222176327242492005461486"}, + {"3.13", 10, "1.1410330046"}, + {"3.13", 25, "1.1410330045520618486427824"}, + {"3.13", 50, "1.14103300455206184864278239988848193837089629107972"}, + {"3.13", 100, "1.1410330045520618486427823998884819383708962910797239760817078430268177216960996098918971117211892839"}, + {"5.71", 25, "1.7422190236679188486939833"}, + {"5.7185108151957193571930205", 50, "1.74370842450178929149992165925283704012576949094645"}, + {"839101.0351", 25, "13.6400864014410013994397240"}, + {"839101.0351094726488848490572028502", 50, "13.64008640145229044389152437468283605382056561604272"}, + {"5023583755703750094849.03519358513093500275017501750602739169823", 25, "49.9684305274348922267409953"}, + {"5023583755703750094849.03519358513093500275017501750602739169823", -1, "50.0"}, + } { + d, _ := NewFromString(testCase.Dec) + expected, _ := NewFromString(testCase.Expected) + + ln, err := d.Ln(testCase.Precision) + if err != nil { + t.Fatal(err) + } + + if ln.Cmp(expected) != 0 { + t.Errorf("expected %s, got %s, for decimal %s", testCase.Expected, ln.String(), testCase.Dec) + } + } +} + +func TestDecimal_LnZero(t *testing.T) { + d := New(0, 0) + + _, err := d.Ln(5) + + if err == nil { + t.Errorf("expected error, natural logarithm of 0 cannot be represented (-infinity)") + } +} + +func TestDecimal_LnNegative(t *testing.T) { + d := New(-20, 2) + + _, err := d.Ln(5) + + if err == nil { + t.Errorf("expected error, natural logarithm cannot be calculated for nagative decimals") + } +} + func TestDecimal_NumDigits(t *testing.T) { for _, testCase := range []struct { Dec string