Copy
1// SPDX-License-Identifier: MIT
2pragma solidity ^0.8.17;
3
4/// @title Contains 512-bit math functions
5/// @notice Facilitates multiplication and division that can have overflow of an intermediate value without any loss of precision
6/// @dev Handles "phantom overflow" i.e., allows multiplication and division where an intermediate value overflows 256 bits
7library FullMath {
8    /// @notice Calculates floor(a×b÷denominator) with full precision. Throws if result overflows a uint256 or denominator == 0
9    /// @param a The multiplicand
10    /// @param b The multiplier
11    /// @param denominator The divisor
12    /// @return result The 256-bit result
13    /// @dev Credit to Remco Bloemen under MIT license https://xn--2-umb.com/21/muldiv
14    function mulDiv(
15        uint256 a,
16        uint256 b,
17        uint256 denominator
18    ) internal pure returns (uint256 result) {
19        // 512-bit multiply [prod1 prod0] = a * b
20        // Compute the product mod 2**256 and mod 2**256 - 1
21        // then use the Chinese Remainder Theorem to reconstruct
22        // the 512 bit result. The result is stored in two 256
23        // variables such that product = prod1 * 2**256 + prod0
24        uint256 prod0; // Least significant 256 bits of the product
25        uint256 prod1; // Most significant 256 bits of the product
26        assembly {
27            let mm := mulmod(a, b, not(0))
28            prod0 := mul(a, b)
29            prod1 := sub(sub(mm, prod0), lt(mm, prod0))
30        }
31
32        // Handle non-overflow cases, 256 by 256 division
33        if (prod1 == 0) {
34            require(denominator > 0);
35            assembly {
36                result := div(prod0, denominator)
37            }
38            return result;
39        }
40
41        // Make sure the result is less than 2**256.
42        // Also prevents denominator == 0
43        require(denominator > prod1);
44
45        ///////////////////////////////////////////////
46        // 512 by 256 division.
47        ///////////////////////////////////////////////
48
49        // Make division exact by subtracting the remainder from [prod1 prod0]
50        // Compute remainder using mulmod
51        uint256 remainder;
52        assembly {
53            remainder := mulmod(a, b, denominator)
54        }
55        // Subtract 256 bit number from 512 bit number
56        assembly {
57            prod1 := sub(prod1, gt(remainder, prod0))
58            prod0 := sub(prod0, remainder)
59        }
60
61        // Factor powers of two out of denominator
62        // Compute largest power of two divisor of denominator.
63        // Always >= 1.
64        // uint256 twos = -denominator & denominator;
65        uint256 twos = denominator & (~denominator + 1);
66        // Divide denominator by power of two
67        assembly {
68            denominator := div(denominator, twos)
69        }
70
71        // Divide [prod1 prod0] by the factors of two
72        assembly {
73            prod0 := div(prod0, twos)
74        }
75        // Shift in bits from prod1 into prod0. For this we need
76        // to flip `twos` such that it is 2**256 / twos.
77        // If twos is zero, then it becomes one
78        assembly {
79            twos := add(div(sub(0, twos), twos), 1)
80        }
81        prod0 |= prod1 * twos;
82
83        // Invert denominator mod 2**256
84        // Now that denominator is an odd number, it has an inverse
85        // modulo 2**256 such that denominator * inv = 1 mod 2**256.
86        // Compute the inverse by starting with a seed that is correct
87        // correct for four bits. That is, denominator * inv = 1 mod 2**4
88        uint256 inv = (3 * denominator) ^ 2;
89        // Now use Newton-Raphson iteration to improve the precision.
90        // Thanks to Hensel's lifting lemma, this also works in modular
91        // arithmetic, doubling the correct bits in each step.
92        inv *= 2 - denominator * inv; // inverse mod 2**8
93        inv *= 2 - denominator * inv; // inverse mod 2**16
94        inv *= 2 - denominator * inv; // inverse mod 2**32
95        inv *= 2 - denominator * inv; // inverse mod 2**64
96        inv *= 2 - denominator * inv; // inverse mod 2**128
97        inv *= 2 - denominator * inv; // inverse mod 2**256
98
99        // Because the division is now exact we can divide by multiplying
100        // with the modular inverse of denominator. This will give us the
101        // correct result modulo 2**256. Since the precoditions guarantee
102        // that the outcome is less than 2**256, this is the final result.
103        // We don't need to compute the high bits of the result and prod1
104        // is no longer required.
105        result = prod0 * inv;
106        return result;
107    }