Add 512bits add and mult operations (OpenZeppelin#5035)

Signed-off-by: Hadrien Croubois <hadrien.croubois@gmail.com>

Author: Amxx

.changeset/blue-nails-give.md

@@ -0,0 +1,5 @@
+---
+'openzeppelin-solidity': minor
+---
+
+`Math`: Add `add512`, `mul512` and `mulShr`.

contracts/utils/math/Math.sol

@@ -17,6 +17,34 @@ library Math {
         Expand // Away from zero
     }
 
+    /**
+     * @dev Return the 512-bit addition of two uint256.
+     *
+     * The result is stored in two 256 variables such that sum = high * 2²⁵⁶ + low.
+     */
+    function add512(uint256 a, uint256 b) internal pure returns (uint256 high, uint256 low) {
+        assembly ("memory-safe") {
+            low := add(a, b)
+            high := lt(low, a)
+        }
+    }
+
+    /**
+     * @dev Return the 512-bit multiplication of two uint256.
+     *
+     * The result is stored in two 256 variables such that product = high * 2²⁵⁶ + low.
+     */
+    function mul512(uint256 a, uint256 b) internal pure returns (uint256 high, uint256 low) {
+        // 512-bit multiply [high low] = x * y. Compute the product mod 2²⁵⁶ and mod 2²⁵⁶ - 1, then use
+        // the Chinese Remainder Theorem to reconstruct the 512 bit result. The result is stored in two 256
+        // variables such that product = high * 2²⁵⁶ + low.
+        assembly ("memory-safe") {
+            let mm := mulmod(a, b, not(0))
+            low := mul(a, b)
+            high := sub(sub(mm, low), lt(mm, low))
+        }
+    }
+
     /**
      * @dev Returns the addition of two unsigned integers, with an success flag (no overflow).
      */
@@ -143,42 +171,34 @@ library Math {
      */
     function mulDiv(uint256 x, uint256 y, uint256 denominator) internal pure returns (uint256 result) {
         unchecked {
-            // 512-bit multiply [prod1 prod0] = x * y. Compute the product mod 2²⁵⁶ and mod 2²⁵⁶ - 1, then use
-            // the Chinese Remainder Theorem to reconstruct the 512 bit result. The result is stored in two 256
-            // variables such that product = prod1 * 2²⁵⁶ + prod0.
-            uint256 prod0 = x * y; // Least significant 256 bits of the product
-            uint256 prod1; // Most significant 256 bits of the product
-            assembly {
-                let mm := mulmod(x, y, not(0))
-                prod1 := sub(sub(mm, prod0), lt(mm, prod0))
-            }
+            (uint256 high, uint256 low) = mul512(x, y);
 
             // Handle non-overflow cases, 256 by 256 division.
-            if (prod1 == 0) {
+            if (high == 0) {
                 // Solidity will revert if denominator == 0, unlike the div opcode on its own.
                 // The surrounding unchecked block does not change this fact.
                 // See https://docs.soliditylang.org/en/latest/control-structures.html#checked-or-unchecked-arithmetic.
-                return prod0 / denominator;
+                return low / denominator;
             }
 
             // Make sure the result is less than 2²⁵⁶. Also prevents denominator == 0.
-            if (denominator <= prod1) {
+            if (denominator <= high) {
                 Panic.panic(ternary(denominator == 0, Panic.DIVISION_BY_ZERO, Panic.UNDER_OVERFLOW));
             }
 
             ///////////////////////////////////////////////
             // 512 by 256 division.
             ///////////////////////////////////////////////
 
-            // Make division exact by subtracting the remainder from [prod1 prod0].
+            // Make division exact by subtracting the remainder from [high low].
             uint256 remainder;
             assembly {
                 // Compute remainder using mulmod.
                 remainder := mulmod(x, y, denominator)
 
                 // Subtract 256 bit number from 512 bit number.
-                prod1 := sub(prod1, gt(remainder, prod0))
-                prod0 := sub(prod0, remainder)
+                high := sub(high, gt(remainder, low))
+                low := sub(low, remainder)
             }
 
             // Factor powers of two out of denominator and compute largest power of two divisor of denominator.
@@ -189,15 +209,15 @@ library Math {
                 // Divide denominator by twos.
                 denominator := div(denominator, twos)
 
-                // Divide [prod1 prod0] by twos.
-                prod0 := div(prod0, twos)
+                // Divide [high low] by twos.
+                low := div(low, twos)
 
                 // Flip twos such that it is 2²⁵⁶ / twos. If twos is zero, then it becomes one.
                 twos := add(div(sub(0, twos), twos), 1)
             }
 
-            // Shift in bits from prod1 into prod0.
-            prod0 |= prod1 * twos;
+            // Shift in bits from high into low.
+            low |= high * twos;
 
             // Invert denominator mod 2²⁵⁶. Now that denominator is an odd number, it has an inverse modulo 2²⁵⁶ such
             // that denominator * inv ≡ 1 mod 2²⁵⁶. Compute the inverse by starting with a seed that is correct for
@@ -215,9 +235,9 @@ library Math {
 
             // Because the division is now exact we can divide by multiplying with the modular inverse of denominator.
             // This will give us the correct result modulo 2²⁵⁶. Since the preconditions guarantee that the outcome is
-            // less than 2²⁵⁶, this is the final result. We don't need to compute the high bits of the result and prod1
+            // less than 2²⁵⁶, this is the final result. We don't need to compute the high bits of the result and high
             // is no longer required.
-            result = prod0 * inverse;
+            result = low * inverse;
             return result;
         }
     }
@@ -229,6 +249,26 @@ library Math {
         return mulDiv(x, y, denominator) + SafeCast.toUint(unsignedRoundsUp(rounding) && mulmod(x, y, denominator) > 0);
     }
 
+    /**
+     * @dev Calculates floor(x * y >> n) with full precision. Throws if result overflows a uint256.
+     */
+    function mulShr(uint256 x, uint256 y, uint8 n) internal pure returns (uint256 result) {
+        unchecked {
+            (uint256 high, uint256 low) = mul512(x, y);
+            if (high >= 1 << n) {
+                Panic.panic(Panic.UNDER_OVERFLOW);
+            }
+            return (high << (256 - n)) | (low >> n);
+        }
+    }
+
+    /**
+     * @dev Calculates x * y >> n with full precision, following the selected rounding direction.
+     */
+    function mulShr(uint256 x, uint256 y, uint8 n, Rounding rounding) internal pure returns (uint256) {
+        return mulShr(x, y, n) + SafeCast.toUint(unsignedRoundsUp(rounding) && mulmod(x, y, 1 << n) > 0);
+    }
+
     /**
      * @dev Calculate the modular multiplicative inverse of a number in Z/nZ.
      *

test/helpers/enums.js

@@ -1,12 +1,14 @@
-function Enum(...options) {
-  return Object.fromEntries(options.map((key, i) => [key, BigInt(i)]));
-}
+const { ethers } = require('ethers');
+
+const Enum = (...options) => Object.fromEntries(options.map((key, i) => [key, BigInt(i)]));
+const EnumTyped = (...options) => Object.fromEntries(options.map((key, i) => [key, ethers.Typed.uint8(i)]));
 
 module.exports = {
   Enum,
+  EnumTyped,
   ProposalState: Enum('Pending', 'Active', 'Canceled', 'Defeated', 'Succeeded', 'Queued', 'Expired', 'Executed'),
   VoteType: Object.assign(Enum('Against', 'For', 'Abstain'), { Parameters: 255n }),
-  Rounding: Enum('Floor', 'Ceil', 'Trunc', 'Expand'),
+  Rounding: EnumTyped('Floor', 'Ceil', 'Trunc', 'Expand'),
   OperationState: Enum('Unset', 'Waiting', 'Ready', 'Done'),
-  RevertType: Enum('None', 'RevertWithoutMessage', 'RevertWithMessage', 'RevertWithCustomError', 'Panic'),
+  RevertType: EnumTyped('None', 'RevertWithoutMessage', 'RevertWithMessage', 'RevertWithCustomError', 'Panic'),
 };

test/utils/math/Math.t.sol

@@ -11,6 +11,48 @@ contract MathTest is Test {
         assertEq(Math.ternary(f, a, b), f ? a : b);
     }
 
+    // ADD512 & MUL512
+    function testAdd512(uint256 a, uint256 b) public pure {
+        (uint256 high, uint256 low) = Math.add512(a, b);
+
+        // test against tryAdd
+        (bool success, uint256 result) = Math.tryAdd(a, b);
+        if (success) {
+            assertEq(high, 0);
+            assertEq(low, result);
+        } else {
+            assertEq(high, 1);
+        }
+
+        // test against unchecked
+        unchecked {
+            assertEq(low, a + b); // unchecked allow overflow
+        }
+    }
+
+    function testMul512(uint256 a, uint256 b) public pure {
+        (uint256 high, uint256 low) = Math.mul512(a, b);
+
+        // test against tryMul
+        (bool success, uint256 result) = Math.tryMul(a, b);
+        if (success) {
+            assertEq(high, 0);
+            assertEq(low, result);
+        } else {
+            assertGt(high, 0);
+        }
+
+        // test against unchecked
+        unchecked {
+            assertEq(low, a * b); // unchecked allow overflow
+        }
+
+        // test against alternative method
+        (uint256 _high, uint256 _low) = _mulKaratsuba(a, b);
+        assertEq(high, _high);
+        assertEq(low, _low);
+    }
+
     // MIN & MAX
     function testSymbolicMinMax(uint256 a, uint256 b) public pure {
         assertEq(Math.min(a, b), a < b ? a : b);
@@ -184,7 +226,7 @@ contract MathTest is Test {
     // MULDIV
     function testMulDiv(uint256 x, uint256 y, uint256 d) public pure {
         // Full precision for x * y
-        (uint256 xyHi, uint256 xyLo) = _mulHighLow(x, y);
+        (uint256 xyHi, uint256 xyLo) = Math.mul512(x, y);
 
         // Assume result won't overflow (see {testMulDivDomain})
         // This also checks that `d` is positive
@@ -194,9 +236,9 @@ contract MathTest is Test {
         uint256 q = Math.mulDiv(x, y, d);
 
         // Full precision for q * d
-        (uint256 qdHi, uint256 qdLo) = _mulHighLow(q, d);
+        (uint256 qdHi, uint256 qdLo) = Math.mul512(q, d);
         // Add remainder of x * y / d (computed as rem = (x * y % d))
-        (uint256 qdRemLo, uint256 c) = _addCarry(qdLo, mulmod(x, y, d));
+        (uint256 c, uint256 qdRemLo) = Math.add512(qdLo, mulmod(x, y, d));
         uint256 qdRemHi = qdHi + c;
 
         // Full precision check that x * y = q * d + rem
@@ -206,7 +248,7 @@ contract MathTest is Test {
 
     /// forge-config: default.allow_internal_expect_revert = true
     function testMulDivDomain(uint256 x, uint256 y, uint256 d) public {
-        (uint256 xyHi, ) = _mulHighLow(x, y);
+        (uint256 xyHi, ) = Math.mul512(x, y);
 
         // Violate {testMulDiv} assumption (covers d is 0 and result overflow)
         vm.assume(xyHi >= d);
@@ -266,26 +308,13 @@ contract MathTest is Test {
         }
     }
 
-    function _nativeModExp(uint256 b, uint256 e, uint256 m) private pure returns (uint256) {
-        if (m == 1) return 0;
-        uint256 r = 1;
-        while (e > 0) {
-            if (e % 2 > 0) {
-                r = mulmod(r, b, m);
-            }
-            b = mulmod(b, b, m);
-            e >>= 1;
-        }
-        return r;
-    }
-
     // Helpers
     function _asRounding(uint8 r) private pure returns (Math.Rounding) {
         vm.assume(r < uint8(type(Math.Rounding).max));
         return Math.Rounding(r);
     }
 
-    function _mulHighLow(uint256 x, uint256 y) private pure returns (uint256 high, uint256 low) {
+    function _mulKaratsuba(uint256 x, uint256 y) private pure returns (uint256 high, uint256 low) {
         (uint256 x0, uint256 x1) = (x & type(uint128).max, x >> 128);
         (uint256 y0, uint256 y1) = (y & type(uint128).max, y >> 128);
 
@@ -305,10 +334,16 @@ contract MathTest is Test {
         }
     }
 
-    function _addCarry(uint256 x, uint256 y) private pure returns (uint256 res, uint256 carry) {
-        unchecked {
-            res = x + y;
+    function _nativeModExp(uint256 b, uint256 e, uint256 m) private pure returns (uint256) {
+        if (m == 1) return 0;
+        uint256 r = 1;
+        while (e > 0) {
+            if (e % 2 > 0) {
+                r = mulmod(r, b, m);
+            }
+            b = mulmod(b, b, m);
+            e >>= 1;
         }
-        carry = res < x ? 1 : 0;
+        return r;
     }
 }