PRBMath.sol

// SPDX-License-Identifier: Unlicense
pragma solidity >=0.8.13;

/// @dev Common mathematical functions used in both SD59x18 and UD60x18. Note that this shared library does not
/// always operate with SD59x18 and UD60x18 numbers.
library PRBMath {
    /// CUSTOM ERRORS ///

    /// @notice Emitted when the ending result in the fixed-point version of `mulDiv` would overflow uint256.
    error PRBMath__MulDivFixedPointOverflow(uint256 x, uint256 y);

    /// @notice Emitted when the ending result in `mulDiv` would overflow uint256.
    error PRBMath__MulDivOverflow(uint256 x, uint256 y, uint256 denominator);

    /// @notice Emitted when attempting to run `mulDiv` with one of the inputs `type(int256).min`.
    error PRBMath__MulDivSignedInputTooSmall();

    /// @notice Emitted when the ending result in the signed version of `mulDiv` would overflow int256.
    error PRBMath__MulDivSignedOverflow(int256 x, int256 y);

    /// CONSTANTS ///

    /// @dev How many trailing decimals can be represented.
    uint256 internal constant SCALE = 1e18;

    /// @dev Largest power of two that is a divisor of `SCALE`.
    uint256 internal constant SCALE_LPOTD = 262144;

    /// @dev The `SCALE` number inverted mod 2^256.
    uint256 internal constant SCALE_INVERSE =
        78156646155174841979727994598816262306175212592076161876661_508869554232690281;

    /// FUNCTIONS ///

    /// @notice Calculates the binary exponent of x using the binary fraction method.
    /// @dev Has to use 192.64-bit fixed-point numbers.
    /// See https://ethereum.stackexchange.com/a/96594/24693.
    /// @param x The exponent as an unsigned 192.64-bit fixed-point number.
    /// @return result The result as an unsigned 60.18-decimal fixed-point number.
    function exp2(uint256 x) internal pure returns (uint256 result) {
        unchecked {
            // Start from 0.5 in the 192.64-bit fixed-point format.
            result = 0x800000000000000000000000000000000000000000000000;

            // Multiply the result by root(2, 2^-i) when the bit at position i is 1. None of the intermediary results overflows
            // because the initial result is 2^191 and all magic factors are less than 2^65.
            if (x & 0xFF00000000000000 > 0) {
                if (x & 0x8000000000000000 > 0) {
                    result = (result * 0x16A09E667F3BCC909) >> 64;
                }
                if (x & 0x4000000000000000 > 0) {
                    result = (result * 0x1306FE0A31B7152DF) >> 64;
                }
                if (x & 0x2000000000000000 > 0) {
                    result = (result * 0x1172B83C7D517ADCE) >> 64;
                }
                if (x & 0x1000000000000000 > 0) {
                    result = (result * 0x10B5586CF9890F62A) >> 64;
                }
                if (x & 0x800000000000000 > 0) {
                    result = (result * 0x1059B0D31585743AE) >> 64;
                }
                if (x & 0x400000000000000 > 0) {
                    result = (result * 0x102C9A3E778060EE7) >> 64;
                }
                if (x & 0x200000000000000 > 0) {
                    result = (result * 0x10163DA9FB33356D8) >> 64;
                }
                if (x & 0x100000000000000 > 0) {
                    result = (result * 0x100B1AFA5ABCBED61) >> 64;
                }
            }

            if (x & 0xFF000000000000 > 0) {
                if (x & 0x80000000000000 > 0) {
                    result = (result * 0x10058C86DA1C09EA2) >> 64;
                }
                if (x & 0x40000000000000 > 0) {
                    result = (result * 0x1002C605E2E8CEC50) >> 64;
                }
                if (x & 0x20000000000000 > 0) {
                    result = (result * 0x100162F3904051FA1) >> 64;
                }
                if (x & 0x10000000000000 > 0) {
                    result = (result * 0x1000B175EFFDC76BA) >> 64;
                }
                if (x & 0x8000000000000 > 0) {
                    result = (result * 0x100058BA01FB9F96D) >> 64;
                }
                if (x & 0x4000000000000 > 0) {
                    result = (result * 0x10002C5CC37DA9492) >> 64;
                }
                if (x & 0x2000000000000 > 0) {
                    result = (result * 0x1000162E525EE0547) >> 64;
                }
                if (x & 0x1000000000000 > 0) {
                    result = (result * 0x10000B17255775C04) >> 64;
                }
            }

            if (x & 0xFF0000000000 > 0) {
                if (x & 0x800000000000 > 0) {
                    result = (result * 0x1000058B91B5BC9AE) >> 64;
                }
                if (x & 0x400000000000 > 0) {
                    result = (result * 0x100002C5C89D5EC6D) >> 64;
                }
                if (x & 0x200000000000 > 0) {
                    result = (result * 0x10000162E43F4F831) >> 64;
                }
                if (x & 0x100000000000 > 0) {
                    result = (result * 0x100000B1721BCFC9A) >> 64;
                }
                if (x & 0x80000000000 > 0) {
                    result = (result * 0x10000058B90CF1E6E) >> 64;
                }
                if (x & 0x40000000000 > 0) {
                    result = (result * 0x1000002C5C863B73F) >> 64;
                }
                if (x & 0x20000000000 > 0) {
                    result = (result * 0x100000162E430E5A2) >> 64;
                }
                if (x & 0x10000000000 > 0) {
                    result = (result * 0x1000000B172183551) >> 64;
                }
            }

            if (x & 0xFF00000000 > 0) {
                if (x & 0x8000000000 > 0) {
                    result = (result * 0x100000058B90C0B49) >> 64;
                }
                if (x & 0x4000000000 > 0) {
                    result = (result * 0x10000002C5C8601CC) >> 64;
                }
                if (x & 0x2000000000 > 0) {
                    result = (result * 0x1000000162E42FFF0) >> 64;
                }
                if (x & 0x1000000000 > 0) {
                    result = (result * 0x10000000B17217FBB) >> 64;
                }
                if (x & 0x800000000 > 0) {
                    result = (result * 0x1000000058B90BFCE) >> 64;
                }
                if (x & 0x400000000 > 0) {
                    result = (result * 0x100000002C5C85FE3) >> 64;
                }
                if (x & 0x200000000 > 0) {
                    result = (result * 0x10000000162E42FF1) >> 64;
                }
                if (x & 0x100000000 > 0) {
                    result = (result * 0x100000000B17217F8) >> 64;
                }
            }

            if (x & 0xFF00000000 > 0) {
                if (x & 0x80000000 > 0) {
                    result = (result * 0x10000000058B90BFC) >> 64;
                }
                if (x & 0x40000000 > 0) {
                    result = (result * 0x1000000002C5C85FE) >> 64;
                }
                if (x & 0x20000000 > 0) {
                    result = (result * 0x100000000162E42FF) >> 64;
                }
                if (x & 0x10000000 > 0) {
                    result = (result * 0x1000000000B17217F) >> 64;
                }
                if (x & 0x8000000 > 0) {
                    result = (result * 0x100000000058B90C0) >> 64;
                }
                if (x & 0x4000000 > 0) {
                    result = (result * 0x10000000002C5C860) >> 64;
                }
                if (x & 0x2000000 > 0) {
                    result = (result * 0x1000000000162E430) >> 64;
                }
                if (x & 0x1000000 > 0) {
                    result = (result * 0x10000000000B17218) >> 64;
                }
            }

            if (x & 0xFF0000 > 0) {
                if (x & 0x800000 > 0) {
                    result = (result * 0x1000000000058B90C) >> 64;
                }
                if (x & 0x400000 > 0) {
                    result = (result * 0x100000000002C5C86) >> 64;
                }
                if (x & 0x200000 > 0) {
                    result = (result * 0x10000000000162E43) >> 64;
                }
                if (x & 0x100000 > 0) {
                    result = (result * 0x100000000000B1721) >> 64;
                }
                if (x & 0x80000 > 0) {
                    result = (result * 0x10000000000058B91) >> 64;
                }
                if (x & 0x40000 > 0) {
                    result = (result * 0x1000000000002C5C8) >> 64;
                }
                if (x & 0x20000 > 0) {
                    result = (result * 0x100000000000162E4) >> 64;
                }
                if (x & 0x10000 > 0) {
                    result = (result * 0x1000000000000B172) >> 64;
                }
            }

            if (x & 0xFF00 > 0) {
                if (x & 0x8000 > 0) {
                    result = (result * 0x100000000000058B9) >> 64;
                }
                if (x & 0x4000 > 0) {
                    result = (result * 0x10000000000002C5D) >> 64;
                }
                if (x & 0x2000 > 0) {
                    result = (result * 0x1000000000000162E) >> 64;
                }
                if (x & 0x1000 > 0) {
                    result = (result * 0x10000000000000B17) >> 64;
                }
                if (x & 0x800 > 0) {
                    result = (result * 0x1000000000000058C) >> 64;
                }
                if (x & 0x400 > 0) {
                    result = (result * 0x100000000000002C6) >> 64;
                }
                if (x & 0x200 > 0) {
                    result = (result * 0x10000000000000163) >> 64;
                }
                if (x & 0x100 > 0) {
                    result = (result * 0x100000000000000B1) >> 64;
                }
            }

            if (x & 0xFF > 0) {
                if (x & 0x80 > 0) {
                    result = (result * 0x10000000000000059) >> 64;
                }
                if (x & 0x40 > 0) {
                    result = (result * 0x1000000000000002C) >> 64;
                }
                if (x & 0x20 > 0) {
                    result = (result * 0x10000000000000016) >> 64;
                }
                if (x & 0x10 > 0) {
                    result = (result * 0x1000000000000000B) >> 64;
                }
                if (x & 0x8 > 0) {
                    result = (result * 0x10000000000000006) >> 64;
                }
                if (x & 0x4 > 0) {
                    result = (result * 0x10000000000000003) >> 64;
                }
                if (x & 0x2 > 0) {
                    result = (result * 0x10000000000000001) >> 64;
                }
                if (x & 0x1 > 0) {
                    result = (result * 0x10000000000000001) >> 64;
                }
            }

            // We're doing two things at the same time:
            //
            //   1. Multiply the result by 2^n + 1, where "2^n" is the integer part and the one is added to account for
            //      the fact that we initially set the result to 0.5. This is accomplished by subtracting from 191
            //      rather than 192.
            //   2. Convert the result to the unsigned 60.18-decimal fixed-point format.
            //
            // This works because 2^(191-ip) = 2^ip / 2^191, where "ip" is the integer part "2^n".
            result *= SCALE;
            result >>= (191 - (x >> 64));
        }
    }

    /// @notice Finds the zero-based index of the first 1 in the binary representation of x.
    /// @dev See the note on msb in the "Find First Set" Wikipedia article https://en.wikipedia.org/wiki/Find_first_set
    /// @param x The uint256 number for which to find the index of the most significant bit.
    /// @return result The index of the most significant bit as an uint256.
    function msb(uint256 x) internal pure returns (uint256 result) {
        unchecked {
            if (x >= 2**128) {
                x >>= 128;
                result += 128;
            }
            if (x >= 2**64) {
                x >>= 64;
                result += 64;
            }
            if (x >= 2**32) {
                x >>= 32;
                result += 32;
            }
            if (x >= 2**16) {
                x >>= 16;
                result += 16;
            }
            if (x >= 2**8) {
                x >>= 8;
                result += 8;
            }
            if (x >= 2**4) {
                x >>= 4;
                result += 4;
            }
            if (x >= 2**2) {
                x >>= 2;
                result += 2;
            }
            // No need to shift x any more.
            if (x >= 2**1) {
                result += 1;
            }
        }
    }

    /// @notice Calculates floor(x*y÷denominator) with full precision.
    ///
    /// @dev Credit to Remco Bloemen under MIT license https://xn--2-umb.com/21/muldiv.
    ///
    /// Requirements:
    /// - The denominator cannot be zero.
    /// - The result must fit within uint256.
    ///
    /// Caveats:
    /// - This function does not work with fixed-point numbers.
    ///
    /// @param x The multiplicand as an uint256.
    /// @param y The multiplier as an uint256.
    /// @param denominator The divisor as an uint256.
    /// @return result The result as an uint256.
    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 {
            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. "prod1 > 0 == false" is equivalent to "prod1 == 0" but faster.
        if (prod1 > 0 == false) {
            unchecked {
                return prod0 / denominator;
            }
        }

        // Make sure the result is less than 2^256. Also prevents denominator == 0.
        if (prod1 >= denominator) {
            revert PRBMath__MulDivOverflow(x, y, denominator);
        }

        ///////////////////////////////////////////////
        // 512 by 256 division.
        ///////////////////////////////////////////////

        // Make division exact by subtracting the remainder from [prod1 prod0].
        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)
        }

        // Factor powers of two out of denominator and compute largest power of two divisor of denominator. Always >= 1.
        // See https://cs.stackexchange.com/q/138556/92363.
        unchecked {
            // Does not overflow because the denominator cannot be zero at this stage in the function.
            uint256 lpotdod = denominator & (~denominator + 1);
            assembly {
                // Divide denominator by lpotdod.
                denominator := div(denominator, lpotdod)

                // Divide [prod1 prod0] by lpotdod.
                prod0 := div(prod0, lpotdod)

                // Flip lpotdod such that it is 2^256 / lpotdod. If lpotdod is zero, then it becomes one.
                lpotdod := add(div(sub(0, lpotdod), lpotdod), 1)
            }

            // Shift in bits from prod1 into prod0.
            prod0 |= prod1 * lpotdod;

            // 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.
            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;
        }
    }

    /// @notice Calculates floor(x*y÷1e18) with full precision.
    ///
    /// @dev Variant of `mulDiv` with constant folding, i.e. in which the denominator is always 1e18. Before returning the
    /// final result, we add 1 if `(x * y) % SCALE >= HALF_SCALE`. Without this adjustment, 6.6e-19 would be truncated to 0
    /// instead of being rounded to 1e-18. See "Listing 6" and text above it at https://accu.org/index.php/journals/1717.
    ///
    /// Requirements:
    /// - The result must fit within uint256.
    ///
    /// Caveats:
    /// - The body is purposely left uncommented; to understand how this works, see the NatSpec comments in `mulDiv`.
    /// - It is assumed that the result can never be `type(uint256).max` when x and y solve the following two equations:
    ///     1. x * y = type(uint256).max * SCALE
    ///     2. (x * y) % SCALE >= SCALE / 2
    ///
    /// @param x The multiplicand as an unsigned 60.18-decimal fixed-point number.
    /// @param y The multiplier as an unsigned 60.18-decimal fixed-point number.
    /// @return result The result as an unsigned 60.18-decimal fixed-point number.
    function mulDiv18(uint256 x, uint256 y) internal pure returns (uint256 result) {
        uint256 prod0;
        uint256 prod1;
        assembly {
            let mm := mulmod(x, y, not(0))
            prod0 := mul(x, y)
            prod1 := sub(sub(mm, prod0), lt(mm, prod0))
        }

        if (prod1 >= SCALE) {
            revert PRBMath__MulDivFixedPointOverflow(x, y);
        }

        uint256 remainder;
        uint256 roundUpUnit;
        assembly {
            remainder := mulmod(x, y, SCALE)
            roundUpUnit := gt(remainder, 499999999999999999)
        }

        if (prod1 > 0 == false) {
            unchecked {
                return (prod0 / SCALE) + roundUpUnit;
            }
        }

        assembly {
            result := add(
                mul(
                    or(
                        div(sub(prod0, remainder), SCALE_LPOTD),
                        mul(sub(prod1, gt(remainder, prod0)), add(div(sub(0, SCALE_LPOTD), SCALE_LPOTD), 1))
                    ),
                    SCALE_INVERSE
                ),
                roundUpUnit
            )
        }
    }

    /// @notice Calculates floor(x*y÷denominator) with full precision.
    ///
    /// @dev An extension of `mulDiv` for signed numbers. Works by computing the signs and the absolute values separately.
    ///
    /// Requirements:
    /// - None of the inputs can be `type(int256).min`.
    /// - The result must fit within int256.
    ///
    /// @param x The multiplicand as an int256.
    /// @param y The multiplier as an int256.
    /// @param denominator The divisor as an int256.
    /// @return result The result as an int256.
    function mulDivSigned(
        int256 x,
        int256 y,
        int256 denominator
    ) internal pure returns (int256 result) {
        if (x == type(int256).min || y == type(int256).min || denominator == type(int256).min) {
            revert PRBMath__MulDivSignedInputTooSmall();
        }

        // Get hold of the absolute values of x, y and the denominator.
        uint256 ax;
        uint256 ay;
        uint256 ad;
        unchecked {
            ax = x < 0 ? uint256(-x) : uint256(x);
            ay = y < 0 ? uint256(-y) : uint256(y);
            ad = denominator < 0 ? uint256(-denominator) : uint256(denominator);
        }

        // Compute the absolute value of (x*y)÷denominator. The result must fit within int256.
        uint256 rAbs = mulDiv(ax, ay, ad);
        if (rAbs > uint256(type(int256).max)) {
            revert PRBMath__MulDivSignedOverflow(x, y);
        }

        // Get the signs of x, y and the denominator.
        uint256 sx;
        uint256 sy;
        uint256 sd;
        assembly {
            // This works thanks to two's complement.
            // "sgt" stands for "signed greater than" and "sub(0,1)" is max uint256.
            sx := sgt(x, sub(0, 1))
            sy := sgt(y, sub(0, 1))
            sd := sgt(denominator, sub(0, 1))
        }

        // XOR over sx, sy and sd. What this does is to check whether there are 1 or 3 negative signs in the inputs.
        // If yes, the result should be negative. Otherwise, it should be positive.
        unchecked {
            result = sx ^ sy ^ sd == 0 ? -int256(rAbs) : int256(rAbs);
        }
    }

    /// @notice Calculates the square root of x, rounding down.
    /// @dev Uses the Babylonian method https://en.wikipedia.org/wiki/Methods_of_computing_square_roots#Babylonian_method.
    ///
    /// Caveats:
    /// - This function does not work with fixed-point numbers.
    ///
    /// @param x The uint256 number for which to calculate the square root.
    /// @return result The result as an uint256.
    function sqrt(uint256 x) internal pure returns (uint256 result) {
        // Equivalent to "x == 0" but faster.
        if (x > 0 == false) {
            return 0;
        }

        // Set the initial guess to the least power of two that is greater than or equal to sqrt(x).
        uint256 xAux = uint256(x);
        result = 1;
        if (xAux >= 2**128) {
            xAux >>= 128;
            result <<= 64;
        }
        if (xAux >= 2**64) {
            xAux >>= 64;
            result <<= 32;
        }
        if (xAux >= 2**32) {
            xAux >>= 32;
            result <<= 16;
        }
        if (xAux >= 2**16) {
            xAux >>= 16;
            result <<= 8;
        }
        if (xAux >= 2**8) {
            xAux >>= 8;
            result <<= 4;
        }
        if (xAux >= 2**4) {
            xAux >>= 4;
            result <<= 2;
        }
        if (xAux >= 2**2) {
            result <<= 1;
        }

        // The operations can never overflow because the result is max 2^127 when it enters this block.
        unchecked {
            // Seven iterations should be enough
            result = (result + x / result) >> 1;
            result = (result + x / result) >> 1;
            result = (result + x / result) >> 1;
            result = (result + x / result) >> 1;
            result = (result + x / result) >> 1;
            result = (result + x / result) >> 1;
            result = (result + x / result) >> 1;

            // Round down the result.
            uint256 roundedDownResult = x / result;
            result >= roundedDownResult ? roundedDownResult : result;
        }
    }
}