-
Notifications
You must be signed in to change notification settings - Fork 12k
/
Copy pathMath.sol
226 lines (201 loc) · 8.62 KB
/
Math.sol
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
// SPDX-License-Identifier: MIT
// OpenZeppelin Contracts (last updated v4.5.0) (utils/math/Math.sol)
pragma solidity ^0.8.0;
/**
* @dev Standard math utilities missing in the Solidity language.
*/
library Math {
enum Rounding {
Down, // Toward negative infinity
Up, // Toward infinity
Zero // Toward zero
}
/**
* @dev Returns the largest of two numbers.
*/
function max(uint256 a, uint256 b) internal pure returns (uint256) {
return a >= b ? a : b;
}
/**
* @dev Returns the smallest of two numbers.
*/
function min(uint256 a, uint256 b) internal pure returns (uint256) {
return a < b ? a : b;
}
/**
* @dev Returns the average of two numbers. The result is rounded towards
* zero.
*/
function average(uint256 a, uint256 b) internal pure returns (uint256) {
// (a + b) / 2 can overflow.
return (a & b) + (a ^ b) / 2;
}
/**
* @dev Returns the ceiling of the division of two numbers.
*
* This differs from standard division with `/` in that it rounds up instead
* of rounding down.
*/
function ceilDiv(uint256 a, uint256 b) internal pure returns (uint256) {
// (a + b - 1) / b can overflow on addition, so we distribute.
return a == 0 ? 0 : (a - 1) / b + 1;
}
/**
* @notice Calculates floor(x * y / denominator) with full precision. Throws if result overflows a uint256 or denominator == 0
* @dev Original credit to Remco Bloemen under MIT license (https://xn--2-umb.com/21/muldiv)
* with further edits by Uniswap Labs also under MIT license.
*/
function mulDiv(
uint256 x,
uint256 y,
uint256 denominator
) internal pure returns (uint256 result) {
unchecked {
// 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.
if (prod1 == 0) {
return prod0 / denominator;
}
// Make sure the result is less than 2^256. Also prevents denominator == 0.
require(denominator > prod1);
///////////////////////////////////////////////
// 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.
// Does not overflow because the denominator cannot be zero at this stage in the function.
uint256 twos = denominator & (~denominator + 1);
assembly {
// Divide denominator by twos.
denominator := div(denominator, twos)
// Divide [prod1 prod0] by twos.
prod0 := div(prod0, twos)
// Flip twos such that it is 2^256 / twos. If twos is zero, then it becomes one.
twos := add(div(sub(0, twos), twos), 1)
}
// Shift in bits from prod1 into prod0.
prod0 |= prod1 * twos;
// 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;
return result;
}
}
/**
* @notice Calculates x * y / denominator with full precision, following the selected rounding direction.
*/
function mulDiv(
uint256 x,
uint256 y,
uint256 denominator,
Rounding rounding
) internal pure returns (uint256) {
uint256 result = mulDiv(x, y, denominator);
if (rounding == Rounding.Up && mulmod(x, y, denominator) > 0) {
result += 1;
}
return result;
}
/**
* @dev Returns the square root of a number. It the number is not a perfect square, the value is rounded down.
*
* Inspired by Henry S. Warren, Jr.'s "Hacker's Delight" (Chapter 11).
*/
function sqrt(uint256 a) internal pure returns (uint256) {
if (a == 0) {
return 0;
}
// For our first guess, we get the biggest power of 2 which is smaller than the square root of the target.
// We know that the "msb" (most significant bit) of our target number `a` is a power of 2 such that we have
// `msb(a) <= a < 2*msb(a)`.
// We also know that `k`, the position of the most significant bit, is such that `msb(a) = 2**k`.
// This gives `2**k < a <= 2**(k+1)` → `2**(k/2) <= sqrt(a) < 2 ** (k/2+1)`.
// Using an algorithm similar to the msb conmputation, we are able to compute `result = 2**(k/2)` which is a
// good first aproximation of `sqrt(a)` with at least 1 correct bit.
uint256 result = 1;
uint256 x = a;
if (x >> 128 > 0) {
x >>= 128;
result <<= 64;
}
if (x >> 64 > 0) {
x >>= 64;
result <<= 32;
}
if (x >> 32 > 0) {
x >>= 32;
result <<= 16;
}
if (x >> 16 > 0) {
x >>= 16;
result <<= 8;
}
if (x >> 8 > 0) {
x >>= 8;
result <<= 4;
}
if (x >> 4 > 0) {
x >>= 4;
result <<= 2;
}
if (x >> 2 > 0) {
result <<= 1;
}
// At this point `result` is an estimation with one bit of precision. We know the true value is a uint128,
// since it is the square root of a uint256. Newton's method converges quadratically (precision doubles at
// every iteration). We thus need at most 7 iteration to turn our partial result with one bit of precision
// into the expected uint128 result.
unchecked {
result = (result + a / result) >> 1;
result = (result + a / result) >> 1;
result = (result + a / result) >> 1;
result = (result + a / result) >> 1;
result = (result + a / result) >> 1;
result = (result + a / result) >> 1;
result = (result + a / result) >> 1;
return min(result, a / result);
}
}
/**
* @notice Calculates sqrt(a), following the selected rounding direction.
*/
function sqrt(uint256 a, Rounding rounding) internal pure returns (uint256) {
uint256 result = sqrt(a);
if (rounding == Rounding.Up && result * result < a) {
result += 1;
}
return result;
}
}