diff --git a/src/uu/factor/build.rs b/src/uu/factor/build.rs index acc4d9a3f..3a1b7df94 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 81f655eb7..c3c8dd674 100644 --- a/src/uu/factor/src/factor.rs +++ b/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, +} + +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 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) { - 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) { - 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!(); } diff --git a/src/uu/factor/src/miller_rabin.rs b/src/uu/factor/src/miller_rabin.rs new file mode 100644 index 000000000..f8ad493dd --- /dev/null +++ b/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(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::(n) + } else { + test::(n) + } + .is_prime() +} diff --git a/src/uu/factor/src/numeric.rs b/src/uu/factor/src/numeric.rs index 474c80613..95b75e986 100644 --- a/src/uu/factor/src/numeric.rs +++ b/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 + } } diff --git a/src/uu/factor/src/rho.rs b/src/uu/factor/src/rho.rs new file mode 100644 index 000000000..e864519c9 --- /dev/null +++ b/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(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(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::(n) + } else { + _factor::(n) + } + }; + + let mut factors = Factors::new(); + if num == 1 { + return factors; + } + + match miller_rabin::test::(num) { + Prime => { + factors.push(num); + return factors; + } + + Composite(d) => { + num /= d; + factors *= _factor(d) + } + + Pseudoprime => {} + }; + + let divisor = find_divisor::(num); + factors *= _factor(divisor); + factors *= _factor(num / divisor); + factors +} + +pub(crate) fn factor(n: u64) -> Factors { + if n < 1 << 63 { + _factor::(n) + } else { + _factor::(n) + } +} diff --git a/src/uu/factor/src/table.rs b/src/uu/factor/src/table.rs new file mode 100644 index 000000000..3a4f44e85 --- /dev/null +++ b/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) +}