crypto: Use fixed-window exponentiation in modexp

lib/evmone_precompiles/modexp.cpp

@@ -280,8 +280,8 @@ class Exponent
 
     [[nodiscard]] size_t bit_width() const noexcept { return bit_width_; }
 
-    /// Returns the bit value of the exponent at the given index, counting from the most significant
-    /// bit (e[0] is the top bit).
+    /// Returns the bit value of the exponent at the given index, counting from the least
+    /// significant bit (e[0] is the bottom bit, e[bit_width() - 1] is the top bit, always set).
     bool operator[](size_t index) const noexcept
     {
         // TODO: Replace this with a custom iterator type.
@@ -368,8 +368,33 @@ template <>
     mul_amm_256(r, x, y, mod, mod_inv);
 }
 
+/// Maximum window width used by the windowed method in modexp_odd.
+constexpr unsigned MAX_WINDOW_WIDTH = 4;
+
+/// Number of precomputed values for the max width windowed method.
+constexpr size_t MAX_PRECOMPUTED = (size_t{1} << MAX_WINDOW_WIDTH) - 1;
+
+/// Selects the fixed-window width from the exponent bit length.
+///
+/// TODO: Switch to a sliding window: the table then holds only odd powers. Measured ~3-7%.
+/// TODO: Tune for the densest exponent instead of the average, because gas is charged on
+///   exponent bit length and ignores Hamming weight. The width then collapses to
+///   min(MAX_WINDOW_WIDTH, (bit_width(exp_bits) + 1) / 2). Measured +10.7% worst case.
+constexpr unsigned window_width(size_t exp_bits) noexcept
+{
+    // Break-even points for a random exponent, where the 2^w extra table multiplies stop
+    // being repaid: 2^w / ((1-2^-w)/w - (1-2^-(w+1))/(w+1)) = 16, 48, 140 (rounded to 144).
+    if (exp_bits <= 16)
+        return 1;
+    if (exp_bits <= 48)
+        return 2;
+    if (exp_bits <= 144)
+        return 3;
+    return MAX_WINDOW_WIDTH;
+}
+
 /// Computes result[] = base[]^exp % mod[] for odd mod[] (mod[0] % 2 != 0).
-/// Scratch space required: 4n + 3*base.size() + 2 words, where n = mod.size().
+/// Scratch space required: (MAX_PRECOMPUTED + 3)*mod.size() + 3*base.size() + 2 words.
 void modexp_odd(std::span<uint64_t> result, std::span<const uint64_t> base, Exponent exp,
     std::span<const uint64_t> mod, std::span<uint64_t> scratch) noexcept
 {
@@ -380,17 +405,24 @@ void modexp_odd(std::span<uint64_t> result, std::span<const uint64_t> base, Expo
 
     const auto n = mod.size();
     const auto mod_inv = -evmmax::modinv(mod[0]);
+    const auto exp_bits = exp.bit_width();
 
-    // Layout: u[n+base.size()] | base_mont[n] | t/rem_scratch[max(n, 2*(n+base.size())+2)]
-    // t and rem_scratch share the same region (exclusive lifetimes).
-    assert(scratch.size() >= 4 * n + 3 * base.size() + 2);
+    const auto w = window_width(exp_bits);
+    const auto table_size = (size_t{1} << w) - 1;
 
-    // Compute base_mont = (base * R) % mod, where R = 2^(n*64).
-    // The numerator u = base << (n*64): base in the upper words, lower n words are zero.
+    // Layout: u[n + base.size()] | table[MAX_PRECOMPUTED*n]
+    //       | rem_scratch[2*n + 2*base.size() + 2].
+    // u and rem_scratch are dead after the to-Montgomery conversion; u's first n words are
+    // then reused as the exponentiation double-buffer.
+    assert(scratch.size() >= (MAX_PRECOMPUTED + 3) * n + 3 * base.size() + 2);
     const auto u = scratch.subspan(0, n + base.size());
-    const auto base_mont = scratch.subspan(n + base.size(), n);
-    const auto rem_scratch = scratch.subspan(2 * n + base.size(), 2 * n + 2 * base.size() + 2);
+    const auto table = scratch.subspan(n + base.size(), MAX_PRECOMPUTED * n);
+    const auto base_mont = table.first(n);
+    const auto rem_scratch =
+        scratch.subspan(n + base.size() + MAX_PRECOMPUTED * n, 2 * n + 2 * base.size() + 2);
 
+    // Compute base_mont = table[0] = (base * R) % mod, where R = 2^(n*64).
+    // The numerator u = base << (n*64): base in the upper words, lower n words are zero.
     std::ranges::fill(u.first(n), uint64_t{0});  // Lower n words of u must be zero.
     std::ranges::copy(base, u.subspan(n).begin());
     rem(base_mont, u, mod, rem_scratch);
@@ -399,28 +431,58 @@ void modexp_odd(std::span<uint64_t> result, std::span<const uint64_t> base, Expo
     const auto exp_loop = [&]<size_t N>() {
         auto r_cur = std::span<uint64_t, N>{result};
         auto r_tmp = std::span<uint64_t, N>{u.first(n)};
-        const auto bm = std::span<const uint64_t, N>{base_mont};
         const auto m = std::span<const uint64_t, N>{mod};
 
-        std::ranges::copy(bm, r_cur.begin());
-        for (auto i = exp.bit_width() - 1; i != 0; --i)
+        // base_mont^j, for j in 1..table_size.
+        const auto precomputed = [table, n](size_t j) noexcept {
+            return std::span<uint64_t, N>{table.subspan((j - 1) * n, n)};
+        };
+
+        // precomputed(1) = base_mont is already set.
+        for (size_t j = 2; j <= table_size; ++j)
+            mul_amm<N>(precomputed(j), precomputed(j - 1), precomputed(1), m, mod_inv);
+
+        // Reads the `width` exponent bits starting at index `lo`.
+        // TODO: A window spans at most two adjacent bytes, so it could be read with one
+        //   two-byte load, a shift and a mask. Est. 1-3%, and only for a 4-word modulus
+        //   with a very long exponent; measure before doing it.
+        const auto window = [&](size_t lo, size_t width) noexcept {
+            size_t v = 0;
+            for (size_t b = 0; b < width; ++b)
+                v |= size_t{exp[lo + b]} << b;
+            return v;
+        };
+
+        // Windows tile from the bottom, so the ragged one is processed first, at the top.
+        // The top bit is always set, so that first window is nonzero.
+        // TODO: Tiling from the top instead would save w - top_width squarings when
+        //   exp_bits % w != 0 (up to 4%), at the cost of a special-cased final iteration.
+        const size_t top_width = (exp_bits - 1) % w + 1;
+        std::ranges::copy(precomputed(window(exp_bits - top_width, top_width)), r_cur.begin());
+
+        for (size_t pos = exp_bits - top_width; pos != 0;)
         {
-            mul_amm<N>(r_tmp, r_cur, r_cur, m, mod_inv);  // Square.
-            if (exp[i - 1])
-                mul_amm<N>(r_cur, r_tmp, bm, m, mod_inv);  // Multiply.
-            else
+            pos -= w;
+            for (unsigned s = 0; s != w; ++s)  // square w times
+            {
+                mul_amm<N>(r_tmp, r_cur, r_cur, m, mod_inv);
+                std::swap(r_cur, r_tmp);
+            }
+            if (const size_t v = window(pos, w); v != 0)  // multiply by base_mont^v
+            {
+                mul_amm<N>(r_tmp, r_cur, precomputed(v), m, mod_inv);
                 std::swap(r_cur, r_tmp);
+            }
         }
 
-        // Convert from Montgomery form: multiply by 1.
+        // Convert from Montgomery form: multiply by 1. Reuses precomputed(1) storage.
         std::ranges::fill(base_mont, uint64_t{0});
         base_mont[0] = 1;
         mul_amm<N>(r_tmp, r_cur, std::span<const uint64_t, N>{base_mont}, m, mod_inv);
-        std::swap(r_cur, r_tmp);
 
         // If the result ended up in scratch, copy to result.
-        if (r_cur.data() != result.data())
-            std::ranges::copy(r_cur, result.begin());
+        if (r_tmp.data() != result.data())
+            std::ranges::copy(r_tmp, result.begin());
     };
 
     if (n == 4)
@@ -531,10 +593,11 @@ void modexp(std::span<const uint8_t> base_bytes, std::span<const uint8_t> exp_by
 
     // Bump allocator for all working memory (values + scratch).
     // Stack buffer covers inputs up to the EIP-7823 limit (1024 bytes).
-    // Capacity: values[b+2m] + op scratch[4m+3b+2] + CRT[m+2] = 4b+7m+4 words.
+    // Capacity: values[b+2m] + op scratch[(MAX_PRECOMPUTED+3)m+3b+2] + CRT[m+2]
+    //           = 4b + (MAX_PRECOMPUTED+6)m + 4 words.
     // The worst case is an even modulus with 1 trailing zero bit (odd_size=m, pow2_size=1).
     static constexpr size_t MAX_SIZE = 1024 / sizeof(uint64_t);  // EIP-7823
-    static constexpr size_t STACK_CAPACITY = 4 * MAX_SIZE + 7 * MAX_SIZE + 4;
+    static constexpr size_t STACK_CAPACITY = 4 * MAX_SIZE + (6 + MAX_PRECOMPUTED) * MAX_SIZE + 4;
     alignas(uint64_t) std::byte stack_buf[STACK_CAPACITY * sizeof(uint64_t)];
     std::pmr::monotonic_buffer_resource pool{stack_buf, sizeof(stack_buf)};
     std::pmr::polymorphic_allocator<uint64_t> alloc{&pool};
@@ -578,7 +641,8 @@ void modexp(std::span<const uint8_t> base_bytes, std::span<const uint8_t> exp_by
         const auto need_crt = !pow2_is_trivial && !odd_is_trivial;
 
         // Allocate operation scratch (dead after each call, reused sequentially).
-        const size_t odd_scratch = !odd_is_trivial ? 4 * odd_size + 3 * base.size() + 2 : 0;
+        const size_t odd_scratch =
+            !odd_is_trivial ? (MAX_PRECOMPUTED + 3) * odd_size + 3 * base.size() + 2 : 0;
         const size_t pow2_scratch = !pow2_is_trivial ? pow2_size : 0;
         const size_t inv_scratch = need_crt ? 2 * pow2_size : 0;
         const size_t op_scratch_size = std::max({odd_scratch, pow2_scratch, inv_scratch});

test/unittests/precompiles_expmod_test.cpp

@@ -290,6 +290,59 @@ TEST_P(expmod, inputs)
         {"02", "80", "0300000000000000000000000000000000", "0100000000000000000000000000000000"},
         // 2^129 mod (7 * 2^128): carry propagates and is absorbed in nonzero word.
         {"02", "0081", "0700000000000000000000000000000000", "0200000000000000000000000000000000"},
+
+        // Fixed-window exponentiation in modexp_odd. One case per window width w, and
+        // per width of the leading partial window ((exp_bits - 1) % w + 1), which is what
+        // aligns the remaining windows. The exponents are picked so that the windows
+        // consumed cover 0 (multiply skipped), 1 and 2^w - 1 (first and last precomputed
+        // power). Modulus is the secp256k1 field prime: odd, 4 words, so these also cover
+        // the mul_amm<4> specialization.
+        // exp_bits=16, w=1: plain binary square-and-multiply, no table.
+        {"03", "8005", "fffffffffffffffffffffffffffffffffffffffffffffffffffffffefffffc2f",
+            "79c4559d064ab3615f6da729a1f67265b88ee2eaba22838109bea30fb7bee31b"},
+        // exp_bits=17, w=2, leading window 1 bit.
+        {"03", "01001b", "fffffffffffffffffffffffffffffffffffffffffffffffffffffffefffffc2f",
+            "a890a61d8d745fae67a345fb031b048c0cf8952b43622263de0fdc4391a6c6a9"},
+        // exp_bits=18, w=2, leading window 2 bits.
+        {"03", "0200c9", "fffffffffffffffffffffffffffffffffffffffffffffffffffffffefffffc2f",
+            "600614416289329cf72ef906cdfc1dea20339051ec80ed3ff692eb14ed33be81"},
+        // exp_bits=48, w=2: last exponent size before w becomes 3.
+        {"03", "80013b71b865", "fffffffffffffffffffffffffffffffffffffffffffffffffffffffefffffc2f",
+            "fd66fdbe1f0c43e6640c121c366b9061c7f13964a572828c8e3968a50dba847f"},
+        // exp_bits=49, w=3, leading window 1 bit.
+        {"03", "0100d2c92fc182", "fffffffffffffffffffffffffffffffffffffffffffffffffffffffefffffc2f",
+            "651aace134976d8456fcc35686a57cf12670b2e596dabecd0ddae9984ced96c4"},
+        // exp_bits=50, w=3, leading window 2 bits.
+        {"03", "0200a6a7ef231d", "fffffffffffffffffffffffffffffffffffffffffffffffffffffffefffffc2f",
+            "f722a91e1faa3b57f0a19af8d4506b395a0a342e9ee2cbe65cd7a63155d38537"},
+        // exp_bits=51, w=3, leading window 3 bits.
+        {"03", "04013929f7999c", "fffffffffffffffffffffffffffffffffffffffffffffffffffffffefffffc2f",
+            "06f41e370c4ef45a2bc5e1ade1504fbe35e5a42a8f8c2b17ad16a6c657900d48"},
+        // exp_bits=144, w=3: last exponent size before w becomes 4.
+        {"03", "8004cb3ff13151bb9f84a488a5d62e79a680",
+            "fffffffffffffffffffffffffffffffffffffffffffffffffffffffefffffc2f",
+            "97265df41405de7f9b35c1037c349ef367cffd34ed6a86cb933fe14f84bb12d1"},
+        // exp_bits=145, w=4, leading window 1 bit.
+        {"03", "010014b0a1922289f0b19f56c6c373b0e5cd4a",
+            "fffffffffffffffffffffffffffffffffffffffffffffffffffffffefffffc2f",
+            "3587c0d41ce1eb59ec2fa686877d8166aa9740f2410f9271592e5f283e3bd738"},
+        // exp_bits=146, w=4, leading window 2 bits.
+        {"03", "02008d61508c16734bdbe4a9578f4c8185d260",
+            "fffffffffffffffffffffffffffffffffffffffffffffffffffffffefffffc2f",
+            "f65d573e0ba5bdc7cc0e31072eb946ffe5138d0cd4bc936cc1a714d17cdaf954"},
+        // exp_bits=147, w=4, leading window 3 bits.
+        {"03", "040160dce60c2531e93ae750b53938d5b04faf",
+            "fffffffffffffffffffffffffffffffffffffffffffffffffffffffefffffc2f",
+            "0648a7caabfd3d4b972c034830faf933179ed038e1e6a6c4c3ad26f330fe1397"},
+        // exp_bits=148, w=4, leading window 4 bits.
+        {"03", "0802ae8d294c48793907af3e71b536ed84fa84",
+            "fffffffffffffffffffffffffffffffffffffffffffffffffffffffefffffc2f",
+            "40c2770e749bcbf7949855252da0258cc5ae80658427a4af8ba3489a81182ee9"},
+        // Same, with a 5-word modulus: the windowed loop above only ever runs through the
+        // mul_amm<4> specialization, this covers the generic instantiation.
+        {"03", "08f83d563ebc382e09e4b8245edebc817af708",
+            "fffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffefffffc2f",
+            "8016137e4c542dd66f4ab5f668fc0ac76d43353a675f3d4616a56f23757e463ca1093164385ef006"},
     };
 
     for (const auto& [base_hex, exp_hex, mod_hex, expected_result_hex] : test_cases)