From a03c1a1b2aaede69187e5631e4af2d33dd646670 Mon Sep 17 00:00:00 2001 From: Alexey Sharov Date: Sat, 1 Aug 2026 19:53:29 +0700 Subject: [PATCH 01/16] crypto: Use fixed-window exponentiation in modexp modexp_odd used binary square-and-multiply: one Montgomery multiply per set exponent bit. For large exponents that is roughly twice the multiplies a windowed method needs. Precompute a small table of base powers (b^1 .. b^(2^w - 1) in Montgomery form) and consume w exponent bits per multiply. The window width scales with the exponent size (w = 1..4) so the table cost stays amortized even for a sparse exponent, and small exponents keep the plain binary path (w = 1). With w = 1 the loop is identical to the previous binary square-and-multiply. Measured ~1.5-1.6x on large-exponent modexp (256-bit modulus, 256-bit exponent: 19.6us -> 12.9us; 4096-bit modulus, 8192-bit exponent: 96.5ms -> 60.6ms on an AMD EPYC 4344P); smaller exponents also improve and none regress. The power table adds MODEXP_TABLE_MAX*n words to the stack scratch buffer. Add expmod.windowing_vs_gmp: a differential test against GMP over many exponent bit-lengths and patterns, for odd and even moduli. --- lib/evmone_precompiles/modexp.cpp | 97 +++++++++++++++++----- test/unittests/precompiles_expmod_test.cpp | 53 ++++++++++++ 2 files changed, 129 insertions(+), 21 deletions(-) diff --git a/lib/evmone_precompiles/modexp.cpp b/lib/evmone_precompiles/modexp.cpp index 9e96ce3c6a..efd9b5b1cd 100644 --- a/lib/evmone_precompiles/modexp.cpp +++ b/lib/evmone_precompiles/modexp.cpp @@ -370,6 +370,11 @@ template <> /// 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(). +/// Maximum fixed-window width used by modexp_odd, and the resulting size of the +/// precomputed power table (b^1 .. b^(2^w - 1)). Bounds the extra scratch space. +constexpr unsigned MODEXP_WINDOW_MAX = 4; +constexpr size_t MODEXP_TABLE_MAX = (size_t{1} << MODEXP_WINDOW_MAX) - 1; + void modexp_odd(std::span result, std::span base, Exponent exp, std::span mod, std::span scratch) noexcept { @@ -380,39 +385,87 @@ void modexp_odd(std::span result, std::span base, Expo const auto n = mod.size(); const auto mod_inv = -evmmax::modinv(mod[0]); - - // 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); - - // 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. + const auto exp_bits = exp.bit_width(); + + // Fixed-window exponentiation width. w == 1 is plain binary square-and-multiply + // (no table); a wider window precomputes b^1..b^(2^w-1) once and then does one + // multiply per w exponent bits instead of one per set bit. The width scales with + // the exponent size so the table (which grows as 2^w) stays amortized even for a + // sparse exponent, and small exponents stay on the plain binary path. + const unsigned w = [exp_bits]() -> unsigned { + if (exp_bits > 144) + return MODEXP_WINDOW_MAX; + if (exp_bits > 48) + return 3; + if (exp_bits > 16) + return 2; + return 1; + }(); + const size_t table_size = (size_t{1} << w) - 1; // entries b^1 .. b^(2^w - 1) + + // Layout: u[n+base.size()] | table[MODEXP_TABLE_MAX*n] | rem_scratch[2n+2b+2]. + // table[0] doubles as base_mont (b^1 in Montgomery form). rem_scratch is only + // live during the initial to-Montgomery conversion. + assert(scratch.size() >= (MODEXP_TABLE_MAX + 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(), MODEXP_TABLE_MAX * n); + const auto base_mont = table.first(n); + const auto rem_scratch = + scratch.subspan(n + base.size() + MODEXP_TABLE_MAX * 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); // Double-buffer exponentiation loop, parameterized by mul_amm size. 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}; + auto r_cur = std::span{result}; + auto r_tmp = std::span{u.first(n)}; + + // Precompute the odd/all-powers table: table[j] = base^(j+1) in Montgomery + // form. table[0] = base_mont is already set; for w == 1 the loop is empty. + for (size_t j = 1; j < table_size; ++j) + { + const auto prev = std::span{table.subspan((j - 1) * n, n)}; + const auto cur = std::span{table.subspan(j * n, n)}; + mul_amm(cur, prev, bm, m, mod_inv); // b^(j+1) = b^j * b + } - std::ranges::copy(bm, r_cur.begin()); - for (auto i = exp.bit_width() - 1; i != 0; --i) + // Reads the w-bit (or fewer) window whose most-significant bit is at index `hi`. + const auto window = [&](size_t hi, size_t width) noexcept { + size_t v = 0; + for (size_t b = 0; b < width; ++b) + v = (v << 1) | (exp[hi - b] ? size_t{1} : size_t{0}); + return v; + }; + + // Process the most-significant (possibly short) window first, then full + // w-bit windows. The top bit is always set, so the first window is nonzero. + const size_t top_width = (exp_bits - 1) % w + 1; + const size_t top_val = window(exp_bits - 1, top_width); + std::ranges::copy(table.subspan((top_val - 1) * n, n), 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 - 1, w); v != 0) // multiply by b^v + { + const auto tv = std::span{table.subspan((v - 1) * n, n)}; + mul_amm(r_tmp, r_cur, tv, m, mod_inv); std::swap(r_cur, r_tmp); + } } - // Convert from Montgomery form: multiply by 1. + // Convert from Montgomery form: multiply by 1. Reuses table[0] 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); @@ -531,10 +584,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[(TABLE_MAX+4)m+3b+2] + CRT[m+2]. + // The op scratch grows by the modexp_odd power table (MODEXP_TABLE_MAX*m 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 + (7 + MODEXP_TABLE_MAX) * 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 +632,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 ? (MODEXP_TABLE_MAX + 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..730cadfb21 100644 --- a/test/unittests/precompiles_expmod_test.cpp +++ b/test/unittests/precompiles_expmod_test.cpp @@ -455,3 +455,56 @@ TEST(expmod, huge_inputs_analysis) EXPECT_EQ(max_output_size, expected_output_size); } } + +#ifdef EVMONE_PRECOMPILES_GMP +namespace +{ +/// Runs modexp through a specific implementation and returns the result bytes. +evmc::bytes run_expmod( + ExpmodExecuteFn fn, const evmc::bytes& base, const evmc::bytes& exp, const evmc::bytes& mod) +{ + evmc::bytes input(3 * 32, 0); + using namespace intx; + be::unsafe::store(&input[0], uint256{base.size()}); + be::unsafe::store(&input[32], uint256{exp.size()}); + be::unsafe::store(&input[64], uint256{mod.size()}); + input += base; + input += exp; + input += mod; + evmc::bytes result(mod.size(), 0xfe); + const auto [status, output_size] = fn(input.data(), input.size(), result.data(), result.size()); + EXPECT_EQ(status, EVMC_SUCCESS); + EXPECT_EQ(output_size, mod.size()); + return result; +} +} // namespace + +// Differential test for the fixed-window exponentiation in modexp_odd: exercises +// many exponent bit-lengths (crossing the binary<->window threshold) and bit +// patterns (all window values) against the GMP reference, for odd and even moduli. +TEST(expmod, windowing_vs_gmp) +{ + const auto base = make_val(32, 0xab, 0xcd, 0xa5); + for (const auto& mod : { + make_val(32, 0xff, 0xff, 0xff), // 2^256-1 (odd: direct window path) + make_val(32, 0xff, 0xfe, 0xff), // even (odd part via CRT still windows) + }) + { + for (size_t size = 2; size <= 40; ++size) // exponent 16..320 bits + { + for (const uint8_t msb : {uint8_t{0x01}, uint8_t{0x80}, uint8_t{0xff}}) + { + for (const uint8_t fill : {uint8_t{0x00}, uint8_t{0xa5}, uint8_t{0xff}}) + { + const auto exp = make_val(size, msb, 0x01, fill); + const auto ev = + run_expmod(&evmone::state::expmod_execute_evmone, base, exp, mod); + const auto gm = run_expmod(&evmone::state::expmod_execute_gmp, base, exp, mod); + EXPECT_EQ(ev, gm) + << "exp_size=" << size << " msb=" << +msb << " fill=" << +fill; + } + } + } + } +} +#endif From c797cdb7f0b3e2122f0cd457429d41139704a56c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Pawe=C5=82=20Bylica?= Date: Sat, 8 Aug 2026 20:06:24 +0200 Subject: [PATCH 02/16] crypto: Fix modexp_odd doc comment and scratch requirement The window constants were inserted between modexp_odd's doc comment and the function, so the whole block documented MODEXP_WINDOW_MAX instead and the function was left undocumented. The stated scratch requirement was also left at the pre-windowing value. It is now (MODEXP_TABLE_MAX + 3)*n + 3*base.size() + 2 words, matching both the assert in modexp_odd and the odd_scratch computation in modexp(): u[n + b] + table[MODEXP_TABLE_MAX*n] + rem_scratch[2n + 2b + 2]. --- lib/evmone_precompiles/modexp.cpp | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/lib/evmone_precompiles/modexp.cpp b/lib/evmone_precompiles/modexp.cpp index efd9b5b1cd..14422dff1f 100644 --- a/lib/evmone_precompiles/modexp.cpp +++ b/lib/evmone_precompiles/modexp.cpp @@ -368,13 +368,14 @@ template <> mul_amm_256(r, x, y, mod, mod_inv); } -/// 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(). /// Maximum fixed-window width used by modexp_odd, and the resulting size of the /// precomputed power table (b^1 .. b^(2^w - 1)). Bounds the extra scratch space. constexpr unsigned MODEXP_WINDOW_MAX = 4; constexpr size_t MODEXP_TABLE_MAX = (size_t{1} << MODEXP_WINDOW_MAX) - 1; +/// Computes result[] = base[]^exp % mod[] for odd mod[] (mod[0] % 2 != 0). +/// Scratch space required: (MODEXP_TABLE_MAX + 3)*n + 3*base.size() + 2 words, +/// where n = mod.size(). void modexp_odd(std::span result, std::span base, Exponent exp, std::span mod, std::span scratch) noexcept { From 729c0e5f194b6eeb2d2df55e94fdb08bb78f59e8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Pawe=C5=82=20Bylica?= Date: Sat, 8 Aug 2026 20:08:54 +0200 Subject: [PATCH 03/16] crypto: Restore exact modexp stack buffer accounting The pre-windowing capacity 4b + 7m + 4 was exactly the sum of the allocations: mod[m] + result[m] + op scratch[4m + 3b + 2] + result_odd[m] + CRT[2], i.e. the m coefficient was 2 + 4 + 1. Windowing raised the op scratch m term from 4 to MODEXP_TABLE_MAX + 3, so the total is 2 + (TABLE_MAX + 3) + 1 = TABLE_MAX + 6, but MODEXP_TABLE_MAX was added to the old 7 instead, which double-counts one m. Use (6 + MODEXP_TABLE_MAX): 3332 -> 3204 words, exactly the worst-case demand (1024-byte base and modulus, even modulus with 1 trailing zero bit). --- lib/evmone_precompiles/modexp.cpp | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/lib/evmone_precompiles/modexp.cpp b/lib/evmone_precompiles/modexp.cpp index 14422dff1f..e4aaad1746 100644 --- a/lib/evmone_precompiles/modexp.cpp +++ b/lib/evmone_precompiles/modexp.cpp @@ -585,11 +585,12 @@ 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[(TABLE_MAX+4)m+3b+2] + CRT[m+2]. + // Capacity: values[b+2m] + op scratch[(TABLE_MAX+3)m+3b+2] + CRT[m+2] + // = 4b + (TABLE_MAX+6)m + 4 words. // The op scratch grows by the modexp_odd power table (MODEXP_TABLE_MAX*m 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 + MODEXP_TABLE_MAX) * MAX_SIZE + 4; + static constexpr size_t STACK_CAPACITY = 4 * MAX_SIZE + (6 + MODEXP_TABLE_MAX) * 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}; From e31f86b44d20fbe17241ea1a654e6bb81216565a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Pawe=C5=82=20Bylica?= Date: Sat, 8 Aug 2026 20:09:38 +0200 Subject: [PATCH 04/16] crypto: Rename the modexp window constants MODEXP_ is redundant for file-local constants in modexp.cpp, and TABLE_MAX did not say what the table holds. Use MAX_WINDOW_WIDTH and MAX_PRECOMPUTED, the latter matching the MAX_SIZE ordering already used in this file, and split the shared doc comment so each constant documents itself. --- lib/evmone_precompiles/modexp.cpp | 32 ++++++++++++++++--------------- 1 file changed, 17 insertions(+), 15 deletions(-) diff --git a/lib/evmone_precompiles/modexp.cpp b/lib/evmone_precompiles/modexp.cpp index e4aaad1746..a0d81b39cc 100644 --- a/lib/evmone_precompiles/modexp.cpp +++ b/lib/evmone_precompiles/modexp.cpp @@ -368,13 +368,15 @@ template <> mul_amm_256(r, x, y, mod, mod_inv); } -/// Maximum fixed-window width used by modexp_odd, and the resulting size of the -/// precomputed power table (b^1 .. b^(2^w - 1)). Bounds the extra scratch space. -constexpr unsigned MODEXP_WINDOW_MAX = 4; -constexpr size_t MODEXP_TABLE_MAX = (size_t{1} << MODEXP_WINDOW_MAX) - 1; +/// Maximum fixed-window width used by modexp_odd. +constexpr unsigned MAX_WINDOW_WIDTH = 4; + +/// Number of base powers b^1 .. b^(2^w - 1) precomputed for the widest window. +/// Bounds the extra scratch space taken by the power table. +constexpr size_t MAX_PRECOMPUTED = (size_t{1} << MAX_WINDOW_WIDTH) - 1; /// Computes result[] = base[]^exp % mod[] for odd mod[] (mod[0] % 2 != 0). -/// Scratch space required: (MODEXP_TABLE_MAX + 3)*n + 3*base.size() + 2 words, +/// Scratch space required: (MAX_PRECOMPUTED + 3)*n + 3*base.size() + 2 words, /// where n = mod.size(). void modexp_odd(std::span result, std::span base, Exponent exp, std::span mod, std::span scratch) noexcept @@ -395,7 +397,7 @@ void modexp_odd(std::span result, std::span base, Expo // sparse exponent, and small exponents stay on the plain binary path. const unsigned w = [exp_bits]() -> unsigned { if (exp_bits > 144) - return MODEXP_WINDOW_MAX; + return MAX_WINDOW_WIDTH; if (exp_bits > 48) return 3; if (exp_bits > 16) @@ -404,15 +406,15 @@ void modexp_odd(std::span result, std::span base, Expo }(); const size_t table_size = (size_t{1} << w) - 1; // entries b^1 .. b^(2^w - 1) - // Layout: u[n+base.size()] | table[MODEXP_TABLE_MAX*n] | rem_scratch[2n+2b+2]. + // Layout: u[n+base.size()] | table[MAX_PRECOMPUTED*n] | rem_scratch[2n+2b+2]. // table[0] doubles as base_mont (b^1 in Montgomery form). rem_scratch is only // live during the initial to-Montgomery conversion. - assert(scratch.size() >= (MODEXP_TABLE_MAX + 3) * n + 3 * base.size() + 2); + assert(scratch.size() >= (MAX_PRECOMPUTED + 3) * n + 3 * base.size() + 2); const auto u = scratch.subspan(0, n + base.size()); - const auto table = scratch.subspan(n + base.size(), MODEXP_TABLE_MAX * n); + 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() + MODEXP_TABLE_MAX * n, 2 * n + 2 * base.size() + 2); + 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. @@ -585,12 +587,12 @@ 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[(TABLE_MAX+3)m+3b+2] + CRT[m+2] - // = 4b + (TABLE_MAX+6)m + 4 words. - // The op scratch grows by the modexp_odd power table (MODEXP_TABLE_MAX*m words). + // Capacity: values[b+2m] + op scratch[(MAX_PRECOMPUTED+3)m+3b+2] + CRT[m+2] + // = 4b + (MAX_PRECOMPUTED+6)m + 4 words. + // The op scratch grows by the modexp_odd power table (MAX_PRECOMPUTED*m 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 + (6 + MODEXP_TABLE_MAX) * 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}; @@ -635,7 +637,7 @@ void modexp(std::span base_bytes, std::span exp_by // Allocate operation scratch (dead after each call, reused sequentially). const size_t odd_scratch = - !odd_is_trivial ? (MODEXP_TABLE_MAX + 3) * odd_size + 3 * base.size() + 2 : 0; + !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}); From de2cfe374b6d53d8a352b4ac58870541629c98b2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Pawe=C5=82=20Bylica?= Date: Sat, 8 Aug 2026 20:31:49 +0200 Subject: [PATCH 05/16] test: Cover the modexp window widths without GMP The windowing was covered by a differential test against GMP, which only builds in the precompiles-gmp CI job; a default build got no coverage of the new code at all. Replace it with vectors in the existing expmod input table. Instrumenting modexp_odd over the previous table shows it only ever reached w=1 (65 calls) and w=4 (1 call): the w=2 and w=3 bands were never executed, and w=4 was sampled at a single exponent size. The new vectors add one case per window width and per width of the leading partial window, with exponents whose windows cover 0 (multiply skipped), 1 and 2^w-1 (first and last precomputed power). Coverage becomes w=1/2/3/4 = 66/3/4/5 calls. Mutation testing the windowing code (15 hand-written mutants) goes from 11 to 13 killed. The two mutants that only alter top_width when the leading window is partial, e.g. top_width = (exp_bits - 1) % w + 1 -> top_width = min(w, exp_bits) are no-ops at exp_bits=256 (the old table's only w>1 case) and survived before. The two remaining survivors change the window width only, which is a performance parameter: every width computes the same result, so no correctness test can kill them. --- test/unittests/precompiles_expmod_test.cpp | 101 ++++++++++----------- 1 file changed, 48 insertions(+), 53 deletions(-) diff --git a/test/unittests/precompiles_expmod_test.cpp b/test/unittests/precompiles_expmod_test.cpp index 730cadfb21..c2dfbb0b60 100644 --- a/test/unittests/precompiles_expmod_test.cpp +++ b/test/unittests/precompiles_expmod_test.cpp @@ -290,6 +290,54 @@ 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"}, }; for (const auto& [base_hex, exp_hex, mod_hex, expected_result_hex] : test_cases) @@ -455,56 +503,3 @@ TEST(expmod, huge_inputs_analysis) EXPECT_EQ(max_output_size, expected_output_size); } } - -#ifdef EVMONE_PRECOMPILES_GMP -namespace -{ -/// Runs modexp through a specific implementation and returns the result bytes. -evmc::bytes run_expmod( - ExpmodExecuteFn fn, const evmc::bytes& base, const evmc::bytes& exp, const evmc::bytes& mod) -{ - evmc::bytes input(3 * 32, 0); - using namespace intx; - be::unsafe::store(&input[0], uint256{base.size()}); - be::unsafe::store(&input[32], uint256{exp.size()}); - be::unsafe::store(&input[64], uint256{mod.size()}); - input += base; - input += exp; - input += mod; - evmc::bytes result(mod.size(), 0xfe); - const auto [status, output_size] = fn(input.data(), input.size(), result.data(), result.size()); - EXPECT_EQ(status, EVMC_SUCCESS); - EXPECT_EQ(output_size, mod.size()); - return result; -} -} // namespace - -// Differential test for the fixed-window exponentiation in modexp_odd: exercises -// many exponent bit-lengths (crossing the binary<->window threshold) and bit -// patterns (all window values) against the GMP reference, for odd and even moduli. -TEST(expmod, windowing_vs_gmp) -{ - const auto base = make_val(32, 0xab, 0xcd, 0xa5); - for (const auto& mod : { - make_val(32, 0xff, 0xff, 0xff), // 2^256-1 (odd: direct window path) - make_val(32, 0xff, 0xfe, 0xff), // even (odd part via CRT still windows) - }) - { - for (size_t size = 2; size <= 40; ++size) // exponent 16..320 bits - { - for (const uint8_t msb : {uint8_t{0x01}, uint8_t{0x80}, uint8_t{0xff}}) - { - for (const uint8_t fill : {uint8_t{0x00}, uint8_t{0xa5}, uint8_t{0xff}}) - { - const auto exp = make_val(size, msb, 0x01, fill); - const auto ev = - run_expmod(&evmone::state::expmod_execute_evmone, base, exp, mod); - const auto gm = run_expmod(&evmone::state::expmod_execute_gmp, base, exp, mod); - EXPECT_EQ(ev, gm) - << "exp_size=" << size << " msb=" << +msb << " fill=" << +fill; - } - } - } - } -} -#endif From 1787aa1771070075adbd3206df2e798849e9a2e6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Pawe=C5=82=20Bylica?= Date: Sat, 8 Aug 2026 21:13:15 +0200 Subject: [PATCH 06/16] crypto: Document where the modexp window thresholds come from 16, 48 and 144 read as arbitrary. They are the break-even points for a random exponent, exp_bits = 2^w / ((1-2^-w)/w - (1-2^-(w+1))/(w+1)) = 16, 48, 140. A random exponent is the worst case for windowing, so each band already errs towards the smaller window: a dense exponent prefers the next width up at every threshold. Verified by benchmarking forced widths 1..5 over exponents of 8..512 bits at mod_len 32 and 256, with a random and an all-ones exponent (27 sizes x 4 sets). The chosen width is within 4% of the best width in 94 of 108 points, and every larger miss is a case where a wider window would have won: up to 11% at exp_bits=48 and 29% for a dense exponent at exp_bits<=16, both given up deliberately. Choosing a width that is too large never costs more than 3.2%. Closed forms were tried and rejected: the thresholds grow by ~2.9x per width, so any bit_width()-based rule lands on 2x spacing. The best of them, clamp(bit_width(exp_bits) >> 1, 1, 4), has a better worst case (9.5% vs 28.9%) but only by preferring wider windows earlier, which loses up to 9.5% on the random exponent this heuristic is tuned for. --- lib/evmone_precompiles/modexp.cpp | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/lib/evmone_precompiles/modexp.cpp b/lib/evmone_precompiles/modexp.cpp index a0d81b39cc..271c37808b 100644 --- a/lib/evmone_precompiles/modexp.cpp +++ b/lib/evmone_precompiles/modexp.cpp @@ -392,9 +392,13 @@ void modexp_odd(std::span result, std::span base, Expo // Fixed-window exponentiation width. w == 1 is plain binary square-and-multiply // (no table); a wider window precomputes b^1..b^(2^w-1) once and then does one - // multiply per w exponent bits instead of one per set bit. The width scales with - // the exponent size so the table (which grows as 2^w) stays amortized even for a - // sparse exponent, and small exponents stay on the plain binary path. + // multiply per w exponent bits instead of one per set bit. + // + // The thresholds are the break-even points for a random exponent, where the 2^w extra + // table multiplies stop being repaid by the sparser multiplies in the loop: + // exp_bits = 2^w / ((1-2^-w)/w - (1-2^-(w+1))/(w+1)), giving 16, 48 and 140. + // A random exponent is the worst case for windowing, so this errs towards the smaller + // window: a dense exponent would prefer the next width up at every threshold. const unsigned w = [exp_bits]() -> unsigned { if (exp_bits > 144) return MAX_WINDOW_WIDTH; From 1ef93a8391b607dc2373dfc8b5d621cbad6fc9d0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Pawe=C5=82=20Bylica?= Date: Sat, 8 Aug 2026 21:13:24 +0200 Subject: [PATCH 07/16] crypto: Fix Exponent::operator[] index documentation The comment said e[0] is the top bit, but the implementation indexes from the least significant bit: byte_index = index / 8 selects data_[exp_size - 1 - ...], i.e. the last byte, so e[0] is the bottom bit and the top bit is at bit_width() - 1. Both the old binary loop and the window reads rely on the actual behaviour; only the comment was wrong. --- lib/evmone_precompiles/modexp.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/lib/evmone_precompiles/modexp.cpp b/lib/evmone_precompiles/modexp.cpp index 271c37808b..02bb49a459 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. From dd7b231e3fd2579046f5bf2f1e45d248b7fe8284 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Pawe=C5=82=20Bylica?= Date: Sun, 9 Aug 2026 01:51:43 +0200 Subject: [PATCH 08/16] crypto: Correct the modexp windowing comments Review follow-ups, comments only: - "odd/all-powers table" described the odd-powers table of a sliding window, which this is not. The table holds every power b^1..b^(2^w - 1). - The threshold derivation gave 140 while the code tests 144, with nothing saying the value had been rounded. - The scratch layout noted the lifetime of rem_scratch, which is never reused, and omitted the reuse that does matter: u's first n words become the exponentiation double-buffer once the to-Montgomery conversion is done. --- lib/evmone_precompiles/modexp.cpp | 15 +++++++++------ 1 file changed, 9 insertions(+), 6 deletions(-) diff --git a/lib/evmone_precompiles/modexp.cpp b/lib/evmone_precompiles/modexp.cpp index 02bb49a459..06287b5fb3 100644 --- a/lib/evmone_precompiles/modexp.cpp +++ b/lib/evmone_precompiles/modexp.cpp @@ -396,7 +396,8 @@ void modexp_odd(std::span result, std::span base, Expo // // The thresholds are the break-even points for a random exponent, where the 2^w extra // table multiplies stop being repaid by the sparser multiplies in the loop: - // exp_bits = 2^w / ((1-2^-w)/w - (1-2^-(w+1))/(w+1)), giving 16, 48 and 140. + // exp_bits = 2^w / ((1-2^-w)/w - (1-2^-(w+1))/(w+1)), giving 16, 48 and 140 (rounded + // up to 144 below; the two widths are within 0.2% of each other over 140..144). // A random exponent is the worst case for windowing, so this errs towards the smaller // window: a dense exponent would prefer the next width up at every threshold. const unsigned w = [exp_bits]() -> unsigned { @@ -410,9 +411,11 @@ void modexp_odd(std::span result, std::span base, Expo }(); const size_t table_size = (size_t{1} << w) - 1; // entries b^1 .. b^(2^w - 1) - // Layout: u[n+base.size()] | table[MAX_PRECOMPUTED*n] | rem_scratch[2n+2b+2]. - // table[0] doubles as base_mont (b^1 in Montgomery form). rem_scratch is only - // live during the initial to-Montgomery conversion. + // Layout: u[n + base.size()] | table[MAX_PRECOMPUTED*n] + // | rem_scratch[2*n + 2*base.size() + 2]. + // table[0] doubles as base_mont (b^1 in Montgomery form). Both u and rem_scratch are + // dead after the to-Montgomery conversion, and 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 table = scratch.subspan(n + base.size(), MAX_PRECOMPUTED * n); @@ -433,8 +436,8 @@ void modexp_odd(std::span result, std::span base, Expo auto r_cur = std::span{result}; auto r_tmp = std::span{u.first(n)}; - // Precompute the odd/all-powers table: table[j] = base^(j+1) in Montgomery - // form. table[0] = base_mont is already set; for w == 1 the loop is empty. + // Precompute the power table: table[j] = base^(j+1) in Montgomery form, i.e. all + // of b^1..b^(2^w-1). table[0] = base_mont is already set; empty for w == 1. for (size_t j = 1; j < table_size; ++j) { const auto prev = std::span{table.subspan((j - 1) * n, n)}; From 2062f6a680dbb334fdcafc463ca595ffd9e577dc Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Pawe=C5=82=20Bylica?= Date: Sun, 9 Aug 2026 12:15:33 +0200 Subject: [PATCH 09/16] crypto: Record the modexp windowing follow-ups as TODOs Four measured improvements that are out of scope for this change, with the magnitude of each so the next reader can judge whether to bother: - Sliding window (odd-powers table only). This is what GMP's mpn_powm and OpenSSL's BN_mod_exp_mont both use for modexp. Half the table for a given width, worth ~3-7% here. - Worst-case width thresholds. Gas is charged on exponent bit length, not Hamming weight, so the densest exponent is the costliest input at a given charge; tuning for it yields a closed form and +10.7% worst-case time per gas. - Ragged window at the bottom rather than the top: up to 4%, but exactly 0 when exp_bits is a multiple of w, which covers the common 256/512/2048/8192 sizes. - Batched window reads: below the noise floor except in one corner. Also note that the ragged-window-at-the-top layout is the standard m-ary one (blst's ec_mult.h computes the same "top excess bits modulo window size"), so the third item is a deviation from common practice, not a fix. --- lib/evmone_precompiles/modexp.cpp | 29 +++++++++++++++++++++++++++++ 1 file changed, 29 insertions(+) diff --git a/lib/evmone_precompiles/modexp.cpp b/lib/evmone_precompiles/modexp.cpp index 06287b5fb3..01cbe0cc40 100644 --- a/lib/evmone_precompiles/modexp.cpp +++ b/lib/evmone_precompiles/modexp.cpp @@ -400,6 +400,19 @@ void modexp_odd(std::span result, std::span base, Expo // up to 144 below; the two widths are within 0.2% of each other over 140..144). // A random exponent is the worst case for windowing, so this errs towards the smaller // window: a dense exponent would prefer the next width up at every threshold. + // + // TODO: Switch to a sliding window, as GMP's mpn_powm and OpenSSL's BN_mod_exp_mont + // both do. Its table holds only the odd powers b^1, b^3 .. b^(2^w-1), so it is half + // the size for a given width, which buys one extra width at the same memory. + // Break-even points become 2^(w-1) / (1/(w+1) - 1/(w+2)) = 6, 24, 80, 240, 672 — + // GMP's win_size() uses exactly these (7, 25, 81, 241, 673). Measured ~3-7% over + // the fixed window here. + // TODO: Tune the thresholds for the worst case rather than the average. Gas is charged + // on exponent bit length, not Hamming weight, so the costliest input at a given + // charge is the densest exponent, whose break-even points are 2^w*w*(w+1) = 4, 24, + // 96, 320. Those have bit widths 3, 5, 7, 9, so the width collapses to a closed form, + // min(MAX_WINDOW_WIDTH, (bit_width(exp_bits) + 1) / 2). Measured +10.7% worst-case + // time per gas over the whole modexp benchmark matrix. const unsigned w = [exp_bits]() -> unsigned { if (exp_bits > 144) return MAX_WINDOW_WIDTH; @@ -446,6 +459,13 @@ void modexp_odd(std::span result, std::span base, Expo } // Reads the w-bit (or fewer) window whose most-significant bit is at index `hi`. + // One Exponent::operator[] per exponent bit, the same count as the binary loop + // this replaced, so windowing added no extraction work per bit. + // + // TODO: A w <= 4 bit field spans at most two adjacent bytes, so a window could be + // read with one two-byte load, a shift and a mask instead of w indexed bit + // reads. Below the noise floor everywhere except the n=4 / 8192-bit corner, + // where it is worth an estimated 1-3% of cycles; measure before doing it. const auto window = [&](size_t hi, size_t width) noexcept { size_t v = 0; for (size_t b = 0; b < width; ++b) @@ -455,6 +475,15 @@ void modexp_odd(std::span result, std::span base, Expo // Process the most-significant (possibly short) window first, then full // w-bit windows. The top bit is always set, so the first window is nonzero. + // This is the usual m-ary layout: windows tile from the bottom, so the ragged + // one lands at the top (blst's ec_mult.h does the same, "top excess bits + // modulo target window size"). + // + // TODO: Putting the ragged window at the bottom instead would use a full-width + // top window and save w - top_width squarings for the same multiply count. + // Worth 0% whenever exp_bits % w == 0 (so nothing at 256/512/2048/8192 bits), + // but 4.0% at exp_bits=17, 2.9% at 49, 2.0% at 33 and 1.6% at 145. Costs a + // special-cased final iteration. const size_t top_width = (exp_bits - 1) % w + 1; const size_t top_val = window(exp_bits - 1, top_width); std::ranges::copy(table.subspan((top_val - 1) * n, n), r_cur.begin()); From b48fbcf6386433fa84fcec41a0ecefbd84530d31 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Pawe=C5=82=20Bylica?= Date: Sun, 9 Aug 2026 14:14:22 +0200 Subject: [PATCH 10/16] crypto: Trim the modexp windowing comments The windowing comments had grown to explain more than the code needed, with the same facts stated in several places and measurement detail that belongs in a commit message. Extract the window-width lambda as window_width(), so the choice is named rather than described, and hang the one non-rederivable fact (where the thresholds come from) off it as a doc comment. Codegen is unchanged. Drop the rest: the random-vs-dense justification, the "table holds b^1..b^(2^w-1)" fact repeated at four sites, the blst precedent note, the comparison against the binary loop this replaced, and a "grows by" remark that only made sense relative to the previous revision of the file. Keep the TODOs but reduce each to its actionable point. modexp.cpp goes from 167 comment lines to 140 (24% -> 20%). --- lib/evmone_precompiles/modexp.cpp | 92 ++++++++++++------------------- 1 file changed, 34 insertions(+), 58 deletions(-) diff --git a/lib/evmone_precompiles/modexp.cpp b/lib/evmone_precompiles/modexp.cpp index 01cbe0cc40..16cfbe354a 100644 --- a/lib/evmone_precompiles/modexp.cpp +++ b/lib/evmone_precompiles/modexp.cpp @@ -375,6 +375,27 @@ constexpr unsigned MAX_WINDOW_WIDTH = 4; /// Bounds the extra scratch space taken by the power table. constexpr size_t MAX_PRECOMPUTED = (size_t{1} << MAX_WINDOW_WIDTH) - 1; +/// Selects the fixed-window width: one multiply per w exponent bits instead of one per set +/// bit, at the cost of a 2^w - 1 entry table. w == 1 is plain square-and-multiply. +/// The thresholds are the break-even points for a random exponent, +/// exp_bits = 2^w / ((1-2^-w)/w - (1-2^-(w+1))/(w+1)) ≈ 16, 48, 140 (rounded up to 144). +/// +/// TODO: Switch to a sliding window, as GMP's mpn_powm and OpenSSL's BN_mod_exp_mont do: +/// the table then holds only odd powers, half the entries per width. 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 +{ + if (exp_bits > 144) + return MAX_WINDOW_WIDTH; + if (exp_bits > 48) + return 3; + if (exp_bits > 16) + return 2; + return 1; +} + /// Computes result[] = base[]^exp % mod[] for odd mod[] (mod[0] % 2 != 0). /// Scratch space required: (MAX_PRECOMPUTED + 3)*n + 3*base.size() + 2 words, /// where n = mod.size(). @@ -390,45 +411,13 @@ void modexp_odd(std::span result, std::span base, Expo const auto mod_inv = -evmmax::modinv(mod[0]); const auto exp_bits = exp.bit_width(); - // Fixed-window exponentiation width. w == 1 is plain binary square-and-multiply - // (no table); a wider window precomputes b^1..b^(2^w-1) once and then does one - // multiply per w exponent bits instead of one per set bit. - // - // The thresholds are the break-even points for a random exponent, where the 2^w extra - // table multiplies stop being repaid by the sparser multiplies in the loop: - // exp_bits = 2^w / ((1-2^-w)/w - (1-2^-(w+1))/(w+1)), giving 16, 48 and 140 (rounded - // up to 144 below; the two widths are within 0.2% of each other over 140..144). - // A random exponent is the worst case for windowing, so this errs towards the smaller - // window: a dense exponent would prefer the next width up at every threshold. - // - // TODO: Switch to a sliding window, as GMP's mpn_powm and OpenSSL's BN_mod_exp_mont - // both do. Its table holds only the odd powers b^1, b^3 .. b^(2^w-1), so it is half - // the size for a given width, which buys one extra width at the same memory. - // Break-even points become 2^(w-1) / (1/(w+1) - 1/(w+2)) = 6, 24, 80, 240, 672 — - // GMP's win_size() uses exactly these (7, 25, 81, 241, 673). Measured ~3-7% over - // the fixed window here. - // TODO: Tune the thresholds for the worst case rather than the average. Gas is charged - // on exponent bit length, not Hamming weight, so the costliest input at a given - // charge is the densest exponent, whose break-even points are 2^w*w*(w+1) = 4, 24, - // 96, 320. Those have bit widths 3, 5, 7, 9, so the width collapses to a closed form, - // min(MAX_WINDOW_WIDTH, (bit_width(exp_bits) + 1) / 2). Measured +10.7% worst-case - // time per gas over the whole modexp benchmark matrix. - const unsigned w = [exp_bits]() -> unsigned { - if (exp_bits > 144) - return MAX_WINDOW_WIDTH; - if (exp_bits > 48) - return 3; - if (exp_bits > 16) - return 2; - return 1; - }(); - const size_t table_size = (size_t{1} << w) - 1; // entries b^1 .. b^(2^w - 1) + const auto w = window_width(exp_bits); + const size_t table_size = (size_t{1} << w) - 1; // Layout: u[n + base.size()] | table[MAX_PRECOMPUTED*n] // | rem_scratch[2*n + 2*base.size() + 2]. - // table[0] doubles as base_mont (b^1 in Montgomery form). Both u and rem_scratch are - // dead after the to-Montgomery conversion, and u's first n words are then reused as - // the exponentiation double-buffer. + // 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 table = scratch.subspan(n + base.size(), MAX_PRECOMPUTED * n); @@ -449,23 +438,18 @@ void modexp_odd(std::span result, std::span base, Expo auto r_cur = std::span{result}; auto r_tmp = std::span{u.first(n)}; - // Precompute the power table: table[j] = base^(j+1) in Montgomery form, i.e. all - // of b^1..b^(2^w-1). table[0] = base_mont is already set; empty for w == 1. + // table[j] = b^(j+1) in Montgomery form; table[0] = base_mont is already set. for (size_t j = 1; j < table_size; ++j) { const auto prev = std::span{table.subspan((j - 1) * n, n)}; const auto cur = std::span{table.subspan(j * n, n)}; - mul_amm(cur, prev, bm, m, mod_inv); // b^(j+1) = b^j * b + mul_amm(cur, prev, bm, m, mod_inv); } // Reads the w-bit (or fewer) window whose most-significant bit is at index `hi`. - // One Exponent::operator[] per exponent bit, the same count as the binary loop - // this replaced, so windowing added no extraction work per bit. - // - // TODO: A w <= 4 bit field spans at most two adjacent bytes, so a window could be - // read with one two-byte load, a shift and a mask instead of w indexed bit - // reads. Below the noise floor everywhere except the n=4 / 8192-bit corner, - // where it is worth an estimated 1-3% of cycles; measure before doing it. + // 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 hi, size_t width) noexcept { size_t v = 0; for (size_t b = 0; b < width; ++b) @@ -473,17 +457,10 @@ void modexp_odd(std::span result, std::span base, Expo return v; }; - // Process the most-significant (possibly short) window first, then full - // w-bit windows. The top bit is always set, so the first window is nonzero. - // This is the usual m-ary layout: windows tile from the bottom, so the ragged - // one lands at the top (blst's ec_mult.h does the same, "top excess bits - // modulo target window size"). - // - // TODO: Putting the ragged window at the bottom instead would use a full-width - // top window and save w - top_width squarings for the same multiply count. - // Worth 0% whenever exp_bits % w == 0 (so nothing at 256/512/2048/8192 bits), - // but 4.0% at exp_bits=17, 2.9% at 49, 2.0% at 33 and 1.6% at 145. Costs a - // special-cased final iteration. + // 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; const size_t top_val = window(exp_bits - 1, top_width); std::ranges::copy(table.subspan((top_val - 1) * n, n), r_cur.begin()); @@ -625,7 +602,6 @@ void modexp(std::span base_bytes, std::span exp_by // Stack buffer covers inputs up to the EIP-7823 limit (1024 bytes). // Capacity: values[b+2m] + op scratch[(MAX_PRECOMPUTED+3)m+3b+2] + CRT[m+2] // = 4b + (MAX_PRECOMPUTED+6)m + 4 words. - // The op scratch grows by the modexp_odd power table (MAX_PRECOMPUTED*m 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 + (6 + MAX_PRECOMPUTED) * MAX_SIZE + 4; From 23ebacc1a5a66273547091212deb1798cf4d223c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Pawe=C5=82=20Bylica?= Date: Sun, 9 Aug 2026 17:54:19 +0200 Subject: [PATCH 11/16] crypto: Drop a reorder and cover the generic mul_amm instantiation MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Restore the original order of the exp_loop locals. Moving r_cur/r_tmp below bm/m changed nothing — the four initializers are independent — and only made the diff four lines longer. Add one vector with a 5-word modulus. Every other windowed case uses the 4-word secp256k1 prime, so instrumenting the dispatch shows the windowed loop was only ever reaching the mul_amm<4> specialization; the generic instantiation ran at w=1 only. Coverage by (instantiation, width) becomes: generic w=1 63 calls mul_amm<4> w=1 3 calls generic w=4 1 call mul_amm<4> w=2 3 calls mul_amm<4> w=3 4 calls mul_amm<4> w=4 5 calls --- lib/evmone_precompiles/modexp.cpp | 4 ++-- test/unittests/precompiles_expmod_test.cpp | 5 +++++ 2 files changed, 7 insertions(+), 2 deletions(-) diff --git a/lib/evmone_precompiles/modexp.cpp b/lib/evmone_precompiles/modexp.cpp index 16cfbe354a..f11fda841c 100644 --- a/lib/evmone_precompiles/modexp.cpp +++ b/lib/evmone_precompiles/modexp.cpp @@ -433,10 +433,10 @@ void modexp_odd(std::span result, std::span base, Expo // Double-buffer exponentiation loop, parameterized by mul_amm size. const auto exp_loop = [&]() { - const auto bm = std::span{base_mont}; - const auto m = std::span{mod}; 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}; // table[j] = b^(j+1) in Montgomery form; table[0] = base_mont is already set. for (size_t j = 1; j < table_size; ++j) diff --git a/test/unittests/precompiles_expmod_test.cpp b/test/unittests/precompiles_expmod_test.cpp index c2dfbb0b60..6726ba77d5 100644 --- a/test/unittests/precompiles_expmod_test.cpp +++ b/test/unittests/precompiles_expmod_test.cpp @@ -338,6 +338,11 @@ TEST_P(expmod, inputs) {"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) From 3c15cf8544fbee4de37033b3fcc8548860d6d873 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Pawe=C5=82=20Bylica?= Date: Sun, 9 Aug 2026 19:39:56 +0200 Subject: [PATCH 12/16] crypto: Address review comments on window_width - Order the window_width conditions by ascending exponent size. - Inline mod.size() in the scratch requirement rather than naming it for a single use; the line still fits. - Deduce table_size. --- lib/evmone_precompiles/modexp.cpp | 17 ++++++++--------- 1 file changed, 8 insertions(+), 9 deletions(-) diff --git a/lib/evmone_precompiles/modexp.cpp b/lib/evmone_precompiles/modexp.cpp index f11fda841c..bc903f4986 100644 --- a/lib/evmone_precompiles/modexp.cpp +++ b/lib/evmone_precompiles/modexp.cpp @@ -387,18 +387,17 @@ constexpr size_t MAX_PRECOMPUTED = (size_t{1} << MAX_WINDOW_WIDTH) - 1; /// min(MAX_WINDOW_WIDTH, (bit_width(exp_bits) + 1) / 2). Measured +10.7% worst case. constexpr unsigned window_width(size_t exp_bits) noexcept { - if (exp_bits > 144) - return MAX_WINDOW_WIDTH; - if (exp_bits > 48) - return 3; - if (exp_bits > 16) + if (exp_bits <= 16) + return 1; + if (exp_bits <= 48) return 2; - return 1; + 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: (MAX_PRECOMPUTED + 3)*n + 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 { @@ -412,7 +411,7 @@ void modexp_odd(std::span result, std::span base, Expo const auto exp_bits = exp.bit_width(); const auto w = window_width(exp_bits); - const size_t table_size = (size_t{1} << w) - 1; + const auto table_size = (size_t{1} << w) - 1; // Layout: u[n + base.size()] | table[MAX_PRECOMPUTED*n] // | rem_scratch[2*n + 2*base.size() + 2]. From 9bdc46cf58fdbece7d1632dc759378774835c7cd Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Pawe=C5=82=20Bylica?= Date: Sun, 9 Aug 2026 23:13:08 +0200 Subject: [PATCH 13/16] cleanups --- lib/evmone_precompiles/modexp.cpp | 19 +++++++------------ 1 file changed, 7 insertions(+), 12 deletions(-) diff --git a/lib/evmone_precompiles/modexp.cpp b/lib/evmone_precompiles/modexp.cpp index bc903f4986..3e781f1503 100644 --- a/lib/evmone_precompiles/modexp.cpp +++ b/lib/evmone_precompiles/modexp.cpp @@ -368,20 +368,15 @@ template <> mul_amm_256(r, x, y, mod, mod_inv); } -/// Maximum fixed-window width used by modexp_odd. +/// Maximum window width used by the windowed method in modexp_odd. constexpr unsigned MAX_WINDOW_WIDTH = 4; -/// Number of base powers b^1 .. b^(2^w - 1) precomputed for the widest window. -/// Bounds the extra scratch space taken by the power table. +/// 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: one multiply per w exponent bits instead of one per set -/// bit, at the cost of a 2^w - 1 entry table. w == 1 is plain square-and-multiply. -/// The thresholds are the break-even points for a random exponent, -/// exp_bits = 2^w / ((1-2^-w)/w - (1-2^-(w+1))/(w+1)) ≈ 16, 48, 140 (rounded up to 144). +/// Selects the fixed-window width off exponent bits. /// -/// TODO: Switch to a sliding window, as GMP's mpn_powm and OpenSSL's BN_mod_exp_mont do: -/// the table then holds only odd powers, half the entries per width. Measured ~3-7%. +/// 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. @@ -437,7 +432,7 @@ void modexp_odd(std::span result, std::span base, Expo const auto bm = std::span{base_mont}; const auto m = std::span{mod}; - // table[j] = b^(j+1) in Montgomery form; table[0] = base_mont is already set. + // table[j] = base_mont^(j+1); table[0] = base_mont is already set. for (size_t j = 1; j < table_size; ++j) { const auto prev = std::span{table.subspan((j - 1) * n, n)}; @@ -472,7 +467,7 @@ void modexp_odd(std::span result, std::span base, Expo 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 - 1, w); v != 0) // multiply by b^v + if (const size_t v = window(pos + w - 1, w); v != 0) // multiply by base^v { const auto tv = std::span{table.subspan((v - 1) * n, n)}; mul_amm(r_tmp, r_cur, tv, m, mod_inv); @@ -600,7 +595,7 @@ 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[(MAX_PRECOMPUTED+3)m+3b+2] + CRT[m+2] - // = 4b + (MAX_PRECOMPUTED+6)m + 4 words. + // = 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 + (6 + MAX_PRECOMPUTED) * MAX_SIZE + 4; From a2fcd4e883661b8aa0e45f486ad0a6a77369289c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Pawe=C5=82=20Bylica?= Date: Sun, 9 Aug 2026 23:43:09 +0200 Subject: [PATCH 14/16] crypto: Introduce precomputed() and fix window_width comments Replace the three repeated span-cast expressions over the power table with a precomputed(j) helper, capturing only what it uses. Saves ~120 instructions of generated code in modexp(); no runtime difference, the offsets are computed once per window against a whole mul_amm. Also: move the threshold derivation into window_width's body (a caller does not need it), fix "off exponent bits", and use base_mont consistently in the two comments naming the table contents. --- lib/evmone_precompiles/modexp.cpp | 24 +++++++++++++----------- 1 file changed, 13 insertions(+), 11 deletions(-) diff --git a/lib/evmone_precompiles/modexp.cpp b/lib/evmone_precompiles/modexp.cpp index 3e781f1503..d7e573b640 100644 --- a/lib/evmone_precompiles/modexp.cpp +++ b/lib/evmone_precompiles/modexp.cpp @@ -374,7 +374,7 @@ 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 off exponent bits. +/// 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 @@ -382,6 +382,8 @@ constexpr size_t MAX_PRECOMPUTED = (size_t{1} << MAX_WINDOW_WIDTH) - 1; /// 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) @@ -432,13 +434,14 @@ void modexp_odd(std::span result, std::span base, Expo const auto bm = std::span{base_mont}; const auto m = std::span{mod}; - // table[j] = base_mont^(j+1); table[0] = base_mont is already set. + // The j-th precomputed value. + const auto precomputed = [table, n](size_t j) noexcept { + return std::span{table.subspan(j * n, n)}; + }; + + // precomputed[j] = base_mont^(j+1); precomputed[0] = base_mont is already set. for (size_t j = 1; j < table_size; ++j) - { - const auto prev = std::span{table.subspan((j - 1) * n, n)}; - const auto cur = std::span{table.subspan(j * n, n)}; - mul_amm(cur, prev, bm, m, mod_inv); - } + mul_amm(precomputed(j), precomputed(j - 1), bm, m, mod_inv); // Reads the w-bit (or fewer) window whose most-significant bit is at index `hi`. // TODO: A window spans at most two adjacent bytes, so it could be read with one @@ -457,7 +460,7 @@ void modexp_odd(std::span result, std::span base, Expo // 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; const size_t top_val = window(exp_bits - 1, top_width); - std::ranges::copy(table.subspan((top_val - 1) * n, n), r_cur.begin()); + std::ranges::copy(precomputed(top_val - 1), r_cur.begin()); for (size_t pos = exp_bits - top_width; pos != 0;) { @@ -467,10 +470,9 @@ void modexp_odd(std::span result, std::span base, Expo 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 - 1, w); v != 0) // multiply by base^v + if (const size_t v = window(pos + w - 1, w); v != 0) // multiply by base_mont^v { - const auto tv = std::span{table.subspan((v - 1) * n, n)}; - mul_amm(r_tmp, r_cur, tv, m, mod_inv); + mul_amm(r_tmp, r_cur, precomputed(v - 1), m, mod_inv); std::swap(r_cur, r_tmp); } } From 9e2ff4a0a1245cb1839dba4c89e78f293c7f857e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Pawe=C5=82=20Bylica?= Date: Mon, 10 Aug 2026 00:44:11 +0200 Subject: [PATCH 15/16] crypto: Remove index compensations in the windowing loop Four bookkeeping simplifications, no change to the multiplies performed or the exponent bits read: - window() now takes the low bit index instead of the high one, so callers pass the position they already have rather than re-deriving the top bit. The second call site was adding back the w that the line above had just subtracted. - precomputed(k) is base_mont^k rather than base_mont^(k+1), which drops the -1 from both value-bearing call sites and makes the code read like its comment. - bm was left over from the pre-windowing loop and had one use; precomputed(1) says the same thing. - The from-Montgomery multiply leaves its result in r_tmp, so test and copy that instead of swapping first. --- lib/evmone_precompiles/modexp.cpp | 33 ++++++++++++++----------------- 1 file changed, 15 insertions(+), 18 deletions(-) diff --git a/lib/evmone_precompiles/modexp.cpp b/lib/evmone_precompiles/modexp.cpp index d7e573b640..2eb9fa8450 100644 --- a/lib/evmone_precompiles/modexp.cpp +++ b/lib/evmone_precompiles/modexp.cpp @@ -431,26 +431,25 @@ 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}; - // The j-th precomputed value. - const auto precomputed = [table, n](size_t j) noexcept { - return std::span{table.subspan(j * n, n)}; + // base_mont^k, for k in 1..table_size. + const auto precomputed = [table, n](size_t k) noexcept { + return std::span{table.subspan((k - 1) * n, n)}; }; - // precomputed[j] = base_mont^(j+1); precomputed[0] = base_mont is already set. - for (size_t j = 1; j < table_size; ++j) - mul_amm(precomputed(j), precomputed(j - 1), bm, m, mod_inv); + // precomputed(1) = base_mont is already set. + for (size_t k = 2; k <= table_size; ++k) + mul_amm(precomputed(k), precomputed(k - 1), precomputed(1), m, mod_inv); - // Reads the w-bit (or fewer) window whose most-significant bit is at index `hi`. + // 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 hi, size_t width) noexcept { + const auto window = [&](size_t lo, size_t width) noexcept { size_t v = 0; for (size_t b = 0; b < width; ++b) - v = (v << 1) | (exp[hi - b] ? size_t{1} : size_t{0}); + v |= size_t{exp[lo + b]} << b; return v; }; @@ -459,8 +458,7 @@ void modexp_odd(std::span result, std::span base, Expo // 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; - const size_t top_val = window(exp_bits - 1, top_width); - std::ranges::copy(precomputed(top_val - 1), r_cur.begin()); + std::ranges::copy(precomputed(window(exp_bits - top_width, top_width)), r_cur.begin()); for (size_t pos = exp_bits - top_width; pos != 0;) { @@ -470,22 +468,21 @@ void modexp_odd(std::span result, std::span base, Expo 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 - 1, w); v != 0) // multiply by base_mont^v + if (const size_t v = window(pos, w); v != 0) // multiply by base_mont^v { - mul_amm(r_tmp, r_cur, precomputed(v - 1), m, mod_inv); + mul_amm(r_tmp, r_cur, precomputed(v), m, mod_inv); std::swap(r_cur, r_tmp); } } - // Convert from Montgomery form: multiply by 1. Reuses table[0] storage. + // 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) From c29da9a42ac9297d6ace040a0faed59f6de60ab2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Pawe=C5=82=20Bylica?= Date: Mon, 10 Aug 2026 00:51:02 +0200 Subject: [PATCH 16/16] crypto: Keep the loop variable named j The 1-based precomputed() indexing did not need a different letter. --- lib/evmone_precompiles/modexp.cpp | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/lib/evmone_precompiles/modexp.cpp b/lib/evmone_precompiles/modexp.cpp index 2eb9fa8450..923600983b 100644 --- a/lib/evmone_precompiles/modexp.cpp +++ b/lib/evmone_precompiles/modexp.cpp @@ -433,14 +433,14 @@ void modexp_odd(std::span result, std::span base, Expo auto r_tmp = std::span{u.first(n)}; const auto m = std::span{mod}; - // base_mont^k, for k in 1..table_size. - const auto precomputed = [table, n](size_t k) noexcept { - return std::span{table.subspan((k - 1) * n, n)}; + // 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 k = 2; k <= table_size; ++k) - mul_amm(precomputed(k), precomputed(k - 1), precomputed(1), m, mod_inv); + 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