// SPDX-License-Identifier: WTFPL pragma solidity >=0.8.0; import "./PRBMathCommon.sol"; /// @title PRBMathSD59x18 /// @author Paul Razvan Berg /// @notice Smart contract library for advanced fixed-point math. It works with int256 numbers considered to have 18 /// trailing decimals. We call this number representation signed 59.18-decimal fixed-point, since the numbers can have /// a sign and there can be up to 59 digits in the integer part and up to 18 decimals in the fractional part. The numbers /// are bound by the minimum and the maximum values permitted by the Solidity type int256. library PRBMathSD59x18 { /// @dev log2(e) as a signed 59.18-decimal fixed-point number. int256 internal constant LOG2_E = 1442695040888963407; /// @dev Half the SCALE number. int256 internal constant HALF_SCALE = 5e17; /// @dev The maximum value a signed 59.18-decimal fixed-point number can have. int256 internal constant MAX_SD59x18 = 57896044618658097711785492504343953926634992332820282019728792003956564819967; /// @dev The maximum whole value a signed 59.18-decimal fixed-point number can have. int256 internal constant MAX_WHOLE_SD59x18 = 57896044618658097711785492504343953926634992332820282019728000000000000000000; /// @dev The minimum value a signed 59.18-decimal fixed-point number can have. int256 internal constant MIN_SD59x18 = -57896044618658097711785492504343953926634992332820282019728792003956564819968; /// @dev The minimum whole value a signed 59.18-decimal fixed-point number can have. int256 internal constant MIN_WHOLE_SD59x18 = -57896044618658097711785492504343953926634992332820282019728000000000000000000; /// @dev How many trailing decimals can be represented. int256 internal constant SCALE = 1e18; /// INTERNAL FUNCTIONS /// /// @notice Calculate the absolute value of x. /// /// @dev Requirements: /// - x must be greater than MIN_SD59x18. /// /// @param x The number to calculate the absolute value for. /// @param result The absolute value of x. function abs(int256 x) internal pure returns (int256 result) { unchecked { require(x > MIN_SD59x18); result = x < 0 ? -x : x; } } /// @notice Calculates arithmetic average of x and y, rounding down. /// @param x The first operand as a signed 59.18-decimal fixed-point number. /// @param y The second operand as a signed 59.18-decimal fixed-point number. /// @return result The arithmetic average as a signed 59.18-decimal fixed-point number. function avg(int256 x, int256 y) internal pure returns (int256 result) { // The operations can never overflow. unchecked { // The last operand checks if both x and y are odd and if that is the case, we add 1 to the result. We need // to do this because if both numbers are odd, the 0.5 remainder gets truncated twice. result = (x >> 1) + (y >> 1) + (x & y & 1); } } /// @notice Yields the least greatest signed 59.18 decimal fixed-point number greater than or equal to x. /// /// @dev Optimised for fractional value inputs, because for every whole value there are (1e18 - 1) fractional counterparts. /// See https://en.wikipedia.org/wiki/Floor_and_ceiling_functions. /// /// Requirements: /// - x must be less than or equal to MAX_WHOLE_SD59x18. /// /// @param x The signed 59.18-decimal fixed-point number to ceil. /// @param result The least integer greater than or equal to x, as a signed 58.18-decimal fixed-point number. function ceil(int256 x) internal pure returns (int256 result) { require(x <= MAX_WHOLE_SD59x18); unchecked { int256 remainder = x % SCALE; if (remainder == 0) { result = x; } else { // Solidity uses C fmod style, which returns a modulus with the same sign as x. result = x - remainder; if (x > 0) { result += SCALE; } } } } /// @notice Divides two signed 59.18-decimal fixed-point numbers, returning a new signed 59.18-decimal fixed-point number. /// /// @dev Variant of "mulDiv" that works with signed numbers. Works by computing the signs and the absolute values separately. /// /// Requirements: /// - All from "PRBMathCommon.mulDiv". /// - None of the inputs can be type(int256).min. /// - y cannot be zero. /// - The result must fit within int256. /// /// Caveats: /// - All from "PRBMathCommon.mulDiv". /// /// @param x The numerator as a signed 59.18-decimal fixed-point number. /// @param y The denominator as a signed 59.18-decimal fixed-point number. /// @param result The quotient as a signed 59.18-decimal fixed-point number. function div(int256 x, int256 y) internal pure returns (int256 result) { require(x > type(int256).min); require(y > type(int256).min); // Get hold of the absolute values of x and y. uint256 ax; uint256 ay; unchecked { ax = x < 0 ? uint256(-x) : uint256(x); ay = y < 0 ? uint256(-y) : uint256(y); } // Compute the absolute value of (x*SCALE)÷y. The result must fit within int256. uint256 resultUnsigned = PRBMathCommon.mulDiv(ax, uint256(SCALE), ay); require(resultUnsigned <= uint256(type(int256).max)); // Get the signs of x and y. uint256 sx; uint256 sy; assembly { sx := sgt(x, sub(0, 1)) sy := sgt(y, sub(0, 1)) } // XOR over sx and sy. This is basically checking whether the inputs have the same sign. If yes, the result // should be positive. Otherwise, it should be negative. result = sx ^ sy == 1 ? -int256(resultUnsigned) : int256(resultUnsigned); } /// @notice Returns Euler's number as a signed 59.18-decimal fixed-point number. /// @dev See https://en.wikipedia.org/wiki/E_(mathematical_constant). function e() internal pure returns (int256 result) { result = 2718281828459045235; } /// @notice Calculates the natural exponent of x. /// /// @dev Based on the insight that e^x = 2^(x * log2(e)). /// /// Requirements: /// - All from "log2". /// - x must be less than 88722839111672999628. /// /// @param x The exponent as a signed 59.18-decimal fixed-point number. /// @return result The result as a signed 59.18-decimal fixed-point number. function exp(int256 x) internal pure returns (int256 result) { // Without this check, the value passed to "exp2" would be less than -59794705707972522261. if (x < -41446531673892822322) { return 0; } // Without this check, the value passed to "exp2" would be greater than 128e18. require(x < 88722839111672999628); // Do the fixed-point multiplication inline to save gas. unchecked { int256 doubleScaleProduct = x * LOG2_E; result = exp2((doubleScaleProduct + HALF_SCALE) / SCALE); } } /// @notice Calculates the binary exponent of x using the binary fraction method. /// /// @dev See https://ethereum.stackexchange.com/q/79903/24693. /// /// Requirements: /// - x must be 128e18 or less. /// - The result must fit within MAX_SD59x18. /// /// Caveats: /// - For any x less than -59794705707972522261, the result is zero. /// /// @param x The exponent as a signed 59.18-decimal fixed-point number. /// @return result The result as a signed 59.18-decimal fixed-point number. function exp2(int256 x) internal pure returns (int256 result) { // This works because 2^(-x) = 1/2^x. if (x < 0) { // 2**59.794705707972522262 is the maximum number whose inverse does not turn into zero. if (x < -59794705707972522261) { return 0; } // Do the fixed-point inversion inline to save gas. The numerator is SCALE * SCALE. unchecked { result = 1e36 / exp2(-x); } } else { // 2**128 doesn't fit within the 128.128-bit fixed-point representation. require(x < 128e18); unchecked { // Convert x to the 128.128-bit fixed-point format. uint256 x128x128 = (uint256(x) << 128) / uint256(SCALE); // Safe to convert the result to int256 directly because the maximum input allowed is 128e18. result = int256(PRBMathCommon.exp2(x128x128)); } } } /// @notice Yields the greatest signed 59.18 decimal fixed-point number less than or equal to x. /// /// @dev Optimised for fractional value inputs, because for every whole value there are (1e18 - 1) fractional counterparts. /// See https://en.wikipedia.org/wiki/Floor_and_ceiling_functions. /// /// Requirements: /// - x must be greater than or equal to MIN_WHOLE_SD59x18. /// /// @param x The signed 59.18-decimal fixed-point number to floor. /// @param result The greatest integer less than or equal to x, as a signed 58.18-decimal fixed-point number. function floor(int256 x) internal pure returns (int256 result) { require(x >= MIN_WHOLE_SD59x18); unchecked { int256 remainder = x % SCALE; if (remainder == 0) { result = x; } else { // Solidity uses C fmod style, which returns a modulus with the same sign as x. result = x - remainder; if (x < 0) { result -= SCALE; } } } } /// @notice Yields the excess beyond the floor of x for positive numbers and the part of the number to the right /// of the radix point for negative numbers. /// @dev Based on the odd function definition. https://en.wikipedia.org/wiki/Fractional_part /// @param x The signed 59.18-decimal fixed-point number to get the fractional part of. /// @param result The fractional part of x as a signed 59.18-decimal fixed-point number. function frac(int256 x) internal pure returns (int256 result) { unchecked { result = x % SCALE; } } /// @notice Calculates geometric mean of x and y, i.e. sqrt(x * y), rounding down. /// /// @dev Requirements: /// - x * y must fit within MAX_SD59x18, lest it overflows. /// - x * y cannot be negative. /// /// @param x The first operand as a signed 59.18-decimal fixed-point number. /// @param y The second operand as a signed 59.18-decimal fixed-point number. /// @return result The result as a signed 59.18-decimal fixed-point number. function gm(int256 x, int256 y) internal pure returns (int256 result) { if (x == 0) { return 0; } unchecked { // Checking for overflow this way is faster than letting Solidity do it. int256 xy = x * y; require(xy / x == y); // The product cannot be negative. require(xy >= 0); // We don't need to multiply by the SCALE here because the x*y product had already picked up a factor of SCALE // during multiplication. See the comments within the "sqrt" function. result = int256(PRBMathCommon.sqrt(uint256(xy))); } } /// @notice Calculates 1 / x, rounding towards zero. /// /// @dev Requirements: /// - x cannot be zero. /// /// @param x The signed 59.18-decimal fixed-point number for which to calculate the inverse. /// @return result The inverse as a signed 59.18-decimal fixed-point number. function inv(int256 x) internal pure returns (int256 result) { unchecked { // 1e36 is SCALE * SCALE. result = 1e36 / x; } } /// @notice Calculates the natural logarithm of x. /// /// @dev Based on the insight that ln(x) = log2(x) / log2(e). /// /// Requirements: /// - All from "log2". /// /// Caveats: /// - All from "log2". /// - This doesn't return exactly 1 for 2718281828459045235, for that we would need more fine-grained precision. /// /// @param x The signed 59.18-decimal fixed-point number for which to calculate the natural logarithm. /// @return result The natural logarithm as a signed 59.18-decimal fixed-point number. function ln(int256 x) internal pure returns (int256 result) { // Do the fixed-point multiplication inline to save gas. This is overflow-safe because the maximum value that log2(x) // can return is 195205294292027477728. unchecked { result = (log2(x) * SCALE) / LOG2_E; } } /// @notice Calculates the common logarithm of x. /// /// @dev First checks if x is an exact power of ten and it stops if yes. If it's not, calculates the common /// logarithm based on the insight that log10(x) = log2(x) / log2(10). /// /// Requirements: /// - All from "log2". /// /// Caveats: /// - All from "log2". /// /// @param x The signed 59.18-decimal fixed-point number for which to calculate the common logarithm. /// @return result The common logarithm as a signed 59.18-decimal fixed-point number. function log10(int256 x) internal pure returns (int256 result) { require(x > 0); // Note that the "mul" in this block is the assembly mul operation, not the "mul" function defined in this contract. // prettier-ignore assembly { switch x case 1 { result := mul(SCALE, sub(0, 18)) } case 10 { result := mul(SCALE, sub(1, 18)) } case 100 { result := mul(SCALE, sub(2, 18)) } case 1000 { result := mul(SCALE, sub(3, 18)) } case 10000 { result := mul(SCALE, sub(4, 18)) } case 100000 { result := mul(SCALE, sub(5, 18)) } case 1000000 { result := mul(SCALE, sub(6, 18)) } case 10000000 { result := mul(SCALE, sub(7, 18)) } case 100000000 { result := mul(SCALE, sub(8, 18)) } case 1000000000 { result := mul(SCALE, sub(9, 18)) } case 10000000000 { result := mul(SCALE, sub(10, 18)) } case 100000000000 { result := mul(SCALE, sub(11, 18)) } case 1000000000000 { result := mul(SCALE, sub(12, 18)) } case 10000000000000 { result := mul(SCALE, sub(13, 18)) } case 100000000000000 { result := mul(SCALE, sub(14, 18)) } case 1000000000000000 { result := mul(SCALE, sub(15, 18)) } case 10000000000000000 { result := mul(SCALE, sub(16, 18)) } case 100000000000000000 { result := mul(SCALE, sub(17, 18)) } case 1000000000000000000 { result := 0 } case 10000000000000000000 { result := SCALE } case 100000000000000000000 { result := mul(SCALE, 2) } case 1000000000000000000000 { result := mul(SCALE, 3) } case 10000000000000000000000 { result := mul(SCALE, 4) } case 100000000000000000000000 { result := mul(SCALE, 5) } case 1000000000000000000000000 { result := mul(SCALE, 6) } case 10000000000000000000000000 { result := mul(SCALE, 7) } case 100000000000000000000000000 { result := mul(SCALE, 8) } case 1000000000000000000000000000 { result := mul(SCALE, 9) } case 10000000000000000000000000000 { result := mul(SCALE, 10) } case 100000000000000000000000000000 { result := mul(SCALE, 11) } case 1000000000000000000000000000000 { result := mul(SCALE, 12) } case 10000000000000000000000000000000 { result := mul(SCALE, 13) } case 100000000000000000000000000000000 { result := mul(SCALE, 14) } case 1000000000000000000000000000000000 { result := mul(SCALE, 15) } case 10000000000000000000000000000000000 { result := mul(SCALE, 16) } case 100000000000000000000000000000000000 { result := mul(SCALE, 17) } case 1000000000000000000000000000000000000 { result := mul(SCALE, 18) } case 10000000000000000000000000000000000000 { result := mul(SCALE, 19) } case 100000000000000000000000000000000000000 { result := mul(SCALE, 20) } case 1000000000000000000000000000000000000000 { result := mul(SCALE, 21) } case 10000000000000000000000000000000000000000 { result := mul(SCALE, 22) } case 100000000000000000000000000000000000000000 { result := mul(SCALE, 23) } case 1000000000000000000000000000000000000000000 { result := mul(SCALE, 24) } case 10000000000000000000000000000000000000000000 { result := mul(SCALE, 25) } case 100000000000000000000000000000000000000000000 { result := mul(SCALE, 26) } case 1000000000000000000000000000000000000000000000 { result := mul(SCALE, 27) } case 10000000000000000000000000000000000000000000000 { result := mul(SCALE, 28) } case 100000000000000000000000000000000000000000000000 { result := mul(SCALE, 29) } case 1000000000000000000000000000000000000000000000000 { result := mul(SCALE, 30) } case 10000000000000000000000000000000000000000000000000 { result := mul(SCALE, 31) } case 100000000000000000000000000000000000000000000000000 { result := mul(SCALE, 32) } case 1000000000000000000000000000000000000000000000000000 { result := mul(SCALE, 33) } case 10000000000000000000000000000000000000000000000000000 { result := mul(SCALE, 34) } case 100000000000000000000000000000000000000000000000000000 { result := mul(SCALE, 35) } case 1000000000000000000000000000000000000000000000000000000 { result := mul(SCALE, 36) } case 10000000000000000000000000000000000000000000000000000000 { result := mul(SCALE, 37) } case 100000000000000000000000000000000000000000000000000000000 { result := mul(SCALE, 38) } case 1000000000000000000000000000000000000000000000000000000000 { result := mul(SCALE, 39) } case 10000000000000000000000000000000000000000000000000000000000 { result := mul(SCALE, 40) } case 100000000000000000000000000000000000000000000000000000000000 { result := mul(SCALE, 41) } case 1000000000000000000000000000000000000000000000000000000000000 { result := mul(SCALE, 42) } case 10000000000000000000000000000000000000000000000000000000000000 { result := mul(SCALE, 43) } case 100000000000000000000000000000000000000000000000000000000000000 { result := mul(SCALE, 44) } case 1000000000000000000000000000000000000000000000000000000000000000 { result := mul(SCALE, 45) } case 10000000000000000000000000000000000000000000000000000000000000000 { result := mul(SCALE, 46) } case 100000000000000000000000000000000000000000000000000000000000000000 { result := mul(SCALE, 47) } case 1000000000000000000000000000000000000000000000000000000000000000000 { result := mul(SCALE, 48) } case 10000000000000000000000000000000000000000000000000000000000000000000 { result := mul(SCALE, 49) } case 100000000000000000000000000000000000000000000000000000000000000000000 { result := mul(SCALE, 50) } case 1000000000000000000000000000000000000000000000000000000000000000000000 { result := mul(SCALE, 51) } case 10000000000000000000000000000000000000000000000000000000000000000000000 { result := mul(SCALE, 52) } case 100000000000000000000000000000000000000000000000000000000000000000000000 { result := mul(SCALE, 53) } case 1000000000000000000000000000000000000000000000000000000000000000000000000 { result := mul(SCALE, 54) } case 10000000000000000000000000000000000000000000000000000000000000000000000000 { result := mul(SCALE, 55) } case 100000000000000000000000000000000000000000000000000000000000000000000000000 { result := mul(SCALE, 56) } case 1000000000000000000000000000000000000000000000000000000000000000000000000000 { result := mul(SCALE, 57) } case 10000000000000000000000000000000000000000000000000000000000000000000000000000 { result := mul(SCALE, 58) } default { result := MAX_SD59x18 } } if (result == MAX_SD59x18) { // Do the fixed-point division inline to save gas. The denominator is log2(10). unchecked { result = (log2(x) * SCALE) / 332192809488736234; } } } /// @notice Calculates the binary logarithm of x. /// /// @dev Based on the iterative approximation algorithm. /// https://en.wikipedia.org/wiki/Binary_logarithm#Iterative_approximation /// /// Requirements: /// - x must be greater than zero. /// /// Caveats: /// - The results are nor perfectly accurate to the last digit, due to the lossy precision of the iterative approximation. /// /// @param x The signed 59.18-decimal fixed-point number for which to calculate the binary logarithm. /// @return result The binary logarithm as a signed 59.18-decimal fixed-point number. function log2(int256 x) internal pure returns (int256 result) { require(x > 0); unchecked { // This works because log2(x) = -log2(1/x). int256 sign; if (x >= SCALE) { sign = 1; } else { sign = -1; // Do the fixed-point inversion inline to save gas. The numerator is SCALE * SCALE. assembly { x := div(1000000000000000000000000000000000000, x) } } // Calculate the integer part of the logarithm and add it to the result and finally calculate y = x * 2^(-n). uint256 n = PRBMathCommon.mostSignificantBit(uint256(x / SCALE)); // The integer part of the logarithm as a signed 59.18-decimal fixed-point number. The operation can't overflow // because n is maximum 255, SCALE is 1e18 and sign is either 1 or -1. result = int256(n) * SCALE; // This is y = x * 2^(-n). int256 y = x >> n; // If y = 1, the fractional part is zero. if (y == SCALE) { return result * sign; } // Calculate the fractional part via the iterative approximation. // The "delta >>= 1" part is equivalent to "delta /= 2", but shifting bits is faster. for (int256 delta = int256(HALF_SCALE); delta > 0; delta >>= 1) { y = (y * y) / SCALE; // Is y^2 > 2 and so in the range [2,4)? if (y >= 2 * SCALE) { // Add the 2^(-m) factor to the logarithm. result += delta; // Corresponds to z/2 on Wikipedia. y >>= 1; } } result *= sign; } } /// @notice Multiplies two signed 59.18-decimal fixed-point numbers together, returning a new signed 59.18-decimal /// fixed-point number. /// /// @dev Variant of "mulDiv" that works with signed numbers and employs constant folding, i.e. the denominator is /// alawys 1e18. /// /// Requirements: /// - All from "PRBMathCommon.mulDivFixedPoint". /// - The result must fit within MAX_SD59x18. /// /// Caveats: /// - The body is purposely left uncommented; see the NatSpec comments in "PRBMathCommon.mulDiv" to understand how this works. /// /// @param x The multiplicand as a signed 59.18-decimal fixed-point number. /// @param y The multiplier as a signed 59.18-decimal fixed-point number. /// @return result The result as a signed 59.18-decimal fixed-point number. function mul(int256 x, int256 y) internal pure returns (int256 result) { require(x > MIN_SD59x18); require(y > MIN_SD59x18); unchecked { uint256 ax; uint256 ay; ax = x < 0 ? uint256(-x) : uint256(x); ay = y < 0 ? uint256(-y) : uint256(y); uint256 resultUnsigned = PRBMathCommon.mulDivFixedPoint(ax, ay); require(resultUnsigned <= uint256(MAX_SD59x18)); uint256 sx; uint256 sy; assembly { sx := sgt(x, sub(0, 1)) sy := sgt(y, sub(0, 1)) } result = sx ^ sy == 1 ? -int256(resultUnsigned) : int256(resultUnsigned); } } /// @notice Retrieves PI as a signed 59.18-decimal fixed-point number. function pi() internal pure returns (int256 result) { result = 3141592653589793238; } /// @notice Raises x (signed 59.18-decimal fixed-point number) to the power of y (basic unsigned integer) using the /// famous algorithm "exponentiation by squaring". /// /// @dev See https://en.wikipedia.org/wiki/Exponentiation_by_squaring /// /// Requirements: /// - All from "abs" and "PRBMathCommon.mulDivFixedPoint". /// - The result must fit within MAX_SD59x18. /// /// Caveats: /// - All from "PRBMathCommon.mulDivFixedPoint". /// - Assumes 0^0 is 1. /// /// @param x The base as a signed 59.18-decimal fixed-point number. /// @param y The exponent as an uint256. /// @return result The result as a signed 59.18-decimal fixed-point number. function pow(int256 x, uint256 y) internal pure returns (int256 result) { uint256 absX = uint256(abs(x)); // Calculate the first iteration of the loop in advance. uint256 absResult = y & 1 > 0 ? absX : uint256(SCALE); // Equivalent to "for(y /= 2; y > 0; y /= 2)" but faster. for (y >>= 1; y > 0; y >>= 1) { absX = PRBMathCommon.mulDivFixedPoint(absX, absX); // Equivalent to "y % 2 == 1" but faster. if (y & 1 > 0) { absResult = PRBMathCommon.mulDivFixedPoint(absResult, absX); } } // The result must fit within the 59.18-decimal fixed-point representation. require(absResult <= uint256(MAX_SD59x18)); // Is the base negative and the exponent an odd number? bool isNegative = x < 0 && y & 1 == 1; result = isNegative ? -int256(absResult) : int256(absResult); } /// @notice Returns 1 as a signed 59.18-decimal fixed-point number. function scale() internal pure returns (int256 result) { result = SCALE; } /// @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. /// /// Requirements: /// - x cannot be negative. /// - x must be less than MAX_SD59x18 / SCALE. /// /// Caveats: /// - The maximum fixed-point number permitted is 57896044618658097711785492504343953926634.992332820282019729. /// /// @param x The signed 59.18-decimal fixed-point number for which to calculate the square root. /// @return result The result as a signed 59.18-decimal fixed-point . function sqrt(int256 x) internal pure returns (int256 result) { require(x >= 0); require(x < 57896044618658097711785492504343953926634992332820282019729); unchecked { // Multiply x by the SCALE to account for the factor of SCALE that is picked up when multiplying two signed // 59.18-decimal fixed-point numbers together (in this case, those two numbers are both the square root). result = int256(PRBMathCommon.sqrt(uint256(x * SCALE))); } } }