Add Montgomery multiexponentiation

4 files changed, 95 insertions, 0 deletions

diff --git a/src/lib/math/numbertheory/monty.cpp b/src/lib/math/numbertheory/monty.cpp
index 41c01abb1..33d15de5b 100644
--- a/src/lib/math/numbertheory/monty.cpp
+++ b/src/lib/math/numbertheory/monty.cpp
@@ -301,6 +301,12 @@ Montgomery_Int Montgomery_Int::operator*(const Montgomery_Int& other) const
    return Montgomery_Int(m_params, m_params->mul(m_v, other.m_v, ws), false);
    }
 
+Montgomery_Int Montgomery_Int::mul(const Montgomery_Int& other,
+                                   secure_vector<word>& ws) const
+   {
+   return Montgomery_Int(m_params, m_params->mul(m_v, other.m_v, ws), false);
+   }
+
 Montgomery_Int& Montgomery_Int::mul_by(const Montgomery_Int& other,
                                        secure_vector<word>& ws)
    {
diff --git a/src/lib/math/numbertheory/monty.h b/src/lib/math/numbertheory/monty.h
index 499a4ac91..9f369f1a5 100644
--- a/src/lib/math/numbertheory/monty.h
+++ b/src/lib/math/numbertheory/monty.h
@@ -75,6 +75,9 @@ class Montgomery_Int final
 
       Montgomery_Int& operator*=(const secure_vector<word>& other);
 
+      Montgomery_Int mul(const Montgomery_Int& other,
+                         secure_vector<word>& ws) const;
+
       Montgomery_Int& mul_by(const Montgomery_Int& other,
                              secure_vector<word>& ws);
 
diff --git a/src/lib/math/numbertheory/monty_exp.cpp b/src/lib/math/numbertheory/monty_exp.cpp
index 657686044..6f70ef7c6 100644
--- a/src/lib/math/numbertheory/monty_exp.cpp
+++ b/src/lib/math/numbertheory/monty_exp.cpp
@@ -8,6 +8,7 @@
 
 #include <botan/internal/monty_exp.h>
 #include <botan/internal/ct_utils.h>
+#include <botan/internal/rounding.h>
 #include <botan/numthry.h>
 #include <botan/reducer.h>
 #include <botan/monty.h>
@@ -125,5 +126,81 @@ BigInt monty_execute(const Montgomery_Exponentation_State& precomputed_state,
    return precomputed_state.exponentiation(k);
    }
 
+BigInt monty_multi_exp(std::shared_ptr<const Montgomery_Params> params_p,
+                       const BigInt& x_bn,
+                       const BigInt& z1,
+                       const BigInt& y_bn,
+                       const BigInt& z2)
+   {
+   if(z1.is_negative() || z2.is_negative())
+      throw Invalid_Argument("multi_exponentiate exponents must be positive");
+
+   const size_t z_bits = round_up(std::max(z1.bits(), z2.bits()), 2);
+
+   secure_vector<word> ws;
+
+   const Montgomery_Int one(params_p, params_p->R1(), false);
+   //const Montgomery_Int one(params_p, 1);
+
+   const Montgomery_Int x1(params_p, x_bn);
+   const Montgomery_Int x2 = x1.square(ws);
+   const Montgomery_Int x3 = x2.mul(x1, ws);
+
+   const Montgomery_Int y1(params_p, y_bn);
+   const Montgomery_Int y2 = y1.square(ws);
+   const Montgomery_Int y3 = y2.mul(y1, ws);
+
+   const Montgomery_Int y1x1 = y1.mul(x1, ws);
+   const Montgomery_Int y1x2 = y1.mul(x2, ws);
+   const Montgomery_Int y1x3 = y1.mul(x3, ws);
+
+   const Montgomery_Int y2x1 = y2.mul(x1, ws);
+   const Montgomery_Int y2x2 = y2.mul(x2, ws);
+   const Montgomery_Int y2x3 = y2.mul(x3, ws);
+
+   const Montgomery_Int y3x1 = y3.mul(x1, ws);
+   const Montgomery_Int y3x2 = y3.mul(x2, ws);
+   const Montgomery_Int y3x3 = y3.mul(x3, ws);
+
+   const Montgomery_Int* M[16] = {
+      &one,
+      &x1,                    // 0001
+      &x2,                    // 0010
+      &x3,                    // 0011
+      &y1,                    // 0100
+      &y1x1,
+      &y1x2,
+      &y1x3,
+      &y2,                    // 1000
+      &y2x1,
+      &y2x2,
+      &y2x3,
+      &y3,                    // 1100
+      &y3x1,
+      &y3x2,
+      &y3x3
+   };
+
+   Montgomery_Int H = one;
+
+   for(size_t i = 0; i != z_bits; i += 2)
+      {
+      if(i > 0)
+         {
+         H.square_this(ws);
+         H.square_this(ws);
+         }
+
+      const uint8_t z1_b = z1.get_substring(z_bits - i - 2, 2);
+      const uint8_t z2_b = z2.get_substring(z_bits - i - 2, 2);
+
+      const uint8_t z12 = (4*z2_b) + z1_b;
+
+      H.mul_by(*M[z12], ws);
+      }
+
+   return H.value();
+   }
+
 }
 
diff --git a/src/lib/math/numbertheory/monty_exp.h b/src/lib/math/numbertheory/monty_exp.h
index 8c644d221..6eeec8bb4 100644
--- a/src/lib/math/numbertheory/monty_exp.h
+++ b/src/lib/math/numbertheory/monty_exp.h
@@ -32,6 +32,15 @@ monty_precompute(std::shared_ptr<const Montgomery_Params> params_p,
 BigInt monty_execute(const Montgomery_Exponentation_State& precomputed_state,
                      const BigInt& k);
 
+/**
+* Return (x^z1 * y^z2) % p
+*/
+BigInt monty_multi_exp(std::shared_ptr<const Montgomery_Params> params_p,
+                       const BigInt& x,
+                       const BigInt& z1,
+                       const BigInt& y,
+                       const BigInt& z2);
+
 }
 
 #endif