feat(scripts): polynomial and rational approximations (backport #3620) (

#3685)

Co-authored-by: Roman <[email protected]>

.gitignore

@@ -233,4 +233,8 @@ blocks.db
 
 # Ignore e2e test artifacts (which clould leak information if commited)
 .ash_history
-.bash_history
+.bash_history
+
+# Python
+**/venv/*
+**/*.pyc

CHANGELOG.md

@@ -44,8 +44,9 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ### Misc Improvements
 
-* [#2788](https://github.com/osmosis-labs/osmosis/pull/2788) Add logarithm base 2 implementation.
 * [#3677](https://github.com/osmosis-labs/osmosis/pull/3677) Add methods for cloning and mutative multiplication on osmomath.BigDec.
+* [#3676](https://github.com/osmosis-labs/osmosis/pull/3676) implement `PowerInteger` function on `osmomath.BigDec` 
+* [#3678](https://github.com/osmosis-labs/osmosis/pull/3678) implement mutative `PowerIntegerMut` function on `osmomath.BigDec`.
 
 ### Features
 

osmomath/decimal.go

@@ -989,3 +989,31 @@ func (x BigDec) CustomBaseLog(base BigDec) BigDec {
 
 	return y
 }
+
+// PowerInteger takes a given decimal to an integer power
+// and returns the result. Non-mutative. Uses square and multiply
+// algorithm for performing the calculation.
+func (d BigDec) PowerInteger(power uint64) BigDec {
+	clone := d.Clone()
+	return clone.PowerIntegerMut(power)
+}
+
+// PowerIntegerMut takes a given decimal to an integer power
+// and returns the result. Mutative. Uses square and multiply
+// algorithm for performing the calculation.
+func (d BigDec) PowerIntegerMut(power uint64) BigDec {
+	if power == 0 {
+		return OneDec()
+	}
+	tmp := OneDec()
+
+	for i := power; i > 1; {
+		if i%2 != 0 {
+			tmp = tmp.MulMut(d)
+		}
+		i /= 2
+		d = d.MulMut(d)
+	}
+
+	return d.MulMut(tmp)
+}

osmomath/decimal_test.go

@@ -24,6 +24,20 @@ func TestDecimalTestSuite(t *testing.T) {
 	suite.Run(t, new(decimalTestSuite))
 }
 
+// assertMutResult given expected value after applying a math operation, a start value,
+// mutative and non mutative results with start values, asserts that mutation are only applied
+// to the mutative versions. Also, asserts that both results match the expected value.
+func (s *decimalTestSuite) assertMutResult(expectedResult, startValue, mutativeResult, nonMutativeResult, mutativeStartValue, nonMutativeStartValue BigDec) {
+	// assert both results are as expected.
+	s.Require().Equal(expectedResult, mutativeResult)
+	s.Require().Equal(expectedResult, nonMutativeResult)
+
+	// assert that mutative method mutated the receiver
+	s.Require().Equal(mutativeStartValue, expectedResult)
+	// assert that non-mutative method did not mutate the receiver
+	s.Require().Equal(nonMutativeStartValue, startValue)
+}
+
 func TestDecApproxEq(t *testing.T) {
 	// d1 = 0.55, d2 = 0.6, tol = 0.1
 	d1 := NewDecWithPrec(55, 2)
@@ -1031,6 +1045,139 @@ func (s *decimalTestSuite) TestCustomBaseLog() {
 	}
 }
 
+func (s *decimalTestSuite) TestPowerInteger() {
+	var expectedErrTolerance = MustNewDecFromStr("0.000000000000000000000000000000100000")
+
+	tests := map[string]struct {
+		base           BigDec
+		exponent       uint64
+		expectedResult BigDec
+
+		expectedToleranceOverwrite BigDec
+	}{
+		"0^2": {
+			base:     ZeroDec(),
+			exponent: 2,
+
+			expectedResult: ZeroDec(),
+		},
+		"1^2": {
+			base:     OneDec(),
+			exponent: 2,
+
+			expectedResult: OneDec(),
+		},
+		"4^4": {
+			base:     MustNewDecFromStr("4"),
+			exponent: 4,
+
+			expectedResult: MustNewDecFromStr("256"),
+		},
+		"5^3": {
+			base:     MustNewDecFromStr("5"),
+			exponent: 4,
+
+			expectedResult: MustNewDecFromStr("625"),
+		},
+		"e^10": {
+			base:     eulersNumber,
+			exponent: 10,
+
+			// https://www.wolframalpha.com/input?i=e%5E10+41+digits
+			expectedResult: MustNewDecFromStr("22026.465794806716516957900645284244366354"),
+		},
+		"geom twap overflow: 2^log_2{max spot price + 1}": {
+			base: twoBigDec,
+			// add 1 for simplicity of calculation to isolate overflow.
+			exponent: uint64(BigDecFromSDKDec(gammtypes.MaxSpotPrice).Add(OneDec()).LogBase2().TruncateInt().Uint64()),
+
+			// https://www.wolframalpha.com/input?i=2%5E%28floor%28+log+base+2+%282%5E128%29%29%29+++39+digits
+			expectedResult: MustNewDecFromStr("340282366920938463463374607431768211456"),
+		},
+		"geom twap overflow: 2^log_2{max spot price}": {
+			base:     twoBigDec,
+			exponent: uint64(BigDecFromSDKDec(gammtypes.MaxSpotPrice).LogBase2().TruncateInt().Uint64()),
+
+			// https://www.wolframalpha.com/input?i=2%5E%28floor%28+log+base+2+%282%5E128+-+1%29%29%29+++39+digits
+			expectedResult: MustNewDecFromStr("170141183460469231731687303715884105728"),
+		},
+		"geom twap overflow: 2^log_2{max spot price / 2 - 2017}": { // 2017 is prime.
+			base:     twoBigDec,
+			exponent: uint64(BigDecFromSDKDec(gammtypes.MaxSpotPrice.Quo(sdk.NewDec(2)).Sub(sdk.NewDec(2017))).LogBase2().TruncateInt().Uint64()),
+
+			// https://www.wolframalpha.com/input?i=e%5E10+41+digits
+			expectedResult: MustNewDecFromStr("85070591730234615865843651857942052864"),
+		},
+
+		// sdk.Dec test vectors copied from osmosis-labs/cosmos-sdk:
+
+		"1.0 ^ (10) => 1.0": {
+			base:     OneDec(),
+			exponent: 10,
+
+			expectedResult: OneDec(),
+		},
+		"0.5 ^ 2 => 0.25": {
+			base:     NewDecWithPrec(5, 1),
+			exponent: 2,
+
+			expectedResult: NewDecWithPrec(25, 2),
+		},
+		"0.2 ^ 2 => 0.04": {
+			base:     NewDecWithPrec(2, 1),
+			exponent: 2,
+
+			expectedResult: NewDecWithPrec(4, 2),
+		},
+		"3 ^ 3 => 27": {
+			base:     NewBigDec(3),
+			exponent: 3,
+
+			expectedResult: NewBigDec(27),
+		},
+		"-3 ^ 4 = 81": {
+			base:     NewBigDec(-3),
+			exponent: 4,
+
+			expectedResult: NewBigDec(81),
+		},
+		"-3 ^ 50 = 717897987691852588770249": {
+			base:     NewBigDec(-3),
+			exponent: 50,
+
+			expectedResult: MustNewDecFromStr("717897987691852588770249"),
+		},
+		"-3 ^ 51 = -2153693963075557766310747": {
+			base:     NewBigDec(-3),
+			exponent: 51,
+
+			expectedResult: MustNewDecFromStr("-2153693963075557766310747"),
+		},
+		"1.414213562373095049 ^ 2 = 2": {
+			base:     NewDecWithPrec(1414213562373095049, 18),
+			exponent: 2,
+
+			expectedResult:             NewBigDec(2),
+			expectedToleranceOverwrite: MustNewDecFromStr("0.0000000000000000006"),
+		},
+	}
+
+	for name, tc := range tests {
+		tc := tc
+		s.Run(name, func() {
+
+			tolerance := expectedErrTolerance
+			if !tc.expectedToleranceOverwrite.IsNil() {
+				tolerance = tc.expectedToleranceOverwrite
+			}
+
+			actualResult := tc.base.PowerInteger(tc.exponent)
+			require.True(DecApproxEq(s.T(), tc.expectedResult, actualResult, tolerance))
+
+		})
+	}
+}
+
 func (s *decimalTestSuite) TestClone() {
 
 	// The value to change the underlying copy's
@@ -1099,16 +1246,51 @@ func (s *decimalTestSuite) TestMul_Mutation() {
 			startNonMut := tc.startValue.Clone()
 
 			resultMut := startMut.MulMut(mulBy)
-			result := startNonMut.Mul(mulBy)
+			resultNonMut := startNonMut.Mul(mulBy)
+
+			s.assertMutResult(tc.expectedMulResult, tc.startValue, resultMut, resultNonMut, startMut, startNonMut)
+		})
+	}
+}
+
+// TestMul_Mutation tests that PowerIntegerMut mutates the receiver
+// while PowerInteger is not.
+func (s *decimalTestSuite) TestPowerInteger_Mutation() {
+
+	exponent := uint64(2)
+
+	tests := map[string]struct {
+		startValue     BigDec
+		expectedResult BigDec
+	}{
+		"1": {
+			startValue:     OneDec(),
+			expectedResult: OneDec(),
+		},
+		"-3": {
+			startValue:     MustNewDecFromStr("-3"),
+			expectedResult: MustNewDecFromStr("9"),
+		},
+		"0": {
+			startValue:     ZeroDec(),
+			expectedResult: ZeroDec(),
+		},
+		"4": {
+			startValue:     MustNewDecFromStr("4.5"),
+			expectedResult: MustNewDecFromStr("20.25"),
+		},
+	}
+
+	for name, tc := range tests {
+		s.Run(name, func() {
+
+			startMut := tc.startValue.Clone()
+			startNonMut := tc.startValue.Clone()
 
-			// assert both results are as expectde.
-			s.Require().Equal(tc.expectedMulResult, resultMut)
-			s.Require().Equal(tc.expectedMulResult, result)
+			resultMut := startMut.PowerIntegerMut(exponent)
+			resultNonMut := startNonMut.PowerInteger(exponent)
 
-			// assert MulMut mutated the receiver
-			s.Require().Equal(tc.expectedMulResult, startMut)
-			// assert Mul did not mutate the receiver
-			s.Require().Equal(tc.startValue, startNonMut)
+			s.assertMutResult(tc.expectedResult, tc.startValue, resultMut, resultNonMut, startMut, startNonMut)
 		})
 	}
 }

osmomath/math.go

@@ -18,6 +18,10 @@ var (
 	one_half sdk.Dec = sdk.MustNewDecFromStr("0.5")
 	one      sdk.Dec = sdk.OneDec()
 	two      sdk.Dec = sdk.MustNewDecFromStr("2")
+
+	// https://www.wolframalpha.com/input?i=2.718281828459045235360287471352662498&assumption=%22ClashPrefs%22+-%3E+%7B%22Math%22%7D
+	// nolint: unused
+	eulersNumber = MustNewDecFromStr("2.718281828459045235360287471352662498")
 )
 
 // Returns the internal "power precision".

scripts/approximations/README.md

@@ -0,0 +1,65 @@
+# Math Operations Approximations
+
+## Context
+
+This is a set of scripts to approximate a mathematical operation using polynomial
+and rational approximations.
+
+The `main` script is in its respective function of `main.py`. It does the following:
+
+1. Computes polynomial and rational approximations of a given function (e^x by default), 
+returning the coefficients.
+
+1. Computes (x,y) coordinates for every kind of approximation given the same x coordinates.
+Plots the results for rough comparison.
+
+1. Plots the results for rough comparison.
+
+2. Computes the max error for every approximation given the same x coordinates.
+
+3. Computes and plots max errors for every approximation with a varying number of parameters.
+
+In other words, this script runs various approximation methods, plots their results and deltas
+from actual function values. It can be configured to print the maximum error.
+The exact behavior is controlled by the global variables at the top of `main.py`.
+
+The following are the resources used to create this:
+- <https://xn--2-umb.com/22/approximation>
+- <https://sites.tufts.edu/atasissa/files/2019/09/remez.pdf>
+
+In `main.py`, there is also an `exponent_approximation_choice` script.
+
+This is a shorter and simpler version of `main` that isolates the 13-parameter
+Chebyshev Rational approximation of e^x. We are planning to use it in production.
+Therefore, we need to perform coefficient truncations to 36 decimal points
+(the max `osmomath` supported precision). This truncation is applied
+to `exponent_approximation_choice` but not `main`.
+
+## Configuration
+
+Several parameters can be changed on the needs basis at the
+top of `main.py`.
+
+Some of the parameters include the function we are approximating, the [x_start, x_end] range of
+the approximation, and the number of terms to be used. For the full parameter list, see `main.py`.
+
+## Usage
+
+Assuming that you are in the root of the repository and have Sympy installed:
+
+```bash
+# Create a virtual environment.
+python3 -m venv ~/approx-venv
+
+# Start the environment
+source ~/approx-venv/bin/activate
+
+# Install dependencies in the virtual environment.
+pip install -r scripts/approximations/requirements.txt
+
+# Run the script.
+python3 scripts/approximations/main.py
+
+# Run tests
+python3 scripts/approximations/rational_test.py
+```