diff --git a/lib/evmone_precompiles/modexp.cpp b/lib/evmone_precompiles/modexp.cpp index 9e96ce3c6a..923600983b 100644 --- a/lib/evmone_precompiles/modexp.cpp +++ b/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 result, std::span base, Exponent exp, std::span mod, std::span scratch) noexcept { @@ -380,17 +405,24 @@ void modexp_odd(std::span result, std::span 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 result, std::span base, Expo const auto exp_loop = [&]() { auto r_cur = std::span{result}; auto r_tmp = std::span{u.first(n)}; - const auto bm = std::span{base_mont}; const auto m = std::span{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{table.subspan((j - 1) * n, n)}; + }; + + // precomputed(1) = base_mont is already set. + for (size_t j = 2; j <= table_size; ++j) + mul_amm(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(r_tmp, r_cur, r_cur, m, mod_inv); // Square. - if (exp[i - 1]) - mul_amm(r_cur, r_tmp, bm, m, mod_inv); // Multiply. - else + pos -= w; + for (unsigned s = 0; s != w; ++s) // square w times + { + mul_amm(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(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(r_tmp, r_cur, std::span{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 base_bytes, std::span 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 alloc{&pool}; @@ -578,7 +641,8 @@ void modexp(std::span base_bytes, std::span 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}); diff --git a/test/unittests/precompiles_expmod_test.cpp b/test/unittests/precompiles_expmod_test.cpp index 1f487853a0..6726ba77d5 100644 --- a/test/unittests/precompiles_expmod_test.cpp +++ b/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)