Add specialized reduction for P-521 along with 9x9 Comba routines.

Roughly 35-50% faster on my laptop (depending on if mlock is enabled,
the overhead in that allocator is becoming much more of a hotspot).

6 files changed, 229 insertions, 3 deletions

diff --git a/src/lib/math/ec_gfp/curve_gfp.cpp b/src/lib/math/ec_gfp/curve_gfp.cpp
index 2fa9f391e..2ee4ae41e 100644
--- a/src/lib/math/ec_gfp/curve_gfp.cpp
+++ b/src/lib/math/ec_gfp/curve_gfp.cpp
@@ -6,6 +6,7 @@
 */
 
 #include <botan/curve_gfp.h>
+#include <botan/internal/curve_nistp.h>
 #include <botan/internal/mp_core.h>
 
 namespace Botan {
@@ -37,6 +38,8 @@ class CurveGFp_Montgomery : public CurveGFp_Repr
 
       const BigInt& get_b_rep() const override { return m_b_r; }
 
+      size_t get_p_words() const override { return m_p_words; }
+
       void to_curve_rep(BigInt& x, secure_vector<word>& ws) const override;
 
       void from_curve_rep(BigInt& x, secure_vector<word>& ws) const override;
@@ -113,9 +116,38 @@ void CurveGFp_Montgomery::curve_sqr(BigInt& z, const BigInt& x,
 
 }
 
+// Default implementation
+void CurveGFp_Repr::normalize(BigInt& x, secure_vector<word>& ws, size_t /*bound*/) const
+   {
+   const BigInt& p = get_p();
+
+   while(x.is_negative())
+      x += p;
+
+   const size_t p_words = get_p_words();
+   const word* prime = p.data();
+
+   x.grow_to(p_words + 1);
+
+   if(ws.size() < p_words + 1)
+      ws.resize(p_words + 1);
+
+   //FIXME: take into account bound if > 0
+   while(true)
+      {
+      if(bigint_sub3(&ws[0], x.data(), p_words, prime, p_words)) // borrow?
+         break;
+
+      x.swap_reg(ws);
+      }
+   }
+
 std::shared_ptr<CurveGFp_Repr>
 CurveGFp::choose_repr(const BigInt& p, const BigInt& a, const BigInt& b)
    {
+   if(p == CurveGFp_P521::prime())
+      return std::shared_ptr<CurveGFp_Repr>(new CurveGFp_P521(a, b));
+
    return std::shared_ptr<CurveGFp_Repr>(new CurveGFp_Montgomery(p, a, b));
    }
 
diff --git a/src/lib/math/ec_gfp/curve_gfp.h b/src/lib/math/ec_gfp/curve_gfp.h
index 67fb3d2cf..59639d537 100644
--- a/src/lib/math/ec_gfp/curve_gfp.h
+++ b/src/lib/math/ec_gfp/curve_gfp.h
@@ -24,6 +24,8 @@ class CurveGFp_Repr
       virtual const BigInt& get_a() const = 0;
       virtual const BigInt& get_b() const = 0;
 
+      virtual size_t get_p_words() const = 0;
+
       /*
       * Returns to_curve_rep(get_a())
       */
@@ -43,6 +45,10 @@ class CurveGFp_Repr
 
       virtual void curve_sqr(BigInt& z, const BigInt& x,
                              secure_vector<word>& ws) const = 0;
+
+      virtual void normalize(BigInt& x,
+                             secure_vector<word>& ws,
+                             size_t bound) const;
    };
 
 /**
@@ -109,6 +115,8 @@ class BOTAN_DLL CurveGFp
          return xt;
          }
 
+      // TODO: from_rep taking && ref
+
       void mul(BigInt& z, const BigInt& x, const BigInt& y, secure_vector<word>& ws) const
          {
          m_repr->curve_mul(z, x, y, ws);
@@ -133,6 +141,16 @@ class BOTAN_DLL CurveGFp
          return z;
          }
 
+      /**
+      * Adjust x to be in [0,p)
+      * @param bound if greater than zero, assume that no more than bound
+      *        additions or subtractions are required to move x into range.
+      */
+      void normalize(BigInt& x, secure_vector<word>& ws, size_t bound = 0) const
+         {
+         m_repr->normalize(x, ws, bound);
+         }
+
       void swap(CurveGFp& other)
          {
          std::swap(m_repr, other.m_repr);
diff --git a/src/lib/math/ec_gfp/curve_nistp.cpp b/src/lib/math/ec_gfp/curve_nistp.cpp
new file mode 100644
index 000000000..1aac01f56
--- /dev/null
+++ b/src/lib/math/ec_gfp/curve_nistp.cpp
@@ -0,0 +1,95 @@
+/*
+* NIST curve reduction
+* (C) 2014 Jack Lloyd
+*
+* Distributed under the terms of the Botan license
+*/
+
+#include <botan/internal/curve_nistp.h>
+#include <botan/internal/mp_core.h>
+#include <botan/internal/rounding.h>
+#include <botan/hex.h>
+
+namespace Botan {
+
+void CurveGFp_NIST::curve_mul(BigInt& z, const BigInt& x, const BigInt& y,
+                              secure_vector<word>& ws) const
+   {
+   if(x.is_zero() || y.is_zero())
+      {
+      z = 0;
+      return;
+      }
+
+   const size_t p_words = get_p_words();
+   const size_t output_size = 2*p_words + 1;
+   ws.resize(2*(p_words+2));
+
+   z.grow_to(output_size);
+   z.clear();
+
+   bigint_mul(z.mutable_data(), output_size, &ws[0],
+              x.data(), x.size(), x.sig_words(),
+              y.data(), y.size(), y.sig_words());
+
+   this->redc(z, ws);
+   }
+
+void CurveGFp_NIST::curve_sqr(BigInt& z, const BigInt& x,
+                              secure_vector<word>& ws) const
+   {
+   if(x.is_zero())
+      {
+      z = 0;
+      return;
+      }
+
+   const size_t p_words = get_p_words();
+   const size_t output_size = 2*p_words + 1;
+
+   ws.resize(2*(p_words+2));
+
+   z.grow_to(output_size);
+   z.clear();
+
+   bigint_sqr(z.mutable_data(), output_size, &ws[0],
+              x.data(), x.size(), x.sig_words());
+
+   this->redc(z, ws);
+   }
+
+//static
+const BigInt& CurveGFp_P521::prime()
+   {
+   static const BigInt p521("0x1FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF"
+                               "FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF");
+
+   return p521;
+   }
+
+void CurveGFp_P521::redc(BigInt& x, secure_vector<word>& ws) const
+   {
+   const size_t p_words = get_p_words();
+
+   const size_t shift_words = 521 / MP_WORD_BITS,
+                shift_bits  = 521 % MP_WORD_BITS;
+
+   const size_t x_sw = x.sig_words();
+
+   if(x_sw < p_words)
+      return; // already smaller
+
+   if(ws.size() < p_words + 1)
+      ws.resize(p_words + 1);
+
+   clear_mem(&ws[0], ws.size());
+   bigint_shr2(&ws[0], x.data(), x_sw, shift_words, shift_bits);
+
+   x.mask_bits(521);
+
+   bigint_add3(x.mutable_data(), x.data(), p_words, &ws[0], p_words);
+
+   normalize(x, ws, max_redc_subtractions());
+   }
+
+}
diff --git a/src/lib/math/ec_gfp/curve_nistp.h b/src/lib/math/ec_gfp/curve_nistp.h
new file mode 100644
index 000000000..ffa32d377
--- /dev/null
+++ b/src/lib/math/ec_gfp/curve_nistp.h
@@ -0,0 +1,75 @@
+/*
+* NIST elliptic curves over GF(p)
+* (C) 2014 Jack Lloyd
+*
+* Distributed under the terms of the Botan license
+*/
+
+#ifndef BOTAN_GFP_CURVE_NIST_H__
+#define BOTAN_GFP_CURVE_NIST_H__
+
+#include <botan/curve_gfp.h>
+#include <memory>
+
+namespace Botan {
+
+class CurveGFp_NIST : public CurveGFp_Repr
+   {
+   public:
+      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)
+         {
+         }
+
+      size_t get_p_words() const override { return m_p_words; }
+
+      const BigInt& get_a() const override { return m_a; }
+
+      const BigInt& get_b() const override { return m_b; }
+
+      const BigInt& get_a_rep() const override { return m_a; }
+
+      const BigInt& get_b_rep() const override { return m_b; }
+
+      void to_curve_rep(BigInt& x, secure_vector<word>& ws) const override
+         { redc(x, ws); }
+
+      void from_curve_rep(BigInt& x, secure_vector<word>& ws) const override
+         { redc(x, ws); }
+
+      void curve_mul(BigInt& z, const BigInt& x, const BigInt& y,
+                     secure_vector<word>& ws) const override;
+
+      void curve_sqr(BigInt& z, const BigInt& x,
+                     secure_vector<word>& ws) const override;
+   private:
+      virtual void redc(BigInt& x, secure_vector<word>& ws) const = 0;
+
+      virtual size_t max_redc_subtractions() const = 0;
+
+      // Curve parameters
+      BigInt m_a, m_b;
+      size_t m_p_words; // cache of m_p.sig_words()
+   };
+
+/**
+* The NIST P-521 curve
+*/
+class CurveGFp_P521 : public CurveGFp_NIST
+   {
+   public:
+      CurveGFp_P521(const BigInt& a, const BigInt& b) : CurveGFp_NIST(521, a, b) {}
+
+      static const BigInt& prime();
+
+      const BigInt& get_p() const override { return CurveGFp_P521::prime(); }
+
+   private:
+      void redc(BigInt& x, secure_vector<word>& ws) const override;
+
+      size_t max_redc_subtractions() const override { return 1; }
+   };
+
+}
+
+#endif
diff --git a/src/lib/math/ec_gfp/info.txt b/src/lib/math/ec_gfp/info.txt
index 2e67f608e..9af40c752 100644
--- a/src/lib/math/ec_gfp/info.txt
+++ b/src/lib/math/ec_gfp/info.txt
@@ -7,6 +7,10 @@ curve_gfp.h
 point_gfp.h
 </header:public>
 
+<header:internal>
+curve_nistp.h
+</header:internal>
+
 <requires>
 numbertheory
 </requires>
diff --git a/src/lib/math/ec_gfp/point_gfp.cpp b/src/lib/math/ec_gfp/point_gfp.cpp
index 3d244d0f0..6bae35e5f 100644
--- a/src/lib/math/ec_gfp/point_gfp.cpp
+++ b/src/lib/math/ec_gfp/point_gfp.cpp
@@ -479,18 +479,20 @@ BigInt decompress_point(bool yMod2,
    {
    BigInt xpow3 = x * x * x;
 
+   const BigInt& p = curve.get_p();
+
    BigInt g = curve.get_a() * x;
    g += xpow3;
    g += curve.get_b();
-   g = g % curve.get_p();
+   g = g % p;
 
-   BigInt z = ressol(g, curve.get_p());
+   BigInt z = ressol(g, p);
 
    if(z < 0)
       throw Illegal_Point("error during EC point decompression");
 
    if(z.get_bit(0) != yMod2)
-      z = curve.get_p() - z;
+      z = p - z;
 
    return z;
    }