2025 08 16 no normal mul

src/lib/implementation/LibDecimalFloatImplementation.sol

@@ -94,42 +94,106 @@ library LibDecimalFloatImplementation {
         }
     }
 
-    /// Stack only implementation of `mul`.
-    function mul(int256 signedCoefficientA, int256 exponentA, int256 signedCoefficientB, int256 exponentB)
+    function absUnsignedSignedCoefficient(int256 signedCoefficient) internal pure returns (uint256) {
+        unchecked {
+            if (signedCoefficient < 0) {
+                if (signedCoefficient == type(int256).min) {
+                    return uint256(type(int256).max) + 1;
+                } else {
+                    return uint256(-signedCoefficient);
+                }
+            } else {
+                return uint256(signedCoefficient);
+            }
+        }
+    }
+
+    function unabsUnsignedMulOrDivLossy(int256 a, int256 b, uint256 signedCoefficientAbs, int256 exponent)
         internal
         pure
         returns (int256, int256)
     {
         unchecked {
-            // Unchecked mul the coefficients and add the exponents.
-            int256 signedCoefficient = signedCoefficientA * signedCoefficientB;
-
-            // Need to return early if the result is zero to avoid divide by
-            // zero in the overflow check.
-            if (signedCoefficient == 0) {
-                return (NORMALIZED_ZERO_SIGNED_COEFFICIENT, NORMALIZED_ZERO_EXPONENT);
+            // Need to minus the coefficient because a and b had different signs.
+            if ((a ^ b) < 0) {
+                if (signedCoefficientAbs > uint256(type(int256).max)) {
+                    if (signedCoefficientAbs == uint256(type(int256).max) + 1) {
+                        // Edge case where the absolute value is exactly
+                        // type(int256).min.
+                        return (type(int256).min, exponent);
+                    } else {
+                        return (-int256(signedCoefficientAbs / 10), exponent + 1);
+                    }
+                } else {
+                    return (-int256(signedCoefficientAbs), exponent);
+                }
+            } else {
+                if (signedCoefficientAbs > uint256(type(int256).max)) {
+                    return (int256(signedCoefficientAbs / 10), exponent + 1);
+                } else {
+                    return (int256(signedCoefficientAbs), exponent);
+                }
             }
+        }
+    }
 
-            int256 exponent = exponentA + exponentB;
+    /// Stack only implementation of `mul`.
+    function mul(int256 signedCoefficientA, int256 exponentA, int256 signedCoefficientB, int256 exponentB)
+        internal
+        pure
+        returns (int256 signedCoefficient, int256 exponent)
+    {
+        bool isZero;
+        assembly ("memory-safe") {
+            isZero := or(iszero(signedCoefficientA), iszero(signedCoefficientB))
+        }
+        if (isZero) {
+            // These sets are redundant as both are zero but this makes it
+            // clearer and more explicit.
+            signedCoefficient = MAXIMIZED_ZERO_SIGNED_COEFFICIENT;
+            exponent = MAXIMIZED_ZERO_EXPONENT;
+        } else {
+            exponent = exponentA + exponentB;
 
-            // No jumps to see if we overflowed.
-            bool didOverflow;
-            assembly ("memory-safe") {
-                didOverflow :=
-                    or(
-                        iszero(eq(sdiv(signedCoefficient, signedCoefficientA), signedCoefficientB)),
-                        iszero(eq(sub(exponent, exponentA), exponentB))
-                    )
-            }
-            // If we did overflow, normalize and try again. Normalized values
-            // cannot overflow, so this will always succeed, provided the
-            // exponents are not out of bounds.
-            if (didOverflow) {
-                (signedCoefficientA, exponentA) = normalize(signedCoefficientA, exponentA);
-                (signedCoefficientB, exponentB) = normalize(signedCoefficientB, exponentB);
-                return mul(signedCoefficientA, exponentA, signedCoefficientB, exponentB);
+            // mulDiv only works with unsigned integers, so get the absolute
+            // values of the coefficients.
+            uint256 signedCoefficientAAbs = absUnsignedSignedCoefficient(signedCoefficientA);
+            uint256 signedCoefficientBAbs = absUnsignedSignedCoefficient(signedCoefficientB);
+
+            (uint256 prod1,) = mul512(signedCoefficientAAbs, signedCoefficientBAbs);
+
+            uint256 adjustExponent = 0;
+            unchecked {
+                if (prod1 > 1e37) {
+                    prod1 /= 1e37;
+                    adjustExponent += 37;
+                }
+                if (prod1 > 1e18) {
+                    prod1 /= 1e18;
+                    adjustExponent += 18;
+                }
+                if (prod1 > 1e9) {
+                    prod1 /= 1e9;
+                    adjustExponent += 9;
+                }
+                if (prod1 > 1e4) {
+                    prod1 /= 1e4;
+                    adjustExponent += 4;
+                }
+                while (prod1 > 0) {
+                    prod1 /= 10;
+                    adjustExponent++;
+                }
             }
-            return (signedCoefficient, exponent);
+
+            exponent += int256(adjustExponent);
+
+            (signedCoefficient, exponent) = unabsUnsignedMulOrDivLossy(
+                signedCoefficientA,
+                signedCoefficientB,
+                mulDiv(signedCoefficientAAbs, signedCoefficientBAbs, uint256(10) ** adjustExponent),
+                exponent
+            );
         }
     }
 
@@ -186,42 +250,24 @@ library LibDecimalFloatImplementation {
     function div(int256 signedCoefficientA, int256 exponentA, int256 signedCoefficientB, int256 exponentB)
         internal
         pure
-        returns (int256, int256)
+        returns (int256 signedCoefficient, int256 exponent)
     {
-        uint256 scale = 1e76;
-        int256 adjustExponent = 76;
-        int256 signedCoefficient;
-
-        unchecked {
+        if (signedCoefficientA == 0) {
+            signedCoefficient = MAXIMIZED_ZERO_SIGNED_COEFFICIENT;
+            exponent = MAXIMIZED_ZERO_EXPONENT;
+        } else {
             // Move both coefficients into the e75/e76 range, so that the result
             // of division will not cause a mulDiv overflow.
             (signedCoefficientA, exponentA) = maximize(signedCoefficientA, exponentA);
             (signedCoefficientB, exponentB) = maximize(signedCoefficientB, exponentB);
 
             // mulDiv only works with unsigned integers, so get the absolute
             // values of the coefficients.
-            uint256 signedCoefficientAAbs;
-            if (signedCoefficientA > 0) {
-                signedCoefficientAAbs = uint256(signedCoefficientA);
-            } else if (signedCoefficientA < 0) {
-                if (signedCoefficientA == type(int256).min) {
-                    signedCoefficientAAbs = uint256(type(int256).max) + 1;
-                } else {
-                    signedCoefficientAAbs = uint256(-signedCoefficientA);
-                }
-            } else {
-                return (MAXIMIZED_ZERO_SIGNED_COEFFICIENT, MAXIMIZED_ZERO_EXPONENT);
-            }
-            uint256 signedCoefficientBAbs;
-            if (signedCoefficientB < 0) {
-                if (signedCoefficientB == type(int256).min) {
-                    signedCoefficientBAbs = uint256(type(int256).max) + 1;
-                } else {
-                    signedCoefficientBAbs = uint256(-signedCoefficientB);
-                }
-            } else {
-                signedCoefficientBAbs = uint256(signedCoefficientB);
-            }
+            uint256 signedCoefficientAAbs = absUnsignedSignedCoefficient(signedCoefficientA);
+            uint256 signedCoefficientBAbs = absUnsignedSignedCoefficient(signedCoefficientB);
+
+            uint256 scale = 1e76;
+            int256 adjustExponent = 76;
 
             // We are going to scale the numerator up by the largest power of ten
             // that is smaller than the denominator. This will always overflow
@@ -233,31 +279,35 @@ library LibDecimalFloatImplementation {
                 scale = 1e75;
                 adjustExponent = 75;
             }
-            uint256 signedCoefficientAbs = mulDiv(signedCoefficientAAbs, scale, signedCoefficientBAbs);
-            signedCoefficient = (signedCoefficientA ^ signedCoefficientB) < 0
-                ? -int256(signedCoefficientAbs)
-                : int256(signedCoefficientAbs);
+            exponent = exponentA - exponentB - adjustExponent;
+
+            (signedCoefficient, exponent) = unabsUnsignedMulOrDivLossy(
+                signedCoefficientA,
+                signedCoefficientB,
+                mulDiv(signedCoefficientAAbs, scale, signedCoefficientBAbs),
+                exponent
+            );
         }
+    }
 
-        // Keep the exponent calculation outside the unchecked block so that we
-        // don't silently under/overflow.
-        int256 exponent = exponentA - exponentB - adjustExponent;
-        return (signedCoefficient, exponent);
+    /// mul512 from Open Zeppelin.
+    /// Simply part of the original mulDiv function abstracted out for reuse
+    /// elsewhere.
+    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))
+        }
     }
 
     /// mulDiv as seen in Open Zeppelin, PRB Math, Solady, and other libraries.
     /// Credit to Remco Bloemen under MIT license: https://2π.com/21/muldiv
     function mulDiv(uint256 x, uint256 y, uint256 denominator) internal pure returns (uint256 result) {
-        // 512-bit multiply [prod1 prod0] = x * y. Compute the product mod 2^256 and mod 2^256 - 1, then use
-        // use the Chinese Remainder Theorem to reconstruct the 512-bit result. The result is stored in two 256
-        // variables such that product = prod1 * 2^256 + prod0.
-        uint256 prod0; // Least significant 256 bits of the product
-        uint256 prod1; // Most significant 256 bits of the product
-        assembly ("memory-safe") {
-            let mm := mulmod(x, y, not(0))
-            prod0 := mul(x, y)
-            prod1 := sub(sub(mm, prod0), lt(mm, prod0))
-        }
+        (uint256 prod1, uint256 prod0) = mul512(x, y);
 
         // Handle non-overflow cases, 256 by 256 division.
         if (prod1 == 0) {

test/lib/LibDecimalFloatSlow.sol

@@ -13,29 +13,37 @@ library LibDecimalFloatSlow {
         returns (int256, int256)
     {
         unchecked {
-            int256 signedCoefficient = signedCoefficientA * signedCoefficientB;
-            int256 exponent = exponentA + exponentB;
-
             // If the expected signed coefficient is 0 then everything is just
             // normalized 0.
-            if (signedCoefficient == 0) {
+            if (signedCoefficientA == 0 || signedCoefficientB == 0) {
                 return (0, 0);
             }
-            // If nothing overflowed then our expected outcome is correct.
-            else if (signedCoefficient / signedCoefficientA == signedCoefficientB && exponent - exponentA == exponentB)
-            {
-                return (signedCoefficient, exponent);
-            }
-            // If something overflowed then we have to normalize and try again.
-            else {
-                (signedCoefficientA, exponentA) = LibDecimalFloatImplementation.normalize(signedCoefficientA, exponentA);
-                (signedCoefficientB, exponentB) = LibDecimalFloatImplementation.normalize(signedCoefficientB, exponentB);
 
-                signedCoefficient = signedCoefficientA * signedCoefficientB;
-                exponent = exponentA + exponentB;
+            int256 exponent = exponentA + exponentB;
 
-                return (signedCoefficient, exponent);
+            uint256 signedCoefficientAAbs =
+                LibDecimalFloatImplementation.absUnsignedSignedCoefficient(signedCoefficientA);
+            uint256 signedCoefficientBAbs =
+                LibDecimalFloatImplementation.absUnsignedSignedCoefficient(signedCoefficientB);
+
+            (uint256 prod1,) = LibDecimalFloatImplementation.mul512(signedCoefficientAAbs, signedCoefficientBAbs);
+
+            uint256 adjustExponent = 0;
+            while (prod1 > 0) {
+                prod1 /= 10;
+                adjustExponent++;
             }
+
+            uint256 signedCoefficientAbs = LibDecimalFloatImplementation.mulDiv(
+                signedCoefficientAAbs, signedCoefficientBAbs, uint256(10) ** adjustExponent
+            );
+
+            exponent += int256(adjustExponent);
+            int256 signedCoefficient;
+            (signedCoefficient, exponent) = LibDecimalFloatImplementation.unabsUnsignedMulOrDivLossy(
+                signedCoefficientA, signedCoefficientB, signedCoefficientAbs, exponent
+            );
+            return (signedCoefficient, exponent);
         }
     }
 
@@ -183,4 +191,19 @@ library LibDecimalFloatSlow {
 
         return ltSlow(signedCoefficientA, exponentA, signedCoefficientB, exponentB);
     }
+
+    function maximizeSlow(int256 signedCoefficient, int256 exponent) internal pure returns (int256, int256) {
+        unchecked {
+            if (signedCoefficient == 0) {
+                return (0, 0);
+            }
+            int256 trySignedCoefficient = signedCoefficient * 10;
+            while (trySignedCoefficient / 10 == signedCoefficient) {
+                signedCoefficient = trySignedCoefficient;
+                exponent--;
+                trySignedCoefficient *= 10;
+            }
+            return (signedCoefficient, exponent);
+        }
+    }
 }

test/src/lib/LibDecimalFloat.pow.t.sol

@@ -47,6 +47,8 @@ contract LibDecimalFloatPowTest is LogTest {
             5e37, -38, 6e37, -36, 8.710801393728222996515679442508710801393728222996515679442508710801e66, -66 - 19
         );
         // Issues found in fuzzing from here.
+        // 99999 ^ 12182 = 8.853071703048649170130397094169464632911643045383977634639832230468640539353...e60910
+        // 8.853071703048649170130397094169464632911643045383977634639832230468640539353e75 e60910
         checkPow(99999, 0, 12182, 0, 1000, 60907);
         checkPow(1785215562, 0, 18, 0, 3388, 163);
     }

test/src/lib/implementation/LibDecimalFloatImplementation.maximize.t.sol

@@ -10,6 +10,7 @@ import {
     MAXIMIZED_ZERO_EXPONENT,
     MAXIMIZED_ZERO_SIGNED_COEFFICIENT
 } from "src/lib/implementation/LibDecimalFloatImplementation.sol";
+import {LibDecimalFloatSlow} from "test/lib/LibDecimalFloatSlow.sol";
 
 contract LibDecimalFloatImplementationMaximizeTest is Test {
     function isMaximized(int256 signedCoefficient, int256 exponent) internal pure returns (bool) {
@@ -79,4 +80,15 @@ contract LibDecimalFloatImplementationMaximizeTest is Test {
         assertEq(actualSignedCoefficient, maximizedSignedCoefficient);
         assertEq(actualExponent, maximizedExponent);
     }
+
+    /// Maximization against reference.
+    function testMaximizedReference(int256 signedCoefficient, int256 exponent) external pure {
+        exponent = bound(exponent, EXPONENT_MIN, EXPONENT_MAX);
+        (int256 actualSignedCoefficient, int256 actualExponent) =
+            LibDecimalFloatImplementation.maximize(signedCoefficient, exponent);
+        (int256 expectedSignedCoefficient, int256 expectedExponent) =
+            LibDecimalFloatSlow.maximizeSlow(signedCoefficient, exponent);
+        assertEq(actualSignedCoefficient, expectedSignedCoefficient);
+        assertEq(actualExponent, expectedExponent);
+    }
 }