1
Fork 0
mirror of https://github.com/RGBCube/uutils-coreutils synced 2025-08-04 23:17:46 +00:00

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.
This commit is contained in:
nicoo 2020-05-24 15:41:23 +02:00
parent d9095a2539
commit a1b2522750
4 changed files with 122 additions and 108 deletions

View file

@ -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!();
}

View file

@ -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;

54
src/uu/factor/src/rho.rs Normal file
View file

@ -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
}

View file

@ -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)
//}
}