From 43ee92c40f57caed18838ebb5d6bfc9b67578eb1 Mon Sep 17 00:00:00 2001 From: nicoo Date: Tue, 23 Jun 2020 12:09:07 +0200 Subject: [PATCH] factor::numeric: Generalise modular inverse computation --- Cargo.lock | 3 ++ src/uu/factor/Cargo.toml | 4 +++ src/uu/factor/build.rs | 4 +-- src/uu/factor/src/numeric.rs | 57 ++++++++++++++++++++++++++---------- 4 files changed, 50 insertions(+), 18 deletions(-) diff --git a/Cargo.lock b/Cargo.lock index 25997280c..c6dc3bbfa 100644 --- a/Cargo.lock +++ b/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)", diff --git a/src/uu/factor/Cargo.toml b/src/uu/factor/Cargo.toml index 7589fe1bb..5a3cc6077 100644 --- a/src/uu/factor/Cargo.toml +++ b/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" } diff --git a/src/uu/factor/build.rs b/src/uu/factor/build.rs index 719ca63bb..642b646c6 100644 --- a/src/uu/factor/build.rs +++ b/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(); diff --git a/src/uu/factor/src/numeric.rs b/src/uu/factor/src/numeric.rs index 29babb281..f1b447583 100644 --- a/src/uu/factor/src/numeric.rs +++ b/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(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(".wrapping_mul(&newt)); t = newt; newt = newtp; - let newrp = r.wrapping_sub(quot.wrapping_mul(newr)); + let newrp = r.wrapping_sub(".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() { // 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::() + } + + #[test] + fn test_inverter_u64() { + test_inverter::() } #[test]