From 2c0b81280e1bb373073acaf22ad4d5bd99ec8be0 Mon Sep 17 00:00:00 2001
From: Duncan Townsend <git@duncancmt.com>
Date: Sun, 22 Feb 2026 15:03:07 +0100
Subject: [PATCH 01/17] add alt 512-bit sqrt path and fuzz gas comparison

---
 src/utils/512Math.sol     | 77 +++++++++++++++++++++++++++++++++++++++
 test/0.8.25/512Math.t.sol | 42 +++++++++++++++++++++
 2 files changed, 119 insertions(+)

diff --git a/src/utils/512Math.sol b/src/utils/512Math.sol
index 883dcffef..711779038 100644
--- a/src/utils/512Math.sol
+++ b/src/utils/512Math.sol
@@ -1846,6 +1846,50 @@ library Lib512MathArithmetic {
         }
     }
 
+    // Alternative square root implementation using the Karatsuba-style split
+    // from SimonSuckut/Solidity_Uint512 (`Uint512.sqrt512`) adapted to this
+    // library's helpers and conventions.
+    function _sqrtAlt(uint256 x_hi, uint256 x_lo) private pure returns (uint256 r) {
+        unchecked {
+            // Normalize so the top limb has its MSB in bit 255 or 254.
+            uint256 shift = x_hi.clz() & 0xfe;
+            (, x_hi, x_lo) = _shl256(x_hi, x_lo, shift);
+
+            // Split √x into high/low 128-bit limbs:
+            //   x = [x_hi x_lo], r = (sp << 128) + q.
+            uint256 sp = x_hi.sqrt();
+            uint256 rp = x_hi - sp * sp;
+
+            uint256 nom;
+            uint256 denom;
+            uint256 u;
+            uint256 q;
+            assembly ("memory-safe") {
+                nom := or(shl(0x80, rp), shr(0x80, x_lo))
+                denom := shl(0x01, sp)
+                q := div(nom, denom)
+                u := mod(nom, denom)
+
+                // `nom` can be 257 bits. Handle the carry bit separately to
+                // avoid using 512/256 division.
+                let carry := shr(0x80, rp)
+                let x := mul(carry, not(0x00))
+                q := add(q, div(x, denom))
+                u := add(u, add(carry, mod(x, denom)))
+                q := add(q, div(u, denom))
+                u := mod(u, denom)
+            }
+
+            r = (sp << 128) + q;
+
+            uint256 rl = (u << 128) | (x_lo & 0xffffffffffffffffffffffffffffffff);
+            uint256 rr = q * q;
+            r = r.unsafeDec((q >> 128) > (u >> 128) || (((q >> 128) == (u >> 128)) && rl < rr));
+
+            return r >> (shift >> 1);
+        }
+    }
+
     function sqrt(uint512 x) internal pure returns (uint256) {
         (uint256 x_hi, uint256 x_lo) = x.into();
 
@@ -1862,6 +1906,20 @@ library Lib512MathArithmetic {
         return r.unsafeDec(_gt(r2_hi, r2_lo, x_hi, x_lo));
     }
 
+    function sqrtAlt(uint512 x) internal pure returns (uint256) {
+        (uint256 x_hi, uint256 x_lo) = x.into();
+
+        if (x_hi == 0) {
+            return x_lo.sqrt();
+        }
+
+        uint256 r = _sqrtAlt(x_hi, x_lo);
+
+        // Clamp to floor if the Karatsuba estimate overshot by 1.
+        (uint256 r2_hi, uint256 r2_lo) = _mul(r, r);
+        return r.unsafeDec(_gt(r2_hi, r2_lo, x_hi, x_lo));
+    }
+
     function osqrtUp(uint512 r, uint512 x) internal pure returns (uint512) {
         (uint256 x_hi, uint256 x_lo) = x.into();
 
@@ -1879,10 +1937,29 @@ library Lib512MathArithmetic {
         return r.from(r_hi, r_lo);
     }
 
+    function osqrtUpAlt(uint512 r, uint512 x) internal pure returns (uint512) {
+        (uint256 x_hi, uint256 x_lo) = x.into();
+
+        if (x_hi == 0) {
+            return r.from(0, x_lo.sqrtUp());
+        }
+
+        uint256 r_lo = _sqrtAlt(x_hi, x_lo);
+
+        (uint256 r2_hi, uint256 r2_lo) = _mul(r_lo, r_lo);
+        uint256 r_hi;
+        (r_hi, r_lo) = _add(0, r_lo, _gt(x_hi, x_lo, r2_hi, r2_lo).toUint());
+        return r.from(r_hi, r_lo);
+    }
+
     function isqrtUp(uint512 r) internal pure returns (uint512) {
         return osqrtUp(r, r);
     }
 
+    function isqrtUpAlt(uint512 r) internal pure returns (uint512) {
+        return osqrtUpAlt(r, r);
+    }
+
     function oshr(uint512 r, uint512 x, uint256 s) internal pure returns (uint512) {
         (uint256 x_hi, uint256 x_lo) = x.into();
         (uint256 r_hi, uint256 r_lo) = _shr(x_hi, x_lo, s);
diff --git a/test/0.8.25/512Math.t.sol b/test/0.8.25/512Math.t.sol
index d202c4518..d3f8a6ea8 100644
--- a/test/0.8.25/512Math.t.sol
+++ b/test/0.8.25/512Math.t.sol
@@ -239,6 +239,26 @@ contract Lib512MathTest is Test {
         }
     }
 
+    function test512Math_sqrtAlt(uint256 x_hi, uint256 x_lo) external pure {
+        uint512 x = alloc().from(x_hi, x_lo);
+        uint256 r = x.sqrtAlt();
+
+        (uint256 r2_lo, uint256 r2_hi) = SlowMath.fullMul(r, r);
+        assertTrue((r2_hi < x_hi) || (r2_hi == x_hi && r2_lo <= x_lo), "sqrtAlt too high");
+
+        if (r == type(uint256).max) {
+            assertTrue(
+                x_hi > 0xfffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffe
+                    || (x_hi == 0xfffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffe && x_lo != 0),
+                "sqrtAlt too low (overflow)"
+            );
+        } else {
+            r++;
+            (r2_lo, r2_hi) = SlowMath.fullMul(r, r);
+            assertTrue((r2_hi > x_hi) || (r2_hi == x_hi && r2_lo > x_lo), "sqrtAlt too low");
+        }
+    }
+
     function test512Math_divUpAlt(uint256 x_hi, uint256 x_lo, uint256 y_hi, uint256 y_lo) external view {
         vm.assume(y_hi != 0);
 
@@ -422,6 +442,28 @@ contract Lib512MathTest is Test {
         }
     }
 
+    function test512Math_osqrtUpAlt(uint256 x_hi, uint256 x_lo) external pure {
+        uint512 x = alloc().from(x_hi, x_lo);
+        (uint256 r_hi, uint256 r_lo) = alloc().osqrtUpAlt(x).into();
+
+        if (r_hi == 0 && r_lo == 0) {
+            assertTrue(x_hi == 0 && x_lo == 0, "sqrtUpAlt of nonzero is zero");
+        } else if (r_hi != 0) {
+            assertTrue(r_hi == 1 && r_lo == 0, "overflow result must be exactly 2^256");
+            (uint256 r_dec2_lo, uint256 r_dec2_hi) = SlowMath.fullMul(type(uint256).max, type(uint256).max);
+            assertTrue((r_dec2_hi < x_hi) || (r_dec2_hi == x_hi && r_dec2_lo < x_lo), "sqrtUpAlt too high");
+        } else {
+            (uint256 r2_lo, uint256 r2_hi) = SlowMath.fullMul(r_lo, r_lo);
+            assertTrue((r2_hi > x_hi) || (r2_hi == x_hi && r2_lo >= x_lo), "sqrtUpAlt too low");
+
+            if (r_lo != 1) {
+                uint256 r_dec_lo = r_lo - 1;
+                (r2_lo, r2_hi) = SlowMath.fullMul(r_dec_lo, r_dec_lo);
+                assertTrue((r2_hi < x_hi) || (r2_hi == x_hi && r2_lo < x_lo), "sqrtUpAlt too high");
+            }
+        }
+    }
+
     function test512Math_oshrUp(uint256 x_hi, uint256 x_lo, uint256 s) external pure {
         s = bound(s, 0, 512);
 

From a6bd3ca37718257b6527ddcf0dce3dff8275c465 Mon Sep 17 00:00:00 2001
From: Duncan Townsend <git@duncancmt.com>
Date: Sun, 22 Feb 2026 15:39:40 +0100
Subject: [PATCH 02/17] Comments, and a bit of gas golfing

---
 src/utils/512Math.sol | 67 ++++++++++++++++++++++---------------------
 1 file changed, 35 insertions(+), 32 deletions(-)

diff --git a/src/utils/512Math.sol b/src/utils/512Math.sol
index 711779038..37e152793 100644
--- a/src/utils/512Math.sol
+++ b/src/utils/512Math.sol
@@ -1846,46 +1846,49 @@ library Lib512MathArithmetic {
         }
     }
 
-    // Alternative square root implementation using the Karatsuba-style split
-    // from SimonSuckut/Solidity_Uint512 (`Uint512.sqrt512`) adapted to this
-    // library's helpers and conventions.
     function _sqrtAlt(uint256 x_hi, uint256 x_lo) private pure returns (uint256 r) {
         unchecked {
-            // Normalize so the top limb has its MSB in bit 255 or 254.
+            /// Our general approach is to apply Zimmerman's "Karatsuba Square Root" algorithm
+            /// https://inria.hal.science/inria-00072854/document with the helpers from Solady and
+            /// 512Math. This approach is inspired by
+            /// https://github.com/SimonSuckut/Solidity_Uint512/
+
+            // Normalize `x` so the top limb has its MSB in bit 255 or 254.
+            //   x ≥ 2²⁵³
             uint256 shift = x_hi.clz() & 0xfe;
             (, x_hi, x_lo) = _shl256(x_hi, x_lo, shift);
 
-            // Split √x into high/low 128-bit limbs:
-            //   x = [x_hi x_lo], r = (sp << 128) + q.
-            uint256 sp = x_hi.sqrt();
-            uint256 rp = x_hi - sp * sp;
+            // We treat `r` as a ≤2-limb bigint where each limb is half a machine word (128 bits).
+            // Spliting √x in this way lets us apply the ordinary 256-bit `sqrt` to the top
+            // limb. Then we can recover the bottom limb without 512-bit division.
+            uint256 r_hi = x_hi.sqrt();
+            uint256 res = x_hi - r_hi * r_hi;
 
-            uint256 nom;
-            uint256 denom;
-            uint256 u;
-            uint256 q;
+            uint256 r_lo;
             assembly ("memory-safe") {
-                nom := or(shl(0x80, rp), shr(0x80, x_lo))
-                denom := shl(0x01, sp)
-                q := div(nom, denom)
-                u := mod(nom, denom)
-
-                // `nom` can be 257 bits. Handle the carry bit separately to
-                // avoid using 512/256 division.
-                let carry := shr(0x80, rp)
-                let x := mul(carry, not(0x00))
-                q := add(q, div(x, denom))
-                u := add(u, add(carry, mod(x, denom)))
-                q := add(q, div(u, denom))
-                u := mod(u, denom)
+                let n := or(shl(0x80, res), shr(0x80, x_lo))
+                let d := shl(0x01, r_hi)
+                r_lo := div(n, d)
+
+                // It's possible that `n` was 257 bits and overflowed. Explicitly handling the carry
+                // avoids 512-bit division.
+                let c := shr(0x80, res)
+                let c_ := mul(c, not(0x00))
+                res := mod(n, d)
+                r_lo := add(r_lo, div(c_, d))
+                res := add(res, add(c, mod(c_, d)))
+                r_lo := add(r_lo, div(res, d))
+                res := mod(res, d)
             }
-
-            r = (sp << 128) + q;
-
-            uint256 rl = (u << 128) | (x_lo & 0xffffffffffffffffffffffffffffffff);
-            uint256 rr = q * q;
-            r = r.unsafeDec((q >> 128) > (u >> 128) || (((q >> 128) == (u >> 128)) && rl < rr));
-
+            r = (r_hi << 128) + r_lo;
+
+            r = r.unsafeDec(
+                ((r_lo >> 128) > (res >> 128))
+                .or(
+                    ((r_lo >> 128) == (res >> 128))
+                    .and((res << 128) | (x_lo & 0xffffffffffffffffffffffffffffffff) < r_lo * r_lo)
+                )
+            );
             return r >> (shift >> 1);
         }
     }

From 6eacdc7467e02d0fe6f8227f8afd6deb708e625c Mon Sep 17 00:00:00 2001
From: Duncan Townsend <git@duncancmt.com>
Date: Sun, 22 Feb 2026 15:58:47 +0100
Subject: [PATCH 03/17] Golf

---
 src/utils/512Math.sol | 6 +-----
 1 file changed, 1 insertion(+), 5 deletions(-)

diff --git a/src/utils/512Math.sol b/src/utils/512Math.sol
index 37e152793..363478f8b 100644
--- a/src/utils/512Math.sol
+++ b/src/utils/512Math.sol
@@ -1916,11 +1916,7 @@ library Lib512MathArithmetic {
             return x_lo.sqrt();
         }
 
-        uint256 r = _sqrtAlt(x_hi, x_lo);
-
-        // Clamp to floor if the Karatsuba estimate overshot by 1.
-        (uint256 r2_hi, uint256 r2_lo) = _mul(r, r);
-        return r.unsafeDec(_gt(r2_hi, r2_lo, x_hi, x_lo));
+        return _sqrtAlt(x_hi, x_lo);
     }
 
     function osqrtUp(uint512 r, uint512 x) internal pure returns (uint512) {

From 2665d20a6f5c9fcde9e89398f4144c14eaa1ab53 Mon Sep 17 00:00:00 2001
From: Duncan Townsend <git@duncancmt.com>
Date: Sun, 22 Feb 2026 16:09:48 +0100
Subject: [PATCH 04/17] Pedantry

---
 src/utils/512Math.sol | 5 ++++-
 1 file changed, 4 insertions(+), 1 deletion(-)

diff --git a/src/utils/512Math.sol b/src/utils/512Math.sol
index 363478f8b..56ddb1081 100644
--- a/src/utils/512Math.sol
+++ b/src/utils/512Math.sol
@@ -1873,7 +1873,7 @@ library Lib512MathArithmetic {
                 // It's possible that `n` was 257 bits and overflowed. Explicitly handling the carry
                 // avoids 512-bit division.
                 let c := shr(0x80, res)
-                let c_ := mul(c, not(0x00))
+                let c_ := mul(not(0x00), c)
                 res := mod(n, d)
                 r_lo := add(r_lo, div(c_, d))
                 res := add(res, add(c, mod(c_, d)))
@@ -1882,6 +1882,7 @@ library Lib512MathArithmetic {
             }
             r = (r_hi << 128) + r_lo;
 
+            // Handle negative residue
             r = r.unsafeDec(
                 ((r_lo >> 128) > (res >> 128))
                 .or(
@@ -1889,6 +1890,8 @@ library Lib512MathArithmetic {
                     .and((res << 128) | (x_lo & 0xffffffffffffffffffffffffffffffff) < r_lo * r_lo)
                 )
             );
+
+            // Un-normalize
             return r >> (shift >> 1);
         }
     }

From d0936be3c79a940082d9f64380d9c2bdc2e50eb2 Mon Sep 17 00:00:00 2001
From: Duncan Townsend <git@duncancmt.com>
Date: Sun, 22 Feb 2026 16:14:33 +0100
Subject: [PATCH 05/17] Comment

---
 src/utils/512Math.sol | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/utils/512Math.sol b/src/utils/512Math.sol
index 56ddb1081..248f1aaa8 100644
--- a/src/utils/512Math.sol
+++ b/src/utils/512Math.sol
@@ -1859,8 +1859,8 @@ library Lib512MathArithmetic {
             (, x_hi, x_lo) = _shl256(x_hi, x_lo, shift);
 
             // We treat `r` as a ≤2-limb bigint where each limb is half a machine word (128 bits).
-            // Spliting √x in this way lets us apply the ordinary 256-bit `sqrt` to the top
-            // limb. Then we can recover the bottom limb without 512-bit division.
+            // Spliting √x in this way lets us apply the ordinary 256-bit `sqrt` to the top limb of
+            // `x`. Then we can recover the bottom limb or `r` without 512-bit division.
             uint256 r_hi = x_hi.sqrt();
             uint256 res = x_hi - r_hi * r_hi;
 

From d3dff10eb373ed14336014b6e4190724d85a9345 Mon Sep 17 00:00:00 2001
From: Duncan Townsend <git@duncancmt.com>
Date: Sun, 22 Feb 2026 16:33:04 +0100
Subject: [PATCH 06/17] Golf

---
 src/utils/512Math.sol | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/src/utils/512Math.sol b/src/utils/512Math.sol
index 248f1aaa8..cf4836c32 100644
--- a/src/utils/512Math.sol
+++ b/src/utils/512Math.sol
@@ -1873,10 +1873,10 @@ library Lib512MathArithmetic {
                 // It's possible that `n` was 257 bits and overflowed. Explicitly handling the carry
                 // avoids 512-bit division.
                 let c := shr(0x80, res)
-                let c_ := mul(not(0x00), c)
+                let neg_c := sub(0x00, c)
                 res := mod(n, d)
-                r_lo := add(r_lo, div(c_, d))
-                res := add(res, add(c, mod(c_, d)))
+                r_lo := add(r_lo, div(neg_c, d))
+                res := add(res, add(c, mod(neg_c, d)))
                 r_lo := add(r_lo, div(res, d))
                 res := mod(res, d)
             }

From e5d43470d01923b36d91b56bae377e6792cc8816 Mon Sep 17 00:00:00 2001
From: Duncan Townsend <git@duncancmt.com>
Date: Sun, 22 Feb 2026 16:45:49 +0100
Subject: [PATCH 07/17] Golf

---
 src/utils/512Math.sol | 34 +++++++++++++++++++++++++++++++---
 1 file changed, 31 insertions(+), 3 deletions(-)

diff --git a/src/utils/512Math.sol b/src/utils/512Math.sol
index cf4836c32..4ed405faf 100644
--- a/src/utils/512Math.sol
+++ b/src/utils/512Math.sol
@@ -1859,12 +1859,40 @@ library Lib512MathArithmetic {
             (, x_hi, x_lo) = _shl256(x_hi, x_lo, shift);
 
             // We treat `r` as a ≤2-limb bigint where each limb is half a machine word (128 bits).
-            // Spliting √x in this way lets us apply the ordinary 256-bit `sqrt` to the top limb of
+            // Spliting √x in this way lets us apply "ordinary" 256-bit `sqrt` to the top word of
             // `x`. Then we can recover the bottom limb or `r` without 512-bit division.
-            uint256 r_hi = x_hi.sqrt();
+            //
+            // Implementing this as:
+            //   uint256 r_hi = x_hi.sqrt();
+            // is correct, but duplicates the normalization that we just did above and performs a
+            // more-costly initialization step. solc is not smart enough to optimize this away, so
+            // we inline and do it ourselves.
+            uint256 r_hi;
+            assembly ("memory-safe") {
+                // Initialization requires that the first bit of `r_hi` be in the correct
+                // position. This is correct from the normalization above.
+                r_hi := 0x80000000000000000000000000000000
+
+                // Seven Babylonian steps is sufficient for convergence.
+                r_hi := shr(0x01, add(r_hi, div(x_hi, r_hi)))
+                r_hi := shr(0x01, add(r_hi, div(x_hi, r_hi)))
+                r_hi := shr(0x01, add(r_hi, div(x_hi, r_hi)))
+                r_hi := shr(0x01, add(r_hi, div(x_hi, r_hi)))
+                r_hi := shr(0x01, add(r_hi, div(x_hi, r_hi)))
+                r_hi := shr(0x01, add(r_hi, div(x_hi, r_hi)))
+                r_hi := shr(0x01, add(r_hi, div(x_hi, r_hi)))
+
+                // The Babylonian step can oscillate between ⌊√x_hi⌋ and ⌈√x_hi⌉. Clean that up.
+                r_hi := sub(r_hi, lt(div(x_hi, r_hi), r_hi))
+            }
             uint256 res = x_hi - r_hi * r_hi;
 
             uint256 r_lo;
+            // `res` is (almost) a single limb. Create a new machine word `w` with `res` as the
+            // upper limb and shifting in the next limb of `x` (namely `x_lo >> 128`) as the lower
+            // limb. The next step of Zimmerman's algorithm is:
+            //   r_lo = w / (2 · r_hi)
+            //   res = w % (2 · r_hi)
             assembly ("memory-safe") {
                 let n := or(shl(0x80, res), shr(0x80, x_lo))
                 let d := shl(0x01, r_hi)
@@ -1882,7 +1910,7 @@ library Lib512MathArithmetic {
             }
             r = (r_hi << 128) + r_lo;
 
-            // Handle negative residue
+            // Then, if res · 2¹²⁸ + x_lo % 2¹²⁸ < r_lo², decrement `r`
             r = r.unsafeDec(
                 ((r_lo >> 128) > (res >> 128))
                 .or(

From bf12bb86d47d5c6c53ab97f5e9b84da611d51b4f Mon Sep 17 00:00:00 2001
From: Duncan Townsend <git@duncancmt.com>
Date: Sun, 22 Feb 2026 16:50:35 +0100
Subject: [PATCH 08/17] Comments

---
 src/utils/512Math.sol | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/src/utils/512Math.sol b/src/utils/512Math.sol
index 4ed405faf..269e38795 100644
--- a/src/utils/512Math.sol
+++ b/src/utils/512Math.sol
@@ -1860,7 +1860,7 @@ library Lib512MathArithmetic {
 
             // We treat `r` as a ≤2-limb bigint where each limb is half a machine word (128 bits).
             // Spliting √x in this way lets us apply "ordinary" 256-bit `sqrt` to the top word of
-            // `x`. Then we can recover the bottom limb or `r` without 512-bit division.
+            // `x`. Then we can recover the bottom limb of `r` without 512-bit division.
             //
             // Implementing this as:
             //   uint256 r_hi = x_hi.sqrt();
@@ -1898,8 +1898,8 @@ library Lib512MathArithmetic {
                 let d := shl(0x01, r_hi)
                 r_lo := div(n, d)
 
-                // It's possible that `n` was 257 bits and overflowed. Explicitly handling the carry
-                // avoids 512-bit division.
+                // It's possible that `n` was 257 bits and overflowed (`res` was not just a single
+                // limb). Explicitly handling the carry avoids 512-bit division.
                 let c := shr(0x80, res)
                 let neg_c := sub(0x00, c)
                 res := mod(n, d)

From aa5b92dcc0327c4990d8f4d5b603ffd181f6a93e Mon Sep 17 00:00:00 2001
From: Duncan Townsend <git@duncancmt.com>
Date: Sun, 22 Feb 2026 16:54:04 +0100
Subject: [PATCH 09/17] Remove the old 512-bit `sqrt` functions, replace them
 with the new, optimized ones

---
 src/utils/512Math.sol     | 196 +-------------------------------------
 test/0.8.25/512Math.t.sol |  42 --------
 2 files changed, 1 insertion(+), 237 deletions(-)

diff --git a/src/utils/512Math.sol b/src/utils/512Math.sol
index 269e38795..b01063552 100644
--- a/src/utils/512Math.sol
+++ b/src/utils/512Math.sol
@@ -1689,164 +1689,7 @@ library Lib512MathArithmetic {
         return omodAlt(r, y, r);
     }
 
-    // hi ≈ x · y / 2²⁵⁶ (±1)
-    function _inaccurateMulHi(uint256 x, uint256 y) private pure returns (uint256 hi) {
-        assembly ("memory-safe") {
-            hi := sub(mulmod(x, y, not(0x00)), mul(x, y))
-        }
-    }
-
-    // gas benchmark 2025/09/20: ~1425 gas
     function _sqrt(uint256 x_hi, uint256 x_lo) private pure returns (uint256 r) {
-        /// Our general approach here is to compute the inverse of the square root of the argument
-        /// using Newton-Raphson iterations. Then we combine (multiply) this inverse square root
-        /// approximation with the argument to approximate the square root of the argument. After
-        /// that, a final fixup step is applied to get the exact result. We compute the inverse of
-        /// the square root rather than the square root directly because then our Newton-Raphson
-        /// iteration can avoid the extremely expensive 512-bit division subroutine.
-        unchecked {
-            /// First, we normalize `x` by separating it into a mantissa and exponent. We use
-            /// even-exponent normalization.
-
-            // `e` is half the exponent of `x`
-            // e    = ⌊bitlength(x)/2⌋
-            // invE = 256 - e
-            uint256 invE = (x_hi.clz() + 1) >> 1; // invE ∈ [0, 128]
-
-            // Extract mantissa M by shifting x right by 2·e - 255 bits
-            // `M` is the mantissa of `x` as a Q1.255; M ∈ [½, 2)
-            (, uint256 M) = _shr(x_hi, x_lo, 257 - (invE << 1)); // scale: 2⁽²⁵⁵⁻²ᵉ⁾
-
-            /// Pick an initial estimate (seed) for Y using a lookup table. Even-exponent
-            /// normalization means our mantissa is geometrically symmetric around 1, leading to 16
-            /// buckets on the low side and 32 buckets on the high side.
-            // `Y` _ultimately_ approximates the inverse square root of fixnum `M` as a
-            // Q3.253. However, as a gas optimization, the number of fractional bits in `Y` rises
-            // through the steps, giving an inhomogeneous fixed-point representation. Y ≈∈ [√½, √2]
-            uint256 Y; // scale: 2⁽²⁵³⁺ᵉ⁾
-            uint256 Mbucket;
-            assembly ("memory-safe") {
-                // Extract the upper 6 bits of `M` to be used as a table index. `M >> 250 < 16` is
-                // invalid (that would imply M<½), so our lookup table only needs to handle only 16
-                // through 63.
-                Mbucket := shr(0xfa, M)
-                // We can't fit 48 seeds into a single word, so we split the table in 2 and use `c`
-                // to select which table we index.
-                let c := lt(0x27, Mbucket)
-
-                // Each entry is 10 bits and the entries are ordered from lowest `i` to highest. The
-                // seed is the value for `Y` for the midpoint of the bucket, rounded to 10
-                // significant bits. That is, Y ≈ 1/√(2·M_mid), as a Q247.9. The 2 comes from the
-                // half-scale difference between Y and √M. The optimality of this choice was
-                // verified by fuzzing.
-                let table_hi := 0x71dc26f1b76c9ad6a5a46819c661946418c621856057e5ed775d1715b96b
-                let table_lo := 0xb26b4a8690a027198e559263e8ce2887e15832047f1f47b5e677dd974dcd
-                let table := xor(table_lo, mul(xor(table_hi, table_lo), c))
-
-                // Index the table to obtain the initial seed of `Y`.
-                let shift := add(0x186, mul(0x0a, sub(mul(0x18, c), Mbucket)))
-                // We begin the Newton-Raphson iterations with `Y` in Q247.9 format.
-                Y := and(0x3ff, shr(shift, table))
-
-                // The worst-case seed for `Y` occurs when `Mbucket = 16`. For monotone quadratic
-                // convergence, we desire that 1/√3 < Y·√M < √(5/3). At the boundaries (worst case)
-                // of the `Mbucket = 16` range, we are 0.407351 (41.3680%) from the lower bound and
-                // 0.275987 (27.1906%) from the higher bound.
-            }
-
-            /// Perform 5 Newton-Raphson iterations. 5 is enough iterations for sufficient
-            /// convergence that our final fixup step produces an exact result.
-            // The Newton-Raphson iteration for 1/√M is:
-            //     Y ≈ Y · (3 - M · Y²) / 2
-            // The implementation of this iteration is deliberately imprecise. No matter how many
-            // times you run it, you won't converge `Y` on the closest Q3.253 to √M. However, this
-            // is acceptable because the cleanup step applied after the final call is very tolerant
-            // of error in the low bits of `Y`.
-
-            // `M` is Q1.255
-            // `Y` is Q247.9
-            {
-                uint256 Y2 = Y * Y;                    // scale: 2¹⁸
-                // Because `M` is Q1.255, multiplying `Y2` by `M` and taking the high word
-                // implicitly divides `MY2` by 2. We move the division by 2 inside the subtraction
-                // from 3 by adjusting the minuend.
-                uint256 MY2 = _inaccurateMulHi(M, Y2); // scale: 2¹⁸
-                uint256 T = 1.5 * 2 ** 18 - MY2;       // scale: 2¹⁸
-                Y *= T;                                // scale: 2²⁷
-            }
-            // `Y` is Q229.27
-            {
-                uint256 Y2 = Y * Y;                    // scale: 2⁵⁴
-                uint256 MY2 = _inaccurateMulHi(M, Y2); // scale: 2⁵⁴
-                uint256 T = 1.5 * 2 ** 54 - MY2;       // scale: 2⁵⁴
-                Y *= T;                                // scale: 2⁸¹
-            }
-            // `Y` is Q175.81
-            {
-                uint256 Y2 = Y * Y;                    // scale: 2¹⁶²
-                uint256 MY2 = _inaccurateMulHi(M, Y2); // scale: 2¹⁶²
-                uint256 T = 1.5 * 2 ** 162 - MY2;      // scale: 2¹⁶²
-                Y = Y * T >> 116;                      // scale: 2¹²⁷
-            }
-            // `Y` is Q129.127
-            if (invE < 95 - Mbucket) {
-                // Generally speaking, for relatively smaller `e` (lower values of `x`) and for
-                // relatively larger `M`, we can skip the 5th N-R iteration. The constant `95` is
-                // derived by extensive fuzzing. Attempting a higher-order approximation of the
-                // relationship between `M` and `invE` consumes, on average, more gas. When this
-                // branch is not taken, the correct bits that this iteration would obtain are
-                // shifted away during the denormalization step. This branch is net gas-optimizing.
-                uint256 Y2 = Y * Y;                    // scale: 2²⁵⁴
-                uint256 MY2 = _inaccurateMulHi(M, Y2); // scale: 2²⁵⁴
-                uint256 T = 1.5 * 2 ** 254 - MY2;      // scale: 2²⁵⁴
-                Y = _inaccurateMulHi(Y << 2, T);       // scale: 2¹²⁷
-            }
-            // `Y` is Q129.127
-            {
-                uint256 Y2 = Y * Y;                    // scale: 2²⁵⁴
-                uint256 MY2 = _inaccurateMulHi(M, Y2); // scale: 2²⁵⁴
-                uint256 T = 1.5 * 2 ** 254 - MY2;      // scale: 2²⁵⁴
-                Y = _inaccurateMulHi(Y << 128, T);     // scale: 2²⁵³
-            }
-            // `Y` is Q3.253
-
-            /// When we combine `Y` with `M` to form our approximation of the square root, we have
-            /// to un-normalize by the half-scale value. This is where even-exponent normalization
-            /// comes in because the half-scale is integral.
-            ///     M   = ⌊x · 2⁽²⁵⁵⁻²ᵉ⁾⌋
-            ///     Y   ≈ 2²⁵³ / √(M / 2²⁵⁵)
-            ///     Y   ≈ 2³⁸¹ / √(2·M)
-            ///     M·Y ≈ 2³⁸¹ · √(M/2)
-            ///     M·Y ≈ 2⁽⁵⁰⁸⁻ᵉ⁾ · √x
-            ///     r0  ≈ M·Y / 2⁽⁵⁰⁸⁻ᵉ⁾ ≈ ⌊√x⌋
-            // We shift right by `508 - e` to account for both the Q3.253 scaling and
-            // denormalization. We don't care about accuracy in the low bits of `r0`, so we can cut
-            // some corners.
-            (, uint256 r0) = _shr(_inaccurateMulHi(M, Y), 0, 252 + invE);
-
-            /// `r0` is only an approximation of √x, so we perform a single Babylonian step to fully
-            /// converge on ⌊√x⌋ or ⌈√x⌉.  The Babylonian step is:
-            ///     r = ⌊(r0 + ⌊x/r0⌋) / 2⌋
-            // Rather than use the more-expensive division routine that returns a 512-bit result,
-            // because the value the upper word of the quotient can take is highly constrained, we
-            // can compute the quotient mod 2²⁵⁶ and recover the high word separately. Although
-            // `_div` does an expensive Newton-Raphson-Hensel modular inversion:
-            //     ⌊x/r0⌋ ≡ ⌊x/2ⁿ⌋·⌊r0/2ⁿ⌋⁻¹ mod 2²⁵⁶ (for r % 2⁽ⁿ⁺¹⁾ = 2ⁿ)
-            // and we already have a pretty good estimate for r0⁻¹, namely `Y`, refining `Y` into
-            // the appropriate inverse requires a series of 768-bit multiplications that take more
-            // gas.
-            uint256 q_lo = _div(x_hi, x_lo, r0);
-            uint256 q_hi = (r0 <= x_hi).toUint();
-            (uint256 s_hi, uint256 s_lo) = _add(q_hi, q_lo, r0);
-            // `oflo` here is either 0 or 1. When `oflo == 1`, `r == 0`, and the correct value for
-            // `r` is `type(uint256).max`.
-            uint256 oflo;
-            (oflo, r) = _shr256(s_hi, s_lo, 1);
-            r -= oflo; // underflow is desired
-        }
-    }
-
-    function _sqrtAlt(uint256 x_hi, uint256 x_lo) private pure returns (uint256 r) {
         unchecked {
             /// Our general approach is to apply Zimmerman's "Karatsuba Square Root" algorithm
             /// https://inria.hal.science/inria-00072854/document with the helpers from Solady and
@@ -1931,23 +1774,7 @@ library Lib512MathArithmetic {
             return x_lo.sqrt();
         }
 
-        uint256 r = _sqrt(x_hi, x_lo);
-
-        // Because the Babylonian step can give ⌈√x⌉ if x+1 is a perfect square, we have to
-        // check whether we've overstepped by 1 and clamp as appropriate. ref:
-        // https://en.wikipedia.org/wiki/Integer_square_root#Using_only_integer_division
-        (uint256 r2_hi, uint256 r2_lo) = _mul(r, r);
-        return r.unsafeDec(_gt(r2_hi, r2_lo, x_hi, x_lo));
-    }
-
-    function sqrtAlt(uint512 x) internal pure returns (uint256) {
-        (uint256 x_hi, uint256 x_lo) = x.into();
-
-        if (x_hi == 0) {
-            return x_lo.sqrt();
-        }
-
-        return _sqrtAlt(x_hi, x_lo);
+        return _sqrt(x_hi, x_lo);
     }
 
     function osqrtUp(uint512 r, uint512 x) internal pure returns (uint512) {
@@ -1959,23 +1786,6 @@ library Lib512MathArithmetic {
 
         uint256 r_lo = _sqrt(x_hi, x_lo);
 
-        // The Babylonian step can give ⌈√x⌉ if x+1 is a perfect square. This is
-        // fine. If the Babylonian step gave ⌊√x⌋ ≠ √x, we have to round up.
-        (uint256 r2_hi, uint256 r2_lo) = _mul(r_lo, r_lo);
-        uint256 r_hi;
-        (r_hi, r_lo) = _add(0, r_lo, _gt(x_hi, x_lo, r2_hi, r2_lo).toUint());
-        return r.from(r_hi, r_lo);
-    }
-
-    function osqrtUpAlt(uint512 r, uint512 x) internal pure returns (uint512) {
-        (uint256 x_hi, uint256 x_lo) = x.into();
-
-        if (x_hi == 0) {
-            return r.from(0, x_lo.sqrtUp());
-        }
-
-        uint256 r_lo = _sqrtAlt(x_hi, x_lo);
-
         (uint256 r2_hi, uint256 r2_lo) = _mul(r_lo, r_lo);
         uint256 r_hi;
         (r_hi, r_lo) = _add(0, r_lo, _gt(x_hi, x_lo, r2_hi, r2_lo).toUint());
@@ -1986,10 +1796,6 @@ library Lib512MathArithmetic {
         return osqrtUp(r, r);
     }
 
-    function isqrtUpAlt(uint512 r) internal pure returns (uint512) {
-        return osqrtUpAlt(r, r);
-    }
-
     function oshr(uint512 r, uint512 x, uint256 s) internal pure returns (uint512) {
         (uint256 x_hi, uint256 x_lo) = x.into();
         (uint256 r_hi, uint256 r_lo) = _shr(x_hi, x_lo, s);
diff --git a/test/0.8.25/512Math.t.sol b/test/0.8.25/512Math.t.sol
index d3f8a6ea8..d202c4518 100644
--- a/test/0.8.25/512Math.t.sol
+++ b/test/0.8.25/512Math.t.sol
@@ -239,26 +239,6 @@ contract Lib512MathTest is Test {
         }
     }
 
-    function test512Math_sqrtAlt(uint256 x_hi, uint256 x_lo) external pure {
-        uint512 x = alloc().from(x_hi, x_lo);
-        uint256 r = x.sqrtAlt();
-
-        (uint256 r2_lo, uint256 r2_hi) = SlowMath.fullMul(r, r);
-        assertTrue((r2_hi < x_hi) || (r2_hi == x_hi && r2_lo <= x_lo), "sqrtAlt too high");
-
-        if (r == type(uint256).max) {
-            assertTrue(
-                x_hi > 0xfffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffe
-                    || (x_hi == 0xfffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffe && x_lo != 0),
-                "sqrtAlt too low (overflow)"
-            );
-        } else {
-            r++;
-            (r2_lo, r2_hi) = SlowMath.fullMul(r, r);
-            assertTrue((r2_hi > x_hi) || (r2_hi == x_hi && r2_lo > x_lo), "sqrtAlt too low");
-        }
-    }
-
     function test512Math_divUpAlt(uint256 x_hi, uint256 x_lo, uint256 y_hi, uint256 y_lo) external view {
         vm.assume(y_hi != 0);
 
@@ -442,28 +422,6 @@ contract Lib512MathTest is Test {
         }
     }
 
-    function test512Math_osqrtUpAlt(uint256 x_hi, uint256 x_lo) external pure {
-        uint512 x = alloc().from(x_hi, x_lo);
-        (uint256 r_hi, uint256 r_lo) = alloc().osqrtUpAlt(x).into();
-
-        if (r_hi == 0 && r_lo == 0) {
-            assertTrue(x_hi == 0 && x_lo == 0, "sqrtUpAlt of nonzero is zero");
-        } else if (r_hi != 0) {
-            assertTrue(r_hi == 1 && r_lo == 0, "overflow result must be exactly 2^256");
-            (uint256 r_dec2_lo, uint256 r_dec2_hi) = SlowMath.fullMul(type(uint256).max, type(uint256).max);
-            assertTrue((r_dec2_hi < x_hi) || (r_dec2_hi == x_hi && r_dec2_lo < x_lo), "sqrtUpAlt too high");
-        } else {
-            (uint256 r2_lo, uint256 r2_hi) = SlowMath.fullMul(r_lo, r_lo);
-            assertTrue((r2_hi > x_hi) || (r2_hi == x_hi && r2_lo >= x_lo), "sqrtUpAlt too low");
-
-            if (r_lo != 1) {
-                uint256 r_dec_lo = r_lo - 1;
-                (r2_lo, r2_hi) = SlowMath.fullMul(r_dec_lo, r_dec_lo);
-                assertTrue((r2_hi < x_hi) || (r2_hi == x_hi && r2_lo < x_lo), "sqrtUpAlt too high");
-            }
-        }
-    }
-
     function test512Math_oshrUp(uint256 x_hi, uint256 x_lo, uint256 s) external pure {
         s = bound(s, 0, 512);
 

From 2a8753fcbe7dd2e92a031f46b590dc327338d655 Mon Sep 17 00:00:00 2001
From: Duncan Townsend <git@duncancmt.com>
Date: Mon, 23 Feb 2026 12:37:10 +0100
Subject: [PATCH 10/17] Typo

---
 src/utils/512Math.sol | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/utils/512Math.sol b/src/utils/512Math.sol
index b01063552..cf948f997 100644
--- a/src/utils/512Math.sol
+++ b/src/utils/512Math.sol
@@ -1696,8 +1696,8 @@ library Lib512MathArithmetic {
             /// 512Math. This approach is inspired by
             /// https://github.com/SimonSuckut/Solidity_Uint512/
 
-            // Normalize `x` so the top limb has its MSB in bit 255 or 254.
-            //   x ≥ 2²⁵³
+            // Normalize `x` so the top word has its MSB in bit 255 or 254.
+            //   x ≥ 2⁵¹⁰
             uint256 shift = x_hi.clz() & 0xfe;
             (, x_hi, x_lo) = _shl256(x_hi, x_lo, shift);
 

From d85bc80a2b08d423de979ef8d6730445e75cea8b Mon Sep 17 00:00:00 2001
From: Duncan Townsend <git@duncancmt.com>
Date: Mon, 23 Feb 2026 16:55:13 +0100
Subject: [PATCH 11/17] Comment

---
 src/utils/512Math.sol | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/src/utils/512Math.sol b/src/utils/512Math.sol
index cf948f997..acd012727 100644
--- a/src/utils/512Math.sol
+++ b/src/utils/512Math.sol
@@ -1753,7 +1753,9 @@ library Lib512MathArithmetic {
             }
             r = (r_hi << 128) + r_lo;
 
-            // Then, if res · 2¹²⁸ + x_lo % 2¹²⁸ < r_lo², decrement `r`
+            // Then, if res · 2¹²⁸ + x_lo % 2¹²⁸ < r_lo², decrement `r`. We have to do this in a
+            // complicated manner because both `res` and `r_lo` can be _slightly_ longer than 1 limb
+            // (128 bits). This is more efficient than performing the full 257-bit comparison.
             r = r.unsafeDec(
                 ((r_lo >> 128) > (res >> 128))
                 .or(

From cb4bd31850d7938640192000a52bea6863fc128c Mon Sep 17 00:00:00 2001
From: Duncan Townsend <git@duncancmt.com>
Date: Mon, 23 Feb 2026 18:08:33 +0100
Subject: [PATCH 12/17] Homogenize

---
 src/utils/512Math.sol | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/utils/512Math.sol b/src/utils/512Math.sol
index acd012727..cfa758167 100644
--- a/src/utils/512Math.sol
+++ b/src/utils/512Math.sol
@@ -1757,9 +1757,9 @@ library Lib512MathArithmetic {
             // complicated manner because both `res` and `r_lo` can be _slightly_ longer than 1 limb
             // (128 bits). This is more efficient than performing the full 257-bit comparison.
             r = r.unsafeDec(
-                ((r_lo >> 128) > (res >> 128))
+                ((res >> 128) < (r_lo >> 128))
                 .or(
-                    ((r_lo >> 128) == (res >> 128))
+                    ((res >> 128) == (r_lo >> 128))
                     .and((res << 128) | (x_lo & 0xffffffffffffffffffffffffffffffff) < r_lo * r_lo)
                 )
             );

From 46c567e54f19fb77cf9e72a81168aa8161ca4fd3 Mon Sep 17 00:00:00 2001
From: Duncan Townsend <git@duncancmt.com>
Date: Mon, 23 Feb 2026 19:12:42 +0100
Subject: [PATCH 13/17] Comment

---
 src/utils/512Math.sol | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/src/utils/512Math.sol b/src/utils/512Math.sol
index cfa758167..ebb500531 100644
--- a/src/utils/512Math.sol
+++ b/src/utils/512Math.sol
@@ -1731,11 +1731,11 @@ library Lib512MathArithmetic {
             uint256 res = x_hi - r_hi * r_hi;
 
             uint256 r_lo;
-            // `res` is (almost) a single limb. Create a new machine word `w` with `res` as the
-            // upper limb and shifting in the next limb of `x` (namely `x_lo >> 128`) as the lower
-            // limb. The next step of Zimmerman's algorithm is:
-            //   r_lo = w / (2 · r_hi)
-            //   res = w % (2 · r_hi)
+            // `res` is (almost) a single limb. Create a new (almost) machine word `n` with `res` as
+            // the upper limb and shifting in the next limb of `x` (namely `x_lo >> 128`) as the
+            // lower limb. The next step of Zimmerman's algorithm is:
+            //   r_lo = n / (2 · r_hi)
+            //   res = n % (2 · r_hi)
             assembly ("memory-safe") {
                 let n := or(shl(0x80, res), shr(0x80, x_lo))
                 let d := shl(0x01, r_hi)

From 678f5a5de2d8b1426699463193c3e68dac6ebf6a Mon Sep 17 00:00:00 2001
From: Duncan Townsend <git@duncancmt.com>
Date: Mon, 23 Feb 2026 20:56:06 +0100
Subject: [PATCH 14/17] Golf

---
 src/utils/512Math.sol | 17 ++++++++++-------
 1 file changed, 10 insertions(+), 7 deletions(-)

diff --git a/src/utils/512Math.sol b/src/utils/512Math.sol
index ebb500531..77fa1fa9e 100644
--- a/src/utils/512Math.sol
+++ b/src/utils/512Math.sol
@@ -1741,15 +1741,18 @@ library Lib512MathArithmetic {
                 let d := shl(0x01, r_hi)
                 r_lo := div(n, d)
 
-                // It's possible that `n` was 257 bits and overflowed (`res` was not just a single
-                // limb). Explicitly handling the carry avoids 512-bit division.
                 let c := shr(0x80, res)
-                let neg_c := sub(0x00, c)
                 res := mod(n, d)
-                r_lo := add(r_lo, div(neg_c, d))
-                res := add(res, add(c, mod(neg_c, d)))
-                r_lo := add(r_lo, div(res, d))
-                res := mod(res, d)
+
+                // It's possible that `n` was 257 bits and overflowed (`res` was not just a single
+                // limb). Explicitly handling the carry avoids 512-bit division.
+                if c {
+                    let neg_c := not(0x00)
+                    r_lo := add(r_lo, div(neg_c, d))
+                    res := add(res, add(0x01, mod(neg_c, d)))
+                    r_lo := add(r_lo, div(res, d))
+                    res := mod(res, d)
+                }
             }
             r = (r_hi << 128) + r_lo;
 

From 22904c1527e19b8142883ab4a09abfac799bcd95 Mon Sep 17 00:00:00 2001
From: Duncan Townsend <git@duncancmt.com>
Date: Mon, 23 Feb 2026 21:02:33 +0100
Subject: [PATCH 15/17] Golf

---
 src/utils/512Math.sol | 7 ++++---
 1 file changed, 4 insertions(+), 3 deletions(-)

diff --git a/src/utils/512Math.sol b/src/utils/512Math.sol
index 77fa1fa9e..191be1d73 100644
--- a/src/utils/512Math.sol
+++ b/src/utils/512Math.sol
@@ -1698,8 +1698,9 @@ library Lib512MathArithmetic {
 
             // Normalize `x` so the top word has its MSB in bit 255 or 254.
             //   x ≥ 2⁵¹⁰
-            uint256 shift = x_hi.clz() & 0xfe;
-            (, x_hi, x_lo) = _shl256(x_hi, x_lo, shift);
+            uint256 shift = x_hi.clz();
+            (, x_hi, x_lo) = _shl256(x_hi, x_lo, shift & 0xfe);
+            shift >>= 1;
 
             // We treat `r` as a ≤2-limb bigint where each limb is half a machine word (128 bits).
             // Spliting √x in this way lets us apply "ordinary" 256-bit `sqrt` to the top word of
@@ -1768,7 +1769,7 @@ library Lib512MathArithmetic {
             );
 
             // Un-normalize
-            return r >> (shift >> 1);
+            return r >> shift;
         }
     }
 

From 041cf8e3eeb5660e651218ddf85b1c6a91a5223e Mon Sep 17 00:00:00 2001
From: Duncan Townsend <git@duncancmt.com>
Date: Mon, 23 Feb 2026 21:08:10 +0100
Subject: [PATCH 16/17] Golf

---
 src/utils/512Math.sol | 9 ++++++++-
 1 file changed, 8 insertions(+), 1 deletion(-)

diff --git a/src/utils/512Math.sol b/src/utils/512Math.sol
index 191be1d73..c63961635 100644
--- a/src/utils/512Math.sol
+++ b/src/utils/512Math.sol
@@ -1729,7 +1729,14 @@ library Lib512MathArithmetic {
                 // The Babylonian step can oscillate between ⌊√x_hi⌋ and ⌈√x_hi⌉. Clean that up.
                 r_hi := sub(r_hi, lt(div(x_hi, r_hi), r_hi))
             }
-            uint256 res = x_hi - r_hi * r_hi;
+
+            // This is cheaper than
+            //   uint256 res = x_hi - r_hi * r_hi;
+            // for no clear reason
+            uint256 res;
+            assembly ("memory-safe") {
+                res := sub(x_hi, mul(r_hi, r_hi))
+            }
 
             uint256 r_lo;
             // `res` is (almost) a single limb. Create a new (almost) machine word `n` with `res` as

From ab9c4b1e36da45fdf94de7a5ba07046f7098e3fa Mon Sep 17 00:00:00 2001
From: Duncan Townsend <git@duncancmt.com>
Date: Tue, 24 Feb 2026 12:15:23 +0100
Subject: [PATCH 17/17] Add a fuzz test to explicitly exercise a known
 weak-point of some `sqrt` implementations

---
 test/0.8.25/512Math.t.sol | 12 ++++++++++++
 1 file changed, 12 insertions(+)

diff --git a/test/0.8.25/512Math.t.sol b/test/0.8.25/512Math.t.sol
index d202c4518..727817109 100644
--- a/test/0.8.25/512Math.t.sol
+++ b/test/0.8.25/512Math.t.sol
@@ -422,6 +422,18 @@ contract Lib512MathTest is Test {
         }
     }
 
+    function test512Math_sqrt_perfectSquare(uint256 r) external pure {
+        uint512 x = alloc().omul(r, r);
+        assertEq(x.sqrt(), r);
+    }
+
+    function test512Math_osqrtUp_perfectSquare(uint256 r) external pure {
+        uint512 x = alloc().omul(r, r);
+        (uint256 r_hi, uint256 r_lo) = alloc().osqrtUp(x).into();
+        assertEq(r_hi, 0);
+        assertEq(r_lo, r);
+    }
+
     function test512Math_oshrUp(uint256 x_hi, uint256 x_lo, uint256 s) external pure {
         s = bound(s, 0, 512);