Multi-exponentiation

benches/bench.rs

@@ -101,6 +101,35 @@ fn bench_montgomery_ops<M: Measurement>(group: &mut BenchmarkGroup<'_, M>) {
             BatchSize::SmallInput,
         )
     });
+
+    for i in [1, 2, 3, 4, 10, 100] {
+        group.bench_function(
+            format!("multi_exponentiate for {i} bases, U256^U256"),
+            |b| {
+                b.iter_batched(
+                    || {
+                        let bases_and_exponents: Vec<(DynResidue<{ U256::LIMBS }>, U256)> = (1..=i)
+                            .map(|_| {
+                                let x = U256::random(&mut OsRng);
+                                let x_m = DynResidue::new(&x, params);
+                                let p = U256::random(&mut OsRng) | (U256::ONE << (U256::BITS - 1));
+                                (x_m, p)
+                            })
+                            .collect();
+
+                        bases_and_exponents
+                    },
+                    |bases_and_exponents| {
+                        DynResidue::<{ U256::LIMBS }>::multi_exponentiate(
+                            bases_and_exponents,
+                            params,
+                        )
+                    },
+                    BatchSize::SmallInput,
+                )
+            },
+        );
+    }
 }
 
 fn bench_montgomery_conversion<M: Measurement>(group: &mut BenchmarkGroup<'_, M>) {

src/uint/modular/constant_mod/const_pow.rs

@@ -1,6 +1,11 @@
 use crate::{modular::pow::pow_montgomery_form, PowBoundedExp, Uint};
-
+#[cfg(feature = "alloc")]
+extern crate alloc;
 use super::{Residue, ResidueParams};
+#[cfg(feature = "alloc")]
+use crate::modular::pow::multi_exponentiate_montgomery_form;
+#[cfg(feature = "alloc")]
+use alloc::vec::Vec;
 
 impl<MOD: ResidueParams<LIMBS>, const LIMBS: usize> Residue<MOD, LIMBS> {
     /// Raises to the `exponent` power.
@@ -33,6 +38,40 @@ impl<MOD: ResidueParams<LIMBS>, const LIMBS: usize> Residue<MOD, LIMBS> {
             phantom: core::marker::PhantomData,
         }
     }
+
+    /// Raises each `base` to its `exponent` power,
+    #[cfg(feature = "alloc")]
+    pub fn multi_exponentiate(
+        bases_and_exponents: Vec<(Self, Uint<LIMBS>)>,
+    ) -> Residue<MOD, LIMBS> {
+        Self::multi_exponentiate_bounded_exp(bases_and_exponents, Uint::<LIMBS>::BITS)
+    }
+
+    /// Raises each `base` to its `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.
+    #[cfg(feature = "alloc")]
+    pub fn multi_exponentiate_bounded_exp(
+        bases_and_exponents: Vec<(Self, Uint<LIMBS>)>,
+        exponent_bits: usize,
+    ) -> Residue<MOD, LIMBS> {
+        let bases_and_exponents: Vec<(Uint<LIMBS>, Uint<LIMBS>)> = bases_and_exponents
+            .iter()
+            .map(|(base, exp)| (base.montgomery_form, *exp))
+            .collect();
+        Self {
+            montgomery_form: multi_exponentiate_montgomery_form(
+                &bases_and_exponents,
+                exponent_bits,
+                &MOD::MODULUS,
+                &MOD::R,
+                MOD::MOD_NEG_INV,
+            ),
+            phantom: core::marker::PhantomData,
+        }
+    }
 }
 
 impl<MOD: ResidueParams<LIMBS>, const LIMBS: usize> PowBoundedExp<Uint<LIMBS>>
@@ -98,4 +137,39 @@ mod tests {
             U256::from_be_hex("3681BC0FEA2E5D394EB178155A127B0FD2EF405486D354251C385BDD51B9D421");
         assert_eq!(res.retrieve(), expected);
     }
+
+    #[cfg(feature = "alloc")]
+    #[test]
+    fn test_multi_exp() {
+        let base = U256::from(2u8);
+        let base_mod = const_residue!(base, Modulus);
+
+        let exponent = U256::from(33u8);
+
+        let res =
+            crate::modular::constant_mod::Residue::<Modulus, { U256::LIMBS }>::multi_exponentiate(
+                alloc::vec![(base_mod, exponent)],
+            );
+
+        let expected =
+            U256::from_be_hex("0000000000000000000000000000000000000000000000000000000200000000");
+
+        assert_eq!(res.retrieve(), expected);
+
+        let base2 =
+            U256::from_be_hex("3435D18AA8313EBBE4D20002922225B53F75DC4453BB3EEC0378646F79B524A4");
+        let base2_mod = const_residue!(base2, Modulus);
+
+        let exponent2 =
+            U256::from_be_hex("77117F1273373C26C700D076B3F780074D03339F56DD0EFB60E7F58441FD3685");
+
+        let expected = base_mod.pow(&exponent) * base2_mod.pow(&exponent2);
+
+        let res =
+            crate::modular::constant_mod::Residue::<Modulus, { U256::LIMBS }>::multi_exponentiate(
+                alloc::vec![(base_mod, exponent), (base2_mod, exponent2),],
+            );
+
+        assert_eq!(res, expected);
+    }
 }

src/uint/modular/pow.rs

@@ -2,6 +2,9 @@ use crate::{Limb, Uint, Word};
 
 use super::mul::{mul_montgomery_form, square_montgomery_form};
 
+const WINDOW: usize = 4;
+const WINDOW_MASK: Word = (1 << WINDOW) - 1;
+
 /// Performs modular exponentiation using Montgomery's ladder.
 /// `exponent_bits` represents the number of bits to take into account for the exponent.
 ///
@@ -18,9 +21,6 @@ pub const fn pow_montgomery_form<const LIMBS: usize, const RHS_LIMBS: usize>(
         return *r; // 1 in Montgomery form
     }
 
-    const WINDOW: usize = 4;
-    const WINDOW_MASK: Word = (1 << WINDOW) - 1;
-
     // powers[i] contains x^i
     let mut powers = [*r; 1 << WINDOW];
     powers[1] = *x;
@@ -30,6 +30,65 @@ pub const fn pow_montgomery_form<const LIMBS: usize, const RHS_LIMBS: usize>(
         i += 1;
     }
 
+    multi_exponentiate_montgomery_form_internal(
+        &[(powers, *exponent)],
+        exponent_bits,
+        modulus,
+        r,
+        mod_neg_inv,
+    )
+}
+
+/// Performs modular multi-exponentiation using Montgomery's ladder.
+/// `exponent_bits` represents the number of bits to take into account for the exponent.
+///
+/// See: Straus, E. G. Problems and solutions: Addition chains of vectors. American Mathematical Monthly 71 (1964), 806–808.
+///
+/// NOTE: this value is leaked in the time pattern.
+#[cfg(feature = "alloc")]
+pub fn multi_exponentiate_montgomery_form<const LIMBS: usize, const RHS_LIMBS: usize>(
+    bases_and_exponents: &[(Uint<LIMBS>, Uint<RHS_LIMBS>)],
+    exponent_bits: usize,
+    modulus: &Uint<LIMBS>,
+    r: &Uint<LIMBS>,
+    mod_neg_inv: Limb,
+) -> Uint<LIMBS> {
+    if exponent_bits == 0 {
+        return *r; // 1 in Montgomery form
+    }
+
+    let powers_and_exponents: alloc::vec::Vec<([Uint<LIMBS>; 1 << WINDOW], Uint<RHS_LIMBS>)> =
+        bases_and_exponents
+            .iter()
+            .map(|(base, exponent)| {
+                // powers[i] contains x^i
+                let mut powers = [*r; 1 << WINDOW];
+                powers[1] = *base;
+                let mut i = 2;
+                while i < powers.len() {
+                    powers[i] = mul_montgomery_form(&powers[i - 1], base, modulus, mod_neg_inv);
+                    i += 1;
+                }
+                (powers, *exponent)
+            })
+            .collect();
+
+    multi_exponentiate_montgomery_form_internal(
+        powers_and_exponents.as_slice(),
+        exponent_bits,
+        modulus,
+        r,
+        mod_neg_inv,
+    )
+}
+
+const fn multi_exponentiate_montgomery_form_internal<const LIMBS: usize, const RHS_LIMBS: usize>(
+    powers_and_exponents: &[([Uint<LIMBS>; 1 << WINDOW], Uint<RHS_LIMBS>)],
+    exponent_bits: usize,
+    modulus: &Uint<LIMBS>,
+    r: &Uint<LIMBS>,
+    mod_neg_inv: Limb,
+) -> Uint<LIMBS> {
     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;
@@ -40,7 +99,6 @@ pub const fn pow_montgomery_form<const LIMBS: usize, const RHS_LIMBS: usize>(
     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
@@ -50,28 +108,36 @@ pub const fn pow_montgomery_form<const LIMBS: usize, const RHS_LIMBS: usize>(
         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 {
+            if limb_num != starting_limb || window_num != starting_window {
                 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];
-            let mut i = 1;
-            while i < 1 << WINDOW {
-                let choice = Limb::ct_eq(Limb(i as Word), Limb(idx));
-                power = Uint::<LIMBS>::ct_select(&power, &powers[i], choice);
+            let mut i = 0;
+            while i < powers_and_exponents.len() {
+                let (powers, exponent) = powers_and_exponents[i];
+                let w = exponent.as_limbs()[limb_num].0;
+                let mut idx = (w >> (window_num * WINDOW)) & WINDOW_MASK;
+
+                if limb_num == starting_limb && window_num == starting_window {
+                    idx &= starting_window_mask;
+                }
+
+                // Constant-time lookup in the array of powers
+                let mut power = powers[0];
+                let mut j = 1;
+                while j < 1 << WINDOW {
+                    let choice = Limb::ct_eq(Limb(j as Word), Limb(idx));
+                    power = Uint::<LIMBS>::ct_select(&power, &powers[j], choice);
+                    j += 1;
+                }
+
+                z = mul_montgomery_form(&z, &power, modulus, mod_neg_inv);
                 i += 1;
             }
-
-            z = mul_montgomery_form(&z, &power, modulus, mod_neg_inv);
         }
     }
 

src/uint/modular/runtime_mod/runtime_pow.rs

@@ -1,6 +1,7 @@
-use crate::{modular::pow::pow_montgomery_form, PowBoundedExp, Uint};
-
 use super::DynResidue;
+#[cfg(feature = "alloc")]
+use crate::modular::{pow::multi_exponentiate_montgomery_form, runtime_mod::DynResidueParams};
+use crate::{modular::pow::pow_montgomery_form, PowBoundedExp, Uint};
 
 impl<const LIMBS: usize> DynResidue<LIMBS> {
     /// Raises to the `exponent` power.
@@ -33,6 +34,46 @@ impl<const LIMBS: usize> DynResidue<LIMBS> {
             residue_params: self.residue_params,
         }
     }
+
+    /// Raises each `base` to its `exponent` power.
+    #[cfg(feature = "alloc")]
+    pub fn multi_exponentiate(
+        bases_and_exponents: alloc::vec::Vec<(Self, Uint<LIMBS>)>,
+        residue_params: DynResidueParams<LIMBS>,
+    ) -> DynResidue<LIMBS> {
+        Self::multi_exponentiate_bounded_exp(
+            bases_and_exponents,
+            Uint::<LIMBS>::BITS,
+            residue_params,
+        )
+    }
+
+    /// Raises each `base` to its `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.
+    #[cfg(feature = "alloc")]
+    pub fn multi_exponentiate_bounded_exp(
+        bases_and_exponents: alloc::vec::Vec<(Self, Uint<LIMBS>)>,
+        exponent_bits: usize,
+        residue_params: DynResidueParams<LIMBS>,
+    ) -> Self {
+        let bases_and_exponents: alloc::vec::Vec<(Uint<LIMBS>, Uint<LIMBS>)> = bases_and_exponents
+            .iter()
+            .map(|(base, exp)| (base.montgomery_form, *exp))
+            .collect();
+        Self {
+            montgomery_form: multi_exponentiate_montgomery_form(
+                &bases_and_exponents,
+                exponent_bits,
+                &residue_params.modulus,
+                &residue_params.r,
+                residue_params.mod_neg_inv,
+            ),
+            residue_params,
+        }
+    }
 }
 
 impl<const LIMBS: usize> PowBoundedExp<Uint<LIMBS>> for DynResidue<LIMBS> {