-
-
Notifications
You must be signed in to change notification settings - Fork 122
/
Common.sol
647 lines (596 loc) · 23.7 KB
/
Common.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
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
// SPDX-License-Identifier: MIT
pragma solidity >=0.8.13;
/// Common mathematical functions used in both SD59x18 and UD60x18. Note that these global functions do not
/// always operate with SD59x18 and UD60x18 numbers.
/*//////////////////////////////////////////////////////////////////////////
CUSTOM ERRORS
//////////////////////////////////////////////////////////////////////////*/
/// @notice Emitted when the ending result in the fixed-point version of `mulDiv` would overflow uint256.
error PRBMath_MulDiv18_Overflow(uint256 x, uint256 y);
/// @notice Emitted when the ending result in `mulDiv` would overflow uint256.
error PRBMath_MulDiv_Overflow(uint256 x, uint256 y, uint256 denominator);
/// @notice Emitted when attempting to run `mulDiv` with one of the inputs `type(int256).min`.
error PRBMath_MulDivSigned_InputTooSmall();
/// @notice Emitted when the ending result in the signed version of `mulDiv` would overflow int256.
error PRBMath_MulDivSigned_Overflow(int256 x, int256 y);
/*//////////////////////////////////////////////////////////////////////////
CONSTANTS
//////////////////////////////////////////////////////////////////////////*/
/// @dev The maximum value an uint128 number can have.
uint128 constant MAX_UINT128 = type(uint128).max;
/// @dev The maximum value an uint40 number can have.
uint40 constant MAX_UINT40 = type(uint40).max;
/// @dev How many trailing decimals can be represented.
uint256 constant UNIT = 1e18;
/// @dev Largest power of two that is a divisor of `UNIT`.
uint256 constant UNIT_LPOTD = 262144;
/// @dev The `UNIT` number inverted mod 2^256.
uint256 constant UNIT_INVERSE = 78156646155174841979727994598816262306175212592076161876661_508869554232690281;
/*//////////////////////////////////////////////////////////////////////////
FUNCTIONS
//////////////////////////////////////////////////////////////////////////*/
/// @notice Finds the zero-based index of the first one in the binary representation of x.
/// @dev See the note on msb in the "Find First Set" Wikipedia article https://en.wikipedia.org/wiki/Find_first_set
///
/// Each of the steps in this implementation is equivalent to this high-level code:
///
/// ```solidity
/// if (x >= 2 ** 128) {
/// x >>= 128;
/// result += 128;
/// }
/// ```
///
/// Where 128 is swapped with each respective power of two factor. See the full high-level implementation here:
/// https://gist.github.com/PaulRBerg/f932f8693f2733e30c4d479e8e980948
///
/// A list of the Yul instructions used below:
/// - "gt" is "greater than"
/// - "or" is the OR bitwise operator
/// - "shl" is "shift left"
/// - "shr" is "shift right"
///
/// @param x The uint256 number for which to find the index of the most significant bit.
/// @return result The index of the most significant bit as an uint256.
function msb(uint256 x) pure returns (uint256 result) {
// 2^128
assembly ("memory-safe") {
let factor := shl(7, gt(x, 0xFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF))
x := shr(factor, x)
result := or(result, factor)
}
// 2^64
assembly ("memory-safe") {
let factor := shl(6, gt(x, 0xFFFFFFFFFFFFFFFF))
x := shr(factor, x)
result := or(result, factor)
}
// 2^32
assembly ("memory-safe") {
let factor := shl(5, gt(x, 0xFFFFFFFF))
x := shr(factor, x)
result := or(result, factor)
}
// 2^16
assembly ("memory-safe") {
let factor := shl(4, gt(x, 0xFFFF))
x := shr(factor, x)
result := or(result, factor)
}
// 2^8
assembly ("memory-safe") {
let factor := shl(3, gt(x, 0xFF))
x := shr(factor, x)
result := or(result, factor)
}
// 2^4
assembly ("memory-safe") {
let factor := shl(2, gt(x, 0xF))
x := shr(factor, x)
result := or(result, factor)
}
// 2^2
assembly ("memory-safe") {
let factor := shl(1, gt(x, 0x3))
x := shr(factor, x)
result := or(result, factor)
}
// 2^1
// No need to shift x any more.
assembly ("memory-safe") {
let factor := gt(x, 0x1)
result := or(result, factor)
}
}
/// @notice Calculates floor(x*y÷denominator) with full precision.
///
/// @dev Credits to Remco Bloemen under MIT license https://xn--2-umb.com/21/muldiv.
///
/// Requirements:
/// - The denominator cannot be zero.
/// - The result must fit within uint256.
///
/// Caveats:
/// - This function does not work with fixed-point numbers.
///
/// @param x The multiplicand as an uint256.
/// @param y The multiplier as an uint256.
/// @param denominator The divisor as an uint256.
/// @return result The result as an uint256.
function mulDiv(uint256 x, uint256 y, uint256 denominator) pure returns (uint256 result) {
// 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 ("memory-safe") {
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) {
unchecked {
return prod0 / denominator;
}
}
// Make sure the result is less than 2^256. Also prevents denominator == 0.
if (prod1 >= denominator) {
revert PRBMath_MulDiv_Overflow(x, y, denominator);
}
///////////////////////////////////////////////
// 512 by 256 division.
///////////////////////////////////////////////
// Make division exact by subtracting the remainder from [prod1 prod0].
uint256 remainder;
assembly ("memory-safe") {
// Compute remainder using the mulmod Yul instruction.
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.
unchecked {
// Does not overflow because the denominator cannot be zero at this stage in the function.
uint256 lpotdod = denominator & (~denominator + 1);
assembly ("memory-safe") {
// Divide denominator by lpotdod.
denominator := div(denominator, lpotdod)
// Divide [prod1 prod0] by lpotdod.
prod0 := div(prod0, lpotdod)
// Flip lpotdod such that it is 2^256 / lpotdod. If lpotdod is zero, then it becomes one.
lpotdod := add(div(sub(0, lpotdod), lpotdod), 1)
}
// Shift in bits from prod1 into prod0.
prod0 |= prod1 * lpotdod;
// 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;
}
}
/// @notice Calculates floor(x*y÷1e18) with full precision.
///
/// @dev Variant of `mulDiv` with constant folding, i.e. in which the denominator is always 1e18. Before returning the
/// final result, we add 1 if `(x * y) % UNIT >= HALF_UNIT`. Without this adjustment, 6.6e-19 would be truncated to 0
/// instead of being rounded to 1e-18. See "Listing 6" and text above it at https://accu.org/index.php/journals/1717.
///
/// Requirements:
/// - The result must fit within uint256.
///
/// Caveats:
/// - The body is purposely left uncommented; to understand how this works, see the NatSpec comments in `mulDiv`.
/// - It is assumed that the result can never be `type(uint256).max` when x and y solve the following two equations:
/// 1. x * y = type(uint256).max * UNIT
/// 2. (x * y) % UNIT >= UNIT / 2
///
/// @param x The multiplicand as an unsigned 60.18-decimal fixed-point number.
/// @param y The multiplier as an unsigned 60.18-decimal fixed-point number.
/// @return result The result as an unsigned 60.18-decimal fixed-point number.
function mulDiv18(uint256 x, uint256 y) pure returns (uint256 result) {
uint256 prod0;
uint256 prod1;
assembly ("memory-safe") {
let mm := mulmod(x, y, not(0))
prod0 := mul(x, y)
prod1 := sub(sub(mm, prod0), lt(mm, prod0))
}
if (prod1 >= UNIT) {
revert PRBMath_MulDiv18_Overflow(x, y);
}
uint256 remainder;
assembly ("memory-safe") {
remainder := mulmod(x, y, UNIT)
}
if (prod1 == 0) {
unchecked {
return prod0 / UNIT;
}
}
assembly ("memory-safe") {
result := mul(
or(
div(sub(prod0, remainder), UNIT_LPOTD),
mul(sub(prod1, gt(remainder, prod0)), add(div(sub(0, UNIT_LPOTD), UNIT_LPOTD), 1))
),
UNIT_INVERSE
)
}
}
/// @notice Calculates floor(x*y÷denominator) with full precision.
///
/// @dev An extension of `mulDiv` for signed numbers. Works by computing the signs and the absolute values separately.
///
/// Requirements:
/// - None of the inputs can be `type(int256).min`.
/// - The result must fit within int256.
///
/// @param x The multiplicand as an int256.
/// @param y The multiplier as an int256.
/// @param denominator The divisor as an int256.
/// @return result The result as an int256.
function mulDivSigned(int256 x, int256 y, int256 denominator) pure returns (int256 result) {
if (x == type(int256).min || y == type(int256).min || denominator == type(int256).min) {
revert PRBMath_MulDivSigned_InputTooSmall();
}
// Get hold of the absolute values of x, y and the denominator.
uint256 absX;
uint256 absY;
uint256 absD;
unchecked {
absX = x < 0 ? uint256(-x) : uint256(x);
absY = y < 0 ? uint256(-y) : uint256(y);
absD = denominator < 0 ? uint256(-denominator) : uint256(denominator);
}
// Compute the absolute value of (x*y)÷denominator. The result must fit within int256.
uint256 rAbs = mulDiv(absX, absY, absD);
if (rAbs > uint256(type(int256).max)) {
revert PRBMath_MulDivSigned_Overflow(x, y);
}
// Get the signs of x, y and the denominator.
uint256 sx;
uint256 sy;
uint256 sd;
assembly ("memory-safe") {
// This works thanks to two's complement.
// "sgt" stands for "signed greater than" and "sub(0,1)" is max uint256.
sx := sgt(x, sub(0, 1))
sy := sgt(y, sub(0, 1))
sd := sgt(denominator, sub(0, 1))
}
// XOR over sx, sy and sd. What this does is to check whether there are 1 or 3 negative signs in the inputs.
// If there are, the result should be negative. Otherwise, it should be positive.
unchecked {
result = sx ^ sy ^ sd == 0 ? -int256(rAbs) : int256(rAbs);
}
}
/// @notice Calculates the binary exponent of x using the binary fraction method.
/// @dev Has to use 192.64-bit fixed-point numbers.
/// See https://ethereum.stackexchange.com/a/96594/24693.
/// @param x The exponent as an unsigned 192.64-bit fixed-point number.
/// @return result The result as an unsigned 60.18-decimal fixed-point number.
function prbExp2(uint256 x) pure returns (uint256 result) {
unchecked {
// Start from 0.5 in the 192.64-bit fixed-point format.
result = 0x800000000000000000000000000000000000000000000000;
// Multiply the result by root(2, 2^-i) when the bit at position i is 1. None of the intermediary results overflows
// because the initial result is 2^191 and all magic factors are less than 2^65.
if (x & 0xFF00000000000000 > 0) {
if (x & 0x8000000000000000 > 0) {
result = (result * 0x16A09E667F3BCC909) >> 64;
}
if (x & 0x4000000000000000 > 0) {
result = (result * 0x1306FE0A31B7152DF) >> 64;
}
if (x & 0x2000000000000000 > 0) {
result = (result * 0x1172B83C7D517ADCE) >> 64;
}
if (x & 0x1000000000000000 > 0) {
result = (result * 0x10B5586CF9890F62A) >> 64;
}
if (x & 0x800000000000000 > 0) {
result = (result * 0x1059B0D31585743AE) >> 64;
}
if (x & 0x400000000000000 > 0) {
result = (result * 0x102C9A3E778060EE7) >> 64;
}
if (x & 0x200000000000000 > 0) {
result = (result * 0x10163DA9FB33356D8) >> 64;
}
if (x & 0x100000000000000 > 0) {
result = (result * 0x100B1AFA5ABCBED61) >> 64;
}
}
if (x & 0xFF000000000000 > 0) {
if (x & 0x80000000000000 > 0) {
result = (result * 0x10058C86DA1C09EA2) >> 64;
}
if (x & 0x40000000000000 > 0) {
result = (result * 0x1002C605E2E8CEC50) >> 64;
}
if (x & 0x20000000000000 > 0) {
result = (result * 0x100162F3904051FA1) >> 64;
}
if (x & 0x10000000000000 > 0) {
result = (result * 0x1000B175EFFDC76BA) >> 64;
}
if (x & 0x8000000000000 > 0) {
result = (result * 0x100058BA01FB9F96D) >> 64;
}
if (x & 0x4000000000000 > 0) {
result = (result * 0x10002C5CC37DA9492) >> 64;
}
if (x & 0x2000000000000 > 0) {
result = (result * 0x1000162E525EE0547) >> 64;
}
if (x & 0x1000000000000 > 0) {
result = (result * 0x10000B17255775C04) >> 64;
}
}
if (x & 0xFF0000000000 > 0) {
if (x & 0x800000000000 > 0) {
result = (result * 0x1000058B91B5BC9AE) >> 64;
}
if (x & 0x400000000000 > 0) {
result = (result * 0x100002C5C89D5EC6D) >> 64;
}
if (x & 0x200000000000 > 0) {
result = (result * 0x10000162E43F4F831) >> 64;
}
if (x & 0x100000000000 > 0) {
result = (result * 0x100000B1721BCFC9A) >> 64;
}
if (x & 0x80000000000 > 0) {
result = (result * 0x10000058B90CF1E6E) >> 64;
}
if (x & 0x40000000000 > 0) {
result = (result * 0x1000002C5C863B73F) >> 64;
}
if (x & 0x20000000000 > 0) {
result = (result * 0x100000162E430E5A2) >> 64;
}
if (x & 0x10000000000 > 0) {
result = (result * 0x1000000B172183551) >> 64;
}
}
if (x & 0xFF00000000 > 0) {
if (x & 0x8000000000 > 0) {
result = (result * 0x100000058B90C0B49) >> 64;
}
if (x & 0x4000000000 > 0) {
result = (result * 0x10000002C5C8601CC) >> 64;
}
if (x & 0x2000000000 > 0) {
result = (result * 0x1000000162E42FFF0) >> 64;
}
if (x & 0x1000000000 > 0) {
result = (result * 0x10000000B17217FBB) >> 64;
}
if (x & 0x800000000 > 0) {
result = (result * 0x1000000058B90BFCE) >> 64;
}
if (x & 0x400000000 > 0) {
result = (result * 0x100000002C5C85FE3) >> 64;
}
if (x & 0x200000000 > 0) {
result = (result * 0x10000000162E42FF1) >> 64;
}
if (x & 0x100000000 > 0) {
result = (result * 0x100000000B17217F8) >> 64;
}
}
if (x & 0xFF00000000 > 0) {
if (x & 0x80000000 > 0) {
result = (result * 0x10000000058B90BFC) >> 64;
}
if (x & 0x40000000 > 0) {
result = (result * 0x1000000002C5C85FE) >> 64;
}
if (x & 0x20000000 > 0) {
result = (result * 0x100000000162E42FF) >> 64;
}
if (x & 0x10000000 > 0) {
result = (result * 0x1000000000B17217F) >> 64;
}
if (x & 0x8000000 > 0) {
result = (result * 0x100000000058B90C0) >> 64;
}
if (x & 0x4000000 > 0) {
result = (result * 0x10000000002C5C860) >> 64;
}
if (x & 0x2000000 > 0) {
result = (result * 0x1000000000162E430) >> 64;
}
if (x & 0x1000000 > 0) {
result = (result * 0x10000000000B17218) >> 64;
}
}
if (x & 0xFF0000 > 0) {
if (x & 0x800000 > 0) {
result = (result * 0x1000000000058B90C) >> 64;
}
if (x & 0x400000 > 0) {
result = (result * 0x100000000002C5C86) >> 64;
}
if (x & 0x200000 > 0) {
result = (result * 0x10000000000162E43) >> 64;
}
if (x & 0x100000 > 0) {
result = (result * 0x100000000000B1721) >> 64;
}
if (x & 0x80000 > 0) {
result = (result * 0x10000000000058B91) >> 64;
}
if (x & 0x40000 > 0) {
result = (result * 0x1000000000002C5C8) >> 64;
}
if (x & 0x20000 > 0) {
result = (result * 0x100000000000162E4) >> 64;
}
if (x & 0x10000 > 0) {
result = (result * 0x1000000000000B172) >> 64;
}
}
if (x & 0xFF00 > 0) {
if (x & 0x8000 > 0) {
result = (result * 0x100000000000058B9) >> 64;
}
if (x & 0x4000 > 0) {
result = (result * 0x10000000000002C5D) >> 64;
}
if (x & 0x2000 > 0) {
result = (result * 0x1000000000000162E) >> 64;
}
if (x & 0x1000 > 0) {
result = (result * 0x10000000000000B17) >> 64;
}
if (x & 0x800 > 0) {
result = (result * 0x1000000000000058C) >> 64;
}
if (x & 0x400 > 0) {
result = (result * 0x100000000000002C6) >> 64;
}
if (x & 0x200 > 0) {
result = (result * 0x10000000000000163) >> 64;
}
if (x & 0x100 > 0) {
result = (result * 0x100000000000000B1) >> 64;
}
}
if (x & 0xFF > 0) {
if (x & 0x80 > 0) {
result = (result * 0x10000000000000059) >> 64;
}
if (x & 0x40 > 0) {
result = (result * 0x1000000000000002C) >> 64;
}
if (x & 0x20 > 0) {
result = (result * 0x10000000000000016) >> 64;
}
if (x & 0x10 > 0) {
result = (result * 0x1000000000000000B) >> 64;
}
if (x & 0x8 > 0) {
result = (result * 0x10000000000000006) >> 64;
}
if (x & 0x4 > 0) {
result = (result * 0x10000000000000003) >> 64;
}
if (x & 0x2 > 0) {
result = (result * 0x10000000000000001) >> 64;
}
if (x & 0x1 > 0) {
result = (result * 0x10000000000000001) >> 64;
}
}
// We're doing two things at the same time:
//
// 1. Multiply the result by 2^n + 1, where "2^n" is the integer part and the one is added to account for
// the fact that we initially set the result to 0.5. This is accomplished by subtracting from 191
// rather than 192.
// 2. Convert the result to the unsigned 60.18-decimal fixed-point format.
//
// This works because 2^(191-ip) = 2^ip / 2^191, where "ip" is the integer part "2^n".
result *= UNIT;
result >>= (191 - (x >> 64));
}
}
/// @notice Calculates the square root of x, rounding down if x is not a perfect square.
/// @dev Uses the Babylonian method https://en.wikipedia.org/wiki/Methods_of_computing_square_roots#Babylonian_method.
/// Credits to OpenZeppelin for the explanations in code comments below.
///
/// Caveats:
/// - This function does not work with fixed-point numbers.
///
/// @param x The uint256 number for which to calculate the square root.
/// @return result The result as an uint256.
function prbSqrt(uint256 x) pure returns (uint256 result) {
if (x == 0) {
return 0;
}
// For our first guess, we get the biggest power of 2 which is smaller than the square root of x.
//
// We know that the "msb" (most significant bit) of x is a power of 2 such that we have:
//
// $$
// msb(x) <= x <= 2*msb(x)$
// $$
//
// We write $msb(x)$ as $2^k$ and we get:
//
// $$
// k = log_2(x)
// $$
//
// Thus we can write the initial inequality as:
//
// $$
// 2^{log_2(x)} <= x <= 2*2^{log_2(x)+1} \\
// sqrt(2^k) <= sqrt(x) < sqrt(2^{k+1}) \\
// 2^{k/2} <= sqrt(x) < 2^{(k+1)/2} <= 2^{(k/2)+1}
// $$
//
// Consequently, $2^{log_2(x) /2}` is a good first approximation of sqrt(x) with at least one correct bit.
uint256 xAux = uint256(x);
result = 1;
if (xAux >= 2 ** 128) {
xAux >>= 128;
result <<= 64;
}
if (xAux >= 2 ** 64) {
xAux >>= 64;
result <<= 32;
}
if (xAux >= 2 ** 32) {
xAux >>= 32;
result <<= 16;
}
if (xAux >= 2 ** 16) {
xAux >>= 16;
result <<= 8;
}
if (xAux >= 2 ** 8) {
xAux >>= 8;
result <<= 4;
}
if (xAux >= 2 ** 4) {
xAux >>= 4;
result <<= 2;
}
if (xAux >= 2 ** 2) {
result <<= 1;
}
// At this point, `result` is an estimation with at least one bit of precision. We know the true value has at
// most 128 bits, 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 + x / result) >> 1;
result = (result + x / result) >> 1;
result = (result + x / result) >> 1;
result = (result + x / result) >> 1;
result = (result + x / result) >> 1;
result = (result + x / result) >> 1;
result = (result + x / result) >> 1;
// Round down the result in case x is not a perfect square.
uint256 roundedDownResult = x / result;
if (result >= roundedDownResult) {
result = roundedDownResult;
}
}
}