From 8b9ac0c7c3231ebde87b875e95dc67d02e9e753a Mon Sep 17 00:00:00 2001 From: nicoo Date: Sat, 20 Mar 2021 11:46:58 +0100 Subject: [PATCH] =?UTF-8?q?Revert=20#1571=20=E2=80=9Cperf/factor=20~=20ded?= =?UTF-8?q?uplicate=20divisors=E2=80=9D=20(#1842)?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit It was a draft PR, not ready for merging, and its premature inclusion caused repeated issues, see 368f47381b7441ddb23bb7460417e0d1ce33e330 & friends. Close #1841. This reverts commits 3743a3e1e7a87bdd0ea7b923dc42386bfc14ccd5, ce218e01b6fcfadf4a9e55348d98f36dafff6c6b, and b7b0c76b8ea60297b43777d17d7322d1d45a9e91. --- src/uu/factor/src/factor.rs | 143 +++++++++++------------------------- src/uu/factor/src/table.rs | 5 +- 2 files changed, 44 insertions(+), 104 deletions(-) diff --git a/src/uu/factor/src/factor.rs b/src/uu/factor/src/factor.rs index 2715bad71c8..7d2e16a1162 100644 --- a/src/uu/factor/src/factor.rs +++ b/src/uu/factor/src/factor.rs @@ -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(num: u64) -> Option { +fn _factor(num: u64, f: Factors) -> Factors { use miller_rabin::Result::*; - let n = A::new(num); - match miller_rabin::test::(n) { - Prime => None, - Composite(d) => Some(d), - Pseudoprime => Some(rho::find_divisor::(n)), - } -} + // Shadow the name, so the recursion automatically goes from “Big” arithmetic to small. + let _factor = |n, f| { + if n < (1 << 32) { + _factor::>(n, f) + } else { + _factor::(n, f) + } + }; -fn find_factor(num: u64) -> Option { - if num < (1 << 32) { - _find_factor::>(num) - } else { - _find_factor::>(num) + if num == 1 { + return f; } + + let n = A::new(num); + let divisor = match miller_rabin::test::(n) { + Prime => { + let mut r = f; + r.push(num); + return r; + } + + Composite(d) => d, + Pseudoprime => rho::find_divisor::(n), + }; + + let f = _factor(divisor, f); + _factor(num / divisor, f) } -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()); if n == 1 { return 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]; + let (factors, n) = table::factor(n, factors); - 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); - } + if n < (1 << 32) { + _factor::>(n, factors) + } else { + _factor::>(n, factors) } - - factors } #[cfg(test)] diff --git a/src/uu/factor/src/table.rs b/src/uu/factor/src/table.rs index b62e801cbc7..d6ef796fc03 100644 --- a/src/uu/factor/src/table.rs +++ b/src/uu/factor/src/table.rs @@ -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) }