diff --git a/src/uu/factor/src/factor.rs b/src/uu/factor/src/factor.rs index 334fe644a..f01cae76d 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, } @@ -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 474c80613..d9e6f9441 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 000000000..949497737 --- /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 000000000..d8cfa1236 --- /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) + //} +}