Repository: osmosis-labs/osmosis
Pull request #3685 (2022-12-10 20:08:10)

feat(scripts): polynomial and rational approximations (backport #3620)
This is an automatic backport of pull request #3620 done by [Mergify](https://mergify.com).
Cherry-pick of 5f55d1beaf379ccd0d5d31274fdd054887c652d6 has failed:
```
On branch mergify/bp/v13.x/pr-3620
Your branch is up to date with 'origin/v13.x'.

You are currently cherry-picking commit 5f55d1be.
  (fix conflicts and run "git cherry-pick --continue")
  (use "git cherry-pick --skip" to skip this patch)
  (use "git cherry-pick --abort" to cancel the cherry-pick operation)

Changes to be committed:
	modified:   osmomath/decimal.go
	modified:   osmomath/decimal_test.go
	modified:   osmomath/math.go
	new file:   scripts/approximations/README.md
	new file:   scripts/approximations/approximations.py
	new file:   scripts/approximations/chebyshev.py
	new file:   scripts/approximations/main.py
	new file:   scripts/approximations/polynomial.py
	new file:   scripts/approximations/rational.py
	new file:   scripts/approximations/rational_test.py
	new file:   scripts/approximations/requirements.txt

Unmerged paths:
  (use "git add <file>..." to mark resolution)
	both modified:   .gitignore
	both modified:   CHANGELOG.md

```


To fix up this pull request, you can check it out locally. See documentation: https://docs.github.com/en/github/collaborating-with-pull-requests/reviewing-changes-in-pull-requests/checking-out-pull-requests-locally

---


<details>
<summary>Mergify commands and options</summary>

<br />

More conditions and actions can be found in the [documentation](https://docs.mergify.com/).

You can also trigger Mergify actions by commenting on this pull request:

- `@Mergifyio refresh` will re-evaluate the rules
- `@Mergifyio rebase` will rebase this PR on its base branch
- `@Mergifyio update` will merge the base branch into this PR
- `@Mergifyio backport <destination>` will backport this PR on `<destination>` branch

Additionally, on Mergify [dashboard](https://dashboard.mergify.com/) you can:

- look at your merge queues
- generate the Mergify configuration with the config editor.

Finally, you can contact us on https://mergify.com
</details>


Patch:
diff --git a/.gitignore b/.gitignore
index ff99e2886..188682326 100644
--- a/.gitignore
+++ b/.gitignore
@@ -233,4 +233,8 @@ blocks.db
 
 # Ignore e2e test artifacts (which clould leak information if commited)
 .ash_history
-.bash_history
\ No newline at end of file
+.bash_history
+
+# Python
+**/venv/*
+**/*.pyc
diff --git a/CHANGELOG.md b/CHANGELOG.md
index d69f82f68..191a1fbc0 100644
--- a/CHANGELOG.md
+++ b/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
 
diff --git a/osmomath/decimal.go b/osmomath/decimal.go
index 55dbd399f..cba13ca44 100644
--- a/osmomath/decimal.go
+++ b/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)
+}
diff --git a/osmomath/math.go b/osmomath/math.go
index d8f341bb6..0e07cac82 100644
--- a/osmomath/math.go
+++ b/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".
diff --git a/scripts/approximations/README.md b/scripts/approximations/README.md
new file mode 100644
index 000000000..ebadb8819
--- /dev/null
+++ b/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
+```
diff --git a/scripts/approximations/approximations.py b/scripts/approximations/approximations.py
new file mode 100644
index 000000000..198496980
--- /dev/null
+++ b/scripts/approximations/approximations.py
@@ -0,0 +1,191 @@
+import sympy as sp
+
+import rational
+import polynomial
+import chebyshev
+
+def linspace(start: sp.Float, end: sp.Float, num_points: int) -> list[sp.Float]:
+    """ Given [start, end] on the x-axis, creates the desired number of
+    equispaced points and returns them in increasing order as a list. 
+    """
+    if num_points == 1:
+        return [(end - start) / 2]
+
+    result: list[sp.Float] = []
+    for i in range(num_points):
+        cur_coord = sp.Float(start + i * (end - start) / (num_points - 1), 200)
+
+        if cur_coord is sp.nan:
+            raise ValueError("cur_coord is nan")
+
+        result.append(cur_coord)
+
+    return result
+
+def approx_and_eval_all(approximated_fn, num_parameters: int, x_coordinates) -> tuple:
+    """ Performs all supported approximations of the given function, evaluates
+    each wih the given x-coordinates.
+
+    Returns y-coordinates as a tuple in the following order:
+    - Evaluated Equispaced Polynomial
+    - Evaluated Chebyshev Polynomial
+    - Evaluated Chebyshev Rational
+    - True Y-coordinates as determined by Sympy's implementation of the function.
+    """
+    x_start = x_coordinates[0]
+    x_end = x_coordinates[-1]
+    
+    # Equispaced Polynomial Approximation
+    coefficients_equispaced_poly = equispaced_poly_approx(approximated_fn, x_start, x_end, num_parameters)
+    y_eqispaced_poly = polynomial.evaluate(x_coordinates, coefficients_equispaced_poly)
+
+    # Chebyshev Polynomial Approximation
+    coefficients_chebyshev_poly = chebyshev_poly_approx(approximated_fn, x_start, x_end, num_parameters)
+    y_chebyshev_poly = polynomial.evaluate(x_coordinates, coefficients_chebyshev_poly)
+
+    # Chebyshev Rational Approximation
+    numerator_coefficients_chebyshev_rational, denominator_coefficients_chebyshev_rational = chebyshev_rational_approx(approximated_fn, x_start, x_end, num_parameters)
+    y_chebyshev_rational = rational.evaluate(x_coordinates, numerator_coefficients_chebyshev_rational, denominator_coefficients_chebyshev_rational)
+
+    # Actual
+    y_actual = get_y_actual(approximated_fn, x_coordinates)
+
+    return (y_eqispaced_poly, y_chebyshev_poly, y_chebyshev_rational, y_actual)
+
+def get_y_actual(approximated_fn, x_coordinates) -> list[sp.Float]:
+    """ Given x coordinates and a function, returns a list of y coordinates, one for each x coordinate
+    that represent the true (x,y) coordinates of the function, as calculated by sp.
+    """
+    y_actual = []
+    for x in x_coordinates:
+        y_actual.append(approximated_fn(x))
+    return y_actual
+
+def compute_max_error(y_approximation, y_actual) -> sp.Float:
+    """ Given an approximated list of y values and actual y values, computes and returns
+    the maximum error delta between them.
+
+    CONTRACT:
+    - y_approximation and y_actual must be the same length
+    - for every i in range(len(y_approximation)), y_approximation[i] and y_actual[i] must correspond to the
+    same x coordinate
+    """
+    if len(y_approximation) != len(y_actual):
+        raise ValueError(F"y_approximation ({len(y_approximation)}) and y_actual ({len(y_actual)}) must be the same length.")
+
+    max: sp.Float = None
+    for i in range(len(y_approximation)):
+        cur_abs = sp.Abs(y_approximation[i] - y_actual[i])
+        if cur_abs is sp.nan:
+            raise ValueError(F"cur_abs is nan. y_approximation[i] ({y_approximation[i]}) and y_actual[i] ({y_actual[i]})")
+        if max is None:
+            max = cur_abs
+        else:
+            max = sp.Max(max, cur_abs)
+    return max
+
+def compute_error_range(y_approximation: list, y_actual: list) -> list[sp.Float]:
+    """ Given an approximated list of y values and actual y values, computes and returns
+    error deltas between them.
+
+    CONTRACT:
+    - y_approximation and y_actual must be the same length
+    - for every i in range(len(y_approximation)), y_approximation[i] and y_actual[i] must correspond to the
+    same x coordinate
+    """
+    result = []
+    for i in range(len(y_approximation)):
+        cur_abs = sp.Abs(y_approximation[i] - y_actual[i])
+        if cur_abs is sp.nan:
+            raise ValueError(F"cur_abs is nan. y_approximation[i] ({y_approximation[i]}) and y_actual[i] ({y_actual[i]})")
+        result.append(cur_abs)
+    return result
+
+def equispaced_poly_approx(fn, x_start: sp.Float, x_end: sp.Float, num_terms: int):
+    """ Returns the coefficients for an equispaced polynomial between x_start and x_end with num_terms terms.
+
+    The return value is a list of num_terms polynomial coefficients needed to get the returned y coordinates from returned x coordinates.
+    """
+    # Compute equispaced coordinates.
+    equispaced_nodes_x = linspace(x_start, x_end, num_terms)
+    y_nodes = sp.Matrix([fn(x) for x in equispaced_nodes_x])
+
+    # Construct a system of linear equations.
+    vandermonde_matrix = polynomial.construct_vandermonde_matrix(equispaced_nodes_x)
+
+    coef = vandermonde_matrix.solve(y_nodes)
+
+    return coef
+
+def chebyshev_poly_approx(fn, x_start: int, x_end: int, num_terms: int):
+    """ Returns the coefficients for a polynomial constructed from Chebyshev nodes between x_start and x_end with num_terms terms.
+
+    Equation for Chebyshev nodes:
+    x_i = (x_start + x_end)/2 + (x_end - x_start)/2 * cos((2*i + 1)/(2*num_terms) * pi)
+
+    The return value is a list of num_terms polynomial coefficients needed to get the returned y coordinates from returned x coordinates.
+    """
+    # Compute Chebyshev coordinates.
+    x_estimated, y_estimated = chebyshev.get_nodes(fn, x_start, x_end, num_terms)
+
+    # Construct a system of linear equations.
+    vandermonde_matrix = polynomial.construct_vandermonde_matrix(x_estimated)
+
+    # Solve the matrix to get the coefficients used in the final approximation polynomial.
+    coef = vandermonde_matrix.solve(sp.Matrix(y_estimated))
+
+    return coef
+
+def chebyshev_rational_approx(fn, x_start: int, x_end: int, num_parameters: int):
+    """ Returns a rational approximation between x_start and x_end with num_terms terms
+    using Chebyshev nodes.
+
+    Equation for Chebyshev nodes:
+    x_i = (x_start + x_end)/2 + (x_end - x_start)/2 * cos((2*i + 1)/(2*num_terms) * pi)
+
+    We approximate h(x) = p(x) / q(x)
+
+    Assume num_terms_numerator = 3.
+    Then, we have
+
+    h(x) = (p_0 + p_1 x + p_2 x^2) / (1 + q_1 x + q_2 x^2) 
+
+    The return value is a list with a 2 items where:
+    - item 1: num_terms equispaced x coordinates between x_start and x_end
+    - item 2: num_terms y coordinates for the equispaced x coordinates
+    """
+    if num_parameters % 2 == 0:
+        # if num_parameters is 6, we want (3, 4) terms
+        # assume h(x) = p(x) / q(x)
+        # we always set the first term of q(x) to 1 for ease of calculation.
+        # so, we want p(x) to have 3 terms and q(x) to have 4 terms.
+        num_terms_numerator = num_parameters // 2
+        num_terms_denominator = num_parameters // 2 + 1
+    else:
+        # if num_parameters is 5, we want (3, 3) terms
+        # assume h(x) = p(x) / q(x)
+        # we always set the first term of q(x) to 1 for ease of calculation.
+        # so, we want p(x) to have 3 terms and q(x) to have 4 terms.
+        num_terms_numerator = num_parameters // 2 + 1
+        num_terms_denominator = num_parameters // 2 + 1
+
+    # Compute Chebyshev coordinates.
+    x_chebyshev, y_chebyshev = chebyshev.get_nodes(fn, x_start, x_end, num_parameters)
+
+    # Construct a system of linear equations.
+    matrix = rational.construct_rational_eval_matrix(x_chebyshev, y_chebyshev, num_terms_numerator, num_terms_denominator)
+
+    # Solve the matrix to get the coefficients used in the final approximation polynomial.
+    # coef = np.linalg.solve(np.array(matrix), y_chebyshev)
+    coef = matrix.solve(sp.Matrix(y_chebyshev))
+
+    # first num_terms_numerator values are the numerator coefficients
+    # next num_terms_numerator - 1 values are the denominator coefficients
+    coef_numerator = coef[:num_terms_numerator]
+    coef_denominator = coef[num_terms_numerator:]
+
+    # h(x) = (p_0 + p_1 x + p_2 x^2) / (1 + q_1 x + q_2 x^2)
+    # Therefore, we insert 1 here.
+    coef_denominator = [1] + coef_denominator
+
+    return [coef_numerator, coef_denominator]
diff --git a/scripts/approximations/chebyshev.py b/scripts/approximations/chebyshev.py
new file mode 100644
index 000000000..0c2986d6f
--- /dev/null
+++ b/scripts/approximations/chebyshev.py
@@ -0,0 +1,29 @@
+import sympy as sp
+from typing import Tuple
+
+# N.B. We must evaluate x values early.
+# Otherwise, solving for the coefficients will take a long time.
+_x_eval_precision = 200
+_x_eval_chop = 1e-200
+
+def get_nodes(fn, x_start: sp.Float, x_end: sp.Float, num_terms: int) -> Tuple[list, list]:
+    """ Returns Chebyshev nodes between x_start and x_end with num_terms terms
+    and the given function fn.
+
+    Equation for Chebyshev nodes:
+    x_i = (x_start + x_end)/2 + (x_end - x_start)/2 * cos((2*i + 1)/(2*num_terms) * pi)
+    
+    The first returned value is a list of x coordinates for the Chebyshev nodes.
+    The second returned value is a list of y coordinates for the Chebyshev nodes.
+    """
+    x_estimated = []
+    y_estimated = []
+
+    for i in range (num_terms):
+        x_i = ((x_start + x_end) / 2 + (x_end - x_start) / 2 * sp.cos((2*sp.Float(i,_x_eval_precision) + 1).evalf(chop=_x_eval_chop) * sp.pi.evalf(chop=_x_eval_chop) / (2 * sp.Float(num_terms, _x_eval_precision))))
+        y_i = fn(x_i)
+
+        x_estimated.append(x_i)
+        y_estimated.append(y_i)
+
+    return x_estimated, y_estimated
diff --git a/scripts/approximations/main.py b/scripts/approximations/main.py
new file mode 100644
index 000000000..b91433855
--- /dev/null
+++ b/scripts/approximations/main.py
@@ -0,0 +1,194 @@
+import matplotlib.pyplot as plt
+import sympy as sp
+
+import approximations
+import rational
+
+##########################
+# Configuration Parameters
+
+# start of the interval to calculate the approximation on
+x_start = 0
+
+# end of the interval to calculate the approximation on
+x_end = 1
+
+# number of paramters to use for the approximations.
+num_parameters = 13
+
+# number of paramters to use for plotting error deltas.
+num_parameters_errors = 30
+
+# number of (x,y) coordinates used to plot the resulting approximation.
+num_points_plot = 100000
+
+# function to approximate
+approximated_fn = lambda x: sp.Pow(sp.E, x)
+
+# fixed point precision used in Osmosis `osmomath` package.
+osmomath_precision = 36
+
+# flag controlling whether to plot each approximation.
+# Plots if true.
+shouldPlotApproximations = True
+
+# flag controlling whether to compute max error for each approximation
+# given the equally spaced x coordinates.
+# Computes if true.
+shouldComputeErrorDelta = True
+
+# flag controlling whether to plot errors over a range.
+# Currently, does so only for Chebyshev Rational Approximation.
+# Computes if true.
+shouldPlotErrorRange  = True
+
+# flag controlling whether to plot max error for every approximation
+# with a varying number of parameters. This is useful to find the most
+# optimal number of parameters to use for each kind of approximation.
+# Plots if true.
+shouldPlotMaxError = True
+
+def plot_error_range(x_coordinates, y_approximation, y_actual):
+    """ Given x coordinates that correspond to approximated y coordinates and actual y coordinates,
+    compute the deltas between y approximated and y actual and plot them in log scale on y.
+    """
+    error_deltas = approximations.compute_error_range(y_approximation, y_actual)
+
+    plt.semilogy(x_coordinates, error_deltas)
+
+    plt.grid(True)
+    plt.title(f"Chebyshev Rational e^x Errors on [{x_start}, {x_end}]. {num_parameters} params, {num_points_plot} points")
+    plt.show()
+
+# This script does the following:
+# - Computes polynomial and rational approximations of a given function (e^x by default).
+# - Computes (x,y) coordinates for every approximation given the same x coordinates.
+# - Plots the results for rough comparison.
+# - Computes the max error for every approximation given the same x coordinates.
+# - Computes and plots max errors for every approximation with a varying number of parameters.
+# This script runs various approximation methods, plots their results and deltas
+# from actual function values. The script can also be configured to print the maximum error.
+# The exact behavior is controlled by the global
+# variables at the top of the file. 
+# The following are the resources used:
+# https://xn--2-umb.com/22/approximation/
+# https://sites.tufts.edu/atasissa/files/2019/09/remez.pdf
+def main():
+    # Equispaced x coordinates to be used for plotting every approximation.
+    x_coordinates = approximations.linspace(x_start, x_end, num_points_plot)
+
+    if shouldComputeErrorDelta or shouldPlotApproximations or shouldPlotErrorRange:
+        ###############################################
+        # Approximation With Given Number of Parameters
+        y_eqispaced_poly, y_chebyshev_poly, y_chebyshev_rational, y_actual = approximations.approx_and_eval_all(approximated_fn, num_parameters, x_coordinates)
+
+        ################
+        # Compute Errors
+        if shouldComputeErrorDelta:
+            print(F"\n\nMax Error on [{x_start}, {x_end}]")
+            print(F"{num_points_plot} coordinates equally spaced on the X axis")
+            print(F"{num_parameters} parameters used\n\n")
+
+            # Equispaced Polynomial Approximation
+            max_error_equispaced_poly = approximations.compute_max_error(y_eqispaced_poly, y_actual)
+            print(F"Equispaced Poly: {max_error_equispaced_poly.evalf(chop=1e-100)}")
+
+            # Chebyshev Polynomial Approximation
+            max_error_chebyshev_poly = approximations.compute_max_error(y_chebyshev_poly, y_actual)
+            print(F"Chebyshev Poly: {max_error_chebyshev_poly.evalf(chop=1e-100)}")
+
+            # Chebyshev Rational Approximation
+            max_error_chebyshev_rational = approximations.compute_max_error(y_chebyshev_rational, y_actual)
+            print(F"Chebyshev Rational: {max_error_chebyshev_rational.evalf(chop=1e-100)}")
+
+        if shouldPlotErrorRange:
+            plot_error_range(x_coordinates, y_chebyshev_rational, y_actual)
+
+        ###############################
+        # Plot Every Approximation Kind
+        if shouldPlotApproximations:
+            # Equispaced Polynomial Approximation
+            plt.plot(x_coordinates, y_eqispaced_poly, label="Equispaced Poly")
+
+            # Chebyshev Polynomial Approximation
+            plt.plot(x_coordinates, y_chebyshev_poly, label="Chebyshev Poly")
+
+            # Chebyshev Rational Approximation
+            plt.plot(x_coordinates, y_chebyshev_rational, label="Chebyshev Rational")
+
+            # Actual With Large Number of Coordinates (evenly spaced on the X-axis)
+            plt.plot(x_coordinates, y_actual, label=F"Actual")
+
+            plt.legend(loc="upper left")
+            plt.grid(True)
+            plt.title(f"Appproximation of e^x on [{x_start}, {x_end}] with {num_parameters} parameters")
+            plt.show()
+
+    #####################################################
+    # Calculate Errors Given Varying Number of Parameters
+    if shouldPlotMaxError:
+        x_axis = []
+
+        deltas_eqispaced_poly = []
+        deltas_chebyshev_poly = []
+        deltas_chebyshev_rational = []
+
+        ################
+        # Compute Deltas
+        # The deltas are taken from actual function values for different number of parameters
+        # This is needed to find the most optimal number of parameters to use.
+        for num_parameters_current in range(1, num_parameters_errors + 1):
+            x_axis.append(int(num_parameters_current))
+            y_eqispaced_poly, y_chebyshev_poly, y_chebyshev_rational, y_actual = approximations.approx_and_eval_all(approximated_fn, num_parameters_current, x_coordinates)
+
+            deltas_eqispaced_poly.append(approximations.compute_max_error(y_eqispaced_poly, y_actual))
+            deltas_chebyshev_poly.append(approximations.compute_max_error(y_chebyshev_poly, y_actual))
+            deltas_chebyshev_rational.append(approximations.compute_max_error(y_chebyshev_rational, y_actual))
+
+        ##################
+        # Plot the results
+
+        # Equispaced Polynomial Approximation
+        plt.semilogy(x_axis, deltas_eqispaced_poly, label="Equispaced Poly")
+
+        # Chebyshev Polynomial Approximation
+        plt.semilogy(x_axis, deltas_chebyshev_poly, label="Chebyshev Poly")
+
+        # Chebyshev Rational Approximation
+        plt.semilogy(x_axis, deltas_chebyshev_rational, label="Chebyshev Rational")
+
+        plt.legend(loc="upper left")
+        plt.grid(True)
+        plt.title(f"Approximation Errors on [{x_start}, {x_end}]")
+        plt.gca().invert_yaxis()
+        plt.xlabel('Number of Parameters')
+        plt.ylabel(F"-log_10{{ max | f'(x) - f(x) | }} where x is in [{x_start}, {x_end}]")
+        plt.show()
+
+# This script isolates the 13-parameter Chebyshev Rational approximation of e^x
+# We are planning to use it in production. Therefore, we need to peform coefficient
+# truncations to 36 decimal points (the max osmomath supported precision).
+def exponent_approximation_choice():
+    # Equispaced x coordinates to be used for plotting every approximation.
+    x_coordinates = approximations.linspace(x_start, x_end, num_points_plot)
+    x_coordinates = [sp.Float(sp.N(coef, osmomath_precision + 1), osmomath_precision + 1) for coef in x_coordinates]
+
+    # Chebyshev Rational Approximation to get the coefficients.
+    coef_numerator, coef_denominator = approximations.chebyshev_rational_approx(approximated_fn, x_start, x_end, num_parameters)
+
+    # Truncate the coefficients to osmomath precision.
+    coef_numerator = [sp.Float(sp.N(coef, osmomath_precision + 1), osmomath_precision + 1) for coef in coef_numerator]
+    coef_denominator = [sp.Float(sp.N(coef, osmomath_precision + 1), osmomath_precision + 1) for coef in coef_denominator]
+
+    # Evaluate approximation.
+    y_chebyshev_rational = rational.evaluate(x_coordinates, coef_numerator, coef_denominator)
+
+    # Compute Actual Values
+    y_actual = approximations.get_y_actual(approximated_fn, x_coordinates)
+
+    plot_error_range(x_coordinates, y_chebyshev_rational, y_actual)
+
+if __name__ == "__main__":
+    # Uncomment to run the main script.
+    #main()
+    exponent_approximation_choice()
diff --git a/scripts/approximations/polynomial.py b/scripts/approximations/polynomial.py
new file mode 100644
index 000000000..5c6af2043
--- /dev/null
+++ b/scripts/approximations/polynomial.py
@@ -0,0 +1,41 @@
+import sympy as sp
+
+def construct_vandermonde_matrix(x_list: list[sp.Float]) -> sp.Matrix:
+    """ Constructs a Vandermonde matrix for a polynomial approximation.
+    from the list of x values given.
+    
+    len(x_list) * len(x_list)
+    x_list is the list of all x values to construct the matrix from.
+    
+    V = [1 x_0 x_0^2 ... x_0^{n-1}]
+    [1 x_1 x_2^1 ... x_1^{n-1}]
+    ...
+    [1 x_n x_n^2 ... x_n^{n-1}]
+
+    Vandermonde matrix is a matrix with the terms of a geometric progression in each row.
+    We use Vandermonde matrix to convert the system of linear equations into a linear algebra problem
+    that we can solve.
+    """
+    num_terms = len(x_list)
+
+    matrix = []
+
+    for i in range(num_terms):
+        row = []
+        for j in range(num_terms):
+            row.append(sp.Pow(x_list[i], j))
+        matrix.append(row)
+
+    return sp.Matrix(matrix)
+
+def evaluate(x, coeffs):
+    """ Evaluates the polynomial. Given a list of x coordinates and a list of coefficients, returns a list of
+    y coordinates, one for each x coordinate. The coefficients must be in ascending order.
+    """
+    y = []
+    for x_i in x:
+        y_i = 0
+        for i in range(len(coeffs)):
+            y_i += coeffs[i]*sp.Pow(x_i, i)
+        y.append(y_i)
+    return y
diff --git a/scripts/approximations/rational.py b/scripts/approximations/rational.py
new file mode 100644
index 000000000..4adea2d91
--- /dev/null
+++ b/scripts/approximations/rational.py
@@ -0,0 +1,44 @@
+import sympy as sp
+
+import polynomial
+
+def construct_rational_eval_matrix(x_list: list, y_list: list, num_terms_numerator: int, num_terms_denominator) -> sp.Matrix:
+    """ Constructs a matrix to use for computing coefficients for a rational approximation
+        from the list of x and y values given.
+        len(x_list) * len(x_list)
+        x_list is the list of all x values to construct the matrix from.
+        V = [1 x_0 x_0^2 ... x_0^{n-1} - y_0*x_0 - y_0*x_0^2 ... - y_0*x_0^{n-1}]
+            [1 x_1 x_2^1 ... x_1^{n-1} - y_1*x_1 - y_1*x_1^2 ... - y_1*x_1^{n-1}]
+            ...
+            [1 x_n x_n^2 ... x_n^{n-1} - y_n*x_n - y_n*x_n^2 ... - y_n*x_n^{n-1}]
+    """
+    matrix = []
+
+    for i in range(num_terms_numerator + num_terms_denominator - 1):
+        row = []
+        for j in range(num_terms_numerator):
+            row.append(sp.Pow(x_list[i], j))
+
+        for j in range(num_terms_denominator):
+            # denominator terms
+            if j > 0:
+                row.append(-1 * sp.Pow(x_list[i], j) * y_list[i])
+
+        matrix.append(row)
+
+    return sp.Matrix(matrix)
+
+def evaluate(x: list, coefficients_numerator: list, coefficients_denominator: list):
+    """ Evaluates the rational function. Assume rational h(x) = p(x) / q(x)
+    Given a list of x coordinates, a list of coefficients of p(x) - coefficients_numerator, and a list of
+    coefficients of q(x) - coefficients_denominator, returns a list of y coordinates, one for each x coordinate.
+    """
+    p = polynomial.evaluate(x, coefficients_numerator)
+    q = polynomial.evaluate(x, coefficients_denominator)
+
+    result = []
+
+    for i in range(len(x)):
+        result.append(p[i] / q[i])
+
+    return result
diff --git a/scripts/approximations/requirements.txt b/scripts/approximations/requirements.txt
new file mode 100644
index 000000000..ce1140c3e
--- /dev/null
+++ b/scripts/approximations/requirements.txt
@@ -0,0 +1,16 @@
+contourpy==1.0.6
+cycler==0.11.0
+fonttools==4.38.0
+kiwisolver==1.4.4
+matplotlib==3.6.2
+mpmath==1.2.1
+numpy==1.23.5
+packaging==22.0
+Pillow==9.3.0
+pyparsing==3.0.9
+PyQt5==5.15.7
+PyQt5-Qt5==5.15.2
+PyQt5-sip==12.11.0
+python-dateutil==2.8.2
+six==1.16.0
+sympy==1.11.1


Test patch:
diff --git a/osmomath/decimal_test.go b/osmomath/decimal_test.go
index 1b09c9f0e..d881f86ba 100644
--- a/osmomath/decimal_test.go
+++ b/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)
 		})
 	}
 }
diff --git a/scripts/approximations/rational_test.py b/scripts/approximations/rational_test.py
new file mode 100644
index 000000000..099b15d5a
--- /dev/null
+++ b/scripts/approximations/rational_test.py
@@ -0,0 +1,36 @@
+import unittest
+import math
+import sympy as sp
+
+import rational
+
+class TestChebyshevRational(unittest.TestCase):
+
+    # assume (3,3) rational 
+    # h(x) = p(x) / q(x) = (p_0 + p_1 * x + p_2 * x^2) / (1 + q_1 * x + q_2 * x^2)
+    #  p_0 + p_1 * x + p_2 * x^2 - q_1 * y * x - q_2 * y * x^2 = y
+    # construct 5 x 5 matrix solving for p and q given x_0 to x_4 and y_0 to y_4
+    # assume function is y = e**x
+    def test_construct_matrix_3_3(self):
+        fn = lambda x: math.e**x
+        x = [1, 2, 3, 4, 5]
+        y = list(map(fn, x))
+
+        coeffs = rational.construct_rational_eval_matrix(x, y, 3, 3)
+        
+        # correct matrix size
+        self.assertEqual(len(x) * len(x), len(coeffs))
+
+        # first row is correct
+        for i in range(len(x)):
+            x_i = x[i]
+            y_i = y[i]
+
+            expected_row = [sp.Pow(x_i, 0), sp.Pow(x_i, 1), sp.Pow(x_i, 2), -1 * sp.Pow(x_i, 1) * y_i, -1 * sp.Pow(x_i, 2) * y_i]
+
+            actual_row = coeffs.row(i).tolist()[0]
+
+            self.assertEqual(expected_row, actual_row)
+
+if __name__ == '__main__':
+    unittest.main()