Fix rounding of 2CLP pools

pkg/pool-gyro/contracts/Gyro2CLPPool.sol

@@ -90,9 +90,10 @@ contract Gyro2CLPPool is IGyro2CLPPool, BalancerPoolToken {
         // New invariant
         invariant = invariant.mulUp(invariantRatio);
         uint256 squareNewInv = invariant * invariant;
-        // L / sqrt(beta)
+
+        // L / sqrt(beta), rounded down to maximize newBalance.
         uint256 a = invariant.divDown(sqrtBeta);
-        // L * sqrt(alpha)
+        // L * sqrt(alpha), rounded down to maximize newBalance (b is in the denominator).
         uint256 b = invariant.mulDown(sqrtAlpha);
 
         if (tokenInIndex == 0) {
@@ -167,7 +168,7 @@ contract Gyro2CLPPool is IGyro2CLPPool, BalancerPoolToken {
 
         uint256 currentInvariant = Gyro2CLPMath.calculateInvariant(balances, sqrtAlpha, sqrtBeta, rounding);
 
-        uint256[2] memory virtualOffsets = _calculateVirtualOffsets(currentInvariant, sqrtAlpha, sqrtBeta);
+        uint256[2] memory virtualOffsets = _calculateVirtualOffsets(currentInvariant, sqrtAlpha, sqrtBeta, rounding);
 
         virtualBalanceIn = tokenInIsToken0 ? virtualOffsets[0] : virtualOffsets[1];
         virtualBalanceOut = tokenInIsToken0 ? virtualOffsets[1] : virtualOffsets[0];
@@ -177,10 +178,11 @@ contract Gyro2CLPPool is IGyro2CLPPool, BalancerPoolToken {
     function _calculateVirtualOffsets(
         uint256 invariant,
         uint256 sqrtAlpha,
-        uint256 sqrtBeta
+        uint256 sqrtBeta,
+        Rounding rounding
     ) internal view virtual returns (uint256[2] memory virtualOffsets) {
-        virtualOffsets[0] = Gyro2CLPMath.calculateVirtualParameter0(invariant, sqrtBeta);
-        virtualOffsets[1] = Gyro2CLPMath.calculateVirtualParameter1(invariant, sqrtAlpha);
+        virtualOffsets[0] = Gyro2CLPMath.calculateVirtualParameter0(invariant, sqrtBeta, rounding);
+        virtualOffsets[1] = Gyro2CLPMath.calculateVirtualParameter1(invariant, sqrtAlpha, rounding);
         return virtualOffsets;
     }
 

pkg/pool-gyro/contracts/lib/Gyro2CLPMath.sol

@@ -45,7 +45,7 @@ library Gyro2CLPMath {
             rounding
         );
 
-        return calculateQuadratic(a, mb, bSquare, mc);
+        return calculateQuadratic(a, mb, bSquare, mc, rounding);
     }
 
     /**
@@ -116,17 +116,26 @@ library Gyro2CLPMath {
         uint256 a,
         uint256 mb,
         uint256 bSquare, // b^2 can be calculated separately with more precision
-        uint256 mc
+        uint256 mc,
+        Rounding rounding
     ) internal pure returns (uint256 invariant) {
-        uint256 denominator = a.mulUp(2 * FixedPoint.ONE);
+        function(uint256, uint256) pure returns (uint256) _mulUpOrDown = rounding == Rounding.ROUND_DOWN
+            ? FixedPoint.mulDown
+            : FixedPoint.mulUp;
+
+        function(uint256, uint256) pure returns (uint256) _divUpOrDown = rounding == Rounding.ROUND_DOWN
+            ? FixedPoint.divDown
+            : FixedPoint.divUp;
+
+        uint256 denominator = 2 * a;
         // Order multiplications for fixed point precision.
-        uint256 addTerm = (mc.mulDown(4 * FixedPoint.ONE)).mulDown(a);
+        uint256 addTerm = _mulUpOrDown(4 * mc, a);
         // The minus sign in the radicand cancels out in this special case.
         uint256 radicand = bSquare + addTerm;
         uint256 sqrResult = GyroPoolMath.sqrt(radicand, 5);
         // The minus sign in the numerator cancels out in this special case.
         uint256 numerator = mb + sqrResult;
-        invariant = numerator.divDown(denominator);
+        invariant = _divUpOrDown(numerator, denominator);
     }
 
     /**
@@ -224,13 +233,21 @@ library Gyro2CLPMath {
     }
 
     /// @dev Calculate the virtual offset `a` for reserves `x`, as in (x+a)*(y+b)=L^2.
-    function calculateVirtualParameter0(uint256 invariant, uint256 _sqrtBeta) internal pure returns (uint256) {
-        return invariant.divDown(_sqrtBeta);
+    function calculateVirtualParameter0(
+        uint256 invariant,
+        uint256 _sqrtBeta,
+        Rounding rounding
+    ) internal pure returns (uint256) {
+        return rounding == Rounding.ROUND_DOWN ? invariant.divDown(_sqrtBeta) : invariant.divUp(_sqrtBeta);
     }
 
     /// @dev Calculate the virtual offset `b` for reserves `y`, as in (x+a)*(y+b)=L^2.
-    function calculateVirtualParameter1(uint256 invariant, uint256 _sqrtAlpha) internal pure returns (uint256) {
-        return invariant.mulDown(_sqrtAlpha);
+    function calculateVirtualParameter1(
+        uint256 invariant,
+        uint256 _sqrtAlpha,
+        Rounding rounding
+    ) internal pure returns (uint256) {
+        return rounding == Rounding.ROUND_DOWN ? invariant.mulDown(_sqrtAlpha) : invariant.mulUp(_sqrtAlpha);
     }
 
     /**

pkg/pool-gyro/test/foundry/Gyro2CLPMathRounding.sol

@@ -0,0 +1,100 @@
+// SPDX-License-Identifier: GPL-3.0-or-later
+
+pragma solidity ^0.8.24;
+
+import "forge-std/Test.sol";
+
+import { Rounding } from "@balancer-labs/v3-interfaces/contracts/vault/VaultTypes.sol";
+
+import { ArrayHelpers } from "@balancer-labs/v3-solidity-utils/contracts/test/ArrayHelpers.sol";
+
+import "../../contracts/lib/Gyro2CLPMath.sol";
+
+struct QuadraticTerms {
+    uint256 a;
+    uint256 mb;
+    uint256 bSquare;
+    uint256 mc;
+}
+
+contract Gyro2CLPMathRoundingTest is Test {
+    using ArrayHelpers for *;
+
+    uint256 internal constant MIN_DIFF_ALPHA_BETA = 100;
+
+    uint256 internal constant MIN_SQRT_ALPHA = 0.01e18;
+    uint256 internal constant MAX_SQRT_ALPHA = 10e18 - MIN_DIFF_ALPHA_BETA;
+    // Make sqrtBeta 0.5% higher than sqrtAlpha
+    uint256 internal constant MIN_SQRT_BETA = MIN_SQRT_ALPHA;
+    uint256 internal constant MAX_SQRT_BETA = MAX_SQRT_ALPHA + MIN_DIFF_ALPHA_BETA;
+
+    function testCalculateQuadraticTermsRounding__Fuzz(
+        uint256[2] memory balances,
+        uint256 sqrtAlpha,
+        uint256 sqrtBeta
+    ) public pure {
+        balances[0] = bound(balances[0], 1e16, 1e8 * 1e18);
+        balances[1] = bound(balances[1], 1e16, 1e8 * 1e18);
+        sqrtAlpha = bound(sqrtAlpha, MIN_SQRT_ALPHA, MAX_SQRT_ALPHA);
+        sqrtBeta = bound(sqrtBeta, sqrtAlpha + MIN_DIFF_ALPHA_BETA, MAX_SQRT_BETA);
+
+        QuadraticTerms memory qTermsDown;
+        QuadraticTerms memory qTermsUp;
+
+        (qTermsDown.a, qTermsDown.mb, qTermsDown.bSquare, qTermsDown.mc) = Gyro2CLPMath.calculateQuadraticTerms(
+            balances.toMemoryArray(),
+            sqrtAlpha,
+            sqrtBeta,
+            Rounding.ROUND_DOWN
+        );
+        (qTermsUp.a, qTermsUp.mb, qTermsUp.bSquare, qTermsUp.mc) = Gyro2CLPMath.calculateQuadraticTerms(
+            balances.toMemoryArray(),
+            sqrtAlpha,
+            sqrtBeta,
+            Rounding.ROUND_UP
+        );
+
+        assertLe(qTermsUp.a, qTermsDown.a, "Wrong rounding result (a)");
+        assertGe(qTermsUp.mb, qTermsDown.mb, "Wrong rounding result (mb)");
+        assertGe(qTermsUp.bSquare, qTermsDown.bSquare, "Wrong rounding result (bSquare)");
+        assertGe(qTermsUp.mc, qTermsDown.mc, "Wrong rounding result (mc)");
+    }
+
+    function testCalculateQuadraticRounding__Fuzz(uint256 a, uint256 mb, uint256 mc) public pure {
+        a = bound(a, 1, 1e18); // 0 < a < FP(1)
+        mb = bound(mb, 1, 1e8 * 1e18);
+        uint256 bSquare = mb * mb; // This is an approximation just to unit fuzz calculate quadratic.
+        mc = bound(mc, 1, 1e8 * 1e18);
+
+        uint256 quadraticRoundDown = Gyro2CLPMath.calculateQuadratic(a, mb, bSquare, mc, Rounding.ROUND_DOWN);
+        uint256 quadraticRoundUp = Gyro2CLPMath.calculateQuadratic(a, mb, bSquare, mc, Rounding.ROUND_UP);
+
+        assertGe(quadraticRoundUp, quadraticRoundDown, "Wrong rounding result");
+    }
+
+    function testComputeInvariantRounding__Fuzz(
+        uint256[2] memory balances,
+        uint256 sqrtAlpha,
+        uint256 sqrtBeta
+    ) public pure {
+        balances[0] = bound(balances[0], 1e16, 1e8 * 1e18);
+        balances[1] = bound(balances[1], 1e16, 1e8 * 1e18);
+        sqrtAlpha = bound(sqrtAlpha, MIN_SQRT_ALPHA, MAX_SQRT_ALPHA);
+        sqrtBeta = bound(sqrtBeta, sqrtAlpha + MIN_DIFF_ALPHA_BETA, MAX_SQRT_BETA);
+
+        uint256 invariantDown = Gyro2CLPMath.calculateInvariant(
+            balances.toMemoryArray(),
+            sqrtAlpha,
+            sqrtBeta,
+            Rounding.ROUND_DOWN
+        );
+        uint256 invariantUp = Gyro2CLPMath.calculateInvariant(
+            balances.toMemoryArray(),
+            sqrtAlpha,
+            sqrtBeta,
+            Rounding.ROUND_UP
+        );
+
+        assertGe(invariantUp, invariantDown, "Wrong rounding result");
+    }
+}