ipsilon
diff --git a/‎lib/evmone_precompiles/ecc.hpp‎
Lines changed: 101 additions & 21 deletions b/‎lib/evmone_precompiles/ecc.hpp‎
Lines changed: 101 additions & 21 deletions
diff --git a/‎test/unittests/CMakeLists.txt‎
Lines changed: 1 addition & 0 deletions b/‎test/unittests/CMakeLists.txt‎
Lines changed: 1 addition & 0 deletions
diff --git a/‎test/unittests/crypto_wnaf.cpp‎
Lines changed: 107 additions & 0 deletions b/‎test/unittests/crypto_wnaf.cpp‎
Lines changed: 107 additions & 0 deletions
@@ -492,40 +492,120 @@ ProjPoint<Curve> mul(const AffinePoint<Curve>& p, typename Curve::uint_type c) n
 
 /// Computes multi-scalar multiplication of u×P ⊕ v×Q.
 ///
-/// The implementation uses the "Straus-Shamir trick": https://eprint.iacr.org/2003/257.pdf#page=7.
-template <typename Curve>
-ProjPoint<Curve> msm(const typename Curve::uint_type& u, const AffinePoint<Curve>& p,
-    const typename Curve::uint_type& v, const AffinePoint<Curve>& q)
+/// The implementation uses windowed NAF (wNAF) recoding.
+template <typename UIntT>
+class NAF
 {
-    ProjPoint<Curve> r;
+public:
+    using digit_type = int8_t;
 
-    const auto w = u | v;
-    const auto bit_width = sizeof(w) * 8 - intx::clz(w);
-    if (bit_width == 0)
-        return r;
+private:
+    /// The storage for the NAF digits, starting from the least significant one.
+    /// For a k-bit scalar, there can be at most k+1 digits.
+    std::array<digit_type, sizeof(UIntT) * 8 + 1> digits_{};
 
-    // Precompute affine P + Q. Works correctly if P == Q.
-    const auto h = add_affine(p, q);
+    /// The number of non-zero digits used in the representation.
+    size_t width_ = 0;
 
-    // Create lookup table for points. The index 0 is unused.
-    // TODO: Put 0 at index 0 and use it in the loop to avoid the branch.
-    const AffinePoint<Curve>* const points[]{nullptr, &p, &q, &h};
+public:
+    size_t width() const noexcept { return width_; }
 
-    for (auto i = bit_width; i != 0; --i)
+    digit_type operator[](size_t i) const noexcept { return digits_[i]; }
+
+    void set(size_t i, digit_type d) noexcept
+    {
+        assert(d != 0);
+        assert(i >= width_);
+        digits_[i] = d;
+        width_ = i + 1;
+    }
+};
+
+/// Convert an unsigned scalar value to its windowed Non-adjacent Form (wNAF).
+///
+/// See
+/// https://en.wikipedia.org/wiki/Elliptic_curve_point_multiplication#w-ary_non-adjacent_form_(wNAF)_method.
+template <unsigned W, typename UIntT>
+constexpr NAF<UIntT> to_wnaf(UIntT k) noexcept
+{
+    using digit_type = NAF<UIntT>::digit_type;
+    static_assert(W >= 2);
+    static_assert(W <= sizeof(digit_type) * 8);
+    constexpr unsigned RADIX = 1 << W;
+
+    NAF<UIntT> naf;
+    for (size_t i = 0; k != 0; ++i, k >>= 1)
+    {
+        const auto r = static_cast<unsigned>(k) % RADIX;
+        if (r % 2 != 0)
+        {
+            const auto d_sign = r > RADIX / 2;
+            const auto d_abs = d_sign ? RADIX - r : r;
+            const auto d = d_sign ? -d_abs : d_abs;
+            naf.set(i, static_cast<digit_type>(d));
+            k -= d_sign ? -UIntT{d_abs} : UIntT{d_abs};  // intx lacks sign extending conversion.
+        }
+    }
+    return naf;
+}
+
+template <unsigned W, typename Curve>
+void precompute_wnaf_table(
+    std::span<ProjPoint<Curve>, 1 << (W - 2)> table, const AffinePoint<Curve>& p) noexcept
+{
+    table[0] = ProjPoint{p};           // 1P.
+    const auto two_p = dbl(table[0]);  // 2P.
+
+    for (size_t i = 1; i < table.size(); ++i)
+        table[i] = add(table[i - 1], two_p);  // (2i+3)P = (2i+1)P + 2P.
+}
+
+template <unsigned W, size_t S, typename Curve>
+ProjPoint<Curve> msm_wnaf(std::span<const AffinePoint<Curve>* const, S> points,
+    std::span<const typename Curve::uint_type* const, S> scalars) noexcept
+{
+    static constexpr size_t TABLE_SIZE = 1 << (W - 2);
+
+    std::array<NAF<typename Curve::uint_type>, S> nafs;
+    std::array<ProjPoint<Curve>, S * TABLE_SIZE> joint_table;
+
+    for (size_t s = 0; s < S; ++s)
+    {
+        nafs[s] = to_wnaf<W>(*scalars[s]);
+        precompute_wnaf_table<W>(
+            std::span<ProjPoint<Curve>, TABLE_SIZE>{&joint_table[s * TABLE_SIZE], TABLE_SIZE},
+            *points[s]);
+    }
+
+    ProjPoint<Curve> r;
+    const auto max_width =
+        std::ranges::max(nafs, {}, &NAF<typename Curve::uint_type>::width).width();
+    for (size_t i = max_width; i != 0; --i)
     {
         r = dbl(r);
 
-        const auto u_bit = bit_test(u, i - 1);
-        const auto v_bit = bit_test(v, i - 1);
-        const auto idx = 2 * size_t{v_bit} + size_t{u_bit};
-        if (idx == 0)
-            continue;
-        r = add(r, *points[idx]);
+        for (size_t s = 0; s < S; ++s)
+        {
+            const auto d = nafs[s][i - 1];
+            if (d == 0)  // TODO: likely
+                continue;
+
+            const auto* table = &joint_table[s * TABLE_SIZE];
+            const auto& pt = table[(static_cast<unsigned>(std::abs(d)) - 1) / 2];
+            r = add(r, d >= 0 ? pt : -pt);
+        }
     }
 
     return r;
 }
 
+template <typename Curve>
+ProjPoint<Curve> msm(const typename Curve::uint_type& u, const AffinePoint<Curve>& p,
+    const typename Curve::uint_type& v, const AffinePoint<Curve>& q)
+{
+    return msm_wnaf<4, 2, Curve>(std::array{&p, &q}, std::array{&u, &v});
+}
+
 template <typename UIntT>
 struct SignedScalar
 {
 
@@ -10,6 +10,7 @@ target_sources(
     baseline_analysis_test.cpp
     blockchaintest_loader_test.cpp
     bytecode_test.cpp
+    crypto_wnaf.cpp
     evm_fixture.cpp
     evm_fixture.hpp
     evm_test.cpp
 
@@ -0,0 +1,107 @@
+// evmone: Fast Ethereum Virtual Machine implementation
+// Copyright 2023 The evmone Authors.
+// SPDX-License-Identifier: Apache-2.0
+
+#include <evmone_precompiles/ecc.hpp>
+#include <gtest/gtest.h>
+#include <intx/intx.hpp>
+#include <random>
+
+using namespace evmmax::ecc;
+
+namespace
+{
+template <typename UIntT>
+UIntT evaluate(NAF<UIntT> naf)
+{
+    UIntT result = 0;
+    UIntT base = 1;
+    for (size_t i = 0; i < naf.width(); ++i)
+    {
+        const auto d = naf[i];
+        const auto d_abs = static_cast<unsigned>(std::abs(d));
+        const auto d_sign = d < 0;
+        const auto r_abs = UIntT{d_abs};
+        const auto r = d_sign ? -r_abs : r_abs;
+        result += r * base;
+        base <<= 1;
+    }
+
+    if (naf.width() == 0)
+    {
+        // NAF == 0 <=> result == 0.
+        EXPECT_EQ(result, 0);
+    }
+    else
+    {
+        // The most significant digit must be non-zero.
+        EXPECT_NE(naf[naf.width() - 1], 0);
+    }
+    return result;
+}
+}  // namespace
+
+TEST(crypto_wnaf, example1)
+{
+    const auto naf = to_wnaf<3>(uint32_t{21});
+    EXPECT_EQ(naf.width(), 4u);
+    EXPECT_EQ(naf[0], -3);
+    EXPECT_EQ(naf[1], 0);
+    EXPECT_EQ(naf[2], 0);
+    EXPECT_EQ(naf[3], 3);
+    EXPECT_EQ(naf[4], 0);
+    EXPECT_EQ(evaluate(naf), 21u);
+}
+
+TEST(crypto_wnaf, max_width)
+{
+    const auto x = uint32_t{0xfffffffe};
+    const auto naf = to_wnaf<4>(x);
+    EXPECT_EQ(naf.width(), 33u);
+    EXPECT_EQ(naf[0], 0);
+    EXPECT_EQ(naf[1], -1);
+    for (size_t i = 2; i <= 31; ++i)
+        EXPECT_EQ(naf[i], 0);
+    EXPECT_EQ(naf[32], 1);
+    EXPECT_EQ(evaluate(naf), x);
+}
+
+TEST(crypto_wnaf, max_digit)
+{
+    const auto x = uint32_t{0xfffffcfe};
+    const auto naf = to_wnaf<8>(x);
+    EXPECT_EQ(naf.width(), 33u);
+    EXPECT_EQ(naf[1], 127);
+    EXPECT_EQ(evaluate(naf), x);
+}
+
+TEST(crypto_wnaf, min_digit)
+{
+    const auto x = uint32_t{0x102};
+    const auto naf = to_wnaf<8>(x);
+    EXPECT_EQ(naf.width(), 10u);
+    EXPECT_EQ(naf[1], -127);
+    EXPECT_EQ(evaluate(naf), x);
+}
+
+TEST(crypto_wnaf, uint256_fuzz)
+{
+    std::mt19937_64 rng{std::random_device{}()};
+    std::uniform_int_distribution<uint64_t> dist{};
+    const intx::uint256 start{dist(rng), dist(rng), dist(rng), dist(rng)};
+
+    for (size_t i = 0; i < 100; ++i)
+    {
+        const auto x = start + i;
+        const auto naf2 = to_wnaf<2>(x);
+        ASSERT_EQ(evaluate(naf2), x);
+        const auto naf3 = to_wnaf<3>(x);
+        ASSERT_EQ(evaluate(naf3), x);
+        const auto naf4 = to_wnaf<4>(x);
+        ASSERT_EQ(evaluate(naf4), x);
+        const auto naf5 = to_wnaf<5>(x);
+        ASSERT_EQ(evaluate(naf5), x);
+        const auto naf8 = to_wnaf<8>(x);
+        ASSERT_EQ(evaluate(naf8), x);
+    }
+}