From 171e5c2037c11a09eebeb6db011ba2035539baf4 Mon Sep 17 00:00:00 2001 From: Chris Gorman Date: Thu, 29 Jun 2023 15:43:43 -0700 Subject: [PATCH 01/23] Improved integer square root. --- contracts/utils/math/Math.sol | 39 +++++++++++++---------------------- 1 file changed, 14 insertions(+), 25 deletions(-) diff --git a/contracts/utils/math/Math.sol b/contracts/utils/math/Math.sol index d372295d7a8..fc9d0f2a9ad 100644 --- a/contracts/utils/math/Math.sol +++ b/contracts/utils/math/Math.sol @@ -224,39 +224,28 @@ library Math { /** * @dev Returns the square root of a number. If the number is not a perfect square, the value is rounded down. - * - * Inspired by Henry S. Warren, Jr.'s "Hacker's Delight" (Chapter 11). */ - function sqrt(uint256 a) internal pure returns (uint256) { - if (a == 0) { - return 0; - } - - // For our first guess, we get the biggest power of 2 which is smaller than the square root of the target. - // - // We know that the "msb" (most significant bit) of our target number `a` is a power of 2 such that we have - // `msb(a) <= a < 2*msb(a)`. This value can be written `msb(a)=2**k` with `k=log2(a)`. - // - // This can be rewritten `2**log2(a) <= a < 2**(log2(a) + 1)` - // → `sqrt(2**k) <= sqrt(a) < sqrt(2**(k+1))` - // → `2**(k/2) <= sqrt(a) < 2**((k+1)/2) <= 2**(k/2 + 1)` - // - // Consequently, `2**(log2(a) / 2)` is a good first approximation of `sqrt(a)` with at least 1 correct bit. - uint256 result = 1 << (log2(a) >> 1); - - // At this point `result` is an estimation with one bit of precision. We know the true value is a uint128, - // since it is the square root of a uint256. Newton's method converges quadratically (precision doubles at - // every iteration). We thus need at most 7 iteration to turn our partial result with one bit of precision - // into the expected uint128 result. + function sqrt(uint256 a) internal pure returns (uint256 result) { unchecked { + if (a <= 1) { return a; } + if (a >= ((1 << 128) - 1)**2) { return (1 << 128) - 1; } + uint256 aAux = a; + result = 1; + if (aAux >= (1 << 128)) { aAux >>= 128; result <<= 64; } + if (aAux >= (1 << 64 )) { aAux >>= 64; result <<= 32; } + if (aAux >= (1 << 32 )) { aAux >>= 32; result <<= 16; } + if (aAux >= (1 << 16 )) { aAux >>= 16; result <<= 8; } + if (aAux >= (1 << 8 )) { aAux >>= 8; result <<= 4; } + if (aAux >= (1 << 4 )) { aAux >>= 4; result <<= 2; } + if (aAux >= (1 << 2 )) { result <<= 1; } + result += (result >> 1); result = (result + a / result) >> 1; result = (result + a / result) >> 1; result = (result + a / result) >> 1; result = (result + a / result) >> 1; result = (result + a / result) >> 1; result = (result + a / result) >> 1; - result = (result + a / result) >> 1; - return min(result, a / result); + return result * result <= a ? result : (result - 1); } } From beec0572a5a0590cc2d657ddcae9eb2eda1aa0a7 Mon Sep 17 00:00:00 2001 From: Chris Gorman Date: Sat, 22 Jul 2023 14:12:05 -0700 Subject: [PATCH 02/23] Small updates to sqrt function. --- contracts/utils/math/Math.sol | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/contracts/utils/math/Math.sol b/contracts/utils/math/Math.sol index fc9d0f2a9ad..4bdc04de3c9 100644 --- a/contracts/utils/math/Math.sol +++ b/contracts/utils/math/Math.sol @@ -225,12 +225,12 @@ library Math { /** * @dev Returns the square root of a number. If the number is not a perfect square, the value is rounded down. */ - function sqrt(uint256 a) internal pure returns (uint256 result) { + function sqrt(uint256 a) internal pure returns (uint256) { unchecked { if (a <= 1) { return a; } if (a >= ((1 << 128) - 1)**2) { return (1 << 128) - 1; } uint256 aAux = a; - result = 1; + uint256 result = 1; if (aAux >= (1 << 128)) { aAux >>= 128; result <<= 64; } if (aAux >= (1 << 64 )) { aAux >>= 64; result <<= 32; } if (aAux >= (1 << 32 )) { aAux >>= 32; result <<= 16; } @@ -245,7 +245,10 @@ library Math { result = (result + a / result) >> 1; result = (result + a / result) >> 1; result = (result + a / result) >> 1; - return result * result <= a ? result : (result - 1); + if (result * result <= x) { + return result; + } + return result-1; } } From b6c29b166e5ddf08bfd43515556aee93e3dcf4df Mon Sep 17 00:00:00 2001 From: Chris Gorman Date: Mon, 24 Jul 2023 09:48:03 -0700 Subject: [PATCH 03/23] Fixed variable naming. --- contracts/utils/math/Math.sol | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/contracts/utils/math/Math.sol b/contracts/utils/math/Math.sol index 75fdcd8b5f4..84b0fb5e420 100644 --- a/contracts/utils/math/Math.sol +++ b/contracts/utils/math/Math.sol @@ -239,7 +239,7 @@ library Math { result = (result + a / result) >> 1; result = (result + a / result) >> 1; result = (result + a / result) >> 1; - if (result * result <= x) { + if (result * result <= a) { return result; } return result-1; From 0d2bbb3815aec20524dc8e518f90849168c9e43c Mon Sep 17 00:00:00 2001 From: Chris Gorman Date: Mon, 24 Jul 2023 14:54:20 -0700 Subject: [PATCH 04/23] Added some comments to help clarify the algorithm. --- contracts/utils/math/Math.sol | 22 ++++++++++++++++++++-- 1 file changed, 20 insertions(+), 2 deletions(-) diff --git a/contracts/utils/math/Math.sol b/contracts/utils/math/Math.sol index 84b0fb5e420..2c2db530778 100644 --- a/contracts/utils/math/Math.sol +++ b/contracts/utils/math/Math.sol @@ -216,13 +216,25 @@ library Math { /** * @dev Returns the square root of a number. If the number is not a perfect square, the value is rounded * towards zero. - * - * Inspired by Henry S. Warren, Jr.'s "Hacker's Delight" (Chapter 11). */ function sqrt(uint256 a) internal pure returns (uint256) { unchecked { + // Take care of easy edge cases if (a <= 1) { return a; } + // This check ensures no overflow if (a >= ((1 << 128) - 1)**2) { return (1 << 128) - 1; } + + // If we have + // + // 2^{e-1} <= sqrt(x) < 2^{e}, + // + // then at the end of initialization, we will have + // + // result == 2^{e-1} + 2^{e-2}. + // + // This ensures that + // + // abs(sqrt(x) - result) <= 2^{e-2}. uint256 aAux = a; uint256 result = 1; if (aAux >= (1 << 128)) { aAux >>= 128; result <<= 64; } @@ -233,12 +245,18 @@ library Math { if (aAux >= (1 << 4 )) { aAux >>= 4; result <<= 2; } if (aAux >= (1 << 2 )) { result <<= 1; } result += (result >> 1); + + // Perform the 6 required Newton iteration result = (result + a / result) >> 1; result = (result + a / result) >> 1; result = (result + a / result) >> 1; result = (result + a / result) >> 1; result = (result + a / result) >> 1; result = (result + a / result) >> 1; + + // We either have + // + // Isqrt(x) == result or Isqrt(x) == result-1. if (result * result <= a) { return result; } From 1c8e8abe47d2286b25173a1ae3c8e1d3d60e5d48 Mon Sep 17 00:00:00 2001 From: Chris Gorman Date: Tue, 25 Jul 2023 09:55:33 -0700 Subject: [PATCH 05/23] Updated edge case for clarity. --- contracts/utils/math/Math.sol | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/contracts/utils/math/Math.sol b/contracts/utils/math/Math.sol index 2c2db530778..edebe3e4202 100644 --- a/contracts/utils/math/Math.sol +++ b/contracts/utils/math/Math.sol @@ -220,9 +220,13 @@ library Math { function sqrt(uint256 a) internal pure returns (uint256) { unchecked { // Take care of easy edge cases - if (a <= 1) { return a; } - // This check ensures no overflow - if (a >= ((1 << 128) - 1)**2) { return (1 << 128) - 1; } + if (a <= 1) { + return a; + } + // This check ensures no overflow at return + if (a >= uint256(type(uint128).max)**2) { + return type(uint128).max; + } // If we have // @@ -246,7 +250,7 @@ library Math { if (aAux >= (1 << 2 )) { result <<= 1; } result += (result >> 1); - // Perform the 6 required Newton iteration + // Perform the 6 required Newton iterations result = (result + a / result) >> 1; result = (result + a / result) >> 1; result = (result + a / result) >> 1; From 8e2ab78af47bca93a2c953e3629e5d1332ceb66a Mon Sep 17 00:00:00 2001 From: Chris Gorman Date: Wed, 30 Aug 2023 22:46:57 -0700 Subject: [PATCH 06/23] Made small improvements. --- contracts/utils/math/Math.sol | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/contracts/utils/math/Math.sol b/contracts/utils/math/Math.sol index edebe3e4202..9f49fe544b7 100644 --- a/contracts/utils/math/Math.sol +++ b/contracts/utils/math/Math.sol @@ -241,14 +241,14 @@ library Math { // abs(sqrt(x) - result) <= 2^{e-2}. uint256 aAux = a; uint256 result = 1; - if (aAux >= (1 << 128)) { aAux >>= 128; result <<= 64; } + if (aAux >= (1 << 128)) { aAux >>= 128; result = 1 << 64; } if (aAux >= (1 << 64 )) { aAux >>= 64; result <<= 32; } if (aAux >= (1 << 32 )) { aAux >>= 32; result <<= 16; } if (aAux >= (1 << 16 )) { aAux >>= 16; result <<= 8; } if (aAux >= (1 << 8 )) { aAux >>= 8; result <<= 4; } if (aAux >= (1 << 4 )) { aAux >>= 4; result <<= 2; } if (aAux >= (1 << 2 )) { result <<= 1; } - result += (result >> 1); + result = (3 * result) >> 1; // Perform the 6 required Newton iterations result = (result + a / result) >> 1; From 8b3f15f9633d822bf17244f917dde017f484b2c7 Mon Sep 17 00:00:00 2001 From: Chris Gorman Date: Mon, 12 Feb 2024 18:01:12 -0700 Subject: [PATCH 07/23] Added additional comments to include associated proof. --- contracts/utils/math/Math.sol | 76 ++++++++++++++++++++++++++++++++--- 1 file changed, 71 insertions(+), 5 deletions(-) diff --git a/contracts/utils/math/Math.sol b/contracts/utils/math/Math.sol index 3a8857e49b8..f7cc5a088b1 100644 --- a/contracts/utils/math/Math.sol +++ b/contracts/utils/math/Math.sol @@ -221,21 +221,27 @@ library Math { /** * @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 + // Take care of easy edge cases when a == 0 or a == 1 if (a <= 1) { return a; } - // This check ensures no overflow at return + // This check ensures no overflow at return; see below. if (a >= uint256(type(uint128).max)**2) { return type(uint128).max; } // If we have // - // 2^{e-1} <= sqrt(x) < 2^{e}, + // 2^{2e-2} <= a < 2^{2e} + // + // or equivalently + // + // 2^{e-1} <= sqrt(a) < 2^{e}, // // then at the end of initialization, we will have // @@ -243,7 +249,10 @@ library Math { // // This ensures that // - // abs(sqrt(x) - result) <= 2^{e-2}. + // abs(result - sqrt(a)) <= 2^{e-2}. + // + // Using this choice for initialization ensures that only + // 6 Newton iterations are required. uint256 aAux = a; uint256 result = 1; if (aAux >= (1 << 128)) { aAux >>= 128; result = 1 << 64; } @@ -253,19 +262,76 @@ library Math { if (aAux >= (1 << 8 )) { aAux >>= 8; result <<= 4; } if (aAux >= (1 << 4 )) { aAux >>= 4; result <<= 2; } if (aAux >= (1 << 2 )) { result <<= 1; } + // At this point, we have + // result == 2^{e-1}. result = (3 * result) >> 1; + // After this update, we have + // result == 2^{e-1} + 2^{e-2}. // Perform the 6 required Newton iterations + // + // We let + // + // f(x) = (x + a/x)/2. + // + // In this case, if + // + // eps = x - sqrt(a), + // + // then + // + // f(x) - sqrt(a) = (eps^2)/(2x) + // + // This allows for the computation of error bounds + // when Newton's method is used with real numbers. + // When using real numbers, we will always have + // + // f(x) >= sqrt(a); + // + // In the situation where these operations are all integer operations, + // we only have + // + // floor(f(x)) >= Isqrt(a), + // + // but this is not a problem + // as we are actually trying to compute Isqrt(a). + // + // From the initialization, we start with + // + // abs(error) := abs(result - sqrt(a)) + // <= 2^{e-2} result = (result + a / result) >> 1; + // error := result - sqrt(a) + // <= 2^{e-4.5} result = (result + a / result) >> 1; + // error := result - sqrt(a) + // <= 2^{e-9} result = (result + a / result) >> 1; + // error := result - sqrt(a) + // <= 2^{e-18} result = (result + a / result) >> 1; + // error := result - sqrt(a) + // <= 2^{e-36} result = (result + a / result) >> 1; + // error := result - sqrt(a) + // <= 2^{e-72} result = (result + a / result) >> 1; + // error := result - sqrt(a) + // <= 2^{e-144} // We either have // - // Isqrt(x) == result or Isqrt(x) == result-1. + // result == Isqrt(a) or result == Isqrt(a) + 1. + // + // The computation of + // + // result * result + // + // could overflow if result == 2^128 + // (only possible when a >= (uint128.max)^2 and + // result == uint128.max + 1), + // but this is not possible because that edge case + // was taken care of by previous logic. if (result * result <= a) { return result; } From bd9d969e0313602b554933fcb482cc58a0e93f62 Mon Sep 17 00:00:00 2001 From: ernestognw Date: Tue, 13 Feb 2024 21:20:59 -0600 Subject: [PATCH 08/23] Reduce sqrt comments --- contracts/utils/math/Math.sol | 124 ++++++++-------------------------- 1 file changed, 29 insertions(+), 95 deletions(-) diff --git a/contracts/utils/math/Math.sol b/contracts/utils/math/Math.sol index eced28e222d..f726b45ad00 100644 --- a/contracts/utils/math/Math.sol +++ b/contracts/utils/math/Math.sol @@ -338,8 +338,9 @@ library Math { /** * @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. + * + * 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 { @@ -352,26 +353,13 @@ library Math { return type(uint128).max; } - // If we have - // - // 2^{2e-2} <= a < 2^{2e} - // - // or equivalently - // - // 2^{e-1} <= sqrt(a) < 2^{e}, - // - // then at the end of initialization, we will have - // - // result == 2^{e-1} + 2^{e-2}. - // - // This ensures that - // - // abs(result - sqrt(a)) <= 2^{e-2}. - // - // Using this choice for initialization ensures that only - // 6 Newton iterations are required. uint256 aAux = a; - uint256 result = 1; + uint256 result = 1; // Initial approximate bit length. + + // For the first guess of `result` (e), we get the biggest power of 2 which is smaller than sqrt(a) + // (i.e. 2^e <= sqrt(a)). We know that e is at most 127 given (2^128)^2 overflows an uint256. + // Thus, we approximate e by iterating 2^{i/2} where i starts at 128, and applying the exponent + // to e if the result is still smaller than a (up to e == 127). if (aAux >= (1 << 128)) { aAux >>= 128; result = 1 << 64; } if (aAux >= (1 << 64 )) { aAux >>= 64; result <<= 32; } if (aAux >= (1 << 32 )) { aAux >>= 32; result <<= 16; } @@ -379,80 +367,26 @@ library Math { if (aAux >= (1 << 8 )) { aAux >>= 8; result <<= 4; } if (aAux >= (1 << 4 )) { aAux >>= 4; result <<= 2; } if (aAux >= (1 << 2 )) { result <<= 1; } - // At this point, we have - // result == 2^{e-1}. - result = (3 * result) >> 1; - // After this update, we have - // result == 2^{e-1} + 2^{e-2}. - - // Perform the 6 required Newton iterations - // - // We let - // - // f(x) = (x + a/x)/2. - // - // In this case, if - // - // eps = x - sqrt(a), - // - // then - // - // f(x) - sqrt(a) = (eps^2)/(2x) - // - // This allows for the computation of error bounds - // when Newton's method is used with real numbers. - // When using real numbers, we will always have - // - // f(x) >= sqrt(a); - // - // In the situation where these operations are all integer operations, - // we only have - // - // floor(f(x)) >= Isqrt(a), - // - // but this is not a problem - // as we are actually trying to compute Isqrt(a). - // - // From the initialization, we start with - // - // abs(error) := abs(result - sqrt(a)) - // <= 2^{e-2} - result = (result + a / result) >> 1; - // error := result - sqrt(a) - // <= 2^{e-4.5} - result = (result + a / result) >> 1; - // error := result - sqrt(a) - // <= 2^{e-9} - result = (result + a / result) >> 1; - // error := result - sqrt(a) - // <= 2^{e-18} - result = (result + a / result) >> 1; - // error := result - sqrt(a) - // <= 2^{e-36} - result = (result + a / result) >> 1; - // error := result - sqrt(a) - // <= 2^{e-72} - result = (result + a / result) >> 1; - // error := result - sqrt(a) - // <= 2^{e-144} - - // We either have - // - // result == Isqrt(a) or result == Isqrt(a) + 1. - // - // The computation of - // - // result * result - // - // could overflow if result == 2^128 - // (only possible when a >= (uint128.max)^2 and - // result == uint128.max + 1), - // but this is not possible because that edge case - // was taken care of by previous logic. - if (result * result <= a) { - return result; - } - return result-1; + + // 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 (result * 2), ensuring that result - sqrt(a) <= 2^{e-2}. + result = (3 * result) >> 1; // (result + result << 1) / 2 + + // Each Newton iteration will have f(x) = (x + a / x) / 2. + // Given the error (ε) is defined by x - sqrt(a), then we know that + // ε + 1 == ε^2 / 2x <= ε^2 / 2 * sqrt(a). + result = (result + a / result) >> 1; // err := result - sqrt(a) <= 2^{e-4.5} + result = (result + a / result) >> 1; // err := result - sqrt(a) <= 2^{e-9} + result = (result + a / result) >> 1; // err := result - sqrt(a) <= 2^{e-18} + result = (result + a / result) >> 1; // err := result - sqrt(a) <= 2^{e-36} + result = (result + a / result) >> 1; // err := result - sqrt(a)<= 2^{e-72} + result = (result + a / result) >> 1; // err := result - sqrt(a) <= 2^{e-144} + + // After 6 iterations, no more precision can be obtained since the max result is 127. + // Squaring result could overflow if a >= type(uint128).max^2, case discarded at the start. + // result is either sqrt(a) or sqrt(a) + 1. + return result * result <= a ? result : result - 1; } } From dfe2fbcd4a8c03f5c6e99f785314b8a95ea7a60d Mon Sep 17 00:00:00 2001 From: ernestognw Date: Tue, 13 Feb 2024 21:21:25 -0600 Subject: [PATCH 09/23] Lint --- contracts/utils/math/Math.sol | 40 ++++++++++++++++++++++++++--------- 1 file changed, 30 insertions(+), 10 deletions(-) diff --git a/contracts/utils/math/Math.sol b/contracts/utils/math/Math.sol index f726b45ad00..effb6bc2073 100644 --- a/contracts/utils/math/Math.sol +++ b/contracts/utils/math/Math.sol @@ -338,7 +338,7 @@ library Math { /** * @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. */ @@ -349,7 +349,7 @@ library Math { return a; } // This check ensures no overflow at return; see below. - if (a >= uint256(type(uint128).max)**2) { + if (a >= uint256(type(uint128).max) ** 2) { return type(uint128).max; } @@ -360,13 +360,33 @@ library Math { // (i.e. 2^e <= sqrt(a)). We know that e is at most 127 given (2^128)^2 overflows an uint256. // Thus, we approximate e by iterating 2^{i/2} where i starts at 128, and applying the exponent // to e if the result is still smaller than a (up to e == 127). - if (aAux >= (1 << 128)) { aAux >>= 128; result = 1 << 64; } - if (aAux >= (1 << 64 )) { aAux >>= 64; result <<= 32; } - if (aAux >= (1 << 32 )) { aAux >>= 32; result <<= 16; } - if (aAux >= (1 << 16 )) { aAux >>= 16; result <<= 8; } - if (aAux >= (1 << 8 )) { aAux >>= 8; result <<= 4; } - if (aAux >= (1 << 4 )) { aAux >>= 4; result <<= 2; } - if (aAux >= (1 << 2 )) { result <<= 1; } + if (aAux >= (1 << 128)) { + aAux >>= 128; + result = 1 << 64; + } + if (aAux >= (1 << 64)) { + aAux >>= 64; + result <<= 32; + } + if (aAux >= (1 << 32)) { + aAux >>= 32; + result <<= 16; + } + if (aAux >= (1 << 16)) { + aAux >>= 16; + result <<= 8; + } + if (aAux >= (1 << 8)) { + aAux >>= 8; + result <<= 4; + } + if (aAux >= (1 << 4)) { + aAux >>= 4; + result <<= 2; + } + if (aAux >= (1 << 2)) { + result <<= 1; + } // We can use the fact that 2^e <= sqrt(a) to improve the estimation // by computing the arithmetic mean between the current estimation and @@ -382,7 +402,7 @@ library Math { result = (result + a / result) >> 1; // err := result - sqrt(a) <= 2^{e-36} result = (result + a / result) >> 1; // err := result - sqrt(a)<= 2^{e-72} result = (result + a / result) >> 1; // err := result - sqrt(a) <= 2^{e-144} - + // After 6 iterations, no more precision can be obtained since the max result is 127. // Squaring result could overflow if a >= type(uint128).max^2, case discarded at the start. // result is either sqrt(a) or sqrt(a) + 1. From a53d378cee1bce90bc3d46408ad2084f120c09c9 Mon Sep 17 00:00:00 2001 From: Chris Gorman Date: Tue, 13 Feb 2024 21:25:02 -0700 Subject: [PATCH 10/23] Correcting return logic --- contracts/utils/math/Math.sol | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/contracts/utils/math/Math.sol b/contracts/utils/math/Math.sol index effb6bc2073..7aae82d8b24 100644 --- a/contracts/utils/math/Math.sol +++ b/contracts/utils/math/Math.sol @@ -406,7 +406,10 @@ library Math { // After 6 iterations, no more precision can be obtained since the max result is 127. // Squaring result could overflow if a >= type(uint128).max^2, case discarded at the start. // result is either sqrt(a) or sqrt(a) + 1. - return result * result <= a ? result : result - 1; + if (result * result <= a) { + return result; + } + return result-1; } } From 6c7401231fd312c37454259c4c7351108a602e66 Mon Sep 17 00:00:00 2001 From: Hadrien Croubois Date: Wed, 14 Feb 2024 16:44:11 +0100 Subject: [PATCH 11/23] remove overflow check (early) in favor of overflow save math at the end. This optimized gas cost for 'small' values that are not catched by the initial check. --- contracts/utils/math/Math.sol | 16 ++++------------ 1 file changed, 4 insertions(+), 12 deletions(-) diff --git a/contracts/utils/math/Math.sol b/contracts/utils/math/Math.sol index 7aae82d8b24..1ce44cd939b 100644 --- a/contracts/utils/math/Math.sol +++ b/contracts/utils/math/Math.sol @@ -348,13 +348,9 @@ library Math { if (a <= 1) { return a; } - // This check ensures no overflow at return; see below. - if (a >= uint256(type(uint128).max) ** 2) { - return type(uint128).max; - } uint256 aAux = a; - uint256 result = 1; // Initial approximate bit length. + uint256 result = 1; // For the first guess of `result` (e), we get the biggest power of 2 which is smaller than sqrt(a) // (i.e. 2^e <= sqrt(a)). We know that e is at most 127 given (2^128)^2 overflows an uint256. @@ -362,7 +358,7 @@ library Math { // to e if the result is still smaller than a (up to e == 127). if (aAux >= (1 << 128)) { aAux >>= 128; - result = 1 << 64; + result <<= 64; } if (aAux >= (1 << 64)) { aAux >>= 64; @@ -391,7 +387,7 @@ library Math { // 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 (result * 2), ensuring that result - sqrt(a) <= 2^{e-2}. - result = (3 * result) >> 1; // (result + result << 1) / 2 + result = (3 * result) >> 1; // Each Newton iteration will have f(x) = (x + a / x) / 2. // Given the error (ε) is defined by x - sqrt(a), then we know that @@ -404,12 +400,8 @@ library Math { result = (result + a / result) >> 1; // err := result - sqrt(a) <= 2^{e-144} // After 6 iterations, no more precision can be obtained since the max result is 127. - // Squaring result could overflow if a >= type(uint128).max^2, case discarded at the start. // result is either sqrt(a) or sqrt(a) + 1. - if (result * result <= a) { - return result; - } - return result-1; + return result - SafeCast.toUint(result > a / result); } } From fe8abc3faeece8c55a3ebe27595471fb8dccb8bd Mon Sep 17 00:00:00 2001 From: Hadrien Croubois Date: Wed, 14 Feb 2024 16:57:08 +0100 Subject: [PATCH 12/23] format --- contracts/utils/math/Math.sol | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/contracts/utils/math/Math.sol b/contracts/utils/math/Math.sol index 1ce44cd939b..5ae79be098a 100644 --- a/contracts/utils/math/Math.sol +++ b/contracts/utils/math/Math.sol @@ -396,7 +396,7 @@ library Math { result = (result + a / result) >> 1; // err := result - sqrt(a) <= 2^{e-9} result = (result + a / result) >> 1; // err := result - sqrt(a) <= 2^{e-18} result = (result + a / result) >> 1; // err := result - sqrt(a) <= 2^{e-36} - result = (result + a / result) >> 1; // err := result - sqrt(a)<= 2^{e-72} + result = (result + a / result) >> 1; // err := result - sqrt(a) <= 2^{e-72} result = (result + a / result) >> 1; // err := result - sqrt(a) <= 2^{e-144} // After 6 iterations, no more precision can be obtained since the max result is 127. From 382ab6fa1b631ba97c91e62edc79c2361eb99289 Mon Sep 17 00:00:00 2001 From: ernestognw Date: Wed, 14 Feb 2024 18:09:06 -0600 Subject: [PATCH 13/23] Iterate on comments --- contracts/utils/math/Math.sol | 52 +++++++++++++++++++---------------- 1 file changed, 28 insertions(+), 24 deletions(-) diff --git a/contracts/utils/math/Math.sol b/contracts/utils/math/Math.sol index 5ae79be098a..05c4cd0baa7 100644 --- a/contracts/utils/math/Math.sol +++ b/contracts/utils/math/Math.sol @@ -352,55 +352,59 @@ library Math { uint256 aAux = a; uint256 result = 1; - // For the first guess of `result` (e), we get the biggest power of 2 which is smaller than sqrt(a) - // (i.e. 2^e <= sqrt(a)). We know that e is at most 127 given (2^128)^2 overflows an uint256. - // Thus, we approximate e by iterating 2^{i/2} where i starts at 128, and applying the exponent - // to e if the result is still smaller than a (up to e == 127). + // 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 approximate 2**n by adding 2**e to our result and substracting + // it from the target if 2**e is still less than its square root + // (i.e. 2**e <= sqrt(aAux)) for all e/2 from 128 to 0. + // We know that e is at most 127 because (2**128)**2 = 2**256is bigger than any uint256. if (aAux >= (1 << 128)) { aAux >>= 128; - result <<= 64; + result <<= 64; // e/2 = 64 } if (aAux >= (1 << 64)) { aAux >>= 64; - result <<= 32; + result <<= 32; // e/2 = 32 } if (aAux >= (1 << 32)) { aAux >>= 32; - result <<= 16; + result <<= 16; // e/2 = 16 } if (aAux >= (1 << 16)) { aAux >>= 16; - result <<= 8; + result <<= 8; // e/2 = 8 } if (aAux >= (1 << 8)) { aAux >>= 8; - result <<= 4; + result <<= 4; // e/2 = 4 } if (aAux >= (1 << 4)) { aAux >>= 4; - result <<= 2; + result <<= 2; // e/2 = 2 } if (aAux >= (1 << 2)) { result <<= 1; } - // We can use the fact that 2^e <= sqrt(a) to improve the estimation + // 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 (result * 2), ensuring that result - sqrt(a) <= 2^{e-2}. + // the next one (result * 2), ensuring that result - sqrt(a) <= 2**(e-2). result = (3 * result) >> 1; - // Each Newton iteration will have f(x) = (x + a / x) / 2. - // Given the error (ε) is defined by x - sqrt(a), then we know that - // ε + 1 == ε^2 / 2x <= ε^2 / 2 * sqrt(a). - result = (result + a / result) >> 1; // err := result - sqrt(a) <= 2^{e-4.5} - result = (result + a / result) >> 1; // err := result - sqrt(a) <= 2^{e-9} - result = (result + a / result) >> 1; // err := result - sqrt(a) <= 2^{e-18} - result = (result + a / result) >> 1; // err := result - sqrt(a) <= 2^{e-36} - result = (result + a / result) >> 1; // err := result - sqrt(a) <= 2^{e-72} - result = (result + a / result) >> 1; // err := result - sqrt(a) <= 2^{e-144} - - // After 6 iterations, no more precision can be obtained since the max result is 127. - // result is either sqrt(a) or sqrt(a) + 1. + // We define the error as ε = result - sqrt(a) and we've shown that ε <= 2**(e-2). + // Then, we know that ε + 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); } } From c60542b18ba75b672807fade0b5dc631c20268df Mon Sep 17 00:00:00 2001 From: ernestognw Date: Wed, 14 Feb 2024 19:06:46 -0600 Subject: [PATCH 14/23] Iterate on comments --- contracts/utils/math/Math.sol | 13 ++++++++----- 1 file changed, 8 insertions(+), 5 deletions(-) diff --git a/contracts/utils/math/Math.sol b/contracts/utils/math/Math.sol index 05c4cd0baa7..81b36a9fb03 100644 --- a/contracts/utils/math/Math.sol +++ b/contracts/utils/math/Math.sol @@ -357,7 +357,7 @@ library Math { // We approximate 2**n by adding 2**e to our result and substracting // it from the target if 2**e is still less than its square root // (i.e. 2**e <= sqrt(aAux)) for all e/2 from 128 to 0. - // We know that e is at most 127 because (2**128)**2 = 2**256is bigger than any uint256. + // We know that e is at most 127 because (2**128)**2 = 2**256 is bigger than any uint256. if (aAux >= (1 << 128)) { aAux >>= 128; result <<= 64; // e/2 = 64 @@ -388,13 +388,16 @@ library Math { // 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 (result * 2), ensuring that result - sqrt(a) <= 2**(e-2). + // the next one (result * 2). result = (3 * result) >> 1; - // We define the error as ε = result - sqrt(a) and we've shown that ε <= 2**(e-2). - // Then, we know that ε + 1 == ε**2 / 2x <= ε**2 / 2 * sqrt(a) as shown in - // Walter Rudin. Principles of Mathematical Analysis. + // 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) + + // Then, we know that ε + 1 == ε**2 / 2x <= ε**2 / 2 * sqrt(a) as shown in 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) From ddea29227a7c9f2d060dfe2c7c983b02c079a92f Mon Sep 17 00:00:00 2001 From: ernestognw Date: Wed, 14 Feb 2024 19:10:31 -0600 Subject: [PATCH 15/23] Lint --- contracts/utils/math/Math.sol | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/contracts/utils/math/Math.sol b/contracts/utils/math/Math.sol index 81b36a9fb03..49cb1504513 100644 --- a/contracts/utils/math/Math.sol +++ b/contracts/utils/math/Math.sol @@ -393,7 +393,7 @@ library Math { // 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) + // 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) From 360f65d0b88453cfb4a581c4271eb3b55b693311 Mon Sep 17 00:00:00 2001 From: ernestognw Date: Wed, 14 Feb 2024 19:15:09 -0600 Subject: [PATCH 16/23] Codespell --- contracts/utils/math/Math.sol | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/contracts/utils/math/Math.sol b/contracts/utils/math/Math.sol index 49cb1504513..bebdd203a73 100644 --- a/contracts/utils/math/Math.sol +++ b/contracts/utils/math/Math.sol @@ -354,7 +354,7 @@ library Math { // 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 approximate 2**n by adding 2**e to our result and substracting + // We approximate 2**n by adding 2**e to our result and subtracting // it from the target if 2**e is still less than its square root // (i.e. 2**e <= sqrt(aAux)) for all e/2 from 128 to 0. // We know that e is at most 127 because (2**128)**2 = 2**256 is bigger than any uint256. From 4288c743760c7403114b0233d45b4610b2121549 Mon Sep 17 00:00:00 2001 From: ernestognw Date: Wed, 14 Feb 2024 20:15:59 -0600 Subject: [PATCH 17/23] Iterate on comments --- contracts/utils/math/Math.sol | 43 +++++++++++++++++------------------ 1 file changed, 21 insertions(+), 22 deletions(-) diff --git a/contracts/utils/math/Math.sol b/contracts/utils/math/Math.sol index bebdd203a73..33784704bc2 100644 --- a/contracts/utils/math/Math.sol +++ b/contracts/utils/math/Math.sol @@ -353,57 +353,56 @@ library Math { 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 approximate 2**n by adding 2**e to our result and subtracting - // it from the target if 2**e is still less than its square root - // (i.e. 2**e <= sqrt(aAux)) for all e/2 from 128 to 0. - // We know that e is at most 127 because (2**128)**2 = 2**256 is bigger than any uint256. + // 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; // e/2 = 64 + result <<= 64; // sqrt(a) >= 2**(e/2) } if (aAux >= (1 << 64)) { aAux >>= 64; - result <<= 32; // e/2 = 32 + result <<= 32; // sqrt(a) >= 2**(e/2) } if (aAux >= (1 << 32)) { aAux >>= 32; - result <<= 16; // e/2 = 16 + result <<= 16; // sqrt(a) >= 2**(e/2) } if (aAux >= (1 << 16)) { aAux >>= 16; - result <<= 8; // e/2 = 8 + result <<= 8; // sqrt(a) >= 2**(e/2) } if (aAux >= (1 << 8)) { aAux >>= 8; - result <<= 4; // e/2 = 4 + result <<= 4; // sqrt(a) >= 2**(e/2) } if (aAux >= (1 << 4)) { aAux >>= 4; - result <<= 2; // e/2 = 2 + result <<= 2; // sqrt(a) >= 2**(e/2) } if (aAux >= (1 << 2)) { result <<= 1; } - // 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 (result * 2). + // 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) + // 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) - // Then, we know that ε + 1 == ε**2 / 2x <= ε**2 / 2 * sqrt(a) as shown in - 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) + 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 From 256189af89d43c6bcc54c0dab6d3a5f418ecd930 Mon Sep 17 00:00:00 2001 From: Hadrien Croubois Date: Thu, 15 Feb 2024 11:08:05 +0100 Subject: [PATCH 18/23] rewrite most comments of the sqrt function --- contracts/utils/math/Math.sol | 128 ++++++++++++++++++++-------------- 1 file changed, 77 insertions(+), 51 deletions(-) diff --git a/contracts/utils/math/Math.sol b/contracts/utils/math/Math.sol index 33784704bc2..2f9a0ef7e3a 100644 --- a/contracts/utils/math/Math.sol +++ b/contracts/utils/math/Math.sol @@ -349,65 +349,91 @@ library Math { 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) + // In this function, we use Newton's method to get a root of `f(x) := x² - a`. This is also known as + // Heron's method. This involves building a sequence x_n that converges toward sqrt(a). For each + // iteration x_n, we define ε_n = | x_n - sqrt(a) |. This represents the error between the current value + // and the result. + // + // For our first estimation, we get the biggest power of 2 which is smaller than the square root of the + // target. (i.e. x_n = 2**e ≤ sqrt(a) < 2**(e+1)). We know that e cannot be greater than 127 because + // (2¹²⁸)² = 2²⁵⁶ is bigger than any uint256. + // + // By noticing that + // `2**e ≤ sqrt(a) < 2**(e+1) → (2**e)² ≤ a < 2**(e+1) → 2**(2*e) ≤ a < 2**(2*e+2)` + // we can deduce that the `e` we are looking for is `log2(a) / 2`, and can be computed using a method + // similar to the msb function. + uint256 aa = a; + uint256 xn = 1; + + if (aa >= (1 << 128)) { + aa >>= 128; + xn <<= 64; } - if (aAux >= (1 << 64)) { - aAux >>= 64; - result <<= 32; // sqrt(a) >= 2**(e/2) + if (aa >= (1 << 64)) { + aa >>= 64; + xn <<= 32; } - if (aAux >= (1 << 32)) { - aAux >>= 32; - result <<= 16; // sqrt(a) >= 2**(e/2) + if (aa >= (1 << 32)) { + aa >>= 32; + xn <<= 16; } - if (aAux >= (1 << 16)) { - aAux >>= 16; - result <<= 8; // sqrt(a) >= 2**(e/2) + if (aa >= (1 << 16)) { + aa >>= 16; + xn <<= 8; } - if (aAux >= (1 << 8)) { - aAux >>= 8; - result <<= 4; // sqrt(a) >= 2**(e/2) + if (aa >= (1 << 8)) { + aa >>= 8; + xn <<= 4; } - if (aAux >= (1 << 4)) { - aAux >>= 4; - result <<= 2; // sqrt(a) >= 2**(e/2) + if (aa >= (1 << 4)) { + aa >>= 4; + xn <<= 2; } - if (aAux >= (1 << 2)) { - result <<= 1; + if (aa >= (1 << 2)) { + xn <<= 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); + // We now have x_n such that `x_n = 2**e <= sqrt(a) < 2**(e+1) = 2 * x_n`. This implies ε_n < 2**e. + // + // We can refine our estimation by noticing that the the middle of that interval minimizes the error. + // If we move x_n to equal 2**e + 2**(e-1), then we reduce the error to ε_n < 2**(e-1). + // This is going to be or x_0 (and ε_0) + xn = (3 * xn) >> 1; // ε_0 := | x_0 - sqrt(a) | < 2**(e-1) + + // From here, we iterate using Heron's method + // x_{n+1} = (x_n + a / x_n) / 2 + // + // One should note that: + // x_{n+1}² - a = ((x_n + a / x_n) / 2)² - a + // = ((x_n² + a) / (2 * x_n))² - a + // = (x_n⁴ + 2*a*x_n² + a²) / (4 * x_n²) - a + // = (x_n⁴ + 2*a*x_n² + a² - 4*a*x_n²) / (4 * x_n²) + // = (x_n⁴ - 2*a*x_n² a²) / (4 * x_n²) + // = (x_n² - a)² / (2 * x_n)² + // = ((x_n² - a) / (2 * x_n))² + // ≥ 0 + // Which proves that for all n ≥ 1, x_n ≥ sqrt(a) + // + // This gives us the proof of quadratic convergence of the sequence: + // ε_{n+1} = | x_{n+1} - sqrt(a) | + // = | (x_n + a / x_n) / 2 - sqrt(a) | + // = | (x_n² + a - 2*x_n*sqrt(a)) / (2 * x_n) | + // = | (x_n - sqrt(a))² / (2 * x_n) | + // = | ε_n² / (2 * x_n) | + // = ε_n² / | (2 * x_n) | + // ≤ ε_n² + // That last inequality holds because x_0 ≥ 1 by construction, and x_n ≥ sqrt(a) ≥ 1 for n ≥ 1. + xn = (xn + a / xn) >> 1; // ε_1 := | x_1 - sqrt(a) | < 2**(e-4.5) ← todo 4.5? how? + xn = (xn + a / xn) >> 1; // ε_2 := | x_2 - sqrt(a) | < 2**(e-9) + xn = (xn + a / xn) >> 1; // ε_3 := | x_3 - sqrt(a) | < 2**(e-18) + xn = (xn + a / xn) >> 1; // ε_4 := | x_4 - sqrt(a) | < 2**(e-36) + xn = (xn + a / xn) >> 1; // ε_5 := | x_5 - sqrt(a) | < 2**(e-72) + xn = (xn + a / xn) >> 1; // ε_6 := | x_6 - sqrt(a) | < 2**(e-144) + + // Because e < 128 (as discussed during the first estimation phase), we know have reached a precision + // ε_6 < 2**(e-144) < 1. Given we're operating on integers, then we can ensure that xn is now either + // sqrt(a) or sqrt(a) + 1. + return xn - SafeCast.toUint(xn > a / xn); } } From 6d00bbda079eac6886d268539242f6546969364c Mon Sep 17 00:00:00 2001 From: Hadrien Croubois Date: Thu, 15 Feb 2024 11:36:57 +0100 Subject: [PATCH 19/23] improve --- contracts/utils/math/Math.sol | 28 ++++++++++++++-------------- 1 file changed, 14 insertions(+), 14 deletions(-) diff --git a/contracts/utils/math/Math.sol b/contracts/utils/math/Math.sol index 2f9a0ef7e3a..1055b464f7f 100644 --- a/contracts/utils/math/Math.sol +++ b/contracts/utils/math/Math.sol @@ -393,12 +393,12 @@ library Math { xn <<= 1; } - // We now have x_n such that `x_n = 2**e <= sqrt(a) < 2**(e+1) = 2 * x_n`. This implies ε_n < 2**e. + // We now have x_n such that `x_n = 2**e ≤ sqrt(a) < 2**(e+1) = 2 * x_n`. This implies ε_n ≤ 2**e. // // We can refine our estimation by noticing that the the middle of that interval minimizes the error. - // If we move x_n to equal 2**e + 2**(e-1), then we reduce the error to ε_n < 2**(e-1). - // This is going to be or x_0 (and ε_0) - xn = (3 * xn) >> 1; // ε_0 := | x_0 - sqrt(a) | < 2**(e-1) + // If we move x_n to equal 2**e + 2**(e-1), then we reduce the error to ε_n ≤ 2**(e-1). + // This is going to be our x_0 (and ε_0) + xn = (3 * xn) >> 1; // ε_0 := | x_0 - sqrt(a) | ≤ 2**(e-1) // From here, we iterate using Heron's method // x_{n+1} = (x_n + a / x_n) / 2 @@ -406,9 +406,9 @@ library Math { // One should note that: // x_{n+1}² - a = ((x_n + a / x_n) / 2)² - a // = ((x_n² + a) / (2 * x_n))² - a - // = (x_n⁴ + 2*a*x_n² + a²) / (4 * x_n²) - a - // = (x_n⁴ + 2*a*x_n² + a² - 4*a*x_n²) / (4 * x_n²) - // = (x_n⁴ - 2*a*x_n² a²) / (4 * x_n²) + // = (x_n⁴ + 2 * a * x_n² + a²) / (4 * x_n²) - a + // = (x_n⁴ + 2 * a * x_n² + a² - 4 * a * x_n²) / (4 * x_n²) + // = (x_n⁴ - 2 * a * x_n² + a²) / (4 * x_n²) // = (x_n² - a)² / (2 * x_n)² // = ((x_n² - a) / (2 * x_n))² // ≥ 0 @@ -423,15 +423,15 @@ library Math { // = ε_n² / | (2 * x_n) | // ≤ ε_n² // That last inequality holds because x_0 ≥ 1 by construction, and x_n ≥ sqrt(a) ≥ 1 for n ≥ 1. - xn = (xn + a / xn) >> 1; // ε_1 := | x_1 - sqrt(a) | < 2**(e-4.5) ← todo 4.5? how? - xn = (xn + a / xn) >> 1; // ε_2 := | x_2 - sqrt(a) | < 2**(e-9) - xn = (xn + a / xn) >> 1; // ε_3 := | x_3 - sqrt(a) | < 2**(e-18) - xn = (xn + a / xn) >> 1; // ε_4 := | x_4 - sqrt(a) | < 2**(e-36) - xn = (xn + a / xn) >> 1; // ε_5 := | x_5 - sqrt(a) | < 2**(e-72) - xn = (xn + a / xn) >> 1; // ε_6 := | x_6 - sqrt(a) | < 2**(e-144) + xn = (xn + a / xn) >> 1; // ε_1 := | x_1 - sqrt(a) | ≤ 2**(e-4.5) ← todo 4.5? how? + xn = (xn + a / xn) >> 1; // ε_2 := | x_2 - sqrt(a) | ≤ 2**(e-9) + xn = (xn + a / xn) >> 1; // ε_3 := | x_3 - sqrt(a) | ≤ 2**(e-18) + xn = (xn + a / xn) >> 1; // ε_4 := | x_4 - sqrt(a) | ≤ 2**(e-36) + xn = (xn + a / xn) >> 1; // ε_5 := | x_5 - sqrt(a) | ≤ 2**(e-72) + xn = (xn + a / xn) >> 1; // ε_6 := | x_6 - sqrt(a) | ≤ 2**(e-144) // Because e < 128 (as discussed during the first estimation phase), we know have reached a precision - // ε_6 < 2**(e-144) < 1. Given we're operating on integers, then we can ensure that xn is now either + // ε_6 ≤ 2**(e-144) < 1. Given we're operating on integers, then we can ensure that xn is now either // sqrt(a) or sqrt(a) + 1. return xn - SafeCast.toUint(xn > a / xn); } From 6c8be261f8e79b861e2b3f644d2e26b1ca402138 Mon Sep 17 00:00:00 2001 From: Hadrien Croubois Date: Thu, 15 Feb 2024 11:54:48 +0100 Subject: [PATCH 20/23] =?UTF-8?q?shift=20notation=20e+1=20=E2=86=92=20e?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- contracts/utils/math/Math.sol | 22 +++++++++++----------- 1 file changed, 11 insertions(+), 11 deletions(-) diff --git a/contracts/utils/math/Math.sol b/contracts/utils/math/Math.sol index 1055b464f7f..a92b3cd6c7a 100644 --- a/contracts/utils/math/Math.sol +++ b/contracts/utils/math/Math.sol @@ -351,17 +351,17 @@ library Math { // In this function, we use Newton's method to get a root of `f(x) := x² - a`. This is also known as // Heron's method. This involves building a sequence x_n that converges toward sqrt(a). For each - // iteration x_n, we define ε_n = | x_n - sqrt(a) |. This represents the error between the current value + // iteration x_n, we define `ε_n = | x_n - sqrt(a) |`. This represents the error between the current value // and the result. // - // For our first estimation, we get the biggest power of 2 which is smaller than the square root of the - // target. (i.e. x_n = 2**e ≤ sqrt(a) < 2**(e+1)). We know that e cannot be greater than 127 because - // (2¹²⁸)² = 2²⁵⁶ is bigger than any uint256. + // For our first estimation, we consider `e` the smallest power of 2 which is bigger than the square root + // of the target. (i.e. `2**(e-1) ≤ sqrt(a) < 2**e`). We know that `e ≤ 128` because `(2¹²⁸)² = 2²⁵⁶` is + // bigger than any uint256. // // By noticing that - // `2**e ≤ sqrt(a) < 2**(e+1) → (2**e)² ≤ a < 2**(e+1) → 2**(2*e) ≤ a < 2**(2*e+2)` - // we can deduce that the `e` we are looking for is `log2(a) / 2`, and can be computed using a method - // similar to the msb function. + // `2**(e-1) ≤ sqrt(a) < 2**e → (2**(e-1))² ≤ a < (2**e)² → 2**(2*e-2) ≤ a < 2**(2*e)` + // we can deduce that `e - 1` is `log2(a) / 2`. We can thus compute `x_n = 2**(e-1)` using a method similar + // to the msb function. uint256 aa = a; uint256 xn = 1; @@ -393,12 +393,12 @@ library Math { xn <<= 1; } - // We now have x_n such that `x_n = 2**e ≤ sqrt(a) < 2**(e+1) = 2 * x_n`. This implies ε_n ≤ 2**e. + // 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). // // We can refine our estimation by noticing that the the middle of that interval minimizes the error. - // If we move x_n to equal 2**e + 2**(e-1), then we reduce the error to ε_n ≤ 2**(e-1). + // If we move x_n to equal 2**(e-1) + 2**(e-2), then we reduce the error to ε_n ≤ 2**(e-2). // This is going to be our x_0 (and ε_0) - xn = (3 * xn) >> 1; // ε_0 := | x_0 - sqrt(a) | ≤ 2**(e-1) + xn = (3 * xn) >> 1; // ε_0 := | x_0 - sqrt(a) | ≤ 2**(e-2) // From here, we iterate using Heron's method // x_{n+1} = (x_n + a / x_n) / 2 @@ -430,7 +430,7 @@ library Math { xn = (xn + a / xn) >> 1; // ε_5 := | x_5 - sqrt(a) | ≤ 2**(e-72) xn = (xn + a / xn) >> 1; // ε_6 := | x_6 - sqrt(a) | ≤ 2**(e-144) - // Because e < 128 (as discussed during the first estimation phase), we know have reached a precision + // Because e ≤ 128 (as discussed during the first estimation phase), we know have reached a precision // ε_6 ≤ 2**(e-144) < 1. Given we're operating on integers, then we can ensure that xn is now either // sqrt(a) or sqrt(a) + 1. return xn - SafeCast.toUint(xn > a / xn); From 67d032858706c4fb964834ebcb2c3f889ae60bf7 Mon Sep 17 00:00:00 2001 From: Hadrien Croubois Date: Thu, 15 Feb 2024 13:13:00 +0100 Subject: [PATCH 21/23] complete proof --- contracts/utils/math/Math.sol | 29 ++++++++++++++++++++--------- 1 file changed, 20 insertions(+), 9 deletions(-) diff --git a/contracts/utils/math/Math.sol b/contracts/utils/math/Math.sol index a92b3cd6c7a..4de27652d3e 100644 --- a/contracts/utils/math/Math.sol +++ b/contracts/utils/math/Math.sol @@ -412,7 +412,7 @@ library Math { // = (x_n² - a)² / (2 * x_n)² // = ((x_n² - a) / (2 * x_n))² // ≥ 0 - // Which proves that for all n ≥ 1, x_n ≥ sqrt(a) + // Which proves that for all n ≥ 1, sqrt(a) ≤ x_n // // This gives us the proof of quadratic convergence of the sequence: // ε_{n+1} = | x_{n+1} - sqrt(a) | @@ -421,14 +421,25 @@ library Math { // = | (x_n - sqrt(a))² / (2 * x_n) | // = | ε_n² / (2 * x_n) | // = ε_n² / | (2 * x_n) | - // ≤ ε_n² - // That last inequality holds because x_0 ≥ 1 by construction, and x_n ≥ sqrt(a) ≥ 1 for n ≥ 1. - xn = (xn + a / xn) >> 1; // ε_1 := | x_1 - sqrt(a) | ≤ 2**(e-4.5) ← todo 4.5? how? - xn = (xn + a / xn) >> 1; // ε_2 := | x_2 - sqrt(a) | ≤ 2**(e-9) - xn = (xn + a / xn) >> 1; // ε_3 := | x_3 - sqrt(a) | ≤ 2**(e-18) - xn = (xn + a / xn) >> 1; // ε_4 := | x_4 - sqrt(a) | ≤ 2**(e-36) - xn = (xn + a / xn) >> 1; // ε_5 := | x_5 - sqrt(a) | ≤ 2**(e-72) - xn = (xn + a / xn) >> 1; // ε_6 := | x_6 - sqrt(a) | ≤ 2**(e-144) + // + // For the first iteration, we have a special case where x_0 is known: + // ε_1 = ε_0² / | (2 * x_0) | + // ≤ (2**(e-2))² / (2 * (2**(e-1) + 2**(e-2))) + // ≤ 2**(2*e-4) / (2 * 3 * 2**(e-2)) + // ≤ 2**(e-3) / 3 + // ≤ 2**(e-4.5) + // + // For the following iterations, we use the fact that, 2**(e-1) ≤ sqrt(a) ≤ x_n + // ε_{n+1} = ε_n² / | (2 * x_n) | + // ≤ (2**(e-k))² / (2 * 2**(e-1)) + // ≤ 2**(2*e-2*k) / 2**e + // ≤ 2**(e-2*k) + xn = (xn + a / xn) >> 1; // ε_1 := | x_1 - sqrt(a) | ≤ 2**(e-4.5) -- special case, see above + xn = (xn + a / xn) >> 1; // ε_2 := | x_2 - sqrt(a) | ≤ 2**(e-9) -- general case with k = 4.5 + xn = (xn + a / xn) >> 1; // ε_3 := | x_3 - sqrt(a) | ≤ 2**(e-18) -- general case with k = 9 + xn = (xn + a / xn) >> 1; // ε_4 := | x_4 - sqrt(a) | ≤ 2**(e-36) -- general case with k = 18 + xn = (xn + a / xn) >> 1; // ε_5 := | x_5 - sqrt(a) | ≤ 2**(e-72) -- general case with k = 36 + xn = (xn + a / xn) >> 1; // ε_6 := | x_6 - sqrt(a) | ≤ 2**(e-144) -- general case with k = 72 // Because e ≤ 128 (as discussed during the first estimation phase), we know have reached a precision // ε_6 ≤ 2**(e-144) < 1. Given we're operating on integers, then we can ensure that xn is now either From 85a56d65ab079d9f35810975f83e57d4af97a6d7 Mon Sep 17 00:00:00 2001 From: Hadrien Croubois Date: Thu, 15 Feb 2024 13:16:38 +0100 Subject: [PATCH 22/23] improvement --- contracts/utils/math/Math.sol | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/contracts/utils/math/Math.sol b/contracts/utils/math/Math.sol index 4de27652d3e..405f881e06d 100644 --- a/contracts/utils/math/Math.sol +++ b/contracts/utils/math/Math.sol @@ -425,8 +425,9 @@ library Math { // For the first iteration, we have a special case where x_0 is known: // ε_1 = ε_0² / | (2 * x_0) | // ≤ (2**(e-2))² / (2 * (2**(e-1) + 2**(e-2))) - // ≤ 2**(2*e-4) / (2 * 3 * 2**(e-2)) + // ≤ 2**(2*e-4) / (3 * 2**(e-1)) // ≤ 2**(e-3) / 3 + // ≤ 2**(e-3-log2(3)) // ≤ 2**(e-4.5) // // For the following iterations, we use the fact that, 2**(e-1) ≤ sqrt(a) ≤ x_n From e748bfcdc8f6246293ec53e52abe0340158b7412 Mon Sep 17 00:00:00 2001 From: Hadrien Croubois Date: Fri, 16 Feb 2024 10:14:35 +0100 Subject: [PATCH 23/23] remove reference to Heron's method to avoid confusion --- contracts/utils/math/Math.sol | 11 +++++------ 1 file changed, 5 insertions(+), 6 deletions(-) diff --git a/contracts/utils/math/Math.sol b/contracts/utils/math/Math.sol index 405f881e06d..1571429900c 100644 --- a/contracts/utils/math/Math.sol +++ b/contracts/utils/math/Math.sol @@ -349,10 +349,9 @@ library Math { return a; } - // In this function, we use Newton's method to get a root of `f(x) := x² - a`. This is also known as - // Heron's method. This involves building a sequence x_n that converges toward sqrt(a). For each - // iteration x_n, we define `ε_n = | x_n - sqrt(a) |`. This represents the error between the current value - // and the result. + // In this function, we use Newton's method to get a root of `f(x) := x² - a`. It involves building a + // sequence x_n that converges toward sqrt(a). For each iteration x_n, we also define the error between + // the current value as `ε_n = | x_n - sqrt(a) |`. // // For our first estimation, we consider `e` the smallest power of 2 which is bigger than the square root // of the target. (i.e. `2**(e-1) ≤ sqrt(a) < 2**e`). We know that `e ≤ 128` because `(2¹²⁸)² = 2²⁵⁶` is @@ -400,7 +399,7 @@ library Math { // This is going to be our x_0 (and ε_0) xn = (3 * xn) >> 1; // ε_0 := | x_0 - sqrt(a) | ≤ 2**(e-2) - // From here, we iterate using Heron's method + // From here, Newton's method give us: // x_{n+1} = (x_n + a / x_n) / 2 // // One should note that: @@ -430,7 +429,7 @@ library Math { // ≤ 2**(e-3-log2(3)) // ≤ 2**(e-4.5) // - // For the following iterations, we use the fact that, 2**(e-1) ≤ sqrt(a) ≤ x_n + // For the following iterations, we use the fact that, 2**(e-1) ≤ sqrt(a) ≤ x_n: // ε_{n+1} = ε_n² / | (2 * x_n) | // ≤ (2**(e-k))² / (2 * 2**(e-1)) // ≤ 2**(2*e-2*k) / 2**e