Merge pull request #80 from cuviper/karatsuba-sqrt

implement Karatsuba Square Root

Author: cuviper

Cargo.toml

@@ -8,7 +8,7 @@ categories = ["algorithms", "science", "no-std"]
 license = "MIT OR Apache-2.0"
 repository = "https://github.com/rust-num/num-integer"
 name = "num-integer"
-version = "0.1.46"
+version = "0.1.47"
 readme = "README.md"
 exclude = ["/ci/*", "/.github/*"]
 edition = "2018"

RELEASES.md

@@ -1,3 +1,12 @@
+# Release 0.1.47 (2025-08-11)
+
+- [Implement Karatsuba Square Root for 128-bit `Roots::sqrt`][80]
+- Miscellaneous other improvements to docs and testing.
+
+**Contributors**: @durin42, @cuviper, @mikem8891
+
+[80]: https://github.com/rust-num/num-integer/pull/80
+
 # Release 0.1.46 (2024-02-07)
 
 - [Upgrade to 2018 edition, **MSRV 1.31**][51]

src/roots.rs

@@ -196,6 +196,44 @@ fn log2<T: PrimInt>(x: T) -> u32 {
     bits::<T>() - 1 - x.leading_zeros()
 }
 
+/// 128-bit Karatsuba Square Root, using b = 2³²
+///
+/// Reference:
+/// Paul Zimmermann. Karatsuba Square Root. [Research Report] RR-3805, INRIA. 1999, pp.8.
+/// <https://inria.hal.science/inria-00072854/en/>
+#[inline]
+fn karatsuba_sqrt(n: u128) -> u128 {
+    // Algorithm SqrtRem(n = a₃b³ + a₂b² + a₁b + a₀)
+    // Input: 0 ≤ aᵢ < b with a₃ ≥ b/4
+    // Output: (s,r) such that s² ≤ n = s² + r < (s+1)²
+    debug_assert!(n.leading_zeros() < 2);
+    let a0 = n as u32 as u128;
+    let a1 = (n >> 32) as u32 as u128;
+    let a23 = n >> 64;
+
+    // (s',r') ← SqrtRem(a₃b + a₂)
+    let s1 = (a23 as u64).sqrt() as u128;
+    let r1 = a23 - s1 * s1;
+
+    // (q,u) ← DivRem(r'b + a₁, 2s')
+    let (q, u) = ((r1 << 32) | a1).div_rem(&(2 * s1));
+
+    // s ← s'b + q
+    let mut s = (s1 << 32) + q;
+
+    // r ← ub + a₀ - q²
+    // if r < 0 then
+    //   r ← r + 2s - 1
+    //   s ← s - 1
+    //
+    // but to avoid negatives, we compare and adjust before subtraction,
+    // and in this case we don't care about the actual remainder.
+    if ((u << 32) | a0) < q * q {
+        s -= 1;
+    }
+    s
+}
+
 macro_rules! unsigned_roots {
     ($T:ident) => {
         impl Roots for $T {
@@ -275,17 +313,14 @@ macro_rules! unsigned_roots {
             fn sqrt(&self) -> Self {
                 fn go(a: $T) -> $T {
                     if bits::<$T>() > 64 {
-                        // 128-bit division is slow, so do a bitwise `sqrt` until it's small enough.
+                        // 128-bit division is slow in the Babylonian method,
+                        // so use 64-bit sqrt and Karatsuba if needed.
                         return if a <= core::u64::MAX as $T {
                             (a as u64).sqrt() as $T
                         } else {
-                            let lo = (a >> 2u32).sqrt() << 1;
-                            let hi = lo + 1;
-                            if hi * hi <= a {
-                                hi
-                            } else {
-                                lo
-                            }
+                            let shift = a.leading_zeros() / 2;
+                            let n = a << (shift * 2);
+                            (karatsuba_sqrt(n as u128) >> shift) as $T
                         };
                     }