diff --git a/src/uu/factor/src/miller_rabin.rs b/src/uu/factor/src/miller_rabin.rs
index 9f9dfc913..cdc1722fd 100644
--- a/src/uu/factor/src/miller_rabin.rs
+++ b/src/uu/factor/src/miller_rabin.rs
@@ -20,7 +20,7 @@ impl Result {
 
 // Deterministic Miller-Rabin primality-checking algorithm, adapted to extract
 // (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::*;
 
     if n < 2 {
@@ -32,28 +32,22 @@ pub(crate) fn test(n: u64) -> Result {
 
     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() {
         let mut x = a % n;
         if x == 0 {
             break;
         }
 
-        if pow(x, n - 1, n, mul) != 1 {
+        if A::pow(x, n - 1, n) != 1 {
             return Pseudoprime;
         }
-        x = pow(x, r, n, mul);
+        x = A::pow(x, r, n);
         if x == 1 || x == n - 1 {
             break;
         }
 
         loop {
-            let y = mul(x, x, n);
+            let y = A::mul(x, x, n);
             if y == 1 {
                 return Composite(gcd(x - 1, n));
             }
@@ -72,5 +66,10 @@ pub(crate) fn test(n: u64) -> Result {
 // Used by build.rs' tests
 #[allow(dead_code)]
 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()
 }
diff --git a/src/uu/factor/src/numeric.rs b/src/uu/factor/src/numeric.rs
index 7dfba41fb..95b75e986 100644
--- a/src/uu/factor/src/numeric.rs
+++ b/src/uu/factor/src/numeric.rs
@@ -21,71 +21,87 @@ pub fn gcd(mut a: u64, mut b: u64) -> u64 {
     a
 }
 
-pub fn big_add(a: u64, b: u64, m: u64) -> u64 {
-    let Wrapping(msb_mod_m) = Wrapping(MAX_U64) - Wrapping(m) + Wrapping(1);
-    let msb_mod_m = msb_mod_m % m;
+pub(crate) trait Arithmetic {
+    fn add(a: u64, b: u64, modulus: u64) -> u64;
+    fn mul(a: u64, b: u64, modulus: u64) -> u64;
 
-    let Wrapping(res) = Wrapping(a) + Wrapping(b);
-    if b <= MAX_U64 - a {
-        res
-    } else {
-        (res + msb_mod_m) % m
+    fn pow(mut a: u64, mut b: u64, m: u64) -> u64 {
+        let mut result = 1;
+        while b > 0 {
+            if b & 1 != 0 {
+                result = Self::mul(result, a, m);
+            }
+            a = Self::mul(a, a, m);
+            b >>= 1;
+        }
+        result
     }
 }
 
-// computes (a + b) % m using the russian peasant algorithm
-// 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
-}
+pub(crate) struct Big {}
 
-// computes (a + b) % m using the russian peasant algorithm
-// Only necessary when m >= 2^63; otherwise, just wastes time.
-pub fn big_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;
+impl Arithmetic for Big {
+    fn add(a: u64, b: u64, m: u64) -> u64 {
+        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
+        let Wrapping(res) = Wrapping(a) + Wrapping(b);
+        if b <= MAX_U64 - a {
+            res
         } else {
-            (next_a + msb_mod_m) % m
-        };
-        b >>= 1;
+            (res + msb_mod_m) % m
+        }
+    }
+
+    // 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) fn pow(mut a: u64, mut b: u64, m: u64, mul: fn(u64, u64, u64) -> u64) -> u64 {
-    let mut result = 1;
-    while b > 0 {
-        if b & 1 != 0 {
-            result = mul(result, a, m);
+pub(crate) struct Small {}
+
+impl Arithmetic for Small {
+    // computes (a + b) % m using the russian peasant algorithm
+    // CAUTION: Will overflow if m >= 2^63
+    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);
-        b >>= 1;
+        result
+    }
+
+    fn add(a: u64, b: u64, m: u64) -> u64 {
+        (a + b) % m
     }
-    result
 }
diff --git a/src/uu/factor/src/rho.rs b/src/uu/factor/src/rho.rs
index acd07e2da..f8adb4193 100644
--- a/src/uu/factor/src/rho.rs
+++ b/src/uu/factor/src/rho.rs
@@ -6,36 +6,31 @@ use rand::rngs::SmallRng;
 use rand::{thread_rng, SeedableRng};
 use std::cmp::{max, min};
 
-fn pseudorandom_function(x: u64, a: u64, b: u64, num: 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 {
+fn find_divisor<A: Arithmetic>(n: u64) -> u64 {
     #![allow(clippy::many_single_char_names)]
-    let range = Uniform::new(1, num);
-    let mut rng = SmallRng::from_rng(&mut thread_rng()).unwrap();
-    let mut x = range.sample(&mut rng);
-    let mut y = x;
-    let mut a = range.sample(&mut rng);
-    let mut b = range.sample(&mut rng);
+    let mut rand = {
+        let range = Uniform::new(1, n);
+        let mut rng = SmallRng::from_rng(&mut thread_rng()).unwrap();
+        move || range.sample(&mut rng)
+    };
+
+    let quadratic = |a, b| move |x| A::add(A::mul(a, A::mul(x, x, n), n), b, n);
 
     loop {
-        x = pseudorandom_function(x, a, b, num);
-        y = pseudorandom_function(y, a, b, num);
-        y = pseudorandom_function(y, a, b, num);
-        let d = gcd(num, max(x, y) - min(x, y));
-        if d == num {
-            // Failure, retry with different function
-            x = range.sample(&mut rng);
-            y = x;
-            a = range.sample(&mut rng);
-            b = range.sample(&mut rng);
-        } else if d > 1 {
-            return d;
+        let f = quadratic(rand(), rand());
+        let mut x = rand();
+        let mut y = x;
+
+        loop {
+            x = f(x);
+            y = f(f(y));
+            let d = gcd(n, max(x, y) - min(x, y));
+            if d == n {
+                // Failure, retry with a different quadratic
+                break;
+            } else if d > 1 {
+                return d;
+            }
         }
     }
 }
@@ -46,7 +41,11 @@ pub(crate) fn factor(mut num: u64) -> 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 => {
             factors.push(num);
             return factors;
@@ -60,7 +59,11 @@ pub(crate) fn factor(mut num: u64) -> Factors {
         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(num / divisor);
     factors