Repository: osmosis-labs/osmosis
Pull request #2847 (2022-09-24 17:51:49)

[stableswap]: Combine `solveCfmm` and `solveCfmmMulti`
<!-- < < < < < < < < < < < < < < < < < < < < < < < < < < < < < < < < < ☺
v                               ✰  Thanks for creating a PR! ✰    
v    Before smashing the submit button please review the checkboxes.
v    If a checkbox is n/a - please still include it but + a little note why
v    If your PR doesn't close an issue, that's OK!  Just remove the Closes: #XXX line!
☺ > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > >  -->

Closes: #2845

## What is the purpose of the change

This PR merges `solveCfmmMulti` into `solveCfmm` and updates tests and pool logic accordingly.

Note: majority of diff is from #2843 and test refactors (one side benefit of this change is that it allows our tests to be _drastically_ cleaner and more readable (no more manual `uReserve` and `wSumSquares` calculations!))

Context from issue:

> All pool-level CFMM calls currently only use solveCfmm, which is only connected to our two-asset CFMM solver. Many of these functions also hard code the assumption that there are only two assets (e.g. by always specifically taking the [first two](https://github.com/osmosis-labs/osmosis/blob/39b32018d736f93fbb4ecdc85d53dae457eec5a1/x/gamm/pool-models/stableswap/amm.go#L379) assets in pool reserves).
> 
> We should combine these two solver functions and properly direct inputs towards the correct solver depending on the number of assets they have.

## Brief Changelog

- Merge `solveCfmmMulti` into `solveCfmm`
- Change `solveCfmm` to take a single `sdk.Coins` or `PoolAssets` input (latter might be safer if possible)
- Add logic in `solveCfmm` that routes to two-asset solver if only two assets and to multi-asset solver otherwise. This logic should calculate `uReserve` and `wSumSquares` internally from the input tokens to prevent unintuitive calculations elsewhere in the codebase and to prevent user error around these calcs in the future
- Update pool function logic to accommodate multi-asset pools (previously would just take first two tokens in token reserves)

## Testing and Verifying

All tests previously relating to `solveCfmm` and `solveCfmmMulti` should now be properly refactored to use the new `solveCfmm`

## Documentation and Release Note

  - Does this pull request introduce a new feature or user-facing behavior changes? (no)
  - Is a relevant changelog entry added to the `Unreleased` section in `CHANGELOG.md`? (no)
  - How is the feature or change documented? (not documented)

ValarDragon: Do we not want to keep the direct solvers at all?

Patch:
diff --git a/osmomath/decimal.go b/osmomath/decimal.go
index 3125d64fa..1e82faf19 100644
--- a/osmomath/decimal.go
+++ b/osmomath/decimal.go
@@ -516,12 +516,22 @@ func (d BigDec) SDKDec() sdk.Dec {
 	return truncatedDec
 }
 
-// BigDecFromSdkDec returns the BigDec representation of an SdkDec.
+// BigDecFromSdkDec returns the BigDec representation of an SDKDec.
 // Values in any additional decimal places are truncated.
 func BigDecFromSDKDec(d sdk.Dec) BigDec {
 	return NewDecFromBigIntWithPrec(d.BigInt(), sdk.Precision)
 }
 
+// BigDecFromSdkDecSlice returns the []BigDec representation of an []SDKDec.
+// Values in any additional decimal places are truncated.
+func BigDecFromSDKDecSlice(ds []sdk.Dec) []BigDec {
+	result := make([]BigDec, len(ds))
+	for i, d := range ds {
+		result[i] = NewDecFromBigIntWithPrec(d.BigInt(), sdk.Precision)
+	}
+	return result
+}
+
 //     ____
 //  __|    |__   "chop 'em
 //       ` \     round!"
diff --git a/x/gamm/pool-models/stableswap/amm.go b/x/gamm/pool-models/stableswap/amm.go
index 6256a0d94..10ad1929c 100644
--- a/x/gamm/pool-models/stableswap/amm.go
+++ b/x/gamm/pool-models/stableswap/amm.go
@@ -43,238 +43,25 @@ func cfmmConstantMulti(xReserve, yReserve, uReserve, vSumSquares osmomath.BigDec
 	return xyu.Mul(x2.Add(y2).Add(vSumSquares))
 }
 
-// solidly CFMM is xy(x^2 + y^2) = k
+// solidly CFMM is xy(x^2 + y^2) = k, and our multi-asset CFMM is xyz(x^2 + y^2 + w) = k
 // So we want to solve for a given addition of `b` units of y into the pool,
 // how many units `a` of x do we get out.
-// So we solve the following expression for `a`
+// So we solve the following expression for `a` in two-asset pools:
 // xy(x^2 + y^2) = (x - a)(y + b)((x - a)^2 + (y + b)^2)
-func solveCfmm(xReserve, yReserve, yIn osmomath.BigDec) osmomath.BigDec {
-	if !xReserve.IsPositive() || !yReserve.IsPositive() || !yIn.IsPositive() {
-		panic("invalid input: reserves and input must be positive")
-	}
-
-	// use the following wolfram alpha link to solve the equation
-	// https://www.wolframalpha.com/input?i=solve+for+a%2C+xy%28x%5E2+%2B+y%5E2%29+%3D+%28x+-+a%29%28y+%2B+b%29%28%28x+-+a%29%5E2+%2B+%28y+%2Bb%29%5E2%29+
-	// This returns (copied from wolfram):
-	// assuming (correctly) that b + y!=0
-	// a = (-27 b^2 x^3 y - 27 b^2 x y^3 + sqrt((-27 b^2 x^3 y - 27 b^2 x y^3 - 54 b x^3 y^2 - 54 b x y^4 - 27 x^3 y^3 - 27 x y^5)^2 + 4 (3 b^4 + 12 b^3 y + 18 b^2 y^2 + 12 b y^3 + 3 y^4)^3) - 54 b x^3 y^2 - 54 b x y^4 - 27 x^3 y^3 - 27 x y^5)^(1/3)/(3 2^(1/3) (b + y)) - (2^(1/3) (3 b^4 + 12 b^3 y + 18 b^2 y^2 + 12 b y^3 + 3 y^4))/(3 (b + y) (-27 b^2 x^3 y - 27 b^2 x y^3 + sqrt((-27 b^2 x^3 y - 27 b^2 x y^3 - 54 b x^3 y^2 - 54 b x y^4 - 27 x^3 y^3 - 27 x y^5)^2 + 4 (3 b^4 + 12 b^3 y + 18 b^2 y^2 + 12 b y^3 + 3 y^4)^3) - 54 b x^3 y^2 - 54 b x y^4 - 27 x^3 y^3 - 27 x y^5)^(1/3)) + (b x + x y)/(b + y) and b + y!=0
-	// We simplify and separate out terms to get that its the following:
-	// The key substitutions are that 3(b+y)^4 = 3 b^4 + 12 b^3 y + 18 b^2 y^2 + 12 b y^3 + 3 y^4
-	// and -27 x y (b + y)^2 (x^2 + y^2) = -27 b^2 x^3 y - 27 b^2 x y^3 - 54 b x^3 y^2 - 54 b x y^4 - 27 x^3 y^3 - 27 x y^5
-	// I added {} myself, making better distinctions between entirely distinct terms.
-	// a = {(-27 x y (b + y)^2 (x^2 + y^2)
-	//			+ sqrt(
-	//				(-27 x y (b + y)^2 (x^2 + y^2))^2
-	//				+ 108 ((b+y)^4)^3
-	// 			)^(1/3)
-	// 		  / (3 2^(1/3) (b + y))}
-	//		- {(2^(1/3) (3 (b + y)^4))
-	//		  /(3 (b + y)
-	// 			(-27 x y (b + y)^2 (x^2 + y^2)
-	// 				+ sqrt(
-	//					(-27 x y (b + y)^2 (x^2 + y^2))^2
-	//					+ 108 ((b+y)^4)^3)
-	// 			)^(1/3))}
-	//      + {(b x + x y)/(b + y)}
-	// we further simplify, and call:
-	// foo = (-27 x y (b + y)^2 (x^2 + y^2)
-	// 			+ sqrt(
-	//				(-27 x y (b + y)^2 (x^2 + y^2))^2
-	//				+ 108 ((b+y)^4)^3)
-	//		 )^(1/3)
-	// Thus, a is then:
-	// a = {foo / (3 2^(1/3) (b + y))}
-	//		- {(3 * 2^(1/3) (b+y)^4)
-	//		  /(3 (b + y) foo)}
-	//      + {(b x + x y)/(b + y)}
-	// Let:
-	// term1 = {foo / (3 2^(1/3) (b + y))}
-	// term2 = {(3 * 2^(1/3) (b+y)^4) /(3 (b + y) foo)} =  2^(1/3) (b+y)^3 / foo
-	// term3 = {(b x + x y)/(b + y)}
-
-	// prelude, compute all the xy cross terms. Consider keeping these precomputed in the struct,
-	// and maybe in state.
-	x := xReserve
-	y := yReserve
-	x2py2 := x.Mul(x).Add(y.Mul(y))
-
-	xy := x.Mul(y)
-
-	b := yIn
-
-	bpy := b.Add(y)
-	bpy2 := bpy.Mul(bpy)
-	bpy3 := bpy2.Mul(bpy)
-	bpy4 := bpy2.Mul(bpy2)
-
-	// TODO: Come back and optimize alot of the calculations
-
-	// Now we compute foo
-	// foo = (-27 x y (b + y)^2 (x^2 + y^2)
-	// 			+ sqrt(
-	//				(-27 x y (b + y)^2 (x^2 + y^2))^2
-	//				+ 108 ((b+y)^4)^3)
-	//		 )^(1/3)
-	// This has a y^12 term in it, which is unappealing, so we spend some energy reducing this max bitlen.
-	// foo = (-27 x y (b + y)^2 (x^2 + y^2)
-	// 			+ (b + y)^2 sqrt(
-	//				729 (x y (x^2 + y^2))^2
-	//				+ 108 (b+y)^8)
-	//		 )^(1/3)
-	// let e = x y (x^2 + y^2))
-	// foo = (-27 (b + y)^2 e
-	// 			+ (b + y)^2 sqrt(
-	//				729 e^2 + 108 (b+y)^8)
-	//		 )^(1/3)
-
-	e := xy.Mul(x2py2) // xy(x^2 + y^2)
-
-	// t1 = -27 (b + y)^2 e
-	t1 := e.Mul(bpy2).MulInt64(-27)
-
-	// compute d = (b + y)^2 sqrt(729 e^2 + 108 (b+y)^8)
-	bpy8 := bpy4.Mul(bpy4)
-	sqrt_inner := e.Mul(e).MulInt64(729).Add(bpy8.MulInt64(108)) // 729 e^2 + 108 (b+y)^8
-	sqrt, err := sqrt_inner.ApproxSqrt()
-	if err != nil {
-		panic(err)
-	}
-	d := sqrt.Mul(bpy2)
-
-	// foo = (t1 + d)^(1/3)
-	foo3 := t1.Add(d)
-	foo, _ := foo3.ApproxRoot(3)
-
-	// a = {foo / (3 2^(1/3) (b + y))}
-	//		- {(2^(1/3) banana) / (3 (b + y) foo}
-	//      + {(b x + x y)/(b + y)}
-
-	// term1 := {foo / (3 2^(1/3) (b + y))}
-	term1Denominator := threeCubeRootTwo.Mul(bpy)
-	term1 := foo.Quo(term1Denominator)
-	// term2 := {(2^(1/3) (b+y)^3) / (foo}
-	term2 := cubeRootTwo.Mul(bpy3)
-	term2 = term2.Quo(foo)
-	// term3 := {(b x + x y)/(b + y)}
-	term3Numerator := b.Mul(x).Add(xy)
-	term3 := term3Numerator.Quo(bpy)
-
-	a := term1.Sub(term2).Add(term3)
-
-	if a.GTE(xReserve) {
-		panic("invalid output: greater than full pool reserves")
-	}
-
-	return a
-}
-
-// Our multi-asset CFMM is xyz(x^2 + y^2 + w) = k
-// So we want to solve for a given addition of `b` units of y into the pool,
-// how many units `a` of x do we get out.
-// So we solve the following expression for `a`
+// and the following expression for `a` in multi-asset pools:
 // xyz(x^2 + y^2 + w) = (x - a)(y + b)z((x - a)^2 + (y + b)^2 + w)
-func solveCfmmMulti(xReserve, yReserve, wSumSquares, yIn osmomath.BigDec) osmomath.BigDec {
-	if !xReserve.IsPositive() || !yReserve.IsPositive() || !yIn.IsPositive() {
-		panic("invalid input: reserves and input must be positive")
-	}
-
-	// Use the following wolfram alpha link to solve the equation
-	// https://www.wolframalpha.com/input?i=solve+for+a%2C+xyz%28x%5E2+%2B+y%5E2+%2B+w%29+%3D+%28x+-+a%29%28y+%2B+b%29z%28%28x+-+a%29%5E2+%2B+%28y+%2Bb%29%5E2+%2B+w%29
-	// This returns (copied from wolfram):
-	// assuming (correctly) that b + y!=0
-	// a = (-27 b^2 w x y - 27 b^2 x^3 y - 27 b^2 x y^3 + sqrt((-27 b^2 w x y - 27 b^2 x^3 y - 27 b^2 x y^3 - 54 b w x y^2 - 54 b x^3 y^2 - 54 b x y^4 - 27 w x y^3 - 27 x^3 y^3 - 27 x y^5)^2 + 4 (3 b^4 + 12 b^3 y + 3 b^2 w + 18 b^2 y^2 + 6 b w y + 12 b y^3 + 3 w y^2 + 3 y^4)^3) - 54 b w x y^2 - 54 b x^3 y^2 - 54 b x y^4 - 27 w x y^3 - 27 x^3 y^3 - 27 x y^5)^(1/3)/(3 2^(1/3) (b + y)) - (2^(1/3) (3 b^4 + 12 b^3 y + 3 b^2 w + 18 b^2 y^2 + 6 b w y + 12 b y^3 + 3 w y^2 + 3 y^4))/(3 (b + y) (-27 b^2 w x y - 27 b^2 x^3 y - 27 b^2 x y^3 + sqrt((-27 b^2 w x y - 27 b^2 x^3 y - 27 b^2 x y^3 - 54 b w x y^2 - 54 b x^3 y^2 - 54 b x y^4 - 27 w x y^3 - 27 x^3 y^3 - 27 x y^5)^2 + 4 (3 b^4 + 12 b^3 y + 3 b^2 w + 18 b^2 y^2 + 6 b w y + 12 b y^3 + 3 w y^2 + 3 y^4)^3) - 54 b w x y^2 - 54 b x^3 y^2 - 54 b x y^4 - 27 w x y^3 - 27 x^3 y^3 - 27 x y^5)^(1/3)) + (b x + x y)/(b + y) and b + y!=0
-	//
-	// The key substitutions are (where w represents the sum of the squares as represented in the multi-asset CFMM function):
-	// 1. S1: 3 (b + y)^2 (b^2 + 2 b y + y^2 + w) = 3 b^4 + 12 b^3 y + 3b^2w + 18 b^2 y^2 + 6bwy + 12 b y^3 + 3wy^2 + 3 y^4
-	// 2. S2: -27 x y (b + y)^2 (x^2 + y^2 + w) = -27b^2xyw - 27b^2x^3y - 27b^2xy^3 - 54bxy^2w - 54bx^3y^2 - 54bxy^4 - 27xy^3w - 27x^3y^3 - 27xy^5
-	//
-	// This is the simplified version using the substitutions above, to be expanded afterwards below:
-	// a = (1 / [3 * 2^(1/3) * (b + y)])
-	// 			* [S2 + sqrt(S2^2 + 4*(S1^3))]^(1/3)
-	//			- {2^(1/3) * S1 /
-	// 				[3 * (b + y) * (S2 + sqrt(S2^2 + 4*(S1^3)))^(1/3)]}
-	// 			+ [(bx + xy) / (b + y)]
-	//
-	// To further simplify, let:
-	// foo = (S2 + sqrt(S2^2 + 4*(S1^3)))^(1/3)
-	// bpy = (b + y)
-	//
-	// Thus, a further simplifies to:
-	// a = (foo / (3 * 2^(1/3) * bpy))
-	//			- (2^(1/3) * S1 / (3 * bpy * foo))
-	// 			+ ((bx + xy) / bpy)
-	//
-	// Finally, let:
-	// term1 = (foo / (3 * 2^(1/3) * bpy))
-	// term2 = (2^(1/3) * S1 / (3 * bpy * foo))
-	// term3 = ((bx + xy) / bpy)
-	//
-	// The final result should be:
-	// a = term1 - term2 + term3
-
-	// Prelude, compute all the xy cross terms. Consider keeping these precomputed in the struct,
-	// and maybe in state.
-
-	x := xReserve
-	x2 := x.Mul(x)
-	y := yReserve
-	y2 := y.Mul(y)
-	w := wSumSquares
-	b := yIn
-
-	xy := x.Mul(y)
-
-	bpy := b.Add(y)
-	bpy2 := bpy.Mul(bpy)
-	// bpy3 := bpy2.Mul(bpy)
-	// bpy4 := bpy2.Mul(bpy2)
-
-	// S1 = 3 (b + y)^2 (b^2 + 2 b y + y^2 + w)
-	s1_inner := bpy2.Add(w) // (b^2 + 2 b y + y^2 + w) = (b+y)^2 + w
-	s1 := bpy2.MulInt64(3).Mul(s1_inner)
-	// S2 = -27 x y (b + y)^2 (x^2 + y^2 + w)
-	s2_inner := x2.Add(y2).Add(w) // (x^2 + y^2 + w)
-	s2 := bpy2.MulInt64(-27).Mul(xy).Mul(s2_inner)
-
-	// Calculating foo directly causes an integer overflow due to having a y^12 term, so we expand and factor to reduce its max bitlen:
-	// Original foo = (S2 + sqrt(S2^2 + 4*(S1^3)))^(1/3)
-	// Expand:
-	// foo = (S2 + sqrt((-27 x y (b + y)^2 (x^2 + y^2 + w))^2 +
-	// 					4*(3 (b + y)^2 (b^2 + 2 b y + y^2 + w))^3)
-	// 		 )^(1/3)
-	// Factor (b+y)^2 out from within the square root:
-	// foo = (S2 + (b + y)^2
-	// 			* sqrt((-27 x y (x^2 + y^2 + w))^2 +
-	// 					4 * ((b + y)^2) * (3 (b^2 + 2 b y + y^2 + w))^3)
-	// 		 )^(1/3)
-	sqrt_inner_term1 := xy.MulInt64(-27).Mul(s2_inner)                // (-27 x y (x^2 + y^2 + w))
-	sqrt_inner_term1_2 := sqrt_inner_term1.Mul(sqrt_inner_term1)      // sqrt_inner_term1^2
-	s1_inner_3 := (s1_inner).Mul(s1_inner).Mul(s1_inner).MulInt64(27) // (3 (b^2 + 2 b y + y^2 + w))^3 = 27 s1_inner^3
-	sqrt_inner_term2 := s1_inner_3.MulInt64(4).Mul(bpy2)              // 4 * ((b + y)^2) * s1_inner_3
-	sqrt_inner := sqrt_inner_term1_2.Add(sqrt_inner_term2)            // sqrt_inner_term1 + sqrt_inner_term2
-	sqrt, err := sqrt_inner.ApproxSqrt()
-	if err != nil {
-		panic(err)
+func solveCfmm(xReserve, yReserve osmomath.BigDec, remReserves []osmomath.BigDec, yIn osmomath.BigDec) osmomath.BigDec {
+	if len(remReserves) == 0 {
+		return solveCFMMBinarySearch(cfmmConstant)(xReserve, yReserve, yIn)
 	}
-	foo3 := s2.Add(bpy2.Mul(sqrt))
-	foo, _ := foo3.ApproxRoot(3)
-
-	// term1 = (foo / (3 * 2^(1/3) * bar))
-	term1Denominator := threeCubeRootTwo.Mul(bpy) // 3 * 2^(1/3) * (b + y)
-	term1 := foo.Quo(term1Denominator)
-	// term2 = (2^(1/3) * S1 / (3 * bar * foo))
-	term2 := (cubeRootTwo.Mul(s1)).Quo(foo.Mul(bpy).MulInt64(3))
-	// term3 = ((bx + xy) / bpy)
-	term3Numerator := b.Mul(x).Add(xy)
-	term3 := term3Numerator.Quo(bpy)
 
-	a := term1.Sub(term2).Add(term3)
-
-	if a.GTE(xReserve) {
-		panic("invalid output: greater than full pool reserves")
+	uReserve := osmomath.OneDec()
+	wSumSquares := osmomath.ZeroDec()
+	for _, assetReserve := range remReserves {
+		uReserve = uReserve.Mul(assetReserve)
+		wSumSquares = wSumSquares.Add(assetReserve.Mul(assetReserve))
 	}
-
-	return a
+	return solveCFMMBinarySearchMulti(cfmmConstantMulti)(xReserve, yReserve, uReserve, wSumSquares, yIn)
 }
 
 func approxDecEqual(a, b, tol osmomath.BigDec) bool {
@@ -324,25 +111,29 @@ func solveCFMMBinarySearch(constantFunction func(osmomath.BigDec, osmomath.BigDe
 // added for future extension
 func solveCFMMBinarySearchMulti(constantFunction func(osmomath.BigDec, osmomath.BigDec, osmomath.BigDec, osmomath.BigDec) osmomath.BigDec) func(osmomath.BigDec, osmomath.BigDec, osmomath.BigDec, osmomath.BigDec, osmomath.BigDec) osmomath.BigDec {
 	return func(xReserve, yReserve, uReserve, wSumSquares, yIn osmomath.BigDec) osmomath.BigDec {
-		if !xReserve.IsPositive() || !yReserve.IsPositive() || wSumSquares.IsNegative() || !yIn.IsPositive() {
+		if !xReserve.IsPositive() || !yReserve.IsPositive() || !uReserve.IsPositive() || wSumSquares.IsNegative() || !yIn.IsPositive() {
 			panic("invalid input: reserves and input must be positive")
+		} else if yIn.GTE(yReserve) {
+			panic("cannot input more than pool reserves")
 		}
 		k := constantFunction(xReserve, yReserve, uReserve, wSumSquares)
-		yf := yReserve.Add(yIn)
-		x_low_est := osmomath.ZeroDec()
-		x_high_est := xReserve
-		x_est := (x_high_est.Add(x_low_est)).Quo(twodec)
-		cur_k := constantFunction(x_est, yf, uReserve, wSumSquares)
-		for i := 0; !approxDecEqual(cur_k, k, k_threshold) && i < 256; i++ {
-			if cur_k.GT(k) {
-				x_high_est = x_est
-			} else if cur_k.LT(k) {
-				x_low_est = x_est
-			}
-			x_est = (x_high_est.Add(x_low_est)).Quo(twodec)
-			cur_k = constantFunction(x_est, yf, uReserve, wSumSquares)
+		yFinal := yReserve.Add(yIn)
+		xLowEst := osmomath.ZeroDec()
+		xHighEst := xReserve
+		maxIterations := 256
+		errTolerance := osmoutils.ErrTolerance{AdditiveTolerance: sdk.OneInt(), MultiplicativeTolerance: sdk.Dec{}}
+
+		// create single-input CFMM to pass into binary search
+		calcXEst := func(xEst osmomath.BigDec) (osmomath.BigDec, error) {
+			return constantFunction(xEst, yFinal, uReserve, wSumSquares), nil
 		}
-		xOut := xReserve.Sub(x_est)
+
+		xEst, err := osmoutils.BinarySearchBigDec(calcXEst, xLowEst, xHighEst, k, errTolerance, maxIterations)
+		if err != nil {
+			panic(err)
+		}
+
+		xOut := xReserve.Sub(xEst)
 		if xOut.GTE(xReserve) {
 			panic("invalid output: greater than full pool reserves")
 		}
@@ -367,7 +158,7 @@ func spotPrice(baseReserve, quoteReserve osmomath.BigDec) osmomath.BigDec {
 	// xReserve & yReserve.
 	a := osmomath.OneDec()
 	// no need to divide by a, since a = 1.
-	return solveCfmm(baseReserve, quoteReserve, a)
+	return solveCfmm(baseReserve, quoteReserve, []osmomath.BigDec{}, a)
 }
 
 // returns outAmt as a decimal
@@ -377,8 +168,9 @@ func (p *Pool) calcOutAmtGivenIn(tokenIn sdk.Coin, tokenOutDenom string, swapFee
 		return sdk.Dec{}, err
 	}
 	tokenInSupply, tokenOutSupply := reserves[0], reserves[1]
+	remReserves := osmomath.BigDecFromSDKDecSlice(reserves[2:])
 	// We are solving for the amount of token out, hence x = tokenOutSupply, y = tokenInSupply
-	cfmmOut := solveCfmm(osmomath.BigDecFromSDKDec(tokenOutSupply), osmomath.BigDecFromSDKDec(tokenInSupply), osmomath.BigDecFromSDKDec(tokenIn.Amount.ToDec()))
+	cfmmOut := solveCfmm(osmomath.BigDecFromSDKDec(tokenOutSupply), osmomath.BigDecFromSDKDec(tokenInSupply), remReserves, osmomath.BigDecFromSDKDec(tokenIn.Amount.ToDec()))
 	outAmt := p.getDescaledPoolAmt(tokenOutDenom, cfmmOut)
 	return outAmt.SDKDec(), nil
 }
@@ -390,9 +182,10 @@ func (p *Pool) calcInAmtGivenOut(tokenOut sdk.Coin, tokenInDenom string, swapFee
 		return sdk.Dec{}, err
 	}
 	tokenInSupply, tokenOutSupply := reserves[0], reserves[1]
+	remReserves := osmomath.BigDecFromSDKDecSlice(reserves[2:])
 	// We are solving for the amount of token in, cfmm(x,y) = cfmm(x + x_in, y - y_out)
 	// x = tokenInSupply, y = tokenOutSupply, yIn = -tokenOutAmount
-	cfmmIn := solveCfmm(osmomath.BigDecFromSDKDec(tokenInSupply), osmomath.BigDecFromSDKDec(tokenOutSupply), osmomath.BigDecFromSDKDec(tokenOut.Amount.ToDec().Neg()))
+	cfmmIn := solveCfmm(osmomath.BigDecFromSDKDec(tokenInSupply), osmomath.BigDecFromSDKDec(tokenOutSupply), remReserves, osmomath.BigDecFromSDKDec(tokenOut.Amount.ToDec().Neg()))
 	inAmt := p.getDescaledPoolAmt(tokenInDenom, cfmmIn.Neg())
 	return inAmt.SDKDec(), nil
 }


Test patch:
diff --git a/osmomath/decimal_test.go b/osmomath/decimal_test.go
index 09a419e41..1650186dd 100644
--- a/osmomath/decimal_test.go
+++ b/osmomath/decimal_test.go
@@ -200,6 +200,32 @@ func (s *decimalTestSuite) TestBigDecFromSdkDec() {
 	}
 }
 
+func (s *decimalTestSuite) TestBigDecFromSdkDecSlice() {
+	tests := []struct {
+		d    []sdk.Dec
+		want []BigDec
+		expPanic bool
+	}{
+		{[]sdk.Dec{sdk.MustNewDecFromStr("0.000000000000000000")}, []BigDec{NewBigDec(0)}, false},
+		{[]sdk.Dec{sdk.MustNewDecFromStr("0.000000000000000000"), sdk.MustNewDecFromStr("1.000000000000000000")}, []BigDec{NewBigDec(0), NewBigDec(1)}, false},
+		{[]sdk.Dec{sdk.MustNewDecFromStr("1.000000000000000000"), sdk.MustNewDecFromStr("0.000000000000000000"), sdk.MustNewDecFromStr("0.000123400000000000")}, []BigDec{NewBigDec(1), NewBigDec(0), NewDecWithPrec(12340, 8)}, false},
+		{[]sdk.Dec{sdk.MustNewDecFromStr("10.000000000000000000")}, []BigDec{NewBigDec(10)}, false},
+		{[]sdk.Dec{sdk.MustNewDecFromStr("12340.000000000000000000")}, []BigDec{NewBigDec(12340)}, false},
+		{[]sdk.Dec{sdk.MustNewDecFromStr("1.234000000000000000"), sdk.MustNewDecFromStr("12340.000000000000000000")}, []BigDec{NewDecWithPrec(12340, 4), NewBigDec(12340)}, false},
+		{[]sdk.Dec{sdk.MustNewDecFromStr("0.123400000000000000"), sdk.MustNewDecFromStr("12340.000000000000000000")}, []BigDec{NewDecWithPrec(12340, 5), NewBigDec(12340)}, false},
+		{[]sdk.Dec{sdk.MustNewDecFromStr("0.000123400000000000"), sdk.MustNewDecFromStr("10.090090090090090090")}, []BigDec{NewDecWithPrec(12340, 8), NewDecWithPrec(1009009009009009009, 17)}, false},
+		{[]sdk.Dec{sdk.MustNewDecFromStr("10.090090090090090090"), sdk.MustNewDecFromStr("10.090090090090090090")}, []BigDec{NewDecWithPrec(1009009009009009009, 17), NewDecWithPrec(1009009009009009009, 17)}, false},
+	}
+	for tcIndex, tc := range tests {
+		if tc.expPanic {
+			s.Require().Panics(func() { BigDecFromSDKDecSlice(tc.d) })
+		} else {
+			value := BigDecFromSDKDecSlice(tc.d)
+			s.Require().Equal(tc.want, value, "bad BigDecFromSdkDec(), index: %v", tcIndex)
+		}
+	}
+}
+
 func (s *decimalTestSuite) TestEqualities() {
 	tests := []struct {
 		d1, d2     BigDec
diff --git a/x/gamm/pool-models/stableswap/amm_bench_test.go b/x/gamm/pool-models/stableswap/amm_bench_test.go
index a868b8f86..37eec728f 100644
--- a/x/gamm/pool-models/stableswap/amm_bench_test.go
+++ b/x/gamm/pool-models/stableswap/amm_bench_test.go
@@ -10,19 +10,43 @@ import (
 func BenchmarkCFMM(b *testing.B) {
 	// Uses solveCfmm
 	for i := 0; i < b.N; i++ {
-		runCalc(solveCfmm)
+		runCalcCFMM(solveCfmm)
 	}
 }
 
-func BenchmarkBinarySearch(b *testing.B) {
+func BenchmarkBinarySearchTwoAsset(b *testing.B) {
 	for i := 0; i < b.N; i++ {
-		runCalc(solveCFMMBinarySearch(cfmmConstant))
+		runCalcTwoAsset(solveCFMMBinarySearch(cfmmConstant))
 	}
 }
 
-func runCalc(solve func(osmomath.BigDec, osmomath.BigDec, osmomath.BigDec) osmomath.BigDec) {
+func BenchmarkBinarySearchMultiAsset(b *testing.B) {
+	for i := 0; i < b.N; i++ {
+		runCalcMultiAsset(solveCFMMBinarySearchMulti(cfmmConstantMulti))
+	}
+}
+
+func runCalcCFMM(solve func(osmomath.BigDec, osmomath.BigDec, []osmomath.BigDec, osmomath.BigDec) osmomath.BigDec) {
+	xReserve := osmomath.NewBigDec(rand.Int63n(100000) + 50000)
+	yReserve := osmomath.NewBigDec(rand.Int63n(100000) + 50000)
+	yIn := osmomath.NewBigDec(rand.Int63n(100000))
+	solve(xReserve, yReserve, []osmomath.BigDec{}, yIn)
+}
+
+func runCalcTwoAsset(solve func(osmomath.BigDec, osmomath.BigDec, osmomath.BigDec) osmomath.BigDec) {
 	xReserve := osmomath.NewBigDec(rand.Int63n(100000) + 50000)
 	yReserve := osmomath.NewBigDec(rand.Int63n(100000) + 50000)
 	yIn := osmomath.NewBigDec(rand.Int63n(100000))
 	solve(xReserve, yReserve, yIn)
 }
+
+func runCalcMultiAsset(solve func(osmomath.BigDec, osmomath.BigDec, osmomath.BigDec, osmomath.BigDec, osmomath.BigDec) osmomath.BigDec) {
+	xReserve := osmomath.NewBigDec(rand.Int63n(100000) + 50000)
+	yReserve := osmomath.NewBigDec(rand.Int63n(100000) + 50000)
+	mReserve := osmomath.NewBigDec(rand.Int63n(100000) + 50000)
+	nReserve := osmomath.NewBigDec(rand.Int63n(100000) + 50000)
+	u := mReserve.Mul(nReserve)
+	w := mReserve.Mul(mReserve).Add(nReserve.Mul(nReserve))
+	yIn := osmomath.NewBigDec(rand.Int63n(100000))
+	solve(xReserve, yReserve, u, w, yIn)
+}
diff --git a/x/gamm/pool-models/stableswap/amm_test.go b/x/gamm/pool-models/stableswap/amm_test.go
index 019bdcf62..7b5923d81 100644
--- a/x/gamm/pool-models/stableswap/amm_test.go
+++ b/x/gamm/pool-models/stableswap/amm_test.go
@@ -13,162 +13,386 @@ import (
 	"github.com/osmosis-labs/osmosis/v12/x/gamm/pool-models/internal/test_helpers"
 )
 
-// twoAssetCFMMTestCase defines a testcase for TestCFMMInvariantTwoAssets
-// and TestCFMMInvariantTwoAssetsDirect
-type twoAssetCFMMTestCase struct {
+// CFMMTestCase defines a testcase for stableswap pools
+type CFMMTestCase struct {
 	xReserve    osmomath.BigDec
 	yReserve    osmomath.BigDec
+	remReserves []osmomath.BigDec
 	yIn         osmomath.BigDec
 	expectPanic bool
 }
 
-var twoAssetCFMMTestCases = map[string]twoAssetCFMMTestCase{
-	// sanity checks
-	"small pool small input": {
-		xReserve:    osmomath.NewBigDec(100),
-		yReserve:    osmomath.NewBigDec(100),
-		yIn:         osmomath.NewBigDec(1),
-		expectPanic: false,
-	},
-	"small pool large input": {
-		xReserve:    osmomath.NewBigDec(100),
-		yReserve:    osmomath.NewBigDec(100),
-		yIn:         osmomath.NewBigDec(99),
-		expectPanic: false,
-	},
-	"medium pool medium join": {
-		xReserve:    osmomath.NewBigDec(100000),
-		yReserve:    osmomath.NewBigDec(100000),
-		yIn:         osmomath.NewBigDec(10000),
-		expectPanic: false,
-	},
-	"large pool medium join": {
-		xReserve:    osmomath.NewBigDec(10000000),
-		yReserve:    osmomath.NewBigDec(10000000),
-		yIn:         osmomath.NewBigDec(10000),
-		expectPanic: false,
-	},
-	"large pool large join": {
-		xReserve:    osmomath.NewBigDec(10000000),
-		yReserve:    osmomath.NewBigDec(10000000),
-		yIn:         osmomath.NewBigDec(1000000),
-		expectPanic: false,
-	},
-	"very large pool medium join": {
-		xReserve:    osmomath.NewBigDec(1000000000),
-		yReserve:    osmomath.NewBigDec(1000000000),
-		yIn:         osmomath.NewBigDec(100000),
-		expectPanic: false,
-	},
-	"billion token pool hundred million token join": {
-		xReserve:    osmomath.NewBigDec(1000000000),
-		yReserve:    osmomath.NewBigDec(1000000000),
-		yIn:         osmomath.NewBigDec(100000000),
-		expectPanic: false,
-	},
-
-	// uneven reserves
-	"xReserve double yReserve (small)": {
-		xReserve:    osmomath.NewBigDec(100),
-		yReserve:    osmomath.NewBigDec(50),
-		yIn:         osmomath.NewBigDec(1),
-		expectPanic: false,
-	},
-	"yReserve double xReserve (small)": {
-		xReserve:    osmomath.NewBigDec(50),
-		yReserve:    osmomath.NewBigDec(100),
-		yIn:         osmomath.NewBigDec(1),
-		expectPanic: false,
-	},
-	"xReserve double yReserve (large)": {
-		xReserve:    osmomath.NewBigDec(13789470),
-		yReserve:    osmomath.NewBigDec(59087324),
-		yIn:         osmomath.NewBigDec(1047829),
-		expectPanic: false,
-	},
-	"yReserve double xReserve (large)": {
-		xReserve:    osmomath.NewBigDec(50000000),
-		yReserve:    osmomath.NewBigDec(100000000),
-		yIn:         osmomath.NewBigDec(1000000),
-		expectPanic: false,
-	},
-	"uneven medium pool medium join": {
-		xReserve:    osmomath.NewBigDec(123456),
-		yReserve:    osmomath.NewBigDec(434245),
-		yIn:         osmomath.NewBigDec(23314),
-		expectPanic: false,
-	},
-	"uneven large pool medium join": {
-		xReserve:    osmomath.NewBigDec(11023432),
-		yReserve:    osmomath.NewBigDec(17432897),
-		yIn:         osmomath.NewBigDec(89734),
-		expectPanic: false,
-	},
-	"uneven large pool large join": {
-		xReserve:    osmomath.NewBigDec(38987364),
-		yReserve:    osmomath.NewBigDec(52893462),
-		yIn:         osmomath.NewBigDec(9819874),
-		expectPanic: false,
-	},
-	"uneven very large pool medium join": {
-		xReserve:    osmomath.NewBigDec(1473891748),
-		yReserve:    osmomath.NewBigDec(7438971234),
-		yIn:         osmomath.NewBigDec(100000),
-		expectPanic: false,
-	},
-	"uneven billion token pool billion token join": {
-		xReserve:    osmomath.NewBigDec(2678238934),
-		yReserve:    osmomath.NewBigDec(1573917894),
-		yIn:         osmomath.NewBigDec(5378748),
-		expectPanic: false,
-	},
-
-	// panic catching
-	"yIn greater than pool reserves": {
-		xReserve:    osmomath.NewBigDec(100),
-		yReserve:    osmomath.NewBigDec(100),
-		yIn:         osmomath.NewBigDec(1000),
-		expectPanic: true,
-	},
-	"xReserve negative": {
-		xReserve:    osmomath.NewBigDec(-100),
-		yReserve:    osmomath.NewBigDec(100),
-		yIn:         osmomath.NewBigDec(1),
-		expectPanic: true,
-	},
-	"yReserve negative": {
-		xReserve:    osmomath.NewBigDec(100),
-		yReserve:    osmomath.NewBigDec(-100),
-		yIn:         osmomath.NewBigDec(1),
-		expectPanic: true,
-	},
-	"yIn negative": {
-		xReserve:    osmomath.NewBigDec(100),
-		yReserve:    osmomath.NewBigDec(100),
-		yIn:         osmomath.NewBigDec(-1),
-		expectPanic: true,
-	},
-
-	// overflows
-	"xReserve near max bitlen": {
-		xReserve:    osmomath.NewDecFromBigInt(new(big.Int).Sub(new(big.Int).Exp(big.NewInt(2), big.NewInt(1024), nil), big.NewInt(1))),
-		yReserve:    osmomath.NewBigDec(100),
-		yIn:         osmomath.NewBigDec(1),
-		expectPanic: true,
-	},
-	"yReserve near max bitlen": {
-		xReserve:    osmomath.NewBigDec(100),
-		yReserve:    osmomath.NewDecFromBigInt(new(big.Int).Sub(new(big.Int).Exp(big.NewInt(2), big.NewInt(1024), nil), big.NewInt(1))),
-		yIn:         osmomath.NewBigDec(1),
-		expectPanic: true,
-	},
-	"both assets near max bitlen": {
-		xReserve:    osmomath.NewDecFromBigInt(new(big.Int).Sub(new(big.Int).Exp(big.NewInt(2), big.NewInt(1024), nil), big.NewInt(1))),
-		yReserve:    osmomath.NewDecFromBigInt(new(big.Int).Sub(new(big.Int).Exp(big.NewInt(2), big.NewInt(1024), nil), big.NewInt(1))),
-		yIn:         osmomath.NewBigDec(1),
-		expectPanic: true,
-	},
-}
+var (
+	overflowDec           = osmomath.NewDecFromBigInt(new(big.Int).Sub(new(big.Int).Exp(big.NewInt(2), big.NewInt(1024), nil), big.NewInt(1)))
+	twoAssetCFMMTestCases = map[string]CFMMTestCase{
+		// sanity checks
+		"small pool small input": {
+			xReserve:    osmomath.NewBigDec(100),
+			yReserve:    osmomath.NewBigDec(100),
+			remReserves: []osmomath.BigDec{},
+			yIn:         osmomath.NewBigDec(1),
+			expectPanic: false,
+		},
+		"small pool large input": {
+			xReserve:    osmomath.NewBigDec(100),
+			yReserve:    osmomath.NewBigDec(100),
+			remReserves: []osmomath.BigDec{},
+			yIn:         osmomath.NewBigDec(99),
+			expectPanic: false,
+		},
+		"medium pool medium join": {
+			xReserve:    osmomath.NewBigDec(100000),
+			yReserve:    osmomath.NewBigDec(100000),
+			remReserves: []osmomath.BigDec{},
+			yIn:         osmomath.NewBigDec(10000),
+			expectPanic: false,
+		},
+		"large pool medium join": {
+			xReserve:    osmomath.NewBigDec(10000000),
+			yReserve:    osmomath.NewBigDec(10000000),
+			remReserves: []osmomath.BigDec{},
+			yIn:         osmomath.NewBigDec(10000),
+			expectPanic: false,
+		},
+		"large pool large join": {
+			xReserve:    osmomath.NewBigDec(10000000),
+			yReserve:    osmomath.NewBigDec(10000000),
+			remReserves: []osmomath.BigDec{},
+			yIn:         osmomath.NewBigDec(1000000),
+			expectPanic: false,
+		},
+		"very large pool medium join": {
+			xReserve:    osmomath.NewBigDec(1000000000),
+			yReserve:    osmomath.NewBigDec(1000000000),
+			remReserves: []osmomath.BigDec{},
+			yIn:         osmomath.NewBigDec(100000),
+			expectPanic: false,
+		},
+		"billion token pool hundred million token join": {
+			xReserve:    osmomath.NewBigDec(1000000000),
+			yReserve:    osmomath.NewBigDec(1000000000),
+			remReserves: []osmomath.BigDec{},
+			yIn:         osmomath.NewBigDec(100000000),
+			expectPanic: false,
+		},
+
+		// uneven reserves
+		"xReserve double yReserve (small)": {
+			xReserve:    osmomath.NewBigDec(100),
+			yReserve:    osmomath.NewBigDec(50),
+			remReserves: []osmomath.BigDec{},
+			yIn:         osmomath.NewBigDec(1),
+			expectPanic: false,
+		},
+		"yReserve double xReserve (small)": {
+			xReserve:    osmomath.NewBigDec(50),
+			yReserve:    osmomath.NewBigDec(100),
+			remReserves: []osmomath.BigDec{},
+			yIn:         osmomath.NewBigDec(1),
+			expectPanic: false,
+		},
+		"xReserve double yReserve (large)": {
+			xReserve:    osmomath.NewBigDec(13789470),
+			yReserve:    osmomath.NewBigDec(59087324),
+			remReserves: []osmomath.BigDec{},
+			yIn:         osmomath.NewBigDec(1047829),
+			expectPanic: false,
+		},
+		"yReserve double xReserve (large)": {
+			xReserve:    osmomath.NewBigDec(50000000),
+			yReserve:    osmomath.NewBigDec(100000000),
+			remReserves: []osmomath.BigDec{},
+			yIn:         osmomath.NewBigDec(1000000),
+			expectPanic: false,
+		},
+		"uneven medium pool medium join": {
+			xReserve:    osmomath.NewBigDec(123456),
+			yReserve:    osmomath.NewBigDec(434245),
+			remReserves: []osmomath.BigDec{},
+			yIn:         osmomath.NewBigDec(23314),
+			expectPanic: false,
+		},
+		"uneven large pool medium join": {
+			xReserve:    osmomath.NewBigDec(11023432),
+			yReserve:    osmomath.NewBigDec(17432897),
+			remReserves: []osmomath.BigDec{},
+			yIn:         osmomath.NewBigDec(89734),
+			expectPanic: false,
+		},
+		"uneven large pool large join": {
+			xReserve:    osmomath.NewBigDec(38987364),
+			yReserve:    osmomath.NewBigDec(52893462),
+			remReserves: []osmomath.BigDec{},
+			yIn:         osmomath.NewBigDec(9819874),
+			expectPanic: false,
+		},
+		"uneven very large pool medium join": {
+			xReserve:    osmomath.NewBigDec(1473891748),
+			yReserve:    osmomath.NewBigDec(7438971234),
+			remReserves: []osmomath.BigDec{},
+			yIn:         osmomath.NewBigDec(100000),
+			expectPanic: false,
+		},
+		"uneven billion token pool billion token join": {
+			xReserve:    osmomath.NewBigDec(2678238934),
+			yReserve:    osmomath.NewBigDec(1573917894),
+			remReserves: []osmomath.BigDec{},
+			yIn:         osmomath.NewBigDec(5378748),
+			expectPanic: false,
+		},
+
+		// panic catching
+		"yIn greater than pool reserves": {
+			xReserve:    osmomath.NewBigDec(100),
+			yReserve:    osmomath.NewBigDec(100),
+			remReserves: []osmomath.BigDec{},
+			yIn:         osmomath.NewBigDec(1000),
+			expectPanic: true,
+		},
+		"xReserve negative": {
+			xReserve:    osmomath.NewBigDec(-100),
+			yReserve:    osmomath.NewBigDec(100),
+			remReserves: []osmomath.BigDec{},
+			yIn:         osmomath.NewBigDec(1),
+			expectPanic: true,
+		},
+		"yReserve negative": {
+			xReserve:    osmomath.NewBigDec(100),
+			yReserve:    osmomath.NewBigDec(-100),
+			remReserves: []osmomath.BigDec{},
+			yIn:         osmomath.NewBigDec(1),
+			expectPanic: true,
+		},
+		"yIn negative": {
+			xReserve:    osmomath.NewBigDec(100),
+			yReserve:    osmomath.NewBigDec(100),
+			remReserves: []osmomath.BigDec{},
+			yIn:         osmomath.NewBigDec(-1),
+			expectPanic: true,
+		},
+
+		// overflows
+		"xReserve near max bitlen": {
+			xReserve:    overflowDec,
+			yReserve:    osmomath.NewBigDec(100),
+			remReserves: []osmomath.BigDec{},
+			yIn:         osmomath.NewBigDec(1),
+			expectPanic: true,
+		},
+		"yReserve near max bitlen": {
+			xReserve:    osmomath.NewBigDec(100),
+			yReserve:    overflowDec,
+			remReserves: []osmomath.BigDec{},
+			yIn:         osmomath.NewBigDec(1),
+			expectPanic: true,
+		},
+		"both assets near max bitlen": {
+			xReserve:    overflowDec,
+			yReserve:    overflowDec,
+			remReserves: []osmomath.BigDec{},
+			yIn:         osmomath.NewBigDec(1),
+			expectPanic: true,
+		},
+	}
+
+	multiAssetCFMMTestCases = map[string]CFMMTestCase{
+		// sanity checks
+		"even 3-asset small pool, small input": {
+			xReserve: osmomath.NewBigDec(100),
+			yReserve: osmomath.NewBigDec(100),
+			// represents a 3-asset pool with 100 in each reserve
+			remReserves: []osmomath.BigDec{osmomath.NewBigDec(100)},
+			yIn:         osmomath.NewBigDec(1),
+			expectPanic: false,
+		},
+		"even 3-asset medium pool, small input": {
+			xReserve: osmomath.NewBigDec(100000),
+			yReserve: osmomath.NewBigDec(100000),
+			// represents a 3-asset pool with 100,000 in each reserve
+			remReserves: []osmomath.BigDec{osmomath.NewBigDec(100000)},
+			yIn:         osmomath.NewBigDec(100),
+			expectPanic: false,
+		},
+		"even 4-asset small pool, small input": {
+			xReserve: osmomath.NewBigDec(100),
+			yReserve: osmomath.NewBigDec(100),
+			// represents a 4-asset pool with 100 in each reserve
+			remReserves: []osmomath.BigDec{osmomath.NewBigDec(100), osmomath.NewBigDec(100)},
+			yIn:         osmomath.NewBigDec(1),
+			expectPanic: false,
+		},
+		"even 4-asset medium pool, small input": {
+			xReserve: osmomath.NewBigDec(100000),
+			yReserve: osmomath.NewBigDec(100000),
+			// represents a 4-asset pool with 100,000 in each reserve
+			remReserves: []osmomath.BigDec{osmomath.NewBigDec(100000), osmomath.NewBigDec(100000)},
+			yIn:         osmomath.NewBigDec(1),
+			expectPanic: false,
+		},
+		/* TODO: increase BigDec precision (36 -> 72) to be able to accommodate this
+		"even 4-asset large pool, small input": {
+			xReserve: osmomath.NewBigDec(100000000),
+			yReserve: osmomath.NewBigDec(100000000),
+			// represents a 4-asset pool with 100M in each reserve
+			remReserves: []osmomath.BigDec{osmomath.NewBigDec(100000000), osmomath.NewBigDec(100000000)},
+			yIn: osmomath.NewBigDec(100),
+			expectPanic: false,
+		},
+		*/
+
+		// uneven pools
+		"uneven 3-asset pool, even swap assets as pool minority": {
+			xReserve: osmomath.NewBigDec(100),
+			yReserve: osmomath.NewBigDec(100),
+			// the asset not being swapped has 100,000 token reserves (swap assets in pool minority)
+			remReserves: []osmomath.BigDec{osmomath.NewBigDec(100000)},
+			yIn:         osmomath.NewBigDec(10),
+			expectPanic: false,
+		},
+		"uneven 3-asset pool, uneven swap assets as pool minority, y > x": {
+			xReserve: osmomath.NewBigDec(100),
+			yReserve: osmomath.NewBigDec(200),
+			// the asset not being swapped has 100,000 token reserves (swap assets in pool minority)
+			remReserves: []osmomath.BigDec{osmomath.NewBigDec(100000)},
+			yIn:         osmomath.NewBigDec(10),
+			expectPanic: false,
+		},
+		"uneven 3-asset pool, uneven swap assets as pool minority, x > y": {
+			xReserve: osmomath.NewBigDec(200),
+			yReserve: osmomath.NewBigDec(100),
+			// the asset not being swapped has 100,000 token reserves (swap assets in pool minority)
+			remReserves: []osmomath.BigDec{osmomath.NewBigDec(100000)},
+			yIn:         osmomath.NewBigDec(10),
+			expectPanic: false,
+		},
+		"uneven 3-asset pool, no round numbers": {
+			xReserve: osmomath.NewBigDec(1178349),
+			yReserve: osmomath.NewBigDec(8329743),
+			// the asset not being swapped has 329,847 token reserves (swap assets in pool minority)
+			remReserves: []osmomath.BigDec{osmomath.NewBigDec(329847)},
+			yIn:         osmomath.NewBigDec(10),
+			expectPanic: false,
+		},
+		"uneven 4-asset pool, small input and swap assets in pool minority": {
+			xReserve: osmomath.NewBigDec(100),
+			yReserve: osmomath.NewBigDec(100),
+			// the assets not being swapped have 100,000 token reserves each (swap assets in pool minority)
+			remReserves: []osmomath.BigDec{osmomath.NewBigDec(100000), osmomath.NewBigDec(100000)},
+			yIn:         osmomath.NewBigDec(10),
+			expectPanic: false,
+		},
+		"uneven 4-asset pool, even swap assets in pool majority": {
+			xReserve: osmomath.NewBigDec(100000),
+			yReserve: osmomath.NewBigDec(100000),
+			// the assets not being swapped have 100 token reserves each (swap assets in pool majority)
+			remReserves: []osmomath.BigDec{osmomath.NewBigDec(100), osmomath.NewBigDec(100)},
+			yIn:         osmomath.NewBigDec(10),
+			expectPanic: false,
+		},
+		"uneven 4-asset pool, uneven swap assets in pool majority, y > x": {
+			xReserve: osmomath.NewBigDec(100000),
+			yReserve: osmomath.NewBigDec(200000),
+			// the assets not being swapped have 100 token reserves each (swap assets in pool majority)
+			remReserves: []osmomath.BigDec{osmomath.NewBigDec(100), osmomath.NewBigDec(100)},
+			yIn:         osmomath.NewBigDec(10),
+			expectPanic: false,
+		},
+		"uneven 4-asset pool, uneven swap assets in pool majority, y < x": {
+			xReserve: osmomath.NewBigDec(200000),
+			yReserve: osmomath.NewBigDec(100000),
+			// the assets not being swapped have 100 token reserves each (swap assets in pool majority)
+			remReserves: []osmomath.BigDec{osmomath.NewBigDec(100), osmomath.NewBigDec(100)},
+			yIn:         osmomath.NewBigDec(10),
+			expectPanic: false,
+		},
+		"uneven 4-asset pool, no round numbers": {
+			xReserve: osmomath.NewBigDec(1178349),
+			yReserve: osmomath.NewBigDec(8329743),
+			// the assets not being swapped have 329,847 tokens and 4,372,897 respectively
+			remReserves: []osmomath.BigDec{osmomath.NewBigDec(329847), osmomath.NewBigDec(4372897)},
+			yIn:         osmomath.NewBigDec(10),
+			expectPanic: false,
+		},
+
+		// panic catching
+		"negative xReserve": {
+			xReserve: osmomath.NewBigDec(-100),
+			yReserve: osmomath.NewBigDec(100),
+			// represents a 4-asset pool with 100 in each reserve
+			remReserves: []osmomath.BigDec{osmomath.NewBigDec(100), osmomath.NewBigDec(100)},
+			yIn:         osmomath.NewBigDec(1),
+			expectPanic: true,
+		},
+		"negative yReserve": {
+			xReserve: osmomath.NewBigDec(100),
+			yReserve: osmomath.NewBigDec(-100),
+			// represents a 4-asset pool with 100 in each reserve
+			remReserves: []osmomath.BigDec{osmomath.NewBigDec(100), osmomath.NewBigDec(100)},
+			yIn:         osmomath.NewBigDec(1),
+			expectPanic: true,
+		},
+		"negative remReserve": {
+			xReserve: osmomath.NewBigDec(100),
+			yReserve: osmomath.NewBigDec(100),
+			// represents a 4-asset pool with 100 in each reserve
+			remReserves: []osmomath.BigDec{osmomath.NewBigDec(-100), osmomath.NewBigDec(100)},
+			yIn:         osmomath.NewBigDec(1),
+			expectPanic: true,
+		},
+		"negative yIn": {
+			xReserve: osmomath.NewBigDec(100),
+			yReserve: osmomath.NewBigDec(100),
+			// represents a 4-asset pool with 100 in each reserve
+			remReserves: []osmomath.BigDec{},
+			yIn:         osmomath.NewBigDec(-1),
+			expectPanic: true,
+		},
+		"input greater than pool reserves (even 4-asset pool)": {
+			xReserve:    osmomath.NewBigDec(100),
+			yReserve:    osmomath.NewBigDec(100),
+			remReserves: []osmomath.BigDec{osmomath.NewBigDec(100), osmomath.NewBigDec(100)},
+			yIn:         osmomath.NewBigDec(1000),
+			expectPanic: true,
+		},
+
+		// overflows
+		"xReserve overflows in 4-asset pool": {
+			xReserve:    overflowDec,
+			yReserve:    osmomath.NewBigDec(100),
+			remReserves: []osmomath.BigDec{osmomath.NewBigDec(100), osmomath.NewBigDec(100)},
+			yIn:         osmomath.NewBigDec(1),
+			expectPanic: true,
+		},
+		"yReserve overflows in 4-asset pool": {
+			xReserve:    osmomath.NewBigDec(100),
+			yReserve:    overflowDec,
+			remReserves: []osmomath.BigDec{osmomath.NewBigDec(100), osmomath.NewBigDec(100)},
+			yIn:         osmomath.NewBigDec(1),
+			expectPanic: true,
+		},
+		"remReserve overflows in 3-asset pool": {
+			xReserve:    osmomath.NewBigDec(100),
+			yReserve:    osmomath.NewBigDec(100),
+			remReserves: []osmomath.BigDec{overflowDec},
+			yIn:         osmomath.NewBigDec(1),
+			expectPanic: true,
+		},
+		"remReserve overflows in 4-asset pool": {
+			xReserve:    osmomath.NewBigDec(100),
+			yReserve:    osmomath.NewBigDec(100),
+			remReserves: []osmomath.BigDec{osmomath.NewBigDec(100), overflowDec},
+			yIn:         osmomath.NewBigDec(1),
+			expectPanic: true,
+		},
+		"yIn overflows in 4-asset pool": {
+			xReserve:    osmomath.NewBigDec(100),
+			yReserve:    osmomath.NewBigDec(100),
+			remReserves: []osmomath.BigDec{osmomath.NewBigDec(100), osmomath.NewBigDec(100)},
+			yIn:         overflowDec,
+			expectPanic: true,
+		},
+	}
+)
 
 type StableSwapTestSuite struct {
 	test_helpers.CfmmCommonTestSuite
@@ -177,26 +401,21 @@ type StableSwapTestSuite struct {
 func TestCFMMInvariantTwoAssets(t *testing.T) {
 	kErrTolerance := osmomath.OneDec()
 
-	// TODO: switch solveCfmm to binary search and replace this with test case suite
-	tests := map[string]twoAssetCFMMTestCase{}
+	tests := twoAssetCFMMTestCases
 
 	for name, test := range tests {
 		t.Run(name, func(t *testing.T) {
 			// system under test
 			sut := func() {
+				// ensure there are only two assets
+				require.True(t, len(test.remReserves) == 0)
+
 				// using two-asset cfmm
 				k0 := cfmmConstant(test.xReserve, test.yReserve)
-				xOut := solveCfmm(test.xReserve, test.yReserve, test.yIn)
+				xOut := solveCfmm(test.xReserve, test.yReserve, test.remReserves, test.yIn)
 
 				k1 := cfmmConstant(test.xReserve.Sub(xOut), test.yReserve.Add(test.yIn))
 				osmomath.DecApproxEq(t, k0, k1, kErrTolerance)
-
-				// using multi-asset cfmm (should be equivalent with u = 1, w = 0)
-				k2 := cfmmConstantMulti(test.xReserve, test.yReserve, osmomath.OneDec(), osmomath.ZeroDec())
-				osmomath.DecApproxEq(t, k2, k0, kErrTolerance)
-				xOut2 := solveCfmmMulti(test.xReserve, test.yReserve, osmomath.ZeroDec(), test.yIn)
-				k3 := cfmmConstantMulti(test.xReserve.Sub(xOut2), test.yReserve.Add(test.yIn), osmomath.OneDec(), osmomath.ZeroDec())
-				osmomath.DecApproxEq(t, k2, k3, kErrTolerance)
 			}
 
 			osmoassert.ConditionalPanic(t, test.expectPanic, sut)
@@ -229,94 +448,43 @@ func TestCFMMInvariantTwoAssetsBinarySearch(t *testing.T) {
 func TestCFMMInvariantMultiAssets(t *testing.T) {
 	kErrTolerance := osmomath.OneDec()
 
-	tests := map[string]struct {
-		xReserve    osmomath.BigDec
-		yReserve    osmomath.BigDec
-		uReserve    osmomath.BigDec
-		wSumSquares osmomath.BigDec
-		yIn         osmomath.BigDec
-		expectPanic bool
-	}{
-		"4-asset pool, small input": {
-			osmomath.NewBigDec(100),
-			osmomath.NewBigDec(100),
-			// represents a 4-asset pool with 100 in each reserve
-			osmomath.NewBigDec(200),
-			osmomath.NewBigDec(20000),
-			osmomath.NewBigDec(1),
-			false,
-		},
-		"4-asset pool, large input": {
-			osmomath.NewBigDec(100),
-			osmomath.NewBigDec(100),
-			osmomath.NewBigDec(200),
-			osmomath.NewBigDec(20000),
-			osmomath.NewBigDec(1000),
-			false,
-		},
-		// This test fails due to a bug in our original solver
-		// "large pool, large input": {
-		// 	sdk.NewDec(100000),
-		// 	sdk.NewDec(100000),
-		// 	sdk.NewDec(10000),
-		// },
+	tests := multiAssetCFMMTestCases
 
-		// panic catching
-		"negative xReserve": {
-			osmomath.NewBigDec(-100),
-			osmomath.NewBigDec(100),
-			// represents a 4-asset pool with 100 in each reserve
-			osmomath.NewBigDec(200),
-			osmomath.NewBigDec(20000),
-			osmomath.NewBigDec(1),
-			true,
-		},
-		"negative yReserve": {
-			osmomath.NewBigDec(100),
-			osmomath.NewBigDec(-100),
-			// represents a 4-asset pool with 100 in each reserve
-			osmomath.NewBigDec(200),
-			osmomath.NewBigDec(20000),
-			osmomath.NewBigDec(1),
-			true,
-		},
-		"negative uReserve": {
-			osmomath.NewBigDec(100),
-			osmomath.NewBigDec(100),
-			// represents a 4-asset pool with 100 in each reserve
-			osmomath.NewBigDec(-200),
-			osmomath.NewBigDec(20000),
-			osmomath.NewBigDec(1),
-			true,
-		},
-		"negative sumSquares": {
-			osmomath.NewBigDec(100),
-			osmomath.NewBigDec(100),
-			// represents a 4-asset pool with 100 in each reserve
-			osmomath.NewBigDec(200),
-			osmomath.NewBigDec(-20000),
-			osmomath.NewBigDec(1),
-			true,
-		},
-		"negative yIn": {
-			osmomath.NewBigDec(100),
-			osmomath.NewBigDec(100),
-			// represents a 4-asset pool with 100 in each reserve
-			osmomath.NewBigDec(200),
-			osmomath.NewBigDec(-20000),
-			osmomath.NewBigDec(1),
-			true,
-		},
+	for name, test := range tests {
+		t.Run(name, func(t *testing.T) {
+			// system under test
+			sut := func() {
+				uReserve := calcUReserve(test.remReserves)
+				wSumSquares := calcWSumSquares(test.remReserves)
+
+				// using multi-asset cfmm
+				k2 := cfmmConstantMulti(test.xReserve, test.yReserve, uReserve, wSumSquares)
+				xOut2 := solveCfmm(test.xReserve, test.yReserve, test.remReserves, test.yIn)
+				k3 := cfmmConstantMulti(test.xReserve.Sub(xOut2), test.yReserve.Add(test.yIn), uReserve, wSumSquares)
+				osmomath.DecApproxEq(t, k2, k3, kErrTolerance)
+			}
+
+			osmoassert.ConditionalPanic(t, test.expectPanic, sut)
+		})
 	}
+}
+
+func TestCFMMInvariantMultiAssetsBinarySearch(t *testing.T) {
+	kErrTolerance := osmomath.OneDec()
+
+	tests := multiAssetCFMMTestCases
 
 	for name, test := range tests {
 		t.Run(name, func(t *testing.T) {
 			// system under test
 			sut := func() {
+				uReserve := calcUReserve(test.remReserves)
+				wSumSquares := calcWSumSquares(test.remReserves)
+
 				// using multi-asset cfmm
-				k2 := cfmmConstantMulti(test.xReserve, test.yReserve, test.uReserve, test.wSumSquares)
-				xOut2 := solveCfmmMulti(test.xReserve, test.yReserve, test.wSumSquares, test.yIn)
-				k3 := cfmmConstantMulti(test.xReserve.Sub(xOut2), test.yReserve.Add(test.yIn), test.uReserve, test.wSumSquares)
+				k2 := cfmmConstantMulti(test.xReserve, test.yReserve, uReserve, wSumSquares)
+				xOut2 := solveCFMMBinarySearchMulti(cfmmConstantMulti)(test.xReserve, test.yReserve, uReserve, wSumSquares, test.yIn)
+				k3 := cfmmConstantMulti(test.xReserve.Sub(xOut2), test.yReserve.Add(test.yIn), uReserve, wSumSquares)
 				osmomath.DecApproxEq(t, k2, k3, kErrTolerance)
 			}
 
@@ -425,3 +593,19 @@ func (suite *StableSwapTestSuite) Test_StableSwap_CalculateAmountOutAndIn_Invers
 		}
 	}
 }
+
+func calcUReserve(remReserves []osmomath.BigDec) osmomath.BigDec {
+	uReserve := osmomath.OneDec()
+	for _, assetReserve := range remReserves {
+		uReserve = uReserve.Mul(assetReserve)
+	}
+	return uReserve
+}
+
+func calcWSumSquares(remReserves []osmomath.BigDec) osmomath.BigDec {
+	wSumSquares := osmomath.ZeroDec()
+	for _, assetReserve := range remReserves {
+		wSumSquares = wSumSquares.Add(assetReserve.Mul(assetReserve))
+	}
+	return wSumSquares
+}