diff --git a/examples/primes.rs b/examples/primes.rs new file mode 100644 index 00000000..717dac05 --- /dev/null +++ b/examples/primes.rs @@ -0,0 +1,152 @@ +//! Example generating prime numbers in a few large bit sizes. + +#![cfg_attr(not(feature = "rand"), allow(unused))] + +extern crate num_bigint; +extern crate num_integer; +extern crate num_traits; +extern crate rand; + +use num_bigint::{BigUint, RandBigInt}; +use num_integer::Integer; +use num_traits::{One, ToPrimitive, Zero}; +use rand::{Rng, SeedableRng, XorShiftRng}; +use std::env; +use std::u64; + +#[cfg(not(feature = "rand"))] +fn main() { + println!("This example requires `--features rand`"); +} + +#[cfg(feature = "rand")] +fn main() { + // Use a seeded Rng for benchmarking purposes. + let seeded = env::var_os("SEEDED_PRIMES").is_some(); + let (mut thread_rng, mut seeded_rng); + let mut rng: &mut Rng = if seeded { + seeded_rng = XorShiftRng::from_seed([4, 3, 2, 1]); + &mut seeded_rng + } else { + thread_rng = rand::thread_rng(); + &mut thread_rng + }; + + for &bits in &[256, 512, 1024, 2048] { + let min = BigUint::one() << (bits - 1); + for i in 1usize.. { + let mut n = &min + rng.gen_biguint_below(&min); + if n.is_even() { + n += 1u32; + } + if is_prime(&n) { + println!("Found a {}-bit prime in {} attempts:", bits, i); + println!("{}", n); + println!(""); + break; + } + } + } +} + +/// Attempt to determine whether n is prime. This algorithm is exact for +/// numbers less than 2^64, and highly probabilistic for other numbers. +#[cfg(feature = "rand")] +fn is_prime(n: &BigUint) -> bool { + // guarantees of primality given by Sloane's [A014233](https://oeis.org/A014233) + const A014233: [(u8, u64); 12] = [ + (2, 2_047), + (3, 1_373_653), + (5, 25_326_001), + (7, 3_215_031_751), + (11, 2_152_302_898_747), + (13, 3_474_749_660_383), + (17, 341_550_071_728_321), + (19, 341_550_071_728_321), + (23, 3_825_123_056_546_413_051), + (29, 3_825_123_056_546_413_051), + (31, 3_825_123_056_546_413_051), + (37, u64::MAX), + ]; + + let n64 = n.to_u64(); + if n.is_even() { + return n64 == Some(2); + } else if n64 == Some(1) { + return false; + } + + // try some simple divisors + for &(p, _) in &A014233[1..] { + if n64 == Some(u64::from(p)) { + return true; + } + if (n % p).is_zero() { + return false; + } + } + + // try series of primes as witnesses + for &(p, u) in &A014233 { + if witness(&BigUint::from(p), n) { + return false; + } + if let Some(n) = n64 { + if n < u || u == u64::MAX { + return true; + } + } + } + + // Now we're into the "probably prime" arena... + // Generate a few witnesses in the range `[2, n - 2]` + let mut rng = rand::thread_rng(); + let max = n - 3u32; + for _ in 0..10 { + let w = rng.gen_biguint_below(&max) + 2u32; + if witness(&w, n) { + return false; + } + } + true +} + +/// Test a possible witness to the compositeness of n. +/// +/// Computes a**(n-1) (mod n), and checks if the result gives evidence that +/// n is composite. Returning false means that n is likely prime, but not +/// necessarily. +#[cfg(feature = "rand")] +fn witness(a: &BigUint, n: &BigUint) -> bool { + // let n-1 = (2**t)*u, where t ≥ 1 and u is odd + // TODO `trailing_zeros` would make this trivial + let n_m1 = n - 1u32; + let (t, u) = if n_m1.is_even() { + let mut t = 1usize; + let mut u = n_m1.clone() >> 1; + while u.is_even() { + t += 1; + u >>= 1; + } + (t, u) + } else { + (0, n_m1.clone()) + }; + + let mut x = a.modpow(&u, n); + if x.is_one() || x == n_m1 { + return false; + } + + for _ in 1..t { + x = &x * &x % n; + if x.is_one() { + return true; + } + if x == n_m1 { + return false; + } + } + + true +} diff --git a/src/biguint.rs b/src/biguint.rs index 83102476..bbeeb1a4 100644 --- a/src/biguint.rs +++ b/src/biguint.rs @@ -1653,6 +1653,9 @@ impl BigUint { /// Panics if the modulus is zero. pub fn modpow(&self, exponent: &Self, modulus: &Self) -> Self { assert!(!modulus.is_zero(), "divide by zero!"); + if modulus.is_one() { return BigUint::zero(); } + if exponent.is_zero() { return BigUint::one(); } + if self.is_zero() || self.is_one() { return self.clone(); } // For an odd modulus, we can use Montgomery multiplication in base 2^32. if modulus.is_odd() { @@ -1660,19 +1663,16 @@ impl BigUint { } // Otherwise do basically the same as `num::pow`, but with a modulus. - let one = BigUint::one(); - if exponent.is_zero() { return one; } - let mut base = self % modulus; let mut exp = exponent.clone(); while exp.is_even() { base = &base * &base % modulus; exp >>= 1; } - if exp == one { return base } + if exp.is_one() { return base } let mut acc = base.clone(); - while exp > one { + while !exp.is_one() { exp >>= 1; base = &base * &base % modulus; if exp.is_odd() { diff --git a/src/monty.rs b/src/monty.rs index 3cf526b8..29204b3f 100644 --- a/src/monty.rs +++ b/src/monty.rs @@ -1,7 +1,8 @@ use integer::Integer; -use traits::Zero; +use traits::{One, Zero}; use biguint::BigUint; +use big_digit::BITS; struct MontyReducer<'a> { n: &'a BigUint, @@ -49,79 +50,96 @@ impl<'a> MontyReducer<'a> { let n0inv = inv_mod_u32(n.data[0]); MontyReducer { n: n, n0inv: n0inv } } -} -// Montgomery Reduction -// -// Reference: -// Brent & Zimmermann, Modern Computer Arithmetic, v0.5.9, Algorithm 2.6 -fn monty_redc(a: BigUint, mr: &MontyReducer) -> BigUint { - let mut c = a.data; - let n = &mr.n.data; - let n_size = n.len(); - - // Allocate sufficient work space - c.resize(2 * n_size + 2, 0); - - // β is the size of a word, in this case 32 bits. So "a mod β" is - // equivalent to masking a to 32 bits. - // mu <- -N^(-1) mod β - let mu = 0u32.wrapping_sub(mr.n0inv); - - // 1: for i = 0 to (n-1) - for i in 0..n_size { - // 2: q_i <- mu*c_i mod β - let q_i = c[i].wrapping_mul(mu); - - // 3: C <- C + q_i * N * β^i - super::algorithms::mac_digit(&mut c[i..], n, q_i); + /// Map a number to the Montgomery domain + fn map(&self, x: &BigUint) -> BigUint { + let shift = self.n.data.len() * BITS; + (x << shift) % self.n } - // 4: R <- C * β^(-n) - // This is an n-word bitshift, equivalent to skipping n words. - let ret = BigUint::new(c[n_size..].to_vec()); + /// Montgomery Reduction + /// + /// Reference: + /// Brent & Zimmermann, Modern Computer Arithmetic, v0.5.9, Algorithm 2.6 + fn redc(&self, a: BigUint) -> BigUint { + let mut c = a.data; + let n = &self.n.data; + let n_size = n.len(); + + // Allocate sufficient work space + c.resize(2 * n_size + 2, 0); + + // β is the size of a word, in this case 32 bits. So "a mod β" is + // equivalent to masking a to 32 bits. + // mu <- -N^(-1) mod β + let mu = 0u32.wrapping_sub(self.n0inv); + + // 1: for i = 0 to (n-1) + for i in 0..n_size { + // 2: q_i <- mu*c_i mod β + let q_i = c[i].wrapping_mul(mu); + + // 3: C <- C + q_i * N * β^i + super::algorithms::mac_digit(&mut c[i..], n, q_i); + } + + // 4: R <- C * β^(-n) + // This is an n-word bitshift, equivalent to skipping n words. + let ret = BigUint::new(c[n_size..].to_vec()); - // 5: if R >= β^n then return R-N else return R. - if &ret < mr.n { - ret - } else { - ret - mr.n + // 5: if R >= β^n then return R-N else return R. + if &ret < self.n { + ret + } else { + ret - self.n + } } -} -// Montgomery Multiplication -fn monty_mult(a: BigUint, b: &BigUint, mr: &MontyReducer) -> BigUint { - monty_redc(a * b, mr) -} + /// Montgomery Multiplication + fn mul(&self, a: BigUint, b: &BigUint) -> BigUint { + self.redc(a * b) + } + + /// Montgomery Squaring + fn square(&self, a: BigUint) -> BigUint { + // TODO: Replace with an optimised squaring function + self.redc(&a * &a) + } + + /// Montgomery Exponentiation + fn pow(&self, mut base: BigUint, exp: &BigUint) -> BigUint { + debug_assert!(!exp.is_zero()); + + // Binary exponentiation + let mut exp = exp.clone(); + while exp.is_even() { + base = self.square(base); + exp >>= 1; + } + if exp.is_one() { return base; } + + let mut acc = base.clone(); + while !exp.is_one() { + exp >>= 1; + base = self.square(base); + if exp.is_odd() { + acc = self.mul(acc, &base); + } + } -// Montgomery Squaring -fn monty_sqr(a: BigUint, mr: &MontyReducer) -> BigUint { - // TODO: Replace with an optimised squaring function - monty_redc(&a * &a, mr) + acc + } } pub fn monty_modpow(a: &BigUint, exp: &BigUint, modulus: &BigUint) -> BigUint{ let mr = MontyReducer::new(modulus); - // Calculate the Montgomery parameter - let mut v = vec![0; modulus.data.len()]; - v.push(1); - let r = BigUint::new(v); - // Map the base to the Montgomery domain - let mut apri = a * &r % modulus; - - // Binary exponentiation - let mut ans = &r % modulus; - let mut e = exp.clone(); - while !e.is_zero() { - if e.is_odd() { - ans = monty_mult(ans, &apri, &mr); - } - apri = monty_sqr(apri, &mr); - e = e >> 1; - } + let base = mr.map(a); + + // Do the computation + let answer = mr.pow(base, exp); // Map the result back to the residues domain - monty_redc(ans, &mr) + mr.redc(answer) }