diff --git a/benches/bigint.rs b/benches/bigint.rs index 0eafb93..23932e6 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,40 @@ 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); + 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/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/biguint.rs b/bigint/src/biguint.rs index 28d1e90..439332f 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::monty_modpow; use UsizePromotion; @@ -233,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() { @@ -240,6 +250,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 { @@ -1611,6 +1622,38 @@ impl BigUint { self.normalize(); self } + + /// Returns `(self ^ exponent) % modulus`. + pub fn modpow(&self, exponent: &Self, modulus: &Self) -> Self { + 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 + } } #[cfg(feature = "serde")] diff --git a/bigint/src/monty.rs b/bigint/src/monty.rs new file mode 100644 index 0000000..3cf526b --- /dev/null +++ b/bigint/src/monty.rs @@ -0,0 +1,127 @@ +use integer::Integer; +use traits::Zero; + +use biguint::BigUint; + +struct MontyReducer<'a> { + n: &'a BigUint, + 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) -> 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; + + // 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); + // Downcasting acts like a mod 2^32 too. + u as u32 +} + +impl<'a> MontyReducer<'a> { + fn new(n: &'a BigUint) -> Self { + let n0inv = inv_mod_u32(n.data[0]); + MontyReducer { n: n, n0inv: n0inv } + } +} + +// Montgomery Reduction +// +// Reference: +// Brent & Zimmermann, Modern Computer Arithmetic, v0.5.9, Algorithm 2.6 +fn monty_redc(a: BigUint, mr: &MontyReducer) -> BigUint { + let mut c = a.data; + let n = &mr.n.data; + let n_size = n.len(); + + // Allocate sufficient work space + c.resize(2 * n_size + 2, 0); + + // β is the size of a word, in this case 32 bits. So "a mod β" is + // equivalent to masking a to 32 bits. + // mu <- -N^(-1) mod β + let mu = 0u32.wrapping_sub(mr.n0inv); + + // 1: for i = 0 to (n-1) + for i in 0..n_size { + // 2: q_i <- mu*c_i mod β + let q_i = c[i].wrapping_mul(mu); + + // 3: C <- C + q_i * N * β^i + super::algorithms::mac_digit(&mut c[i..], n, q_i); + } + + // 4: R <- C * β^(-n) + // This is an n-word bitshift, equivalent to skipping n words. + let ret = BigUint::new(c[n_size..].to_vec()); + + // 5: if R >= β^n then return R-N else return R. + if &ret < mr.n { + ret + } else { + ret - mr.n + } +} + +// 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, modulus: &BigUint) -> BigUint{ + let mr = MontyReducer::new(modulus); + + // Calculate the Montgomery parameter + let mut v = vec![0; modulus.data.len()]; + v.push(1); + let r = BigUint::new(v); + + // Map the base to the Montgomery domain + let mut apri = a * &r % modulus; + + // Binary exponentiation + let mut ans = &r % modulus; + let mut e = exp.clone(); + while !e.is_zero() { + if e.is_odd() { + ans = monty_mult(ans, &apri, &mr); + } + apri = monty_sqr(apri, &mr); + e = e >> 1; + } + + // 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..e27cc17 100644 --- a/bigint/src/tests/biguint.rs +++ b/bigint/src/tests/biguint.rs @@ -1089,6 +1089,89 @@ 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::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); + + 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); + 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); + + 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)>)> { let bits = big_digit::BITS; vec![(Zero::zero(), @@ -1468,6 +1551,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(); @@ -1476,6 +1561,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]