Optimize EC point doubling for a == 0 and a == -3

3 files changed, 61 insertions, 9 deletions

diff --git a/src/lib/pubkey/ec_group/curve_gfp.cpp b/src/lib/pubkey/ec_group/curve_gfp.cpp
index 8953edba4..9f614ac61 100644
--- a/src/lib/pubkey/ec_group/curve_gfp.cpp
+++ b/src/lib/pubkey/ec_group/curve_gfp.cpp
@@ -34,8 +34,14 @@ class CurveGFp_Montgomery final : public CurveGFp_Repr
          m_r3  = mod_p.multiply(m_r, m_r2);
          m_a_r = mod_p.multiply(m_r, m_a);
          m_b_r = mod_p.multiply(m_r, m_b);
+
+         m_a_is_zero = m_a.is_zero();
+         m_a_is_minus_3 = (m_a + 3 == m_p);
          }
 
+      bool a_is_zero() const override { return m_a_is_zero; }
+      bool a_is_minus_3() const override { return m_a_is_minus_3; }
+
       const BigInt& get_a() const override { return m_a; }
 
       const BigInt& get_b() const override { return m_b; }
@@ -78,6 +84,9 @@ class CurveGFp_Montgomery final : public CurveGFp_Repr
       // Montgomery parameters
       BigInt m_r, m_r2, m_r3;
       word m_p_dash;
+
+      bool m_a_is_zero;
+      bool m_a_is_minus_3;
    };
 
 BigInt CurveGFp_Montgomery::invert_element(const BigInt& x, secure_vector<word>& ws) const
@@ -190,8 +199,12 @@ class CurveGFp_NIST : public CurveGFp_Repr
       CurveGFp_NIST(size_t p_bits, const BigInt& a, const BigInt& b) :
          m_a(a), m_b(b), m_p_words((p_bits + BOTAN_MP_WORD_BITS - 1) / BOTAN_MP_WORD_BITS)
          {
+         // All Solinas prime curves are assumed a == -3
          }
 
+      bool a_is_zero() const override { return false; }
+      bool a_is_minus_3() const override { return true; }
+
       const BigInt& get_a() const override { return m_a; }
 
       const BigInt& get_b() const override { return m_b; }
@@ -233,7 +246,6 @@ class CurveGFp_NIST : public CurveGFp_Repr
       size_t m_p_words; // cache of m_p.sig_words()
    };
 
-
 BigInt CurveGFp_NIST::invert_element(const BigInt& x, secure_vector<word>& ws) const
    {
    BOTAN_UNUSED(ws);
diff --git a/src/lib/pubkey/ec_group/curve_gfp.h b/src/lib/pubkey/ec_group/curve_gfp.h
index 8a83a9c41..3a17c6678 100644
--- a/src/lib/pubkey/ec_group/curve_gfp.h
+++ b/src/lib/pubkey/ec_group/curve_gfp.h
@@ -30,6 +30,10 @@ class BOTAN_UNSTABLE_API CurveGFp_Repr
 
       virtual bool is_one(const BigInt& x) const = 0;
 
+      virtual bool a_is_zero() const = 0;
+
+      virtual bool a_is_minus_3() const = 0;
+
       /*
       * Returns to_curve_rep(get_a())
       */
@@ -116,6 +120,9 @@ class BOTAN_UNSTABLE_API CurveGFp final
 
       const BigInt& get_b_rep() const { return m_repr->get_b_rep(); }
 
+      bool a_is_minus_3() const { return m_repr->a_is_minus_3(); }
+      bool a_is_zero() const { return m_repr->a_is_zero(); }
+
       bool is_one(const BigInt& x) const { return m_repr->is_one(x); }
 
       BigInt invert_element(const BigInt& x, secure_vector<word>& ws) const
diff --git a/src/lib/pubkey/ec_group/point_gfp.cpp b/src/lib/pubkey/ec_group/point_gfp.cpp
index acbda95d4..860459252 100644
--- a/src/lib/pubkey/ec_group/point_gfp.cpp
+++ b/src/lib/pubkey/ec_group/point_gfp.cpp
@@ -316,14 +316,47 @@ void PointGFp::mult2(std::vector<BigInt>& ws_bn)
    T1 <<= 2; // * 4
    T1.reduce_below(p, T3.get_word_vector());
 
-   m_curve.sqr(T3, m_coord_z, ws); // z^2
-   m_curve.sqr(T4, T3, ws); // z^4
-   m_curve.mul(T3, m_curve.get_a_rep(), T4, ws);
-
-   m_curve.sqr(T4, m_coord_x, ws);
-   T4 *= 3;
-   T4 += T3;
-   T4.reduce_below(p, T3.get_word_vector());
+   if(m_curve.a_is_zero())
+      {
+      // if a == 0 then 3*x^2 + a*z^4 is just 3*x^2
+      m_curve.sqr(T4, m_coord_x, ws); // x^2
+      T4 *= 3; // 3*x^2
+      T4.reduce_below(p, T3.get_word_vector());
+      }
+   else if(m_curve.a_is_minus_3())
+      {
+      /*
+      if a == -3 then
+        3*x^2 + a*z^4 == 3*x^2 - 3*z^4 == 3*(x^2-z^4) == 3*(x-z^2)*(x+z^2)
+      */
+      m_curve.sqr(T3, m_coord_z, ws); // z^2
+
+      // (x-z^2)
+      T2 = m_coord_x;
+      T2 -= T3;
+      if(T2.is_negative())
+         T2 += p;
+
+      // (x+z^2)
+      T3 += m_coord_x;
+      T3.reduce_below(p, T4.get_word_vector());
+
+      m_curve.mul(T4, T2, T3, ws); // (x-z^2)*(x+z^2)
+
+      T4 *= 3; // 3*(x-z^2)*(x+z^2)
+      T4.reduce_below(p, T3.get_word_vector());
+      }
+   else
+      {
+      m_curve.sqr(T3, m_coord_z, ws); // z^2
+      m_curve.sqr(T4, T3, ws); // z^4
+      m_curve.mul(T3, m_curve.get_a_rep(), T4, ws); // a*z^4
+
+      m_curve.sqr(T4, m_coord_x, ws); // x^2
+      T4 *= 3; // 3*x^2
+      T4 += T3; // 3*x^2 + a*z^4
+      T4.reduce_below(p, T3.get_word_vector());
+      }
 
    m_curve.sqr(T2, T4, ws);
    T2 -= T1;