115 lines
5.1 KiB
Solidity
115 lines
5.1 KiB
Solidity
// SPDX-License-Identifier: MIT
|
|
pragma solidity >=0.8.4;
|
|
|
|
/// @notice Arithmetic library with operations for fixed-point numbers.
|
|
/// @author Solmate (https://github.com/Rari-Capital/solmate/blob/main/src/utils/FixedPointMathLib.sol)
|
|
library FixedPointMathLib {
|
|
|
|
/*//////////////////////////////////////////////////////////////
|
|
SIMPLIFIED FIXED POINT OPERATIONS
|
|
//////////////////////////////////////////////////////////////*/
|
|
|
|
/// @dev The scalar of ETH and most ERC20s.
|
|
uint256 internal constant WAD = 1e18;
|
|
|
|
function mulWadDown(uint256 x, uint256 y) internal pure returns (uint256) {
|
|
// Equivalent to (x * y) / WAD rounded down.
|
|
return mulDivDown(x, y, WAD);
|
|
}
|
|
|
|
function divWadDown(uint256 x, uint256 y) internal pure returns (uint256) {
|
|
// Equivalent to (x * WAD) / y rounded down.
|
|
return mulDivDown(x, WAD, y);
|
|
}
|
|
|
|
/*//////////////////////////////////////////////////////////////
|
|
LOW LEVEL FIXED POINT OPERATIONS
|
|
//////////////////////////////////////////////////////////////*/
|
|
|
|
function mulDivDown(
|
|
uint256 x,
|
|
uint256 y,
|
|
uint256 denominator
|
|
) internal pure returns (uint256 z) {
|
|
assembly {
|
|
// Store x * y in z for now.
|
|
z := mul(x, y)
|
|
|
|
// Equivalent to require(denominator != 0 && (x == 0 || (x * y) / x == y))
|
|
if iszero(and(iszero(iszero(denominator)), or(iszero(x), eq(div(z, x), y)))) {
|
|
revert(0, 0)
|
|
}
|
|
|
|
// Divide z by the denominator.
|
|
z := div(z, denominator)
|
|
}
|
|
}
|
|
|
|
/*//////////////////////////////////////////////////////////////
|
|
GENERAL NUMBER UTILITIES
|
|
//////////////////////////////////////////////////////////////*/
|
|
|
|
function sqrt(uint256 x) internal pure returns (uint256 z) {
|
|
assembly {
|
|
let y := x // We start y at x, which will help us make our initial estimate.
|
|
|
|
z := 181 // The "correct" value is 1, but this saves a multiplication later.
|
|
|
|
// This segment is to get a reasonable initial estimate for the Babylonian method. With a bad
|
|
// start, the correct # of bits increases ~linearly each iteration instead of ~quadratically.
|
|
|
|
// We check y >= 2^(k + 8) but shift right by k bits
|
|
// each branch to ensure that if x >= 256, then y >= 256.
|
|
if iszero(lt(y, 0x10000000000000000000000000000000000)) {
|
|
y := shr(128, y)
|
|
z := shl(64, z)
|
|
}
|
|
if iszero(lt(y, 0x1000000000000000000)) {
|
|
y := shr(64, y)
|
|
z := shl(32, z)
|
|
}
|
|
if iszero(lt(y, 0x10000000000)) {
|
|
y := shr(32, y)
|
|
z := shl(16, z)
|
|
}
|
|
if iszero(lt(y, 0x1000000)) {
|
|
y := shr(16, y)
|
|
z := shl(8, z)
|
|
}
|
|
|
|
// Goal was to get z*z*y within a small factor of x. More iterations could
|
|
// get y in a tighter range. Currently, we will have y in [256, 256*2^16).
|
|
// We ensured y >= 256 so that the relative difference between y and y+1 is small.
|
|
// That's not possible if x < 256 but we can just verify those cases exhaustively.
|
|
|
|
// Now, z*z*y <= x < z*z*(y+1), and y <= 2^(16+8), and either y >= 256, or x < 256.
|
|
// Correctness can be checked exhaustively for x < 256, so we assume y >= 256.
|
|
// Then z*sqrt(y) is within sqrt(257)/sqrt(256) of sqrt(x), or about 20bps.
|
|
|
|
// For s in the range [1/256, 256], the estimate f(s) = (181/1024) * (s+1) is in the range
|
|
// (1/2.84 * sqrt(s), 2.84 * sqrt(s)), with largest error when s = 1 and when s = 256 or 1/256.
|
|
|
|
// Since y is in [256, 256*2^16), let a = y/65536, so that a is in [1/256, 256). Then we can estimate
|
|
// sqrt(y) using sqrt(65536) * 181/1024 * (a + 1) = 181/4 * (y + 65536)/65536 = 181 * (y + 65536)/2^18.
|
|
|
|
// There is no overflow risk here since y < 2^136 after the first branch above.
|
|
z := shr(18, mul(z, add(y, 65536))) // A mul() is saved from starting z at 181.
|
|
|
|
// Given the worst case multiplicative error of 2.84 above, 7 iterations should be enough.
|
|
z := shr(1, add(z, div(x, z)))
|
|
z := shr(1, add(z, div(x, z)))
|
|
z := shr(1, add(z, div(x, z)))
|
|
z := shr(1, add(z, div(x, z)))
|
|
z := shr(1, add(z, div(x, z)))
|
|
z := shr(1, add(z, div(x, z)))
|
|
z := shr(1, add(z, div(x, z)))
|
|
|
|
// If x+1 is a perfect square, the Babylonian method cycles between
|
|
// floor(sqrt(x)) and ceil(sqrt(x)). This statement ensures we return floor.
|
|
// See: https://en.wikipedia.org/wiki/Integer_square_root#Using_only_integer_division
|
|
// Since the ceil is rare, we save gas on the assignment and repeat division in the rare case.
|
|
// If you don't care whether the floor or ceil square root is returned, you can remove this statement.
|
|
z := sub(z, lt(div(x, z), z))
|
|
}
|
|
}
|
|
} |