feat(osmomath): log2 approximation (#2788)

* feat: osmomath log2 approximation

* lint

* fix comment

* comment

* Update osmomath/decimal.go

* improve accuracy with narrower range

* comment

* implement and test precise log for x >= 1

* tests for negative values

* make non-mutative and test

* bench

* changelog

* remove redundant assignments

* improve comments

CHANGELOG.md

@@ -44,6 +44,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ### Features
 
+* [#2788](https://github.com/osmosis-labs/osmosis/pull/2788) Add logarithm base 2 implementation.
 * [#2739](https://github.com/osmosis-labs/osmosis/pull/2739) Add pool type query
 
 ### Bug fixes

osmomath/decimal.go

@@ -30,6 +30,9 @@ const (
 
 	// max number of iterations in ApproxRoot function
 	maxApproxRootIterations = 100
+
+	// max number of iterations in Log2 function
+	maxLog2Iterations = 300
 )
 
 var (
@@ -39,6 +42,10 @@ var (
 	zeroInt              = big.NewInt(0)
 	oneInt               = big.NewInt(1)
 	tenInt               = big.NewInt(10)
+
+	// initialized in init() since requires
+	// precision to be defined.
+	twoBigDec BigDec
 )
 
 // Decimal errors
@@ -54,6 +61,8 @@ func init() {
 	for i := 0; i <= Precision; i++ {
 		precisionMultipliers[i] = calcPrecisionMultiplier(int64(i))
 	}
+
+	twoBigDec = NewBigDec(2)
 }
 
 func precisionInt() *big.Int {
@@ -211,7 +220,8 @@ func (d BigDec) GTE(d2 BigDec) bool   { return (d.i).Cmp(d2.i) >= 0 }          /
 func (d BigDec) LT(d2 BigDec) bool    { return (d.i).Cmp(d2.i) < 0 }           // less than
 func (d BigDec) LTE(d2 BigDec) bool   { return (d.i).Cmp(d2.i) <= 0 }          // less than or equal
 func (d BigDec) Neg() BigDec          { return BigDec{new(big.Int).Neg(d.i)} } // reverse the decimal sign
-func (d BigDec) Abs() BigDec          { return BigDec{new(big.Int).Abs(d.i)} } // absolute value
+// nolint: stylecheck
+func (d BigDec) Abs() BigDec { return BigDec{new(big.Int).Abs(d.i)} } // absolute value
 
 // BigInt returns a copy of the underlying big.Int.
 func (d BigDec) BigInt() *big.Int {
@@ -856,3 +866,54 @@ func DecApproxEq(t *testing.T, d1 BigDec, d2 BigDec, tol BigDec) (*testing.T, bo
 	diff := d1.Sub(d2).Abs()
 	return t, diff.LTE(tol), "expected |d1 - d2| <:\t%v\ngot |d1 - d2| = \t\t%v", tol.String(), diff.String()
 }
+
+// LogBase2 returns log_2 {x}.
+// Rounds down by truncations during division and right shifting.
+// Accurate up to 32 precision digits.
+// Implementation is based on:
+// https://stm32duinoforum.com/forum/dsp/BinaryLogarithm.pdf
+func (x BigDec) LogBase2() BigDec {
+	// create a new decimal to avoid mutating
+	// the receiver's int buffer.
+	xCopy := ZeroDec()
+	xCopy.i = new(big.Int).Set(x.i)
+	if xCopy.LTE(ZeroDec()) {
+		panic(fmt.Sprintf("log is not defined at <= 0, given (%s)", xCopy))
+	}
+
+	// Normalize x to be 1 <= x < 2.
+
+	// y is the exponent that results in a whole multiple of 2.
+	y := ZeroDec()
+
+	// repeat until: x >= 1.
+	for xCopy.LT(OneDec()) {
+		xCopy.i.Lsh(xCopy.i, 1)
+		y = y.Sub(OneDec())
+	}
+
+	// repeat until: x < 2.
+	for xCopy.GTE(twoBigDec) {
+		xCopy.i.Rsh(xCopy.i, 1)
+		y = y.Add(OneDec())
+	}
+
+	b := OneDec().Quo(twoBigDec)
+
+	// N.B. At this point x is a positive real number representing
+	// mantissa of the log. We estimate it using the following
+	// algorithm:
+	// https://stm32duinoforum.com/forum/dsp/BinaryLogarithm.pdf
+	// This has shown precision of 32 digits relative
+	// to Wolfram Alpha in tests.
+	for i := 0; i < maxLog2Iterations; i++ {
+		xCopy = xCopy.Mul(xCopy)
+		if xCopy.GTE(twoBigDec) {
+			xCopy.i.Rsh(xCopy.i, 1)
+			y = y.Add(b)
+		}
+		b.i.Rsh(b.i, 1)
+	}
+
+	return y
+}

osmomath/decimal_test.go

@@ -11,6 +11,8 @@ import (
 	"github.com/stretchr/testify/require"
 	"github.com/stretchr/testify/suite"
 	"gopkg.in/yaml.v2"
+
+	"github.com/osmosis-labs/osmosis/v12/app/apptesting/osmoassert"
 )
 
 type decimalTestSuite struct {
@@ -152,8 +154,8 @@ func (s *decimalTestSuite) TestDecFloat64() {
 
 func (s *decimalTestSuite) TestSdkDec() {
 	tests := []struct {
-		d    BigDec
-		want sdk.Dec
+		d        BigDec
+		want     sdk.Dec
 		expPanic bool
 	}{
 		{NewBigDec(0), sdk.MustNewDecFromStr("0.000000000000000000"), false},
@@ -177,8 +179,8 @@ func (s *decimalTestSuite) TestSdkDec() {
 
 func (s *decimalTestSuite) TestBigDecFromSdkDec() {
 	tests := []struct {
-		d    sdk.Dec
-		want BigDec
+		d        sdk.Dec
+		want     BigDec
 		expPanic bool
 	}{
 		{sdk.MustNewDecFromStr("0.000000000000000000"), NewBigDec(0), false},
@@ -202,8 +204,8 @@ func (s *decimalTestSuite) TestBigDecFromSdkDec() {
 
 func (s *decimalTestSuite) TestBigDecFromSdkDecSlice() {
 	tests := []struct {
-		d    []sdk.Dec
-		want []BigDec
+		d        []sdk.Dec
+		want     []BigDec
 		expPanic bool
 	}{
 		{[]sdk.Dec{sdk.MustNewDecFromStr("0.000000000000000000")}, []BigDec{NewBigDec(0)}, false},
@@ -441,14 +443,14 @@ func (s *decimalTestSuite) TestDecCeil() {
 		input    BigDec
 		expected BigDec
 	}{
-		{MustNewDecFromStr("0.001"), NewBigDec(1)}, // 0.001 => 1.0
-		{MustNewDecFromStr("-0.001"), ZeroDec()},   // -0.001 => 0.0
-		{ZeroDec(), ZeroDec()}, // 0.0 => 0.0
-		{MustNewDecFromStr("0.9"), NewBigDec(1)},    // 0.9 => 1.0
+		{MustNewDecFromStr("0.001"), NewBigDec(1)},   // 0.001 => 1.0
+		{MustNewDecFromStr("-0.001"), ZeroDec()},     // -0.001 => 0.0
+		{ZeroDec(), ZeroDec()},                       // 0.0 => 0.0
+		{MustNewDecFromStr("0.9"), NewBigDec(1)},     // 0.9 => 1.0
 		{MustNewDecFromStr("4.001"), NewBigDec(5)},   // 4.001 => 5.0
 		{MustNewDecFromStr("-4.001"), NewBigDec(-4)}, // -4.001 => -4.0
-		{MustNewDecFromStr("4.7"), NewBigDec(5)},   // 4.7 => 5.0
-		{MustNewDecFromStr("-4.7"), NewBigDec(-4)}, // -4.7 => -4.0
+		{MustNewDecFromStr("4.7"), NewBigDec(5)},     // 4.7 => 5.0
+		{MustNewDecFromStr("-4.7"), NewBigDec(-4)},   // -4.7 => -4.0
 	}
 
 	for i, tc := range testCases {
@@ -463,11 +465,11 @@ func (s *decimalTestSuite) TestPower() {
 		power    uint64
 		expected BigDec
 	}{
-		{OneDec(), 10, OneDec()},                                               // 1.0 ^ (10) => 1.0
-		{NewDecWithPrec(5, 1), 2, NewDecWithPrec(25, 2)},                       // 0.5 ^ 2 => 0.25
-		{NewDecWithPrec(2, 1), 2, NewDecWithPrec(4, 2)},                        // 0.2 ^ 2 => 0.04
-		{NewDecFromInt(NewInt(3)), 3, NewDecFromInt(NewInt(27))},               // 3 ^ 3 => 27
-		{NewDecFromInt(NewInt(-3)), 4, NewDecFromInt(NewInt(81))},              // -3 ^ 4 = 81
+		{OneDec(), 10, OneDec()},                                                                   // 1.0 ^ (10) => 1.0
+		{NewDecWithPrec(5, 1), 2, NewDecWithPrec(25, 2)},                                           // 0.5 ^ 2 => 0.25
+		{NewDecWithPrec(2, 1), 2, NewDecWithPrec(4, 2)},                                            // 0.2 ^ 2 => 0.04
+		{NewDecFromInt(NewInt(3)), 3, NewDecFromInt(NewInt(27))},                                   // 3 ^ 3 => 27
+		{NewDecFromInt(NewInt(-3)), 4, NewDecFromInt(NewInt(81))},                                  // -3 ^ 4 = 81
 		{MustNewDecFromStr("1.414213562373095048801688724209698079"), 2, NewDecFromInt(NewInt(2))}, // 1.414213562373095048801688724209698079 ^ 2 = 2
 	}
 
@@ -483,14 +485,14 @@ func (s *decimalTestSuite) TestApproxRoot() {
 		root     uint64
 		expected BigDec
 	}{
-		{OneDec(), 10, OneDec()},                                                       // 1.0 ^ (0.1) => 1.0
-		{NewDecWithPrec(25, 2), 2, NewDecWithPrec(5, 1)},                               // 0.25 ^ (0.5) => 0.5
-		{NewDecWithPrec(4, 2), 2, NewDecWithPrec(2, 1)},                                // 0.04 ^ (0.5) => 0.2
-		{NewDecFromInt(NewInt(27)), 3, NewDecFromInt(NewInt(3))},                       // 27 ^ (1/3) => 3
-		{NewDecFromInt(NewInt(-81)), 4, NewDecFromInt(NewInt(-3))},                     // -81 ^ (0.25) => -3
-		{NewDecFromInt(NewInt(2)), 2, MustNewDecFromStr("1.414213562373095048801688724209698079")},         // 2 ^ (0.5) => 1.414213562373095048801688724209698079
+		{OneDec(), 10, OneDec()},                                                                         // 1.0 ^ (0.1) => 1.0
+		{NewDecWithPrec(25, 2), 2, NewDecWithPrec(5, 1)},                                                 // 0.25 ^ (0.5) => 0.5
+		{NewDecWithPrec(4, 2), 2, NewDecWithPrec(2, 1)},                                                  // 0.04 ^ (0.5) => 0.2
+		{NewDecFromInt(NewInt(27)), 3, NewDecFromInt(NewInt(3))},                                         // 27 ^ (1/3) => 3
+		{NewDecFromInt(NewInt(-81)), 4, NewDecFromInt(NewInt(-3))},                                       // -81 ^ (0.25) => -3
+		{NewDecFromInt(NewInt(2)), 2, MustNewDecFromStr("1.414213562373095048801688724209698079")},       // 2 ^ (0.5) => 1.414213562373095048801688724209698079
 		{NewDecWithPrec(1005, 3), 31536000, MustNewDecFromStr("1.000000000158153903837946258002096839")}, // 1.005 ^ (1/31536000) ≈ 1.000000000158153903837946258002096839
-		{SmallestDec(), 2, NewDecWithPrec(1, 18)},                                       // 1e-36 ^ (0.5) => 1e-18
+		{SmallestDec(), 2, NewDecWithPrec(1, 18)},                                                        // 1e-36 ^ (0.5) => 1e-18
 		{SmallestDec(), 3, MustNewDecFromStr("0.000000000001000000000000000002431786")},                  // 1e-36 ^ (1/3) => 1e-12
 		{NewDecWithPrec(1, 8), 3, MustNewDecFromStr("0.002154434690031883721759293566519280")},           // 1e-8 ^ (1/3) ≈ 0.002154434690031883721759293566519
 	}
@@ -511,11 +513,11 @@ func (s *decimalTestSuite) TestApproxSqrt() {
 		input    BigDec
 		expected BigDec
 	}{
-		{OneDec(), OneDec()},                                                // 1.0 => 1.0
-		{NewDecWithPrec(25, 2), NewDecWithPrec(5, 1)},                       // 0.25 => 0.5
-		{NewDecWithPrec(4, 2), NewDecWithPrec(2, 1)},                        // 0.09 => 0.3
-		{NewDecFromInt(NewInt(9)), NewDecFromInt(NewInt(3))},                // 9 => 3
-		{NewDecFromInt(NewInt(-9)), NewDecFromInt(NewInt(-3))},              // -9 => -3
+		{OneDec(), OneDec()},                                                                    // 1.0 => 1.0
+		{NewDecWithPrec(25, 2), NewDecWithPrec(5, 1)},                                           // 0.25 => 0.5
+		{NewDecWithPrec(4, 2), NewDecWithPrec(2, 1)},                                            // 0.09 => 0.3
+		{NewDecFromInt(NewInt(9)), NewDecFromInt(NewInt(3))},                                    // 9 => 3
+		{NewDecFromInt(NewInt(-9)), NewDecFromInt(NewInt(-3))},                                  // -9 => -3
 		{NewDecFromInt(NewInt(2)), MustNewDecFromStr("1.414213562373095048801688724209698079")}, // 2 => 1.414213562373095048801688724209698079
 	}
 
@@ -650,3 +652,90 @@ func BenchmarkMarshalTo(b *testing.B) {
 		}
 	}
 }
+
+func (s *decimalTestSuite) TestLog2() {
+	var expectedErrTolerance = MustNewDecFromStr("0.000000000000000000000000000000000100")
+
+	tests := map[string]struct {
+		initialValue BigDec
+		expected     BigDec
+
+		expectedPanic bool
+	}{
+		"log_2{-1}; invalid; panic": {
+			initialValue:  OneDec().Neg(),
+			expectedPanic: true,
+		},
+		"log_2{0}; invalid; panic": {
+			initialValue:  ZeroDec(),
+			expectedPanic: true,
+		},
+		"log_2{0.001} = -9.965784284662087043610958288468170528": {
+			initialValue: MustNewDecFromStr("0.001"),
+			// From: https://www.wolframalpha.com/input?i=log+base+2+of+0.999912345+with+33+digits
+			expected: MustNewDecFromStr("-9.965784284662087043610958288468170528"),
+		},
+		"log_2{0.56171821941421412902170941} = -0.832081497183140708984033250637831402": {
+			initialValue: MustNewDecFromStr("0.56171821941421412902170941"),
+			// From: https://www.wolframalpha.com/input?i=log+base+2+of+0.56171821941421412902170941+with+36+digits
+			expected: MustNewDecFromStr("-0.832081497183140708984033250637831402"),
+		},
+		"log_2{0.999912345} = -0.000126464976533858080645902722235833": {
+			initialValue: MustNewDecFromStr("0.999912345"),
+			// From: https://www.wolframalpha.com/input?i=log+base+2+of+0.999912345+with+37+digits
+			expected: MustNewDecFromStr("-0.000126464976533858080645902722235833"),
+		},
+		"log_2{1} = 0": {
+			initialValue: NewBigDec(1),
+			expected:     NewBigDec(0),
+		},
+		"log_2{2} = 1": {
+			initialValue: NewBigDec(2),
+			expected:     NewBigDec(1),
+		},
+		"log_2{7} = 2.807354922057604107441969317231830809": {
+			initialValue: NewBigDec(7),
+			// From: https://www.wolframalpha.com/input?i=log+base+2+of+7+37+digits
+			expected: MustNewDecFromStr("2.807354922057604107441969317231830809"),
+		},
+		"log_2{512} = 9": {
+			initialValue: NewBigDec(512),
+			expected:     NewBigDec(9),
+		},
+		"log_2{580} = 9.179909090014934468590092754117374938": {
+			initialValue: NewBigDec(580),
+			// From: https://www.wolframalpha.com/input?i=log+base+2+of+600+37+digits
+			expected: MustNewDecFromStr("9.179909090014934468590092754117374938"),
+		},
+		"log_2{1024} = 10": {
+			initialValue: NewBigDec(1024),
+			expected:     NewBigDec(10),
+		},
+		"log_2{1024.987654321} = 10.001390817654141324352719749259888355": {
+			initialValue: NewDecWithPrec(1024987654321, 9),
+			// From: https://www.wolframalpha.com/input?i=log+base+2+of+1024.987654321+38+digits
+			expected: MustNewDecFromStr("10.001390817654141324352719749259888355"),
+		},
+		"log_2{912648174127941279170121098210.92821920190204131121} = 99.525973560175362367047484597337715868": {
+			initialValue: MustNewDecFromStr("912648174127941279170121098210.92821920190204131121"),
+			// From: https://www.wolframalpha.com/input?i=log+base+2+of+912648174127941279170121098210.92821920190204131121+38+digits
+			expected: MustNewDecFromStr("99.525973560175362367047484597337715868"),
+		},
+	}
+
+	for name, tc := range tests {
+		s.Run(name, func() {
+			osmoassert.ConditionalPanic(s.T(), tc.expectedPanic, func() {
+				// Create a copy to test that the original was not modified.
+				// That is, that LogbBase2() is non-mutative.
+				initialCopy := ZeroDec()
+				initialCopy.i.Set(tc.initialValue.i)
+
+				// system under test.
+				res := tc.initialValue.LogBase2()
+				require.True(DecApproxEq(s.T(), tc.expected, res, expectedErrTolerance))
+				require.Equal(s.T(), initialCopy, tc.initialValue)
+			})
+		})
+	}
+}

osmomath/log2_bench_test.go

+package osmomath
+
+import (
+	"math/rand"
+	"testing"
+)
+
+func BenchmarkLog2(b *testing.B) {
+	tests := []BigDec{
+		MustNewDecFromStr("1.2"),
+		MustNewDecFromStr("1.234"),
+		MustNewDecFromStr("1024"),
+		NewBigDec(2048 * 2048 * 2048 * 2048 * 2048),
+		MustNewDecFromStr("999999999999999999999999999999999999999999999999999999.9122181273612911"),
+		MustNewDecFromStr("0.563289239121902491248219047129047129"),
+	}
+
+	for i := 0; i < b.N; i++ {
+		b.StopTimer()
+		test := tests[rand.Int63n(int64(len(tests)))]
+		b.StartTimer()
+		_ = test.LogBase2()
+	}
+}

osmomath/rounding_direction_test.go

@@ -42,7 +42,6 @@ func TestDivIntByU64ToBigDec(t *testing.T) {
 	addTCForAllRoundingModes("odd divided by 2", sdk.NewInt(5), 2, NewDecWithPrec(25, 1))
 
 	for name, tt := range tests {
-		fmt.Println("start")
 		t.Run(name, func(t *testing.T) {
 			got, err := DivIntByU64ToBigDec(tt.i, tt.u, tt.round)
 			require.Equal(t, tt.want, got)