feat: implement prime generation related features

also switches the benchmarks to criterion
2019-01-15 17:02:12 +03:00 · 2019-01-15 17:02:12 +03:00 · d4c0b721c9
commit d4c0b721c9
parent b3f10f4788
22 changed files with 1705 additions and 549 deletions
--- a/Cargo.toml
+++ b/Cargo.toml
@ -1,8 +1,8 @@
 [package]
 authors = ["The Rust Project Developers"]
 description = "Big integer implementation for Rust"
-documentation = "https://docs.rs/num-bigint"
-homepage = "https://github.com/rust-num/num-bigint"
+documentation = "https://docs.rs/num-bigint-dig"
+homepage = "https://github.com/dignifiedquire/num-bigint"
 keywords = ["mathematics", "numerics", "bignum"]
 categories = [ "algorithms", "data-structures", "science" ]
 license = "MIT/Apache-2.0"
@ -11,22 +11,15 @@ repository = "https://github.com/dignifiedquier/num-bigint"
 version = "0.2.1"
 readme = "README.md"
 build = "build.rs"
+autobenches = false

 [package.metadata.docs.rs]
-features = ["std", "serde", "rand"]
-
-[[bench]]
-name = "bigint"
-
-[[bench]]
-name = "factorial"
-
-[[bench]]
-name = "gcd"
+features = ["std", "serde", "rand", "prime"]

 [[bench]]
 harness = false
-name = "shootout-pidigits"
+name = "bench_main"
+required-features = ["prime"]

 [dependencies.smallvec]
 version = "0.6.7"
@ -40,6 +33,11 @@ default-features = false
 version = "0.2.4"
 default-features = false

+[dependencies.num-iter]
+version = "0.1.37"
+default-features = false
+
+
 [dependencies.rand]
 optional = true
 version = "0.5"
@ -55,8 +53,19 @@ features = ["std"]
 [dev-dependencies.serde_test]
 version = "1.0"

+[dependencies.lazy_static]
+version = "1.2.0"
+
+[dependencies.byteorder]
+version = "1.2.7"
+
+[dev-dependencies]
+criterion = "0.2"
+rand_chacha = "0.1"
+
 [features]
 default = ["std", "i128", "u64_digit"]
 i128 = ["num-integer/i128", "num-traits/i128"]
 std = ["num-integer/std", "num-traits/std", "smallvec/std"]
 u64_digit = []
+prime = ["rand"]
--- a/README.md
+++ b/README.md
@ -35,6 +35,8 @@ Implementations for `i128` and `u128` are only available with Rust 1.26 and
 later.  The build script automatically detects this, but you can make it
 mandatory by enabling the `i128` crate feature.

+The `prime` feature gate enables algorithms and support for dealing with large primes.
+
 ## Releases

 Release notes are available in [RELEASES.md](RELEASES.md).
@ -63,3 +65,9 @@ table offers a brief comparison to a few alternatives.
 [`rust-gmp`]: https://crates.io/crates/rust-gmp
 [`ramp`]: https://crates.io/crates/ramp
 [`apint`]: https://crates.io/crates/apint
+
+## Benchmarks
+
+```
+cargo bench --features prime
+```
--- a/benches/bench_main.rs
+++ b/benches/bench_main.rs
@ -0,0 +1,16 @@
+#[macro_use]
+extern crate criterion;
+extern crate num_bigint_dig as num_bigint;
+extern crate num_integer;
+extern crate num_traits;
+extern crate rand;
+extern crate rand_chacha;
+
+mod benchmarks;
+
+criterion_main! {
+    benchmarks::prime_benches::benches,
+    benchmarks::gcd::benches,
+    benchmarks::factorial::benches,
+    benchmarks::bigint::benches,
+}
--- a/benches/benchmarks/bigint.rs
+++ b/benches/benchmarks/bigint.rs
@ -0,0 +1,381 @@
+use criterion::Criterion;
+use num_bigint::{BigInt, BigUint, RandBigInt};
+use num_traits::{FromPrimitive, Num, One, Pow, Zero};
+use rand::{SeedableRng, StdRng};
+use std::mem::replace;
+
+fn get_rng() -> StdRng {
+    let mut seed = [0; 32];
+    for i in 1..32 {
+        seed[usize::from(i)] = i;
+    }
+    SeedableRng::from_seed(seed)
+}
+
+fn multiply_bench(c: &mut Criterion, name: String, xbits: usize, ybits: usize) {
+    let mut rng = get_rng();
+    let x = rng.gen_bigint(xbits);
+    let y = rng.gen_bigint(ybits);
+
+    c.bench_function(&name, move |b| b.iter(|| &x * &y));
+}
+
+fn divide_bench(c: &mut Criterion, name: String, xbits: usize, ybits: usize) {
+    let mut rng = get_rng();
+    let x = rng.gen_bigint(xbits);
+    let y = rng.gen_bigint(ybits);
+
+    c.bench_function(&name, move |b| b.iter(|| &x / &y));
+}
+
+fn factorial(n: usize) -> BigUint {
+    let mut f: BigUint = One::one();
+    for i in 1..(n + 1) {
+        let bu: BigUint = FromPrimitive::from_usize(i).unwrap();
+        f = f * bu;
+    }
+    f
+}
+
+/// Compute Fibonacci numbers
+fn fib(n: usize) -> BigUint {
+    let mut f0: BigUint = Zero::zero();
+    let mut f1: BigUint = One::one();
+    for _ in 0..n {
+        let f2 = f0 + &f1;
+        f0 = replace(&mut f1, f2);
+    }
+    f0
+}
+
+/// Compute Fibonacci numbers with two ops per iteration
+/// (add and subtract, like issue #200)
+fn fib2(n: usize) -> BigUint {
+    let mut f0: BigUint = Zero::zero();
+    let mut f1: BigUint = One::one();
+    for _ in 0..n {
+        f1 = f1 + &f0;
+        f0 = &f1 - f0;
+    }
+    f0
+}
+
+fn multiply_0(c: &mut Criterion) {
+    multiply_bench(c, "multiply_0".to_string(), 1 << 8, 1 << 8);
+}
+
+fn multiply_1(c: &mut Criterion) {
+    multiply_bench(c, "multiply_1".to_string(), 1 << 8, 1 << 16);
+}
+
+fn multiply_2(c: &mut Criterion) {
+    multiply_bench(c, "multiply_2".to_string(), 1 << 16, 1 << 16);
+}
+
+fn multiply_3(c: &mut Criterion) {
+    multiply_bench(c, "multiply_3".to_string(), 1 << 16, 1 << 17);
+}
+
+fn divide_0(c: &mut Criterion) {
+    divide_bench(c, "divide_0".to_string(), 1 << 8, 1 << 6);
+}
+
+fn divide_1(c: &mut Criterion) {
+    divide_bench(c, "divide_1".to_string(), 1 << 12, 1 << 8);
+}
+
+fn divide_2(c: &mut Criterion) {
+    divide_bench(c, "divide_2".to_string(), 1 << 16, 1 << 12);
+}
+
+fn factorial_100(c: &mut Criterion) {
+    c.bench_function("factorial_100", move |b| b.iter(|| factorial(100)));
+}
+
+fn fib_100(c: &mut Criterion) {
+    c.bench_function("fib_100", move |b| b.iter(|| fib(100)));
+}
+
+fn fib_1000(c: &mut Criterion) {
+    c.bench_function("fib_1000", move |b| b.iter(|| fib(1000)));
+}
+
+fn fib_10000(c: &mut Criterion) {
+    c.bench_function("fib_10000", move |b| b.iter(|| fib(10000)));
+}
+
+fn fib2_100(c: &mut Criterion) {
+    c.bench_function("fib2_100", move |b| b.iter(|| fib2(100)));
+}
+
+fn fib2_1000(c: &mut Criterion) {
+    c.bench_function("fib2_1000", move |b| b.iter(|| fib2(1000)));
+}
+
+fn fib2_10000(c: &mut Criterion) {
+    c.bench_function("fib2_10000", move |b| b.iter(|| fib2(10000)));
+}
+
+fn fac_to_string(c: &mut Criterion) {
+    let fac = factorial(100);
+    c.bench_function("fac_to_string", move |b| b.iter(|| fac.to_string()));
+}
+
+fn fib_to_string(c: &mut Criterion) {
+    let fib = fib(100);
+    c.bench_function("fib_to_string", move |b| b.iter(|| fib.to_string()));
+}
+
+fn to_str_radix_bench(c: &mut Criterion, radix: u32) {
+    let mut rng = get_rng();
+    let x = rng.gen_bigint(1009);
+    c.bench_function(&format!("to_str_radix_bench_{:?}", radix), move |b| {
+        b.iter(|| x.to_str_radix(radix))
+    });
+}
+
+fn to_str_radix_02(c: &mut Criterion) {
+    to_str_radix_bench(c, 2);
+}
+
+fn to_str_radix_08(c: &mut Criterion) {
+    to_str_radix_bench(c, 8);
+}
+
+fn to_str_radix_10(c: &mut Criterion) {
+    to_str_radix_bench(c, 10);
+}
+
+fn to_str_radix_16(c: &mut Criterion) {
+    to_str_radix_bench(c, 16);
+}
+
+fn to_str_radix_36(c: &mut Criterion) {
+    to_str_radix_bench(c, 36);
+}
+
+fn from_str_radix_bench(c: &mut Criterion, radix: u32) {
+    let mut rng = get_rng();
+    let x = rng.gen_bigint(1009);
+    let s = x.to_str_radix(radix);
+    assert_eq!(x, BigInt::from_str_radix(&s, radix).unwrap());
+
+    c.bench_function(&format!("from_str_radix_bench{:?}", radix), move |b| {
+        b.iter(|| BigInt::from_str_radix(&s, radix))
+    });
+}
+
+fn from_str_radix_02(c: &mut Criterion) {
+    from_str_radix_bench(c, 2);
+}
+
+fn from_str_radix_08(c: &mut Criterion) {
+    from_str_radix_bench(c, 8);
+}
+
+fn from_str_radix_10(c: &mut Criterion) {
+    from_str_radix_bench(c, 10);
+}
+
+fn from_str_radix_16(c: &mut Criterion) {
+    from_str_radix_bench(c, 16);
+}
+
+fn from_str_radix_36(c: &mut Criterion) {
+    from_str_radix_bench(c, 36);
+}
+
+fn rand_bench(c: &mut Criterion, bits: usize) {
+    let mut rng = get_rng();
+    c.bench_function(&format!("rand_bench_{:?}", bits), move |b| {
+        b.iter(|| rng.gen_bigint(bits))
+    });
+}
+
+fn rand_64(c: &mut Criterion) {
+    rand_bench(c, 1 << 6);
+}
+
+fn rand_256(c: &mut Criterion) {
+    rand_bench(c, 1 << 8);
+}
+
+fn rand_1009(c: &mut Criterion) {
+    rand_bench(c, 1009);
+}
+
+fn rand_2048(c: &mut Criterion) {
+    rand_bench(c, 1 << 11);
+}
+
+fn rand_4096(c: &mut Criterion) {
+    rand_bench(c, 1 << 12);
+}
+
+fn rand_8192(c: &mut Criterion) {
+    rand_bench(c, 1 << 13);
+}
+
+fn rand_65536(c: &mut Criterion) {
+    rand_bench(c, 1 << 16);
+}
+
+fn rand_131072(c: &mut Criterion) {
+    rand_bench(c, 1 << 17);
+}
+
+fn shl(c: &mut Criterion) {
+    let n = BigUint::one() << 1000;
+
+    c.bench_function("shl", move |b| {
+        b.iter(|| {
+            let mut m = n.clone();
+            for i in 0..50 {
+                m = m << i;
+            }
+        })
+    });
+}
+
+fn shr(c: &mut Criterion) {
+    let n = BigUint::one() << 2000;
+
+    c.bench_function("shr", move |b| {
+        b.iter(|| {
+            let mut m = n.clone();
+            for i in 0..50 {
+                m = m << i;
+            }
+        })
+    });
+}
+
+fn hash(c: &mut Criterion) {
+    use std::collections::HashSet;
+    let mut rng = get_rng();
+    let v: Vec<BigInt> = (1000..2000).map(|bits| rng.gen_bigint(bits)).collect();
+    c.bench_function("hash", move |b| {
+        b.iter(|| {
+            let h: HashSet<&BigInt> = v.iter().collect();
+            assert_eq!(h.len(), v.len());
+        })
+    });
+}
+
+fn pow_bench(c: &mut Criterion) {
+    c.bench_function("pow_bench", move |b| {
+        b.iter(|| {
+            let upper = 100_usize;
+            for i in 2..upper + 1 {
+                for j in 2..upper + 1 {
+                    let i_big = BigUint::from_usize(i).unwrap();
+                    i_big.pow(j);
+                }
+            }
+        })
+    });
+}
+
+/// This modulus is the prime from the 2048-bit MODP DH group:
+/// https://tools.ietf.org/html/rfc3526#section-3
+const RFC3526_2048BIT_MODP_GROUP: &'static str =
+    "\
+     FFFFFFFF_FFFFFFFF_C90FDAA2_2168C234_C4C6628B_80DC1CD1\
+     29024E08_8A67CC74_020BBEA6_3B139B22_514A0879_8E3404DD\
+     EF9519B3_CD3A431B_302B0A6D_F25F1437_4FE1356D_6D51C245\
+     E485B576_625E7EC6_F44C42E9_A637ED6B_0BFF5CB6_F406B7ED\
+     EE386BFB_5A899FA5_AE9F2411_7C4B1FE6_49286651_ECE45B3D\
+     C2007CB8_A163BF05_98DA4836_1C55D39A_69163FA8_FD24CF5F\
+     83655D23_DCA3AD96_1C62F356_208552BB_9ED52907_7096966D\
+     670C354E_4ABC9804_F1746C08_CA18217C_32905E46_2E36CE3B\
+     E39E772C_180E8603_9B2783A2_EC07A28F_B5C55DF0_6F4C52C9\
+     DE2BCBF6_95581718_3995497C_EA956AE5_15D22618_98FA0510\
+     15728E5A_8AACAA68_FFFFFFFF_FFFFFFFF";
+
+fn modpow(c: &mut Criterion) {
+    let mut rng = get_rng();
+    let base = rng.gen_biguint(2048);
+    let e = rng.gen_biguint(2048);
+    let m = BigUint::from_str_radix(RFC3526_2048BIT_MODP_GROUP, 16).unwrap();
+
+    c.bench_function("modpow", move |b| b.iter(|| base.modpow(&e, &m)));
+}
+
+fn modpow_even(c: &mut Criterion) {
+    let mut rng = get_rng();
+    let base = rng.gen_biguint(2048);
+    let e = rng.gen_biguint(2048);
+    // Make the modulus even, so monty (base-2^32) doesn't apply.
+    let m = BigUint::from_str_radix(RFC3526_2048BIT_MODP_GROUP, 16).unwrap() - 1u32;
+
+    c.bench_function("modpow_even", move |b| {
+        b.iter(|| base.modpow(&e, &m));
+    });
+}
+
+fn roots_sqrt(c: &mut Criterion) {
+    let mut rng = get_rng();
+    let x = rng.gen_biguint(2048);
+    c.bench_function("roots_sqrt", move |b| b.iter(|| x.sqrt()));
+}
+
+fn roots_cbrt(c: &mut Criterion) {
+    let mut rng = get_rng();
+    let x = rng.gen_biguint(2048);
+    c.bench_function("roots_cbrt", move |b| b.iter(|| x.cbrt()));
+}
+
+fn roots_nth_100(c: &mut Criterion) {
+    let mut rng = get_rng();
+    let x = rng.gen_biguint(2048);
+    c.bench_function("roots_nth_100", move |b| b.iter(|| x.nth_root(100)));
+}
+
+criterion_group! {
+    name = benches;
+    config = Criterion::default();
+    targets =
+        multiply_0,
+        multiply_1,
+        multiply_2,
+        multiply_3,
+        divide_0,
+        divide_1,
+        divide_2,
+        factorial_100,
+        fib_100,
+        fib_1000,
+        fib_10000,
+        fib2_100,
+        fib2_1000,
+        fib2_10000,
+        fac_to_string,
+        fib_to_string,
+        to_str_radix_02,
+        to_str_radix_08,
+        to_str_radix_10,
+        to_str_radix_16,
+        to_str_radix_36,
+        from_str_radix_02,
+        from_str_radix_08,
+        from_str_radix_10,
+        from_str_radix_16,
+        from_str_radix_36,
+        rand_64,
+        rand_256,
+        rand_1009,
+        rand_2048,
+        rand_4096,
+        rand_8192,
+        rand_65536,
+        rand_131072,
+        shl,
+        shr,
+        hash,
+        pow_bench,
+        modpow,
+        modpow_even,
+        roots_sqrt,
+        roots_cbrt,
+        roots_nth_100,
+}
--- a/benches/benchmarks/factorial.rs
+++ b/benches/benchmarks/factorial.rs
@ -0,0 +1,54 @@
+use criterion::Criterion;
+use num_bigint::BigUint;
+use num_traits::One;
+use std::ops::{Div, Mul};
+
+fn factorial_mul_biguint(c: &mut Criterion) {
+    c.bench_function("factorial_mul_biguint", move |b| {
+        b.iter(|| {
+            (1u32..1000)
+                .map(BigUint::from)
+                .fold(BigUint::one(), Mul::mul)
+        })
+    });
+}
+
+fn factorial_mul_u32(c: &mut Criterion) {
+    c.bench_function("factorial_mul_u32", move |b| {
+        b.iter(|| (1u32..1000).fold(BigUint::one(), Mul::mul))
+    });
+}
+
+// The division test is inspired by this blog comparison:
+// <https://tiehuis.github.io/big-integers-in-zig#division-test-single-limb>
+
+fn factorial_div_biguint(c: &mut Criterion) {
+    let n: BigUint = (1u32..1000).fold(BigUint::one(), Mul::mul);
+
+    c.bench_function("factorial_div_biguint", move |b| {
+        b.iter(|| {
+            (1u32..1000)
+                .rev()
+                .map(BigUint::from)
+                .fold(n.clone(), Div::div)
+        })
+    });
+}
+
+fn factorial_div_u32(c: &mut Criterion) {
+    let n: BigUint = (1u32..1000).fold(BigUint::one(), Mul::mul);
+
+    c.bench_function("factorial_div_u32", move |b| {
+        b.iter(|| (1u32..1000).rev().fold(n.clone(), Div::div))
+    });
+}
+
+criterion_group! {
+    name = benches;
+    config = Criterion::default();
+    targets =
+        factorial_mul_biguint,
+        factorial_mul_u32,
+        factorial_div_biguint,
+        factorial_div_u32,
+}
--- a/benches/benchmarks/gcd.rs
+++ b/benches/benchmarks/gcd.rs
@ -0,0 +1,83 @@
+use criterion::Criterion;
+use num_bigint::{BigUint, RandBigInt};
+use num_integer::Integer;
+use num_traits::Zero;
+use rand::{SeedableRng, StdRng};
+
+fn get_rng() -> StdRng {
+    let mut seed = [0; 32];
+    for i in 1..32 {
+        seed[usize::from(i)] = i;
+    }
+    SeedableRng::from_seed(seed)
+}
+
+fn bench(c: &mut Criterion, name: String, bits: usize, gcd: fn(&BigUint, &BigUint) -> BigUint) {
+    let mut rng = get_rng();
+    let x = rng.gen_biguint(bits);
+    let y = rng.gen_biguint(bits);
+
+    assert_eq!(euclid(&x, &y), x.gcd(&y));
+
+    c.bench_function(&name, move |b| b.iter(|| gcd(&x, &y)));
+}
+
+fn euclid(x: &BigUint, y: &BigUint) -> BigUint {
+    // Use Euclid's algorithm
+    let mut m = x.clone();
+    let mut n = y.clone();
+    while !m.is_zero() {
+        let temp = m;
+        m = n % &temp;
+        n = temp;
+    }
+    return n;
+}
+
+fn gcd_euclid_0064(c: &mut Criterion) {
+    bench(c, "gcd_euclid_0064".to_string(), 64, euclid);
+}
+
+fn gcd_euclid_0256(c: &mut Criterion) {
+    bench(c, "gcd_euclid_0256".to_string(), 256, euclid);
+}
+
+fn gcd_euclid_1024(c: &mut Criterion) {
+    bench(c, "gcd_euclid_1024".to_string(), 1024, euclid);
+}
+
+fn gcd_euclid_4096(c: &mut Criterion) {
+    bench(c, "gcd_euclid_4096".to_string(), 4096, euclid);
+}
+
+// Integer for BigUint now uses Stein for gcd
+
+fn gcd_stein_0064(c: &mut Criterion) {
+    bench(c, "gcd_stein_0064".to_string(), 64, BigUint::gcd);
+}
+
+fn gcd_stein_0256(c: &mut Criterion) {
+    bench(c, "gcd_stein_0256".to_string(), 256, BigUint::gcd);
+}
+
+fn gcd_stein_1024(c: &mut Criterion) {
+    bench(c, "gcd_stein_1024".to_string(), 1024, BigUint::gcd);
+}
+
+fn gcd_stein_4096(c: &mut Criterion) {
+    bench(c, "gcd_stein_4096".to_string(), 4096, BigUint::gcd);
+}
+
+criterion_group! {
+    name = benches;
+    config = Criterion::default();
+    targets =
+        gcd_euclid_0064,
+        gcd_euclid_0256,
+        gcd_euclid_1024,
+        gcd_euclid_4096,
+        gcd_stein_0064,
+        gcd_stein_0256,
+        gcd_stein_1024,
+        gcd_stein_4096,
+}
--- a/benches/benchmarks/mod.rs
+++ b/benches/benchmarks/mod.rs
@ -0,0 +1,4 @@
+pub mod bigint;
+pub mod factorial;
+pub mod gcd;
+pub mod prime_benches;
--- a/benches/benchmarks/prime_benches.rs
+++ b/benches/benchmarks/prime_benches.rs
@ -0,0 +1,85 @@
+use criterion::Criterion;
+use num_bigint::prime;
+use num_bigint::BigUint;
+use num_bigint::RandPrime;
+use rand::SeedableRng;
+use rand_chacha::ChaChaRng;
+
+const NUM: &'static str = "203956878356401977405765866929034577280193993314348263094772646453283062722701277632936616063144088173312372882677123879538709400158306567338328279154499698366071906766440037074217117805690872792848149112022286332144876183376326512083574821647933992961249917319836219304274280243803104015000563790123";
+
+fn probably_prime_0(c: &mut Criterion) {
+    let x = BigUint::parse_bytes(NUM.as_bytes(), 10).unwrap();
+
+    c.bench_function("probably_prime_0", move |b| {
+        b.iter(|| prime::probably_prime(&x, 0))
+    });
+}
+
+fn probably_prime_1(c: &mut Criterion) {
+    let x = BigUint::parse_bytes(NUM.as_bytes(), 10).unwrap();
+
+    c.bench_function("probably_prime_1", move |b| {
+        b.iter(|| prime::probably_prime(&x, 1))
+    });
+}
+
+fn probably_prime_5(c: &mut Criterion) {
+    let x = BigUint::parse_bytes(NUM.as_bytes(), 10).unwrap();
+
+    c.bench_function("probably_prime_5", move |b| {
+        b.iter(|| prime::probably_prime(&x, 5))
+    });
+}
+
+fn probably_prime_10(c: &mut Criterion) {
+    let x = BigUint::parse_bytes(NUM.as_bytes(), 10).unwrap();
+
+    c.bench_function("probably_prime_10", move |b| {
+        b.iter(|| prime::probably_prime(&x, 10))
+    });
+}
+
+fn probably_prime_20(c: &mut Criterion) {
+    let x = BigUint::parse_bytes(NUM.as_bytes(), 10).unwrap();
+
+    c.bench_function("probably_prime_20", move |b| {
+        b.iter(|| prime::probably_prime(&x, 20))
+    });
+}
+
+fn bench_prime_lucas(c: &mut Criterion) {
+    let x = BigUint::parse_bytes(NUM.as_bytes(), 10).unwrap();
+
+    c.bench_function("bench_prime_lucas", move |b| {
+        b.iter(|| prime::probably_prime_lucas(&x))
+    });
+}
+
+fn bench_prime_miller_rabin(c: &mut Criterion) {
+    let x = BigUint::parse_bytes(NUM.as_bytes(), 10).unwrap();
+
+    c.bench_function("bench_prime_miller_rabin", move |b| {
+        b.iter(|| prime::probably_prime_miller_rabin(&x, 1, true))
+    });
+}
+
+fn bench_gen_prime(c: &mut Criterion) {
+    c.bench_function("bench_gen_prime", move |b| {
+        let rng = &mut ChaChaRng::from_seed([0u8; 32]);
+        b.iter(|| rng.gen_prime(1024))
+    });
+}
+
+criterion_group! {
+    name = benches;
+    config = Criterion::default();
+    targets =
+        probably_prime_0,
+        probably_prime_1,
+        probably_prime_5,
+        probably_prime_10,
+        probably_prime_20,
+        bench_prime_lucas,
+        bench_prime_miller_rabin,
+        bench_gen_prime,
+}
--- a/benches/bigint.rs
+++ b/benches/bigint.rs
@ -1,369 +0,0 @@
-#![feature(test)]
-#![cfg(feature = "rand")]
-
-extern crate num_bigint_dig as num_bigint;
-extern crate num_integer;
-extern crate num_traits;
-extern crate rand;
-extern crate test;
-
-use num_bigint::{BigInt, BigUint, RandBigInt};
-use num_traits::{FromPrimitive, Num, One, Pow, Zero};
-use rand::{SeedableRng, StdRng};
-use std::mem::replace;
-use test::Bencher;
-
-fn get_rng() -> StdRng {
-    let mut seed = [0; 32];
-    for i in 1..32 {
-        seed[usize::from(i)] = i;
-    }
-    SeedableRng::from_seed(seed)
-}
-
-fn multiply_bench(b: &mut Bencher, xbits: usize, ybits: usize) {
-    let mut rng = get_rng();
-    let x = rng.gen_bigint(xbits);
-    let y = rng.gen_bigint(ybits);
-
-    b.iter(|| &x * &y);
-}
-
-fn divide_bench(b: &mut Bencher, xbits: usize, ybits: usize) {
-    let mut rng = get_rng();
-    let x = rng.gen_bigint(xbits);
-    let y = rng.gen_bigint(ybits);
-
-    b.iter(|| &x / &y);
-}
-
-fn factorial(n: usize) -> BigUint {
-    let mut f: BigUint = One::one();
-    for i in 1..(n + 1) {
-        let bu: BigUint = FromPrimitive::from_usize(i).unwrap();
-        f = f * bu;
-    }
-    f
-}
-
-/// Compute Fibonacci numbers
-fn fib(n: usize) -> BigUint {
-    let mut f0: BigUint = Zero::zero();
-    let mut f1: BigUint = One::one();
-    for _ in 0..n {
-        let f2 = f0 + &f1;
-        f0 = replace(&mut f1, f2);
-    }
-    f0
-}
-
-/// Compute Fibonacci numbers with two ops per iteration
-/// (add and subtract, like issue #200)
-fn fib2(n: usize) -> BigUint {
-    let mut f0: BigUint = Zero::zero();
-    let mut f1: BigUint = One::one();
-    for _ in 0..n {
-        f1 = f1 + &f0;
-        f0 = &f1 - f0;
-    }
-    f0
-}
-
-#[bench]
-fn multiply_0(b: &mut Bencher) {
-    multiply_bench(b, 1 << 8, 1 << 8);
-}
-
-#[bench]
-fn multiply_1(b: &mut Bencher) {
-    multiply_bench(b, 1 << 8, 1 << 16);
-}
-
-#[bench]
-fn multiply_2(b: &mut Bencher) {
-    multiply_bench(b, 1 << 16, 1 << 16);
-}
-
-#[bench]
-fn multiply_3(b: &mut Bencher) {
-    multiply_bench(b, 1 << 16, 1 << 17);
-}
-
-#[bench]
-fn divide_0(b: &mut Bencher) {
-    divide_bench(b, 1 << 8, 1 << 6);
-}
-
-#[bench]
-fn divide_1(b: &mut Bencher) {
-    divide_bench(b, 1 << 12, 1 << 8);
-}
-
-#[bench]
-fn divide_2(b: &mut Bencher) {
-    divide_bench(b, 1 << 16, 1 << 12);
-}
-
-#[bench]
-fn factorial_100(b: &mut Bencher) {
-    b.iter(|| factorial(100));
-}
-
-#[bench]
-fn fib_100(b: &mut Bencher) {
-    b.iter(|| fib(100));
-}
-
-#[bench]
-fn fib_1000(b: &mut Bencher) {
-    b.iter(|| fib(1000));
-}
-
-#[bench]
-fn fib_10000(b: &mut Bencher) {
-    b.iter(|| fib(10000));
-}
-
-#[bench]
-fn fib2_100(b: &mut Bencher) {
-    b.iter(|| fib2(100));
-}
-
-#[bench]
-fn fib2_1000(b: &mut Bencher) {
-    b.iter(|| fib2(1000));
-}
-
-#[bench]
-fn fib2_10000(b: &mut Bencher) {
-    b.iter(|| fib2(10000));
-}
-
-#[bench]
-fn fac_to_string(b: &mut Bencher) {
-    let fac = factorial(100);
-    b.iter(|| fac.to_string());
-}
-
-#[bench]
-fn fib_to_string(b: &mut Bencher) {
-    let fib = fib(100);
-    b.iter(|| fib.to_string());
-}
-
-fn to_str_radix_bench(b: &mut Bencher, radix: u32) {
-    let mut rng = get_rng();
-    let x = rng.gen_bigint(1009);
-    b.iter(|| x.to_str_radix(radix));
-}
-
-#[bench]
-fn to_str_radix_02(b: &mut Bencher) {
-    to_str_radix_bench(b, 2);
-}
-
-#[bench]
-fn to_str_radix_08(b: &mut Bencher) {
-    to_str_radix_bench(b, 8);
-}
-
-#[bench]
-fn to_str_radix_10(b: &mut Bencher) {
-    to_str_radix_bench(b, 10);
-}
-
-#[bench]
-fn to_str_radix_16(b: &mut Bencher) {
-    to_str_radix_bench(b, 16);
-}
-
-#[bench]
-fn to_str_radix_36(b: &mut Bencher) {
-    to_str_radix_bench(b, 36);
-}
-
-fn from_str_radix_bench(b: &mut Bencher, radix: u32) {
-    let mut rng = get_rng();
-    let x = rng.gen_bigint(1009);
-    let s = x.to_str_radix(radix);
-    assert_eq!(x, BigInt::from_str_radix(&s, radix).unwrap());
-    b.iter(|| BigInt::from_str_radix(&s, radix));
-}
-
-#[bench]
-fn from_str_radix_02(b: &mut Bencher) {
-    from_str_radix_bench(b, 2);
-}
-
-#[bench]
-fn from_str_radix_08(b: &mut Bencher) {
-    from_str_radix_bench(b, 8);
-}
-
-#[bench]
-fn from_str_radix_10(b: &mut Bencher) {
-    from_str_radix_bench(b, 10);
-}
-
-#[bench]
-fn from_str_radix_16(b: &mut Bencher) {
-    from_str_radix_bench(b, 16);
-}
-
-#[bench]
-fn from_str_radix_36(b: &mut Bencher) {
-    from_str_radix_bench(b, 36);
-}
-
-fn rand_bench(b: &mut Bencher, bits: usize) {
-    let mut rng = get_rng();
-
-    b.iter(|| rng.gen_bigint(bits));
-}
-
-#[bench]
-fn rand_64(b: &mut Bencher) {
-    rand_bench(b, 1 << 6);
-}
-
-#[bench]
-fn rand_256(b: &mut Bencher) {
-    rand_bench(b, 1 << 8);
-}
-
-#[bench]
-fn rand_1009(b: &mut Bencher) {
-    rand_bench(b, 1009);
-}
-
-#[bench]
-fn rand_2048(b: &mut Bencher) {
-    rand_bench(b, 1 << 11);
-}
-
-#[bench]
-fn rand_4096(b: &mut Bencher) {
-    rand_bench(b, 1 << 12);
-}
-
-#[bench]
-fn rand_8192(b: &mut Bencher) {
-    rand_bench(b, 1 << 13);
-}
-
-#[bench]
-fn rand_65536(b: &mut Bencher) {
-    rand_bench(b, 1 << 16);
-}
-
-#[bench]
-fn rand_131072(b: &mut Bencher) {
-    rand_bench(b, 1 << 17);
-}
-
-#[bench]
-fn shl(b: &mut Bencher) {
-    let n = BigUint::one() << 1000;
-    b.iter(|| {
-        let mut m = n.clone();
-        for i in 0..50 {
-            m = m << i;
-        }
-    })
-}
-
-#[bench]
-fn shr(b: &mut Bencher) {
-    let n = BigUint::one() << 2000;
-    b.iter(|| {
-        let mut m = n.clone();
-        for i in 0..50 {
-            m = m >> i;
-        }
-    })
-}
-
-#[bench]
-fn hash(b: &mut Bencher) {
-    use std::collections::HashSet;
-    let mut rng = get_rng();
-    let v: Vec<BigInt> = (1000..2000).map(|bits| rng.gen_bigint(bits)).collect();
-    b.iter(|| {
-        let h: HashSet<&BigInt> = v.iter().collect();
-        assert_eq!(h.len(), v.len());
-    });
-}
-
-#[bench]
-fn pow_bench(b: &mut Bencher) {
-    b.iter(|| {
-        let upper = 100_usize;
-        for i in 2..upper + 1 {
-            for j in 2..upper + 1 {
-                let i_big = BigUint::from_usize(i).unwrap();
-                i_big.pow(j);
-            }
-        }
-    });
-}
-
-/// This modulus is the prime from the 2048-bit MODP DH group:
-/// https://tools.ietf.org/html/rfc3526#section-3
-const RFC3526_2048BIT_MODP_GROUP: &'static str =
-    "\
-     FFFFFFFF_FFFFFFFF_C90FDAA2_2168C234_C4C6628B_80DC1CD1\
-     29024E08_8A67CC74_020BBEA6_3B139B22_514A0879_8E3404DD\
-     EF9519B3_CD3A431B_302B0A6D_F25F1437_4FE1356D_6D51C245\
-     E485B576_625E7EC6_F44C42E9_A637ED6B_0BFF5CB6_F406B7ED\
-     EE386BFB_5A899FA5_AE9F2411_7C4B1FE6_49286651_ECE45B3D\
-     C2007CB8_A163BF05_98DA4836_1C55D39A_69163FA8_FD24CF5F\
-     83655D23_DCA3AD96_1C62F356_208552BB_9ED52907_7096966D\
-     670C354E_4ABC9804_F1746C08_CA18217C_32905E46_2E36CE3B\
-     E39E772C_180E8603_9B2783A2_EC07A28F_B5C55DF0_6F4C52C9\
-     DE2BCBF6_95581718_3995497C_EA956AE5_15D22618_98FA0510\
-     15728E5A_8AACAA68_FFFFFFFF_FFFFFFFF";
-
-#[bench]
-fn modpow(b: &mut Bencher) {
-    let mut rng = get_rng();
-    let base = rng.gen_biguint(2048);
-    let e = rng.gen_biguint(2048);
-    let m = BigUint::from_str_radix(RFC3526_2048BIT_MODP_GROUP, 16).unwrap();
-
-    b.iter(|| base.modpow(&e, &m));
-}
-
-#[bench]
-fn modpow_even(b: &mut Bencher) {
-    let mut rng = get_rng();
-    let base = rng.gen_biguint(2048);
-    let e = rng.gen_biguint(2048);
-    // Make the modulus even, so monty (base-2^32) doesn't apply.
-    let m = BigUint::from_str_radix(RFC3526_2048BIT_MODP_GROUP, 16).unwrap() - 1u32;
-
-    b.iter(|| base.modpow(&e, &m));
-}
-
-#[bench]
-fn roots_sqrt(b: &mut Bencher) {
-    let mut rng = get_rng();
-    let x = rng.gen_biguint(2048);
-
-    b.iter(|| x.sqrt());
-}
-
-#[bench]
-fn roots_cbrt(b: &mut Bencher) {
-    let mut rng = get_rng();
-    let x = rng.gen_biguint(2048);
-
-    b.iter(|| x.cbrt());
-}
-
-#[bench]
-fn roots_nth_100(b: &mut Bencher) {
-    let mut rng = get_rng();
-    let x = rng.gen_biguint(2048);
-
-    b.iter(|| x.nth_root(100));
-}
--- a/benches/factorial.rs
+++ b/benches/factorial.rs
@ -1,44 +0,0 @@
-#![feature(test)]
-
-extern crate num_bigint_dig as num_bigint;
-extern crate num_traits;
-extern crate test;
-
-use num_bigint::BigUint;
-use num_traits::One;
-use std::ops::{Div, Mul};
-use test::Bencher;
-
-#[bench]
-fn factorial_mul_biguint(b: &mut Bencher) {
-    b.iter(|| {
-        (1u32..1000)
-            .map(BigUint::from)
-            .fold(BigUint::one(), Mul::mul)
-    });
-}
-
-#[bench]
-fn factorial_mul_u32(b: &mut Bencher) {
-    b.iter(|| (1u32..1000).fold(BigUint::one(), Mul::mul));
-}
-
-// The division test is inspired by this blog comparison:
-// <https://tiehuis.github.io/big-integers-in-zig#division-test-single-limb>
-
-#[bench]
-fn factorial_div_biguint(b: &mut Bencher) {
-    let n: BigUint = (1u32..1000).fold(BigUint::one(), Mul::mul);
-    b.iter(|| {
-        (1u32..1000)
-            .rev()
-            .map(BigUint::from)
-            .fold(n.clone(), Div::div)
-    });
-}
-
-#[bench]
-fn factorial_div_u32(b: &mut Bencher) {
-    let n: BigUint = (1u32..1000).fold(BigUint::one(), Mul::mul);
-    b.iter(|| (1u32..1000).rev().fold(n.clone(), Div::div));
-}
--- a/benches/gcd.rs
+++ b/benches/gcd.rs
@ -1,86 +0,0 @@
-#![feature(test)]
-#![cfg(feature = "rand")]
-
-extern crate num_bigint_dig as num_bigint;
-extern crate num_integer;
-extern crate num_traits;
-extern crate rand;
-extern crate test;
-
-use num_bigint::{BigUint, RandBigInt};
-use num_integer::Integer;
-use num_traits::Zero;
-use rand::{SeedableRng, StdRng};
-use test::Bencher;
-
-fn get_rng() -> StdRng {
-    let mut seed = [0; 32];
-    for i in 1..32 {
-        seed[usize::from(i)] = i;
-    }
-    SeedableRng::from_seed(seed)
-}
-
-fn bench(b: &mut Bencher, bits: usize, gcd: fn(&BigUint, &BigUint) -> BigUint) {
-    let mut rng = get_rng();
-    let x = rng.gen_biguint(bits);
-    let y = rng.gen_biguint(bits);
-
-    assert_eq!(euclid(&x, &y), x.gcd(&y));
-
-    b.iter(|| gcd(&x, &y));
-}
-
-fn euclid(x: &BigUint, y: &BigUint) -> BigUint {
-    // Use Euclid's algorithm
-    let mut m = x.clone();
-    let mut n = y.clone();
-    while !m.is_zero() {
-        let temp = m;
-        m = n % &temp;
-        n = temp;
-    }
-    return n;
-}
-
-#[bench]
-fn gcd_euclid_0064(b: &mut Bencher) {
-    bench(b, 64, euclid);
-}
-
-#[bench]
-fn gcd_euclid_0256(b: &mut Bencher) {
-    bench(b, 256, euclid);
-}
-
-#[bench]
-fn gcd_euclid_1024(b: &mut Bencher) {
-    bench(b, 1024, euclid);
-}
-
-#[bench]
-fn gcd_euclid_4096(b: &mut Bencher) {
-    bench(b, 4096, euclid);
-}
-
-// Integer for BigUint now uses Stein for gcd
-
-#[bench]
-fn gcd_stein_0064(b: &mut Bencher) {
-    bench(b, 64, BigUint::gcd);
-}
-
-#[bench]
-fn gcd_stein_0256(b: &mut Bencher) {
-    bench(b, 256, BigUint::gcd);
-}
-
-#[bench]
-fn gcd_stein_1024(b: &mut Bencher) {
-    bench(b, 1024, BigUint::gcd);
-}
-
-#[bench]
-fn gcd_stein_4096(b: &mut Bencher) {
-    bench(b, 4096, BigUint::gcd);
-}
--- a/src/algorithms.rs
+++ b/src/algorithms.rs
@ -1,21 +1,18 @@
+use crate::VEC_SIZE;
+use big_digit::{self, BigDigit, DoubleBigDigit, SignedDoubleBigDigit};
+use bigint::BigInt;
+use bigint::Sign;
+use bigint::Sign::{Minus, NoSign, Plus};
+use biguint::BigUint;
+use integer::Integer;
+use num_traits;
+use num_traits::{One, Signed, Zero};
+use smallvec::SmallVec;
 use std::borrow::Cow;
 use std::cmp;
 use std::cmp::Ordering::{self, Equal, Greater, Less};
 use std::iter::repeat;
 use std::mem;
-use traits;
-use traits::{One, Zero};
-
-use biguint::BigUint;
-
-use bigint::BigInt;
-use bigint::Sign;
-use bigint::Sign::{Minus, NoSign, Plus};
-
-use big_digit::{self, BigDigit, DoubleBigDigit, SignedDoubleBigDigit};
-
-use super::super::VEC_SIZE;
-use smallvec::SmallVec;

 // Generic functions for add/subtract/multiply with carry/borrow:

@ -606,11 +603,11 @@ pub fn div_rem(u: &BigUint, d: &BigUint) -> (BigUint, BigUint) {

 /// Find last set bit
 /// fls(0) == 0, fls(u32::MAX) == 32
-pub fn fls<T: traits::PrimInt>(v: T) -> usize {
+pub fn fls<T: num_traits::PrimInt>(v: T) -> usize {
    mem::size_of::<T>() * 8 - v.leading_zeros() as usize
 }

-pub fn ilog2<T: traits::PrimInt>(v: T) -> usize {
+pub fn ilog2<T: num_traits::PrimInt>(v: T) -> usize {
    fls(v) - 1
 }

@ -691,10 +688,142 @@ pub fn cmp_slice(a: &[BigDigit], b: &[BigDigit]) -> Ordering {
    return Equal;
 }

+// Few Functions taken from
+// https://github.com/RustCrypto/RSA/blob/master/src/math.rs
+
+/// Jacobi returns the Jacobi symbol (x/y), either +1, -1, or 0.
+/// The y argument must be an odd integer.
+pub fn jacobi(x: &BigInt, y: &BigInt) -> isize {
+    if !y.is_odd() {
+        panic!(
+            "invalid arguments, y must be an odd integer,but got {:?}",
+            y
+        );
+    }
+
+    let mut a = x.clone();
+    let mut b = y.clone();
+    let mut j = 1;
+
+    if b.is_negative() {
+        if a.is_negative() {
+            j = -1;
+        }
+        b = -b;
+    }
+
+    loop {
+        if b.is_one() {
+            return j;
+        }
+        if a.is_zero() {
+            return 0;
+        }
+
+        a = a.mod_floor(&b);
+        if a.is_zero() {
+            return 0;
+        }
+
+        // a > 0
+
+        // handle factors of 2 in a
+        let s = a.trailing_zeros().unwrap();
+        if s & 1 != 0 {
+            let bmod8 = b.get_limb(0) & 7;
+            if bmod8 == 3 || bmod8 == 5 {
+                j = -j;
+            }
+        }
+
+        let c = &a >> s; // a = 2^s*c
+
+        // swap numerator and denominator
+        if b.get_limb(0) & 3 == 3 && c.get_limb(0) & 3 == 3 {
+            j = -j
+        }
+
+        a = b;
+        b = c.clone();
+    }
+}
+
+/// Calculates the extended eucledian algorithm.
+/// See https://en.wikipedia.org/wiki/Extended_Euclidean_algorithm for details.
+/// The returned values are
+/// - greatest common divisor (1)
+/// - Bezout coefficients (2)
+// TODO: implement optimized variants
+pub fn extended_gcd(a: &BigUint, b: &BigUint) -> (BigInt, BigInt, BigInt) {
+    let mut a = BigInt::from_biguint(Plus, a.clone());
+    let mut b = BigInt::from_biguint(Plus, b.clone());
+
+    let mut ua = BigInt::one();
+    let mut va = BigInt::zero();
+
+    let mut ub = BigInt::zero();
+    let mut vb = BigInt::one();
+
+    let mut q;
+    let mut tmp;
+    let mut r;
+
+    while !b.is_zero() {
+        q = &a / &b;
+        r = &a % &b;
+
+        a = b;
+        b = r;
+
+        tmp = ua;
+        ua = ub.clone();
+        ub = tmp - &q * &ub;
+
+        tmp = va;
+        va = vb.clone();
+        vb = tmp - &q * &vb;
+    }
+
+    (a, ua, va)
+}
+
+/// Calculate the modular inverse of `a`.
+/// Implemenation is based on the naive version from wikipedia.
+#[inline]
+pub fn mod_inverse(g: Cow<BigInt>, n: &BigInt) -> Option<BigInt> {
+    assert!(g.as_ref() != n, "g must not be equal to n");
+    assert!(!n.is_negative(), "negative modulus not supported");
+
+    let n = n.abs();
+    let g = if g.is_negative() {
+        g.mod_floor(&n).to_biguint().unwrap()
+    } else {
+        g.to_biguint().unwrap()
+    };
+
+    let (d, x, _) = extended_gcd(&g, &n.to_biguint().unwrap());
+
+    if !d.is_one() {
+        return None;
+    }
+
+    if x.is_negative() {
+        Some(x + n)
+    } else {
+        Some(x)
+    }
+}
+
 #[cfg(test)]
 mod algorithm_tests {
+    // extern crate rand;
+
+    use algorithms::{extended_gcd, jacobi};
    use big_digit::BigDigit;
-    use traits::Num;
+    use integer::Integer;
+    use num_traits::Num;
+    use num_traits::{FromPrimitive, One};
+    use traits::ModInverse;
    use Sign::Plus;
    use {BigInt, BigUint};

@ -731,4 +860,132 @@ mod algorithm_tests {
        assert!(&a != &b);
        assert_ne!(&a, &b);
    }
+
+    #[test]
+    fn test_extended_gcd_example() {
+        // simple example for wikipedia
+        let a = BigUint::from_u32(240).unwrap();
+        let b = BigUint::from_u32(46).unwrap();
+        let (q, s_k, t_k) = extended_gcd(&a, &b);
+
+        assert_eq!(q, BigInt::from_i32(2).unwrap());
+        assert_eq!(s_k, BigInt::from_i32(-9).unwrap());
+        assert_eq!(t_k, BigInt::from_i32(47).unwrap());
+    }
+
+    #[test]
+    fn test_jacobi() {
+        let cases = [
+            [0, 1, 1],
+            [0, -1, 1],
+            [1, 1, 1],
+            [1, -1, 1],
+            [0, 5, 0],
+            [1, 5, 1],
+            [2, 5, -1],
+            [-2, 5, -1],
+            [2, -5, -1],
+            [-2, -5, 1],
+            [3, 5, -1],
+            [5, 5, 0],
+            [-5, 5, 0],
+            [6, 5, 1],
+            [6, -5, 1],
+            [-6, 5, 1],
+            [-6, -5, -1],
+        ];
+
+        for case in cases.iter() {
+            let x = BigInt::from_i64(case[0]).unwrap();
+            let y = BigInt::from_i64(case[1]).unwrap();
+
+            assert_eq!(case[2] as isize, jacobi(&x, &y), "jacobi({}, {})", x, y);
+        }
+    }
+
+    #[test]
+    fn test_mod_inverse() {
+        let tests = [
+            ["1234567", "458948883992"],
+	    ["239487239847", "2410312426921032588552076022197566074856950548502459942654116941958108831682612228890093858261341614673227141477904012196503648957050582631942730706805009223062734745341073406696246014589361659774041027169249453200378729434170325843778659198143763193776859869524088940195577346119843545301547043747207749969763750084308926339295559968882457872412993810129130294592999947926365264059284647209730384947211681434464714438488520940127459844288859336526896320919633919"],
+	    ["-10", "13"],
+            ["-6193420858199668535", "2881"],
+        ];
+
+        for test in &tests {
+            let element = BigInt::parse_bytes(test[0].as_bytes(), 10).unwrap();
+            let modulus = BigInt::parse_bytes(test[1].as_bytes(), 10).unwrap();
+
+            println!("{} modinv {}", element, modulus);
+            let inverse = element.clone().mod_inverse(&modulus).unwrap();
+            println!("inverse: {}", &inverse);
+            let cmp = (inverse * &element).mod_floor(&modulus);
+
+            assert_eq!(
+                cmp,
+                BigInt::one(),
+                "mod_inverse({}, {}) * {} % {} = {}, not 1",
+                &element,
+                &modulus,
+                &element,
+                &modulus,
+                &cmp
+            );
+        }
+
+        // exhaustive tests for small numbers
+        for n in 2..100 {
+            let modulus = BigInt::from_u64(n).unwrap();
+            for x in 1..n {
+                for sign in vec![1i64, -1i64] {
+                    let element = BigInt::from_i64(sign * x as i64).unwrap();
+                    let gcd = element.gcd(&modulus);
+
+                    if !gcd.is_one() {
+                        continue;
+                    }
+
+                    let inverse = element.clone().mod_inverse(&modulus).unwrap();
+                    let cmp = (&inverse * &element).mod_floor(&modulus);
+                    println!("inverse: {}", &inverse);
+                    assert_eq!(
+                        cmp,
+                        BigInt::one(),
+                        "mod_inverse({}, {}) * {} % {} = {}, not 1",
+                        &element,
+                        &modulus,
+                        &element,
+                        &modulus,
+                        &cmp
+                    );
+                }
+            }
+        }
+    }
+
+}
+
+#[cfg(feature = "prime")]
+mod random_prime_tests {
+    use algorithms::extended_gcd;
+    use bigrand::RandBigInt;
+    use rand::thread_rng;
+    use BigInt;
+    use Sign::Plus;
+
+    #[test]
+    fn test_extended_gcd_assumptions() {
+        let mut rng = thread_rng();
+
+        for i in 1..100 {
+            let a = rng.gen_biguint(i * 128);
+            let b = rng.gen_biguint(i * 128);
+            let (q, s_k, t_k) = extended_gcd(&a, &b);
+
+            let lhs = BigInt::from_biguint(Plus, a) * &s_k;
+            let rhs = BigInt::from_biguint(Plus, b) * &t_k;
+            assert_eq!(q, lhs + &rhs);
+        }
+    }
+
 }
--- a/src/bigint.rs
+++ b/src/bigint.rs
@ -18,22 +18,22 @@ use std::{i64, u64};
 use serde;

 use integer::{Integer, Roots};
-use traits::{
+use num_traits::{
    CheckedAdd, CheckedDiv, CheckedMul, CheckedSub, FromPrimitive, Num, One, Pow, Signed,
    ToPrimitive, Zero,
 };

-use smallvec::SmallVec;
-
 use self::Sign::{Minus, NoSign, Plus};
-use super::VEC_SIZE;
-
 use super::ParseBigIntError;
+use super::VEC_SIZE;
+use algorithms::mod_inverse;
 use big_digit::{self, BigDigit, DoubleBigDigit};
 use biguint;
 use biguint::to_str_radix_reversed;
 use biguint::{BigUint, IntDigits};
-
+use smallvec::SmallVec;
+use std::borrow::Cow;
+use traits::ModInverse;
 use IsizePromotion;
 use UsizePromotion;

@ -3027,6 +3027,20 @@ where
    }
 }

+// Mod Inverse
+
+impl<'a> ModInverse<&'a BigInt> for BigInt {
+    fn mod_inverse(self, m: &'a BigInt) -> Option<BigInt> {
+        mod_inverse(Cow::Owned(self), m)
+    }
+}
+
+impl ModInverse<BigInt> for BigInt {
+    fn mod_inverse(self, m: BigInt) -> Option<BigInt> {
+        mod_inverse(Cow::Owned(self), &m)
+    }
+}
+
 #[test]
 fn test_from_biguint() {
    fn check(inp_s: Sign, inp_n: usize, ans_s: Sign, ans_n: usize) {
--- a/src/bigrand.rs
+++ b/src/bigrand.rs
@ -3,6 +3,7 @@
 use rand::distributions::uniform::{SampleUniform, UniformSampler};
 use rand::prelude::*;
 use rand::AsByteSliceMut;
+use rand::Rng;

 use BigInt;
 use BigUint;
@ -10,9 +11,12 @@ use Sign::*;

 use big_digit::BigDigit;
 use bigint::{into_magnitude, magnitude};
-
 use integer::Integer;
-use traits::Zero;
+use num_iter::range_step;
+use num_traits::Zero;
+use num_traits::{FromPrimitive, ToPrimitive};
+
+use crate::prime::probably_prime;

 pub trait RandBigInt {
    /// Generate a random `BigUint` of the given bit size.
@ -216,3 +220,106 @@ impl Distribution<BigInt> for RandomBits {
        rng.gen_bigint(self.bits)
    }
 }
+
+/// A generic trait for generating random primes.
+///
+/// *Warning*: This is highly dependend on the provided random number generator,
+/// to provide actually random primes.
+///
+/// # Example
+/// ```
+/// extern crate rand;
+/// extern crate num_bigint_dig as num_bigint;
+///
+/// use rand::thread_rng;
+/// use num_bigint::RandPrime;
+///
+/// let mut rng = thread_rng();
+/// let p = rng.gen_prime(1024);
+/// assert_eq!(p.bits(), 1024);
+/// ```
+///
+pub trait RandPrime {
+    /// Generate a random prime number with as many bits as given.
+    fn gen_prime(&mut self, bits: usize) -> BigUint;
+}
+
+/// A list of small, prime numbers that allows us to rapidly
+/// exclude some fraction of composite candidates when searching for a random
+/// prime. This list is truncated at the point where smallPrimesProduct exceeds
+/// a u64. It does not include two because we ensure that the candidates are
+/// odd by construction.
+const SMALL_PRIMES: [u8; 15] = [3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53];
+
+lazy_static! {
+    /// The product of the values in SMALL_PRIMES and allows us
+    /// to reduce a candidate prime by this number and then determine whether it's
+    /// coprime to all the elements of SMALL_PRIMES without further BigUint
+    /// operations.
+    static ref SMALL_PRIMES_PRODUCT: BigUint = BigUint::from_u64(16_294_579_238_595_022_365).unwrap();
+}
+
+impl<R: Rng + ?Sized> RandPrime for R {
+    fn gen_prime(&mut self, bit_size: usize) -> BigUint {
+        if bit_size < 2 {
+            panic!("prime size must be at least 2-bit");
+        }
+
+        let mut b = bit_size % 8;
+        if b == 0 {
+            b = 8;
+        }
+
+        let bytes_len = (bit_size + 7) / 8;
+        let mut bytes = vec![0u8; bytes_len];
+
+        loop {
+            self.fill_bytes(&mut bytes);
+            // Clear bits in the first byte to make sure the candidate has a size <= bits.
+            bytes[0] &= ((1u32 << (b as u32)) - 1) as u8;
+
+            // Don't let the value be too small, i.e, set the most significant two bits.
+            // Setting the top two bits, rather than just the top bit,
+            // means that when two of these values are multiplied together,
+            // the result isn't ever one bit short.
+            if b >= 2 {
+                bytes[0] |= 3u8.wrapping_shl(b as u32 - 2);
+            } else {
+                // Here b==1, because b cannot be zero.
+                bytes[0] |= 1;
+                if bytes_len > 1 {
+                    bytes[1] |= 0x80;
+                }
+            }
+
+            // Make the value odd since an even number this large certainly isn't prime.
+            bytes[bytes_len - 1] |= 1u8;
+
+            let mut p = BigUint::from_bytes_be(&bytes);
+            // must always be a u64, as the SMALL_PRIMES_PRODUCT is a u64
+            let rem = (&p % &*SMALL_PRIMES_PRODUCT).to_u64().unwrap();
+
+            'next: for delta in range_step(0, 1 << 20, 2) {
+                let m = rem + delta;
+
+                for prime in &SMALL_PRIMES {
+                    if m % u64::from(*prime) == 0 && (bit_size > 6 || m != u64::from(*prime)) {
+                        continue 'next;
+                    }
+                }
+
+                if delta > 0 {
+                    p += BigUint::from_u64(delta).unwrap();
+                }
+
+                break;
+            }
+
+            // There is a tiny possibility that, by adding delta, we caused
+            // the number to be one bit too long. Thus we check bit length here.
+            if p.bits() == bit_size && probably_prime(&p, 20) {
+                return p;
+            }
+        }
+    }
+}
--- a/src/biguint.rs
+++ b/src/biguint.rs
@ -19,26 +19,27 @@ use std::{u64, u8};
 use serde;

 use integer::{Integer, Roots};
-use traits::{
+use num_traits::{
    CheckedAdd, CheckedDiv, CheckedMul, CheckedSub, Float, FromPrimitive, Num, One, Pow,
    ToPrimitive, Unsigned, Zero,
 };
+use BigInt;

 use big_digit::{self, BigDigit};

+use bigint::Sign::Plus;
 use smallvec::SmallVec;
+use traits::ModInverse;

-#[path = "algorithms.rs"]
-mod algorithms;
 #[path = "monty.rs"]
 mod monty;

-use self::algorithms::{__add2, __sub2rev, add2, sub2, sub2rev};
-use self::algorithms::{biguint_shl, biguint_shr};
-use self::algorithms::{cmp_slice, fls, ilog2};
-use self::algorithms::{div_rem, div_rem_digit, mac_with_carry, mul3, scalar_mul};
 use self::monty::monty_modpow;
 use super::VEC_SIZE;
+use crate::algorithms::{__add2, __sub2rev, add2, mod_inverse, sub2, sub2rev};
+use crate::algorithms::{biguint_shl, biguint_shr};
+use crate::algorithms::{cmp_slice, fls, ilog2};
+use crate::algorithms::{div_rem, div_rem_digit, mac_with_carry, mul3, scalar_mul};

 use UsizePromotion;

@ -2263,7 +2264,7 @@ impl BigUint {
    /// Strips off trailing zero bigdigits - comparisons require the last element in the vector to
    /// be nonzero.
    #[inline]
-    fn normalize(&mut self) {
+    pub(crate) fn normalize(&mut self) {
        while let Some(&0) = self.data.last() {
            self.data.pop();
        }
@ -2271,7 +2272,7 @@ impl BigUint {

    /// Returns a normalized `BigUint`.
    #[inline]
-    fn normalized(mut self) -> BigUint {
+    pub(crate) fn normalized(mut self) -> BigUint {
        self.normalize();
        self
    }
@ -3140,3 +3141,29 @@ fn test_u128_u32_roundtrip() {
        assert_eq!(u32_to_u128(a, b, c, d), *val);
    }
 }
+
+// Mod Inverse
+
+impl<'a> ModInverse<&'a BigUint> for BigUint {
+    fn mod_inverse(self, m: &'a BigUint) -> Option<BigUint> {
+        match mod_inverse(
+            Cow::Owned(BigInt::from_biguint(Plus, self)),
+            &BigInt::from_biguint(Plus, m.clone()),
+        ) {
+            Some(res) => res.to_biguint(),
+            None => None,
+        }
+    }
+}
+
+impl ModInverse<BigUint> for BigUint {
+    fn mod_inverse(self, m: BigUint) -> Option<BigUint> {
+        match mod_inverse(
+            Cow::Owned(BigInt::from_biguint(Plus, self)),
+            &BigInt::from_biguint(Plus, m),
+        ) {
+            Some(res) => res.to_biguint(),
+            None => None,
+        }
+    }
+}
--- a/src/lib.rs
+++ b/src/lib.rs
@ -1,3 +1,7 @@
+// Copyright 2018 Stichting Organism
+//
+// Copyright 2018 Friedel Ziegelmayer
+//
 // Copyright 2013-2014 The Rust Project Developers. See the COPYRIGHT
 // file at the top-level directory of this distribution and at
 // http://rust-lang.org/COPYRIGHT.
@ -82,14 +86,23 @@

 #[cfg(feature = "rand")]
 extern crate rand;
+
 #[cfg(feature = "serde")]
 extern crate serde;

 #[macro_use]
 extern crate smallvec;

+#[cfg(feature = "prime")]
+#[macro_use]
+extern crate lazy_static;
+
 extern crate num_integer as integer;
-extern crate num_traits as traits;
+extern crate num_iter;
+extern crate num_traits;
+
+#[cfg(feature = "prime")]
+extern crate byteorder;

 use std::error::Error;
 use std::fmt;
@ -100,6 +113,12 @@ mod macros;
 mod bigint;
 mod biguint;

+#[cfg(feature = "prime")]
+pub mod prime;
+
+pub mod algorithms;
+pub mod traits;
+
 #[cfg(feature = "rand")]
 mod bigrand;

@ -166,7 +185,7 @@ pub use bigint::Sign;
 pub use bigint::ToBigInt;

 #[cfg(feature = "rand")]
-pub use bigrand::{RandBigInt, RandomBits, UniformBigInt, UniformBigUint};
+pub use bigrand::{RandBigInt, RandPrime, RandomBits, UniformBigInt, UniformBigUint};

 #[cfg(not(feature = "u64_digit"))]
 pub const VEC_SIZE: usize = 8;
--- a/src/macro.rs
+++ b/src/macro.rs
--- a/src/monty.rs
+++ b/src/monty.rs
@ -1,5 +1,5 @@
+use num_traits::{One, Zero};
 use std::ops::Shl;
-use traits::{One, Zero};

 use big_digit::{self, BigDigit, DoubleBigDigit, SignedDoubleBigDigit};
 use biguint::BigUint;
--- a/src/prime.rs
+++ b/src/prime.rs
@ -0,0 +1,558 @@
+//https://github.com/RustCrypto/RSA/blob/master/src/prime.rs
+//! Implements probabilistic prime checkers.
+
+use crate::algorithms::jacobi;
+use bigrand::RandBigInt;
+use byteorder::{BigEndian, ByteOrder};
+use integer::Integer;
+use num_traits::{FromPrimitive, One, ToPrimitive, Zero};
+use rand::rngs::StdRng;
+use rand::SeedableRng;
+use BigInt;
+use BigUint;
+use Sign::Plus;
+
+lazy_static! {
+    pub(crate) static ref BIG_1: BigUint = BigUint::one();
+    pub(crate) static ref BIG_2: BigUint = BigUint::from_u64(2).unwrap();
+    pub(crate) static ref BIG_3: BigUint = BigUint::from_u64(3).unwrap();
+    pub(crate) static ref BIG_64: BigUint = BigUint::from_u64(64).unwrap();
+}
+
+const PRIMES_A: u64 = 3 * 5 * 7 * 11 * 13 * 17 * 19 * 23 * 37;
+const PRIMES_B: u64 = 29 * 31 * 41 * 43 * 47 * 53;
+
+/// Records the primes < 64.
+const PRIME_BIT_MASK: u64 = 1 << 2
+    | 1 << 3
+    | 1 << 5
+    | 1 << 7
+    | 1 << 11
+    | 1 << 13
+    | 1 << 17
+    | 1 << 19
+    | 1 << 23
+    | 1 << 29
+    | 1 << 31
+    | 1 << 37
+    | 1 << 41
+    | 1 << 43
+    | 1 << 47
+    | 1 << 53
+    | 1 << 59
+    | 1 << 61;
+
+/// ProbablyPrime reports whether x is probably prime,
+/// applying the Miller-Rabin test with n pseudorandomly chosen bases
+/// as well as a Baillie-PSW test.
+//
+/// If x is prime, ProbablyPrime returns true.
+/// If x is chosen randomly and not prime, ProbablyPrime probably returns false.
+/// The probability of returning true for a randomly chosen non-prime is at most ¼ⁿ.
+///
+/// ProbablyPrime is 100% accurate for inputs less than 2⁶⁴.
+/// See Menezes et al., Handbook of Applied Cryptography, 1997, pp. 145-149,
+/// and FIPS 186-4 Appendix F for further discussion of the error probabilities.
+///
+/// ProbablyPrime is not suitable for judging primes that an adversary may
+/// have crafted to fool the test.
+///
+/// This is a port of `ProbablyPrime` from the go std lib.
+pub fn probably_prime(x: &BigUint, n: usize) -> bool {
+    if x.is_zero() {
+        return false;
+    }
+
+    if x < &*BIG_64 {
+        return (PRIME_BIT_MASK & (1 << x.to_u64().unwrap())) != 0;
+    }
+
+    if x.is_even() {
+        return false;
+    }
+
+    let r_a = &(x % PRIMES_A);
+    let r_b = &(x % PRIMES_B);
+
+    if (r_a % 3u32).is_zero()
+        || (r_a % 5u32).is_zero()
+        || (r_a % 7u32).is_zero()
+        || (r_a % 11u32).is_zero()
+        || (r_a % 13u32).is_zero()
+        || (r_a % 17u32).is_zero()
+        || (r_a % 19u32).is_zero()
+        || (r_a % 23u32).is_zero()
+        || (r_a % 37u32).is_zero()
+        || (r_b % 29u32).is_zero()
+        || (r_b % 31u32).is_zero()
+        || (r_b % 41u32).is_zero()
+        || (r_b % 43u32).is_zero()
+        || (r_b % 47u32).is_zero()
+        || (r_b % 53u32).is_zero()
+    {
+        return false;
+    }
+
+    probably_prime_miller_rabin(x, n + 1, true) && probably_prime_lucas(x)
+}
+
+/// Reports whether n passes reps rounds of the Miller-Rabin primality test, using pseudo-randomly chosen bases.
+/// If `force2` is true, one of the rounds is forced to use base 2.
+///
+/// See Handbook of Applied Cryptography, p. 139, Algorithm 4.24.
+pub fn probably_prime_miller_rabin(n: &BigUint, reps: usize, force2: bool) -> bool {
+    // println!("miller-rabin: {}", n);
+    let nm1 = n - &*BIG_1;
+    // determine q, k such that nm1 = q << k
+    let k = nm1.trailing_zeros().unwrap();
+    let q = &nm1 >> k;
+    let nm3 = n - &*BIG_3;
+
+    let mut seed_vec = vec![0u8; 8];
+    BigEndian::write_u64(seed_vec.as_mut_slice(), n.get_limb(0));
+    let mut seed = [0u8; 32];
+    seed[0..8].copy_from_slice(&seed_vec[..]);
+    let mut rng = StdRng::from_seed(seed);
+
+    'nextrandom: for i in 0..reps {
+        let x = if i == reps - 1 && force2 {
+            BIG_2.clone()
+        } else {
+            rng.gen_biguint_below(&nm3) + &*BIG_2
+        };
+
+        let mut y = x.modpow(&q, n);
+        if y.is_one() || y == nm1 {
+            continue;
+        }
+
+        for _ in 1..k {
+            y = y.modpow(&*BIG_2, n);
+            if y == nm1 {
+                break 'nextrandom;
+            }
+            if y.is_one() {
+                return false;
+            }
+        }
+        return false;
+    }
+
+    true
+}
+
+/// Reports whether n passes the "almost extra strong" Lucas probable prime test,
+/// using Baillie-OEIS parameter selection. This corresponds to "AESLPSP" on Jacobsen's tables (link below).
+/// The combination of this test and a Miller-Rabin/Fermat test with base 2 gives a Baillie-PSW test.
+///
+///
+/// References:
+///
+/// Baillie and Wagstaff, "Lucas Pseudoprimes", Mathematics of Computation 35(152),
+/// October 1980, pp. 1391-1417, especially page 1401.
+/// http://www.ams.org/journals/mcom/1980-35-152/S0025-5718-1980-0583518-6/S0025-5718-1980-0583518-6.pdf
+///
+/// Grantham, "Frobenius Pseudoprimes", Mathematics of Computation 70(234),
+/// March 2000, pp. 873-891.
+/// http://www.ams.org/journals/mcom/2001-70-234/S0025-5718-00-01197-2/S0025-5718-00-01197-2.pdf
+///
+/// Baillie, "Extra strong Lucas pseudoprimes", OEIS A217719, https://oeis.org/A217719.
+///
+/// Jacobsen, "Pseudoprime Statistics, Tables, and Data", http://ntheory.org/pseudoprimes.html.
+///
+/// Nicely, "The Baillie-PSW Primality Test", http://www.trnicely.net/misc/bpsw.html.
+/// (Note that Nicely's definition of the "extra strong" test gives the wrong Jacobi condition,
+/// as pointed out by Jacobsen.)
+///
+/// Crandall and Pomerance, Prime Numbers: A Computational Perspective, 2nd ed.
+/// Springer, 2005.
+pub fn probably_prime_lucas(n: &BigUint) -> bool {
+    // println!("lucas: {}", n);
+    // Discard 0, 1.
+    if n.is_zero() || n.is_one() {
+        return false;
+    }
+
+    // Two is the only even prime.
+    if n.to_u64() == Some(2) {
+        return false;
+    }
+
+    // Baillie-OEIS "method C" for choosing D, P, Q,
+    // as in https://oeis.org/A217719/a217719.txt:
+    // try increasing P ≥ 3 such that D = P² - 4 (so Q = 1)
+    // until Jacobi(D, n) = -1.
+    // The search is expected to succeed for non-square n after just a few trials.
+    // After more than expected failures, check whether n is square
+    // (which would cause Jacobi(D, n) = 1 for all D not dividing n).
+    let mut p = 3u64;
+    let n_int = BigInt::from_biguint(Plus, n.clone());
+
+    loop {
+        if p > 10000 {
+            // This is widely believed to be impossible.
+            // If we get a report, we'll want the exact number n.
+            panic!("internal error: cannot find (D/n) = -1 for {:?}", n)
+        }
+
+        let j = jacobi(&BigInt::from_u64(p * p - 4).unwrap(), &n_int);
+
+        if j == -1 {
+            break;
+        }
+        if j == 0 {
+            // d = p²-4 = (p-2)(p+2).
+            // If (d/n) == 0 then d shares a prime factor with n.
+            // Since the loop proceeds in increasing p and starts with p-2==1,
+            // the shared prime factor must be p+2.
+            // If p+2 == n, then n is prime; otherwise p+2 is a proper factor of n.
+            return n_int.to_i64() == Some(p as i64 + 2);
+        }
+
+        // We'll never find (d/n) = -1 if n is a square.
+        // If n is a non-square we expect to find a d in just a few attempts on average.
+        // After 40 attempts, take a moment to check if n is indeed a square.
+        if p == 40 && (&n_int * &n_int).sqrt() == n_int {
+            return false;
+        }
+
+        p += 1;
+    }
+
+    // Grantham definition of "extra strong Lucas pseudoprime", after Thm 2.3 on p. 876
+    // (D, P, Q above have become Δ, b, 1):
+    //
+    // Let U_n = U_n(b, 1), V_n = V_n(b, 1), and Δ = b²-4.
+    // An extra strong Lucas pseudoprime to base b is a composite n = 2^r s + Jacobi(Δ, n),
+    // where s is odd and gcd(n, 2*Δ) = 1, such that either (i) U_s ≡ 0 mod n and V_s ≡ ±2 mod n,
+    // or (ii) V_{2^t s} ≡ 0 mod n for some 0 ≤ t < r-1.
+    //
+    // We know gcd(n, Δ) = 1 or else we'd have found Jacobi(d, n) == 0 above.
+    // We know gcd(n, 2) = 1 because n is odd.
+    //
+    // Arrange s = (n - Jacobi(Δ, n)) / 2^r = (n+1) / 2^r.
+    let mut s = n + &*BIG_1;
+    let r = s.trailing_zeros().unwrap();
+    s = &s >> r;
+    let nm2 = n - &*BIG_2; // n - 2
+
+    // We apply the "almost extra strong" test, which checks the above conditions
+    // except for U_s ≡ 0 mod n, which allows us to avoid computing any U_k values.
+    // Jacobsen points out that maybe we should just do the full extra strong test:
+    // "It is also possible to recover U_n using Crandall and Pomerance equation 3.13:
+    // U_n = D^-1 (2V_{n+1} - PV_n) allowing us to run the full extra-strong test
+    // at the cost of a single modular inversion. This computation is easy and fast in GMP,
+    // so we can get the full extra-strong test at essentially the same performance as the
+    // almost extra strong test."
+
+    // Compute Lucas sequence V_s(b, 1), where:
+    //
+    //	V(0) = 2
+    //	V(1) = P
+    //	V(k) = P V(k-1) - Q V(k-2).
+    //
+    // (Remember that due to method C above, P = b, Q = 1.)
+    //
+    // In general V(k) = α^k + β^k, where α and β are roots of x² - Px + Q.
+    // Crandall and Pomerance (p.147) observe that for 0 ≤ j ≤ k,
+    //
+    //	V(j+k) = V(j)V(k) - V(k-j).
+    //
+    // So in particular, to quickly double the subscript:
+    //
+    //	V(2k) = V(k)² - 2
+    //	V(2k+1) = V(k) V(k+1) - P
+    //
+    // We can therefore start with k=0 and build up to k=s in log₂(s) steps.
+    let mut vk = BIG_2.clone();
+    let mut vk1 = BigUint::from_u64(p).unwrap();
+
+    for i in (0..s.bits()).rev() {
+        if is_bit_set(&s, i) {
+            // k' = 2k+1
+            // V(k') = V(2k+1) = V(k) V(k+1) - P
+            let t1 = (&vk * &vk1) + n - p;
+            vk = &t1 % n;
+            // V(k'+1) = V(2k+2) = V(k+1)² - 2
+            let t1 = (&vk1 * &vk1) + &nm2;
+            vk1 = &t1 % n;
+        } else {
+            // k' = 2k
+            // V(k'+1) = V(2k+1) = V(k) V(k+1) - P
+            let t1 = (&vk * &vk1) + n - p;
+            vk1 = &t1 % n;
+            // V(k') = V(2k) = V(k)² - 2
+            let t1 = (&vk * &vk) + &nm2;
+            vk = &t1 % n;
+        }
+    }
+
+    // Now k=s, so vk = V(s). Check V(s) ≡ ±2 (mod n).
+    if vk.to_u64() == Some(2) || vk == nm2 {
+        // Check U(s) ≡ 0.
+        // As suggested by Jacobsen, apply Crandall and Pomerance equation 3.13:
+        //
+        //	U(k) = D⁻¹ (2 V(k+1) - P V(k))
+        //
+        // Since we are checking for U(k) == 0 it suffices to check 2 V(k+1) == P V(k) mod n,
+        // or P V(k) - 2 V(k+1) == 0 mod n.
+        let mut t1 = &vk * p;
+        let mut t2 = &vk1 << 1;
+
+        if t1 < t2 {
+            ::std::mem::swap(&mut t1, &mut t2);
+        }
+
+        t1 -= t2;
+
+        if (t1 % n).is_zero() {
+            return true;
+        }
+    }
+
+    // Check V(2^t s) ≡ 0 mod n for some 0 ≤ t < r-1.
+    for _ in 0..r - 1 {
+        if vk.is_zero() {
+            return true;
+        }
+
+        // Optimization: V(k) = 2 is a fixed point for V(k') = V(k)² - 2,
+        // so if V(k) = 2, we can stop: we will never find a future V(k) == 0.
+        if vk.to_u64() == Some(2) {
+            return false;
+        }
+
+        // k' = 2k
+        // V(k') = V(2k) = V(k)² - 2
+        let t1 = (&vk * &vk) - &*BIG_2;
+        vk = &t1 % n;
+    }
+
+    false
+}
+
+/// Checks if the i-th bit is set
+#[inline]
+fn is_bit_set(x: &BigUint, i: usize) -> bool {
+    get_bit(x, i) == 1
+}
+
+/// Returns the i-th bit.
+#[inline]
+fn get_bit(x: &BigUint, i: usize) -> u8 {
+    let j = i / 64;
+    // if is out of range of the set words, it is always false.
+    if i >= x.bits() {
+        return 0;
+    }
+
+    (x.get_limb(j) >> (i % 64) & 1) as u8
+}
+
+// pub fn big_prime(size: uint) -> BigUint {
+//   let one: BigUint = One::one();
+//   let two = one + one;
+
+//   let mut rng = task_rng();
+//   let mut candidate = rng.gen_biguint(size);
+//   if candidate.is_even() {
+//     candidate = candidate + one;
+//   }
+//   while !is_prime(&candidate) {
+//     candidate = candidate + two;
+//   }
+//   candidate
+// }
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+    // use RandBigInt;
+
+    lazy_static! {
+        static ref PRIMES: Vec<&'static str> = vec![
+        "2",
+        "3",
+        "5",
+        "7",
+        "11",
+
+        "13756265695458089029",
+        "13496181268022124907",
+        "10953742525620032441",
+        "17908251027575790097",
+
+        // https://golang.org/issue/638
+        "18699199384836356663",
+
+        "98920366548084643601728869055592650835572950932266967461790948584315647051443",
+        "94560208308847015747498523884063394671606671904944666360068158221458669711639",
+
+        // http://primes.utm.edu/lists/small/small3.html
+        "449417999055441493994709297093108513015373787049558499205492347871729927573118262811508386655998299074566974373711472560655026288668094291699357843464363003144674940345912431129144354948751003607115263071543163",
+        "230975859993204150666423538988557839555560243929065415434980904258310530753006723857139742334640122533598517597674807096648905501653461687601339782814316124971547968912893214002992086353183070342498989426570593",
+        "5521712099665906221540423207019333379125265462121169655563495403888449493493629943498064604536961775110765377745550377067893607246020694972959780839151452457728855382113555867743022746090187341871655890805971735385789993",
+        "203956878356401977405765866929034577280193993314348263094772646453283062722701277632936616063144088173312372882677123879538709400158306567338328279154499698366071906766440037074217117805690872792848149112022286332144876183376326512083574821647933992961249917319836219304274280243803104015000563790123",
+        // ECC primes: http://tools.ietf.org/html/draft-ladd-safecurves-02
+        "3618502788666131106986593281521497120414687020801267626233049500247285301239",                                                                                  // Curve1174: 2^251-9
+        "57896044618658097711785492504343953926634992332820282019728792003956564819949",                                                                                 // Curve25519: 2^255-19
+        "9850501549098619803069760025035903451269934817616361666987073351061430442874302652853566563721228910201656997576599",                                           // E-382: 2^382-105
+        "42307582002575910332922579714097346549017899709713998034217522897561970639123926132812109468141778230245837569601494931472367",                                 // Curve41417: 2^414-17
+        "6864797660130609714981900799081393217269435300143305409394463459185543183397656052122559640661454554977296311391480858037121987999716643812574028291115057151", // E-521: 2^521-1
+        ];
+
+        static ref COMPOSITES: Vec<&'static str> = vec![
+            "0",
+            "1",
+
+            "21284175091214687912771199898307297748211672914763848041968395774954376176754",
+            "6084766654921918907427900243509372380954290099172559290432744450051395395951",
+            "84594350493221918389213352992032324280367711247940675652888030554255915464401",
+            "82793403787388584738507275144194252681",
+
+            // Arnault, "Rabin-Miller Primality Test: Composite Numbers Which Pass It",
+            // Mathematics of Computation, 64(209) (January 1995), pp. 335-361.
+            "1195068768795265792518361315725116351898245581", // strong pseudoprime to prime bases 2 through 29
+            // strong pseudoprime to all prime bases up to 200
+            "8038374574536394912570796143419421081388376882875581458374889175222974273765333652186502336163960045457915042023603208766569966760987284043965408232928738791850869166857328267761771029389697739470167082304286871099974399765441448453411558724506334092790222752962294149842306881685404326457534018329786111298960644845216191652872597534901",
+
+            // Extra-strong Lucas pseudoprimes. https://oeis.org/A217719
+            "989",
+            "3239",
+            "5777",
+            "10877",
+            "27971",
+            "29681",
+            "30739",
+            "31631",
+            "39059",
+            "72389",
+            "73919",
+            "75077",
+            "100127",
+            "113573",
+            "125249",
+            "137549",
+            "137801",
+            "153931",
+            "155819",
+            "161027",
+            "162133",
+            "189419",
+            "218321",
+            "231703",
+            "249331",
+            "370229",
+            "429479",
+            "430127",
+            "459191",
+            "473891",
+            "480689",
+            "600059",
+            "621781",
+            "632249",
+            "635627",
+
+            "3673744903",
+            "3281593591",
+            "2385076987",
+            "2738053141",
+            "2009621503",
+            "1502682721",
+            "255866131",
+            "117987841",
+            "587861",
+
+            "6368689",
+            "8725753",
+            "80579735209",
+            "105919633",
+        ];
+    }
+
+    #[test]
+    fn test_primes() {
+        for prime in PRIMES.iter() {
+            let p = BigUint::parse_bytes(prime.as_bytes(), 10).unwrap();
+            for i in [0, 1, 20].iter() {
+                assert!(
+                    probably_prime(&p, *i as usize),
+                    "{} is a prime ({})",
+                    prime,
+                    i,
+                );
+            }
+        }
+    }
+
+    #[test]
+    fn test_composites() {
+        for comp in COMPOSITES.iter() {
+            let p = BigUint::parse_bytes(comp.as_bytes(), 10).unwrap();
+            for i in [0, 1, 20].iter() {
+                assert!(
+                    !probably_prime(&p, *i as usize),
+                    "{} is a composite ({})",
+                    comp,
+                    i,
+                );
+            }
+        }
+    }
+
+    macro_rules! test_pseudo_primes {
+        ($name:ident, $cond:expr, $want:expr) => {
+            #[test]
+            fn $name() {
+                let mut i = 3;
+                let mut want = $want;
+                while i < 100000 {
+                    let n = BigUint::from_u64(i).unwrap();
+                    let pseudo = $cond(&n);
+                    if pseudo && (want.is_empty() || i != want[0]) {
+                        panic!("cond({}) = true, want false", i);
+                    } else if !pseudo && !want.is_empty() && i == want[0] {
+                        panic!("cond({}) = false, want true", i);
+                    }
+                    if !want.is_empty() && i == want[0] {
+                        want = want[1..].to_vec();
+                    }
+                    i += 2;
+                }
+
+                if !want.is_empty() {
+                    panic!("forgot to test: {:?}", want);
+                }
+            }
+        };
+    }
+
+    test_pseudo_primes!(
+        test_probably_prime_miller_rabin,
+        |n| probably_prime_miller_rabin(n, 1, true) && !probably_prime_lucas(n),
+        vec![
+            2047, 3277, 4033, 4681, 8321, 15841, 29341, 42799, 49141, 52633, 65281, 74665, 80581,
+            85489, 88357, 90751,
+        ]
+    );
+
+    test_pseudo_primes!(
+        test_probably_prime_lucas,
+        |n| probably_prime_lucas(n) && !probably_prime_miller_rabin(n, 1, true),
+        vec![989, 3239, 5777, 10877, 27971, 29681, 30739, 31631, 39059, 72389, 73919, 75077,]
+    );
+
+    #[test]
+    fn test_bit_set() {
+        let v = &vec![0b10101001];
+        let num = BigUint::from_slice(&v);
+        assert!(is_bit_set(&num, 0));
+        assert!(!is_bit_set(&num, 1));
+        assert!(!is_bit_set(&num, 2));
+        assert!(is_bit_set(&num, 3));
+        assert!(!is_bit_set(&num, 4));
+        assert!(is_bit_set(&num, 5));
+        assert!(!is_bit_set(&num, 6));
+        assert!(is_bit_set(&num, 7));
+    }
+}
--- a/src/traits.rs
+++ b/src/traits.rs
@ -0,0 +1,10 @@
+/// Generic trait to implement modular inverse
+pub trait ModInverse<R: Sized>: Sized {
+    /// Function to calculate the [modular multiplicative
+    /// inverse](https://en.wikipedia.org/wiki/Modular_multiplicative_inverse) of an integer *a* modulo *m*.
+    ///
+    /// TODO: references
+    /// Returns the modular inverse of `self`.
+    /// If none exists it returns `None`.
+    fn mod_inverse(self, m: R) -> Option<Self>;
+}
--- a/tests/biguint.rs
+++ b/tests/biguint.rs
@ -1,7 +1,6 @@
 extern crate num_bigint_dig as num_bigint;
 extern crate num_integer;
 extern crate num_traits;
-#[macro_use]
 extern crate smallvec;

 use num_bigint::Sign::Plus;
--- a/tests/rand.rs
+++ b/tests/rand.rs
@ -419,3 +419,27 @@ mod bigint {
        seeded_value_stability::<XorShiftRng>(EXPECTED_XOR);
    }
 }
+
+mod prime {
+    use num_bigint::prime::probably_prime;
+    use num_bigint::RandPrime;
+    use rand::prelude::*;
+
+    #[test]
+    fn test_prime_small() {
+        let mut rng = StdRng::from_seed([0u8; 32]);
+        for n in 2..10 {
+            let p = rng.gen_prime(n);
+
+            assert_eq!(p.bits(), n);
+            assert!(probably_prime(&p, 32));
+        }
+    }
+
+    #[test]
+    fn test_gen_prime_1024() {
+        let mut rng = StdRng::from_seed([0u8; 32]);
+        let p = rng.gen_prime(1024);
+        assert_eq!(p.bits(), 1024);
+    }
+}