Use new derivation for f

Author: hoytech

docs/f.md

@@ -0,0 +1,42 @@
+## Implementation of `f`
+
+### Derivation
+
+Formula 15 from the whitepaper:
+
+    y0 + (px / py) * (x0 - x) * (c + (1 - c) * (x0 / x))
+
+Multiply second term by `x/x`:
+
+    y0 + (px / py) * (x0 - x) * ((c * x) + (1 - c) * x0) / x
+
+`c` is scaled by `1e18`:
+
+    y0 + (px / py) * (x0 - x) * ((c * x) + (1e18 - c) * x0) / (x * 1e18)
+
+Re-order division by `py`:
+
+    y0 + px * (x0 - x) * ((c * x) + (1e18 - c) * x0) / (x * 1e18) / py
+
+Use `mulDiv` to avoid intermediate overflow:
+
+    y0 + Math.mulDiv(px * (x0 - x), c * x + (1e18 - c) * x0, x * 1e18) / py
+
+Round up for both divisions (operation is distributive):
+
+    y0 + (Math.mulDiv(px * (x0 - x), c * x + (1e18 - c) * x0, x * 1e18, Math.Rounding.Ceil) + (py-1)) / py
+
+### Boundary Analysis
+
+Pre-conditions: x <= x0, 1 <= {px,py} <= 1e36, {x0,y0} <= type(uint112).max, c <= 1e18
+
+None of the computations for the arguments to `mulDiv` can overflow:
+
+* Arg 1: `px * (x0 - x)`
+  * Upper-bound: `1e36*(2**112 - 1) =~ 232 bits`
+* Arg 2: `c * x + (1e18 - c) * x0`
+  * Upper-bound: `1e18*(2**112 - 1)*2 =~ 173 bits`
+* Arg 3: `x * 1e18`
+  * Upper-bound: `1e18*(2**112 - 1) =~ 172 bits`
+
+If amounts/prices are large, and we travel too far down the curve, then `mulDiv` (or the subsequent `y0` addition) could overflow because its output value cannot be represented as a `uint256`. However, these output values would never be valid anyway, because they exceed `type(uint112).max`.

src/EulerSwap.sol

@@ -7,6 +7,7 @@ import {IUniswapV2Callee} from "./interfaces/IUniswapV2Callee.sol";
 import {IEulerSwap} from "./interfaces/IEulerSwap.sol";
 import {IAllowanceTransfer} from "permit2/src/interfaces/IAllowanceTransfer.sol";
 import {EVCUtil} from "evc/utils/EVCUtil.sol";
+import {Math} from "openzeppelin-contracts/utils/math/Math.sol";
 
 contract EulerSwap is IEulerSwap, EVCUtil {
     bytes32 public constant curve = keccak256("EulerSwap v1");
@@ -236,7 +237,8 @@ contract EulerSwap is IEulerSwap, EVCUtil {
     }
 
     /// @dev EulerSwap curve definition
-    function f(uint256 xt, uint256 px, uint256 py, uint256 x0, uint256 y0, uint256 c) internal pure returns (uint256) {
-        return y0 + px * 1e18 / py * (c * (2 * x0 - xt) / 1e18 + (1e18 - c) * x0 / 1e18 * x0 / xt - x0) / 1e18;
+    /// Pre-conditions: x <= x0, 1 <= {px,py} <= 1e36, {x0,y0} <= type(uint112).max, c <= 1e18
+    function f(uint256 x, uint256 px, uint256 py, uint256 x0, uint256 y0, uint256 c) public pure returns (uint256) {
+        return y0 + (Math.mulDiv(px * (x0 - x), c * x + (1e18 - c) * x0, x * 1e18, Math.Rounding.Ceil) + (py - 1)) / py;
     }
 }

test/EulerSwapTest.t.sol

@@ -175,4 +175,15 @@ contract EulerSwapTest is EulerSwapTestBase {
             assertGe(getHolderNAV(), origNAV);
         }
     }
+
+    function test_fFuncOverflow(uint256 xt, uint256 px, uint256 py, uint256 x0, uint256 y0, uint256 c) public view {
+        x0 = bound(x0, 1, type(uint112).max);
+        y0 = bound(y0, 0, type(uint112).max);
+        xt = bound(xt, 1 + x0 / 1e6, x0); // from 1 millionth of the reserves left up
+        px = bound(px, 1, 1e36);
+        py = bound(py, 1, 1e36);
+        c = bound(c, 1, 1e18);
+
+        eulerSwap.f(xt, px, py, x0, y0, c);
+    }
 }