1
Fork 0
mirror of https://github.com/RGBCube/uutils-coreutils synced 2025-07-28 11:37:44 +00:00

Revert #1571 “perf/factor ~ deduplicate divisors” (#1842)

It was a draft PR, not ready for merging, and its premature inclusion
caused repeated issues, see 368f47381b & friends.

Close #1841.

This reverts commits 3743a3e1e7,
                     ce218e01b6, and
                     b7b0c76b8e.
This commit is contained in:
nicoo 2021-03-20 11:46:58 +01:00 committed by GitHub
parent e9adc5067b
commit 8b9ac0c7c3
No known key found for this signature in database
GPG key ID: 4AEE18F83AFDEB23
2 changed files with 45 additions and 105 deletions

View file

@ -9,7 +9,7 @@ use smallvec::SmallVec;
use std::cell::RefCell;
use std::fmt;
use crate::numeric::{gcd, Arithmetic, Montgomery};
use crate::numeric::{Arithmetic, Montgomery};
use crate::{miller_rabin, rho, table};
type Exponent = u8;
@ -29,20 +29,15 @@ impl Decomposition {
fn add(&mut self, factor: u64, exp: Exponent) {
debug_assert!(exp > 0);
// Assert the factor doesn't already exist in the Decomposition object
debug_assert_eq!(self.0.iter_mut().find(|(f, _)| *f == factor), None);
self.0.push((factor, exp))
}
fn is_one(&self) -> bool {
self.0.is_empty()
}
fn pop(&mut self) -> Option<(u64, Exponent)> {
self.0.pop()
if let Some((_, e)) = self.0.iter_mut().find(|(f, _)| *f == factor) {
*e += exp;
} else {
self.0.push((factor, exp))
}
}
#[cfg(test)]
fn product(&self) -> u64 {
self.0
.iter()
@ -86,11 +81,11 @@ impl Factors {
self.0.borrow_mut().add(prime, exp)
}
#[cfg(test)]
pub fn push(&mut self, prime: u64) {
self.add(prime, 1)
}
#[cfg(test)]
fn product(&self) -> u64 {
self.0.borrow().product()
}
@ -111,116 +106,62 @@ impl fmt::Display for Factors {
}
}
fn _find_factor<A: Arithmetic + miller_rabin::Basis>(num: u64) -> Option<u64> {
fn _factor<A: Arithmetic + miller_rabin::Basis>(num: u64, f: Factors) -> Factors {
use miller_rabin::Result::*;
// Shadow the name, so the recursion automatically goes from “Big” arithmetic to small.
let _factor = |n, f| {
if n < (1 << 32) {
_factor::<Montgomery<u32>>(n, f)
} else {
_factor::<A>(n, f)
}
};
if num == 1 {
return f;
}
let n = A::new(num);
match miller_rabin::test::<A>(n) {
Prime => None,
Composite(d) => Some(d),
Pseudoprime => Some(rho::find_divisor::<A>(n)),
}
let divisor = match miller_rabin::test::<A>(n) {
Prime => {
let mut r = f;
r.push(num);
return r;
}
Composite(d) => d,
Pseudoprime => rho::find_divisor::<A>(n),
};
let f = _factor(divisor, f);
_factor(num / divisor, f)
}
fn find_factor(num: u64) -> Option<u64> {
if num < (1 << 32) {
_find_factor::<Montgomery<u32>>(num)
} else {
_find_factor::<Montgomery<u64>>(num)
}
}
pub fn factor(num: u64) -> Factors {
pub fn factor(mut n: u64) -> Factors {
let mut factors = Factors::one();
if num < 2 {
if n < 2 {
return factors;
}
let mut n = num;
let n_zeros = num.trailing_zeros();
let n_zeros = n.trailing_zeros();
if n_zeros > 0 {
factors.add(2, n_zeros as Exponent);
n >>= n_zeros;
}
debug_assert_eq!(num, n * factors.product());
if n == 1 {
return factors;
}
table::factor(&mut n, &mut factors);
debug_assert_eq!(num, n * factors.product());
let (factors, n) = table::factor(n, factors);
if n == 1 {
return factors;
if n < (1 << 32) {
_factor::<Montgomery<u32>>(n, factors)
} else {
_factor::<Montgomery<u64>>(n, factors)
}
let mut dec = Decomposition::one();
dec.add(n, 1);
while !dec.is_one() {
// Check correctness invariant
debug_assert_eq!(num, factors.product() * dec.product());
let (factor, exp) = dec.pop().unwrap();
if let Some(divisor) = find_factor(factor) {
let mut gcd_queue = Decomposition::one();
let quotient = factor / divisor;
let mut trivial_gcd = quotient == divisor;
if trivial_gcd {
gcd_queue.add(divisor, exp + 1);
} else {
gcd_queue.add(divisor, exp);
gcd_queue.add(quotient, exp);
}
while !trivial_gcd {
debug_assert_eq!(factor, gcd_queue.product());
let mut tmp = Decomposition::one();
trivial_gcd = true;
for i in 0..gcd_queue.0.len() - 1 {
let (mut a, exp_a) = gcd_queue.0[i];
let (mut b, exp_b) = gcd_queue.0[i + 1];
if a == 1 {
continue;
}
let g = gcd(a, b);
if g != 1 {
trivial_gcd = false;
a /= g;
b /= g;
}
if a != 1 {
tmp.add(a, exp_a);
}
if g != 1 {
tmp.add(g, exp_a + exp_b);
}
if i + 1 != gcd_queue.0.len() - 1 {
gcd_queue.0[i + 1].0 = b;
} else if b != 1 {
tmp.add(b, exp_b);
}
}
gcd_queue = tmp;
}
debug_assert_eq!(factor, gcd_queue.product());
dec.0.extend(gcd_queue.0);
} else {
// factor is prime
factors.add(factor, exp);
}
}
factors
}
#[cfg(test)]

View file

@ -14,8 +14,7 @@ use crate::Factors;
include!(concat!(env!("OUT_DIR"), "/prime_table.rs"));
pub(crate) fn factor(n: &mut u64, factors: &mut Factors) {
let mut num = *n;
pub(crate) fn factor(mut num: u64, mut factors: Factors) -> (Factors, u64) {
for &(prime, inv, ceil) in P_INVS_U64 {
if num == 1 {
break;
@ -43,5 +42,5 @@ pub(crate) fn factor(n: &mut u64, factors: &mut Factors) {
}
}
*n = num;
(factors, num)
}