Math.sol 28 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685
  1. // SPDX-License-Identifier: MIT
  2. // OpenZeppelin Contracts (last updated v5.0.0) (utils/math/Math.sol)
  3. pragma solidity ^0.8.20;
  4. import {Panic} from "../Panic.sol";
  5. import {SafeCast} from "./SafeCast.sol";
  6. /**
  7. * @dev Standard math utilities missing in the Solidity language.
  8. */
  9. library Math {
  10. enum Rounding {
  11. Floor, // Toward negative infinity
  12. Ceil, // Toward positive infinity
  13. Trunc, // Toward zero
  14. Expand // Away from zero
  15. }
  16. /**
  17. * @dev Returns the addition of two unsigned integers, with an success flag (no overflow).
  18. */
  19. function tryAdd(uint256 a, uint256 b) internal pure returns (bool success, uint256 result) {
  20. unchecked {
  21. uint256 c = a + b;
  22. if (c < a) return (false, 0);
  23. return (true, c);
  24. }
  25. }
  26. /**
  27. * @dev Returns the subtraction of two unsigned integers, with an success flag (no overflow).
  28. */
  29. function trySub(uint256 a, uint256 b) internal pure returns (bool success, uint256 result) {
  30. unchecked {
  31. if (b > a) return (false, 0);
  32. return (true, a - b);
  33. }
  34. }
  35. /**
  36. * @dev Returns the multiplication of two unsigned integers, with an success flag (no overflow).
  37. */
  38. function tryMul(uint256 a, uint256 b) internal pure returns (bool success, uint256 result) {
  39. unchecked {
  40. // Gas optimization: this is cheaper than requiring 'a' not being zero, but the
  41. // benefit is lost if 'b' is also tested.
  42. // See: https://github.com/OpenZeppelin/openzeppelin-contracts/pull/522
  43. if (a == 0) return (true, 0);
  44. uint256 c = a * b;
  45. if (c / a != b) return (false, 0);
  46. return (true, c);
  47. }
  48. }
  49. /**
  50. * @dev Returns the division of two unsigned integers, with a success flag (no division by zero).
  51. */
  52. function tryDiv(uint256 a, uint256 b) internal pure returns (bool success, uint256 result) {
  53. unchecked {
  54. if (b == 0) return (false, 0);
  55. return (true, a / b);
  56. }
  57. }
  58. /**
  59. * @dev Returns the remainder of dividing two unsigned integers, with a success flag (no division by zero).
  60. */
  61. function tryMod(uint256 a, uint256 b) internal pure returns (bool success, uint256 result) {
  62. unchecked {
  63. if (b == 0) return (false, 0);
  64. return (true, a % b);
  65. }
  66. }
  67. /**
  68. * @dev Branchless ternary evaluation for `a ? b : c`. Gas costs are constant.
  69. *
  70. * IMPORTANT: This function may reduce bytecode size and consume less gas when used standalone.
  71. * However, the compiler may optimize Solidity ternary operations (i.e. `a ? b : c`) to only compute
  72. * one branch when needed, making this function more expensive.
  73. */
  74. function ternary(bool condition, uint256 a, uint256 b) internal pure returns (uint256) {
  75. unchecked {
  76. // branchless ternary works because:
  77. // b ^ (a ^ b) == a
  78. // b ^ 0 == b
  79. return b ^ ((a ^ b) * SafeCast.toUint(condition));
  80. }
  81. }
  82. /**
  83. * @dev Returns the largest of two numbers.
  84. */
  85. function max(uint256 a, uint256 b) internal pure returns (uint256) {
  86. return ternary(a > b, a, b);
  87. }
  88. /**
  89. * @dev Returns the smallest of two numbers.
  90. */
  91. function min(uint256 a, uint256 b) internal pure returns (uint256) {
  92. return ternary(a < b, a, b);
  93. }
  94. /**
  95. * @dev Returns the average of two numbers. The result is rounded towards
  96. * zero.
  97. */
  98. function average(uint256 a, uint256 b) internal pure returns (uint256) {
  99. // (a + b) / 2 can overflow.
  100. return (a & b) + (a ^ b) / 2;
  101. }
  102. /**
  103. * @dev Returns the ceiling of the division of two numbers.
  104. *
  105. * This differs from standard division with `/` in that it rounds towards infinity instead
  106. * of rounding towards zero.
  107. */
  108. function ceilDiv(uint256 a, uint256 b) internal pure returns (uint256) {
  109. if (b == 0) {
  110. // Guarantee the same behavior as in a regular Solidity division.
  111. Panic.panic(Panic.DIVISION_BY_ZERO);
  112. }
  113. // The following calculation ensures accurate ceiling division without overflow.
  114. // Since a is non-zero, (a - 1) / b will not overflow.
  115. // The largest possible result occurs when (a - 1) / b is type(uint256).max,
  116. // but the largest value we can obtain is type(uint256).max - 1, which happens
  117. // when a = type(uint256).max and b = 1.
  118. unchecked {
  119. return SafeCast.toUint(a > 0) * ((a - 1) / b + 1);
  120. }
  121. }
  122. /**
  123. * @dev Calculates floor(x * y / denominator) with full precision. Throws if result overflows a uint256 or
  124. * denominator == 0.
  125. *
  126. * Original credit to Remco Bloemen under MIT license (https://xn--2-umb.com/21/muldiv) with further edits by
  127. * Uniswap Labs also under MIT license.
  128. */
  129. function mulDiv(uint256 x, uint256 y, uint256 denominator) internal pure returns (uint256 result) {
  130. unchecked {
  131. // 512-bit multiply [prod1 prod0] = x * y. Compute the product mod 2²⁵⁶ and mod 2²⁵⁶ - 1, then use
  132. // the Chinese Remainder Theorem to reconstruct the 512 bit result. The result is stored in two 256
  133. // variables such that product = prod1 * 2²⁵⁶ + prod0.
  134. uint256 prod0 = x * y; // Least significant 256 bits of the product
  135. uint256 prod1; // Most significant 256 bits of the product
  136. assembly {
  137. let mm := mulmod(x, y, not(0))
  138. prod1 := sub(sub(mm, prod0), lt(mm, prod0))
  139. }
  140. // Handle non-overflow cases, 256 by 256 division.
  141. if (prod1 == 0) {
  142. // Solidity will revert if denominator == 0, unlike the div opcode on its own.
  143. // The surrounding unchecked block does not change this fact.
  144. // See https://docs.soliditylang.org/en/latest/control-structures.html#checked-or-unchecked-arithmetic.
  145. return prod0 / denominator;
  146. }
  147. // Make sure the result is less than 2²⁵⁶. Also prevents denominator == 0.
  148. if (denominator <= prod1) {
  149. Panic.panic(ternary(denominator == 0, Panic.DIVISION_BY_ZERO, Panic.UNDER_OVERFLOW));
  150. }
  151. ///////////////////////////////////////////////
  152. // 512 by 256 division.
  153. ///////////////////////////////////////////////
  154. // Make division exact by subtracting the remainder from [prod1 prod0].
  155. uint256 remainder;
  156. assembly {
  157. // Compute remainder using mulmod.
  158. remainder := mulmod(x, y, denominator)
  159. // Subtract 256 bit number from 512 bit number.
  160. prod1 := sub(prod1, gt(remainder, prod0))
  161. prod0 := sub(prod0, remainder)
  162. }
  163. // Factor powers of two out of denominator and compute largest power of two divisor of denominator.
  164. // Always >= 1. See https://cs.stackexchange.com/q/138556/92363.
  165. uint256 twos = denominator & (0 - denominator);
  166. assembly {
  167. // Divide denominator by twos.
  168. denominator := div(denominator, twos)
  169. // Divide [prod1 prod0] by twos.
  170. prod0 := div(prod0, twos)
  171. // Flip twos such that it is 2²⁵⁶ / twos. If twos is zero, then it becomes one.
  172. twos := add(div(sub(0, twos), twos), 1)
  173. }
  174. // Shift in bits from prod1 into prod0.
  175. prod0 |= prod1 * twos;
  176. // Invert denominator mod 2²⁵⁶. Now that denominator is an odd number, it has an inverse modulo 2²⁵⁶ such
  177. // that denominator * inv ≡ 1 mod 2²⁵⁶. Compute the inverse by starting with a seed that is correct for
  178. // four bits. That is, denominator * inv ≡ 1 mod 2⁴.
  179. uint256 inverse = (3 * denominator) ^ 2;
  180. // Use the Newton-Raphson iteration to improve the precision. Thanks to Hensel's lifting lemma, this also
  181. // works in modular arithmetic, doubling the correct bits in each step.
  182. inverse *= 2 - denominator * inverse; // inverse mod 2⁸
  183. inverse *= 2 - denominator * inverse; // inverse mod 2¹⁶
  184. inverse *= 2 - denominator * inverse; // inverse mod 2³²
  185. inverse *= 2 - denominator * inverse; // inverse mod 2⁶⁴
  186. inverse *= 2 - denominator * inverse; // inverse mod 2¹²⁸
  187. inverse *= 2 - denominator * inverse; // inverse mod 2²⁵⁶
  188. // Because the division is now exact we can divide by multiplying with the modular inverse of denominator.
  189. // This will give us the correct result modulo 2²⁵⁶. Since the preconditions guarantee that the outcome is
  190. // less than 2²⁵⁶, this is the final result. We don't need to compute the high bits of the result and prod1
  191. // is no longer required.
  192. result = prod0 * inverse;
  193. return result;
  194. }
  195. }
  196. /**
  197. * @dev Calculates x * y / denominator with full precision, following the selected rounding direction.
  198. */
  199. function mulDiv(uint256 x, uint256 y, uint256 denominator, Rounding rounding) internal pure returns (uint256) {
  200. return mulDiv(x, y, denominator) + SafeCast.toUint(unsignedRoundsUp(rounding) && mulmod(x, y, denominator) > 0);
  201. }
  202. /**
  203. * @dev Calculate the modular multiplicative inverse of a number in Z/nZ.
  204. *
  205. * If n is a prime, then Z/nZ is a field. In that case all elements are inversible, except 0.
  206. * If n is not a prime, then Z/nZ is not a field, and some elements might not be inversible.
  207. *
  208. * If the input value is not inversible, 0 is returned.
  209. *
  210. * NOTE: If you know for sure that n is (big) a prime, it may be cheaper to use Fermat's little theorem and get the
  211. * inverse using `Math.modExp(a, n - 2, n)`. See {invModPrime}.
  212. */
  213. function invMod(uint256 a, uint256 n) internal pure returns (uint256) {
  214. unchecked {
  215. if (n == 0) return 0;
  216. // The inverse modulo is calculated using the Extended Euclidean Algorithm (iterative version)
  217. // Used to compute integers x and y such that: ax + ny = gcd(a, n).
  218. // When the gcd is 1, then the inverse of a modulo n exists and it's x.
  219. // ax + ny = 1
  220. // ax = 1 + (-y)n
  221. // ax ≡ 1 (mod n) # x is the inverse of a modulo n
  222. // If the remainder is 0 the gcd is n right away.
  223. uint256 remainder = a % n;
  224. uint256 gcd = n;
  225. // Therefore the initial coefficients are:
  226. // ax + ny = gcd(a, n) = n
  227. // 0a + 1n = n
  228. int256 x = 0;
  229. int256 y = 1;
  230. while (remainder != 0) {
  231. uint256 quotient = gcd / remainder;
  232. (gcd, remainder) = (
  233. // The old remainder is the next gcd to try.
  234. remainder,
  235. // Compute the next remainder.
  236. // Can't overflow given that (a % gcd) * (gcd // (a % gcd)) <= gcd
  237. // where gcd is at most n (capped to type(uint256).max)
  238. gcd - remainder * quotient
  239. );
  240. (x, y) = (
  241. // Increment the coefficient of a.
  242. y,
  243. // Decrement the coefficient of n.
  244. // Can overflow, but the result is casted to uint256 so that the
  245. // next value of y is "wrapped around" to a value between 0 and n - 1.
  246. x - y * int256(quotient)
  247. );
  248. }
  249. if (gcd != 1) return 0; // No inverse exists.
  250. return ternary(x < 0, n - uint256(-x), uint256(x)); // Wrap the result if it's negative.
  251. }
  252. }
  253. /**
  254. * @dev Variant of {invMod}. More efficient, but only works if `p` is known to be a prime greater than `2`.
  255. *
  256. * From https://en.wikipedia.org/wiki/Fermat%27s_little_theorem[Fermat's little theorem], we know that if p is
  257. * prime, then `a**(p-1) ≡ 1 mod p`. As a consequence, we have `a * a**(p-2) ≡ 1 mod p`, which means that
  258. * `a**(p-2)` is the modular multiplicative inverse of a in Fp.
  259. *
  260. * NOTE: this function does NOT check that `p` is a prime greater than `2`.
  261. */
  262. function invModPrime(uint256 a, uint256 p) internal view returns (uint256) {
  263. unchecked {
  264. return Math.modExp(a, p - 2, p);
  265. }
  266. }
  267. /**
  268. * @dev Returns the modular exponentiation of the specified base, exponent and modulus (b ** e % m)
  269. *
  270. * Requirements:
  271. * - modulus can't be zero
  272. * - underlying staticcall to precompile must succeed
  273. *
  274. * IMPORTANT: The result is only valid if the underlying call succeeds. When using this function, make
  275. * sure the chain you're using it on supports the precompiled contract for modular exponentiation
  276. * at address 0x05 as specified in https://eips.ethereum.org/EIPS/eip-198[EIP-198]. Otherwise,
  277. * the underlying function will succeed given the lack of a revert, but the result may be incorrectly
  278. * interpreted as 0.
  279. */
  280. function modExp(uint256 b, uint256 e, uint256 m) internal view returns (uint256) {
  281. (bool success, uint256 result) = tryModExp(b, e, m);
  282. if (!success) {
  283. Panic.panic(Panic.DIVISION_BY_ZERO);
  284. }
  285. return result;
  286. }
  287. /**
  288. * @dev Returns the modular exponentiation of the specified base, exponent and modulus (b ** e % m).
  289. * It includes a success flag indicating if the operation succeeded. Operation will be marked as failed if trying
  290. * to operate modulo 0 or if the underlying precompile reverted.
  291. *
  292. * IMPORTANT: The result is only valid if the success flag is true. When using this function, make sure the chain
  293. * you're using it on supports the precompiled contract for modular exponentiation at address 0x05 as specified in
  294. * https://eips.ethereum.org/EIPS/eip-198[EIP-198]. Otherwise, the underlying function will succeed given the lack
  295. * of a revert, but the result may be incorrectly interpreted as 0.
  296. */
  297. function tryModExp(uint256 b, uint256 e, uint256 m) internal view returns (bool success, uint256 result) {
  298. if (m == 0) return (false, 0);
  299. assembly ("memory-safe") {
  300. let ptr := mload(0x40)
  301. // | Offset | Content | Content (Hex) |
  302. // |-----------|------------|--------------------------------------------------------------------|
  303. // | 0x00:0x1f | size of b | 0x0000000000000000000000000000000000000000000000000000000000000020 |
  304. // | 0x20:0x3f | size of e | 0x0000000000000000000000000000000000000000000000000000000000000020 |
  305. // | 0x40:0x5f | size of m | 0x0000000000000000000000000000000000000000000000000000000000000020 |
  306. // | 0x60:0x7f | value of b | 0x<.............................................................b> |
  307. // | 0x80:0x9f | value of e | 0x<.............................................................e> |
  308. // | 0xa0:0xbf | value of m | 0x<.............................................................m> |
  309. mstore(ptr, 0x20)
  310. mstore(add(ptr, 0x20), 0x20)
  311. mstore(add(ptr, 0x40), 0x20)
  312. mstore(add(ptr, 0x60), b)
  313. mstore(add(ptr, 0x80), e)
  314. mstore(add(ptr, 0xa0), m)
  315. // Given the result < m, it's guaranteed to fit in 32 bytes,
  316. // so we can use the memory scratch space located at offset 0.
  317. success := staticcall(gas(), 0x05, ptr, 0xc0, 0x00, 0x20)
  318. result := mload(0x00)
  319. }
  320. }
  321. /**
  322. * @dev Variant of {modExp} that supports inputs of arbitrary length.
  323. */
  324. function modExp(bytes memory b, bytes memory e, bytes memory m) internal view returns (bytes memory) {
  325. (bool success, bytes memory result) = tryModExp(b, e, m);
  326. if (!success) {
  327. Panic.panic(Panic.DIVISION_BY_ZERO);
  328. }
  329. return result;
  330. }
  331. /**
  332. * @dev Variant of {tryModExp} that supports inputs of arbitrary length.
  333. */
  334. function tryModExp(
  335. bytes memory b,
  336. bytes memory e,
  337. bytes memory m
  338. ) internal view returns (bool success, bytes memory result) {
  339. if (_zeroBytes(m)) return (false, new bytes(0));
  340. uint256 mLen = m.length;
  341. // Encode call args in result and move the free memory pointer
  342. result = abi.encodePacked(b.length, e.length, mLen, b, e, m);
  343. assembly ("memory-safe") {
  344. let dataPtr := add(result, 0x20)
  345. // Write result on top of args to avoid allocating extra memory.
  346. success := staticcall(gas(), 0x05, dataPtr, mload(result), dataPtr, mLen)
  347. // Overwrite the length.
  348. // result.length > returndatasize() is guaranteed because returndatasize() == m.length
  349. mstore(result, mLen)
  350. // Set the memory pointer after the returned data.
  351. mstore(0x40, add(dataPtr, mLen))
  352. }
  353. }
  354. /**
  355. * @dev Returns whether the provided byte array is zero.
  356. */
  357. function _zeroBytes(bytes memory byteArray) private pure returns (bool) {
  358. for (uint256 i = 0; i < byteArray.length; ++i) {
  359. if (byteArray[i] != 0) {
  360. return false;
  361. }
  362. }
  363. return true;
  364. }
  365. /**
  366. * @dev Returns the square root of a number. If the number is not a perfect square, the value is rounded
  367. * towards zero.
  368. *
  369. * This method is based on Newton's method for computing square roots; the algorithm is restricted to only
  370. * using integer operations.
  371. */
  372. function sqrt(uint256 a) internal pure returns (uint256) {
  373. unchecked {
  374. // Take care of easy edge cases when a == 0 or a == 1
  375. if (a <= 1) {
  376. return a;
  377. }
  378. // In this function, we use Newton's method to get a root of `f(x) := x² - a`. It involves building a
  379. // sequence x_n that converges toward sqrt(a). For each iteration x_n, we also define the error between
  380. // the current value as `ε_n = | x_n - sqrt(a) |`.
  381. //
  382. // For our first estimation, we consider `e` the smallest power of 2 which is bigger than the square root
  383. // of the target. (i.e. `2**(e-1) ≤ sqrt(a) < 2**e`). We know that `e ≤ 128` because `(2¹²⁸)² = 2²⁵⁶` is
  384. // bigger than any uint256.
  385. //
  386. // By noticing that
  387. // `2**(e-1) ≤ sqrt(a) < 2**e → (2**(e-1))² ≤ a < (2**e)² → 2**(2*e-2) ≤ a < 2**(2*e)`
  388. // we can deduce that `e - 1` is `log2(a) / 2`. We can thus compute `x_n = 2**(e-1)` using a method similar
  389. // to the msb function.
  390. uint256 aa = a;
  391. uint256 xn = 1;
  392. if (aa >= (1 << 128)) {
  393. aa >>= 128;
  394. xn <<= 64;
  395. }
  396. if (aa >= (1 << 64)) {
  397. aa >>= 64;
  398. xn <<= 32;
  399. }
  400. if (aa >= (1 << 32)) {
  401. aa >>= 32;
  402. xn <<= 16;
  403. }
  404. if (aa >= (1 << 16)) {
  405. aa >>= 16;
  406. xn <<= 8;
  407. }
  408. if (aa >= (1 << 8)) {
  409. aa >>= 8;
  410. xn <<= 4;
  411. }
  412. if (aa >= (1 << 4)) {
  413. aa >>= 4;
  414. xn <<= 2;
  415. }
  416. if (aa >= (1 << 2)) {
  417. xn <<= 1;
  418. }
  419. // We now have x_n such that `x_n = 2**(e-1) ≤ sqrt(a) < 2**e = 2 * x_n`. This implies ε_n ≤ 2**(e-1).
  420. //
  421. // We can refine our estimation by noticing that the middle of that interval minimizes the error.
  422. // If we move x_n to equal 2**(e-1) + 2**(e-2), then we reduce the error to ε_n ≤ 2**(e-2).
  423. // This is going to be our x_0 (and ε_0)
  424. xn = (3 * xn) >> 1; // ε_0 := | x_0 - sqrt(a) | ≤ 2**(e-2)
  425. // From here, Newton's method give us:
  426. // x_{n+1} = (x_n + a / x_n) / 2
  427. //
  428. // One should note that:
  429. // x_{n+1}² - a = ((x_n + a / x_n) / 2)² - a
  430. // = ((x_n² + a) / (2 * x_n))² - a
  431. // = (x_n⁴ + 2 * a * x_n² + a²) / (4 * x_n²) - a
  432. // = (x_n⁴ + 2 * a * x_n² + a² - 4 * a * x_n²) / (4 * x_n²)
  433. // = (x_n⁴ - 2 * a * x_n² + a²) / (4 * x_n²)
  434. // = (x_n² - a)² / (2 * x_n)²
  435. // = ((x_n² - a) / (2 * x_n))²
  436. // ≥ 0
  437. // Which proves that for all n ≥ 1, sqrt(a) ≤ x_n
  438. //
  439. // This gives us the proof of quadratic convergence of the sequence:
  440. // ε_{n+1} = | x_{n+1} - sqrt(a) |
  441. // = | (x_n + a / x_n) / 2 - sqrt(a) |
  442. // = | (x_n² + a - 2*x_n*sqrt(a)) / (2 * x_n) |
  443. // = | (x_n - sqrt(a))² / (2 * x_n) |
  444. // = | ε_n² / (2 * x_n) |
  445. // = ε_n² / | (2 * x_n) |
  446. //
  447. // For the first iteration, we have a special case where x_0 is known:
  448. // ε_1 = ε_0² / | (2 * x_0) |
  449. // ≤ (2**(e-2))² / (2 * (2**(e-1) + 2**(e-2)))
  450. // ≤ 2**(2*e-4) / (3 * 2**(e-1))
  451. // ≤ 2**(e-3) / 3
  452. // ≤ 2**(e-3-log2(3))
  453. // ≤ 2**(e-4.5)
  454. //
  455. // For the following iterations, we use the fact that, 2**(e-1) ≤ sqrt(a) ≤ x_n:
  456. // ε_{n+1} = ε_n² / | (2 * x_n) |
  457. // ≤ (2**(e-k))² / (2 * 2**(e-1))
  458. // ≤ 2**(2*e-2*k) / 2**e
  459. // ≤ 2**(e-2*k)
  460. xn = (xn + a / xn) >> 1; // ε_1 := | x_1 - sqrt(a) | ≤ 2**(e-4.5) -- special case, see above
  461. xn = (xn + a / xn) >> 1; // ε_2 := | x_2 - sqrt(a) | ≤ 2**(e-9) -- general case with k = 4.5
  462. xn = (xn + a / xn) >> 1; // ε_3 := | x_3 - sqrt(a) | ≤ 2**(e-18) -- general case with k = 9
  463. xn = (xn + a / xn) >> 1; // ε_4 := | x_4 - sqrt(a) | ≤ 2**(e-36) -- general case with k = 18
  464. xn = (xn + a / xn) >> 1; // ε_5 := | x_5 - sqrt(a) | ≤ 2**(e-72) -- general case with k = 36
  465. xn = (xn + a / xn) >> 1; // ε_6 := | x_6 - sqrt(a) | ≤ 2**(e-144) -- general case with k = 72
  466. // Because e ≤ 128 (as discussed during the first estimation phase), we know have reached a precision
  467. // ε_6 ≤ 2**(e-144) < 1. Given we're operating on integers, then we can ensure that xn is now either
  468. // sqrt(a) or sqrt(a) + 1.
  469. return xn - SafeCast.toUint(xn > a / xn);
  470. }
  471. }
  472. /**
  473. * @dev Calculates sqrt(a), following the selected rounding direction.
  474. */
  475. function sqrt(uint256 a, Rounding rounding) internal pure returns (uint256) {
  476. unchecked {
  477. uint256 result = sqrt(a);
  478. return result + SafeCast.toUint(unsignedRoundsUp(rounding) && result * result < a);
  479. }
  480. }
  481. /**
  482. * @dev Return the log in base 2 of a positive value rounded towards zero.
  483. * Returns 0 if given 0.
  484. */
  485. function log2(uint256 value) internal pure returns (uint256) {
  486. uint256 result = 0;
  487. uint256 exp;
  488. unchecked {
  489. exp = 128 * SafeCast.toUint(value > (1 << 128) - 1);
  490. value >>= exp;
  491. result += exp;
  492. exp = 64 * SafeCast.toUint(value > (1 << 64) - 1);
  493. value >>= exp;
  494. result += exp;
  495. exp = 32 * SafeCast.toUint(value > (1 << 32) - 1);
  496. value >>= exp;
  497. result += exp;
  498. exp = 16 * SafeCast.toUint(value > (1 << 16) - 1);
  499. value >>= exp;
  500. result += exp;
  501. exp = 8 * SafeCast.toUint(value > (1 << 8) - 1);
  502. value >>= exp;
  503. result += exp;
  504. exp = 4 * SafeCast.toUint(value > (1 << 4) - 1);
  505. value >>= exp;
  506. result += exp;
  507. exp = 2 * SafeCast.toUint(value > (1 << 2) - 1);
  508. value >>= exp;
  509. result += exp;
  510. result += SafeCast.toUint(value > 1);
  511. }
  512. return result;
  513. }
  514. /**
  515. * @dev Return the log in base 2, following the selected rounding direction, of a positive value.
  516. * Returns 0 if given 0.
  517. */
  518. function log2(uint256 value, Rounding rounding) internal pure returns (uint256) {
  519. unchecked {
  520. uint256 result = log2(value);
  521. return result + SafeCast.toUint(unsignedRoundsUp(rounding) && 1 << result < value);
  522. }
  523. }
  524. /**
  525. * @dev Return the log in base 10 of a positive value rounded towards zero.
  526. * Returns 0 if given 0.
  527. */
  528. function log10(uint256 value) internal pure returns (uint256) {
  529. uint256 result = 0;
  530. unchecked {
  531. if (value >= 10 ** 64) {
  532. value /= 10 ** 64;
  533. result += 64;
  534. }
  535. if (value >= 10 ** 32) {
  536. value /= 10 ** 32;
  537. result += 32;
  538. }
  539. if (value >= 10 ** 16) {
  540. value /= 10 ** 16;
  541. result += 16;
  542. }
  543. if (value >= 10 ** 8) {
  544. value /= 10 ** 8;
  545. result += 8;
  546. }
  547. if (value >= 10 ** 4) {
  548. value /= 10 ** 4;
  549. result += 4;
  550. }
  551. if (value >= 10 ** 2) {
  552. value /= 10 ** 2;
  553. result += 2;
  554. }
  555. if (value >= 10 ** 1) {
  556. result += 1;
  557. }
  558. }
  559. return result;
  560. }
  561. /**
  562. * @dev Return the log in base 10, following the selected rounding direction, of a positive value.
  563. * Returns 0 if given 0.
  564. */
  565. function log10(uint256 value, Rounding rounding) internal pure returns (uint256) {
  566. unchecked {
  567. uint256 result = log10(value);
  568. return result + SafeCast.toUint(unsignedRoundsUp(rounding) && 10 ** result < value);
  569. }
  570. }
  571. /**
  572. * @dev Return the log in base 256 of a positive value rounded towards zero.
  573. * Returns 0 if given 0.
  574. *
  575. * Adding one to the result gives the number of pairs of hex symbols needed to represent `value` as a hex string.
  576. */
  577. function log256(uint256 value) internal pure returns (uint256) {
  578. uint256 result = 0;
  579. uint256 isGt;
  580. unchecked {
  581. isGt = SafeCast.toUint(value > (1 << 128) - 1);
  582. value >>= isGt * 128;
  583. result += isGt * 16;
  584. isGt = SafeCast.toUint(value > (1 << 64) - 1);
  585. value >>= isGt * 64;
  586. result += isGt * 8;
  587. isGt = SafeCast.toUint(value > (1 << 32) - 1);
  588. value >>= isGt * 32;
  589. result += isGt * 4;
  590. isGt = SafeCast.toUint(value > (1 << 16) - 1);
  591. value >>= isGt * 16;
  592. result += isGt * 2;
  593. result += SafeCast.toUint(value > (1 << 8) - 1);
  594. }
  595. return result;
  596. }
  597. /**
  598. * @dev Return the log in base 256, following the selected rounding direction, of a positive value.
  599. * Returns 0 if given 0.
  600. */
  601. function log256(uint256 value, Rounding rounding) internal pure returns (uint256) {
  602. unchecked {
  603. uint256 result = log256(value);
  604. return result + SafeCast.toUint(unsignedRoundsUp(rounding) && 1 << (result << 3) < value);
  605. }
  606. }
  607. /**
  608. * @dev Returns whether a provided rounding mode is considered rounding up for unsigned integers.
  609. */
  610. function unsignedRoundsUp(Rounding rounding) internal pure returns (bool) {
  611. return uint8(rounding) % 2 == 1;
  612. }
  613. }