Add fast iterated point doubling

Improves ECDH performance for curves with A == -3 by 5-9%
randombit · Jul 17, 2024 · b5a922f · b5a922f
1 parent 5477167
commit b5a922f
Show file tree

Hide file tree

Showing 5 changed files with 60 additions and 14 deletions.
diff --git a/src/lib/math/numbertheory/nistp_redc.cpp b/src/lib/math/numbertheory/nistp_redc.cpp
@@ -75,7 +75,7 @@ namespace {
 * Treating this MPI as a sequence of 32-bit words in big-endian
 * order, return word i. The array is assumed to be large enough.
 */
-inline uint32_t get_uint32(const word xw[], size_t i) {
+constexpr uint32_t get_uint32(const word xw[], size_t i) {
 #if(BOTAN_MP_WORD_BITS == 32)
    return xw[i];
 #else

diff --git a/src/lib/math/pcurves/pcurves_impl/pcurves_impl.h b/src/lib/math/pcurves/pcurves_impl/pcurves_impl.h
@@ -176,7 +176,7 @@ class IntMod final {
          return Self(Rep::redc(z));
       }
 
-      void square_n(size_t n) {
+      constexpr void square_n(size_t n) {
          std::array<W, 2 * N> z;
          for(size_t i = 0; i != n; ++i) {
             comba_sqr<N>(z.data(), this->data());
@@ -698,14 +698,52 @@ class ProjectiveCurvePoint {
       }
 
       constexpr Self dbl_n(size_t n) const {
-         // TODO it is possible to optimize this by carrying over values from
-         // the previous iteration into the next
+         if constexpr(Self::A_is_minus_3) {
+            /*
+            Repeated doubling using an adaptation of Algorithm 3.23 in
+            "Guide To Elliptic Curve Cryptography"
+            Hankerson, Menezes, Vanstone
 
-         Self pt = (*this);
-         for(size_t i = 0; i != n; ++i) {
-            pt = pt.dbl();
+            For A == -3
+            Cost: 2S + 1*2 + n*(4S + 4M + 2*2 + 1*3 + 4A)
+            Naive doubling
+            Cost: n*(4S + 4M + 2*2 + 1*3 + 5A + 1*4 + 1*8)
+
+            TODO adapt this for A == 0 and/or generic A cases
+            */
+
+            // The inverse of 2 modulo P is (P/2)+1; this avoids a constexpr time
+            // general inversion, which some compilers can't handle
+            constexpr auto INV_2 = FieldElement::from_words(p_div_2_plus_1(Params::PW));
+
+            auto nx = x();
+            auto ny = y();
+            auto nz = z();
+            ny = ny.mul2();
+            auto w = nz.square().square();
+
+            while(n > 0) {
+               const auto ny2 = ny.square();
+               const auto ny4 = ny2.square();
+               const auto t1 = (nx.square() - w).mul3();
+               const auto t2 = nx * ny2;
+               nx = t1.square() - t2.mul2();
+               nz *= ny;
+               ny = t1 * (t2 - nx).mul2() - ny4;
+               n--;
+               if(n > 0) {
+                  w *= ny4;
+               }
+            }
+            ny *= INV_2;
+            return Self(nx, ny, nz);
+         } else {
+            Self pt = (*this);
+            for(size_t i = 0; i != n; ++i) {
+               pt = pt.dbl();
+            }
+            return pt;
          }
-         return pt;
       }
 
       constexpr Self dbl() const {
@@ -718,7 +756,7 @@ class ProjectiveCurvePoint {
             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)
 
-            Cost: 2M + 2A + 1*3
+            Cost: 1M + 1S + 2A + 1*3
             */
             const auto z2 = z().square();
             m = (x() - z2).mul3() * (x() + z2);
@@ -731,7 +769,6 @@ class ProjectiveCurvePoint {
             const auto z2 = z().square();
             m = x().square().mul3() + A * z2.square();
          }
-
          const auto y2 = y().square();
          const auto s = x().mul4() * y2;
          const auto nx = m.square() - s.mul2();

diff --git a/src/lib/math/pcurves/pcurves_impl/pcurves_solinas.h b/src/lib/math/pcurves/pcurves_impl/pcurves_solinas.h
@@ -47,9 +47,9 @@ class SolinasAccum {
 
       static constexpr size_t N32 = N * (WordInfo<W>::bits / 32);
 
-      SolinasAccum(std::array<W, N>& r) : m_r(r), m_S(0), m_idx(0) {}
+      constexpr SolinasAccum(std::array<W, N>& r) : m_r(r), m_S(0), m_idx(0) {}
 
-      void accum(int64_t v) {
+      constexpr void accum(int64_t v) {
          BOTAN_DEBUG_ASSERT(m_idx < N32);
 
          m_S += v;
@@ -65,7 +65,7 @@ class SolinasAccum {
          m_idx += 1;
       }
 
-      W final_carry(int64_t C) {
+      constexpr W final_carry(int64_t C) {
          BOTAN_DEBUG_ASSERT(m_idx == N32);
          m_S += C;
          BOTAN_DEBUG_ASSERT(m_S >= 0);

diff --git a/src/lib/math/pcurves/pcurves_impl/pcurves_util.h b/src/lib/math/pcurves/pcurves_impl/pcurves_util.h
@@ -181,6 +181,15 @@ inline consteval std::array<W, N> p_minus_1_over_2(const std::array<W, N>& p) {
    return r;
 }
 
+template <WordType W, size_t N>
+inline consteval std::array<W, N> p_div_2_plus_1(const std::array<W, N>& p) {
+   const W one = 1;
+   std::array<W, N> r = p;
+   shift_right<1>(r);
+   bigint_add2_nc(r.data(), N, &one, 1);
+   return r;
+}
+
 template <WordType W, size_t N>
 inline consteval size_t count_bits(const std::array<W, N>& p) {
    auto get_bit = [&](size_t i) {

diff --git a/src/tests/test_pcurves.cpp b/src/tests/test_pcurves.cpp
@@ -44,7 +44,7 @@ class Pcurve_Basemul_Tests final : public Text_Based_Test {
                result.test_eq("mul correct", pt4, P_bytes);
 
                auto pt5 = curve->mul(g, scalar.value(), null_rng).to_affine().serialize();
-               result.test_eq("mul correct", pt5, P_bytes);
+               result.test_eq("mul correct (Null_RNG)", pt5, P_bytes);
             } else {
                result.test_failure("Curve rejected scalar input");
             }