From f523b9c359b0b4e8718176aba693c6e061cd0353 Mon Sep 17 00:00:00 2001
From: str4d <str4d@mail.i2p>
Date: Sat, 7 Oct 2017 22:29:19 +0100
Subject: [PATCH 01/14] Implement modpow() for BigUint backed by Montgomery
 Multiplication

Based on this Gist: https://gist.github.com/yshui/027eecdf95248ea69606

Closes #136
---
 benches/bigint.rs           |  25 +++++-
 bigint/src/biguint.rs       |   9 +++
 bigint/src/monty.rs         | 155 ++++++++++++++++++++++++++++++++++++
 bigint/src/tests/biguint.rs |  73 +++++++++++++++++
 4 files changed, 261 insertions(+), 1 deletion(-)
 create mode 100644 bigint/src/monty.rs

diff --git a/benches/bigint.rs b/benches/bigint.rs
index 0eafb93..9e18282 100644
--- a/benches/bigint.rs
+++ b/benches/bigint.rs
@@ -6,7 +6,7 @@ extern crate rand;
 
 use std::mem::replace;
 use test::Bencher;
-use num::{BigInt, BigUint, Zero, One, FromPrimitive};
+use num::{BigInt, BigUint, Zero, One, FromPrimitive, Num};
 use num::bigint::RandBigInt;
 use rand::{SeedableRng, StdRng};
 
@@ -255,3 +255,26 @@ fn pow_bench(b: &mut Bencher) {
         }
     });
 }
+
+#[bench]
+fn modpow(b: &mut Bencher) {
+    let mut rng = get_rng();
+    let base = rng.gen_biguint(2048);
+    let e = rng.gen_biguint(2048);
+    // This modulus is the prime from the 2048-bit MODP DH group:
+    // https://tools.ietf.org/html/rfc3526#section-3
+    let m = BigUint::from_str_radix("\
+        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", 16).unwrap();
+
+    b.iter(|| base.modpow(&e, &m));
+}
diff --git a/bigint/src/biguint.rs b/bigint/src/biguint.rs
index 367e66e..f3a1c23 100644
--- a/bigint/src/biguint.rs
+++ b/bigint/src/biguint.rs
@@ -21,6 +21,8 @@ use traits::{ToPrimitive, FromPrimitive, Float, Num, Unsigned, CheckedAdd, Check
 
 #[path = "algorithms.rs"]
 mod algorithms;
+#[path = "monty.rs"]
+mod monty;
 pub use self::algorithms::big_digit;
 pub use self::big_digit::{BigDigit, DoubleBigDigit, ZERO_BIG_DIGIT};
 
@@ -28,6 +30,7 @@ use self::algorithms::{mac_with_carry, mul3, scalar_mul, div_rem, div_rem_digit}
 use self::algorithms::{__add2, add2, sub2, sub2rev};
 use self::algorithms::{biguint_shl, biguint_shr};
 use self::algorithms::{cmp_slice, fls, ilog2};
+use self::monty::{MontyReducer, monty_modpow};
 
 use UsizePromotion;
 
@@ -1611,6 +1614,12 @@ impl BigUint {
         self.normalize();
         self
     }
+
+    /// Returns `(self ^ exponent) % modulus`.
+    pub fn modpow(&self, exponent: &Self, modulus: &Self) -> Self {
+        let mr = MontyReducer::new(modulus);
+        monty_modpow(self, exponent, &mr)
+    }
 }
 
 #[cfg(feature = "serde")]
diff --git a/bigint/src/monty.rs b/bigint/src/monty.rs
new file mode 100644
index 0000000..422201e
--- /dev/null
+++ b/bigint/src/monty.rs
@@ -0,0 +1,155 @@
+use std::iter::repeat;
+use integer::Integer;
+use traits::{Zero, One};
+
+use biguint::BigUint;
+
+pub struct MontyReducer<'a> {
+    p: &'a BigUint,
+    n: Vec<u32>,
+    n0inv: u64
+}
+
+// Calculate the modular inverse of `num`, using Extended GCD.
+//
+// Reference:
+// Brent & Zimmermann, Modern Computer Arithmetic, v0.5.9, Algorithm 1.20
+fn inv_mod_u32(num: u32) -> u64 {
+    // num needs to be relatively prime to u32::max_value()
+    assert!(num % 2 != 0);
+
+    let mut a: i64 = num as i64;
+    let mut b: i64 = (u32::max_value() as i64) + 1;
+    let mu = b;
+
+    // ExtendedGcd
+    // Input: positive integers a and b
+    // Output: integers (g, u, v) such that g = gcd(a, b) = ua + vb
+    // As we don't need v for modular inverse, we don't calculate it.
+
+    // 1: (u, w) <- (1, 0)
+    let mut u = 1;
+    let mut w = 0;
+    // 3: while b != 0
+    while b != 0 {
+        // 4: (q, r) <- DivRem(a, b)
+        let q = a / b;
+        let r = a % b;
+        // 5: (a, b) <- (b, r)
+        a = b; b = r;
+        // 6: (u, w) <- (w, u - qw)
+        let m = u - w*q;
+        u = w; w = m;
+    }
+
+    assert!(a == 1);
+    // Ensure returned value is in-range
+    if u < 0 {
+        (u + mu) as u64
+    } else {
+        u as u64
+    }
+}
+
+impl<'a> MontyReducer<'a> {
+    pub fn new(p: &'a BigUint) -> Self {
+        let n : Vec<u32> = p.data.clone();
+        let n0inv = inv_mod_u32(n[0]);
+        MontyReducer { p: p, n: n, n0inv: n0inv }
+    }
+}
+
+// Montgomery Reduction
+//
+// Reference:
+// Brent & Zimmermann, Modern Computer Arithmetic, v0.5.9, Algorithm 2.6
+pub fn monty_redc(a: BigUint, mr: &MontyReducer) -> BigUint {
+    let mut c = a.data;
+    let n = &mr.n;
+    let n_size = n.len();
+    let old_size = c.len();
+
+    // Allocate sufficient work space
+    c.reserve(2*n_size+2-old_size);
+    c.extend(repeat(0).take(2*n_size+2-old_size));
+
+    // β is the size of a word, in this case 32 bits. So "a mod β" is
+    // equivalent to masking a to 32 bits.
+    let beta_mask = u32::max_value() as u64;
+    // mu <- -N^(-1) mod β
+    let mu = (beta_mask-mr.n0inv)+1;
+
+    // 1: for i = 0 to (n-1)
+    for i in 0..n_size {
+        // Carry storage
+        let mut carry = 0;
+
+        // 2: q_i <- mu*c_i mod β
+        let q_i = ((c[i] as u64) * mu) & beta_mask;
+
+        // 3: C <- C + q_i * N * β^i
+        // When iterating over each word, this becomes:
+        for j in 0..n_size {
+            // c_(i+j) <- c_(i+j) + q_i * n_j
+            let x = (c[i+j] as u64) + q_i * (n[j] as u64) + carry;
+            c[i+j] = (x & beta_mask) as u32;
+            carry = x >> 32;
+        }
+
+        // Apply the remaining carry to the rest of the work space
+        for j in n_size..2*n_size-i+2 {
+            let x = (c[i+j] as u64) + carry;
+            c[i+j] = (x & beta_mask) as u32;
+            carry = x >> 32;
+        }
+    }
+
+    // 4: R <- C * β^(-n)
+    // This is an n-word bitshift, equivalent to skipping n words.
+    let r : Vec<u32> = c.iter().skip(n_size).cloned().collect();
+    let ret = BigUint::new(r);
+
+    // 5: if R >= β^n then return R-N else return R.
+    if &ret < mr.p {
+        ret
+    } else {
+        &ret-mr.p
+    }
+}
+
+// Montgomery Multiplication
+fn monty_mult(a: BigUint, b: &BigUint, mr: &MontyReducer) -> BigUint {
+    monty_redc(a * b, mr)
+}
+
+// Montgomery Squaring
+fn monty_sqr(a: BigUint, mr: &MontyReducer) -> BigUint {
+    // TODO: Replace with an optimised squaring function
+    monty_redc(&a * &a, mr)
+}
+
+pub fn monty_modpow(a: &BigUint, exp: &BigUint, mr: &MontyReducer) -> BigUint{
+    // Calculate the Montgomery parameter
+    let mut r : BigUint = One::one();
+    while &r < mr.p {
+        r = r << 32;
+    }
+
+    // Map the base to the Montgomery domain
+    let mut apri = a * &r % mr.p;
+
+    // Binary exponentiation
+    let mut ans = &r % mr.p;
+    let mut e = exp.clone();
+    let zero = Zero::zero();
+    while e > zero {
+        if e.is_odd() {
+            ans = monty_mult(ans, &apri, mr);
+        }
+        apri = monty_sqr(apri, mr);
+        e = e >> 1;
+    }
+
+    // Map the result back to the residues domain
+    monty_redc(ans, mr)
+}
diff --git a/bigint/src/tests/biguint.rs b/bigint/src/tests/biguint.rs
index d010bb5..85abb3b 100644
--- a/bigint/src/tests/biguint.rs
+++ b/bigint/src/tests/biguint.rs
@@ -1089,6 +1089,79 @@ fn test_is_even() {
     assert!(((&one << 64) + one).is_odd());
 }
 
+#[test]
+fn test_modpow() {
+    fn check(b: usize, e: usize, m: usize, r: usize) {
+        let big_b: BigUint = FromPrimitive::from_usize(b).unwrap();
+        let big_e: BigUint = FromPrimitive::from_usize(e).unwrap();
+        let big_m: BigUint = FromPrimitive::from_usize(m).unwrap();
+        let big_r: BigUint = FromPrimitive::from_usize(r).unwrap();
+
+        assert_eq!(big_b.modpow(&big_e, &big_m), big_r);
+    }
+
+    check(1, 0, 11, 1);
+    check(0, 15, 11, 0);
+    check(3, 7, 11, 9);
+    check(5, 117, 19, 1);
+}
+
+#[test]
+fn test_modpow_big() {
+    let b = BigUint::from_str_radix("\
+        efac3c0a_0de55551_fee0bfe4_67fa017a_1a898fa1_6ca57cb1\
+        ca9e3248_cacc09a9_b99d6abc_38418d0f_82ae4238_d9a68832\
+        aadec7c1_ac5fed48_7a56a71b_67ac59d5_afb28022_20d9592d\
+        247c4efc_abbd9b75_586088ee_1dc00dc4_232a8e15_6e8191dd\
+        675b6ae0_c80f5164_752940bc_284b7cee_885c1e10_e495345b\
+        8fbe9cfd_e5233fe1_19459d0b_d64be53c_27de5a02_a829976b\
+        33096862_82dad291_bd38b6a9_be396646_ddaf8039_a2573c39\
+        1b14e8bc_2cb53e48_298c047e_d9879e9c_5a521076_f0e27df3\
+        990e1659_d3d8205b_6443ebc0_9918ebee_6764f668_9f2b2be3\
+        b59cbc76_d76d0dfc_d737c3ec_0ccf9c00_ad0554bf_17e776ad\
+        b4edf9cc_6ce540be_76229093_5c53893b", 16).unwrap();
+    let e = BigUint::from_str_radix("\
+        be0e6ea6_08746133_e0fbc1bf_82dba91e_e2b56231_a81888d2\
+        a833a1fc_f7ff002a_3c486a13_4f420bf3_a5435be9_1a5c8391\
+        774d6e6c_085d8357_b0c97d4d_2bb33f7c_34c68059_f78d2541\
+        eacc8832_426f1816_d3be001e_b69f9242_51c7708e_e10efe98\
+        449c9a4a_b55a0f23_9d797410_515da00d_3ea07970_4478a2ca\
+        c3d5043c_bd9be1b4_6dce479d_4302d344_84a939e6_0ab5ada7\
+        12ae34b2_30cc473c_9f8ee69d_2cac5970_29f5bf18_bc8203e4\
+        f3e895a2_13c94f1e_24c73d77_e517e801_53661fdd_a2ce9e47\
+        a73dd7f8_2f2adb1e_3f136bf7_8ae5f3b8_08730de1_a4eff678\
+        e77a06d0_19a522eb_cbefba2a_9caf7736_b157c5c6_2d192591\
+        17946850_2ddb1822_117b68a0_32f7db88", 16).unwrap();
+    // This modulus is the prime from the 2048-bit MODP DH group:
+    // https://tools.ietf.org/html/rfc3526#section-3
+    let m = BigUint::from_str_radix("\
+        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", 16).unwrap();
+    let r = BigUint::from_str_radix("\
+        a1468311_6e56edc9_7a98228b_5e924776_0dd7836e_caabac13\
+        eda5373b_4752aa65_a1454850_40dc770e_30aa8675_6be7d3a8\
+        9d3085e4_da5155cf_b451ef62_54d0da61_cf2b2c87_f495e096\
+        055309f7_77802bbb_37271ba8_1313f1b5_075c75d1_024b6c77\
+        fdb56f17_b05bce61_e527ebfd_2ee86860_e9907066_edd526e7\
+        93d289bf_6726b293_41b0de24_eff82424_8dfd374b_4ec59542\
+        35ced2b2_6b195c90_10042ffb_8f58ce21_bc10ec42_64fda779\
+        d352d234_3d4eaea6_a86111ad_a37e9555_43ca78ce_2885bed7\
+        5a30d182_f1cf6834_dc5b6e27_1a41ac34_a2e91e11_33363ff0\
+        f88a7b04_900227c9_f6e6d06b_7856b4bb_4e354d61_060db6c8\
+        109c4735_6e7db425_7b5d74c7_0b709508", 16).unwrap();
+
+    assert_eq!(b.modpow(&e, &m), r);
+}
+
 fn to_str_pairs() -> Vec<(BigUint, Vec<(u32, String)>)> {
     let bits = big_digit::BITS;
     vec![(Zero::zero(),

From 720893f67bddc1f9c17b682ae3ebdc3e62fd173a Mon Sep 17 00:00:00 2001
From: str4d <str4d@mail.i2p>
Date: Mon, 9 Oct 2017 16:11:18 +0100
Subject: [PATCH 02/14] Add support to BigUint.from_str_radix() for using _ as
 a visual separator

---
 bigint/src/biguint.rs | 1 +
 1 file changed, 1 insertion(+)

diff --git a/bigint/src/biguint.rs b/bigint/src/biguint.rs
index f3a1c23..1c609c2 100644
--- a/bigint/src/biguint.rs
+++ b/bigint/src/biguint.rs
@@ -243,6 +243,7 @@ impl Num for BigUint {
                 b'0'...b'9' => b - b'0',
                 b'a'...b'z' => b - b'a' + 10,
                 b'A'...b'Z' => b - b'A' + 10,
+                b'_' => continue,
                 _ => u8::MAX,
             };
             if d < radix as u8 {

From 2a1fe6e7ef7aaee8039894b34d1bf6850b8f99f9 Mon Sep 17 00:00:00 2001
From: Josh Stone <cuviper@gmail.com>
Date: Sun, 22 Oct 2017 14:57:52 -0700
Subject: [PATCH 03/14] bigint: fix parsing leading _ and test more

---
 bigint/src/biguint.rs       | 7 +++++++
 bigint/src/tests/biguint.rs | 6 ++++++
 2 files changed, 13 insertions(+)

diff --git a/bigint/src/biguint.rs b/bigint/src/biguint.rs
index 1c609c2..82e1483 100644
--- a/bigint/src/biguint.rs
+++ b/bigint/src/biguint.rs
@@ -236,6 +236,13 @@ impl Num for BigUint {
             return Err(e.into());
         }
 
+        if s.starts_with('_') {
+            // Must lead with a real digit!
+            // create ParseIntError::InvalidDigit
+            let e = u64::from_str_radix(s, radix).unwrap_err();
+            return Err(e.into());
+        }
+
         // First normalize all characters to plain digit values
         let mut v = Vec::with_capacity(s.len());
         for b in s.bytes() {
diff --git a/bigint/src/tests/biguint.rs b/bigint/src/tests/biguint.rs
index 85abb3b..89101a6 100644
--- a/bigint/src/tests/biguint.rs
+++ b/bigint/src/tests/biguint.rs
@@ -1541,6 +1541,8 @@ fn test_from_str_radix() {
     assert_eq!(zed, None);
     let blank = BigUint::from_str_radix("_", 2).ok();
     assert_eq!(blank, None);
+    let blank_one = BigUint::from_str_radix("_1", 2).ok();
+    assert_eq!(blank_one, None);
     let plus_one = BigUint::from_str_radix("+1", 10).ok();
     assert_eq!(plus_one, Some(BigUint::from_slice(&[1])));
     let plus_plus_one = BigUint::from_str_radix("++1", 10).ok();
@@ -1549,6 +1551,10 @@ fn test_from_str_radix() {
     assert_eq!(minus_one, None);
     let zero_plus_two = BigUint::from_str_radix("0+2", 10).ok();
     assert_eq!(zero_plus_two, None);
+    let three = BigUint::from_str_radix("1_1", 2).ok();
+    assert_eq!(three, Some(BigUint::from_slice(&[3])));
+    let ff = BigUint::from_str_radix("1111_1111", 2).ok();
+    assert_eq!(ff, Some(BigUint::from_slice(&[0xff])));
 }
 
 #[test]

From c2fba067873c6bfaaade7bf6db446006a391cd43 Mon Sep 17 00:00:00 2001
From: Josh Stone <cuviper@gmail.com>
Date: Sun, 22 Oct 2017 15:05:16 -0700
Subject: [PATCH 04/14] bigint: less pub in monty

---
 bigint/src/biguint.rs |  5 ++---
 bigint/src/monty.rs   | 16 +++++++++-------
 2 files changed, 11 insertions(+), 10 deletions(-)

diff --git a/bigint/src/biguint.rs b/bigint/src/biguint.rs
index 82e1483..0bfc158 100644
--- a/bigint/src/biguint.rs
+++ b/bigint/src/biguint.rs
@@ -30,7 +30,7 @@ use self::algorithms::{mac_with_carry, mul3, scalar_mul, div_rem, div_rem_digit}
 use self::algorithms::{__add2, add2, sub2, sub2rev};
 use self::algorithms::{biguint_shl, biguint_shr};
 use self::algorithms::{cmp_slice, fls, ilog2};
-use self::monty::{MontyReducer, monty_modpow};
+use self::monty::monty_modpow;
 
 use UsizePromotion;
 
@@ -1625,8 +1625,7 @@ impl BigUint {
 
     /// Returns `(self ^ exponent) % modulus`.
     pub fn modpow(&self, exponent: &Self, modulus: &Self) -> Self {
-        let mr = MontyReducer::new(modulus);
-        monty_modpow(self, exponent, &mr)
+        monty_modpow(self, exponent, modulus)
     }
 }
 
diff --git a/bigint/src/monty.rs b/bigint/src/monty.rs
index 422201e..a9ef305 100644
--- a/bigint/src/monty.rs
+++ b/bigint/src/monty.rs
@@ -4,7 +4,7 @@ use traits::{Zero, One};
 
 use biguint::BigUint;
 
-pub struct MontyReducer<'a> {
+struct MontyReducer<'a> {
     p: &'a BigUint,
     n: Vec<u32>,
     n0inv: u64
@@ -52,7 +52,7 @@ fn inv_mod_u32(num: u32) -> u64 {
 }
 
 impl<'a> MontyReducer<'a> {
-    pub fn new(p: &'a BigUint) -> Self {
+    fn new(p: &'a BigUint) -> Self {
         let n : Vec<u32> = p.data.clone();
         let n0inv = inv_mod_u32(n[0]);
         MontyReducer { p: p, n: n, n0inv: n0inv }
@@ -63,7 +63,7 @@ impl<'a> MontyReducer<'a> {
 //
 // Reference:
 // Brent & Zimmermann, Modern Computer Arithmetic, v0.5.9, Algorithm 2.6
-pub fn monty_redc(a: BigUint, mr: &MontyReducer) -> BigUint {
+fn monty_redc(a: BigUint, mr: &MontyReducer) -> BigUint {
     let mut c = a.data;
     let n = &mr.n;
     let n_size = n.len();
@@ -128,7 +128,9 @@ fn monty_sqr(a: BigUint, mr: &MontyReducer) -> BigUint {
     monty_redc(&a * &a, mr)
 }
 
-pub fn monty_modpow(a: &BigUint, exp: &BigUint, mr: &MontyReducer) -> BigUint{
+pub fn monty_modpow(a: &BigUint, exp: &BigUint, modulus: &BigUint) -> BigUint{
+    let mr = MontyReducer::new(modulus);
+
     // Calculate the Montgomery parameter
     let mut r : BigUint = One::one();
     while &r < mr.p {
@@ -144,12 +146,12 @@ pub fn monty_modpow(a: &BigUint, exp: &BigUint, mr: &MontyReducer) -> BigUint{
     let zero = Zero::zero();
     while e > zero {
         if e.is_odd() {
-            ans = monty_mult(ans, &apri, mr);
+            ans = monty_mult(ans, &apri, &mr);
         }
-        apri = monty_sqr(apri, mr);
+        apri = monty_sqr(apri, &mr);
         e = e >> 1;
     }
 
     // Map the result back to the residues domain
-    monty_redc(ans, mr)
+    monty_redc(ans, &mr)
 }

From aea5f8521699e4b89e50f490441ad381f8a888e2 Mon Sep 17 00:00:00 2001
From: Josh Stone <cuviper@gmail.com>
Date: Sun, 22 Oct 2017 15:15:02 -0700
Subject: [PATCH 05/14] bigint::monty: store the inverse as u32

---
 bigint/src/monty.rs | 17 ++++++-----------
 1 file changed, 6 insertions(+), 11 deletions(-)

diff --git a/bigint/src/monty.rs b/bigint/src/monty.rs
index a9ef305..8960c70 100644
--- a/bigint/src/monty.rs
+++ b/bigint/src/monty.rs
@@ -7,20 +7,19 @@ use biguint::BigUint;
 struct MontyReducer<'a> {
     p: &'a BigUint,
     n: Vec<u32>,
-    n0inv: u64
+    n0inv: u32
 }
 
 // Calculate the modular inverse of `num`, using Extended GCD.
 //
 // Reference:
 // Brent & Zimmermann, Modern Computer Arithmetic, v0.5.9, Algorithm 1.20
-fn inv_mod_u32(num: u32) -> u64 {
-    // num needs to be relatively prime to u32::max_value()
+fn inv_mod_u32(num: u32) -> u32 {
+    // num needs to be relatively prime to 2**32 -- i.e. it must be odd.
     assert!(num % 2 != 0);
 
     let mut a: i64 = num as i64;
     let mut b: i64 = (u32::max_value() as i64) + 1;
-    let mu = b;
 
     // ExtendedGcd
     // Input: positive integers a and b
@@ -43,12 +42,8 @@ fn inv_mod_u32(num: u32) -> u64 {
     }
 
     assert!(a == 1);
-    // Ensure returned value is in-range
-    if u < 0 {
-        (u + mu) as u64
-    } else {
-        u as u64
-    }
+    // Downcasting acts like a mod 2^32 too.
+    u as u32
 }
 
 impl<'a> MontyReducer<'a> {
@@ -77,7 +72,7 @@ fn monty_redc(a: BigUint, mr: &MontyReducer) -> BigUint {
     // equivalent to masking a to 32 bits.
     let beta_mask = u32::max_value() as u64;
     // mu <- -N^(-1) mod β
-    let mu = (beta_mask-mr.n0inv)+1;
+    let mu = (beta_mask-mr.n0inv as u64)+1;
 
     // 1: for i = 0 to (n-1)
     for i in 0..n_size {

From 4d358154267ef73dff3ddf159d2a72e3e967ec59 Mon Sep 17 00:00:00 2001
From: Josh Stone <cuviper@gmail.com>
Date: Sun, 22 Oct 2017 15:28:59 -0700
Subject: [PATCH 06/14] bigint::monty: simplify work space allocation

---
 bigint/src/monty.rs | 5 +----
 1 file changed, 1 insertion(+), 4 deletions(-)

diff --git a/bigint/src/monty.rs b/bigint/src/monty.rs
index 8960c70..da352ed 100644
--- a/bigint/src/monty.rs
+++ b/bigint/src/monty.rs
@@ -1,4 +1,3 @@
-use std::iter::repeat;
 use integer::Integer;
 use traits::{Zero, One};
 
@@ -62,11 +61,9 @@ fn monty_redc(a: BigUint, mr: &MontyReducer) -> BigUint {
     let mut c = a.data;
     let n = &mr.n;
     let n_size = n.len();
-    let old_size = c.len();
 
     // Allocate sufficient work space
-    c.reserve(2*n_size+2-old_size);
-    c.extend(repeat(0).take(2*n_size+2-old_size));
+    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.

From 5a0de140c961f2f96b830f0600e5a3d69a1d0a12 Mon Sep 17 00:00:00 2001
From: Josh Stone <cuviper@gmail.com>
Date: Sun, 22 Oct 2017 15:30:17 -0700
Subject: [PATCH 07/14] bigint::monty: use mac_digit

---
 bigint/src/algorithms.rs |  2 +-
 bigint/src/monty.rs      | 18 +-----------------
 2 files changed, 2 insertions(+), 18 deletions(-)

diff --git a/bigint/src/algorithms.rs b/bigint/src/algorithms.rs
index de99ccc..5b74539 100644
--- a/bigint/src/algorithms.rs
+++ b/bigint/src/algorithms.rs
@@ -220,7 +220,7 @@ pub fn sub_sign(a: &[BigDigit], b: &[BigDigit]) -> (Sign, BigUint) {
 
 /// Three argument multiply accumulate:
 /// acc += b * c
-fn mac_digit(acc: &mut [BigDigit], b: &[BigDigit], c: BigDigit) {
+pub fn mac_digit(acc: &mut [BigDigit], b: &[BigDigit], c: BigDigit) {
     if c == 0 {
         return;
     }
diff --git a/bigint/src/monty.rs b/bigint/src/monty.rs
index da352ed..d522df5 100644
--- a/bigint/src/monty.rs
+++ b/bigint/src/monty.rs
@@ -73,27 +73,11 @@ fn monty_redc(a: BigUint, mr: &MontyReducer) -> BigUint {
 
     // 1: for i = 0 to (n-1)
     for i in 0..n_size {
-        // Carry storage
-        let mut carry = 0;
-
         // 2: q_i <- mu*c_i mod β
         let q_i = ((c[i] as u64) * mu) & beta_mask;
 
         // 3: C <- C + q_i * N * β^i
-        // When iterating over each word, this becomes:
-        for j in 0..n_size {
-            // c_(i+j) <- c_(i+j) + q_i * n_j
-            let x = (c[i+j] as u64) + q_i * (n[j] as u64) + carry;
-            c[i+j] = (x & beta_mask) as u32;
-            carry = x >> 32;
-        }
-
-        // Apply the remaining carry to the rest of the work space
-        for j in n_size..2*n_size-i+2 {
-            let x = (c[i+j] as u64) + carry;
-            c[i+j] = (x & beta_mask) as u32;
-            carry = x >> 32;
-        }
+        super::algorithms::mac_digit(&mut c[i..], n, q_i as u32);
     }
 
     // 4: R <- C * β^(-n)

From 5708db0f679636666638305987478d9769eec7e5 Mon Sep 17 00:00:00 2001
From: Josh Stone <cuviper@gmail.com>
Date: Sun, 22 Oct 2017 15:34:46 -0700
Subject: [PATCH 08/14] bigint::monty: simplify redc masks

---
 bigint/src/monty.rs | 7 +++----
 1 file changed, 3 insertions(+), 4 deletions(-)

diff --git a/bigint/src/monty.rs b/bigint/src/monty.rs
index d522df5..59d160f 100644
--- a/bigint/src/monty.rs
+++ b/bigint/src/monty.rs
@@ -67,17 +67,16 @@ fn monty_redc(a: BigUint, mr: &MontyReducer) -> BigUint {
 
     // β is the size of a word, in this case 32 bits. So "a mod β" is
     // equivalent to masking a to 32 bits.
-    let beta_mask = u32::max_value() as u64;
     // mu <- -N^(-1) mod β
-    let mu = (beta_mask-mr.n0inv as u64)+1;
+    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] as u64) * mu) & beta_mask;
+        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 as u32);
+        super::algorithms::mac_digit(&mut c[i..], n, q_i);
     }
 
     // 4: R <- C * β^(-n)

From 7fa27b60073bc5e0fa237ecbf35a27cec93e58fc Mon Sep 17 00:00:00 2001
From: Josh Stone <cuviper@gmail.com>
Date: Sun, 22 Oct 2017 15:36:33 -0700
Subject: [PATCH 09/14] bigint::monty: simplify redc return value

---
 bigint/src/monty.rs | 5 ++---
 1 file changed, 2 insertions(+), 3 deletions(-)

diff --git a/bigint/src/monty.rs b/bigint/src/monty.rs
index 59d160f..2c427b4 100644
--- a/bigint/src/monty.rs
+++ b/bigint/src/monty.rs
@@ -81,14 +81,13 @@ fn monty_redc(a: BigUint, mr: &MontyReducer) -> BigUint {
 
     // 4: R <- C * β^(-n)
     // This is an n-word bitshift, equivalent to skipping n words.
-    let r : Vec<u32> = c.iter().skip(n_size).cloned().collect();
-    let ret = BigUint::new(r);
+    let ret = BigUint::new(c[n_size..].to_vec());
 
     // 5: if R >= β^n then return R-N else return R.
     if &ret < mr.p {
         ret
     } else {
-        &ret-mr.p
+        ret - mr.p
     }
 }
 

From 96c4a26624f0ef733d4dc452371b12634dec9398 Mon Sep 17 00:00:00 2001
From: Josh Stone <cuviper@gmail.com>
Date: Sun, 22 Oct 2017 15:37:48 -0700
Subject: [PATCH 10/14] bigint::monty: simplify modpow parameter init

---
 bigint/src/monty.rs | 7 +++----
 1 file changed, 3 insertions(+), 4 deletions(-)

diff --git a/bigint/src/monty.rs b/bigint/src/monty.rs
index 2c427b4..5d5ec9c 100644
--- a/bigint/src/monty.rs
+++ b/bigint/src/monty.rs
@@ -106,10 +106,9 @@ pub fn monty_modpow(a: &BigUint, exp: &BigUint, modulus: &BigUint) -> BigUint{
     let mr = MontyReducer::new(modulus);
 
     // Calculate the Montgomery parameter
-    let mut r : BigUint = One::one();
-    while &r < mr.p {
-        r = r << 32;
-    }
+    let mut v = vec![0; mr.p.data.len()];
+    v.push(1);
+    let r = BigUint::new(v);
 
     // Map the base to the Montgomery domain
     let mut apri = a * &r % mr.p;

From b380880ed3224d6d6995dfbffc9813122f6a215a Mon Sep 17 00:00:00 2001
From: Josh Stone <cuviper@gmail.com>
Date: Sun, 22 Oct 2017 15:38:50 -0700
Subject: [PATCH 11/14] bigint::monty: simplify modpow zero test

---
 bigint/src/monty.rs | 5 ++---
 1 file changed, 2 insertions(+), 3 deletions(-)

diff --git a/bigint/src/monty.rs b/bigint/src/monty.rs
index 5d5ec9c..dc2ce96 100644
--- a/bigint/src/monty.rs
+++ b/bigint/src/monty.rs
@@ -1,5 +1,5 @@
 use integer::Integer;
-use traits::{Zero, One};
+use traits::Zero;
 
 use biguint::BigUint;
 
@@ -116,8 +116,7 @@ pub fn monty_modpow(a: &BigUint, exp: &BigUint, modulus: &BigUint) -> BigUint{
     // Binary exponentiation
     let mut ans = &r % mr.p;
     let mut e = exp.clone();
-    let zero = Zero::zero();
-    while e > zero {
+    while !e.is_zero() {
         if e.is_odd() {
             ans = monty_mult(ans, &apri, &mr);
         }

From bb0c9324b2a337625baa1eb1755df989734fdaeb Mon Sep 17 00:00:00 2001
From: Josh Stone <cuviper@gmail.com>
Date: Sun, 22 Oct 2017 15:45:01 -0700
Subject: [PATCH 12/14] bigint::monty: deduplicate mr.n and mr.p

---
 bigint/src/monty.rs | 22 ++++++++++------------
 1 file changed, 10 insertions(+), 12 deletions(-)

diff --git a/bigint/src/monty.rs b/bigint/src/monty.rs
index dc2ce96..3cf526b 100644
--- a/bigint/src/monty.rs
+++ b/bigint/src/monty.rs
@@ -4,8 +4,7 @@ use traits::Zero;
 use biguint::BigUint;
 
 struct MontyReducer<'a> {
-    p: &'a BigUint,
-    n: Vec<u32>,
+    n: &'a BigUint,
     n0inv: u32
 }
 
@@ -46,10 +45,9 @@ fn inv_mod_u32(num: u32) -> u32 {
 }
 
 impl<'a> MontyReducer<'a> {
-    fn new(p: &'a BigUint) -> Self {
-        let n : Vec<u32> = p.data.clone();
-        let n0inv = inv_mod_u32(n[0]);
-        MontyReducer { p: p, n: n, n0inv: n0inv }
+    fn new(n: &'a BigUint) -> Self {
+        let n0inv = inv_mod_u32(n.data[0]);
+        MontyReducer { n: n, n0inv: n0inv }
     }
 }
 
@@ -59,7 +57,7 @@ impl<'a> MontyReducer<'a> {
 // 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;
+    let n = &mr.n.data;
     let n_size = n.len();
 
     // Allocate sufficient work space
@@ -84,10 +82,10 @@ fn monty_redc(a: BigUint, mr: &MontyReducer) -> BigUint {
     let ret = BigUint::new(c[n_size..].to_vec());
 
     // 5: if R >= β^n then return R-N else return R.
-    if &ret < mr.p {
+    if &ret < mr.n {
         ret
     } else {
-        ret - mr.p
+        ret - mr.n
     }
 }
 
@@ -106,15 +104,15 @@ 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; mr.p.data.len()];
+    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 % mr.p;
+    let mut apri = a * &r % modulus;
 
     // Binary exponentiation
-    let mut ans = &r % mr.p;
+    let mut ans = &r % modulus;
     let mut e = exp.clone();
     while !e.is_zero() {
         if e.is_odd() {

From 35b7187e8344f4996711862462d38b53aca549e8 Mon Sep 17 00:00:00 2001
From: Josh Stone <cuviper@gmail.com>
Date: Sun, 22 Oct 2017 15:49:51 -0700
Subject: [PATCH 13/14] bigint::monty: use infallible conversions in tests

---
 bigint/src/tests/biguint.rs | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/bigint/src/tests/biguint.rs b/bigint/src/tests/biguint.rs
index 89101a6..f66564d 100644
--- a/bigint/src/tests/biguint.rs
+++ b/bigint/src/tests/biguint.rs
@@ -1092,10 +1092,10 @@ fn test_is_even() {
 #[test]
 fn test_modpow() {
     fn check(b: usize, e: usize, m: usize, r: usize) {
-        let big_b: BigUint = FromPrimitive::from_usize(b).unwrap();
-        let big_e: BigUint = FromPrimitive::from_usize(e).unwrap();
-        let big_m: BigUint = FromPrimitive::from_usize(m).unwrap();
-        let big_r: BigUint = FromPrimitive::from_usize(r).unwrap();
+        let big_b = BigUint::from(b);
+        let big_e = BigUint::from(e);
+        let big_m = BigUint::from(m);
+        let big_r = BigUint::from(r);
 
         assert_eq!(big_b.modpow(&big_e, &big_m), big_r);
     }

From ed10d617b5ce5a25c2ccb553ed2ab1278d02530d Mon Sep 17 00:00:00 2001
From: Josh Stone <cuviper@gmail.com>
Date: Sun, 22 Oct 2017 16:44:05 -0700
Subject: [PATCH 14/14] bigint: Add a modpow fallback for even modulus

---
 benches/bigint.rs           | 42 ++++++++++++++++++++++++-------------
 bigint/src/biguint.rs       | 29 ++++++++++++++++++++++++-
 bigint/src/tests/biguint.rs | 10 +++++++++
 3 files changed, 66 insertions(+), 15 deletions(-)

diff --git a/benches/bigint.rs b/benches/bigint.rs
index 9e18282..23932e6 100644
--- a/benches/bigint.rs
+++ b/benches/bigint.rs
@@ -256,25 +256,39 @@ fn pow_bench(b: &mut Bencher) {
     });
 }
 
+
+/// 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);
-    // This modulus is the prime from the 2048-bit MODP DH group:
-    // https://tools.ietf.org/html/rfc3526#section-3
-    let m = BigUint::from_str_radix("\
-        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", 16).unwrap();
+    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));
 }
diff --git a/bigint/src/biguint.rs b/bigint/src/biguint.rs
index 0bfc158..3e78a43 100644
--- a/bigint/src/biguint.rs
+++ b/bigint/src/biguint.rs
@@ -1625,7 +1625,34 @@ impl BigUint {
 
     /// Returns `(self ^ exponent) % modulus`.
     pub fn modpow(&self, exponent: &Self, modulus: &Self) -> Self {
-        monty_modpow(self, exponent, modulus)
+        assert!(!modulus.is_zero(), "divide by zero!");
+
+        // 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 }
+
+        let mut acc = base.clone();
+        while exp > one {
+            exp >>= 1;
+            base = &base * &base % modulus;
+            if exp.is_odd() {
+                acc = acc * &base % modulus;
+            }
+        }
+        acc
     }
 }
 
diff --git a/bigint/src/tests/biguint.rs b/bigint/src/tests/biguint.rs
index f66564d..e27cc17 100644
--- a/bigint/src/tests/biguint.rs
+++ b/bigint/src/tests/biguint.rs
@@ -1098,6 +1098,11 @@ fn test_modpow() {
         let big_r = BigUint::from(r);
 
         assert_eq!(big_b.modpow(&big_e, &big_m), big_r);
+
+        let even_m = &big_m << 1;
+        let even_modpow = big_b.modpow(&big_e, &even_m);
+        assert!(even_modpow < even_m);
+        assert_eq!(even_modpow % big_m, big_r);
     }
 
     check(1, 0, 11, 1);
@@ -1160,6 +1165,11 @@ fn test_modpow_big() {
         109c4735_6e7db425_7b5d74c7_0b709508", 16).unwrap();
 
     assert_eq!(b.modpow(&e, &m), r);
+
+    let even_m = &m << 1;
+    let even_modpow = b.modpow(&e, &even_m);
+    assert!(even_modpow < even_m);
+    assert_eq!(even_modpow % m, r);
 }
 
 fn to_str_pairs() -> Vec<(BigUint, Vec<(u32, String)>)> {