From d9095a25392e376da06753f219e9976456ecc6fa Mon Sep 17 00:00:00 2001
From: nicoo <nicoo@debian.org>
Date: Sun, 24 May 2020 11:12:39 +0200
Subject: [PATCH 01/15] factor: Refactor (eheh) around a `Factors` datatype

It is clearer to see what is going on, as opposed to passing around an
unmarked `Vec<u64>`, and there is a single place to add invariants checks.

This is also a more compact memory representation: each prime factor is
represented only once, with an additional byte for multiplicity.  The
performance impact is however not significant.
---
 src/uu/factor/src/factor.rs | 97 +++++++++++++++++++++++++++----------
 1 file changed, 72 insertions(+), 25 deletions(-)

diff --git a/src/uu/factor/src/factor.rs b/src/uu/factor/src/factor.rs
index 81f655eb74..334fe644a3 100644
--- a/src/uu/factor/src/factor.rs
+++ b/src/uu/factor/src/factor.rs
@@ -23,9 +23,12 @@ use rand::distributions::{Distribution, Uniform};
 use rand::rngs::SmallRng;
 use rand::{thread_rng, SeedableRng};
 use std::cmp::{max, min};
+use std::collections::HashMap;
+use std::fmt;
 use std::io::{stdin, BufRead};
 use std::mem::swap;
 use std::num::Wrapping;
+use std::ops;
 
 mod numeric;
 
@@ -36,14 +39,6 @@ static SUMMARY: &str = "Print the prime factors of the given number(s).
  If none are specified, read from standard input.";
 static LONG_HELP: &str = "";
 
-fn rho_pollard_pseudorandom_function(x: u64, a: u64, b: u64, num: u64) -> u64 {
-    if num < 1 << 63 {
-        (sm_mul(a, sm_mul(x, x, num), num) + b) % num
-    } else {
-        big_add(big_mul(a, big_mul(x, x, num), num), b, num)
-    }
-}
-
 fn gcd(mut a: u64, mut b: u64) -> u64 {
     while b > 0 {
         a %= b;
@@ -52,6 +47,58 @@ fn gcd(mut a: u64, mut b: u64) -> u64 {
     a
 }
 
+struct Factors {
+    f: HashMap<u64, u8>,
+}
+
+impl Factors {
+    fn new() -> Factors {
+        Factors { f: HashMap::new() }
+    }
+
+    fn add(&mut self, prime: u64, exp: u8) {
+        assert!(exp > 0);
+        self.f.insert(prime, exp + self.f.get(&prime).unwrap_or(&0));
+    }
+
+    fn push(&mut self, prime: u64) {
+        self.add(prime, 1)
+    }
+}
+
+impl ops::MulAssign<Factors> for Factors {
+    fn mul_assign(&mut self, other: Factors) {
+        for (prime, exp) in &other.f {
+            self.add(*prime, *exp);
+        }
+    }
+}
+
+impl fmt::Display for Factors {
+    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
+        // TODO: Use a representation with efficient in-order iteration
+        let mut primes: Vec<&u64> = self.f.keys().collect();
+        primes.sort();
+
+        for p in primes {
+            for _ in 0..self.f[&p] {
+                write!(f, " {}", p)?
+            }
+        }
+
+        Ok(())
+    }
+}
+
+
+fn rho_pollard_pseudorandom_function(x: u64, a: u64, b: u64, num: u64) -> u64 {
+    if num < 1 << 63 {
+        (sm_mul(a, sm_mul(x, x, num), num) + b) % num
+    } else {
+        big_add(big_mul(a, big_mul(x, x, num), num), b, num)
+    }
+}
+
 fn rho_pollard_find_divisor(num: u64) -> u64 {
     #![allow(clippy::many_single_char_names)]
     let range = Uniform::new(1, num);
@@ -78,31 +125,39 @@ fn rho_pollard_find_divisor(num: u64) -> u64 {
     }
 }
 
-fn rho_pollard_factor(num: u64, factors: &mut Vec<u64>) {
+fn rho_pollard_factor(num: u64, factors: &mut Factors) {
     if is_prime(num) {
         factors.push(num);
         return;
     }
+
     let divisor = rho_pollard_find_divisor(num);
     rho_pollard_factor(divisor, factors);
     rho_pollard_factor(num / divisor, factors);
 }
 
-fn table_division(mut num: u64, factors: &mut Vec<u64>) {
+fn table_division(mut num: u64) -> Factors {
+    let mut factors = Factors::new();
+
     if num < 2 {
-        return;
+        factors.push(num);
+        return factors
     }
+
     while num % 2 == 0 {
         num /= 2;
         factors.push(2);
     }
+
     if num == 1 {
-        return;
+        return factors;
     }
+
     if is_prime(num) {
         factors.push(num);
-        return;
+        return factors;
     }
+
     for &(prime, inv, ceil) in P_INVS_U64 {
         if num == 1 {
             break;
@@ -120,7 +175,7 @@ fn table_division(mut num: u64, factors: &mut Vec<u64>) {
                 factors.push(prime);
                 if is_prime(num) {
                     factors.push(num);
-                    return;
+                    return factors;
                 }
             } else {
                 break;
@@ -136,21 +191,13 @@ fn table_division(mut num: u64, factors: &mut Vec<u64>) {
     //trial_division_slow(num, factors);
     //} else if num > 1 {
     // number is still greater than 1, but not so big that we have to worry
-    rho_pollard_factor(num, factors);
+    rho_pollard_factor(num, &mut factors);
+    factors
     //}
 }
 
 fn print_factors(num: u64) {
-    print!("{}:", num);
-
-    let mut factors = Vec::new();
-    // we always start with table division, and go from there
-    table_division(num, &mut factors);
-    factors.sort();
-
-    for fac in &factors {
-        print!(" {}", fac);
-    }
+    print!("{}:{}", num, table_division(num));
     println!();
 }
 

From a1b25227507c36709d7fc2693f36bab05dde5d05 Mon Sep 17 00:00:00 2001
From: nicoo <nicoo@debian.org>
Date: Sun, 24 May 2020 15:41:23 +0200
Subject: [PATCH 02/15] factor: Move each factorisation method to its own
 module

Also decoupled the factorisation methods; now factor::factor contains
the logic that chains the different algorithms and aggregates results.

As a side-effect, rho::factor now performs extraneous allocations (as each
recursive step creates a new `Factors` value, which is then aggregated into
the previous one) but there is no significant performance impact.
---
 src/uu/factor/src/factor.rs  | 123 +++++------------------------------
 src/uu/factor/src/numeric.rs |   9 +++
 src/uu/factor/src/rho.rs     |  54 +++++++++++++++
 src/uu/factor/src/table.rs   |  44 +++++++++++++
 4 files changed, 122 insertions(+), 108 deletions(-)
 create mode 100644 src/uu/factor/src/rho.rs
 create mode 100644 src/uu/factor/src/table.rs

diff --git a/src/uu/factor/src/factor.rs b/src/uu/factor/src/factor.rs
index 334fe644a3..f01cae76d7 100644
--- a/src/uu/factor/src/factor.rs
+++ b/src/uu/factor/src/factor.rs
@@ -18,35 +18,20 @@ extern crate rand;
 #[macro_use]
 extern crate uucore;
 
-use numeric::*;
-use rand::distributions::{Distribution, Uniform};
-use rand::rngs::SmallRng;
-use rand::{thread_rng, SeedableRng};
-use std::cmp::{max, min};
 use std::collections::HashMap;
 use std::fmt;
 use std::io::{stdin, BufRead};
-use std::mem::swap;
-use std::num::Wrapping;
 use std::ops;
 
 mod numeric;
-
-include!(concat!(env!("OUT_DIR"), "/prime_table.rs"));
+mod rho;
+mod table;
 
 static SYNTAX: &str = "[OPTION] [NUMBER]...";
 static SUMMARY: &str = "Print the prime factors of the given number(s).
  If none are specified, read from standard input.";
 static LONG_HELP: &str = "";
 
-fn gcd(mut a: u64, mut b: u64) -> u64 {
-    while b > 0 {
-        a %= b;
-        swap(&mut a, &mut b);
-    }
-    a
-}
-
 struct Factors {
     f: HashMap<u64, u8>,
 }
@@ -90,114 +75,36 @@ impl fmt::Display for Factors {
     }
 }
 
-
-fn rho_pollard_pseudorandom_function(x: u64, a: u64, b: u64, num: u64) -> u64 {
-    if num < 1 << 63 {
-        (sm_mul(a, sm_mul(x, x, num), num) + b) % num
-    } else {
-        big_add(big_mul(a, big_mul(x, x, num), num), b, num)
-    }
-}
-
-fn rho_pollard_find_divisor(num: u64) -> u64 {
-    #![allow(clippy::many_single_char_names)]
-    let range = Uniform::new(1, num);
-    let mut rng = SmallRng::from_rng(&mut thread_rng()).unwrap();
-    let mut x = range.sample(&mut rng);
-    let mut y = x;
-    let mut a = range.sample(&mut rng);
-    let mut b = range.sample(&mut rng);
-
-    loop {
-        x = rho_pollard_pseudorandom_function(x, a, b, num);
-        y = rho_pollard_pseudorandom_function(y, a, b, num);
-        y = rho_pollard_pseudorandom_function(y, a, b, num);
-        let d = gcd(num, max(x, y) - min(x, y));
-        if d == num {
-            // Failure, retry with different function
-            x = range.sample(&mut rng);
-            y = x;
-            a = range.sample(&mut rng);
-            b = range.sample(&mut rng);
-        } else if d > 1 {
-            return d;
-        }
-    }
-}
-
-fn rho_pollard_factor(num: u64, factors: &mut Factors) {
-    if is_prime(num) {
-        factors.push(num);
-        return;
-    }
-
-    let divisor = rho_pollard_find_divisor(num);
-    rho_pollard_factor(divisor, factors);
-    rho_pollard_factor(num / divisor, factors);
-}
-
-fn table_division(mut num: u64) -> Factors {
+fn factor(mut n: u64) -> Factors {
     let mut factors = Factors::new();
 
-    if num < 2 {
-        factors.push(num);
-        return factors
+    if n < 2 {
+        factors.push(n);
+        return factors;
     }
 
-    while num % 2 == 0 {
-        num /= 2;
+    while n % 2 == 0 {
+        n /= 2;
         factors.push(2);
     }
 
-    if num == 1 {
+    if n == 1 {
         return factors;
     }
 
-    if is_prime(num) {
-        factors.push(num);
+    if numeric::is_prime(n) {
+        factors.push(n);
         return factors;
     }
 
-    for &(prime, inv, ceil) in P_INVS_U64 {
-        if num == 1 {
-            break;
-        }
-
-        // inv = prime^-1 mod 2^64
-        // ceil = floor((2^64-1) / prime)
-        // if (num * inv) mod 2^64 <= ceil, then prime divides num
-        // See http://math.stackexchange.com/questions/1251327/
-        // for a nice explanation.
-        loop {
-            let Wrapping(x) = Wrapping(num) * Wrapping(inv); // x = num * inv mod 2^64
-            if x <= ceil {
-                num = x;
-                factors.push(prime);
-                if is_prime(num) {
-                    factors.push(num);
-                    return factors;
-                }
-            } else {
-                break;
-            }
-        }
-    }
-
-    // do we still have more factoring to do?
-    // Decide whether to use Pollard Rho or slow divisibility based on
-    // number's size:
-    //if num >= 1 << 63 {
-    // number is too big to use rho pollard without overflowing
-    //trial_division_slow(num, factors);
-    //} else if num > 1 {
-    // number is still greater than 1, but not so big that we have to worry
-    rho_pollard_factor(num, &mut factors);
+    let (f, n) = table::factor(n);
+    factors *= f;
+    factors *= rho::factor(n);
     factors
-    //}
 }
 
 fn print_factors(num: u64) {
-    print!("{}:{}", num, table_division(num));
+    print!("{}:{}", num, factor(num));
     println!();
 }
 
diff --git a/src/uu/factor/src/numeric.rs b/src/uu/factor/src/numeric.rs
index 474c806131..d9e6f94416 100644
--- a/src/uu/factor/src/numeric.rs
+++ b/src/uu/factor/src/numeric.rs
@@ -9,9 +9,18 @@
 * that was distributed with this source code.
 */
 
+use std::mem::swap;
 use std::num::Wrapping;
 use std::u64::MAX as MAX_U64;
 
+pub fn gcd(mut a: u64, mut b: u64) -> u64 {
+    while b > 0 {
+        a %= b;
+        swap(&mut a, &mut b);
+    }
+    a
+}
+
 pub fn big_add(a: u64, b: u64, m: u64) -> u64 {
     let Wrapping(msb_mod_m) = Wrapping(MAX_U64) - Wrapping(m) + Wrapping(1);
     let msb_mod_m = msb_mod_m % m;
diff --git a/src/uu/factor/src/rho.rs b/src/uu/factor/src/rho.rs
new file mode 100644
index 0000000000..9494977374
--- /dev/null
+++ b/src/uu/factor/src/rho.rs
@@ -0,0 +1,54 @@
+use crate::Factors;
+use numeric::*;
+use rand::distributions::{Distribution, Uniform};
+use rand::rngs::SmallRng;
+use rand::{thread_rng, SeedableRng};
+use std::cmp::{max, min};
+
+fn pseudorandom_function(x: u64, a: u64, b: u64, num: u64) -> u64 {
+    if num < 1 << 63 {
+        (sm_mul(a, sm_mul(x, x, num), num) + b) % num
+    } else {
+        big_add(big_mul(a, big_mul(x, x, num), num), b, num)
+    }
+}
+
+fn find_divisor(num: u64) -> u64 {
+    #![allow(clippy::many_single_char_names)]
+    let range = Uniform::new(1, num);
+    let mut rng = SmallRng::from_rng(&mut thread_rng()).unwrap();
+    let mut x = range.sample(&mut rng);
+    let mut y = x;
+    let mut a = range.sample(&mut rng);
+    let mut b = range.sample(&mut rng);
+
+    loop {
+        x = pseudorandom_function(x, a, b, num);
+        y = pseudorandom_function(y, a, b, num);
+        y = pseudorandom_function(y, a, b, num);
+        let d = gcd(num, max(x, y) - min(x, y));
+        if d == num {
+            // Failure, retry with different function
+            x = range.sample(&mut rng);
+            y = x;
+            a = range.sample(&mut rng);
+            b = range.sample(&mut rng);
+        } else if d > 1 {
+            return d;
+        }
+    }
+}
+
+pub(crate) fn factor(num: u64) -> Factors {
+    let mut factors = Factors::new();
+    if num == 1 { return factors; }
+    if is_prime(num) {
+        factors.push(num);
+        return factors;
+    }
+
+    let divisor = find_divisor(num);
+    factors *= factor(divisor);
+    factors *= factor(num / divisor);
+    factors
+}
diff --git a/src/uu/factor/src/table.rs b/src/uu/factor/src/table.rs
new file mode 100644
index 0000000000..d8cfa1236c
--- /dev/null
+++ b/src/uu/factor/src/table.rs
@@ -0,0 +1,44 @@
+use crate::Factors;
+use numeric::*;
+use std::num::Wrapping;
+
+include!(concat!(env!("OUT_DIR"), "/prime_table.rs"));
+
+pub(crate) fn factor(mut num: u64) -> (Factors, u64) {
+    let mut factors = Factors::new();
+    for &(prime, inv, ceil) in P_INVS_U64 {
+        if num == 1 {
+            break;
+        }
+
+        // inv = prime^-1 mod 2^64
+        // ceil = floor((2^64-1) / prime)
+        // if (num * inv) mod 2^64 <= ceil, then prime divides num
+        // See http://math.stackexchange.com/questions/1251327/
+        // for a nice explanation.
+        loop {
+            let Wrapping(x) = Wrapping(num) * Wrapping(inv); // x = num * inv mod 2^64
+            if x <= ceil {
+                num = x;
+                factors.push(prime);
+                if is_prime(num) {
+                    factors.push(num);
+                    return (factors, 1);
+                }
+            } else {
+                break;
+            }
+        }
+    }
+
+    // do we still have more factoring to do?
+    // Decide whether to use Pollard Rho or slow divisibility based on
+    // number's size:
+    //if num >= 1 << 63 {
+    // number is too big to use rho pollard without overflowing
+    //trial_division_slow(num, factors);
+    //} else if num > 1 {
+    // number is still greater than 1, but not so big that we have to worry
+    (factors, num)
+    //}
+}

From bc11e5796220a69041dc4bbf7cd7db42bbe01cb1 Mon Sep 17 00:00:00 2001
From: nicoo <nicoo@debian.org>
Date: Sun, 24 May 2020 15:41:23 +0200
Subject: [PATCH 03/15] factor::factor: Use u64::trailing_zero instead of
 iterated division

No significant performance impact (most of the time is spent elsewhere),
but an easy and satisfying fix nevertheless.
---
 src/uu/factor/src/factor.rs | 7 ++++---
 1 file changed, 4 insertions(+), 3 deletions(-)

diff --git a/src/uu/factor/src/factor.rs b/src/uu/factor/src/factor.rs
index f01cae76d7..894c61c81a 100644
--- a/src/uu/factor/src/factor.rs
+++ b/src/uu/factor/src/factor.rs
@@ -83,9 +83,10 @@ fn factor(mut n: u64) -> Factors {
         return factors;
     }
 
-    while n % 2 == 0 {
-        n /= 2;
-        factors.push(2);
+    let z = n.trailing_zeros();
+    if z > 0 {
+        factors.add(2, z as u8);
+        n = n >> z;
     }
 
     if n == 1 {

From 418fd6175952c4147de24407d224b140c843b619 Mon Sep 17 00:00:00 2001
From: nicoo <nicoo@debian.org>
Date: Sun, 24 May 2020 15:41:23 +0200
Subject: [PATCH 04/15] factor::factor: Short-circuit the fallback to Pollard's
 rho

When the remainder is smaller than the max. entry in the table,
it is guaranteed to be prime.
---
 src/uu/factor/src/factor.rs | 6 +++++-
 1 file changed, 5 insertions(+), 1 deletion(-)

diff --git a/src/uu/factor/src/factor.rs b/src/uu/factor/src/factor.rs
index 894c61c81a..1e5b3bb8f3 100644
--- a/src/uu/factor/src/factor.rs
+++ b/src/uu/factor/src/factor.rs
@@ -100,7 +100,11 @@ fn factor(mut n: u64) -> Factors {
 
     let (f, n) = table::factor(n);
     factors *= f;
-    factors *= rho::factor(n);
+
+    if n >= table::NEXT_PRIME {
+        factors *= rho::factor(n);
+    }
+
     factors
 }
 

From 169740629bdee46595c6a452a1ee604a9a96b868 Mon Sep 17 00:00:00 2001
From: nicoo <nicoo@debian.org>
Date: Sun, 24 May 2020 15:41:23 +0200
Subject: [PATCH 05/15] factor::table: Remove extraneous calls to the primality
 test
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

50% performance improvement on factoring all numbers between 2 and 10⁶.
---
 src/uu/factor/src/table.rs | 11 ++++-------
 1 file changed, 4 insertions(+), 7 deletions(-)

diff --git a/src/uu/factor/src/table.rs b/src/uu/factor/src/table.rs
index d8cfa1236c..45f5c2d6af 100644
--- a/src/uu/factor/src/table.rs
+++ b/src/uu/factor/src/table.rs
@@ -1,5 +1,4 @@
 use crate::Factors;
-use numeric::*;
 use std::num::Wrapping;
 
 include!(concat!(env!("OUT_DIR"), "/prime_table.rs"));
@@ -14,17 +13,15 @@ pub(crate) fn factor(mut num: u64) -> (Factors, u64) {
         // inv = prime^-1 mod 2^64
         // ceil = floor((2^64-1) / prime)
         // if (num * inv) mod 2^64 <= ceil, then prime divides num
-        // See http://math.stackexchange.com/questions/1251327/
+        // See https://math.stackexchange.com/questions/1251327/
         // for a nice explanation.
         loop {
-            let Wrapping(x) = Wrapping(num) * Wrapping(inv); // x = num * inv mod 2^64
+            let Wrapping(x) = Wrapping(num) * Wrapping(inv);
+
+            // While prime divides num
             if x <= ceil {
                 num = x;
                 factors.push(prime);
-                if is_prime(num) {
-                    factors.push(num);
-                    return (factors, 1);
-                }
             } else {
                 break;
             }

From e1a6dbe6193f97727e0a1493151dd6cb1119a65d Mon Sep 17 00:00:00 2001
From: nicoo <nicoo@debian.org>
Date: Sun, 24 May 2020 15:41:23 +0200
Subject: [PATCH 06/15] factor::table: Remove obsolete, commented code

---
 src/uu/factor/src/table.rs | 9 ---------
 1 file changed, 9 deletions(-)

diff --git a/src/uu/factor/src/table.rs b/src/uu/factor/src/table.rs
index 45f5c2d6af..3a4f44e85c 100644
--- a/src/uu/factor/src/table.rs
+++ b/src/uu/factor/src/table.rs
@@ -28,14 +28,5 @@ pub(crate) fn factor(mut num: u64) -> (Factors, u64) {
         }
     }
 
-    // do we still have more factoring to do?
-    // Decide whether to use Pollard Rho or slow divisibility based on
-    // number's size:
-    //if num >= 1 << 63 {
-    // number is too big to use rho pollard without overflowing
-    //trial_division_slow(num, factors);
-    //} else if num > 1 {
-    // number is still greater than 1, but not so big that we have to worry
     (factors, num)
-    //}
 }

From 74054feb94173fdbc40a9e54e9b37fa579c61e38 Mon Sep 17 00:00:00 2001
From: nicoo <nicoo@debian.org>
Date: Sun, 24 May 2020 15:41:23 +0200
Subject: [PATCH 07/15] factor::factor: Remove extraneous call to the primality
 test

Another 6.97% runtime improvement
---
 src/uu/factor/src/factor.rs | 5 -----
 1 file changed, 5 deletions(-)

diff --git a/src/uu/factor/src/factor.rs b/src/uu/factor/src/factor.rs
index 1e5b3bb8f3..afaf4f8cdb 100644
--- a/src/uu/factor/src/factor.rs
+++ b/src/uu/factor/src/factor.rs
@@ -93,11 +93,6 @@ fn factor(mut n: u64) -> Factors {
         return factors;
     }
 
-    if numeric::is_prime(n) {
-        factors.push(n);
-        return factors;
-    }
-
     let (f, n) = table::factor(n);
     factors *= f;
 

From e3ecc81d974330858e74a8369c2dda89d97eff02 Mon Sep 17 00:00:00 2001
From: nicoo <nicoo@debian.org>
Date: Sun, 24 May 2020 15:41:23 +0200
Subject: [PATCH 08/15] factor: Move the Miller-Rabin primality test to its own
 module.

---
 src/uu/factor/build.rs            |  2 +-
 src/uu/factor/src/factor.rs       |  1 +
 src/uu/factor/src/miller_rabin.rs | 52 +++++++++++++++++++++++++++++++
 src/uu/factor/src/numeric.rs      | 52 +------------------------------
 src/uu/factor/src/rho.rs          |  9 ++++--
 5 files changed, 61 insertions(+), 55 deletions(-)
 create mode 100644 src/uu/factor/src/miller_rabin.rs

diff --git a/src/uu/factor/build.rs b/src/uu/factor/build.rs
index acc4d9a3f4..3a1b7df94e 100644
--- a/src/uu/factor/build.rs
+++ b/src/uu/factor/build.rs
@@ -26,7 +26,7 @@ use std::path::Path;
 use std::u64::MAX as MAX_U64;
 
 #[cfg(test)]
-use numeric::is_prime;
+use miller_rabin::is_prime;
 
 #[cfg(test)]
 #[path = "src/numeric.rs"]
diff --git a/src/uu/factor/src/factor.rs b/src/uu/factor/src/factor.rs
index afaf4f8cdb..5a458a1432 100644
--- a/src/uu/factor/src/factor.rs
+++ b/src/uu/factor/src/factor.rs
@@ -23,6 +23,7 @@ use std::fmt;
 use std::io::{stdin, BufRead};
 use std::ops;
 
+mod miller_rabin;
 mod numeric;
 mod rho;
 mod table;
diff --git a/src/uu/factor/src/miller_rabin.rs b/src/uu/factor/src/miller_rabin.rs
new file mode 100644
index 0000000000..c42aee5006
--- /dev/null
+++ b/src/uu/factor/src/miller_rabin.rs
@@ -0,0 +1,52 @@
+use crate::numeric::*;
+
+fn witness(mut a: u64, exponent: u64, m: u64) -> bool {
+    if a == 0 {
+        return false;
+    }
+
+    let mul = if m < 1 << 63 {
+        sm_mul as fn(u64, u64, u64) -> u64
+    } else {
+        big_mul as fn(u64, u64, u64) -> u64
+    };
+
+    if pow(a, m - 1, m, mul) != 1 {
+        return true;
+    }
+    a = pow(a, exponent, m, mul);
+    if a == 1 {
+        return false;
+    }
+    loop {
+        if a == 1 {
+            return true;
+        }
+        if a == m - 1 {
+            return false;
+        }
+        a = mul(a, a, m);
+    }
+}
+
+// Deterministic Miller-Rabin primality-checking algorithm, adapted to extract
+// (some) dividers; it will fail to factor strong pseudoprimes.
+pub(crate) fn is_prime(num: u64) -> bool {
+    if num < 2 {
+        return false;
+    }
+    if num % 2 == 0 {
+        return num == 2;
+    }
+    let mut exponent = num - 1;
+    while exponent & 1 == 0 {
+        exponent >>= 1;
+    }
+
+    // These witnesses detect all composites up to at least 2^64.
+    // Discovered by Jim Sinclair, according to http://miller-rabin.appspot.com
+    let witnesses = [2, 325, 9_375, 28_178, 450_775, 9_780_504, 1_795_265_022];
+    !witnesses
+        .iter()
+        .any(|&wit| witness(wit % num, exponent, num))
+}
diff --git a/src/uu/factor/src/numeric.rs b/src/uu/factor/src/numeric.rs
index d9e6f94416..7dfba41fbf 100644
--- a/src/uu/factor/src/numeric.rs
+++ b/src/uu/factor/src/numeric.rs
@@ -78,7 +78,7 @@ pub fn big_mul(mut a: u64, mut b: u64, m: u64) -> u64 {
 }
 
 // computes a.pow(b) % m
-fn pow(mut a: u64, mut b: u64, m: u64, mul: fn(u64, u64, u64) -> u64) -> u64 {
+pub(crate) fn pow(mut a: u64, mut b: u64, m: u64, mul: fn(u64, u64, u64) -> u64) -> u64 {
     let mut result = 1;
     while b > 0 {
         if b & 1 != 0 {
@@ -89,53 +89,3 @@ fn pow(mut a: u64, mut b: u64, m: u64, mul: fn(u64, u64, u64) -> u64) -> u64 {
     }
     result
 }
-
-fn witness(mut a: u64, exponent: u64, m: u64) -> bool {
-    if a == 0 {
-        return false;
-    }
-
-    let mul = if m < 1 << 63 {
-        sm_mul as fn(u64, u64, u64) -> u64
-    } else {
-        big_mul as fn(u64, u64, u64) -> u64
-    };
-
-    if pow(a, m - 1, m, mul) != 1 {
-        return true;
-    }
-    a = pow(a, exponent, m, mul);
-    if a == 1 {
-        return false;
-    }
-    loop {
-        if a == 1 {
-            return true;
-        }
-        if a == m - 1 {
-            return false;
-        }
-        a = mul(a, a, m);
-    }
-}
-
-// uses deterministic (i.e., fixed witness set) Miller-Rabin test
-pub fn is_prime(num: u64) -> bool {
-    if num < 2 {
-        return false;
-    }
-    if num % 2 == 0 {
-        return num == 2;
-    }
-    let mut exponent = num - 1;
-    while exponent & 1 == 0 {
-        exponent >>= 1;
-    }
-
-    // These witnesses detect all composites up to at least 2^64.
-    // Discovered by Jim Sinclair, according to http://miller-rabin.appspot.com
-    let witnesses = [2, 325, 9_375, 28_178, 450_775, 9_780_504, 1_795_265_022];
-    !witnesses
-        .iter()
-        .any(|&wit| witness(wit % num, exponent, num))
-}
diff --git a/src/uu/factor/src/rho.rs b/src/uu/factor/src/rho.rs
index 9494977374..947eac4786 100644
--- a/src/uu/factor/src/rho.rs
+++ b/src/uu/factor/src/rho.rs
@@ -1,4 +1,4 @@
-use crate::Factors;
+use crate::{miller_rabin, Factors};
 use numeric::*;
 use rand::distributions::{Distribution, Uniform};
 use rand::rngs::SmallRng;
@@ -41,8 +41,11 @@ fn find_divisor(num: u64) -> u64 {
 
 pub(crate) fn factor(num: u64) -> Factors {
     let mut factors = Factors::new();
-    if num == 1 { return factors; }
-    if is_prime(num) {
+    if num == 1 {
+        return factors;
+    }
+
+    if miller_rabin::is_prime(num) {
         factors.push(num);
         return factors;
     }

From 6b9585b1dc586ec1109145bee183d1cf75d8d933 Mon Sep 17 00:00:00 2001
From: nicoo <nicoo@debian.org>
Date: Sun, 24 May 2020 15:41:23 +0200
Subject: [PATCH 09/15] factor::miller_rabbin: Refactor before extracting
 dividers

Replace iterated division with u64::trailing_zeros, hoist the selection of `mul`
out of the loop, another cool 49.5% runtime improvement.
---
 src/uu/factor/src/miller_rabin.rs | 73 ++++++++++++++++---------------
 1 file changed, 38 insertions(+), 35 deletions(-)

diff --git a/src/uu/factor/src/miller_rabin.rs b/src/uu/factor/src/miller_rabin.rs
index c42aee5006..997cd1cd6e 100644
--- a/src/uu/factor/src/miller_rabin.rs
+++ b/src/uu/factor/src/miller_rabin.rs
@@ -1,52 +1,55 @@
 use crate::numeric::*;
 
-fn witness(mut a: u64, exponent: u64, m: u64) -> bool {
-    if a == 0 {
+// Small set of bases for the Miller-Rabin prime test, valid for all 64b integers;
+//  discovered by Jim Sinclair on 2011-04-20, see miller-rabin.appspot.com
+const BASIS: [u64; 7] = [2, 325, 9375, 28178, 450775, 9780504, 1795265022];
+
+// Deterministic Miller-Rabin primality-checking algorithm, adapted to extract
+// (some) dividers; it will fail to factor strong pseudoprimes.
+pub(crate) fn is_prime(n: u64) -> bool {
+    if n < 2 {
         return false;
     }
+    if n % 2 == 0 {
+        return n == 2;
+    }
+
+    let d = (n - 1).trailing_zeros();
+    let r = (n - 1) >> d;
 
-    let mul = if m < 1 << 63 {
+    let mul = if n < 1 << 63 {
         sm_mul as fn(u64, u64, u64) -> u64
     } else {
         big_mul as fn(u64, u64, u64) -> u64
     };
 
-    if pow(a, m - 1, m, mul) != 1 {
-        return true;
-    }
-    a = pow(a, exponent, m, mul);
-    if a == 1 {
-        return false;
-    }
-    loop {
-        if a == 1 {
-            return true;
+    for a in BASIS.iter() {
+        let mut x = a % n;
+        if x == 0 {
+            break;
         }
-        if a == m - 1 {
+
+        if pow(x, n - 1, n, mul) != 1 {
             return false;
         }
-        a = mul(a, a, m);
-    }
-}
+        x = pow(x, r, n, mul);
+        if x == 1 || x == n - 1 {
+            break;
+        }
 
-// Deterministic Miller-Rabin primality-checking algorithm, adapted to extract
-// (some) dividers; it will fail to factor strong pseudoprimes.
-pub(crate) fn is_prime(num: u64) -> bool {
-    if num < 2 {
-        return false;
-    }
-    if num % 2 == 0 {
-        return num == 2;
-    }
-    let mut exponent = num - 1;
-    while exponent & 1 == 0 {
-        exponent >>= 1;
+        loop {
+            let y = mul(x, x, n);
+            if y == 1 {
+                return false;
+            }
+            if y == n - 1 {
+                // This basis element is not a witness of `n` being composite.
+                // Keep looking.
+                break;
+            }
+            x = y;
+        }
     }
 
-    // These witnesses detect all composites up to at least 2^64.
-    // Discovered by Jim Sinclair, according to http://miller-rabin.appspot.com
-    let witnesses = [2, 325, 9_375, 28_178, 450_775, 9_780_504, 1_795_265_022];
-    !witnesses
-        .iter()
-        .any(|&wit| witness(wit % num, exponent, num))
+    true
 }

From 824103769028240098c0f30b76fb3e067fedea3d Mon Sep 17 00:00:00 2001
From: nicoo <nicoo@debian.org>
Date: Sun, 24 May 2020 15:41:23 +0200
Subject: [PATCH 10/15] factor::miller_rabin: Extract dividers from the
 primality test

Another 36% improvement.
---
 src/uu/factor/src/miller_rabin.rs | 33 +++++++++++++++++++++++++------
 src/uu/factor/src/rho.rs          | 20 ++++++++++++++-----
 2 files changed, 42 insertions(+), 11 deletions(-)

diff --git a/src/uu/factor/src/miller_rabin.rs b/src/uu/factor/src/miller_rabin.rs
index 997cd1cd6e..42c19b24fe 100644
--- a/src/uu/factor/src/miller_rabin.rs
+++ b/src/uu/factor/src/miller_rabin.rs
@@ -4,14 +4,29 @@ use crate::numeric::*;
 //  discovered by Jim Sinclair on 2011-04-20, see miller-rabin.appspot.com
 const BASIS: [u64; 7] = [2, 325, 9375, 28178, 450775, 9780504, 1795265022];
 
+#[derive(Eq, PartialEq)]
+pub(crate) enum Result {
+    Prime,
+    Pseudoprime,
+    Composite(u64),
+}
+
+impl Result {
+    pub(crate) fn is_prime(&self) -> bool {
+        return *self == Result::Prime;
+    }
+}
+
 // Deterministic Miller-Rabin primality-checking algorithm, adapted to extract
 // (some) dividers; it will fail to factor strong pseudoprimes.
-pub(crate) fn is_prime(n: u64) -> bool {
+pub(crate) fn test(n: u64) -> Result {
+    use self::Result::*;
+
     if n < 2 {
-        return false;
+        return Pseudoprime;
     }
     if n % 2 == 0 {
-        return n == 2;
+        return if n == 2 { Prime } else { Composite(2) };
     }
 
     let d = (n - 1).trailing_zeros();
@@ -30,7 +45,7 @@ pub(crate) fn is_prime(n: u64) -> bool {
         }
 
         if pow(x, n - 1, n, mul) != 1 {
-            return false;
+            return Pseudoprime;
         }
         x = pow(x, r, n, mul);
         if x == 1 || x == n - 1 {
@@ -40,7 +55,7 @@ pub(crate) fn is_prime(n: u64) -> bool {
         loop {
             let y = mul(x, x, n);
             if y == 1 {
-                return false;
+                return Composite(gcd(x - 1, n));
             }
             if y == n - 1 {
                 // This basis element is not a witness of `n` being composite.
@@ -51,5 +66,11 @@ pub(crate) fn is_prime(n: u64) -> bool {
         }
     }
 
-    true
+    Prime
+}
+
+// Used by build.rs' tests
+#[allow(dead_code)]
+pub(crate) fn is_prime(n: u64) -> bool {
+    test(n).is_prime()
 }
diff --git a/src/uu/factor/src/rho.rs b/src/uu/factor/src/rho.rs
index 947eac4786..e1a009810a 100644
--- a/src/uu/factor/src/rho.rs
+++ b/src/uu/factor/src/rho.rs
@@ -1,3 +1,4 @@
+use crate::miller_rabin::Result::*;
 use crate::{miller_rabin, Factors};
 use numeric::*;
 use rand::distributions::{Distribution, Uniform};
@@ -39,16 +40,25 @@ fn find_divisor(num: u64) -> u64 {
     }
 }
 
-pub(crate) fn factor(num: u64) -> Factors {
+pub(crate) fn factor(mut num: u64) -> Factors {
     let mut factors = Factors::new();
     if num == 1 {
         return factors;
     }
 
-    if miller_rabin::is_prime(num) {
-        factors.push(num);
-        return factors;
-    }
+    match miller_rabin::test(num) {
+        Prime => {
+            factors.push(num);
+            return factors;
+        }
+
+        Composite(d) => {
+            num = num / d;
+            factors *= factor(d);
+        }
+
+        Pseudoprime => {}
+    };
 
     let divisor = find_divisor(num);
     factors *= factor(divisor);

From 29eb8fd77bb83349fabf47eb1eee1af427e9b80c Mon Sep 17 00:00:00 2001
From: nicoo <nicoo@debian.org>
Date: Sun, 24 May 2020 17:06:19 +0200
Subject: [PATCH 11/15] format: Make clippy happy

---
 src/uu/factor/src/factor.rs       | 2 +-
 src/uu/factor/src/miller_rabin.rs | 6 +++---
 src/uu/factor/src/rho.rs          | 2 +-
 3 files changed, 5 insertions(+), 5 deletions(-)

diff --git a/src/uu/factor/src/factor.rs b/src/uu/factor/src/factor.rs
index 5a458a1432..493d60912e 100644
--- a/src/uu/factor/src/factor.rs
+++ b/src/uu/factor/src/factor.rs
@@ -87,7 +87,7 @@ fn factor(mut n: u64) -> Factors {
     let z = n.trailing_zeros();
     if z > 0 {
         factors.add(2, z as u8);
-        n = n >> z;
+        n >>= z;
     }
 
     if n == 1 {
diff --git a/src/uu/factor/src/miller_rabin.rs b/src/uu/factor/src/miller_rabin.rs
index 42c19b24fe..9f9dfc913c 100644
--- a/src/uu/factor/src/miller_rabin.rs
+++ b/src/uu/factor/src/miller_rabin.rs
@@ -2,6 +2,7 @@ use crate::numeric::*;
 
 // Small set of bases for the Miller-Rabin prime test, valid for all 64b integers;
 //  discovered by Jim Sinclair on 2011-04-20, see miller-rabin.appspot.com
+#[allow(clippy::unreadable_literal)]
 const BASIS: [u64; 7] = [2, 325, 9375, 28178, 450775, 9780504, 1795265022];
 
 #[derive(Eq, PartialEq)]
@@ -13,7 +14,7 @@ pub(crate) enum Result {
 
 impl Result {
     pub(crate) fn is_prime(&self) -> bool {
-        return *self == Result::Prime;
+        *self == Result::Prime
     }
 }
 
@@ -29,8 +30,7 @@ pub(crate) fn test(n: u64) -> Result {
         return if n == 2 { Prime } else { Composite(2) };
     }
 
-    let d = (n - 1).trailing_zeros();
-    let r = (n - 1) >> d;
+    let r = (n - 1) >> (n - 1).trailing_zeros();
 
     let mul = if n < 1 << 63 {
         sm_mul as fn(u64, u64, u64) -> u64
diff --git a/src/uu/factor/src/rho.rs b/src/uu/factor/src/rho.rs
index e1a009810a..acd07e2da7 100644
--- a/src/uu/factor/src/rho.rs
+++ b/src/uu/factor/src/rho.rs
@@ -53,7 +53,7 @@ pub(crate) fn factor(mut num: u64) -> Factors {
         }
 
         Composite(d) => {
-            num = num / d;
+            num /= d;
             factors *= factor(d);
         }
 

From 30fd6a0309694d0acb9ec956f14786fc8b3ae5c2 Mon Sep 17 00:00:00 2001
From: nicoo <nicoo@debian.org>
Date: Sun, 24 May 2020 18:16:21 +0200
Subject: [PATCH 12/15] factor::numeric: Replace lose functions with an
 Arithmetic trait

---
 src/uu/factor/src/miller_rabin.rs |  21 +++---
 src/uu/factor/src/numeric.rs      | 118 +++++++++++++++++-------------
 src/uu/factor/src/rho.rs          |  61 +++++++--------
 3 files changed, 109 insertions(+), 91 deletions(-)

diff --git a/src/uu/factor/src/miller_rabin.rs b/src/uu/factor/src/miller_rabin.rs
index 9f9dfc913c..cdc1722fde 100644
--- a/src/uu/factor/src/miller_rabin.rs
+++ b/src/uu/factor/src/miller_rabin.rs
@@ -20,7 +20,7 @@ impl Result {
 
 // Deterministic Miller-Rabin primality-checking algorithm, adapted to extract
 // (some) dividers; it will fail to factor strong pseudoprimes.
-pub(crate) fn test(n: u64) -> Result {
+pub(crate) fn test<A: Arithmetic>(n: u64) -> Result {
     use self::Result::*;
 
     if n < 2 {
@@ -32,28 +32,22 @@ pub(crate) fn test(n: u64) -> Result {
 
     let r = (n - 1) >> (n - 1).trailing_zeros();
 
-    let mul = if n < 1 << 63 {
-        sm_mul as fn(u64, u64, u64) -> u64
-    } else {
-        big_mul as fn(u64, u64, u64) -> u64
-    };
-
     for a in BASIS.iter() {
         let mut x = a % n;
         if x == 0 {
             break;
         }
 
-        if pow(x, n - 1, n, mul) != 1 {
+        if A::pow(x, n - 1, n) != 1 {
             return Pseudoprime;
         }
-        x = pow(x, r, n, mul);
+        x = A::pow(x, r, n);
         if x == 1 || x == n - 1 {
             break;
         }
 
         loop {
-            let y = mul(x, x, n);
+            let y = A::mul(x, x, n);
             if y == 1 {
                 return Composite(gcd(x - 1, n));
             }
@@ -72,5 +66,10 @@ pub(crate) fn test(n: u64) -> Result {
 // Used by build.rs' tests
 #[allow(dead_code)]
 pub(crate) fn is_prime(n: u64) -> bool {
-    test(n).is_prime()
+    if n < 1 << 63 {
+        test::<Small>(n)
+    } else {
+        test::<Big>(n)
+    }
+    .is_prime()
 }
diff --git a/src/uu/factor/src/numeric.rs b/src/uu/factor/src/numeric.rs
index 7dfba41fbf..95b75e9862 100644
--- a/src/uu/factor/src/numeric.rs
+++ b/src/uu/factor/src/numeric.rs
@@ -21,71 +21,87 @@ pub fn gcd(mut a: u64, mut b: u64) -> u64 {
     a
 }
 
-pub fn big_add(a: u64, b: u64, m: u64) -> u64 {
-    let Wrapping(msb_mod_m) = Wrapping(MAX_U64) - Wrapping(m) + Wrapping(1);
-    let msb_mod_m = msb_mod_m % m;
+pub(crate) trait Arithmetic {
+    fn add(a: u64, b: u64, modulus: u64) -> u64;
+    fn mul(a: u64, b: u64, modulus: u64) -> u64;
 
-    let Wrapping(res) = Wrapping(a) + Wrapping(b);
-    if b <= MAX_U64 - a {
-        res
-    } else {
-        (res + msb_mod_m) % m
+    fn pow(mut a: u64, mut b: u64, m: u64) -> u64 {
+        let mut result = 1;
+        while b > 0 {
+            if b & 1 != 0 {
+                result = Self::mul(result, a, m);
+            }
+            a = Self::mul(a, a, m);
+            b >>= 1;
+        }
+        result
     }
 }
 
-// computes (a + b) % m using the russian peasant algorithm
-// CAUTION: Will overflow if m >= 2^63
-pub fn sm_mul(mut a: u64, mut b: u64, m: u64) -> u64 {
-    let mut result = 0;
-    while b > 0 {
-        if b & 1 != 0 {
-            result = (result + a) % m;
+pub(crate) struct Big {}
+
+impl Arithmetic for Big {
+    fn add(a: u64, b: u64, m: u64) -> u64 {
+        let Wrapping(msb_mod_m) = Wrapping(MAX_U64) - Wrapping(m) + Wrapping(1);
+        let msb_mod_m = msb_mod_m % m;
+
+        let Wrapping(res) = Wrapping(a) + Wrapping(b);
+        if b <= MAX_U64 - a {
+            res
+        } else {
+            (res + msb_mod_m) % m
         }
-        a = (a << 1) % m;
-        b >>= 1;
     }
-    result
-}
 
-// computes (a + b) % m using the russian peasant algorithm
-// Only necessary when m >= 2^63; otherwise, just wastes time.
-pub fn big_mul(mut a: u64, mut b: u64, m: u64) -> u64 {
-    // precompute 2^64 mod m, since we expect to wrap
-    let Wrapping(msb_mod_m) = Wrapping(MAX_U64) - Wrapping(m) + Wrapping(1);
-    let msb_mod_m = msb_mod_m % m;
+    // computes (a + b) % m using the russian peasant algorithm
+    // Only necessary when m >= 2^63; otherwise, just wastes time.
+    fn mul(mut a: u64, mut b: u64, m: u64) -> u64 {
+        // precompute 2^64 mod m, since we expect to wrap
+        let Wrapping(msb_mod_m) = Wrapping(MAX_U64) - Wrapping(m) + Wrapping(1);
+        let msb_mod_m = msb_mod_m % m;
 
-    let mut result = 0;
-    while b > 0 {
-        if b & 1 != 0 {
-            let Wrapping(next_res) = Wrapping(result) + Wrapping(a);
-            let next_res = next_res % m;
-            result = if result <= MAX_U64 - a {
-                next_res
+        let mut result = 0;
+        while b > 0 {
+            if b & 1 != 0 {
+                let Wrapping(next_res) = Wrapping(result) + Wrapping(a);
+                let next_res = next_res % m;
+                result = if result <= MAX_U64 - a {
+                    next_res
+                } else {
+                    (next_res + msb_mod_m) % m
+                };
+            }
+            let Wrapping(next_a) = Wrapping(a) << 1;
+            let next_a = next_a % m;
+            a = if a < 1 << 63 {
+                next_a
             } else {
-                (next_res + msb_mod_m) % m
+                (next_a + msb_mod_m) % m
             };
+            b >>= 1;
         }
-        let Wrapping(next_a) = Wrapping(a) << 1;
-        let next_a = next_a % m;
-        a = if a < 1 << 63 {
-            next_a
-        } else {
-            (next_a + msb_mod_m) % m
-        };
-        b >>= 1;
+        result
     }
-    result
 }
 
-// computes a.pow(b) % m
-pub(crate) fn pow(mut a: u64, mut b: u64, m: u64, mul: fn(u64, u64, u64) -> u64) -> u64 {
-    let mut result = 1;
-    while b > 0 {
-        if b & 1 != 0 {
-            result = mul(result, a, m);
+pub(crate) struct Small {}
+
+impl Arithmetic for Small {
+    // computes (a + b) % m using the russian peasant algorithm
+    // CAUTION: Will overflow if m >= 2^63
+    fn mul(mut a: u64, mut b: u64, m: u64) -> u64 {
+        let mut result = 0;
+        while b > 0 {
+            if b & 1 != 0 {
+                result = (result + a) % m;
+            }
+            a = (a << 1) % m;
+            b >>= 1;
         }
-        a = mul(a, a, m);
-        b >>= 1;
+        result
+    }
+
+    fn add(a: u64, b: u64, m: u64) -> u64 {
+        (a + b) % m
     }
-    result
 }
diff --git a/src/uu/factor/src/rho.rs b/src/uu/factor/src/rho.rs
index acd07e2da7..f8adb41936 100644
--- a/src/uu/factor/src/rho.rs
+++ b/src/uu/factor/src/rho.rs
@@ -6,36 +6,31 @@ use rand::rngs::SmallRng;
 use rand::{thread_rng, SeedableRng};
 use std::cmp::{max, min};
 
-fn pseudorandom_function(x: u64, a: u64, b: u64, num: u64) -> u64 {
-    if num < 1 << 63 {
-        (sm_mul(a, sm_mul(x, x, num), num) + b) % num
-    } else {
-        big_add(big_mul(a, big_mul(x, x, num), num), b, num)
-    }
-}
-
-fn find_divisor(num: u64) -> u64 {
+fn find_divisor<A: Arithmetic>(n: u64) -> u64 {
     #![allow(clippy::many_single_char_names)]
-    let range = Uniform::new(1, num);
-    let mut rng = SmallRng::from_rng(&mut thread_rng()).unwrap();
-    let mut x = range.sample(&mut rng);
-    let mut y = x;
-    let mut a = range.sample(&mut rng);
-    let mut b = range.sample(&mut rng);
+    let mut rand = {
+        let range = Uniform::new(1, n);
+        let mut rng = SmallRng::from_rng(&mut thread_rng()).unwrap();
+        move || range.sample(&mut rng)
+    };
+
+    let quadratic = |a, b| move |x| A::add(A::mul(a, A::mul(x, x, n), n), b, n);
 
     loop {
-        x = pseudorandom_function(x, a, b, num);
-        y = pseudorandom_function(y, a, b, num);
-        y = pseudorandom_function(y, a, b, num);
-        let d = gcd(num, max(x, y) - min(x, y));
-        if d == num {
-            // Failure, retry with different function
-            x = range.sample(&mut rng);
-            y = x;
-            a = range.sample(&mut rng);
-            b = range.sample(&mut rng);
-        } else if d > 1 {
-            return d;
+        let f = quadratic(rand(), rand());
+        let mut x = rand();
+        let mut y = x;
+
+        loop {
+            x = f(x);
+            y = f(f(y));
+            let d = gcd(n, max(x, y) - min(x, y));
+            if d == n {
+                // Failure, retry with a different quadratic
+                break;
+            } else if d > 1 {
+                return d;
+            }
         }
     }
 }
@@ -46,7 +41,11 @@ pub(crate) fn factor(mut num: u64) -> Factors {
         return factors;
     }
 
-    match miller_rabin::test(num) {
+    match if num < 1 << 63 {
+        miller_rabin::test::<Small>(num)
+    } else {
+        miller_rabin::test::<Big>(num)
+    } {
         Prime => {
             factors.push(num);
             return factors;
@@ -60,7 +59,11 @@ pub(crate) fn factor(mut num: u64) -> Factors {
         Pseudoprime => {}
     };
 
-    let divisor = find_divisor(num);
+    let divisor = if num < 1 << 63 {
+        find_divisor::<Small>(num)
+    } else {
+        find_divisor::<Big>(num)
+    };
     factors *= factor(divisor);
     factors *= factor(num / divisor);
     factors

From 543c7b941a37f3d855458121d5b7353c7d6a7cde Mon Sep 17 00:00:00 2001
From: nicoo <nicoo@debian.org>
Date: Sun, 24 May 2020 18:26:26 +0200
Subject: [PATCH 13/15] factor::rho: Small refactor

---
 src/uu/factor/src/rho.rs | 37 +++++++++++++++++++++++--------------
 1 file changed, 23 insertions(+), 14 deletions(-)

diff --git a/src/uu/factor/src/rho.rs b/src/uu/factor/src/rho.rs
index f8adb41936..e864519c9c 100644
--- a/src/uu/factor/src/rho.rs
+++ b/src/uu/factor/src/rho.rs
@@ -35,17 +35,22 @@ fn find_divisor<A: Arithmetic>(n: u64) -> u64 {
     }
 }
 
-pub(crate) fn factor(mut num: u64) -> Factors {
+fn _factor<A: Arithmetic>(mut num: u64) -> Factors {
+    // Shadow the name, so the recursion automatically goes from “Big” arithmetic to small.
+    let _factor = |n| {
+        if n < 1 << 63 {
+            _factor::<Small>(n)
+        } else {
+            _factor::<A>(n)
+        }
+    };
+
     let mut factors = Factors::new();
     if num == 1 {
         return factors;
     }
 
-    match if num < 1 << 63 {
-        miller_rabin::test::<Small>(num)
-    } else {
-        miller_rabin::test::<Big>(num)
-    } {
+    match miller_rabin::test::<A>(num) {
         Prime => {
             factors.push(num);
             return factors;
@@ -53,18 +58,22 @@ pub(crate) fn factor(mut num: u64) -> Factors {
 
         Composite(d) => {
             num /= d;
-            factors *= factor(d);
+            factors *= _factor(d)
         }
 
         Pseudoprime => {}
     };
 
-    let divisor = if num < 1 << 63 {
-        find_divisor::<Small>(num)
-    } else {
-        find_divisor::<Big>(num)
-    };
-    factors *= factor(divisor);
-    factors *= factor(num / divisor);
+    let divisor = find_divisor::<A>(num);
+    factors *= _factor(divisor);
+    factors *= _factor(num / divisor);
     factors
 }
+
+pub(crate) fn factor(n: u64) -> Factors {
+    if n < 1 << 63 {
+        _factor::<Small>(n)
+    } else {
+        _factor::<Big>(n)
+    }
+}

From 36a29489595210c1fbba91c4110f2ab859168bc3 Mon Sep 17 00:00:00 2001
From: nicoo <nicoo@debian.org>
Date: Sun, 24 May 2020 19:10:34 +0200
Subject: [PATCH 14/15] factor::miller_rabin: Avoid unecessary exponentiation
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

Instead of computing a^r and a^(n-1) = a^(r 2ⁱ) separately,
compute the latter by repeatedly squaring the former.

33.6% performance improvement
---
 src/uu/factor/src/miller_rabin.rs | 25 +++++++++++++++++++------
 1 file changed, 19 insertions(+), 6 deletions(-)

diff --git a/src/uu/factor/src/miller_rabin.rs b/src/uu/factor/src/miller_rabin.rs
index cdc1722fde..f8ad493dd1 100644
--- a/src/uu/factor/src/miller_rabin.rs
+++ b/src/uu/factor/src/miller_rabin.rs
@@ -20,6 +20,7 @@ impl Result {
 
 // Deterministic Miller-Rabin primality-checking algorithm, adapted to extract
 // (some) dividers; it will fail to factor strong pseudoprimes.
+#[allow(clippy::many_single_char_names)]
 pub(crate) fn test<A: Arithmetic>(n: u64) -> Result {
     use self::Result::*;
 
@@ -30,18 +31,30 @@ pub(crate) fn test<A: Arithmetic>(n: u64) -> Result {
         return if n == 2 { Prime } else { Composite(2) };
     }
 
-    let r = (n - 1) >> (n - 1).trailing_zeros();
+    // n-1 = r 2ⁱ
+    let i = (n - 1).trailing_zeros();
+    let r = (n - 1) >> i;
 
     for a in BASIS.iter() {
-        let mut x = a % n;
-        if x == 0 {
+        let a = a % n;
+        if a == 0 {
             break;
         }
 
-        if A::pow(x, n - 1, n) != 1 {
-            return Pseudoprime;
+        // x = a^r mod n
+        let mut x = A::pow(a, r, n);
+
+        {
+            // y = ((x²)²...)² i times = x ^ (2ⁱ) = a ^ (r 2ⁱ) = x ^ (n - 1)
+            let mut y = x;
+            for _ in 0..i {
+                y = A::mul(y, y, n)
+            }
+            if y != 1 {
+                return Pseudoprime;
+            };
         }
-        x = A::pow(x, r, n);
+
         if x == 1 || x == n - 1 {
             break;
         }

From 4c3682aec776b396c584d1d782d28e97d8eb486e Mon Sep 17 00:00:00 2001
From: Nicolas Braud-Santoni <nicoo@debian.org>
Date: Sun, 24 May 2020 19:14:37 +0200
Subject: [PATCH 15/15] factor::Factors::add: Split up to work without NLL

Co-authored-by: Roy Ivy III <rivy.dev@gmail.com>
---
 src/uu/factor/src/factor.rs | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/src/uu/factor/src/factor.rs b/src/uu/factor/src/factor.rs
index 493d60912e..c3c8dd6746 100644
--- a/src/uu/factor/src/factor.rs
+++ b/src/uu/factor/src/factor.rs
@@ -44,7 +44,8 @@ impl Factors {
 
     fn add(&mut self, prime: u64, exp: u8) {
         assert!(exp > 0);
-        self.f.insert(prime, exp + self.f.get(&prime).unwrap_or(&0));
+        let n = *self.f.get(&prime).unwrap_or(&0);
+        self.f.insert(prime, exp + n);
     }
 
     fn push(&mut self, prime: u64) {