forked from OpenZeppelin/openzeppelin-contracts
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Math.sol
573 lines (508 loc) · 22.8 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
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
// SPDX-License-Identifier: MIT
// OpenZeppelin Contracts (last updated v5.0.0) (utils/math/Math.sol)
pragma solidity ^0.8.20;
import {Address} from "../Address.sol";
import {Panic} from "../Panic.sol";
import {SafeCast} from "./SafeCast.sol";
/**
* @dev Standard math utilities missing in the Solidity language.
*/
library Math {
enum Rounding {
Floor, // Toward negative infinity
Ceil, // Toward positive infinity
Trunc, // Toward zero
Expand // Away from zero
}
/**
* @dev Returns the addition of two unsigned integers, with an success flag (no overflow).
*/
function tryAdd(uint256 a, uint256 b) internal pure returns (bool success, uint256 result) {
unchecked {
uint256 c = a + b;
if (c < a) return (false, 0);
return (true, c);
}
}
/**
* @dev Returns the subtraction of two unsigned integers, with an success flag (no overflow).
*/
function trySub(uint256 a, uint256 b) internal pure returns (bool success, uint256 result) {
unchecked {
if (b > a) return (false, 0);
return (true, a - b);
}
}
/**
* @dev Returns the multiplication of two unsigned integers, with an success flag (no overflow).
*/
function tryMul(uint256 a, uint256 b) internal pure returns (bool success, uint256 result) {
unchecked {
// Gas optimization: this is cheaper than requiring 'a' not being zero, but the
// benefit is lost if 'b' is also tested.
// See: https://github.com/OpenZeppelin/openzeppelin-contracts/pull/522
if (a == 0) return (true, 0);
uint256 c = a * b;
if (c / a != b) return (false, 0);
return (true, c);
}
}
/**
* @dev Returns the division of two unsigned integers, with a success flag (no division by zero).
*/
function tryDiv(uint256 a, uint256 b) internal pure returns (bool success, uint256 result) {
unchecked {
if (b == 0) return (false, 0);
return (true, a / b);
}
}
/**
* @dev Returns the remainder of dividing two unsigned integers, with a success flag (no division by zero).
*/
function tryMod(uint256 a, uint256 b) internal pure returns (bool success, uint256 result) {
unchecked {
if (b == 0) return (false, 0);
return (true, a % b);
}
}
/**
* @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 towards infinity instead
* of rounding towards zero.
*/
function ceilDiv(uint256 a, uint256 b) internal pure returns (uint256) {
if (b == 0) {
// Guarantee the same behavior as in a regular Solidity division.
Panic.panic(Panic.DIVISION_BY_ZERO);
}
// The following calculation ensures accurate ceiling division without overflow.
// Since a is non-zero, (a - 1) / b will not overflow.
// The largest possible result occurs when (a - 1) / b is type(uint256).max,
// but the largest value we can obtain is type(uint256).max - 1, which happens
// when a = type(uint256).max and b = 1.
unchecked {
return a == 0 ? 0 : (a - 1) / b + 1;
}
}
/**
* @dev Calculates floor(x * y / denominator) with full precision. Throws if result overflows a uint256 or
* denominator == 0.
*
* 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 = x * y; // Least significant 256 bits of the product
uint256 prod1; // Most significant 256 bits of the product
assembly {
let mm := mulmod(x, y, not(0))
prod1 := sub(sub(mm, prod0), lt(mm, prod0))
}
// Handle non-overflow cases, 256 by 256 division.
if (prod1 == 0) {
// Solidity will revert if denominator == 0, unlike the div opcode on its own.
// The surrounding unchecked block does not change this fact.
// See https://docs.soliditylang.org/en/latest/control-structures.html#checked-or-unchecked-arithmetic.
return prod0 / denominator;
}
// Make sure the result is less than 2^256. Also prevents denominator == 0.
if (denominator <= prod1) {
Panic.panic(denominator == 0 ? Panic.DIVISION_BY_ZERO : Panic.UNDER_OVERFLOW);
}
///////////////////////////////////////////////
// 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.
uint256 twos = denominator & (0 - denominator);
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;
}
}
/**
* @dev 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) {
return mulDiv(x, y, denominator) + SafeCast.toUint(unsignedRoundsUp(rounding) && mulmod(x, y, denominator) > 0);
}
/**
* @dev Calculate the modular multiplicative inverse of a number in Z/nZ.
*
* If n is a prime, then Z/nZ is a field. In that case all elements are inversible, expect 0.
* If n is not a prime, then Z/nZ is not a field, and some elements might not be inversible.
*
* If the input value is not inversible, 0 is returned.
*
* NOTE: If you know for sure that n is (big) a prime, it may be cheaper to use Ferma's little theorem and get the
* inverse using `Math.modExp(a, n - 2, n)`.
*/
function invMod(uint256 a, uint256 n) internal pure returns (uint256) {
unchecked {
if (n == 0) return 0;
// The inverse modulo is calculated using the Extended Euclidean Algorithm (iterative version)
// Used to compute integers x and y such that: ax + ny = gcd(a, n).
// When the gcd is 1, then the inverse of a modulo n exists and it's x.
// ax + ny = 1
// ax = 1 + (-y)n
// ax ≡ 1 (mod n) # x is the inverse of a modulo n
// If the remainder is 0 the gcd is n right away.
uint256 remainder = a % n;
uint256 gcd = n;
// Therefore the initial coefficients are:
// ax + ny = gcd(a, n) = n
// 0a + 1n = n
int256 x = 0;
int256 y = 1;
while (remainder != 0) {
uint256 quotient = gcd / remainder;
(gcd, remainder) = (
// The old remainder is the next gcd to try.
remainder,
// Compute the next remainder.
// Can't overflow given that (a % gcd) * (gcd // (a % gcd)) <= gcd
// where gcd is at most n (capped to type(uint256).max)
gcd - remainder * quotient
);
(x, y) = (
// Increment the coefficient of a.
y,
// Decrement the coefficient of n.
// Can overflow, but the result is casted to uint256 so that the
// next value of y is "wrapped around" to a value between 0 and n - 1.
x - y * int256(quotient)
);
}
if (gcd != 1) return 0; // No inverse exists.
return x < 0 ? (n - uint256(-x)) : uint256(x); // Wrap the result if it's negative.
}
}
/**
* @dev Returns the modular exponentiation of the specified base, exponent and modulus (b ** e % m)
*
* Requirements:
* - modulus can't be zero
* - underlying staticcall to precompile must succeed
*
* IMPORTANT: The result is only valid if the underlying call succeeds. When using this function, make
* sure the chain you're using it on supports the precompiled contract for modular exponentiation
* at address 0x05 as specified in https://eips.ethereum.org/EIPS/eip-198[EIP-198]. Otherwise,
* the underlying function will succeed given the lack of a revert, but the result may be incorrectly
* interpreted as 0.
*/
function modExp(uint256 b, uint256 e, uint256 m) internal view returns (uint256) {
(bool success, uint256 result) = tryModExp(b, e, m);
if (!success) {
if (m == 0) {
Panic.panic(Panic.DIVISION_BY_ZERO);
} else {
revert Address.FailedInnerCall();
}
}
return result;
}
/**
* @dev Returns the modular exponentiation of the specified base, exponent and modulus (b ** e % m).
* It includes a success flag indicating if the operation succeeded. Operation will be marked has failed if trying
* to operate modulo 0 or if the underlying precompile reverted.
*
* IMPORTANT: The result is only valid if the success flag is true. When using this function, make sure the chain
* you're using it on supports the precompiled contract for modular exponentiation at address 0x05 as specified in
* https://eips.ethereum.org/EIPS/eip-198[EIP-198]. Otherwise, the underlying function will succeed given the lack
* of a revert, but the result may be incorrectly interpreted as 0.
*/
function tryModExp(uint256 b, uint256 e, uint256 m) internal view returns (bool success, uint256 result) {
if (m == 0) return (false, 0);
/// @solidity memory-safe-assembly
assembly {
let ptr := mload(0x40)
// | Offset | Content | Content (Hex) |
// |-----------|------------|--------------------------------------------------------------------|
// | 0x00:0x1f | size of b | 0x0000000000000000000000000000000000000000000000000000000000000020 |
// | 0x20:0x3f | size of e | 0x0000000000000000000000000000000000000000000000000000000000000020 |
// | 0x40:0x5f | size of m | 0x0000000000000000000000000000000000000000000000000000000000000020 |
// | 0x60:0x7f | value of b | 0x<.............................................................b> |
// | 0x80:0x9f | value of e | 0x<.............................................................e> |
// | 0xa0:0xbf | value of m | 0x<.............................................................m> |
mstore(ptr, 0x20)
mstore(add(ptr, 0x20), 0x20)
mstore(add(ptr, 0x40), 0x20)
mstore(add(ptr, 0x60), b)
mstore(add(ptr, 0x80), e)
mstore(add(ptr, 0xa0), m)
// Given the result < m, it's guaranteed to fit in 32 bytes,
// so we can use the memory scratch space located at offset 0.
success := staticcall(gas(), 0x05, ptr, 0xc0, 0x00, 0x20)
result := mload(0x00)
}
}
/**
* @dev Returns the square root of a number. If the number is not a perfect square, the value is rounded
* towards zero.
*
* This method is based on Newton's method for computing square roots; the algorithm is restricted to only
* using integer operations.
*/
function sqrt(uint256 a) internal pure returns (uint256) {
unchecked {
// Take care of easy edge cases when a == 0 or a == 1
if (a <= 1) {
return a;
}
uint256 aAux = a;
uint256 result = 1;
// For our first guess, we get the biggest power of 2 which is smaller
// than the square root of the target (i.e. result = 2**n <= sqrt(a) < 2**(n+1)).
// We know if result is 2**e + c, then e is bounded to 127 because (2**128)**2 = 2**256,
// which is bigger than any uint256. If result >= 2**e, then sqrt(a) <= 2**e-1.
// We approximate the result by adding 2**e-1 and subtracting 2**e for each e such that
// sqrt(a) < 2**e by cutting e/2 on each iteration until e = 2.
if (aAux >= (1 << 128)) {
aAux >>= 128;
result <<= 64; // sqrt(a) >= 2**(e/2)
}
if (aAux >= (1 << 64)) {
aAux >>= 64;
result <<= 32; // sqrt(a) >= 2**(e/2)
}
if (aAux >= (1 << 32)) {
aAux >>= 32;
result <<= 16; // sqrt(a) >= 2**(e/2)
}
if (aAux >= (1 << 16)) {
aAux >>= 16;
result <<= 8; // sqrt(a) >= 2**(e/2)
}
if (aAux >= (1 << 8)) {
aAux >>= 8;
result <<= 4; // sqrt(a) >= 2**(e/2)
}
if (aAux >= (1 << 4)) {
aAux >>= 4;
result <<= 2; // sqrt(a) >= 2**(e/2)
}
if (aAux >= (1 << 2)) {
result <<= 1;
}
// We know that result <= sqrt(a) < 2 * result. We can use the fact that
// 2**e <= sqrt(a) to improve the estimation by computing the arithmetic
// mean between the current estimation and the next one (e = 1).
result = (3 * result) >> 1;
// We define the error as ε = result - sqrt(a). Then we know that
// result = 2**e−1 + 2**e−2, and therefore ε0 = 2**e−1 + 2**e−2 - sqrt(n),
// leaving ε_0 <= 2**e−2. We also see ε_{+1} == ε**2 / 2x <= ε**2 / 2 * sqrt(a)
// as shown in Walter Rudin. Principles of Mathematical Analysis.
// 3rd ed. McGraw-Hill New York, 1976. Exercise 3.16 (b)
result = (result + a / result) >> 1; // ε_1 := result - sqrt(a) <= 2**(e-4.5)
result = (result + a / result) >> 1; // ε_2 := result - sqrt(a) <= 2**(e-9)
result = (result + a / result) >> 1; // ε_3 := result - sqrt(a) <= 2**(e-18)
result = (result + a / result) >> 1; // ε_4 := result - sqrt(a) <= 2**(e-36)
result = (result + a / result) >> 1; // ε_5 := result - sqrt(a) <= 2**(e-72)
result = (result + a / result) >> 1; // ε_6 := result - sqrt(a) <= 2**(e-144)
// After 6 iterations, the precision of e is already above 128 (i.e. 144). Meaning that
// ε6 <= 1. And given we're operating on integers, then we can ensure that result is
// either sqrt(a) or sqrt(a) + 1.
return result - SafeCast.toUint(result > a / result);
}
}
/**
* @dev Calculates sqrt(a), following the selected rounding direction.
*/
function sqrt(uint256 a, Rounding rounding) internal pure returns (uint256) {
unchecked {
uint256 result = sqrt(a);
return result + SafeCast.toUint(unsignedRoundsUp(rounding) && result * result < a);
}
}
/**
* @dev Return the log in base 2 of a positive value rounded towards zero.
* Returns 0 if given 0.
*/
function log2(uint256 value) internal pure returns (uint256) {
uint256 result = 0;
uint256 exp;
unchecked {
exp = 128 * SafeCast.toUint(value > (1 << 128) - 1);
value >>= exp;
result += exp;
exp = 64 * SafeCast.toUint(value > (1 << 64) - 1);
value >>= exp;
result += exp;
exp = 32 * SafeCast.toUint(value > (1 << 32) - 1);
value >>= exp;
result += exp;
exp = 16 * SafeCast.toUint(value > (1 << 16) - 1);
value >>= exp;
result += exp;
exp = 8 * SafeCast.toUint(value > (1 << 8) - 1);
value >>= exp;
result += exp;
exp = 4 * SafeCast.toUint(value > (1 << 4) - 1);
value >>= exp;
result += exp;
exp = 2 * SafeCast.toUint(value > (1 << 2) - 1);
value >>= exp;
result += exp;
result += SafeCast.toUint(value > 1);
}
return result;
}
/**
* @dev Return the log in base 2, following the selected rounding direction, of a positive value.
* Returns 0 if given 0.
*/
function log2(uint256 value, Rounding rounding) internal pure returns (uint256) {
unchecked {
uint256 result = log2(value);
return result + SafeCast.toUint(unsignedRoundsUp(rounding) && 1 << result < value);
}
}
/**
* @dev Return the log in base 10 of a positive value rounded towards zero.
* Returns 0 if given 0.
*/
function log10(uint256 value) internal pure returns (uint256) {
uint256 result = 0;
unchecked {
if (value >= 10 ** 64) {
value /= 10 ** 64;
result += 64;
}
if (value >= 10 ** 32) {
value /= 10 ** 32;
result += 32;
}
if (value >= 10 ** 16) {
value /= 10 ** 16;
result += 16;
}
if (value >= 10 ** 8) {
value /= 10 ** 8;
result += 8;
}
if (value >= 10 ** 4) {
value /= 10 ** 4;
result += 4;
}
if (value >= 10 ** 2) {
value /= 10 ** 2;
result += 2;
}
if (value >= 10 ** 1) {
result += 1;
}
}
return result;
}
/**
* @dev Return the log in base 10, following the selected rounding direction, of a positive value.
* Returns 0 if given 0.
*/
function log10(uint256 value, Rounding rounding) internal pure returns (uint256) {
unchecked {
uint256 result = log10(value);
return result + SafeCast.toUint(unsignedRoundsUp(rounding) && 10 ** result < value);
}
}
/**
* @dev Return the log in base 256 of a positive value rounded towards zero.
* Returns 0 if given 0.
*
* Adding one to the result gives the number of pairs of hex symbols needed to represent `value` as a hex string.
*/
function log256(uint256 value) internal pure returns (uint256) {
uint256 result = 0;
uint256 isGt;
unchecked {
isGt = SafeCast.toUint(value > (1 << 128) - 1);
value >>= isGt * 128;
result += isGt * 16;
isGt = SafeCast.toUint(value > (1 << 64) - 1);
value >>= isGt * 64;
result += isGt * 8;
isGt = SafeCast.toUint(value > (1 << 32) - 1);
value >>= isGt * 32;
result += isGt * 4;
isGt = SafeCast.toUint(value > (1 << 16) - 1);
value >>= isGt * 16;
result += isGt * 2;
result += SafeCast.toUint(value > (1 << 8) - 1);
}
return result;
}
/**
* @dev Return the log in base 256, following the selected rounding direction, of a positive value.
* Returns 0 if given 0.
*/
function log256(uint256 value, Rounding rounding) internal pure returns (uint256) {
unchecked {
uint256 result = log256(value);
return result + SafeCast.toUint(unsignedRoundsUp(rounding) && 1 << (result << 3) < value);
}
}
/**
* @dev Returns whether a provided rounding mode is considered rounding up for unsigned integers.
*/
function unsignedRoundsUp(Rounding rounding) internal pure returns (bool) {
return uint8(rounding) % 2 == 1;
}
}