From bada7530fb723829230fab35a839ca8bd213018d Mon Sep 17 00:00:00 2001 From: nicoo Date: Mon, 25 May 2020 16:18:16 +0200 Subject: [PATCH 01/14] factor::miller_rabin: Add tests --- src/uu/factor/src/miller_rabin.rs | 18 ++++++++++++++++++ 1 file changed, 18 insertions(+) diff --git a/src/uu/factor/src/miller_rabin.rs b/src/uu/factor/src/miller_rabin.rs index 63ef70d02..26905b996 100644 --- a/src/uu/factor/src/miller_rabin.rs +++ b/src/uu/factor/src/miller_rabin.rs @@ -88,3 +88,21 @@ pub(crate) fn is_prime(n: u64) -> bool { } .is_prime() } + +#[cfg(test)] +mod tests { + use super::is_prime; + const LARGEST_U64_PRIME: u64 = 0xFFFFFFFFFFFFFFC5; + + #[test] + fn largest_prime() { + assert!(is_prime(LARGEST_U64_PRIME)); + } + + #[test] + fn first_primes() { + use crate::table::{NEXT_PRIME, P_INVS_U64}; + assert!(P_INVS_U64.iter().all(|(p, _, _)| is_prime(*p))); + assert!(is_prime(NEXT_PRIME)); + } +} From e91155519a0b3c0b7e0ba5952a5c2d8091a7ea26 Mon Sep 17 00:00:00 2001 From: nicoo Date: Mon, 25 May 2020 16:50:17 +0200 Subject: [PATCH 02/14] factor::factor: Add integration tests --- src/uu/factor/src/factor.rs | 26 ++++++++++++++++++++++++++ 1 file changed, 26 insertions(+) diff --git a/src/uu/factor/src/factor.rs b/src/uu/factor/src/factor.rs index 2fcf66ad4..bd6c26e4a 100644 --- a/src/uu/factor/src/factor.rs +++ b/src/uu/factor/src/factor.rs @@ -47,6 +47,13 @@ impl Factors { fn push(&mut self, prime: u64) { self.add(prime, 1) } + + #[cfg(test)] + fn product(&self) -> u64 { + self.f + .iter() + .fold(1, |acc, (p, exp)| acc * p.pow(*exp as u32)) + } } impl ops::MulAssign for Factors { @@ -132,3 +139,22 @@ pub fn uumain(args: Vec) -> i32 { } 0 } + +#[cfg(test)] +mod tests { + use super::factor; + + #[test] + fn factor_recombines_small() { + assert!((1..10_000) + .map(|i| 2 * i + 1) + .all(|i| factor(i).product() == i)); + } + + #[test] + fn factor_recombines_overflowing() { + assert!((0..250) + .map(|i| 2 * i + 2u64.pow(32) + 1) + .all(|i| factor(i).product() == i)); + } +} From 8a4d0d30ad2f4553efb426007b853a90bf9bf2cd Mon Sep 17 00:00:00 2001 From: nicoo Date: Mon, 25 May 2020 21:44:54 +0200 Subject: [PATCH 03/14] factor::numeric: Implement Montgomery's transform This is a facter way to perform arithmetic mod n, when n is odd and a 64b number. --- src/uu/factor/build.rs | 53 +------- src/uu/factor/src/miller_rabin.rs | 38 +++--- src/uu/factor/src/numeric.rs | 217 +++++++++++++++++++++--------- src/uu/factor/src/rho.rs | 36 +++-- 4 files changed, 195 insertions(+), 149 deletions(-) diff --git a/src/uu/factor/build.rs b/src/uu/factor/build.rs index 77fa3851a..637c8981f 100644 --- a/src/uu/factor/build.rs +++ b/src/uu/factor/build.rs @@ -20,56 +20,19 @@ use std::env::{self, args}; use std::fs::File; use std::io::Write; -use std::num::Wrapping; use std::path::Path; -use std::u64::MAX as MAX_U64; use self::sieve::Sieve; #[cfg(test)] use miller_rabin::is_prime; -#[cfg(test)] #[path = "src/numeric.rs"] mod numeric; +use numeric::inv_mod_u64; mod sieve; -// extended Euclid algorithm -// precondition: a does not divide 2^64 -fn inv_mod_u64(a: u64) -> Option { - let mut t = 0u64; - let mut newt = 1u64; - let mut r = 0u64; - let mut newr = a; - - while newr != 0 { - let quot = if r == 0 { - // 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); - MAX_U64 - } else { - r - } / newr; - - let (tp, Wrapping(newtp)) = (newt, Wrapping(t) - (Wrapping(quot) * Wrapping(newt))); - t = tp; - newt = newtp; - - let (rp, Wrapping(newrp)) = (newr, Wrapping(r) - (Wrapping(quot) * Wrapping(newr))); - r = rp; - newr = newrp; - } - - if r > 1 { - // not invertible - return None; - } - - Some(t) -} - #[cfg_attr(test, allow(dead_code))] fn main() { let out_dir = env::var("OUT_DIR").unwrap(); @@ -95,7 +58,7 @@ fn main() { let mut x = primes.next().unwrap(); for next in primes { // format the table - let outstr = format!("({}, {}, {}),", x, inv_mod_u64(x).unwrap(), MAX_U64 / x); + let outstr = format!("({}, {}, {}),", x, inv_mod_u64(x), u64::MAX / x); if cols + outstr.len() > MAX_WIDTH { write!(file, "\n {}", outstr).unwrap(); cols = 4 + outstr.len(); @@ -116,18 +79,12 @@ fn main() { } #[test] -fn test_inverter() { - let num = 10000; - - let invs = Sieve::odd_primes().map(|x| inv_mod_u64(x).unwrap()); - assert!(Sieve::odd_primes().zip(invs).take(num).all(|(x, y)| { - let Wrapping(z) = Wrapping(x) * Wrapping(y); - is_prime(x) && z == 1 - })); +fn test_generator_isprime() { + assert_eq!(Sieve::odd_primes.take(10_000).all(is_prime)); } #[test] -fn test_generator() { +fn test_generator_10001() { let prime_10001 = Sieve::primes().skip(10_000).next(); assert_eq!(prime_10001, Some(104_743)); } diff --git a/src/uu/factor/src/miller_rabin.rs b/src/uu/factor/src/miller_rabin.rs index 26905b996..978cd6e47 100644 --- a/src/uu/factor/src/miller_rabin.rs +++ b/src/uu/factor/src/miller_rabin.rs @@ -23,9 +23,10 @@ impl Result { // Deterministic Miller-Rabin primality-checking algorithm, adapted to extract // (some) dividers; it will fail to factor strong pseudoprimes. #[allow(clippy::many_single_char_names)] -pub(crate) fn test(n: u64) -> Result { +pub(crate) fn test(m: A) -> Result { use self::Result::*; + let n = m.modulus(); if n < 2 { return Pseudoprime; } @@ -37,36 +38,38 @@ pub(crate) fn test(n: u64) -> Result { let i = (n - 1).trailing_zeros(); let r = (n - 1) >> i; - for a in BASIS.iter() { - let a = a % n; - if a == 0 { + for _a in BASIS.iter() { + let _a = _a % n; + if _a == 0 { break; } + let a = m.from_u64(_a); + // x = a^r mod n - let mut x = A::pow(a, r, n); + let mut x = m.pow(a, r); { // y = ((x²)²...)² i times = x ^ (2ⁱ) = a ^ (r 2ⁱ) = x ^ (n - 1) let mut y = x; for _ in 0..i { - y = A::mul(y, y, n) + y = m.mul(y, y) } - if y != 1 { + if y != m.one() { return Pseudoprime; }; } - if x == 1 || x == n - 1 { + if x == m.one() || x == m.minus_one() { break; } loop { - let y = A::mul(x, x, n); - if y == 1 { - return Composite(gcd(x - 1, n)); + let y = m.mul(x, x); + if y == m.one() { + return Composite(gcd(m.to_u64(x) - 1, m.modulus())); } - if y == n - 1 { + if y == m.minus_one() { // This basis element is not a witness of `n` being composite. // Keep looking. break; @@ -81,12 +84,7 @@ pub(crate) fn test(n: u64) -> Result { // Used by build.rs' tests #[allow(dead_code)] pub(crate) fn is_prime(n: u64) -> bool { - if n < 1 << 63 { - test::(n) - } else { - test::(n) - } - .is_prime() + test::(Montgomery::new(n)).is_prime() } #[cfg(test)] @@ -102,7 +100,9 @@ mod tests { #[test] fn first_primes() { use crate::table::{NEXT_PRIME, P_INVS_U64}; - assert!(P_INVS_U64.iter().all(|(p, _, _)| is_prime(*p))); + for (p, _, _) in P_INVS_U64.iter() { + assert!(is_prime(*p), "{} reported composite", p); + } assert!(is_prime(NEXT_PRIME)); } } diff --git a/src/uu/factor/src/numeric.rs b/src/uu/factor/src/numeric.rs index a1699951a..d9f0ed7c9 100644 --- a/src/uu/factor/src/numeric.rs +++ b/src/uu/factor/src/numeric.rs @@ -9,7 +9,6 @@ use std::mem::swap; use std::num::Wrapping; -use std::u64::MAX as MAX_U64; pub fn gcd(mut a: u64, mut b: u64) -> u64 { while b > 0 { @@ -19,87 +18,179 @@ pub fn gcd(mut a: u64, mut b: u64) -> u64 { a } -pub(crate) trait Arithmetic { - fn add(a: u64, b: u64, modulus: u64) -> u64; - fn mul(a: u64, b: u64, modulus: u64) -> u64; +pub(crate) trait Arithmetic: Copy + Sized { + type I: Copy + Sized + Eq; - fn pow(mut a: u64, mut b: u64, m: u64) -> u64 { - let mut result = 1; + fn new(m: u64) -> Self; + fn modulus(&self) -> u64; + fn from_u64(&self, n: u64) -> Self::I; + fn to_u64(&self, n: Self::I) -> u64; + fn add(&self, a: Self::I, b: Self::I) -> Self::I; + fn mul(&self, a: Self::I, b: Self::I) -> Self::I; + + fn pow(&self, mut a: Self::I, mut b: u64) -> Self::I { + let mut result = self.from_u64(1u64); while b > 0 { if b & 1 != 0 { - result = Self::mul(result, a, m); + result = self.mul(result, a); } - a = Self::mul(a, a, m); + a = self.mul(a, a); b >>= 1; } result } + + fn one(&self) -> Self::I { + self.from_u64(1) + } + fn minus_one(&self) -> Self::I { + self.from_u64(self.modulus() - 1) + } + fn zero(&self) -> Self::I { + self.from_u64(0) + } } -pub(crate) struct Big {} +#[derive(Clone, Copy, Debug)] +pub(crate) struct Montgomery { + a: u64, + n: u64, +} -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 Wrapping(res) = Wrapping(a) + Wrapping(b); - if b <= MAX_U64 - a { - res +impl Montgomery { + /// computes x/R mod n efficiently + fn reduce(&self, x: u64) -> u64 { + // TODO: optimiiiiiiise + let Montgomery { a, n } = self; + let t = x.wrapping_mul(*a); + let nt = (*n as u128) * (t as u128); + let y = ((x as u128 + nt) >> 64) as u64; + if y >= *n { + y - n } else { - (res + msb_mod_m) % m + y } } - - // 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 - } } -pub(crate) struct Small {} +impl Arithmetic for Montgomery { + // Montgomery transform, R=2⁶⁴ + // Provides fast arithmetic mod n (n odd, u64) + type I = Wrapping; -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; + fn new(n: u64) -> Self { + Montgomery { + a: inv_mod_u64(n).wrapping_neg(), + n, } - result } - fn add(a: u64, b: u64, m: u64) -> u64 { - (a + b) % m + fn modulus(&self) -> u64 { + self.n + } + + fn from_u64(&self, x: u64) -> Self::I { + // TODO: optimise! + Wrapping((((x as u128) << 64) % self.n as u128) as u64) + } + + fn to_u64(&self, n: Self::I) -> u64 { + self.reduce(n.0) + } + + fn add(&self, a: Self::I, b: Self::I) -> Self::I { + a + b + } + + fn mul(&self, a: Self::I, b: Self::I) -> Self::I { + Wrapping(self.reduce((a * b).0)) + } +} + +// extended Euclid algorithm +// precondition: a is odd +pub(crate) fn inv_mod_u64(a: u64) -> u64 { + assert!(a % 2 == 1); + let mut t = 0u64; + let mut newt = 1u64; + let mut r = 0u64; + let mut newr = a; + + while newr != 0 { + let quot = if r == 0 { + // 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); + u64::MAX + } else { + r + } / newr; + + let (tp, Wrapping(newtp)) = (newt, Wrapping(t) - (Wrapping(quot) * Wrapping(newt))); + t = tp; + newt = newtp; + + let (rp, Wrapping(newrp)) = (newr, Wrapping(r) - (Wrapping(quot) * Wrapping(newr))); + r = rp; + newr = newrp; + } + + assert_eq!(r, 1); + t +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn test_inverter() { + // All odd integers from 1 to 20 000 + let mut test_values = (0..10_000u64).map(|i| 2 * i + 1); + + assert!(test_values.all(|x| x.wrapping_mul(inv_mod_u64(x)) == 1)); + } + + #[test] + fn test_montgomery_add() { + for n in 0..100 { + let n = 2 * n + 1; + let m = Montgomery::new(n); + for x in 0..n { + let m_x = m.from_u64(x); + for y in 0..=x { + let m_y = m.from_u64(y); + println!("{n:?}, {x:?}, {y:?}", n = n, x = x, y = y); + assert_eq!((x + y) % n, m.to_u64(m.add(m_x, m_y))); + } + } + } + } + + #[test] + fn test_montgomery_mult() { + for n in 0..100 { + let n = 2 * n + 1; + let m = Montgomery::new(n); + for x in 0..n { + let m_x = m.from_u64(x); + for y in 0..=x { + let m_y = m.from_u64(y); + assert_eq!((x * y) % n, m.to_u64(m.mul(m_x, m_y))); + } + } + } + } + + #[test] + fn test_montgomery_roundtrip() { + for n in 0..100 { + let n = 2 * n + 1; + let m = Montgomery::new(n); + for x in 0..n { + let x_ = m.from_u64(x); + assert_eq!(x, m.to_u64(x_)); + } + } } } diff --git a/src/uu/factor/src/rho.rs b/src/uu/factor/src/rho.rs index 6caded033..a6608747a 100644 --- a/src/uu/factor/src/rho.rs +++ b/src/uu/factor/src/rho.rs @@ -3,19 +3,19 @@ use rand::rngs::SmallRng; use rand::{thread_rng, SeedableRng}; use std::cmp::{max, min}; +use crate::{Factors,miller_rabin}; use crate::miller_rabin::Result::*; use crate::numeric::*; -use crate::{miller_rabin, Factors}; -fn find_divisor(n: u64) -> u64 { +fn find_divisor(n: A) -> u64 { #![allow(clippy::many_single_char_names)] let mut rand = { - let range = Uniform::new(1, n); + let range = Uniform::new(1, n.modulus()); let mut rng = SmallRng::from_rng(&mut thread_rng()).unwrap(); - move || range.sample(&mut rng) + move || n.from_u64(range.sample(&mut rng)) }; - let quadratic = |a, b| move |x| A::add(A::mul(a, A::mul(x, x, n), n), b, n); + let quadratic = |a, b| move |x| n.add(n.mul(a, n.mul(x, x)), b); loop { let f = quadratic(rand(), rand()); @@ -25,8 +25,12 @@ fn find_divisor(n: u64) -> u64 { loop { x = f(x); y = f(f(y)); - let d = gcd(n, max(x, y) - min(x, y)); - if d == n { + let d = { + let _x = n.to_u64(x); + let _y = n.to_u64(y); + gcd(n.modulus(), max(_x, _y) - min(_x, _y)) + }; + if d == n.modulus() { // Failure, retry with a different quadratic break; } else if d > 1 { @@ -39,11 +43,8 @@ fn find_divisor(n: u64) -> u64 { fn _factor(mut num: u64) -> Factors { // Shadow the name, so the recursion automatically goes from “Big” arithmetic to small. let _factor = |n| { - if n < 1 << 63 { - _factor::(n) - } else { - _factor::(n) - } + // TODO: Optimise with 32 and 64b versions + _factor::(n) }; let mut factors = Factors::new(); @@ -51,7 +52,8 @@ fn _factor(mut num: u64) -> Factors { return factors; } - match miller_rabin::test::(num) { + let n = A::new(num); + match miller_rabin::test::(n) { Prime => { factors.push(num); return factors; @@ -65,16 +67,12 @@ fn _factor(mut num: u64) -> Factors { Pseudoprime => {} }; - let divisor = find_divisor::(num); + let divisor = find_divisor::(n); factors *= _factor(divisor); factors *= _factor(num / divisor); factors } pub(crate) fn factor(n: u64) -> Factors { - if n < 1 << 63 { - _factor::(n) - } else { - _factor::(n) - } + _factor::(n) } From 33e18b4cd390311782609f3b46d0741430d7bdb4 Mon Sep 17 00:00:00 2001 From: nicoo Date: Sat, 30 May 2020 10:11:05 +0200 Subject: [PATCH 04/14] factor::numeric::Montgomery: Add debug assertions In debug mode, checks that all arithmetic operations coincide with the plain-u64 versions, as long as the latter does not overflow. --- src/uu/factor/src/numeric.rs | 60 +++++++++++++++++++++++++++++++----- 1 file changed, 52 insertions(+), 8 deletions(-) diff --git a/src/uu/factor/src/numeric.rs b/src/uu/factor/src/numeric.rs index d9f0ed7c9..ce25f81bf 100644 --- a/src/uu/factor/src/numeric.rs +++ b/src/uu/factor/src/numeric.rs @@ -29,7 +29,8 @@ pub(crate) trait Arithmetic: Copy + Sized { fn mul(&self, a: Self::I, b: Self::I) -> Self::I; fn pow(&self, mut a: Self::I, mut b: u64) -> Self::I { - let mut result = self.from_u64(1u64); + let (_a, _b) = (a, b); + let mut result = self.one(); while b > 0 { if b & 1 != 0 { result = self.mul(result, a); @@ -37,6 +38,15 @@ pub(crate) trait Arithmetic: Copy + Sized { a = self.mul(a, a); b >>= 1; } + + // Check that r (reduced back to the usual representation) equals + // a^b % n, unless the latter computation overflows + debug_assert!(self + .to_u64(_a) + .checked_pow(_b as u32) + .map(|r| r % self.modulus() == self.to_u64(result)) + .unwrap_or(true)); + result } @@ -79,10 +89,9 @@ impl Arithmetic for Montgomery { type I = Wrapping; fn new(n: u64) -> Self { - Montgomery { - a: inv_mod_u64(n).wrapping_neg(), - n, - } + let a = inv_mod_u64(n).wrapping_neg(); + debug_assert_eq!(n.wrapping_mul(a), 1_u64.wrapping_neg()); + Montgomery { a, n } } fn modulus(&self) -> u64 { @@ -91,7 +100,10 @@ impl Arithmetic for Montgomery { fn from_u64(&self, x: u64) -> Self::I { // TODO: optimise! - Wrapping((((x as u128) << 64) % self.n as u128) as u64) + assert!(x < self.n); + let r = Wrapping((((x as u128) << 64) % self.n as u128) as u64); + debug_assert_eq!(x, self.to_u64(r)); + r } fn to_u64(&self, n: Self::I) -> u64 { @@ -99,11 +111,43 @@ impl Arithmetic for Montgomery { } fn add(&self, a: Self::I, b: Self::I) -> Self::I { - a + b + let r = a + b; + + // Check that r (reduced back to the usual representation) equals + // a+b % n + #[cfg(debug_assertions)] + { + let a_r = self.to_u64(a); + let b_r = self.to_u64(b); + let r_r = self.to_u64(r); + let r_2 = (((a_r as u128) + (b_r as u128)) % (self.n as u128)) as u64; + debug_assert_eq!( + r_r, r_2, + "[{}] = {} ≠ {} = {} + {} = [{}] + [{}] mod {}; a = {}", + r, r_r, r_2, a_r, b_r, a, b, self.n, self.a + ); + } + r } fn mul(&self, a: Self::I, b: Self::I) -> Self::I { - Wrapping(self.reduce((a * b).0)) + let r = Wrapping(self.reduce((a * b).0)); + + // Check that r (reduced back to the usual representation) equals + // a*b % n + #[cfg(debug_assertions)] + { + let a_r = self.to_u64(a); + let b_r = self.to_u64(b); + let r_r = self.to_u64(r); + let r_2 = (((a_r as u128) * (b_r as u128)) % (self.n as u128)) as u64; + debug_assert_eq!( + r_r, r_2, + "[{}] = {} ≠ {} = {} * {} = [{}] * [{}] mod {}; a = {}", + r, r_r, r_2, a_r, b_r, a, b, self.n, self.a + ); + } + r } } From f84d0f93985a6c64cc3838339c4e1c633d676555 Mon Sep 17 00:00:00 2001 From: nicoo Date: Sat, 30 May 2020 11:32:55 +0200 Subject: [PATCH 05/14] factor::Factors::add: Make the precondition check a debug_assert --- src/uu/factor/src/factor.rs | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/uu/factor/src/factor.rs b/src/uu/factor/src/factor.rs index bd6c26e4a..a4ff91ba8 100644 --- a/src/uu/factor/src/factor.rs +++ b/src/uu/factor/src/factor.rs @@ -39,7 +39,7 @@ impl Factors { } fn add(&mut self, prime: u64, exp: u8) { - assert!(exp > 0); + debug_assert!(exp > 0); let n = *self.f.get(&prime).unwrap_or(&0); self.f.insert(prime, exp + n); } From 918035e01eabc51554a97809f86175a6d7fdb773 Mon Sep 17 00:00:00 2001 From: nicoo Date: Sat, 30 May 2020 14:37:09 +0200 Subject: [PATCH 06/14] factor: Fix for old Rust --- src/uu/factor/build.rs | 2 +- src/uu/factor/src/numeric.rs | 17 ++++++++++------- 2 files changed, 11 insertions(+), 8 deletions(-) diff --git a/src/uu/factor/build.rs b/src/uu/factor/build.rs index 637c8981f..1677e44eb 100644 --- a/src/uu/factor/build.rs +++ b/src/uu/factor/build.rs @@ -58,7 +58,7 @@ fn main() { let mut x = primes.next().unwrap(); for next in primes { // format the table - let outstr = format!("({}, {}, {}),", x, inv_mod_u64(x), u64::MAX / x); + let outstr = format!("({}, {}, {}),", x, inv_mod_u64(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 ce25f81bf..38b98b620 100644 --- a/src/uu/factor/src/numeric.rs +++ b/src/uu/factor/src/numeric.rs @@ -40,12 +40,15 @@ pub(crate) trait Arithmetic: Copy + Sized { } // Check that r (reduced back to the usual representation) equals - // a^b % n, unless the latter computation overflows - debug_assert!(self - .to_u64(_a) - .checked_pow(_b as u32) - .map(|r| r % self.modulus() == self.to_u64(result)) - .unwrap_or(true)); + // a^b % n, unless the latter computation overflows + // Temporarily commented-out, as there u64::checked_pow is not available + // on the minimum supported Rust version, nor is an appropriate method + // for compiling the check conditionally. + //debug_assert!(self + // .to_u64(_a) + // .checked_pow(_b as u32) + // .map(|r| r % self.modulus() == self.to_u64(result)) + // .unwrap_or(true)); result } @@ -165,7 +168,7 @@ pub(crate) fn inv_mod_u64(a: u64) -> u64 { // 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); - u64::MAX + std::u64::MAX } else { r } / newr; From 19a0645a0ae946f482060c0395dc67d22738eb67 Mon Sep 17 00:00:00 2001 From: nicoo Date: Sun, 31 May 2020 13:12:30 +0200 Subject: [PATCH 07/14] factor::numeric: Simplify inv_mod_u64 Just call `u64::wrapping_{mul,sub}` instead of (de)constructing Wrapping values. --- src/uu/factor/src/numeric.rs | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/uu/factor/src/numeric.rs b/src/uu/factor/src/numeric.rs index 38b98b620..351cb8a11 100644 --- a/src/uu/factor/src/numeric.rs +++ b/src/uu/factor/src/numeric.rs @@ -173,12 +173,12 @@ pub(crate) fn inv_mod_u64(a: u64) -> u64 { r } / newr; - let (tp, Wrapping(newtp)) = (newt, Wrapping(t) - (Wrapping(quot) * Wrapping(newt))); - t = tp; + let newtp = t.wrapping_sub(quot.wrapping_mul(newt)); + t = newt; newt = newtp; - let (rp, Wrapping(newrp)) = (newr, Wrapping(r) - (Wrapping(quot) * Wrapping(newr))); - r = rp; + let newrp = r.wrapping_sub(quot.wrapping_mul(newr)); + r = newr; newr = newrp; } From 2238065c9d7095dcb6b193474b74ce33334e3351 Mon Sep 17 00:00:00 2001 From: nicoo Date: Sun, 31 May 2020 13:15:53 +0200 Subject: [PATCH 08/14] factor::numeric: Simplify `Montgomery` (remove superfluous Wrapping) --- src/uu/factor/src/numeric.rs | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/src/uu/factor/src/numeric.rs b/src/uu/factor/src/numeric.rs index 351cb8a11..7c6e7ac3e 100644 --- a/src/uu/factor/src/numeric.rs +++ b/src/uu/factor/src/numeric.rs @@ -8,7 +8,6 @@ // * that was distributed with this source code. use std::mem::swap; -use std::num::Wrapping; pub fn gcd(mut a: u64, mut b: u64) -> u64 { while b > 0 { @@ -89,7 +88,7 @@ impl Montgomery { impl Arithmetic for Montgomery { // Montgomery transform, R=2⁶⁴ // Provides fast arithmetic mod n (n odd, u64) - type I = Wrapping; + type I = u64; fn new(n: u64) -> Self { let a = inv_mod_u64(n).wrapping_neg(); @@ -104,13 +103,13 @@ impl Arithmetic for Montgomery { fn from_u64(&self, x: u64) -> Self::I { // TODO: optimise! assert!(x < self.n); - let r = Wrapping((((x as u128) << 64) % self.n as u128) as u64); + let r = (((x as u128) << 64) % self.n as u128) as u64; debug_assert_eq!(x, self.to_u64(r)); r } fn to_u64(&self, n: Self::I) -> u64 { - self.reduce(n.0) + self.reduce(n) } fn add(&self, a: Self::I, b: Self::I) -> Self::I { @@ -134,7 +133,7 @@ impl Arithmetic for Montgomery { } fn mul(&self, a: Self::I, b: Self::I) -> Self::I { - let r = Wrapping(self.reduce((a * b).0)); + let r = self.reduce(a.wrapping_mul(b)); // Check that r (reduced back to the usual representation) equals // a*b % n From cb6051c5804cd377212fb7652f985360b2ba80d8 Mon Sep 17 00:00:00 2001 From: nicoo Date: Mon, 15 Jun 2020 19:00:31 +0200 Subject: [PATCH 09/14] factor::numeric::Montgomery: Fix overflow bug --- src/uu/factor/src/numeric.rs | 27 ++++++++++++++++++--------- 1 file changed, 18 insertions(+), 9 deletions(-) diff --git a/src/uu/factor/src/numeric.rs b/src/uu/factor/src/numeric.rs index 7c6e7ac3e..0525e70b6 100644 --- a/src/uu/factor/src/numeric.rs +++ b/src/uu/factor/src/numeric.rs @@ -71,16 +71,25 @@ pub(crate) struct Montgomery { impl Montgomery { /// computes x/R mod n efficiently - fn reduce(&self, x: u64) -> u64 { + fn reduce(&self, x: u128) -> u64 { + debug_assert!(x < (self.n as u128) << 64); // TODO: optimiiiiiiise let Montgomery { a, n } = self; - let t = x.wrapping_mul(*a); - let nt = (*n as u128) * (t as u128); - let y = ((x as u128 + nt) >> 64) as u64; - if y >= *n { - y - n + let m = (x as u64).wrapping_mul(*a); + let nm = (*n as u128) * (m as u128); + let (xnm, overflow) = (x as u128).overflowing_add(nm); // x + n*m + debug_assert_eq!(xnm % (1 << 64), 0); + + if !overflow { + let y = (xnm >> 64) as u64 + overflow as u64; // (x + n*m) / R + if y >= *n { + y - n + } else { + y + } } else { - y + // y = (2¹²⁸ + xnm)/2⁶⁴ - n = xnm/2⁶⁴ + (2⁶⁴ - n) + (xnm >> 64) as u64 + n.wrapping_neg() } } } @@ -109,7 +118,7 @@ impl Arithmetic for Montgomery { } fn to_u64(&self, n: Self::I) -> u64 { - self.reduce(n) + self.reduce(n as u128) } fn add(&self, a: Self::I, b: Self::I) -> Self::I { @@ -133,7 +142,7 @@ impl Arithmetic for Montgomery { } fn mul(&self, a: Self::I, b: Self::I) -> Self::I { - let r = self.reduce(a.wrapping_mul(b)); + let r = self.reduce((a as u128) * (b as u128)); // Check that r (reduced back to the usual representation) equals // a*b % n From 4851619d622f6dd271b91f19d29f6435db8a8cf2 Mon Sep 17 00:00:00 2001 From: nicoo Date: Mon, 15 Jun 2020 23:04:46 +0200 Subject: [PATCH 10/14] factor::miller_rabin: Avoid repeatedly transforming 1 and -1 Approx. 25% speedup --- src/uu/factor/src/miller_rabin.rs | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/src/uu/factor/src/miller_rabin.rs b/src/uu/factor/src/miller_rabin.rs index 978cd6e47..f000485c7 100644 --- a/src/uu/factor/src/miller_rabin.rs +++ b/src/uu/factor/src/miller_rabin.rs @@ -38,6 +38,9 @@ pub(crate) fn test(m: A) -> Result { let i = (n - 1).trailing_zeros(); let r = (n - 1) >> i; + let one = m.one(); + let minus_one = m.minus_one(); + for _a in BASIS.iter() { let _a = _a % n; if _a == 0 { @@ -55,21 +58,21 @@ pub(crate) fn test(m: A) -> Result { for _ in 0..i { y = m.mul(y, y) } - if y != m.one() { + if y != one { return Pseudoprime; }; } - if x == m.one() || x == m.minus_one() { + if x == one || x == minus_one { break; } loop { let y = m.mul(x, x); - if y == m.one() { + if y == one { return Composite(gcd(m.to_u64(x) - 1, m.modulus())); } - if y == m.minus_one() { + if y == minus_one { // This basis element is not a witness of `n` being composite. // Keep looking. break; From f1788d9e7058c3cb80d54277b3260755a6275965 Mon Sep 17 00:00:00 2001 From: nicoo Date: Tue, 16 Jun 2020 01:17:16 +0200 Subject: [PATCH 11/14] fixup! factor::numeric::Montgomery: Fix overflow bug --- src/uu/factor/src/numeric.rs | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/src/uu/factor/src/numeric.rs b/src/uu/factor/src/numeric.rs index 0525e70b6..2d5d1ec21 100644 --- a/src/uu/factor/src/numeric.rs +++ b/src/uu/factor/src/numeric.rs @@ -80,16 +80,16 @@ impl Montgomery { let (xnm, overflow) = (x as u128).overflowing_add(nm); // x + n*m debug_assert_eq!(xnm % (1 << 64), 0); - if !overflow { - let y = (xnm >> 64) as u64 + overflow as u64; // (x + n*m) / R - if y >= *n { - y - n - } else { - y - } + // (x + n*m) / R + // in case of overflow, this is (2¹²⁸ + xnm)/2⁶⁴ - n = xnm/2⁶⁴ + (2⁶⁴ - n) + let y = + (xnm >> 64) as u64 + + if !overflow { 0 } else { n.wrapping_neg() }; + + if y >= *n { + y - n } else { - // y = (2¹²⁸ + xnm)/2⁶⁴ - n = xnm/2⁶⁴ + (2⁶⁴ - n) - (xnm >> 64) as u64 + n.wrapping_neg() + y } } } From 334e02786dd92e8e4cb0cc5bd2e3985aa829b822 Mon Sep 17 00:00:00 2001 From: nicoo Date: Tue, 16 Jun 2020 15:43:25 +0200 Subject: [PATCH 12/14] factor: Run `cargo fmt` --- src/uu/factor/src/numeric.rs | 4 +--- src/uu/factor/src/rho.rs | 2 +- 2 files changed, 2 insertions(+), 4 deletions(-) diff --git a/src/uu/factor/src/numeric.rs b/src/uu/factor/src/numeric.rs index 2d5d1ec21..f1dc94eed 100644 --- a/src/uu/factor/src/numeric.rs +++ b/src/uu/factor/src/numeric.rs @@ -82,9 +82,7 @@ impl Montgomery { // (x + n*m) / R // in case of overflow, this is (2¹²⁸ + xnm)/2⁶⁴ - n = xnm/2⁶⁴ + (2⁶⁴ - n) - let y = - (xnm >> 64) as u64 + - if !overflow { 0 } else { n.wrapping_neg() }; + let y = (xnm >> 64) as u64 + if !overflow { 0 } else { n.wrapping_neg() }; if y >= *n { y - n diff --git a/src/uu/factor/src/rho.rs b/src/uu/factor/src/rho.rs index a6608747a..5416218e1 100644 --- a/src/uu/factor/src/rho.rs +++ b/src/uu/factor/src/rho.rs @@ -3,9 +3,9 @@ use rand::rngs::SmallRng; use rand::{thread_rng, SeedableRng}; use std::cmp::{max, min}; -use crate::{Factors,miller_rabin}; use crate::miller_rabin::Result::*; use crate::numeric::*; +use crate::{miller_rabin, Factors}; fn find_divisor(n: A) -> u64 { #![allow(clippy::many_single_char_names)] From d1470dadf8d124dabd0fc6a55c402030cb94e7de Mon Sep 17 00:00:00 2001 From: nicoo Date: Tue, 16 Jun 2020 15:45:10 +0200 Subject: [PATCH 13/14] factor::numeric::gcd: Silence the (erroneous) dead code lint --- src/uu/factor/src/numeric.rs | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/src/uu/factor/src/numeric.rs b/src/uu/factor/src/numeric.rs index f1dc94eed..c0f9d9f47 100644 --- a/src/uu/factor/src/numeric.rs +++ b/src/uu/factor/src/numeric.rs @@ -9,7 +9,10 @@ use std::mem::swap; -pub fn gcd(mut a: u64, mut b: u64) -> u64 { +// This is incorrectly reported as dead code, +// presumably when included in build.rs. +#[allow(dead_code)] +pub(crate) fn gcd(mut a: u64, mut b: u64) -> u64 { while b > 0 { a %= b; swap(&mut a, &mut b); From fb08d9ff9e7264d91850a1042ab3d6db39044310 Mon Sep 17 00:00:00 2001 From: nicoo Date: Thu, 18 Jun 2020 13:01:55 +0200 Subject: [PATCH 14/14] factor::numeric::Montgomery::add: Deal with rare overflow case --- src/uu/factor/src/numeric.rs | 12 +++++++++++- 1 file changed, 11 insertions(+), 1 deletion(-) diff --git a/src/uu/factor/src/numeric.rs b/src/uu/factor/src/numeric.rs index c0f9d9f47..4e0cee072 100644 --- a/src/uu/factor/src/numeric.rs +++ b/src/uu/factor/src/numeric.rs @@ -123,7 +123,17 @@ impl Arithmetic for Montgomery { } fn add(&self, a: Self::I, b: Self::I) -> Self::I { - let r = a + b; + let (r, overflow) = a.overflowing_add(b); + + // In case of overflow, a+b = 2⁶⁴ + r = (2⁶⁴ - n) + r (working mod n) + let r = if !overflow { + r + } else { + r + self.n.wrapping_neg() + }; + + // Normalise to [0; n[ + let r = if r < self.n { r } else { r - self.n }; // Check that r (reduced back to the usual representation) equals // a+b % n