Merge pull request #1525 from nbraud/factor/faster

Performance improvements for `factor`

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"]

src/uu/factor/src/factor.rs

@@ -18,139 +18,95 @@ 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 miller_rabin;
 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 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)
+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);
+        let n = *self.f.get(&prime).unwrap_or(&0);
+        self.f.insert(prime, exp + n);
+    }
+
+    fn push(&mut self, prime: u64) {
+        self.add(prime, 1)
     }
 }
 
-fn gcd(mut a: u64, mut b: u64) -> u64 {
-    while b > 0 {
-        a %= b;
-        swap(&mut a, &mut b);
-    }
-    a
-}
-
-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;
+impl ops::MulAssign<Factors> for Factors {
+    fn mul_assign(&mut self, other: Factors) {
+        for (prime, exp) in &other.f {
+            self.add(*prime, *exp);
         }
     }
 }
 
-fn rho_pollard_factor(num: u64, factors: &mut Vec<u64>) {
-    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);
-}
+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();
 
-fn table_division(mut num: u64, factors: &mut Vec<u64>) {
-    if num < 2 {
-        return;
-    }
-    while num % 2 == 0 {
-        num /= 2;
-        factors.push(2);
-    }
-    if num == 1 {
-        return;
-    }
-    if is_prime(num) {
-        factors.push(num);
-        return;
-    }
-    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;
-                }
-            } else {
-                break;
+        for p in primes {
+            for _ in 0..self.f[&p] {
+                write!(f, " {}", p)?
             }
         }
+
+        Ok(())
+    }
+}
+
+fn factor(mut n: u64) -> Factors {
+    let mut factors = Factors::new();
+
+    if n < 2 {
+        factors.push(n);
+        return factors;
     }
 
-    // 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, factors);
-    //}
+    let z = n.trailing_zeros();
+    if z > 0 {
+        factors.add(2, z as u8);
+        n >>= z;
+    }
+
+    if n == 1 {
+        return factors;
+    }
+
+    let (f, n) = table::factor(n);
+    factors *= f;
+
+    if n >= table::NEXT_PRIME {
+        factors *= rho::factor(n);
+    }
+
+    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, factor(num));
     println!();
 }
 

src/uu/factor/src/miller_rabin.rs

@@ -0,0 +1,88 @@
+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)]
+pub(crate) enum Result {
+    Prime,
+    Pseudoprime,
+    Composite(u64),
+}
+
+impl Result {
+    pub(crate) fn is_prime(&self) -> bool {
+        *self == Result::Prime
+    }
+}
+
+// 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::*;
+
+    if n < 2 {
+        return Pseudoprime;
+    }
+    if n % 2 == 0 {
+        return if n == 2 { Prime } else { Composite(2) };
+    }
+
+    // n-1 = r 2ⁱ
+    let i = (n - 1).trailing_zeros();
+    let r = (n - 1) >> i;
+
+    for a in BASIS.iter() {
+        let a = a % n;
+        if a == 0 {
+            break;
+        }
+
+        // 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;
+            };
+        }
+
+        if x == 1 || x == n - 1 {
+            break;
+        }
+
+        loop {
+            let y = A::mul(x, x, n);
+            if y == 1 {
+                return Composite(gcd(x - 1, n));
+            }
+            if y == n - 1 {
+                // This basis element is not a witness of `n` being composite.
+                // Keep looking.
+                break;
+            }
+            x = y;
+        }
+    }
+
+    Prime
+}
+
+// Used by build.rs' tests
+#[allow(dead_code)]
+pub(crate) fn is_prime(n: u64) -> bool {
+    if n < 1 << 63 {
+        test::<Small>(n)
+    } else {
+        test::<Big>(n)
+    }
+    .is_prime()
+}

src/uu/factor/src/numeric.rs

@@ -9,124 +9,99 @@
 * that was distributed with this source code.
 */
 
+use std::mem::swap;
 use std::num::Wrapping;
 use std::u64::MAX as MAX_U64;
 
-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 fn gcd(mut a: u64, mut b: u64) -> u64 {
+    while b > 0 {
+        a %= b;
+        swap(&mut a, &mut b);
+    }
+    a
+}
 
-    let Wrapping(res) = Wrapping(a) + Wrapping(b);
-    if b <= MAX_U64 - a {
-        res
-    } else {
-        (res + 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;
+
+    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;
-        }
-        a = (a << 1) % m;
-        b >>= 1;
-    }
-    result
-}
+pub(crate) struct Big {}
 
-// 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;
+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 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
+        let Wrapping(res) = Wrapping(a) + Wrapping(b);
+        if b <= MAX_U64 - a {
+            res
         } else {
-            (next_a + msb_mod_m) % m
-        };
-        b >>= 1;
-    }
-    result
-}
-
-// computes a.pow(b) % m
-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);
+            (res + msb_mod_m) % m
         }
-        a = mul(a, a, m);
-        b >>= 1;
-    }
-    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
-    };
+    // 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;
 
-    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;
+        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_a + msb_mod_m) % m
+            };
+            b >>= 1;
         }
-        if a == m - 1 {
-            return false;
+        result
+    }
+}
+
+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);
-    }
-}
-
-// 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;
+        result
     }
 
-    // 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))
+    fn add(a: u64, b: u64, m: u64) -> u64 {
+        (a + b) % m
+    }
 }

src/uu/factor/src/rho.rs

@@ -0,0 +1,79 @@
+use crate::miller_rabin::Result::*;
+use crate::{miller_rabin, Factors};
+use numeric::*;
+use rand::distributions::{Distribution, Uniform};
+use rand::rngs::SmallRng;
+use rand::{thread_rng, SeedableRng};
+use std::cmp::{max, min};
+
+fn find_divisor<A: Arithmetic>(n: u64) -> u64 {
+    #![allow(clippy::many_single_char_names)]
+    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 {
+        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;
+            }
+        }
+    }
+}
+
+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 miller_rabin::test::<A>(num) {
+        Prime => {
+            factors.push(num);
+            return factors;
+        }
+
+        Composite(d) => {
+            num /= d;
+            factors *= _factor(d)
+        }
+
+        Pseudoprime => {}
+    };
+
+    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)
+    }
+}

src/uu/factor/src/table.rs

@@ -0,0 +1,32 @@
+use crate::Factors;
+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 https://math.stackexchange.com/questions/1251327/
+        // for a nice explanation.
+        loop {
+            let Wrapping(x) = Wrapping(num) * Wrapping(inv);
+
+            // While prime divides num
+            if x <= ceil {
+                num = x;
+                factors.push(prime);
+            } else {
+                break;
+            }
+        }
+    }
+
+    (factors, num)
+}