Merge remote-tracking branch 'origin/mega-refactor-curve-gas'

Author: hoytech

src/libraries/CurveLib.sol

@@ -39,40 +39,50 @@ library CurveLib {
         }
     }
 
+    /// @dev EulerSwap inverse function definition
+    /// Pre-conditions: 0 < x <= x0, 1 <= {px,py} <= 1e36, {x0,y0} <= type(uint112).max, c <= 1e18
     function fInverse(uint256 y, uint256 px, uint256 py, uint256 x0, uint256 y0, uint256 c)
         internal
         pure
         returns (uint256)
     {
         // components of quadratic equation
-        int256 B = int256((py * (y - y0) + (px - 1)) / px) - (2 * int256(c) - int256(1e18)) * int256(x0) / 1e18;
+        int256 B;
         uint256 C;
         uint256 fourAC;
-        if (x0 < 1e18) {
-            C = ((1e18 - c) * x0 * x0 + (1e18 - 1)) / 1e18; // upper bound of 1e28 for x0 means this is safe
-            fourAC = Math.mulDiv(4 * c, C, 1e18, Math.Rounding.Ceil);
-        } else {
-            C = Math.mulDiv((1e18 - c), x0 * x0, 1e36, Math.Rounding.Ceil); // upper bound of 1e28 for x0 means this is safe
-            fourAC = Math.mulDiv(4 * c, C, 1, Math.Rounding.Ceil);
+        unchecked {
+            B = int256((py * (y - y0) + (px - 1)) / px) - (2 * int256(c) - int256(1e18)) * int256(x0) / 1e18;            
+            if (x0 >= 1e18) {
+                // if x0 >= 1, scale as normal
+                C = Math.mulDiv((1e18 - c), x0 * x0, 1e36, Math.Rounding.Ceil);
+                fourAC = 4 * c * C;
+            } else {
+                // if x0 < 1, then numbers get very small, so decrease scale to 1e18 to increase precision later
+                C = ((1e18 - c) * x0 * x0 + (1e18 - 1)) / 1e18;
+                fourAC = Math.mulDiv(4 * c, C, 1e18, Math.Rounding.Ceil);
+            }
         }
-
-        // solve for the square root
-        uint256 absB = abs(B);
+        
+        uint256 absB = uint256(B >= 0 ? B : -B);
         uint256 squaredB;
         uint256 discriminant;
         uint256 sqrt;
-        if (absB > 1e33) {
+        if (absB < 1e36) {
+            // safe to use naive squaring
+            unchecked {
+                squaredB = absB * absB;
+                discriminant = squaredB + fourAC; // keep in 1e36 scale for increased precision ahead of sqrt
+                sqrt = Math.sqrt(discriminant); // drop back to 1e18 scale
+                sqrt = (sqrt * sqrt < discriminant) ? sqrt + 1 : sqrt;
+            }                    
+        } else {
+            // use scaled, overflow-safe path
             uint256 scale = computeScale(absB);
             squaredB = Math.mulDiv(absB / scale, absB, scale, Math.Rounding.Ceil);
             discriminant = squaredB + fourAC / (scale * scale);
             sqrt = Math.sqrt(discriminant);
             sqrt = (sqrt * sqrt < discriminant) ? sqrt + 1 : sqrt;
             sqrt = sqrt * scale;
-        } else {
-            squaredB = Math.mulDiv(absB, absB, 1, Math.Rounding.Ceil);
-            discriminant = squaredB + fourAC; // keep in 1e36 scale for increased precision ahead of sqrt
-            sqrt = Math.sqrt(discriminant); // drop back to 1e18 scale
-            sqrt = (sqrt * sqrt < discriminant) ? sqrt + 1 : sqrt;
         }
 
         uint256 x;
@@ -90,25 +100,42 @@ library CurveLib {
     }
 
     function computeScale(uint256 x) internal pure returns (uint256 scale) {
+        // calculate number of bits in x
         uint256 bits = 0;
-        uint256 tmp = x;
-
-        while (tmp > 0) {
-            tmp >>= 1;
+        while (x > 0) {
+            x >>= 1;
             bits++;
         }
 
-        // absB * absB must be <= 2^256 ⇒ bits(B) ≤ 128
+        // 2^excessBits is how much we need to scale down to prevent overflow when squaring x
         if (bits > 128) {
-            uint256 excessBits = bits - 128;
-            // 2^excessBits is how much we need to scale down to prevent overflow
+            uint256 excessBits = bits - 128;            
             scale = 1 << excessBits;
         } else {
             scale = 1;
         }
     }
 
-    function abs(int256 x) internal pure returns (uint256) {
-        return uint256(x >= 0 ? x : -x);
+    function binarySearch(IEulerSwap.Params memory p, uint256 newReserve1, uint256 xMin, uint256 xMax)
+        internal
+        pure
+        returns (uint256)
+    {
+        if (xMin < 1) {
+            xMin = 1;
+        }
+        while (xMin < xMax) {
+            uint256 xMid = (xMin + xMax) / 2;
+            uint256 fxMid = f(xMid, p.priceX, p.priceY, p.equilibriumReserve0, p.equilibriumReserve1, p.concentrationX);
+            if (newReserve1 >= fxMid) {
+                xMax = xMid;
+            } else {
+                xMin = xMid + 1;
+            }
+        }
+        if (newReserve1 < f(xMin, p.priceX, p.priceY, p.equilibriumReserve0, p.equilibriumReserve1, p.concentrationX)) {
+            xMin += 1;
+        }
+        return xMin;
     }
 }

test/CurveLib.t.sol

@@ -0,0 +1,89 @@
+// SPDX-License-Identifier: GPL-2.0-or-later
+pragma solidity ^0.8.24;
+
+import "forge-std/Test.sol";
+import "forge-std/console.sol";
+import {EulerSwapTestBase, EulerSwap, TestERC20} from "./EulerSwapTestBase.t.sol";
+import {IEulerSwap} from "../src/interfaces/IEulerSwap.sol";
+import {CurveLib} from "../src/CurveLib.sol";
+
+contract CurveLibTest is EulerSwapTestBase {
+    EulerSwap public eulerSwap;
+
+    function setUp() public virtual override {
+        super.setUp();
+    }
+
+    function testGas_fInverse() public pure {
+        // Set representative values within valid bounds
+        uint256 px = 1e18;
+        uint256 py = 1e18;
+        uint256 x0 = 1e14;
+        uint256 y0 = 1e14;
+        uint256 c = 1e18;
+
+        // Use CurveLib.f to get a valid y
+        uint256 x = 1e12;
+        uint256 y = CurveLib.f(x, px, py, x0, y0, c);
+
+        // Measure gas of fInverse
+        CurveLib.fInverse(y, px, py, x0, y0, c);
+    }
+
+    function test_fuzzfInverse(uint256 x, uint256 px, uint256 py, uint256 x0, uint256 y0, uint256 cx, uint256 cy)
+        public
+        view
+    {
+        // Params
+        px = 1e18;
+        py = bound(py, 1, 1e36);
+        x0 = bound(x0, 1e2, 1e28);
+        y0 = bound(y0, 0, 1e28);
+        cx = bound(cx, 1, 1e18);
+        cy = bound(cy, 1, 1e18);
+        console.log("px", px);
+        console.log("py", py);
+        console.log("x0", x0);
+        console.log("y0", y0);
+        console.log("cx", cx);
+        console.log("cy", cy);
+
+        IEulerSwap.Params memory p = IEulerSwap.Params({
+            vault0: address(0),
+            vault1: address(0),
+            eulerAccount: address(0),
+            equilibriumReserve0: uint112(x0),
+            equilibriumReserve1: uint112(y0),
+            priceX: px,
+            priceY: py,
+            concentrationX: cx,
+            concentrationY: cy,
+            fee: 0,
+            protocolFee: 0,
+            protocolFeeRecipient: address(0)
+        });
+
+        // Note without -2 in the max bound, f() sometimes fails when x gets too close to centre.
+        // Note small x values lead to large y-values, which causes problems for both f() and fInverse(), so we cap it here
+        x = bound(x, 1e2 - 3, x0 - 3);
+
+        uint256 y = CurveLib.f(x, px, py, x0, y0, cx);
+        console.log("y    ", y);
+        uint256 xCalc = CurveLib.fInverse(y, px, py, x0, y0, cx);
+        console.log("xCalc", xCalc);
+        uint256 yCalc = CurveLib.f(xCalc, px, py, x0, y0, cx);
+        uint256 xBin = CurveLib.binarySearch(p, y, 1, x0);
+        uint256 yBin = CurveLib.f(xBin, px, py, x0, y0, cx);
+        console.log("x    ", x);
+        console.log("xCalc", xCalc);
+        console.log("xBin ", xBin);
+        console.log("y    ", y);
+        console.log("yCalc", yCalc);
+        console.log("yBin ", yBin);
+
+        if (x < type(uint112).max && y < type(uint112).max) {
+            assert(CurveLib.verify(p, xCalc, y));
+            assert(int256(xCalc) - int256(xBin) <= 3 || int256(yCalc) - int256(yBin) <= 3); // suspect this is 2 wei error in fInverse() + 1 wei error in f()
+        }
+    }
+}