Multi-exponentiation

benches/bench.rs

@@ -3,7 +3,7 @@ use criterion::{
 };
 use crypto_bigint::{
     modular::runtime_mod::{DynResidue, DynResidueParams},
-    Limb, NonZero, Random, Reciprocal, U128, U2048, U256,
+    Limb, MultiExponentiate, NonZero, Random, Reciprocal, U128, U2048, U256,
 };
 use rand_core::OsRng;
 
@@ -101,6 +101,34 @@ 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.as_slice(),
+                        )
+                    },
+                    BatchSize::SmallInput,
+                )
+            },
+        );
+    }
 }
 
 fn bench_montgomery_conversion<M: Measurement>(group: &mut BenchmarkGroup<'_, M>) {

src/lib.rs

@@ -152,6 +152,8 @@
 //! [`Sub`]: core::ops::Sub
 
 #[cfg(feature = "alloc")]
+#[allow(unused_imports)]
+#[macro_use]
 extern crate alloc;
 
 #[macro_use]

src/traits.rs

@@ -292,6 +292,53 @@ pub trait PowBoundedExp<Exponent> {
     fn pow_bounded_exp(&self, exponent: &Exponent, exponent_bits: usize) -> Self;
 }
 
+/// Linear combination (extended version).
+///
+/// This trait enables providing an optimized implementation of
+/// linear combinations (e.g. Shamir's Trick).
+// TODO(tarcieri): replace the current `LinearCombination` with this in the next release
+
+/// Performs modular multi-exponentiation using Montgomery's ladder.
+///
+/// See: Straus, E. G. Problems and solutions: Addition chains of vectors. American Mathematical Monthly 71 (1964), 806–808.
+pub trait MultiExponentiate<Exponent, BasesAndExponents>: Pow<Exponent> + Sized
+where
+    BasesAndExponents: AsRef<[(Self, Exponent)]> + ?Sized,
+{
+    /// Calculates `x1 ^ k1 * ... * xn ^ kn`.
+    fn multi_exponentiate(bases_and_exponents: &BasesAndExponents) -> Self;
+}
+
+impl<
+        T: MultiExponentiateBoundedExp<Exponent, BasesAndExponents>,
+        Exponent: Bounded,
+        BasesAndExponents,
+    > MultiExponentiate<Exponent, BasesAndExponents> for T
+where
+    BasesAndExponents: AsRef<[(Self, Exponent)]> + ?Sized,
+{
+    fn multi_exponentiate(bases_and_exponents: &BasesAndExponents) -> Self {
+        Self::multi_exponentiate_bounded_exp(bases_and_exponents, Exponent::BITS)
+    }
+}
+
+/// 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.
+pub trait MultiExponentiateBoundedExp<Exponent, BasesAndExponents>: Pow<Exponent> + Sized
+where
+    BasesAndExponents: AsRef<[(Self, Exponent)]> + ?Sized,
+{
+    /// Calculates `x1 ^ k1 * ... * xn ^ kn`.
+    fn multi_exponentiate_bounded_exp(
+        bases_and_exponents: &BasesAndExponents,
+        exponent_bits: usize,
+    ) -> Self;
+}
+
 /// Constant-time inversion.
 pub trait Invert: Sized {
     /// Output of the inversion.

src/uint/modular/constant_mod/const_pow.rs

@@ -1,6 +1,11 @@
-use crate::{modular::pow::pow_montgomery_form, PowBoundedExp, Uint};
+use crate::{modular::pow::pow_montgomery_form, MultiExponentiateBoundedExp, PowBoundedExp, Uint};
 
 use super::{Residue, ResidueParams};
+use crate::modular::pow::multi_exponentiate_montgomery_form_array;
+#[cfg(feature = "alloc")]
+use crate::modular::pow::multi_exponentiate_montgomery_form_slice;
+#[cfg(feature = "alloc")]
+use alloc::vec::Vec;
 
 impl<MOD: ResidueParams<LIMBS>, const LIMBS: usize> Residue<MOD, LIMBS> {
     /// Raises to the `exponent` power.
@@ -35,16 +40,74 @@ impl<MOD: ResidueParams<LIMBS>, const LIMBS: usize> Residue<MOD, LIMBS> {
     }
 }
 
-impl<MOD: ResidueParams<LIMBS>, const LIMBS: usize> PowBoundedExp<Uint<LIMBS>>
-    for Residue<MOD, LIMBS>
+impl<MOD: ResidueParams<LIMBS>, const LIMBS: usize, const RHS_LIMBS: usize>
+    PowBoundedExp<Uint<RHS_LIMBS>> for Residue<MOD, LIMBS>
 {
-    fn pow_bounded_exp(&self, exponent: &Uint<LIMBS>, exponent_bits: usize) -> Self {
+    fn pow_bounded_exp(&self, exponent: &Uint<RHS_LIMBS>, exponent_bits: usize) -> Self {
         self.pow_bounded_exp(exponent, exponent_bits)
     }
 }
 
+impl<const N: usize, MOD: ResidueParams<LIMBS>, const LIMBS: usize, const RHS_LIMBS: usize>
+    MultiExponentiateBoundedExp<Uint<RHS_LIMBS>, [(Self, Uint<RHS_LIMBS>); N]>
+    for Residue<MOD, LIMBS>
+{
+    fn multi_exponentiate_bounded_exp(
+        bases_and_exponents: &[(Self, Uint<RHS_LIMBS>); N],
+        exponent_bits: usize,
+    ) -> Self {
+        let mut bases_and_exponents_montgomery_form =
+            [(Uint::<LIMBS>::ZERO, Uint::<RHS_LIMBS>::ZERO); N];
+
+        let mut i = 0;
+        while i < N {
+            let (base, exponent) = bases_and_exponents[i];
+            bases_and_exponents_montgomery_form[i] = (base.montgomery_form, exponent);
+            i += 1;
+        }
+
+        Self {
+            montgomery_form: multi_exponentiate_montgomery_form_array(
+                &bases_and_exponents_montgomery_form,
+                exponent_bits,
+                &MOD::MODULUS,
+                &MOD::R,
+                MOD::MOD_NEG_INV,
+            ),
+            phantom: core::marker::PhantomData,
+        }
+    }
+}
+
+#[cfg(feature = "alloc")]
+impl<MOD: ResidueParams<LIMBS>, const LIMBS: usize, const RHS_LIMBS: usize>
+    MultiExponentiateBoundedExp<Uint<RHS_LIMBS>, [(Self, Uint<RHS_LIMBS>)]>
+    for Residue<MOD, LIMBS>
+{
+    fn multi_exponentiate_bounded_exp(
+        bases_and_exponents: &[(Self, Uint<RHS_LIMBS>)],
+        exponent_bits: usize,
+    ) -> Self {
+        let bases_and_exponents: Vec<(Uint<LIMBS>, Uint<RHS_LIMBS>)> = bases_and_exponents
+            .iter()
+            .map(|(base, exp)| (base.montgomery_form, *exp))
+            .collect();
+        Self {
+            montgomery_form: multi_exponentiate_montgomery_form_slice(
+                &bases_and_exponents,
+                exponent_bits,
+                &MOD::MODULUS,
+                &MOD::R,
+                MOD::MOD_NEG_INV,
+            ),
+            phantom: core::marker::PhantomData,
+        }
+    }
+}
+
 #[cfg(test)]
 mod tests {
+    use crate::traits::MultiExponentiate;
     use crate::{const_residue, impl_modulus, modular::constant_mod::ResidueParams, U256};
 
     impl_modulus!(
@@ -98,4 +161,73 @@ mod tests {
             U256::from_be_hex("3681BC0FEA2E5D394EB178155A127B0FD2EF405486D354251C385BDD51B9D421");
         assert_eq!(res.retrieve(), expected);
     }
+
+    #[test]
+    fn test_multi_exp_array() {
+        let base = U256::from(2u8);
+        let base_mod = const_residue!(base, Modulus);
+
+        let exponent = U256::from(33u8);
+        let bases_and_exponents = [(base_mod, exponent)];
+        let res =
+            crate::modular::constant_mod::Residue::<Modulus, { U256::LIMBS }>::multi_exponentiate(
+                &bases_and_exponents,
+            );
+
+        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 bases_and_exponents = [(base_mod, exponent), (base2_mod, exponent2)];
+        let res =
+            crate::modular::constant_mod::Residue::<Modulus, { U256::LIMBS }>::multi_exponentiate(
+                &bases_and_exponents,
+            );
+
+        assert_eq!(res, expected);
+    }
+
+    #[cfg(feature = "alloc")]
+    #[test]
+    fn test_multi_exp_slice() {
+        let base = U256::from(2u8);
+        let base_mod = const_residue!(base, Modulus);
+
+        let exponent = U256::from(33u8);
+        let bases_and_exponents = vec![(base_mod, exponent)];
+        let res =
+            crate::modular::constant_mod::Residue::<Modulus, { U256::LIMBS }>::multi_exponentiate(
+                bases_and_exponents.as_slice(),
+            );
+
+        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 bases_and_exponents = vec![(base_mod, exponent), (base2_mod, exponent2)];
+        let res =
+            crate::modular::constant_mod::Residue::<Modulus, { U256::LIMBS }>::multi_exponentiate(
+                bases_and_exponents.as_slice(),
+            );
+
+        assert_eq!(res, expected);
+    }
 }