// SPDX-License-Identifier: AGPL-3.0-only
pragma solidity >=0.8.0;
/// @notice Arithmetic library with operations for fixed-point numbers.
/// @author Solmate (https://github.com/transmissions11/solmate/blob/main/src/utils/FixedPointMathLib.sol)
/// @author Inspired by USM (https://github.com/usmfum/USM/blob/master/contracts/WadMath.sol)
library FixedPointMathLib {
/*//////////////////////////////////////////////////////////////
SIMPLIFIED FIXED POINT OPERATIONS
//////////////////////////////////////////////////////////////*/
uint256 internal constant WAD = 1e18; // The scalar of ETH and most ERC20s.
function mulWadDown(uint256 x, uint256 y) internal pure returns (uint256) {
return mulDivDown(x, y, WAD); // Equivalent to (x * y) / WAD rounded down.
}
function mulWadUp(uint256 x, uint256 y) internal pure returns (uint256) {
return mulDivUp(x, y, WAD); // Equivalent to (x * y) / WAD rounded up.
}
function divWadDown(uint256 x, uint256 y) internal pure returns (uint256) {
return mulDivDown(x, WAD, y); // Equivalent to (x * WAD) / y rounded down.
}
function divWadUp(uint256 x, uint256 y) internal pure returns (uint256) {
return mulDivUp(x, WAD, y); // Equivalent to (x * WAD) / y rounded up.
}
/*//////////////////////////////////////////////////////////////
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)
}
}
function mulDivUp(
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)
}
// First, divide z - 1 by the denominator and add 1.
// We allow z - 1 to underflow if z is 0, because we multiply the
// end result by 0 if z is zero, ensuring we return 0 if z is zero.
z := mul(iszero(iszero(z)), add(div(sub(z, 1), denominator), 1))
}
}
function rpow(
uint256 x,
uint256 n,
uint256 scalar
) internal pure returns (uint256 z) {
assembly {
switch x
case 0 {
switch n
case 0 {
// 0 ** 0 = 1
z := scalar
}
default {
// 0 ** n = 0
z := 0
}
}
default {
switch mod(n, 2)
case 0 {
// If n is even, store scalar in z for now.
z := scalar
}
default {
// If n is odd, store x in z for now.
z := x
}
// Shifting right by 1 is like dividing by 2.
let half := shr(1, scalar)
for {
// Shift n right by 1 before looping to halve it.
n := shr(1, n)
} n {
// Shift n right by 1 each iteration to halve it.
n := shr(1, n)
} {
// Revert immediately if x ** 2 would overflow.
// Equivalent to iszero(eq(div(xx, x), x)) here.
if shr(128, x) {
revert(0, 0)
}
// Store x squared.
let xx := mul(x, x)
// Round to the nearest number.
let xxRound := add(xx, half)
// Revert if xx + half overflowed.
if lt(xxRound, xx) {
revert(0, 0)
}
// Set x to scaled xxRound.
x := div(xxRound, scalar)
// If n is even:
if mod(n, 2) {
// Compute z * x.
let zx := mul(z, x)
// If z * x overflowed:
if iszero(eq(div(zx, x), z)) {
// Revert if x is non-zero.
if iszero(iszero(x)) {
revert(0, 0)
}
}
// Round to the nearest number.
let zxRound := add(zx, half)
// Revert if zx + half overflowed.
if lt(zxRound, zx) {
revert(0, 0)
}
// Return properly scaled zxRound.
z := div(zxRound, scalar)
}
}
}
}
}
/*//////////////////////////////////////////////////////////////
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))
}
}
function unsafeMod(uint256 x, uint256 y) internal pure returns (uint256 z) {
assembly {
// Mod x by y. Note this will return
// 0 instead of reverting if y is zero.
z := mod(x, y)
}
}
function unsafeDiv(uint256 x, uint256 y) internal pure returns (uint256 r) {
assembly {
// Divide x by y. Note this will return
// 0 instead of reverting if y is zero.
r := div(x, y)
}
}
function unsafeDivUp(uint256 x, uint256 y) internal pure returns (uint256 z) {
assembly {
// Add 1 to x * y if x % y > 0. Note this will
// return 0 instead of reverting if y is zero.
z := add(gt(mod(x, y), 0), div(x, y))
}
}
}