⚡️Golf sqrt and cbrt by getting a better initial estimate and removing an iteration

Author: duncancmt

src/utils/clz/FixedPointMathLib.sol

@@ -770,18 +770,13 @@ library FixedPointMathLib {
     function sqrt(uint256 x) internal pure returns (uint256 z) {
         /// @solidity memory-safe-assembly
         assembly {
-            // Step 1: Get the bit position of the most significant bit
-            // n = floor(log2(x))
-            // For x ≈ 2^n, we know sqrt(x) ≈ 2^(n/2)
-            // We use (n+1)/2 instead of n/2 to round up slightly
-            // This gives a better initial approximation
-            //
-            // Formula: z = 2^((n+1)/2) = 2^(floor((n+1)/2))
-            // Implemented as: z = 1 << ((n+1) >> 1)
+            // Initial guess z = 2^⌊(n+1)/2⌋ where n = ⌊log₂(x)⌋. The
+            // alternating-endpoint seed gives ε₁ = 0.0607 after one Babylonian
+            // step for all inputs. With ε_{n+1} ≈ ε²/2, 6 steps yield 2⁻¹⁶¹
+            // relative error (>128 correct bits).
             z := shl(shr(1, sub(256, clz(x))), 1)
 
-            /// (x/z + z) / 2
-            z := shr(1, add(z, div(x, z)))
+            // 6 Babylonian steps; z = (x/z + z) / 2
             z := shr(1, add(z, div(x, z)))
             z := shr(1, add(z, div(x, z)))
             z := shr(1, add(z, div(x, z)))
@@ -799,23 +794,25 @@ library FixedPointMathLib {
     /// @dev Returns the cube root of `x`, rounded down.
     /// Credit to bout3fiddy and pcaversaccio under AGPLv3 license:
     /// https://github.com/pcaversaccio/snekmate/blob/main/src/snekmate/utils/math.vy
-    /// Formally verified by xuwinnie:
-    /// https://github.com/vectorized/solady/blob/main/audits/xuwinnie-solady-cbrt-proof.pdf
     function cbrt(uint256 x) internal pure returns (uint256 z) {
         /// @solidity memory-safe-assembly
         assembly {
-            // Initial guess z = 2^ceil((log2(x) + 2) / 3).
-            // Since log2(x) = 255 - clz(x), the expression shl((257 - clz(x)) / 3, 1)
-            // computes this over-estimate. Guaranteed ≥ cbrt(x) and safe for Newton-Raphson's.
-            z := shl(div(sub(257, clz(x)), 3), 1)
-            // Newton-Raphson's.
-            z := div(add(add(div(x, mul(z, z)), z), z), 3)
+            // Initial guess z ≈ ∛(3/4) · 2𐞥 where q = ⌊(257 − clz(x)) / 3⌋. The
+            // multiplier 233/256 ≈ 0.910 ≈ ∛(3/4) balances the worst-case
+            // over/underestimate across each octave triplet (ε_over ≈ 0.445,
+            // ε_under ≈ −0.278), giving >85 bits of precision after 6 N-R
+            // iterations. The `lt(0, x)` term ensures z ≥ 1 when x > 0 (the
+            // `shr` can produce 0 for small `q`)
+            z := add(shr(8, shl(div(sub(257, clz(x)), 3), 233)), lt(0, x))
+
+            // 6 Newton-Raphson iterations.
             z := div(add(add(div(x, mul(z, z)), z), z), 3)
             z := div(add(add(div(x, mul(z, z)), z), z), 3)
             z := div(add(add(div(x, mul(z, z)), z), z), 3)
             z := div(add(add(div(x, mul(z, z)), z), z), 3)
             z := div(add(add(div(x, mul(z, z)), z), z), 3)
             z := div(add(add(div(x, mul(z, z)), z), z), 3)
+
             // Round down.
             z := sub(z, lt(div(x, mul(z, z)), z))
         }