factor::numeric: Generalise modular inverse computation

Cargo.lock

@@ -1,3 +1,5 @@
+# This file is automatically @generated by Cargo.
+# It is not intended for manual editing.
 [[package]]
 name = "advapi32-sys"
 version = "0.2.0"
@@ -1313,6 +1315,7 @@ dependencies = [
 name = "uu_factor"
 version = "0.0.1"
 dependencies = [
+ "num-traits 0.2.12 (registry+https://github.com/rust-lang/crates.io-index)",
  "rand 0.5.6 (registry+https://github.com/rust-lang/crates.io-index)",
  "uucore 0.0.4 (git+https://github.com/uutils/uucore.git?branch=canary)",
  "uucore_procs 0.0.4 (git+https://github.com/uutils/uucore.git?branch=canary)",

src/uu/factor/Cargo.toml

@@ -14,7 +14,11 @@ edition = "2018"
 [lib]
 path = "src/factor.rs"
 
+[build-dependencies]
+num-traits="0.2"
+
 [dependencies]
+num-traits = "0.2"
 rand = "0.5"
 uucore = { version="0.0.4", package="uucore", git="https://github.com/uutils/uucore.git", branch="canary" }
 uucore_procs = { version="0.0.4", package="uucore_procs", git="https://github.com/uutils/uucore.git", branch="canary" }

src/uu/factor/build.rs

@@ -29,7 +29,7 @@ use miller_rabin::is_prime;
 
 #[path = "src/numeric.rs"]
 mod numeric;
-use numeric::inv_mod_u64;
+use numeric::modular_inverse;
 
 mod sieve;
 
@@ -57,7 +57,7 @@ fn main() {
     let mut x = primes.next().unwrap();
     for next in primes {
         // format the table
-        let outstr = format!("({}, {}, {}),", x, inv_mod_u64(x), std::u64::MAX / x);
+        let outstr = format!("({}, {}, {}),", x, modular_inverse(x), std::u64::MAX / x);
         if cols + outstr.len() > MAX_WIDTH {
             write!(file, "\n    {}", outstr).unwrap();
             cols = 4 + outstr.len();

src/uu/factor/src/numeric.rs

@@ -6,6 +6,11 @@
 // * For the full copyright and license information, please view the LICENSE file
 // * that was distributed with this source code.
 
+use num_traits::{
+    int::PrimInt,
+    ops::wrapping::{WrappingMul, WrappingSub},
+};
+use std::fmt::Debug;
 use std::mem::swap;
 
 // This is incorrectly reported as dead code,
@@ -100,7 +105,7 @@ impl Arithmetic for Montgomery {
     type I = u64;
 
     fn new(n: u64) -> Self {
-        let a = inv_mod_u64(n).wrapping_neg();
+        let a = modular_inverse(n).wrapping_neg();
         debug_assert_eq!(n.wrapping_mul(a), 1_u64.wrapping_neg());
         Montgomery { a, n }
     }
@@ -172,35 +177,42 @@ impl Arithmetic for Montgomery {
     }
 }
 
+pub(crate) trait Int: Debug + PrimInt + WrappingSub + WrappingMul {}
+impl Int for u64 {}
+impl Int for u32 {}
+
 // extended Euclid algorithm
 // precondition: a is odd
-pub(crate) fn inv_mod_u64(a: u64) -> u64 {
-    assert!(a % 2 == 1, "{} is not odd", a);
-    let mut t = 0u64;
-    let mut newt = 1u64;
-    let mut r = 0u64;
+pub(crate) fn modular_inverse<T: Int>(a: T) -> T {
+    let zero = T::zero();
+    let one = T::one();
+    assert!(a % (one + one) == one, "{:?} is not odd", a);
+
+    let mut t = zero;
+    let mut newt = one;
+    let mut r = zero;
     let mut newr = a;
 
-    while newr != 0 {
-        let quot = if r == 0 {
+    while newr != zero {
+        let quot = if r == zero {
             // special case when we're just starting out
             // This works because we know that
             // a does not divide 2^64, so floor(2^64 / a) == floor((2^64-1) / a);
-            std::u64::MAX
+            T::max_value()
         } else {
             r
         } / newr;
 
-        let newtp = t.wrapping_sub(quot.wrapping_mul(newt));
+        let newtp = t.wrapping_sub(&quot.wrapping_mul(&newt));
         t = newt;
         newt = newtp;
 
-        let newrp = r.wrapping_sub(quot.wrapping_mul(newr));
+        let newrp = r.wrapping_sub(&quot.wrapping_mul(&newr));
         r = newr;
         newr = newrp;
     }
 
-    assert_eq!(r, 1);
+    assert_eq!(r, one);
     t
 }
 
@@ -208,12 +220,25 @@ pub(crate) fn inv_mod_u64(a: u64) -> u64 {
 mod tests {
     use super::*;
 
-    #[test]
-    fn test_inverter() {
+    fn test_inverter<T: Int>() {
         // All odd integers from 1 to 20 000
-        let mut test_values = (0..10_000u64).map(|i| 2 * i + 1);
+        let one = T::from(1).unwrap();
+        let two = T::from(2).unwrap();
+        let mut test_values = (0..10_000)
+            .map(|i| T::from(i).unwrap())
+            .map(|i| two * i + one);
 
-        assert!(test_values.all(|x| x.wrapping_mul(inv_mod_u64(x)) == 1));
+        assert!(test_values.all(|x| x.wrapping_mul(&modular_inverse(x)) == one));
+    }
+
+    #[test]
+    fn test_inverter_u32() {
+        test_inverter::<u32>()
+    }
+
+    #[test]
+    fn test_inverter_u64() {
+        test_inverter::<u64>()
     }
 
     #[test]