Math.sol 27 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672
  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. // use 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, expect 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 Ferma's little theorem and get the
  211. * inverse using `Math.modExp(a, n - 2, n)`.
  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 Returns the modular exponentiation of the specified base, exponent and modulus (b ** e % m)
  255. *
  256. * Requirements:
  257. * - modulus can't be zero
  258. * - underlying staticcall to precompile must succeed
  259. *
  260. * IMPORTANT: The result is only valid if the underlying call succeeds. When using this function, make
  261. * sure the chain you're using it on supports the precompiled contract for modular exponentiation
  262. * at address 0x05 as specified in https://eips.ethereum.org/EIPS/eip-198[EIP-198]. Otherwise,
  263. * the underlying function will succeed given the lack of a revert, but the result may be incorrectly
  264. * interpreted as 0.
  265. */
  266. function modExp(uint256 b, uint256 e, uint256 m) internal view returns (uint256) {
  267. (bool success, uint256 result) = tryModExp(b, e, m);
  268. if (!success) {
  269. Panic.panic(Panic.DIVISION_BY_ZERO);
  270. }
  271. return result;
  272. }
  273. /**
  274. * @dev Returns the modular exponentiation of the specified base, exponent and modulus (b ** e % m).
  275. * It includes a success flag indicating if the operation succeeded. Operation will be marked has failed if trying
  276. * to operate modulo 0 or if the underlying precompile reverted.
  277. *
  278. * IMPORTANT: The result is only valid if the success flag is true. When using this function, make sure the chain
  279. * you're using it on supports the precompiled contract for modular exponentiation at address 0x05 as specified in
  280. * https://eips.ethereum.org/EIPS/eip-198[EIP-198]. Otherwise, the underlying function will succeed given the lack
  281. * of a revert, but the result may be incorrectly interpreted as 0.
  282. */
  283. function tryModExp(uint256 b, uint256 e, uint256 m) internal view returns (bool success, uint256 result) {
  284. if (m == 0) return (false, 0);
  285. /// @solidity memory-safe-assembly
  286. assembly {
  287. let ptr := mload(0x40)
  288. // | Offset | Content | Content (Hex) |
  289. // |-----------|------------|--------------------------------------------------------------------|
  290. // | 0x00:0x1f | size of b | 0x0000000000000000000000000000000000000000000000000000000000000020 |
  291. // | 0x20:0x3f | size of e | 0x0000000000000000000000000000000000000000000000000000000000000020 |
  292. // | 0x40:0x5f | size of m | 0x0000000000000000000000000000000000000000000000000000000000000020 |
  293. // | 0x60:0x7f | value of b | 0x<.............................................................b> |
  294. // | 0x80:0x9f | value of e | 0x<.............................................................e> |
  295. // | 0xa0:0xbf | value of m | 0x<.............................................................m> |
  296. mstore(ptr, 0x20)
  297. mstore(add(ptr, 0x20), 0x20)
  298. mstore(add(ptr, 0x40), 0x20)
  299. mstore(add(ptr, 0x60), b)
  300. mstore(add(ptr, 0x80), e)
  301. mstore(add(ptr, 0xa0), m)
  302. // Given the result < m, it's guaranteed to fit in 32 bytes,
  303. // so we can use the memory scratch space located at offset 0.
  304. success := staticcall(gas(), 0x05, ptr, 0xc0, 0x00, 0x20)
  305. result := mload(0x00)
  306. }
  307. }
  308. /**
  309. * @dev Variant of {modExp} that supports inputs of arbitrary length.
  310. */
  311. function modExp(bytes memory b, bytes memory e, bytes memory m) internal view returns (bytes memory) {
  312. (bool success, bytes memory result) = tryModExp(b, e, m);
  313. if (!success) {
  314. Panic.panic(Panic.DIVISION_BY_ZERO);
  315. }
  316. return result;
  317. }
  318. /**
  319. * @dev Variant of {tryModExp} that supports inputs of arbitrary length.
  320. */
  321. function tryModExp(
  322. bytes memory b,
  323. bytes memory e,
  324. bytes memory m
  325. ) internal view returns (bool success, bytes memory result) {
  326. if (_zeroBytes(m)) return (false, new bytes(0));
  327. uint256 mLen = m.length;
  328. // Encode call args in result and move the free memory pointer
  329. result = abi.encodePacked(b.length, e.length, mLen, b, e, m);
  330. /// @solidity memory-safe-assembly
  331. assembly {
  332. let dataPtr := add(result, 0x20)
  333. // Write result on top of args to avoid allocating extra memory.
  334. success := staticcall(gas(), 0x05, dataPtr, mload(result), dataPtr, mLen)
  335. // Overwrite the length.
  336. // result.length > returndatasize() is guaranteed because returndatasize() == m.length
  337. mstore(result, mLen)
  338. // Set the memory pointer after the returned data.
  339. mstore(0x40, add(dataPtr, mLen))
  340. }
  341. }
  342. /**
  343. * @dev Returns whether the provided byte array is zero.
  344. */
  345. function _zeroBytes(bytes memory byteArray) private pure returns (bool) {
  346. for (uint256 i = 0; i < byteArray.length; ++i) {
  347. if (byteArray[i] != 0) {
  348. return false;
  349. }
  350. }
  351. return true;
  352. }
  353. /**
  354. * @dev Returns the square root of a number. If the number is not a perfect square, the value is rounded
  355. * towards zero.
  356. *
  357. * This method is based on Newton's method for computing square roots; the algorithm is restricted to only
  358. * using integer operations.
  359. */
  360. function sqrt(uint256 a) internal pure returns (uint256) {
  361. unchecked {
  362. // Take care of easy edge cases when a == 0 or a == 1
  363. if (a <= 1) {
  364. return a;
  365. }
  366. // In this function, we use Newton's method to get a root of `f(x) := x² - a`. It involves building a
  367. // sequence x_n that converges toward sqrt(a). For each iteration x_n, we also define the error between
  368. // the current value as `ε_n = | x_n - sqrt(a) |`.
  369. //
  370. // For our first estimation, we consider `e` the smallest power of 2 which is bigger than the square root
  371. // of the target. (i.e. `2**(e-1) ≤ sqrt(a) < 2**e`). We know that `e ≤ 128` because `(2¹²⁸)² = 2²⁵⁶` is
  372. // bigger than any uint256.
  373. //
  374. // By noticing that
  375. // `2**(e-1) ≤ sqrt(a) < 2**e → (2**(e-1))² ≤ a < (2**e)² → 2**(2*e-2) ≤ a < 2**(2*e)`
  376. // we can deduce that `e - 1` is `log2(a) / 2`. We can thus compute `x_n = 2**(e-1)` using a method similar
  377. // to the msb function.
  378. uint256 aa = a;
  379. uint256 xn = 1;
  380. if (aa >= (1 << 128)) {
  381. aa >>= 128;
  382. xn <<= 64;
  383. }
  384. if (aa >= (1 << 64)) {
  385. aa >>= 64;
  386. xn <<= 32;
  387. }
  388. if (aa >= (1 << 32)) {
  389. aa >>= 32;
  390. xn <<= 16;
  391. }
  392. if (aa >= (1 << 16)) {
  393. aa >>= 16;
  394. xn <<= 8;
  395. }
  396. if (aa >= (1 << 8)) {
  397. aa >>= 8;
  398. xn <<= 4;
  399. }
  400. if (aa >= (1 << 4)) {
  401. aa >>= 4;
  402. xn <<= 2;
  403. }
  404. if (aa >= (1 << 2)) {
  405. xn <<= 1;
  406. }
  407. // 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).
  408. //
  409. // We can refine our estimation by noticing that the middle of that interval minimizes the error.
  410. // If we move x_n to equal 2**(e-1) + 2**(e-2), then we reduce the error to ε_n ≤ 2**(e-2).
  411. // This is going to be our x_0 (and ε_0)
  412. xn = (3 * xn) >> 1; // ε_0 := | x_0 - sqrt(a) | ≤ 2**(e-2)
  413. // From here, Newton's method give us:
  414. // x_{n+1} = (x_n + a / x_n) / 2
  415. //
  416. // One should note that:
  417. // x_{n+1}² - a = ((x_n + a / x_n) / 2)² - a
  418. // = ((x_n² + a) / (2 * x_n))² - a
  419. // = (x_n⁴ + 2 * a * x_n² + a²) / (4 * x_n²) - a
  420. // = (x_n⁴ + 2 * a * x_n² + a² - 4 * a * x_n²) / (4 * x_n²)
  421. // = (x_n⁴ - 2 * a * x_n² + a²) / (4 * x_n²)
  422. // = (x_n² - a)² / (2 * x_n)²
  423. // = ((x_n² - a) / (2 * x_n))²
  424. // ≥ 0
  425. // Which proves that for all n ≥ 1, sqrt(a) ≤ x_n
  426. //
  427. // This gives us the proof of quadratic convergence of the sequence:
  428. // ε_{n+1} = | x_{n+1} - sqrt(a) |
  429. // = | (x_n + a / x_n) / 2 - sqrt(a) |
  430. // = | (x_n² + a - 2*x_n*sqrt(a)) / (2 * x_n) |
  431. // = | (x_n - sqrt(a))² / (2 * x_n) |
  432. // = | ε_n² / (2 * x_n) |
  433. // = ε_n² / | (2 * x_n) |
  434. //
  435. // For the first iteration, we have a special case where x_0 is known:
  436. // ε_1 = ε_0² / | (2 * x_0) |
  437. // ≤ (2**(e-2))² / (2 * (2**(e-1) + 2**(e-2)))
  438. // ≤ 2**(2*e-4) / (3 * 2**(e-1))
  439. // ≤ 2**(e-3) / 3
  440. // ≤ 2**(e-3-log2(3))
  441. // ≤ 2**(e-4.5)
  442. //
  443. // For the following iterations, we use the fact that, 2**(e-1) ≤ sqrt(a) ≤ x_n:
  444. // ε_{n+1} = ε_n² / | (2 * x_n) |
  445. // ≤ (2**(e-k))² / (2 * 2**(e-1))
  446. // ≤ 2**(2*e-2*k) / 2**e
  447. // ≤ 2**(e-2*k)
  448. xn = (xn + a / xn) >> 1; // ε_1 := | x_1 - sqrt(a) | ≤ 2**(e-4.5) -- special case, see above
  449. xn = (xn + a / xn) >> 1; // ε_2 := | x_2 - sqrt(a) | ≤ 2**(e-9) -- general case with k = 4.5
  450. xn = (xn + a / xn) >> 1; // ε_3 := | x_3 - sqrt(a) | ≤ 2**(e-18) -- general case with k = 9
  451. xn = (xn + a / xn) >> 1; // ε_4 := | x_4 - sqrt(a) | ≤ 2**(e-36) -- general case with k = 18
  452. xn = (xn + a / xn) >> 1; // ε_5 := | x_5 - sqrt(a) | ≤ 2**(e-72) -- general case with k = 36
  453. xn = (xn + a / xn) >> 1; // ε_6 := | x_6 - sqrt(a) | ≤ 2**(e-144) -- general case with k = 72
  454. // Because e ≤ 128 (as discussed during the first estimation phase), we know have reached a precision
  455. // ε_6 ≤ 2**(e-144) < 1. Given we're operating on integers, then we can ensure that xn is now either
  456. // sqrt(a) or sqrt(a) + 1.
  457. return xn - SafeCast.toUint(xn > a / xn);
  458. }
  459. }
  460. /**
  461. * @dev Calculates sqrt(a), following the selected rounding direction.
  462. */
  463. function sqrt(uint256 a, Rounding rounding) internal pure returns (uint256) {
  464. unchecked {
  465. uint256 result = sqrt(a);
  466. return result + SafeCast.toUint(unsignedRoundsUp(rounding) && result * result < a);
  467. }
  468. }
  469. /**
  470. * @dev Return the log in base 2 of a positive value rounded towards zero.
  471. * Returns 0 if given 0.
  472. */
  473. function log2(uint256 value) internal pure returns (uint256) {
  474. uint256 result = 0;
  475. uint256 exp;
  476. unchecked {
  477. exp = 128 * SafeCast.toUint(value > (1 << 128) - 1);
  478. value >>= exp;
  479. result += exp;
  480. exp = 64 * SafeCast.toUint(value > (1 << 64) - 1);
  481. value >>= exp;
  482. result += exp;
  483. exp = 32 * SafeCast.toUint(value > (1 << 32) - 1);
  484. value >>= exp;
  485. result += exp;
  486. exp = 16 * SafeCast.toUint(value > (1 << 16) - 1);
  487. value >>= exp;
  488. result += exp;
  489. exp = 8 * SafeCast.toUint(value > (1 << 8) - 1);
  490. value >>= exp;
  491. result += exp;
  492. exp = 4 * SafeCast.toUint(value > (1 << 4) - 1);
  493. value >>= exp;
  494. result += exp;
  495. exp = 2 * SafeCast.toUint(value > (1 << 2) - 1);
  496. value >>= exp;
  497. result += exp;
  498. result += SafeCast.toUint(value > 1);
  499. }
  500. return result;
  501. }
  502. /**
  503. * @dev Return the log in base 2, following the selected rounding direction, of a positive value.
  504. * Returns 0 if given 0.
  505. */
  506. function log2(uint256 value, Rounding rounding) internal pure returns (uint256) {
  507. unchecked {
  508. uint256 result = log2(value);
  509. return result + SafeCast.toUint(unsignedRoundsUp(rounding) && 1 << result < value);
  510. }
  511. }
  512. /**
  513. * @dev Return the log in base 10 of a positive value rounded towards zero.
  514. * Returns 0 if given 0.
  515. */
  516. function log10(uint256 value) internal pure returns (uint256) {
  517. uint256 result = 0;
  518. unchecked {
  519. if (value >= 10 ** 64) {
  520. value /= 10 ** 64;
  521. result += 64;
  522. }
  523. if (value >= 10 ** 32) {
  524. value /= 10 ** 32;
  525. result += 32;
  526. }
  527. if (value >= 10 ** 16) {
  528. value /= 10 ** 16;
  529. result += 16;
  530. }
  531. if (value >= 10 ** 8) {
  532. value /= 10 ** 8;
  533. result += 8;
  534. }
  535. if (value >= 10 ** 4) {
  536. value /= 10 ** 4;
  537. result += 4;
  538. }
  539. if (value >= 10 ** 2) {
  540. value /= 10 ** 2;
  541. result += 2;
  542. }
  543. if (value >= 10 ** 1) {
  544. result += 1;
  545. }
  546. }
  547. return result;
  548. }
  549. /**
  550. * @dev Return the log in base 10, following the selected rounding direction, of a positive value.
  551. * Returns 0 if given 0.
  552. */
  553. function log10(uint256 value, Rounding rounding) internal pure returns (uint256) {
  554. unchecked {
  555. uint256 result = log10(value);
  556. return result + SafeCast.toUint(unsignedRoundsUp(rounding) && 10 ** result < value);
  557. }
  558. }
  559. /**
  560. * @dev Return the log in base 256 of a positive value rounded towards zero.
  561. * Returns 0 if given 0.
  562. *
  563. * Adding one to the result gives the number of pairs of hex symbols needed to represent `value` as a hex string.
  564. */
  565. function log256(uint256 value) internal pure returns (uint256) {
  566. uint256 result = 0;
  567. uint256 isGt;
  568. unchecked {
  569. isGt = SafeCast.toUint(value > (1 << 128) - 1);
  570. value >>= isGt * 128;
  571. result += isGt * 16;
  572. isGt = SafeCast.toUint(value > (1 << 64) - 1);
  573. value >>= isGt * 64;
  574. result += isGt * 8;
  575. isGt = SafeCast.toUint(value > (1 << 32) - 1);
  576. value >>= isGt * 32;
  577. result += isGt * 4;
  578. isGt = SafeCast.toUint(value > (1 << 16) - 1);
  579. value >>= isGt * 16;
  580. result += isGt * 2;
  581. result += SafeCast.toUint(value > (1 << 8) - 1);
  582. }
  583. return result;
  584. }
  585. /**
  586. * @dev Return the log in base 256, following the selected rounding direction, of a positive value.
  587. * Returns 0 if given 0.
  588. */
  589. function log256(uint256 value, Rounding rounding) internal pure returns (uint256) {
  590. unchecked {
  591. uint256 result = log256(value);
  592. return result + SafeCast.toUint(unsignedRoundsUp(rounding) && 1 << (result << 3) < value);
  593. }
  594. }
  595. /**
  596. * @dev Returns whether a provided rounding mode is considered rounding up for unsigned integers.
  597. */
  598. function unsignedRoundsUp(Rounding rounding) internal pure returns (bool) {
  599. return uint8(rounding) % 2 == 1;
  600. }
  601. }