From fcf5e5a377b1f7dce30389088b4d73283952813f Mon Sep 17 00:00:00 2001
From: Stefan Lindblad <stefan.lindblad@linux.com>
Date: Sun, 11 Feb 2024 14:43:38 +0100
Subject: [PATCH] Use u128 sqrt from the `num` crate

- Remove the local unoptimized square root implementation for unsigned 128 bit integers and import `num` crate instead
- Update to the latest `stderrlog` version
- Update README with up-to-date time estiamtes
---
 Cargo.toml    |  3 ++-
 README.md     |  4 ++--
 src/lib.rs    | 21 ++++-----------------
 tests/test.rs | 25 +------------------------
 4 files changed, 9 insertions(+), 44 deletions(-)

diff --git a/Cargo.toml b/Cargo.toml
index 671d22b..9c7d6c1 100644
--- a/Cargo.toml
+++ b/Cargo.toml
@@ -25,9 +25,10 @@ harness = false
 [dependencies]
 clap = "4.4"
 log = "0.4"
-stderrlog = "0.5"
+stderrlog = "0.6"
 genawaiter = "0.99"
 rayon = "1.8"
+num = "0.4.1"
 
 [dev-dependencies]
 criterion = "0.5"
diff --git a/README.md b/README.md
index 8bee9b1..4c067c6 100644
--- a/README.md
+++ b/README.md
@@ -30,8 +30,8 @@ On an old system (i7-6700):
 - full test suite completes in about 7 minutes
 
 On a modern system (i7-12700):
-- 32-bit, random number in about 7 us and worst case in 66.5 us
-- 64-bit, random number in about 348 ms ([58 .. 765] ms) and worst case in 4.44 s
+- 32-bit, random number in about 6.5 us and worst case in 68 us
+- 64-bit, random number in about 140 ms ([3 .. 340] ms) and worst case in 4.6 s
 - full test suite completes in less than 3 minutes
 
 The above numbers are taken from the included benchmark test, which you can run with the command: `cargo bench`. Note that it will take a few minutes to run the full suite, in which time you need to close all other applications and leave it unattended, to give the benchmark the most processing power possible.
diff --git a/src/lib.rs b/src/lib.rs
index cd15001..f6eb384 100644
--- a/src/lib.rs
+++ b/src/lib.rs
@@ -6,6 +6,7 @@ use std::convert::From;
 use std::fmt;
 use genawaiter::stack::let_gen_using;
 use candidates::prime_wheel_30;
+use num::integer::Roots;
 
 #[derive(Clone, Copy, Debug, Eq, Ord, PartialEq, PartialOrd)]
 pub struct IntFactor {
@@ -97,7 +98,7 @@ impl From<u128> for PrimeFactors {
         let mut pf = PrimeFactors::new();
         if n < 2 { return pf; }
         // A factor of n must have a value less than or equal to sqrt(n)
-        let mut maxf = u128_sqrt(n) + 1;
+        let mut maxf = n.sqrt() + 1;
         let_gen_using!(mpgen, prime_wheel_30);
         let mut x = n;
         for f in mpgen.into_iter() {
@@ -111,7 +112,7 @@ impl From<u128> for PrimeFactors {
             }
             if c > 0 {
                 // A factor of x must have a value less than or equal to sqrt(x)
-                maxf = u128_sqrt(x) + 1;
+                maxf = x.sqrt() + 1;
                 pf.add(f, c);
             }
             if x == 1 {
@@ -164,20 +165,6 @@ impl<'a> IntoIterator for &'a PrimeFactors {
     }
 }
 
-/// Unsigned 128-bit integer square root calculation.
-/// Based on example implementation in C at:
-/// https://en.wikipedia.org/wiki/Integer_square_root
-pub fn u128_sqrt(s: u128) -> u128 {
-    let mut g = s >> 1; // Initial guess
-    if g == 0 { return s; } // sanity check
-    let mut u = (g + s / g) >> 1; // update
-    while u < g { // this also checks for cycle
-        g = u;
-        u = (g + s / g) >> 1;
-    }
-    g
-}
-
 /// Test if the value is a prime number, or not
 pub fn u128_is_prime(n: u128) -> bool {
     if n < 2 { return false; }
@@ -189,7 +176,7 @@ pub fn u128_is_prime(n: u128) -> bool {
         }
     }
     // A factor of n must have a value less than or equal to sqrt(n)
-    let maxf = u128_sqrt(n) + 1;
+    let maxf = n.sqrt() + 1;
     let_gen_using!(mpgen, prime_wheel_30);
     for f in mpgen.into_iter() {
         if f >= maxf {
diff --git a/tests/test.rs b/tests/test.rs
index fdb7a69..def1d28 100644
--- a/tests/test.rs
+++ b/tests/test.rs
@@ -9,8 +9,7 @@ use primefactor::{
     PrimeFactors,
     u128_gcd,
     u128_is_prime,
-    u128_lcm,
-    u128_sqrt};
+    u128_lcm};
 
 #[test]
 fn test_early_prime_wheel_numbers() {
@@ -60,28 +59,6 @@ fn test_prime_wheel_quality() {
     assert!(percent > 25.0);
 }
 
-#[test]
-fn test_int_sqrt_pow_of_2() {
-    let mut rnd = rand::thread_rng();
-    for _ in 1..1000 {
-        let n = rnd.gen_range(1..u128_sqrt(u128::MAX));
-        let sqrt = u128_sqrt(n.pow(2));
-        assert_eq!(sqrt, n);
-    }
-}
-
-#[test]
-fn test_int_sqrt_floor() {
-    let mut rnd = rand::thread_rng();
-    for _ in 1..1000 {
-        // Largest integer in a f64 is 2^53-1 (52 bits mantissa)
-        let n = rnd.gen_range(1..u64::pow(2, 53) as u128);
-        let expt = f64::sqrt(n as f64) as u128;
-        let sqrt = u128_sqrt(n);
-        assert_eq!(sqrt, expt);
-    }
-}
-
 #[test]
 fn test_is_prime() {
     for num in 2..=1000 {