Merge pull request #1549 from nbraud/factor/strong-pseudoprime-bug

factor: Fix very-rare bug in ρ

src/uu/factor/src/factor.rs

@@ -156,4 +156,17 @@ mod tests {
             .map(|i| 2 * i + 2u64.pow(32) + 1)
             .all(|i| factor(i).product() == i));
     }
+
+    #[test]
+    fn factor_recombines_strong_pseudoprime() {
+        // This is a strong pseudoprime (wrt. miller_rabin::BASIS)
+        //  and triggered a bug in rho::factor's codepath handling
+        //  miller_rabbin::Result::Composite
+        let pseudoprime = 17179869183;
+        for _ in 0..20 {
+            // Repeat the test 20 times, as it only fails some fraction
+            // of the time.
+            assert!(factor(pseudoprime).product() == pseudoprime);
+        }
+    }
 }

src/uu/factor/src/numeric.rs

@@ -175,7 +175,7 @@ impl Arithmetic for Montgomery {
 // extended Euclid algorithm
 // precondition: a is odd
 pub(crate) fn inv_mod_u64(a: u64) -> u64 {
-    assert!(a % 2 == 1);
+    assert!(a % 2 == 1, "{} is not odd", a);
     let mut t = 0u64;
     let mut newt = 1u64;
     let mut r = 0u64;

src/uu/factor/src/rho.rs

@@ -48,7 +48,7 @@ fn find_divisor<A: Arithmetic>(n: A) -> u64 {
     }
 }
 
-fn _factor<A: Arithmetic>(mut num: u64) -> Factors {
+fn _factor<A: Arithmetic>(num: u64) -> Factors {
     // Shadow the name, so the recursion automatically goes from “Big” arithmetic to small.
     let _factor = |n| {
         // TODO: Optimise with 32 and 64b versions
@@ -61,6 +61,7 @@ fn _factor<A: Arithmetic>(mut num: u64) -> Factors {
     }
 
     let n = A::new(num);
+    let divisor;
     match miller_rabin::test::<A>(n) {
         Prime => {
             factors.push(num);
@@ -68,14 +69,14 @@ fn _factor<A: Arithmetic>(mut num: u64) -> Factors {
         }
 
         Composite(d) => {
-            num /= d;
+            divisor = d;
-            factors *= _factor(d)
         }
 
-        Pseudoprime => {}
+        Pseudoprime => {
+            divisor = find_divisor::<A>(n);
+        }
     };
 
-    let divisor = find_divisor::<A>(n);
     factors *= _factor(divisor);
     factors *= _factor(num / divisor);
     factors