Revert #1571 “perf/factor ~ deduplicate divisors”

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<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::*;
 
-    let n = A::new(num);
-    match miller_rabin::test::<A>(n) {
-        Prime => None,
-        Composite(d) => Some(d),
-        Pseudoprime => Some(rho::find_divisor::<A>(n)),
-    }
-}
+    // 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)
+        }
+    };
 
-fn find_factor(num: u64) -> Option<u64> {
-    if num < (1 << 32) {
-        _find_factor::<Montgomery<u32>>(num)
-    } else {
-        _find_factor::<Montgomery<u64>>(num)
+    if num == 1 {
+        return f;
     }
+
+    let n = A::new(num);
+    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)
 }
 
-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::<Montgomery<u32>>(n, factors)
+    } else {
+        _factor::<Montgomery<u64>>(n, factors)
     }
-
-    factors
 }
 
 #[cfg(test)]

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