[WIP] Add an example of generating large primes

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
+}

src/biguint.rs

@@ -1653,26 +1653,26 @@ 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() {
             return monty_modpow(self, exponent, modulus);
         }
 
         // 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() {

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)
 }