1
Fork 0
mirror of https://github.com/RGBCube/uutils-coreutils synced 2025-07-29 03:57:44 +00:00

factor::numeric: Replace lose functions with an Arithmetic trait

This commit is contained in:
nicoo 2020-05-24 18:16:21 +02:00
parent 29eb8fd77b
commit 30fd6a0309
3 changed files with 113 additions and 95 deletions

View file

@ -20,7 +20,7 @@ impl Result {
// Deterministic Miller-Rabin primality-checking algorithm, adapted to extract // Deterministic Miller-Rabin primality-checking algorithm, adapted to extract
// (some) dividers; it will fail to factor strong pseudoprimes. // (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::*; use self::Result::*;
if n < 2 { if n < 2 {
@ -32,28 +32,22 @@ pub(crate) fn test(n: u64) -> Result {
let r = (n - 1) >> (n - 1).trailing_zeros(); 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() { for a in BASIS.iter() {
let mut x = a % n; let mut x = a % n;
if x == 0 { if x == 0 {
break; break;
} }
if pow(x, n - 1, n, mul) != 1 { if A::pow(x, n - 1, n) != 1 {
return Pseudoprime; return Pseudoprime;
} }
x = pow(x, r, n, mul); x = A::pow(x, r, n);
if x == 1 || x == n - 1 { if x == 1 || x == n - 1 {
break; break;
} }
loop { loop {
let y = mul(x, x, n); let y = A::mul(x, x, n);
if y == 1 { if y == 1 {
return Composite(gcd(x - 1, n)); return Composite(gcd(x - 1, n));
} }
@ -72,5 +66,10 @@ pub(crate) fn test(n: u64) -> Result {
// Used by build.rs' tests // Used by build.rs' tests
#[allow(dead_code)] #[allow(dead_code)]
pub(crate) fn is_prime(n: u64) -> bool { 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()
} }

View file

@ -21,71 +21,87 @@ pub fn gcd(mut a: u64, mut b: u64) -> u64 {
a a
} }
pub fn big_add(a: u64, b: u64, m: u64) -> u64 { pub(crate) trait Arithmetic {
let Wrapping(msb_mod_m) = Wrapping(MAX_U64) - Wrapping(m) + Wrapping(1); fn add(a: u64, b: u64, modulus: u64) -> u64;
let msb_mod_m = msb_mod_m % m; fn mul(a: u64, b: u64, modulus: u64) -> u64;
let Wrapping(res) = Wrapping(a) + Wrapping(b); fn pow(mut a: u64, mut b: u64, m: u64) -> u64 {
if b <= MAX_U64 - a { let mut result = 1;
res while b > 0 {
} else { if b & 1 != 0 {
(res + msb_mod_m) % m result = Self::mul(result, a, m);
}
a = Self::mul(a, a, m);
b >>= 1;
}
result
} }
} }
// computes (a + b) % m using the russian peasant algorithm pub(crate) struct Big {}
// 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
}
// computes (a + b) % m using the russian peasant algorithm impl Arithmetic for Big {
// Only necessary when m >= 2^63; otherwise, just wastes time. fn add(a: u64, b: u64, m: u64) -> u64 {
pub fn big_mul(mut a: u64, mut b: u64, m: u64) -> u64 { let Wrapping(msb_mod_m) = Wrapping(MAX_U64) - Wrapping(m) + Wrapping(1);
// precompute 2^64 mod m, since we expect to wrap let msb_mod_m = msb_mod_m % m;
let Wrapping(msb_mod_m) = Wrapping(MAX_U64) - Wrapping(m) + Wrapping(1);
let msb_mod_m = msb_mod_m % m;
let mut result = 0; let Wrapping(res) = Wrapping(a) + Wrapping(b);
while b > 0 { if b <= MAX_U64 - a {
if b & 1 != 0 { res
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 { } else {
(next_a + msb_mod_m) % m (res + msb_mod_m) % m
}; }
b >>= 1; }
// 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
} 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;
}
result
} }
result
} }
// computes a.pow(b) % m pub(crate) struct Small {}
pub(crate) fn pow(mut a: u64, mut b: u64, m: u64, mul: fn(u64, u64, u64) -> u64) -> u64 {
let mut result = 1; impl Arithmetic for Small {
while b > 0 { // computes (a + b) % m using the russian peasant algorithm
if b & 1 != 0 { // CAUTION: Will overflow if m >= 2^63
result = mul(result, a, m); 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); result
b >>= 1; }
fn add(a: u64, b: u64, m: u64) -> u64 {
(a + b) % m
} }
result
} }

View file

@ -6,36 +6,31 @@ use rand::rngs::SmallRng;
use rand::{thread_rng, SeedableRng}; use rand::{thread_rng, SeedableRng};
use std::cmp::{max, min}; use std::cmp::{max, min};
fn pseudorandom_function(x: u64, a: u64, b: u64, num: u64) -> u64 { fn find_divisor<A: Arithmetic>(n: 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)] #![allow(clippy::many_single_char_names)]
let range = Uniform::new(1, num); let mut rand = {
let mut rng = SmallRng::from_rng(&mut thread_rng()).unwrap(); let range = Uniform::new(1, n);
let mut x = range.sample(&mut rng); let mut rng = SmallRng::from_rng(&mut thread_rng()).unwrap();
let mut y = x; move || range.sample(&mut rng)
let mut a = range.sample(&mut rng); };
let mut b = range.sample(&mut rng);
let quadratic = |a, b| move |x| A::add(A::mul(a, A::mul(x, x, n), n), b, n);
loop { loop {
x = pseudorandom_function(x, a, b, num); let f = quadratic(rand(), rand());
y = pseudorandom_function(y, a, b, num); let mut x = rand();
y = pseudorandom_function(y, a, b, num); let mut y = x;
let d = gcd(num, max(x, y) - min(x, y));
if d == num { loop {
// Failure, retry with different function x = f(x);
x = range.sample(&mut rng); y = f(f(y));
y = x; let d = gcd(n, max(x, y) - min(x, y));
a = range.sample(&mut rng); if d == n {
b = range.sample(&mut rng); // Failure, retry with a different quadratic
} else if d > 1 { break;
return d; } else if d > 1 {
return d;
}
} }
} }
} }
@ -46,7 +41,11 @@ pub(crate) fn factor(mut num: u64) -> Factors {
return 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 => { Prime => {
factors.push(num); factors.push(num);
return factors; return factors;
@ -60,7 +59,11 @@ pub(crate) fn factor(mut num: u64) -> Factors {
Pseudoprime => {} 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(divisor);
factors *= factor(num / divisor); factors *= factor(num / divisor);
factors factors