Merge GH #1457 Use faster algorithm for ECC multiplication

20 files changed, 536 insertions, 293 deletions

diff --git a/src/cli/speed.cpp b/src/cli/speed.cpp
index 78cbdccf6..465389dc9 100644
--- a/src/cli/speed.cpp
+++ b/src/cli/speed.cpp
@@ -1229,14 +1229,17 @@ class Speed final : public Command
 
             const Botan::BigInt scalar(rng(), group.get_p_bits());
             const Botan::PointGFp& base_point = group.get_base_point();
-            Botan::Blinded_Point_Multiply scalar_mult(base_point, group.get_order(), 4);
+
+            const Botan::PointGFp_Blinded_Multiplier scalar_mult(base_point);
+
+            std::vector<Botan::BigInt> ws;
 
             while(blinded_mult_timer.under(runtime))
                {
                const Botan::PointGFp r1 = mult_timer.run([&]() { return base_point * scalar; });
 
                const Botan::PointGFp r2 = blinded_mult_timer.run(
-               [&]() { return scalar_mult.blinded_multiply(scalar, rng()); });
+                  [&]() { return scalar_mult.mul(scalar, group.get_order(), rng(), ws); });
 
                BOTAN_ASSERT_EQUAL(r1, r2, "Same point computed by both methods");
                }
diff --git a/src/cli/timing_tests.cpp b/src/cli/timing_tests.cpp
index 86c1572ac..3e267829f 100644
--- a/src/cli/timing_tests.cpp
+++ b/src/cli/timing_tests.cpp
@@ -44,7 +44,6 @@
 
 #if defined(BOTAN_HAS_ECDSA)
    #include <botan/ecdsa.h>
-   #include <botan/reducer.h>
    #include <botan/numthry.h>
 #endif
 
@@ -251,19 +250,17 @@ class ECDSA_Timing_Test final : public Timing_Test
       ticks measure_critical_function(std::vector<uint8_t> input) override;
 
    private:
+      const Botan::EC_Group m_group;
       const Botan::ECDSA_PrivateKey m_privkey;
-      const Botan::BigInt m_order;
-      Botan::Blinded_Point_Multiply m_base_point;
-      const Botan::BigInt m_x;
-      const Botan::Modular_Reducer m_mod_order;
+      const Botan::BigInt& m_x;
+      std::vector<Botan::BigInt> m_ws;
    };
 
 ECDSA_Timing_Test::ECDSA_Timing_Test(std::string ecgroup)
-   : m_privkey(Timing_Test::timing_test_rng(), Botan::EC_Group(ecgroup))
-   , m_order(m_privkey.domain().get_order())
-   , m_base_point(m_privkey.domain().get_base_point(), m_order)
+   : m_group(ecgroup)
+   , m_privkey(Timing_Test::timing_test_rng(), m_group)
    , m_x(m_privkey.private_value())
-   , m_mod_order(m_order) {}
+   {}
 
 std::vector<uint8_t> ECDSA_Timing_Test::prepare_input(std::string input)
    {
@@ -274,16 +271,94 @@ std::vector<uint8_t> ECDSA_Timing_Test::prepare_input(std::string input)
 ticks ECDSA_Timing_Test::measure_critical_function(std::vector<uint8_t> input)
    {
    const Botan::BigInt k(input.data(), input.size());
-   const Botan::BigInt msg(Timing_Test::timing_test_rng(), m_order.bits());
+   const Botan::BigInt msg(5); // fixed message to minimize noise
 
    ticks start = get_ticks();
 
    //The following ECDSA operations involve and should not leak any information about k.
-   const Botan::PointGFp k_times_P = m_base_point.blinded_multiply(k, Timing_Test::timing_test_rng());
-   const Botan::BigInt r = m_mod_order.reduce(k_times_P.get_affine_x());
-   const Botan::BigInt s = m_mod_order.multiply(inverse_mod(k, m_order), mul_add(m_x, r, msg));
-   BOTAN_UNUSED(r);
-   BOTAN_UNUSED(s);
+
+   const Botan::BigInt k_inv = Botan::inverse_mod(k, m_group.get_order());
+   const Botan::PointGFp k_times_P = m_group.blinded_base_point_multiply(k, Timing_Test::timing_test_rng(), m_ws);
+   const Botan::BigInt r = m_group.mod_order(k_times_P.get_affine_x());
+   const Botan::BigInt s = m_group.multiply_mod_order(k_inv, mul_add(m_x, r, msg));
+
+   BOTAN_UNUSED(r, s);
+
+   ticks end = get_ticks();
+
+   return (end - start);
+   }
+
+#endif
+
+#if defined(BOTAN_HAS_EC_CURVE_GFP)
+
+class ECC_Mul_Timing_Test final : public Timing_Test
+   {
+   public:
+      ECC_Mul_Timing_Test(std::string ecgroup) :
+         m_group(ecgroup)
+         {}
+
+      std::vector<uint8_t> prepare_input(std::string input) override;
+      ticks measure_critical_function(std::vector<uint8_t> input) override;
+
+   private:
+      const Botan::EC_Group m_group;
+      std::vector<Botan::BigInt> m_ws;
+   };
+
+std::vector<uint8_t> ECC_Mul_Timing_Test::prepare_input(std::string input)
+   {
+   const std::vector<uint8_t> input_vector = Botan::hex_decode(input);
+   return input_vector;
+   }
+
+ticks ECC_Mul_Timing_Test::measure_critical_function(std::vector<uint8_t> input)
+   {
+   const Botan::BigInt k(input.data(), input.size());
+
+   ticks start = get_ticks();
+
+   const Botan::PointGFp k_times_P = m_group.blinded_base_point_multiply(k, Timing_Test::timing_test_rng(), m_ws);
+
+   ticks end = get_ticks();
+
+   return (end - start);
+   }
+
+#endif
+
+#if defined(BOTAN_HAS_NUMBERTHEORY)
+
+class Invmod_Timing_Test final : public Timing_Test
+   {
+   public:
+      Invmod_Timing_Test(size_t p_bits)
+         {
+         m_p = Botan::random_prime(timing_test_rng(), p_bits);
+         }
+
+      std::vector<uint8_t> prepare_input(std::string input) override;
+      ticks measure_critical_function(std::vector<uint8_t> input) override;
+
+   private:
+      Botan::BigInt m_p;
+   };
+
+std::vector<uint8_t> Invmod_Timing_Test::prepare_input(std::string input)
+   {
+   const std::vector<uint8_t> input_vector = Botan::hex_decode(input);
+   return input_vector;
+   }
+
+ticks Invmod_Timing_Test::measure_critical_function(std::vector<uint8_t> input)
+   {
+   const Botan::BigInt k(input.data(), input.size());
+
+   ticks start = get_ticks();
+
+   const Botan::BigInt inv = inverse_mod(k, m_p);
 
    ticks end = get_ticks();
 
@@ -458,6 +533,20 @@ std::unique_ptr<Timing_Test> Timing_Test_Command::lookup_timing_test(const std::
       }
 #endif
 
+#if defined(BOTAN_HAS_EC_CURVE_GFP)
+   if(test_type == "ecc_mul")
+      {
+      return std::unique_ptr<Timing_Test>(new ECC_Mul_Timing_Test("brainpool512r1"));
+      }
+#endif
+
+#if defined(BOTAN_HAS_NUMBERTHEORY)
+   if(test_type == "inverse_mod")
+      {
+      return std::unique_ptr<Timing_Test>(new Invmod_Timing_Test(512));
+      }
+#endif
+
 #if defined(BOTAN_HAS_TLS_CBC)
    if(test_type == "lucky13sha1sec3" || test_type == "lucky13sha1sec4")
       {
diff --git a/src/fuzzer/ecc_helper.h b/src/fuzzer/ecc_helper.h
index b427bc976..cd0ef0ac9 100644
--- a/src/fuzzer/ecc_helper.h
+++ b/src/fuzzer/ecc_helper.h
@@ -8,7 +8,6 @@
 #define ECC_HELPERS_H_
 
 #include "fuzzers.h"
-#include <botan/curve_gfp.h>
 #include <botan/ec_group.h>
 #include <botan/reducer.h>
 
@@ -25,7 +24,10 @@ void check_ecc_math(const Botan::EC_Group& group,
    {
    // These depend only on the group, which is also static
    static const Botan::PointGFp base_point = group.get_base_point();
-   static Botan::Blinded_Point_Multiply blind(base_point, group.get_order(), 4);
+   static Botan::PointGFp_Blinded_Multiplier blind(base_point);
+
+   // This is shared across runs to reduce overhead
+   static std::vector<Botan::BigInt> ws(10);
 
    const size_t hlen = len / 2;
    const Botan::BigInt a = Botan::BigInt::decode(in, hlen);
@@ -42,9 +44,9 @@ void check_ecc_math(const Botan::EC_Group& group,
 
    FUZZER_ASSERT_EQUAL(A1, A2);
 
-   const Botan::PointGFp P1 = blind.blinded_multiply(a, fuzzer_rng());
-   const Botan::PointGFp Q1 = blind.blinded_multiply(b, fuzzer_rng());
-   const Botan::PointGFp R1 = blind.blinded_multiply(c, fuzzer_rng());
+   const Botan::PointGFp P1 = blind.mul(a, group.get_order(), fuzzer_rng(), ws);
+   const Botan::PointGFp Q1 = blind.mul(b, group.get_order(), fuzzer_rng(), ws);
+   const Botan::PointGFp R1 = blind.mul(c, group.get_order(), fuzzer_rng(), ws);
 
    const Botan::PointGFp S1 = P1 + Q1;
    const Botan::PointGFp S2 = Q1 + P1;
diff --git a/src/lib/math/bigint/bigint.cpp b/src/lib/math/bigint/bigint.cpp
index caa78de45..b63bd9be8 100644
--- a/src/lib/math/bigint/bigint.cpp
+++ b/src/lib/math/bigint/bigint.cpp
@@ -119,7 +119,7 @@ int32_t BigInt::cmp(const BigInt& other, bool check_signs) const
 uint32_t BigInt::get_substring(size_t offset, size_t length) const
    {
    if(length > 32)
-      throw Invalid_Argument("BigInt::get_substring: Substring size too big");
+      throw Invalid_Argument("BigInt::get_substring: Substring size " + std::to_string(length) + " too big");
 
    uint64_t piece = 0;
    for(size_t i = 0; i != 8; ++i)
diff --git a/src/lib/math/ec_gfp/point_gfp.cpp b/src/lib/math/ec_gfp/point_gfp.cpp
index f00f030d7..d9599e650 100644
--- a/src/lib/math/ec_gfp/point_gfp.cpp
+++ b/src/lib/math/ec_gfp/point_gfp.cpp
@@ -13,16 +13,16 @@
 
 namespace Botan {
 
-
 PointGFp::PointGFp(const CurveGFp& curve) :
    m_curve(curve),
    m_coord_x(0),
    m_coord_y(1),
    m_coord_z(0)
    {
-   m_curve.to_rep(m_coord_x, m_monty_ws);
-   m_curve.to_rep(m_coord_y, m_monty_ws);
-   m_curve.to_rep(m_coord_z, m_monty_ws);
+   secure_vector<word> monty_ws;
+   m_curve.to_rep(m_coord_x, monty_ws);
+   m_curve.to_rep(m_coord_y, monty_ws);
+   m_curve.to_rep(m_coord_z, monty_ws);
    }
 
 PointGFp::PointGFp(const CurveGFp& curve, const BigInt& x, const BigInt& y) :
@@ -36,9 +36,10 @@ PointGFp::PointGFp(const CurveGFp& curve, const BigInt& x, const BigInt& y) :
    if(y <= 0 || y >= curve.get_p())
       throw Invalid_Argument("Invalid PointGFp affine y");
 
-   m_curve.to_rep(m_coord_x, m_monty_ws);
-   m_curve.to_rep(m_coord_y, m_monty_ws);
-   m_curve.to_rep(m_coord_z, m_monty_ws);
+   secure_vector<word> monty_ws;
+   m_curve.to_rep(m_coord_x, monty_ws);
+   m_curve.to_rep(m_coord_y, monty_ws);
+   m_curve.to_rep(m_coord_z, monty_ws);
    }
 
 void PointGFp::randomize_repr(RandomNumberGenerator& rng)
@@ -49,13 +50,15 @@ void PointGFp::randomize_repr(RandomNumberGenerator& rng)
       while(mask.is_zero())
          mask.randomize(rng, BOTAN_POINTGFP_RANDOMIZE_BLINDING_BITS, false);
 
-      m_curve.to_rep(mask, m_monty_ws);
-      const BigInt mask2 = curve_mult(mask, mask);
-      const BigInt mask3 = curve_mult(mask2, mask);
+      secure_vector<word> monty_ws;
+
+      m_curve.to_rep(mask, monty_ws);
+      const BigInt mask2 = m_curve.mul(mask, mask, monty_ws);
+      const BigInt mask3 = m_curve.mul(mask2, mask, monty_ws);
 
-      m_coord_x = curve_mult(m_coord_x, mask2);
-      m_coord_y = curve_mult(m_coord_y, mask3);
-      m_coord_z = curve_mult(m_coord_z, mask);
+      m_coord_x = m_curve.mul(m_coord_x, mask2, monty_ws);
+      m_coord_y = m_curve.mul(m_coord_y, mask3, monty_ws);
+      m_coord_z = m_curve.mul(m_coord_z, mask, monty_ws);
       }
    }
 
@@ -85,17 +88,23 @@ void PointGFp::add(const PointGFp& rhs, std::vector<BigInt>& ws_bn)
    BigInt& H = ws_bn[6];
    BigInt& r = ws_bn[7];
 
+   secure_vector<word>& monty_ws = ws_bn[8].get_word_vector();
+
    /*
    https://hyperelliptic.org/EFD/g1p/auto-shortw-jacobian-3.html#addition-add-1998-cmo-2
    */
 
-   curve_sqr(rhs_z2, rhs.m_coord_z);
-   curve_mult(U1, m_coord_x, rhs_z2);
-   curve_mult(S1, m_coord_y, curve_mult(rhs.m_coord_z, rhs_z2));
+   m_curve.sqr(rhs_z2, rhs.m_coord_z, monty_ws);
+   m_curve.mul(U1, m_coord_x, rhs_z2, monty_ws);
+   m_curve.mul(S1, m_coord_y,
+               m_curve.mul(rhs.m_coord_z, rhs_z2, monty_ws),
+               monty_ws);
 
-   curve_sqr(lhs_z2, m_coord_z);
-   curve_mult(U2, rhs.m_coord_x, lhs_z2);
-   curve_mult(S2, rhs.m_coord_y, curve_mult(m_coord_z, lhs_z2));
+   m_curve.sqr(lhs_z2, m_coord_z, monty_ws);
+   m_curve.mul(U2, rhs.m_coord_x, lhs_z2, monty_ws);
+   m_curve.mul(S2, rhs.m_coord_y,
+               m_curve.mul(m_coord_z, lhs_z2, monty_ws),
+               monty_ws);
 
    H = U2;
    H -= U1;
@@ -122,13 +131,13 @@ void PointGFp::add(const PointGFp& rhs, std::vector<BigInt>& ws_bn)
       return;
       }
 
-   curve_sqr(U2, H);
+   m_curve.sqr(U2, H, monty_ws);
 
-   curve_mult(S2, U2, H);
+   m_curve.mul(S2, U2, H, monty_ws);
 
-   U2 = curve_mult(U1, U2);
+   U2 = m_curve.mul(U1, U2, monty_ws);
 
-   curve_sqr(m_coord_x, r);
+   m_curve.sqr(m_coord_x, r, monty_ws);
    m_coord_x -= S2;
    m_coord_x -= (U2 << 1);
    while(m_coord_x.is_negative())
@@ -138,12 +147,14 @@ void PointGFp::add(const PointGFp& rhs, std::vector<BigInt>& ws_bn)
    if(U2.is_negative())
       U2 += p;
 
-   curve_mult(m_coord_y, r, U2);
-   m_coord_y -= curve_mult(S1, S2);
+   m_curve.mul(m_coord_y, r, U2, monty_ws);
+   m_coord_y -= m_curve.mul(S1, S2, monty_ws);
    if(m_coord_y.is_negative())
       m_coord_y += p;
 
-   curve_mult(m_coord_z, curve_mult(m_coord_z, rhs.m_coord_z), H);
+   m_curve.mul(m_coord_z,
+               m_curve.mul(m_coord_z, rhs.m_coord_z, monty_ws),
+               H, monty_ws);
    }
 
 // *this *= 2
@@ -151,7 +162,8 @@ void PointGFp::mult2(std::vector<BigInt>& ws_bn)
    {
    if(is_zero())
       return;
-   else if(m_coord_y.is_zero())
+
+   if(m_coord_y.is_zero())
       {
       *this = PointGFp(m_curve); // setting myself to zero
       return;
@@ -173,28 +185,30 @@ void PointGFp::mult2(std::vector<BigInt>& ws_bn)
    BigInt& y = ws_bn[7];
    BigInt& z = ws_bn[8];
 
-   curve_sqr(y_2, m_coord_y);
+   secure_vector<word>& monty_ws = ws_bn[9].get_word_vector();
+
+   m_curve.sqr(y_2, m_coord_y, monty_ws);
 
-   curve_mult(S, m_coord_x, y_2);
+   m_curve.mul(S, m_coord_x, y_2, monty_ws);
    S <<= 2; // * 4
    while(S >= p)
       S -= p;
 
-   curve_sqr(z4, curve_sqr(m_coord_z));
-   curve_mult(a_z4, m_curve.get_a_rep(), z4);
+   m_curve.sqr(z4, m_curve.sqr(m_coord_z, monty_ws), monty_ws);
+   m_curve.mul(a_z4, m_curve.get_a_rep(), z4, monty_ws);
 
-   M = curve_sqr(m_coord_x);
+   M = m_curve.sqr(m_coord_x, monty_ws);
    M *= 3;
    M += a_z4;
    while(M >= p)
       M -= p;
 
-   curve_sqr(x, M);
+   m_curve.sqr(x, M, monty_ws);
    x -= (S << 1);
    while(x.is_negative())
       x += p;
 
-   curve_sqr(U, y_2);
+   m_curve.sqr(U, y_2, monty_ws);
    U <<= 3;
    while(U >= p)
       U -= p;
@@ -203,12 +217,12 @@ void PointGFp::mult2(std::vector<BigInt>& ws_bn)
    while(S.is_negative())
       S += p;
 
-   curve_mult(y, M, S);
+   m_curve.mul(y, M, S, monty_ws);
    y -= U;
    if(y.is_negative())
       y += p;
 
-   curve_mult(z, m_coord_y, m_coord_z);
+   m_curve.mul(z, m_coord_y, m_coord_z, monty_ws);
    z <<= 1;
    if(z >= p)
       z -= p;
@@ -221,7 +235,7 @@ void PointGFp::mult2(std::vector<BigInt>& ws_bn)
 // arithmetic operators
 PointGFp& PointGFp::operator+=(const PointGFp& rhs)
    {
-   std::vector<BigInt> ws(9);
+   std::vector<BigInt> ws(PointGFp::WORKSPACE_SIZE);
    add(rhs, ws);
    return *this;
    }
@@ -247,28 +261,28 @@ PointGFp& PointGFp::operator*=(const BigInt& scalar)
 PointGFp multi_exponentiate(const PointGFp& p1, const BigInt& z1,
                             const PointGFp& p2, const BigInt& z2)
    {
-   const PointGFp p3 = p1 + p2;
+   PointGFp H = p1.zero();
+   const size_t z_bits = std::max(z1.bits(), z2.bits());
 
-   PointGFp H(p1.get_curve()); // create as zero
-   size_t bits_left = std::max(z1.bits(), z2.bits());
+   std::vector<BigInt> ws(PointGFp::WORKSPACE_SIZE);
 
-   std::vector<BigInt> ws(9);
+   PointGFp M[4] = {
+      p1.zero(),
+      p1,
+      p2,
+      p1 + p2,
+   };
 
-   while(bits_left)
+   for(size_t i = 0; i != z_bits; ++i)
       {
       H.mult2(ws);
 
-      const bool z1_b = z1.get_bit(bits_left - 1);
-      const bool z2_b = z2.get_bit(bits_left - 1);
+      const uint8_t z1_b = z1.get_bit(z_bits - i - 1);
+      const uint8_t z2_b = z2.get_bit(z_bits - i - 1);
 
-      if(z1_b == true && z2_b == true)
-         H.add(p3, ws);
-      else if(z1_b)
-         H.add(p1, ws);
-      else if(z2_b)
-         H.add(p2, ws);
+      const uint8_t z12 = (2*z2_b) + z1_b;
 
-      --bits_left;
+      H.add(M[z12], ws);
       }
 
    if(z1.is_negative() != z2.is_negative())
@@ -281,13 +295,11 @@ PointGFp operator*(const BigInt& scalar, const PointGFp& point)
    {
    //BOTAN_ASSERT(point.on_the_curve(), "Input is on the curve");
 
-   const CurveGFp& curve = point.get_curve();
-
    const size_t scalar_bits = scalar.bits();
 
-   std::vector<BigInt> ws(9);
+   std::vector<BigInt> ws(PointGFp::WORKSPACE_SIZE);
 
-   PointGFp R[2] = { PointGFp(curve), point };
+   PointGFp R[2] = { point.zero(), point };
 
    for(size_t i = scalar_bits; i > 0; i--)
       {
@@ -304,98 +316,17 @@ PointGFp operator*(const BigInt& scalar, const PointGFp& point)
    return R[0];
    }
 
-Blinded_Point_Multiply::Blinded_Point_Multiply(const PointGFp& base, const BigInt& order, size_t h) :
-   m_h(h > 0 ? h : 4), m_order(order), m_ws(9)
-   {
-   // Upper bound is a sanity check rather than hard limit
-   if(m_h < 1 || m_h > 8)
-      throw Invalid_Argument("Blinded_Point_Multiply invalid h param");
-
-   const CurveGFp& curve = base.get_curve();
-
-   const PointGFp inv = -base;
-
-   m_U.resize(6*m_h + 3);
-
-   m_U[3*m_h+0] = inv;
-   m_U[3*m_h+1] = PointGFp::zero_of(curve);
-   m_U[3*m_h+2] = base;
-
-   for(size_t i = 1; i <= 3 * m_h + 1; ++i)
-      {
-      m_U[3*m_h+1+i] = m_U[3*m_h+i];
-      m_U[3*m_h+1+i].add(base, m_ws);
-
-      m_U[3*m_h+1-i] = m_U[3*m_h+2-i];
-      m_U[3*m_h+1-i].add(inv, m_ws);
-      }
-   }
-
-PointGFp Blinded_Point_Multiply::blinded_multiply(const BigInt& scalar_in,
-                                                  RandomNumberGenerator& rng)
-   {
-   if(scalar_in.is_negative())
-      throw Invalid_Argument("Blinded_Point_Multiply scalar must be positive");
-
-#if BOTAN_POINTGFP_USE_SCALAR_BLINDING
-   // Choose a small mask m and use k' = k + m*order (Coron's 1st countermeasure)
-   const BigInt mask(rng, (m_order.bits()+1)/2, false);
-   const BigInt scalar = scalar_in + m_order * mask;
-#else
-   const BigInt& scalar = scalar_in;
-#endif
-
-   const size_t scalar_bits = scalar.bits();
-
-   // Randomize each point representation (Coron's 3rd countermeasure)
-   for(size_t i = 0; i != m_U.size(); ++i)
-      m_U[i].randomize_repr(rng);
-
-   PointGFp R = m_U.at(3*m_h + 2); // base point
-   int32_t alpha = 0;
-
-   R.randomize_repr(rng);
-
-   /*
-   Algorithm 7 from "Randomizing the Montgomery Powering Ladder"
-   Duc-Phong Le, Chik How Tan and Michael Tunstall
-   https://eprint.iacr.org/2015/657
-
-   It takes a random walk through (a subset of) the set of addition
-   chains that end in k.
-   */
-   for(size_t i = scalar_bits; i > 0; i--)
-      {
-      const int32_t ki = scalar.get_bit(i);
-
-      // choose gamma from -h,...,h
-      const int32_t gamma = static_cast<int32_t>((rng.next_byte() % (2*m_h))) - m_h;
-      const int32_t l = gamma - 2*alpha + ki - (ki ^ 1);
-
-      R.mult2(m_ws);
-      R.add(m_U.at(3*m_h + 1 + l), m_ws);
-      alpha = gamma;
-      }
-
-   const int32_t k0 = scalar.get_bit(0);
-   R.add(m_U[3*m_h + 1 - alpha - (k0 ^ 1)], m_ws);
-
-
-   //BOTAN_ASSERT(R.on_the_curve(), "Output is on the curve");
-
-   return R;
-   }
-
 BigInt PointGFp::get_affine_x() const
    {
    if(is_zero())
       throw Illegal_Transformation("Cannot convert zero point to affine");
 
-   BigInt z2 = curve_sqr(m_coord_z);
-   m_curve.from_rep(z2, m_monty_ws);
+   secure_vector<word> monty_ws;
+   BigInt z2 = m_curve.sqr(m_coord_z, monty_ws);
+   m_curve.from_rep(z2, monty_ws);
    z2 = inverse_mod(z2, m_curve.get_p());
 
-   return curve_mult(z2, m_coord_x);
+   return m_curve.mul(z2, m_coord_x, monty_ws);
    }
 
 BigInt PointGFp::get_affine_y() const
@@ -403,11 +334,12 @@ BigInt PointGFp::get_affine_y() const
    if(is_zero())
       throw Illegal_Transformation("Cannot convert zero point to affine");
 
-   BigInt z3 = curve_mult(m_coord_z, curve_sqr(m_coord_z));
+   secure_vector<word> monty_ws;
+   BigInt z3 = m_curve.mul(m_coord_z, m_curve.sqr(m_coord_z, monty_ws), monty_ws);
    z3 = inverse_mod(z3, m_curve.get_p());
-   m_curve.to_rep(z3, m_monty_ws);
+   m_curve.to_rep(z3, monty_ws);
 
-   return curve_mult(z3, m_coord_y);
+   return m_curve.mul(z3, m_coord_y, monty_ws);
    }
 
 bool PointGFp::on_the_curve() const
@@ -421,22 +353,24 @@ bool PointGFp::on_the_curve() const
    if(is_zero())
       return true;
 
-   const BigInt y2 = m_curve.from_rep(curve_sqr(m_coord_y), m_monty_ws);
-   const BigInt x3 = curve_mult(m_coord_x, curve_sqr(m_coord_x));
-   const BigInt ax = curve_mult(m_coord_x, m_curve.get_a_rep());
-   const BigInt z2 = curve_sqr(m_coord_z);
+   secure_vector<word> monty_ws;
+
+   const BigInt y2 = m_curve.from_rep(m_curve.sqr(m_coord_y, monty_ws), monty_ws);
+   const BigInt x3 = m_curve.mul(m_coord_x, m_curve.sqr(m_coord_x, monty_ws), monty_ws);
+   const BigInt ax = m_curve.mul(m_coord_x, m_curve.get_a_rep(), monty_ws);
+   const BigInt z2 = m_curve.sqr(m_coord_z, monty_ws);
 
    if(m_coord_z == z2) // Is z equal to 1 (in Montgomery form)?
       {
-      if(y2 != m_curve.from_rep(x3 + ax + m_curve.get_b_rep(), m_monty_ws))
+      if(y2 != m_curve.from_rep(x3 + ax + m_curve.get_b_rep(), monty_ws))
          return false;
       }
 
-   const BigInt z3 = curve_mult(m_coord_z, z2);
-   const BigInt ax_z4 = curve_mult(ax, curve_sqr(z2));
-   const BigInt b_z6 = curve_mult(m_curve.get_b_rep(), curve_sqr(z3));
+   const BigInt z3 = m_curve.mul(m_coord_z, z2, monty_ws);
+   const BigInt ax_z4 = m_curve.mul(ax, m_curve.sqr(z2, monty_ws), monty_ws);
+   const BigInt b_z6 = m_curve.mul(m_curve.get_b_rep(), m_curve.sqr(z3, monty_ws), monty_ws);
 
-   if(y2 != m_curve.from_rep(x3 + ax_z4 + b_z6, m_monty_ws))
+   if(y2 != m_curve.from_rep(x3 + ax_z4 + b_z6, monty_ws))
       return false;
 
    return true;
@@ -449,12 +383,11 @@ void PointGFp::swap(PointGFp& other)
    m_coord_x.swap(other.m_coord_x);
    m_coord_y.swap(other.m_coord_y);
    m_coord_z.swap(other.m_coord_z);
-   m_monty_ws.swap(other.m_monty_ws);
    }
 
 bool PointGFp::operator==(const PointGFp& other) const
    {
-   if(get_curve() != other.get_curve())
+   if(m_curve != other.m_curve)
       return false;
 
    // If this is zero, only equal if other is also zero
diff --git a/src/lib/math/ec_gfp/point_gfp.h b/src/lib/math/ec_gfp/point_gfp.h
index 2e1c626bd..a2c7ae86a 100644
--- a/src/lib/math/ec_gfp/point_gfp.h
+++ b/src/lib/math/ec_gfp/point_gfp.h
@@ -49,6 +49,8 @@ class BOTAN_PUBLIC_API(2,0) PointGFp final
          HYBRID       = 2
       };
 
+      enum { WORKSPACE_SIZE = 10 };
+
       /**
       * Construct an uninitialized PointGFp
       */
@@ -60,11 +62,6 @@ class BOTAN_PUBLIC_API(2,0) PointGFp final
       */
       explicit PointGFp(const CurveGFp& curve);
 
-      static PointGFp zero_of(const CurveGFp& curve)
-         {
-         return PointGFp(curve);
-         }
-
       /**
       * Copy constructor
       */
@@ -120,30 +117,9 @@ class BOTAN_PUBLIC_API(2,0) PointGFp final
       * @param scalar the PointGFp to multiply with *this
       * @result resulting PointGFp
       */
-
       PointGFp& operator*=(const BigInt& scalar);
 
       /**
-      * Multiplication Operator
-      * @param scalar the scalar value
-      * @param point the point value
-      * @return scalar*point on the curve
-      */
-      friend BOTAN_PUBLIC_API(2,0) PointGFp operator*(const BigInt& scalar, const PointGFp& point);
-
-      /**
-      * Multiexponentiation
-      * @param p1 a point
-      * @param z1 a scalar
-      * @param p2 a point
-      * @param z2 a scalar
-      * @result (p1 * z1 + p2 * z2)
-      */
-      friend BOTAN_PUBLIC_API(2,0) PointGFp multi_exponentiate(
-        const PointGFp& p1, const BigInt& z1,
-        const PointGFp& p2, const BigInt& z2);
-
-      /**
       * Negate this point
       * @return *this
       */
@@ -202,50 +178,49 @@ class BOTAN_PUBLIC_API(2,0) PointGFp final
       * Equality operator
       */
       bool operator==(const PointGFp& other) const;
-   private:
-      friend class Blinded_Point_Multiply;
-
-      BigInt curve_mult(const BigInt& x, const BigInt& y) const
-         {
-         BigInt z;
-         m_curve.mul(z, x, y, m_monty_ws);
-         return z;
-         }
-
-      void curve_mult(BigInt& z, const BigInt& x, const BigInt& y) const
-         {
-         m_curve.mul(z, x, y, m_monty_ws);
-         }
-
-      BigInt curve_sqr(const BigInt& x) const
-         {
-         BigInt z;
-         m_curve.sqr(z, x, m_monty_ws);
-         return z;
-         }
-
-      void curve_sqr(BigInt& z, const BigInt& x) const
-         {
-         m_curve.sqr(z, x, m_monty_ws);
-         }
 
       /**
       * Point addition
-      * @param workspace temp space, at least 11 elements
+      * @param workspace temp space, at least 9 elements
       */
       void add(const PointGFp& other, std::vector<BigInt>& workspace);
 
       /**
       * Point doubling
-      * @param workspace temp space, at least 9 elements
+      * @param workspace temp space, at least 10 elements
       */
       void mult2(std::vector<BigInt>& workspace);
 
+      /**
+      * Return the zero (aka infinite) point associated with this curve
+      */
+      PointGFp zero() const { return PointGFp(m_curve); }
+
+   private:
       CurveGFp m_curve;
       BigInt m_coord_x, m_coord_y, m_coord_z;
-      mutable secure_vector<word> m_monty_ws; // workspace for Montgomery
    };
 
+/**
+* Point multiplication operator
+* @param scalar the scalar value
+* @param point the point value
+* @return scalar*point on the curve
+*/
+BOTAN_PUBLIC_API(2,0) PointGFp operator*(const BigInt& scalar, const PointGFp& point);
+
+/**
+* ECC point multiexponentiation
+* @param p1 a point
+* @param z1 a scalar
+* @param p2 a point
+* @param z2 a scalar
+* @result (p1 * z1 + p2 * z2)
+*/
+BOTAN_PUBLIC_API(2,0) PointGFp multi_exponentiate(
+   const PointGFp& p1, const BigInt& z1,
+   const PointGFp& p2, const BigInt& z2);
+
 // relational operators
 inline bool operator!=(const PointGFp& lhs, const PointGFp& rhs)
    {
@@ -299,19 +274,68 @@ PointGFp OS2ECP(const std::vector<uint8_t, Alloc>& data, const CurveGFp& curve)
    { return OS2ECP(data.data(), data.size(), curve); }
 
 /**
-
+* Blinded ECC point multiplication
 */
-class BOTAN_PUBLIC_API(2,0) Blinded_Point_Multiply final
+class BOTAN_PUBLIC_API(2,5) PointGFp_Blinded_Multiplier final
+   {
+   public:
+      /**
+      * @param base_point the point that will be multiplied (eg, the base point of the curve)
+      * @param w the window bits (leave as zero to take a default value)
+      */
+      PointGFp_Blinded_Multiplier(const PointGFp& base_point,
+                                  size_t w = 0);
+
+      /**
+      * @param base_point the point that will be multiplied (eg, the base point of the curve)
+      * @param ws a temporary workspace
+      * @param w the window bits (leave as zero to take a default value)
+      */
+      PointGFp_Blinded_Multiplier(const PointGFp& base_point,
+                                  std::vector<BigInt>& ws,
+                                  size_t w = 0);
+
+      /**
+      * Randomize the internal state. Changing the values may provide
+      * some protection against side channel attacks.
+      * @param rng a random number generator
+      */
+      void randomize(RandomNumberGenerator& rng);
+
+      /**
+      * Perform blinded point multiplication
+      * @param k the scalar
+      * @param group_order the order of the group
+      * @param rng a random number generator
+      * @param ws a temporary workspace
+      * @return base_point*k
+      */
+      PointGFp mul(const BigInt& k,
+                   const BigInt& group_order,
+                   RandomNumberGenerator& rng,
+                   std::vector<BigInt>& ws) const;
+   private:
+      void init(const PointGFp& base_point, size_t w, std::vector<BigInt>& ws);
+
+      std::vector<PointGFp> m_U;
+      size_t m_h;
+   };
+
+class BOTAN_PUBLIC_API(2,0) BOTAN_DEPRECATED("Use PointGFp_Blinded_Multiplier") Blinded_Point_Multiply final
    {
    public:
-      Blinded_Point_Multiply(const PointGFp& base, const BigInt& order, size_t h = 0);
+      Blinded_Point_Multiply(const PointGFp& base, const BigInt& order, size_t h = 0) :
+         m_ws(PointGFp::WORKSPACE_SIZE), m_order(order), m_point_mul(base, m_ws, h) {}
 
-      PointGFp blinded_multiply(const BigInt& scalar, RandomNumberGenerator& rng);
+      PointGFp blinded_multiply(const BigInt& scalar, RandomNumberGenerator& rng)
+         {
+         m_point_mul.randomize(rng);
+         return m_point_mul.mul(scalar, m_order, rng, m_ws);
+         }
    private:
-      const size_t m_h;
-      const BigInt& m_order;
       std::vector<BigInt> m_ws;
-      std::vector<PointGFp> m_U;
+      const BigInt& m_order;
+      PointGFp_Blinded_Multiplier m_point_mul;
    };
 
 }
diff --git a/src/lib/math/ec_gfp/point_mul.cpp b/src/lib/math/ec_gfp/point_mul.cpp
new file mode 100644
index 000000000..7fea02280
--- /dev/null
+++ b/src/lib/math/ec_gfp/point_mul.cpp
@@ -0,0 +1,113 @@
+/*
+* (C) 2015,2018 Jack Lloyd
+*
+* Botan is released under the Simplified BSD License (see license.txt)
+*/
+
+#include <botan/point_gfp.h>
+#include <botan/rng.h>
+#include <botan/internal/rounding.h>
+
+namespace Botan {
+
+PointGFp_Blinded_Multiplier::PointGFp_Blinded_Multiplier(const PointGFp& base,
+                                                         std::vector<BigInt>& ws,
+                                                         size_t w)
+   {
+   init(base, w, ws);
+   }
+
+PointGFp_Blinded_Multiplier::PointGFp_Blinded_Multiplier(const PointGFp& base,
+                                                         size_t w)
+   {
+   std::vector<BigInt> ws(9);
+   init(base, w, ws);
+   }
+
+void PointGFp_Blinded_Multiplier::init(const PointGFp& base,
+                                       size_t w,
+                                       std::vector<BigInt>& ws)
+   {
+   m_h = (w == 0 ? 5 : w);
+
+   if(ws.size() < PointGFp::WORKSPACE_SIZE)
+      ws.resize(PointGFp::WORKSPACE_SIZE);
+
+   // Upper bound is a sanity check rather than hard limit
+   if(m_h < 1 || m_h > 8)
+      throw Invalid_Argument("PointGFp_Blinded_Multiplier invalid w param");
+
+   m_U.resize(1 << m_h);
+   m_U[0] = base.zero();
+   m_U[1] = base;
+
+   for(size_t i = 2; i < m_U.size(); ++i)
+      {
+      m_U[i] = m_U[i-1];
+      m_U[i].add(base, ws);
+      }
+   }
+
+void PointGFp_Blinded_Multiplier::randomize(RandomNumberGenerator& rng)
+   {
+   // Randomize each point representation (Coron's 3rd countermeasure)
+   for(size_t i = 0; i != m_U.size(); ++i)
+      m_U[i].randomize_repr(rng);
+   }
+
+PointGFp PointGFp_Blinded_Multiplier::mul(const BigInt& k,
+                                          const BigInt& group_order,
+                                          RandomNumberGenerator& rng,
+                                          std::vector<BigInt>& ws) const
+   {
+   if(k.is_negative())
+      throw Invalid_Argument("PointGFp_Blinded_Multiplier scalar must be positive");
+
+#if BOTAN_POINTGFP_USE_SCALAR_BLINDING
+   // Choose a small mask m and use k' = k + m*order (Coron's 1st countermeasure)
+   const BigInt mask(rng, group_order.bits() / 4, false);
+   const BigInt scalar = k + group_order * mask;
+#else
+   const BigInt& scalar = k;
+#endif
+
+   if(ws.size() < PointGFp::WORKSPACE_SIZE)
+      ws.resize(PointGFp::WORKSPACE_SIZE);
+
+   const size_t scalar_bits = scalar.bits();
+
+   size_t windows = round_up(scalar_bits, m_h) / m_h;
+
+   PointGFp R = m_U[0];
+
+   if(windows > 0)
+      {
+      windows--;
+      const uint32_t nibble = scalar.get_substring(windows*m_h, m_h);
+      R.add(m_U[nibble], ws);
+
+      /*
+      Randomize after adding the first nibble as before the addition R
+      is zero, and we cannot effectively randomize the point
+      representation of the zero point.
+      */
+      R.randomize_repr(rng);
+
+      while(windows)
+         {
+         for(size_t i = 0; i != m_h; ++i)
+            R.mult2(ws);
+
+         const uint32_t inner_nibble = scalar.get_substring((windows-1)*m_h, m_h);
+         // cache side channel here, we are relying on blinding...
+         R.add(m_U[inner_nibble], ws);
+         windows--;
+         }
+      }
+
+   //BOTAN_ASSERT(R.on_the_curve(), "Output is on the curve");
+
+   return R;
+   }
+
+}
diff --git a/src/lib/pubkey/ec_group/ec_group.cpp b/src/lib/pubkey/ec_group/ec_group.cpp
index 771bd4b0f..102eed6e5 100644
--- a/src/lib/pubkey/ec_group/ec_group.cpp
+++ b/src/lib/pubkey/ec_group/ec_group.cpp
@@ -17,6 +17,10 @@
 #include <botan/mutex.h>
 #include <vector>
 
+#if defined(BOTAN_HAS_SYSTEM_RNG)
+   #include <botan/system_rng.h>
+#endif
+
 namespace Botan {
 
 class EC_Group_Data final
@@ -36,19 +40,28 @@ class EC_Group_Data final
          m_order(order),
          m_cofactor(cofactor),
          m_mod_order(order),
+         m_base_mult(m_base_point, 5),
          m_oid(oid),
          m_p_bits(p.bits()),
-         m_order_bits(order.bits())
+         m_order_bits(order.bits()),
+         m_a_is_minus_3(a == p - 3)
          {
+#if defined(BOTAN_HAS_SYSTEM_RNG)
+         m_base_mult.randomize(system_rng());
+#endif
          }
 
       bool match(const BigInt& p, const BigInt& a, const BigInt& b,
                  const BigInt& g_x, const BigInt& g_y,
                  const BigInt& order, const BigInt& cofactor) const
          {
-         return (this->p() == p && this->a() == a && this->b() == b &&
-                 this->order() == order && this->cofactor() == cofactor &&
-                 this->g_x() == g_x && this->g_y() == g_y);
+         return (this->p() == p &&
+                 this->a() == a &&
+                 this->b() == b &&
+                 this->order() == order &&
+                 this->cofactor() == cofactor &&
+                 this->g_x() == g_x &&
+                 this->g_y() == g_y);
          }
 
       const OID& oid() const { return m_oid; }
@@ -69,6 +82,8 @@ class EC_Group_Data final
       const CurveGFp& curve() const { return m_curve; }
       const PointGFp& base_point() const { return m_base_point; }
 
+      bool a_is_minus_3() const { return m_a_is_minus_3; }
+
       BigInt mod_order(const BigInt& x) const { return m_mod_order.reduce(x); }
 
       BigInt multiply_mod_order(const BigInt& x, const BigInt& y) const
@@ -76,15 +91,24 @@ class EC_Group_Data final
          return m_mod_order.multiply(x, y);
          }
 
+      PointGFp blinded_base_point_multiply(const BigInt& k,
+                                           RandomNumberGenerator& rng,
+                                           std::vector<BigInt>& ws) const
+         {
+         return m_base_mult.mul(k, m_order, rng, ws);
+         }
+
    private:
       CurveGFp m_curve;
       PointGFp m_base_point;
       BigInt m_order;
       BigInt m_cofactor;
       Modular_Reducer m_mod_order;
+      PointGFp_Blinded_Multiplier m_base_mult;
       OID m_oid;
       size_t m_p_bits;
       size_t m_order_bits;
+      bool m_a_is_minus_3;
    };
 
 class EC_Group_Data_Map final
@@ -349,6 +373,11 @@ const CurveGFp& EC_Group::get_curve() const
    return data().curve();
    }
 
+bool EC_Group::a_is_minus_3() const
+   {
+   return data().a_is_minus_3();
+   }
+
 size_t EC_Group::get_p_bits() const
    {
    return data().p_bits();
@@ -421,6 +450,7 @@ PointGFp EC_Group::OS2ECP(const uint8_t bits[], size_t len) const
 
 PointGFp EC_Group::point(const BigInt& x, const BigInt& y) const
    {
+   // TODO: randomize the representation?
    return PointGFp(data().curve(), x, y);
    }
 
@@ -429,6 +459,13 @@ PointGFp EC_Group::point_multiply(const BigInt& x, const PointGFp& pt, const Big
    return multi_exponentiate(get_base_point(), x, pt, y);
    }
 
+PointGFp EC_Group::blinded_base_point_multiply(const BigInt& k,
+                                               RandomNumberGenerator& rng,
+                                               std::vector<BigInt>& ws) const
+   {
+   return data().blinded_base_point_multiply(k, rng, ws);
+   }
+
 PointGFp EC_Group::zero_point() const
    {
    return PointGFp(data().curve());
diff --git a/src/lib/pubkey/ec_group/ec_group.h b/src/lib/pubkey/ec_group/ec_group.h
index a60c71157..2baa2555e 100644
--- a/src/lib/pubkey/ec_group/ec_group.h
+++ b/src/lib/pubkey/ec_group/ec_group.h
@@ -11,7 +11,6 @@
 #define BOTAN_ECC_DOMAIN_PARAMETERS_H_
 
 #include <botan/point_gfp.h>
-#include <botan/curve_gfp.h>
 #include <botan/asn1_oid.h>
 #include <memory>
 #include <set>
@@ -27,6 +26,8 @@ enum EC_Group_Encoding {
    EC_DOMPAR_ENC_OID = 2
 };
 
+class CurveGFp;
+
 class EC_Group_Data;
 class EC_Group_Data_Map;
 
@@ -126,6 +127,11 @@ class BOTAN_PUBLIC_API(2,0) EC_Group final
       BOTAN_DEPRECATED("Avoid CurveGFp") const CurveGFp& get_curve() const;
 
       /**
+      * Return if a == -3 mod p
+      */
+      bool a_is_minus_3() const;
+
+      /**
       * Return the size of p in bits (same as get_p().bits())
       */
       size_t get_p_bits() const;
@@ -206,12 +212,21 @@ class BOTAN_PUBLIC_API(2,0) EC_Group final
       PointGFp point(const BigInt& x, const BigInt& y) const;
 
       /**
-      * Multi exponentiate
+      * Multi exponentiate. Not constant time.
       * @return base_point*x + pt*y
       */
       PointGFp point_multiply(const BigInt& x, const PointGFp& pt, const BigInt& y) const;
 
       /**
+      * Blinded point multiplication, attempts resistance to side channels
+      * @param k the scalar
+      * @param rng a random number generator
+      * @param ws a temp workspace
+      * @return base_point*k
+      */
+      PointGFp blinded_base_point_multiply(const BigInt& k, RandomNumberGenerator& rng, std::vector<BigInt>& ws) const;
+
+      /**
       * Return the zero (or infinite) point on this curve
       */
       PointGFp zero_point() const;
diff --git a/src/lib/pubkey/ecdh/ecdh.cpp b/src/lib/pubkey/ecdh/ecdh.cpp
index 1850696e1..4989fa0a5 100644
--- a/src/lib/pubkey/ecdh/ecdh.cpp
+++ b/src/lib/pubkey/ecdh/ecdh.cpp
@@ -28,26 +28,30 @@ class ECDH_KA_Operation final : public PK_Ops::Key_Agreement_with_KDF
 
       ECDH_KA_Operation(const ECDH_PrivateKey& key, const std::string& kdf, RandomNumberGenerator& rng) :
          PK_Ops::Key_Agreement_with_KDF(kdf),
-         m_domain(key.domain()),
+         m_group(key.domain()),
          m_rng(rng)
          {
-         m_l_times_priv = inverse_mod(m_domain.get_cofactor(), m_domain.get_order()) * key.private_value();
+         m_l_times_priv = inverse_mod(m_group.get_cofactor(), m_group.get_order()) * key.private_value();
          }
 
       secure_vector<uint8_t> raw_agree(const uint8_t w[], size_t w_len) override
          {
-         PointGFp point = m_domain.OS2ECP(w, w_len);
-         PointGFp S = m_domain.get_cofactor() * point;
-         Blinded_Point_Multiply blinder(S, m_domain.get_order());
-         S = blinder.blinded_multiply(m_l_times_priv, m_rng);
-         BOTAN_ASSERT(S.on_the_curve(), "ECDH agreed value was on the curve");
-         return BigInt::encode_1363(S.get_affine_x(), m_domain.get_p_bytes());
+         PointGFp input_point = m_group.get_cofactor() * m_group.OS2ECP(w, w_len);
+         input_point.randomize_repr(m_rng);
+
+         PointGFp_Blinded_Multiplier blinder(input_point, m_ws);
+
+         const PointGFp S = blinder.mul(m_l_times_priv, m_group.get_order(), m_rng, m_ws);
+
+         if(S.on_the_curve() == false)
+            throw Internal_Error("ECDH agreed value was not on the curve");
+         return BigInt::encode_1363(S.get_affine_x(), m_group.get_p_bytes());
          }
    private:
-      const EC_Group& m_domain;
+      const EC_Group m_group;
       BigInt m_l_times_priv;
       RandomNumberGenerator& m_rng;
-
+      std::vector<BigInt> m_ws;
    };
 
 }
diff --git a/src/lib/pubkey/ecdsa/ecdsa.cpp b/src/lib/pubkey/ecdsa/ecdsa.cpp
index 12ccd9608..57bc197c5 100644
--- a/src/lib/pubkey/ecdsa/ecdsa.cpp
+++ b/src/lib/pubkey/ecdsa/ecdsa.cpp
@@ -53,7 +53,6 @@ class ECDSA_Signature_Operation final : public PK_Ops::Signature_with_EMSA
                                 const std::string& emsa) :
          PK_Ops::Signature_with_EMSA(emsa),
          m_group(ecdsa.domain()),
-         m_base_point(m_group.get_base_point(), m_group.get_order()),
          m_x(ecdsa.private_value())
          {
 #if defined(BOTAN_HAS_RFC6979_GENERATOR)
@@ -68,12 +67,13 @@ class ECDSA_Signature_Operation final : public PK_Ops::Signature_with_EMSA
 
    private:
       const EC_Group m_group;
-      Blinded_Point_Multiply m_base_point;
       const BigInt& m_x;
 
 #if defined(BOTAN_HAS_RFC6979_GENERATOR)
       std::string m_rfc6979_hash;
 #endif
+
+      std::vector<BigInt> m_ws;
    };
 
 secure_vector<uint8_t>
@@ -89,7 +89,7 @@ ECDSA_Signature_Operation::raw_sign(const uint8_t msg[], size_t msg_len,
 #endif
 
    const BigInt k_inv = inverse_mod(k, m_group.get_order());
-   const PointGFp k_times_P = m_base_point.blinded_multiply(k, rng);
+   const PointGFp k_times_P = m_group.blinded_base_point_multiply(k, rng, m_ws);
    const BigInt r = m_group.mod_order(k_times_P.get_affine_x());
    const BigInt s = m_group.multiply_mod_order(k_inv, mul_add(m_x, r, m));
 
diff --git a/src/lib/pubkey/ecgdsa/ecgdsa.cpp b/src/lib/pubkey/ecgdsa/ecgdsa.cpp
index f8e5744d9..6cbd3453b 100644
--- a/src/lib/pubkey/ecgdsa/ecgdsa.cpp
+++ b/src/lib/pubkey/ecgdsa/ecgdsa.cpp
@@ -38,7 +38,6 @@ class ECGDSA_Signature_Operation final : public PK_Ops::Signature_with_EMSA
                                 const std::string& emsa) :
          PK_Ops::Signature_with_EMSA(emsa),
          m_group(ecgdsa.domain()),
-         m_base_point(m_group.get_base_point(), m_group.get_order()),
          m_x(ecgdsa.private_value())
          {
          }
@@ -50,8 +49,8 @@ class ECGDSA_Signature_Operation final : public PK_Ops::Signature_with_EMSA
 
    private:
       const EC_Group m_group;
-      Blinded_Point_Multiply m_base_point;
       const BigInt& m_x;
+      std::vector<BigInt> m_ws;
    };
 
 secure_vector<uint8_t>
@@ -62,7 +61,7 @@ ECGDSA_Signature_Operation::raw_sign(const uint8_t msg[], size_t msg_len,
 
    BigInt k = BigInt::random_integer(rng, 1, m_group.get_order());
 
-   const PointGFp k_times_P = m_base_point.blinded_multiply(k, rng);
+   const PointGFp k_times_P = m_group.blinded_base_point_multiply(k, rng, m_ws);
    const BigInt r = m_group.mod_order(k_times_P.get_affine_x());
    const BigInt s = m_group.multiply_mod_order(m_x, mul_sub(k, r, m));
 
diff --git a/src/lib/pubkey/ecies/ecies.cpp b/src/lib/pubkey/ecies/ecies.cpp
index cd09b4c52..1120a850a 100644
--- a/src/lib/pubkey/ecies/ecies.cpp
+++ b/src/lib/pubkey/ecies/ecies.cpp
@@ -66,16 +66,24 @@ class ECIES_ECDH_KA_Operation final : public PK_Ops::Key_Agreement_with_KDF
 
       secure_vector<uint8_t> raw_agree(const uint8_t w[], size_t w_len) override
          {
-         PointGFp point = m_key.domain().OS2ECP(w, w_len);
-         Blinded_Point_Multiply blinder(point, m_key.domain().get_order());
-         PointGFp S = blinder.blinded_multiply(m_key.private_value(), m_rng);
-         BOTAN_ASSERT(S.on_the_curve(), "ECDH agreed value was on the curve");
-         return BigInt::encode_1363(S.get_affine_x(), m_key.domain().get_p_bytes());
+         const EC_Group& group = m_key.domain();
+
+         PointGFp input_point = group.OS2ECP(w, w_len);
+         input_point.randomize_repr(m_rng);
+
+         PointGFp_Blinded_Multiplier blinder(input_point, m_ws);
+
+         const PointGFp S = blinder.mul(m_key.private_value(), group.get_order(), m_rng, m_ws);
+
+         if(S.on_the_curve() == false)
+            throw Internal_Error("ECDH agreed value was not on the curve");
+         return BigInt::encode_1363(S.get_affine_x(), group.get_p_bytes());
          }
 
    private:
       ECIES_PrivateKey m_key;
       RandomNumberGenerator& m_rng;
+      std::vector<BigInt> m_ws;
    };
 
 std::unique_ptr<PK_Ops::Key_Agreement>
diff --git a/src/lib/pubkey/eckcdsa/eckcdsa.cpp b/src/lib/pubkey/eckcdsa/eckcdsa.cpp
index 743d5ab95..be721a6b6 100644
--- a/src/lib/pubkey/eckcdsa/eckcdsa.cpp
+++ b/src/lib/pubkey/eckcdsa/eckcdsa.cpp
@@ -45,7 +45,6 @@ class ECKCDSA_Signature_Operation final : public PK_Ops::Signature_with_EMSA
                                 const std::string& emsa) :
          PK_Ops::Signature_with_EMSA(emsa),
          m_group(eckcdsa.domain()),
-         m_base_point(m_group.get_base_point(), m_group.get_order()),
          m_x(eckcdsa.private_value()),
          m_prefix()
          {
@@ -68,9 +67,9 @@ class ECKCDSA_Signature_Operation final : public PK_Ops::Signature_with_EMSA
 
    private:
       const EC_Group m_group;
-      Blinded_Point_Multiply m_base_point;
       const BigInt& m_x;
       secure_vector<uint8_t> m_prefix;
+      std::vector<BigInt> m_ws;
    };
 
 secure_vector<uint8_t>
@@ -78,7 +77,7 @@ ECKCDSA_Signature_Operation::raw_sign(const uint8_t msg[], size_t,
                                      RandomNumberGenerator& rng)
    {
    const BigInt k = BigInt::random_integer(rng, 1, m_group.get_order());
-   const PointGFp k_times_P = m_base_point.blinded_multiply(k, rng);
+   const PointGFp k_times_P = m_group.blinded_base_point_multiply(k, rng, m_ws);
    const BigInt k_times_P_x = k_times_P.get_affine_x();
 
    secure_vector<uint8_t> to_be_hashed(k_times_P_x.bytes());
diff --git a/src/lib/pubkey/gost_3410/gost_3410.cpp b/src/lib/pubkey/gost_3410/gost_3410.cpp
index 760e667aa..79d3f204d 100644
--- a/src/lib/pubkey/gost_3410/gost_3410.cpp
+++ b/src/lib/pubkey/gost_3410/gost_3410.cpp
@@ -101,7 +101,6 @@ class GOST_3410_Signature_Operation final : public PK_Ops::Signature_with_EMSA
                                     const std::string& emsa) :
          PK_Ops::Signature_with_EMSA(emsa),
          m_group(gost_3410.domain()),
-         m_base_point(m_group.get_base_point(), m_group.get_order()),
          m_x(gost_3410.private_value())
          {}
 
@@ -112,8 +111,8 @@ class GOST_3410_Signature_Operation final : public PK_Ops::Signature_with_EMSA
 
    private:
       const EC_Group m_group;
-      Blinded_Point_Multiply m_base_point;
       const BigInt& m_x;
+      std::vector<BigInt> m_ws;
    };
 
 secure_vector<uint8_t>
@@ -133,7 +132,7 @@ GOST_3410_Signature_Operation::raw_sign(const uint8_t msg[], size_t msg_len,
    if(e == 0)
       e = 1;
 
-   const PointGFp k_times_P = m_base_point.blinded_multiply(k, rng);
+   const PointGFp k_times_P = m_group.blinded_base_point_multiply(k, rng, m_ws);
    BOTAN_ASSERT(k_times_P.on_the_curve(), "GOST 34.10 k*g is on the curve");
 
    const BigInt r = m_group.mod_order(k_times_P.get_affine_x());
diff --git a/src/lib/pubkey/sm2/sm2.cpp b/src/lib/pubkey/sm2/sm2.cpp
index 2af888bbc..9ef30d9bf 100644
--- a/src/lib/pubkey/sm2/sm2.cpp
+++ b/src/lib/pubkey/sm2/sm2.cpp
@@ -83,7 +83,6 @@ class SM2_Signature_Operation final : public PK_Ops::Signature
                               const std::string& ident,
                               const std::string& hash) :
          m_group(sm2.domain()),
-         m_base_point(m_group.get_base_point(), m_group.get_order()),
          m_x(sm2.private_value()),
          m_da_inv(sm2.get_da_inv()),
          m_hash(HashFunction::create_or_throw(hash))
@@ -102,12 +101,12 @@ class SM2_Signature_Operation final : public PK_Ops::Signature
 
    private:
       const EC_Group m_group;
-      Blinded_Point_Multiply m_base_point;
       const BigInt& m_x;
       const BigInt& m_da_inv;
 
       std::vector<uint8_t> m_za;
       std::unique_ptr<HashFunction> m_hash;
+      std::vector<BigInt> m_ws;
    };
 
 secure_vector<uint8_t>
@@ -115,7 +114,7 @@ SM2_Signature_Operation::sign(RandomNumberGenerator& rng)
    {
    const BigInt k = BigInt::random_integer(rng, 1, m_group.get_order());
 
-   const PointGFp k_times_P = m_base_point.blinded_multiply(k, rng);
+   const PointGFp k_times_P = m_group.blinded_base_point_multiply(k, rng, m_ws);
 
    const BigInt e = BigInt::decode(m_hash->final());
    const BigInt r = m_group.mod_order(k_times_P.get_affine_x() + e);
diff --git a/src/lib/pubkey/sm2/sm2_enc.cpp b/src/lib/pubkey/sm2/sm2_enc.cpp
index 462c4b968..b0d773635 100644
--- a/src/lib/pubkey/sm2/sm2_enc.cpp
+++ b/src/lib/pubkey/sm2/sm2_enc.cpp
@@ -47,8 +47,7 @@ class SM2_Encryption_Operation final : public PK_Ops::Encryption
    public:
       SM2_Encryption_Operation(const SM2_Encryption_PublicKey& key, const std::string& kdf_hash) :
          m_group(key.domain()),
-         m_base_point(m_group.get_base_point(), m_group.get_order()),
-         m_public_point(key.public_point(), m_group.get_order()),
+         m_mul_public_point(key.public_point()),
          m_kdf_hash(kdf_hash)
          {}
 
@@ -69,7 +68,7 @@ class SM2_Encryption_Operation final : public PK_Ops::Encryption
 
          const BigInt k = BigInt::random_integer(rng, 1, m_group.get_order());
 
-         const PointGFp C1 = m_base_point.blinded_multiply(k, rng);
+         const PointGFp C1 = m_group.blinded_base_point_multiply(k, rng, m_ws);
          const BigInt x1 = C1.get_affine_x();
          const BigInt y1 = C1.get_affine_y();
          std::vector<uint8_t> x1_bytes(p_bytes);
@@ -77,7 +76,7 @@ class SM2_Encryption_Operation final : public PK_Ops::Encryption
          BigInt::encode_1363(x1_bytes.data(), x1_bytes.size(), x1);
          BigInt::encode_1363(y1_bytes.data(), y1_bytes.size(), y1);
 
-         const PointGFp kPB = m_public_point.blinded_multiply(k, rng);
+         const PointGFp kPB = m_mul_public_point.mul(k, m_group.get_order(), rng, m_ws);
 
          const BigInt x2 = kPB.get_affine_x();
          const BigInt y2 = kPB.get_affine_y();
@@ -114,9 +113,9 @@ class SM2_Encryption_Operation final : public PK_Ops::Encryption
 
    private:
       const EC_Group m_group;
-      Blinded_Point_Multiply m_base_point;
-      Blinded_Point_Multiply m_public_point;
+      PointGFp_Blinded_Multiplier m_mul_public_point;
       const std::string m_kdf_hash;
+      std::vector<BigInt> m_ws;
    };
 
 class SM2_Decryption_Operation final : public PK_Ops::Decryption
@@ -134,8 +133,9 @@ class SM2_Decryption_Operation final : public PK_Ops::Decryption
                                      const uint8_t ciphertext[],
                                      size_t ciphertext_len) override
          {
-         const BigInt& cofactor = m_key.domain().get_cofactor();
-         const size_t p_bytes = m_key.domain().get_p_bytes();
+         const EC_Group& group = m_key.domain();
+         const BigInt& cofactor = group.get_cofactor();
+         const size_t p_bytes = group.get_p_bytes();
 
          valid_mask = 0x00;
 
@@ -160,18 +160,20 @@ class SM2_Decryption_Operation final : public PK_Ops::Decryption
             .end_cons()
             .verify_end();
 
-         const PointGFp C1 = m_key.domain().point(x1, y1);
+         PointGFp C1 = group.point(x1, y1);
          if(!C1.on_the_curve())
             return secure_vector<uint8_t>();
 
-         Blinded_Point_Multiply C1_mul(C1, m_key.domain().get_order());
-
-         if(cofactor > 1 && C1_mul.blinded_multiply(cofactor, m_rng).is_zero())
+         if(cofactor > 1 && (C1 * cofactor).is_zero())
             {
             return secure_vector<uint8_t>();
             }
 
-         const PointGFp dbC1 = C1_mul.blinded_multiply(m_key.private_value(), m_rng);
+         C1.randomize_repr(m_rng);
+
+         PointGFp_Blinded_Multiplier C1_mul(C1);
+
+         const PointGFp dbC1 = C1_mul.mul(m_key.private_value(), group.get_order(), m_rng, m_ws);
 
          const BigInt x2 = dbC1.get_affine_x();
          const BigInt y2 = dbC1.get_affine_y();
@@ -205,6 +207,7 @@ class SM2_Decryption_Operation final : public PK_Ops::Decryption
       const SM2_Encryption_PrivateKey& m_key;
       RandomNumberGenerator& m_rng;
       const std::string m_kdf_hash;
+      std::vector<BigInt> m_ws;
    };
 
 }
diff --git a/src/tests/data/timing/ecc_mul.vec b/src/tests/data/timing/ecc_mul.vec
new file mode 100644
index 000000000..fc703efbf
--- /dev/null
+++ b/src/tests/data/timing/ecc_mul.vec
@@ -0,0 +1,4 @@
+# leading zeros
+01
+# no leading zeros
+FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF
diff --git a/src/tests/data/timing/inverse_mod.vec b/src/tests/data/timing/inverse_mod.vec
new file mode 100644
index 000000000..5ecd03f7b
--- /dev/null
+++ b/src/tests/data/timing/inverse_mod.vec
@@ -0,0 +1,2 @@
+02
+FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF
diff --git a/src/tests/unit_ecc.cpp b/src/tests/unit_ecc.cpp
index bd8295033..166dfcd14 100644
--- a/src/tests/unit_ecc.cpp
+++ b/src/tests/unit_ecc.cpp
@@ -130,13 +130,16 @@ std::vector<Test::Result> ECC_Randomized_Tests::run()
       result.test_eq("infinite order correct", inf.is_zero(), true);
       result.test_eq("infinity on the curve", inf.on_the_curve(), true);
 
+      std::vector<Botan::BigInt> blind_ws;
+
       try
          {
          const size_t trials = (Test::run_long_tests() ? 10 : 3);
          for(size_t i = 0; i < trials; ++i)
             {
-            const size_t h = 1 + (Test::rng().next_byte() % 8);
-            Botan::Blinded_Point_Multiply blind(base_point, group_order, h);
+            const size_t w = 1 + (Test::rng().next_byte() % 8);
+
+            Botan::PointGFp_Blinded_Multiplier blinded(base_point, w);
 
             const Botan::BigInt a = Botan::BigInt::random_integer(Test::rng(), 2, group_order);
             const Botan::BigInt b = Botan::BigInt::random_integer(Test::rng(), 2, group_order);
@@ -146,9 +149,9 @@ std::vector<Test::Result> ECC_Randomized_Tests::run()
             const Botan::PointGFp Q = base_point * b;
             const Botan::PointGFp R = base_point * c;
 
-            const Botan::PointGFp P1 = blind.blinded_multiply(a, Test::rng());
-            const Botan::PointGFp Q1 = blind.blinded_multiply(b, Test::rng());
-            const Botan::PointGFp R1 = blind.blinded_multiply(c, Test::rng());
+            const Botan::PointGFp P1 = blinded.mul(a, group_order, Test::rng(), blind_ws);
+            const Botan::PointGFp Q1 = blinded.mul(b, group_order, Test::rng(), blind_ws);
+            const Botan::PointGFp R1 = blinded.mul(c, group_order, Test::rng(), blind_ws);
 
             const Botan::PointGFp A1 = P + Q;
             const Botan::PointGFp A2 = Q + P;
@@ -278,6 +281,13 @@ Test::Result test_groups()
       result.confirm("EC_Group is known", !group.get_curve_oid().empty());
       result.test_eq("EC_Group has correct bit size", group.get_p().bits(), group.get_p_bits());
       result.test_eq("EC_Group has byte size", group.get_p().bytes(), group.get_p_bytes());
+
+      bool a_is_minus_3 = group.a_is_minus_3();
+
+      if(a_is_minus_3)
+         result.test_eq("Group A equals -3", group.get_a(), group.get_p() - 3);
+      else
+         result.test_ne("Group " + group_name + " A does not equal -3", group.get_a(), group.get_p() - 3);
       }
    return result;
    }