full precision div

crates/float/src/lib.rs

@@ -1368,7 +1368,10 @@ mod tests {
         let zero = Float::parse("0".to_string()).unwrap();
         let err = (one / zero).unwrap_err();
 
-        assert!(matches!(err, FloatError::Revert(_)));
+        assert!(matches!(
+            err,
+            FloatError::DecimalFloat(DecimalFloatErrors::MulDivOverflow(_))
+        ));
     }
 
     #[test]

src/error/ErrDecimalFloat.sol

@@ -30,3 +30,6 @@ error LossyConversionFromFloat(int256 signedCoefficient, int256 exponent);
 
 /// @dev Thrown when attempting to exponentiate 0^b where b is negative.
 error ZeroNegativePower(Float b);
+
+/// @dev Thrown when mulDiv internal to division overflows.
+error MulDivOverflow(uint256 x, uint256 y, uint256 denominator);

src/lib/implementation/LibDecimalFloatImplementation.sol

@@ -1,7 +1,7 @@
 // SPDX-License-Identifier: CAL
 pragma solidity ^0.8.25;
 
-import {ExponentOverflow, Log10Negative, Log10Zero} from "../../error/ErrDecimalFloat.sol";
+import {ExponentOverflow, Log10Negative, Log10Zero, MulDivOverflow} from "../../error/ErrDecimalFloat.sol";
 import {
     LOG_TABLES,
     LOG_TABLES_SMALL,
@@ -188,13 +188,149 @@ library LibDecimalFloatImplementation {
         pure
         returns (int256, int256)
     {
+        uint256 scale = 1e76;
+        int256 adjustExponent = 76;
+        int256 signedCoefficient;
+
         unchecked {
+            // 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) = normalize(signedCoefficientB, exponentB);
+            (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);
+            }
 
-            int256 signedCoefficient = signedCoefficientA / signedCoefficientB;
-            int256 exponent = exponentA - exponentB;
-            return (signedCoefficient, exponent);
+            // We are going to scale the numerator up by the largest power of ten
+            // that is smaller than the denominator. This will always overflow
+            // internally to the mulDiv during the initial multiplication, in
+            // 512 bits, but will subsequently always be reduced back down to
+            // fit in 256 bits by the division of a denominator that is larger
+            // than the scale up.
+            if (signedCoefficientBAbs < scale) {
+                scale = 1e75;
+                adjustExponent = 75;
+            }
+            uint256 signedCoefficientAbs = mulDiv(signedCoefficientAAbs, scale, signedCoefficientBAbs);
+            signedCoefficient = (signedCoefficientA ^ signedCoefficientB) < 0
+                ? -int256(signedCoefficientAbs)
+                : int256(signedCoefficientAbs);
+        }
+
+        // Keep the exponent calculation outside the unchecked block so that we
+        // don't silently under/overflow.
+        int256 exponent = exponentA - exponentB - adjustExponent;
+        return (signedCoefficient, exponent);
+    }
+
+    /// 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))
+        }
+
+        // Handle non-overflow cases, 256 by 256 division.
+        if (prod1 == 0) {
+            unchecked {
+                return prod0 / denominator;
+            }
+        }
+
+        // Make sure the result is less than 2^256. Also prevents denominator == 0.
+        if (prod1 >= denominator) {
+            revert MulDivOverflow(x, y, denominator);
+        }
+
+        ////////////////////////////////////////////////////////////////////////////
+        // 512 by 256 division
+        ////////////////////////////////////////////////////////////////////////////
+
+        // Make division exact by subtracting the remainder from [prod1 prod0].
+        uint256 remainder;
+        assembly ("memory-safe") {
+            // Compute remainder using the mulmod Yul instruction.
+            remainder := mulmod(x, y, denominator)
+
+            // Subtract 256 bit number from 512-bit number.
+            prod1 := sub(prod1, gt(remainder, prod0))
+            prod0 := sub(prod0, remainder)
+        }
+
+        unchecked {
+            // Calculate the largest power of two divisor of the denominator using the unary operator ~. This operation cannot overflow
+            // because the denominator cannot be zero at this point in the function execution. The result is always >= 1.
+            // For more detail, see https://cs.stackexchange.com/q/138556/92363.
+            uint256 lpotdod = denominator & (~denominator + 1);
+            uint256 flippedLpotdod;
+
+            assembly ("memory-safe") {
+                // Factor powers of two out of denominator.
+                // slither-disable-next-line divide-before-multiply
+                denominator := div(denominator, lpotdod)
+
+                // Divide [prod1 prod0] by lpotdod.
+                // slither-disable-next-line divide-before-multiply
+                prod0 := div(prod0, lpotdod)
+
+                // Get the flipped value `2^256 / lpotdod`. If the `lpotdod` is zero, the flipped value is one.
+                // `sub(0, lpotdod)` produces the two's complement version of `lpotdod`, which is equivalent to flipping all the bits.
+                // However, `div` interprets this value as an unsigned value: https://ethereum.stackexchange.com/q/147168/24693
+                flippedLpotdod := add(div(sub(0, lpotdod), lpotdod), 1)
+            }
+
+            // Shift in bits from prod1 into prod0.
+            prod0 |= prod1 * flippedLpotdod;
+
+            // Invert denominator mod 2^256. Now that denominator is an odd number, it has an inverse modulo 2^256 such
+            // that denominator * inv = 1 mod 2^256. Compute the inverse by starting with a seed that is correct for
+            // four bits. That is, denominator * inv = 1 mod 2^4.
+            // slither-disable-next-line incorrect-exp
+            uint256 inverse = (3 * denominator) ^ 2;
+
+            // Use the Newton-Raphson iteration to improve the precision. Thanks to Hensel's lifting lemma, this also works
+            // in modular arithmetic, doubling the correct bits in each step.
+            inverse *= 2 - denominator * inverse; // inverse mod 2^8
+            inverse *= 2 - denominator * inverse; // inverse mod 2^16
+            inverse *= 2 - denominator * inverse; // inverse mod 2^32
+            inverse *= 2 - denominator * inverse; // inverse mod 2^64
+            inverse *= 2 - denominator * inverse; // inverse mod 2^128
+            inverse *= 2 - denominator * inverse; // inverse mod 2^256
+
+            // 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^256. Since the preconditions guarantee that the outcome is
+            // less than 2^256, this is the final result. We don't need to compute the high bits of the result and prod1
+            // is no longer required.
+            result = prod0 * inverse;
         }
     }
 
@@ -367,14 +503,9 @@ library LibDecimalFloatImplementation {
         return signedCoefficientA == signedCoefficientB;
     }
 
-    /// Inverts a float. Equivalent to `1 / x` with modest gas optimizations.
+    /// Inverts a float. Equivalent to `1 / x`.
     function inv(int256 signedCoefficient, int256 exponent) internal pure returns (int256, int256) {
-        (signedCoefficient, exponent) = normalize(signedCoefficient, exponent);
-
-        signedCoefficient = 1e76 / signedCoefficient;
-        exponent = -exponent - 76;
-
-        return (signedCoefficient, exponent);
+        return div(1e76, -76, signedCoefficient, exponent);
     }
 
     /// log10(x) for a float x.

test/lib/LibCommonResults.sol

@@ -1,5 +1,6 @@
 // SPDX-License-Identifier: CAL
 pragma solidity ^0.8.25;
 
-int256 constant ONES = 111111111111111111111111111111111111111;
-int256 constant THREES = 333333333333333333333333333333333333333;
+int256 constant ONES = 1111111111111111111111111111111111111111111111111111111111111111111111111111;
+int256 constant THREES_PACKED = 3333333333333333333333333333333333333333333333333333333333333333333;
+int256 constant THREES = 3333333333333333333333333333333333333333333333333333333333333333333333333333;

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

@@ -2,26 +2,26 @@
 pragma solidity =0.8.25;
 
 import {LibDecimalFloat, Float} from "src/lib/LibDecimalFloat.sol";
-import {THREES, ONES} from "../../lib/LibCommonResults.sol";
+import {THREES_PACKED, ONES} from "../../lib/LibCommonResults.sol";
 
 import {Test} from "forge-std/Test.sol";
 
 contract LibDecimalFloatMixedTest is Test {
     using LibDecimalFloat for Float;
 
     /// (1 / 3) * 555e18
-    function testDiv1Over3() external pure {
+    function testDiv1Over3Mixed() external pure {
         Float a = LibDecimalFloat.packLossless(1, 0);
         Float b = LibDecimalFloat.packLossless(3, 0);
         Float c = a.div(b);
         (int256 signedCoefficientDiv, int256 exponentDiv) = LibDecimalFloat.unpack(c);
-        assertEq(signedCoefficientDiv, THREES, "coefficient");
-        assertEq(exponentDiv, -39, "exponent");
+        assertEq(signedCoefficientDiv, THREES_PACKED, "coefficient");
+        assertEq(exponentDiv, -67, "exponent");
 
         Float d = c.mul(LibDecimalFloat.packLossless(555, 18));
         (int256 signedCoefficientMul, int256 exponentMul) = LibDecimalFloat.unpack(d);
 
-        assertEq(signedCoefficientMul, 184999999999999999999999999999999999999815);
-        assertEq(exponentMul, -21);
+        assertEq(signedCoefficientMul, 1849999999999999999999999999999999999999999999999999999999999999999);
+        assertEq(exponentMul, -46);
     }
 }

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

@@ -38,9 +38,15 @@ contract LibDecimalFloatPowTest is LogTest {
     }
 
     function testPows() external {
-        checkPow(5e37, -38, 3e37, -36, 9.32835820895522388059701492537313432835e38, -48);
-        checkPow(5e37, -38, 6e37, -36, 8.71080139372822299651567944250871080139e38, -57);
-        // // Issues found in fuzzing from here.
+        // 0.5 ^ 30 = 9.3132257e-10
+        checkPow(
+            5e37, -38, 3e37, -36, 9.328358208955223880597014925373134328358208955223880597014925373134e66, -66 - 10
+        );
+        // 0.5 ^ 60 = 8.6736174e-19
+        checkPow(
+            5e37, -38, 6e37, -36, 8.710801393728222996515679442508710801393728222996515679442508710801e66, -66 - 19
+        );
+        // Issues found in fuzzing from here.
         checkPow(99999, 0, 12182, 0, 1000, 60907);
         checkPow(1785215562, 0, 18, 0, 3388, 163);
     }

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

@@ -27,6 +27,12 @@ contract LibDecimalFloatPow10Test is LogTest {
             } else {
                 Float floatPower10 = this.pow10External(float);
                 (int256 signedCoefficientUnpacked, int256 exponentUnpacked) = floatPower10.unpack();
+
+                // Compensate for the implied pack and unpack.
+                (Float resultPacked, bool lossless) = LibDecimalFloat.packLossy(signedCoefficient, exponent);
+                (lossless);
+                (signedCoefficient, exponent) = resultPacked.unpack();
+
                 assertEq(signedCoefficient, signedCoefficientUnpacked);
                 assertEq(exponent, exponentUnpacked);
             }

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

@@ -40,12 +40,12 @@ contract LibDecimalFloatImplementationDivTest is Test {
 
     /// 1 / 3
     function testDiv1Over3() external pure {
-        checkDiv(1, 0, 3, 0, THREES, -39);
+        checkDiv(1, 0, 3, 0, THREES, -76);
     }
 
     /// - 1 / 3
     function testDivNegative1Over3() external pure {
-        checkDiv(-1, 0, 3, 0, -THREES, -39);
+        checkDiv(-1, 0, 3, 0, -THREES, -76);
     }
 
     /// 1 / 3 gas
@@ -56,41 +56,47 @@ contract LibDecimalFloatImplementationDivTest is Test {
 
     /// 1e18 / 3
     function testDiv1e18Over3() external pure {
-        checkDiv(1e18, 0, 3, 0, THREES, -21);
+        checkDiv(1e18, 0, 3, 0, THREES, -58);
     }
 
     /// 10,0 / 1e38,-37 == 1
     function testDivTenOverOOMs() external pure {
-        checkDiv(10, 0, 1e38, -37, 1e39, -39);
+        checkDiv(10, 0, 1e38, -37, 1e76, -76);
     }
 
     /// 1e38,-37 / 2,0 == 5
     function testDivOOMsOverTen() external pure {
-        checkDiv(1e38, -37, 2, 0, 5e38, -38);
+        checkDiv(1e38, -37, 2, 0, 5e75, -75);
     }
 
     /// 5e37,-37 / 2e37,-37 == 2.5
     function testDivOOMs5and2() external pure {
-        checkDiv(5e37, -37, 2e37, -37, 25e38, -39);
+        checkDiv(5e37, -37, 2e37, -37, 2.5e76, -76);
     }
 
     /// (1 / 9) / (1 / 3) == 0.333..
     function testDiv1Over9Over1Over3() external pure {
         // 1 / 9
         (int256 signedCoefficientA, int256 exponentA) = LibDecimalFloatImplementation.div(1, 0, 9, 0);
         assertEq(signedCoefficientA, ONES);
-        assertEq(exponentA, -39);
+        assertEq(exponentA, -76);
 
         // 1 / 3
         (int256 signedCoefficientB, int256 exponentB) = LibDecimalFloatImplementation.div(1, 0, 3, 0);
         assertEq(signedCoefficientB, THREES);
-        assertEq(exponentB, -39);
+        assertEq(exponentB, -76);
 
         // (1 / 9) / (1 / 3)
         (int256 signedCoefficient, int256 exponent) =
             LibDecimalFloatImplementation.div(signedCoefficientA, exponentA, signedCoefficientB, exponentB);
-        assertEq(signedCoefficient, 333333333333333333333333333333333333336);
-        assertEq(exponent, -39);
+        assertEq(signedCoefficient, THREES);
+        assertEq(exponent, -76);
+
+        // (1 / 3) / (1 / 9) == 3
+        (signedCoefficient, exponent) =
+            LibDecimalFloatImplementation.div(signedCoefficientB, exponentB, signedCoefficientA, exponentA);
+        assertEq(signedCoefficient, 3e76);
+        assertEq(exponent, -76);
     }
 
     /// forge-config: default.fuzz.runs = 100
@@ -102,7 +108,7 @@ contract LibDecimalFloatImplementationDivTest is Test {
         int256 di = 0;
         while (true) {
             int256 i = 1;
-            int256 j = -39 - di;
+            int256 j = -76 - di;
             while (true) {
                 // want to see full precision on the THREES regardless of the
                 // scale of the numerator and denominator.

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

@@ -6,10 +6,16 @@ import {LibDecimalFloatSlow} from "test/lib/LibDecimalFloatSlow.sol";
 import {
     LibDecimalFloatImplementation,
     EXPONENT_MIN,
-    EXPONENT_MAX
+    EXPONENT_MAX,
+    MulDivOverflow
 } from "src/lib/implementation/LibDecimalFloatImplementation.sol";
 
 contract LibDecimalFloatImplementationInvTest is Test {
+    function invExternal(int256 signedCoefficient, int256 exponent) external pure returns (int256, int256) {
+        (signedCoefficient, exponent) = LibDecimalFloatImplementation.inv(signedCoefficient, exponent);
+        return (signedCoefficient, exponent);
+    }
+
     /// Compare reference.
     function testInvReference(int256 signedCoefficient, int256 exponent) external pure {
         vm.assume(signedCoefficient != 0);
@@ -33,4 +39,9 @@ contract LibDecimalFloatImplementationInvTest is Test {
         (int256 outputSignedCoefficient, int256 outputExponent) = LibDecimalFloatSlow.invSlow(3e37, -37);
         (outputSignedCoefficient, outputExponent);
     }
+
+    function testInv0() external {
+        vm.expectRevert(abi.encodeWithSelector(MulDivOverflow.selector, 1e76, 1e75, 0));
+        this.invExternal(0, 0);
+    }
 }

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

@@ -32,7 +32,7 @@ contract LibDecimalFloatImplementationPow10Test is LogTest {
         checkPow10(1, 4, 1000, 9997);
     }
 
-    function testExactLookups() external {
+    function testExactLookupsPow10() external {
         // 10^2 = 100
         checkPow10(2, 0, 1000, -1);
         // 10^3 = 1000
@@ -56,7 +56,8 @@ contract LibDecimalFloatImplementationPow10Test is LogTest {
         checkPow10(0.5e37, -37, 3162, -3);
 
         checkPow10(0.3e37, -37, 1995, -3);
-        checkPow10(-0.3e37, -37, 0.501253132832080200501253132832080200501e39, -39);
+        // 10^-0.3 = 0.50118723362
+        checkPow10(-0.3e37, -37, 0.5012531328320802005012531328320802005012531328320802005012531328320802005012e76, -76);
     }
 
     function testInterpolatedLookupsPower() external {