Add BoxedUint::pow

Initial support for modular exponentiation, adapted from the original
implementation of `pow_montgomery_form` this crate used prior to #248:

https://github.com/RustCrypto/crypto-bigint/blob/4838fd96e1bde8b0c5e0ce691c366c7ec930e466/src/uint/modular/pow.rs

Proptested against `num_bigint::BitUint::modpow`.

src/modular/boxed_residue.rs

@@ -2,6 +2,7 @@
 //! is chosen at runtime.
 
 mod mul;
+mod pow;
 
 use super::reduction::montgomery_reduction_boxed;
 use crate::{BoxedUint, Limb, NonZero, Word};

src/modular/boxed_residue/mul.rs

@@ -23,8 +23,14 @@ impl BoxedResidue {
 
     /// Computes the (reduced) square of a residue.
     pub fn square(&self) -> Self {
-        // TODO(tarcieri): optimized implementation
-        self.mul(self)
+        Self {
+            montgomery_form: square_montgomery_form(
+                &self.montgomery_form,
+                &self.residue_params.modulus,
+                self.residue_params.mod_neg_inv,
+            ),
+            residue_params: self.residue_params.clone(),
+        }
     }
 }
 
@@ -83,7 +89,7 @@ impl Square for BoxedResidue {
     }
 }
 
-fn mul_montgomery_form(
+pub(super) fn mul_montgomery_form(
     a: &BoxedUint,
     b: &BoxedUint,
     modulus: &BoxedUint,
@@ -100,3 +106,13 @@ fn mul_montgomery_form(
 
     ret
 }
+
+#[inline]
+pub(super) fn square_montgomery_form(
+    a: &BoxedUint,
+    modulus: &BoxedUint,
+    mod_neg_inv: Limb,
+) -> BoxedUint {
+    // TODO(tarcieri): optimized implementation
+    mul_montgomery_form(a, a, modulus, mod_neg_inv)
+}

src/modular/boxed_residue/pow.rs

@@ -0,0 +1,115 @@
+//! Modular exponentiation support.
+
+use super::{
+    mul::{mul_montgomery_form, square_montgomery_form},
+    BoxedResidue,
+};
+use crate::{BoxedUint, Limb, PowBoundedExp, Word};
+use subtle::ConstantTimeEq;
+
+impl BoxedResidue {
+    /// Raises to the `exponent` power.
+    pub fn pow(&self, exponent: &BoxedUint) -> Self {
+        self.pow_bounded_exp(exponent, exponent.bits_precision())
+    }
+
+    /// Raises to the `exponent` power,
+    /// with `exponent_bits` representing the number of (least significant) bits
+    /// to take into account for the exponent.
+    ///
+    /// NOTE: `exponent_bits` may be leaked in the time pattern.
+    pub fn pow_bounded_exp(&self, exponent: &BoxedUint, exponent_bits: usize) -> Self {
+        Self {
+            montgomery_form: pow_montgomery_form(
+                &self.montgomery_form,
+                exponent,
+                exponent_bits,
+                &self.residue_params.modulus,
+                &self.residue_params.r,
+                self.residue_params.mod_neg_inv,
+            ),
+            residue_params: self.residue_params.clone(),
+        }
+    }
+}
+
+impl PowBoundedExp<BoxedUint> for BoxedResidue {
+    fn pow_bounded_exp(&self, exponent: &BoxedUint, exponent_bits: usize) -> Self {
+        self.pow_bounded_exp(exponent, exponent_bits)
+    }
+}
+
+/// Performs modular exponentiation using Montgomery's ladder.
+/// `exponent_bits` represents the number of bits to take into account for the exponent.
+///
+/// NOTE: this value is leaked in the time pattern.
+fn pow_montgomery_form(
+    x: &BoxedUint,
+    exponent: &BoxedUint,
+    exponent_bits: usize,
+    modulus: &BoxedUint,
+    r: &BoxedUint,
+    mod_neg_inv: Limb,
+) -> BoxedUint {
+    if exponent_bits == 0 {
+        return r.clone(); // 1 in Montgomery form
+    }
+
+    const WINDOW: usize = 4;
+    const WINDOW_MASK: Word = (1 << WINDOW) - 1;
+
+    // powers[i] contains x^i
+    let mut powers = vec![r.clone(); 1 << WINDOW];
+    powers[1] = x.clone();
+    let mut i = 2;
+    while i < powers.len() {
+        powers[i] = mul_montgomery_form(&powers[i - 1], x, modulus, mod_neg_inv);
+        i += 1;
+    }
+
+    let starting_limb = (exponent_bits - 1) / Limb::BITS;
+    let starting_bit_in_limb = (exponent_bits - 1) % Limb::BITS;
+    let starting_window = starting_bit_in_limb / WINDOW;
+    let starting_window_mask = (1 << (starting_bit_in_limb % WINDOW + 1)) - 1;
+
+    let mut z = r.clone(); // 1 in Montgomery form
+
+    let mut limb_num = starting_limb + 1;
+    while limb_num > 0 {
+        limb_num -= 1;
+        let w = exponent.as_limbs()[limb_num].0;
+
+        let mut window_num = if limb_num == starting_limb {
+            starting_window + 1
+        } else {
+            Limb::BITS / WINDOW
+        };
+        while window_num > 0 {
+            window_num -= 1;
+
+            let mut idx = (w >> (window_num * WINDOW)) & WINDOW_MASK;
+
+            if limb_num == starting_limb && window_num == starting_window {
+                idx &= starting_window_mask;
+            } else {
+                let mut i = 0;
+                while i < WINDOW {
+                    i += 1;
+                    z = square_montgomery_form(&z, modulus, mod_neg_inv);
+                }
+            }
+
+            // Constant-time lookup in the array of powers
+            let mut power = powers[0].clone();
+            let mut i = 1;
+            while i < 1 << WINDOW {
+                power = BoxedUint::conditional_select(&power, &powers[i], (i as Word).ct_eq(&idx));
+                i += 1;
+            }
+
+            z = mul_montgomery_form(&z, &power, modulus, mod_neg_inv);
+        }
+    }
+
+    z
+}

tests/boxed_residue_proptests.rs

@@ -18,6 +18,20 @@ fn retrieve_biguint(residue: &BoxedResidue) -> BigUint {
     to_biguint(&residue.retrieve())
 }
 
+fn reduce(n: &BoxedUint, p: BoxedResidueParams) -> BoxedResidue {
+    let bits_precision = p.modulus().bits_precision();
+    let modulus = NonZero::new(p.modulus().clone()).unwrap();
+
+    let n = match n.bits_precision().cmp(&bits_precision) {
+        Ordering::Less => n.widen(bits_precision),
+        Ordering::Equal => n.clone(),
+        Ordering::Greater => n.shorten(bits_precision),
+    };
+
+    let n_reduced = n.rem_vartime(&modulus).widen(p.bits_precision());
+    BoxedResidue::new(&n_reduced, p)
+}
+
 prop_compose! {
     /// Generate a random `BoxedUint`.
     fn uint()(mut bytes in any::<Vec<u8>>()) -> BoxedUint {
@@ -29,47 +43,45 @@ prop_compose! {
 }
 prop_compose! {
     /// Generate a random modulus.
-    fn modulus()(mut a in uint()) -> BoxedResidueParams {
-        if a.is_even().into() {
-            a = a.wrapping_add(&BoxedUint::one());
+    fn modulus()(mut n in uint()) -> BoxedResidueParams {
+        if n.is_even().into() {
+            n = n.wrapping_add(&BoxedUint::one());
         }
 
-        BoxedResidueParams::new(a).expect("modulus should be valid")
+        BoxedResidueParams::new(n).expect("modulus should be valid")
     }
 }
 prop_compose! {
     /// Generate two residues with a common modulus.
-    fn residue_pair()(a in uint(), b in uint(), p in modulus()) -> (BoxedResidue, BoxedResidue) {
-        fn reduce(n: &BoxedUint, p: BoxedResidueParams) -> BoxedResidue {
-            let bits_precision = p.modulus().bits_precision();
-            let modulus = NonZero::new(p.modulus().clone()).unwrap();
-
-            let n = match n.bits_precision().cmp(&bits_precision) {
-                Ordering::Less => n.widen(bits_precision),
-                Ordering::Equal => n.clone(),
-                Ordering::Greater => n.shorten(bits_precision)
-            };
-
-            let n_reduced = n.rem_vartime(&modulus).widen(p.bits_precision());
-            BoxedResidue::new(&n_reduced, p)
-        }
-
-
-        (reduce(&a, p.clone()), reduce(&b, p.clone()))
+    fn residue_pair()(a in uint(), b in uint(), n in modulus()) -> (BoxedResidue, BoxedResidue) {
+        (reduce(&a, n.clone()), reduce(&b, n.clone()))
     }
 }
 
 proptest! {
     #[test]
     fn mul((a, b) in residue_pair()) {
         let p = a.params().modulus();
-        let c = &a * &b;
+        let actual = &a * &b;
 
         let a_bi = retrieve_biguint(&a);
         let b_bi = retrieve_biguint(&b);
         let p_bi = to_biguint(&p);
-        let c_bi = (a_bi * b_bi) % p_bi;
+        let expected = (a_bi * b_bi) % p_bi;
+
+        prop_assert_eq!(retrieve_biguint(&actual), expected);
+    }
+
+    #[test]
+    fn pow(a in uint(), b in uint(), n in modulus()) {
+        let a = reduce(&a, n.clone());
+        let actual = a.pow(&b);
+
+        let a_bi = retrieve_biguint(&a);
+        let b_bi = to_biguint(&b);
+        let n_bi = to_biguint(n.modulus());
+        let expected = a_bi.modpow(&b_bi, &n_bi);
 
-        prop_assert_eq!(retrieve_biguint(&c), c_bi);
+        prop_assert_eq!(retrieve_biguint(&actual), expected);
     }
 }