diff --git a/lib/evmone_precompiles/modexp.cpp b/lib/evmone_precompiles/modexp.cpp index 923600983b..d54b0a2368 100644 --- a/lib/evmone_precompiles/modexp.cpp +++ b/lib/evmone_precompiles/modexp.cpp @@ -292,6 +292,19 @@ class Exponent const auto bit = (byte >> bit_index) & 1; return bit != 0; } + + /// Returns bits [lo, hi] as an integer, the bit at hi being the most significant. + /// The range must span at most 8 bits, so it covers at most two adjacent bytes. + [[nodiscard]] size_t window(size_t lo, size_t hi) const noexcept + { + assert(lo <= hi && hi - lo < 8); + const auto exp_size = (bit_width_ + 7) / 8; + const auto byte_index = exp_size - 1 - lo / 8; + auto bytes = size_t{data_[byte_index]}; + if (byte_index != 0) // Prepend the next more significant byte if there is one. + bytes |= size_t{data_[byte_index - 1]} << 8; + return (bytes >> (lo % 8)) & ((size_t{1} << (hi + 1 - lo)) - 1); + } }; /// Performs the Almost Montgomery Multiplication (AMM). @@ -369,27 +382,26 @@ template <> } /// Maximum window width used by the windowed method in modexp_odd. -constexpr unsigned MAX_WINDOW_WIDTH = 4; +constexpr unsigned MAX_WINDOW_WIDTH = 5; +static_assert(MAX_WINDOW_WIDTH <= 8, "Exponent::window() covers at most two adjacent bytes"); -/// Number of precomputed values for the max width windowed method. -constexpr size_t MAX_PRECOMPUTED = (size_t{1} << MAX_WINDOW_WIDTH) - 1; +/// Number of precomputed base odd powers 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. +/// Selects the sliding-window width from the exponent bit length. 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) + // Break-even points for a random exponent, where the table's extra multiply stops being + // repaid: 2^(w-1) / (1/(w+1) - 1/(w+2)) = 6, 24, 80, 240. Each narrower width is kept + // one bit longer, which measures better on the sparse small exponents seen in practice. + if (exp_bits <= 7) return 1; - if (exp_bits <= 48) + if (exp_bits <= 25) return 2; - if (exp_bits <= 144) + if (exp_bits <= 81) return 3; + if (exp_bits <= 241) + return 4; return MAX_WINDOW_WIDTH; } @@ -408,7 +420,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 auto 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]. @@ -433,46 +445,50 @@ 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^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)}; + // base_mont^v, for odd v. + const auto precomputed = [table, n](size_t v) noexcept { + return std::span{table.subspan((v / 2) * 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; + // Fill the precomputed table (precomputed(1) is already set). + if (table_size > 1) + { + mul_amm(r_tmp, precomputed(1), precomputed(1), m, mod_inv); // r_tmp = base_mont^2. + for (size_t v = 3; v < 2 * table_size; v += 2) + mul_amm(precomputed(v), precomputed(v - 2), r_tmp, m, mod_inv); + } + + // The widest window of at most w bits ending at bit `hi`, which must be set. Trailing + // zero bits are trimmed off, so the value is odd and only odd table entries are + // needed. Returns the value and the index of its lowest bit. + const auto window = [exp, w](size_t hi) noexcept { + const auto lo = hi + 1 >= w ? hi + 1 - w : size_t{0}; + const auto v = exp.window(lo, hi); + const auto tz = static_cast(std::countr_zero(v)); // v != 0: exp[hi] is set. + return std::pair{v >> tz, lo + tz}; }; - // 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()); + // The top bit is always set, so the first window ends there and is loaded directly. + auto [v_top, pos] = window(exp_bits - 1); + std::ranges::copy(precomputed(v_top), r_cur.begin()); - for (size_t pos = exp_bits - top_width; pos != 0;) + while (pos != 0) { - pos -= w; - for (unsigned s = 0; s != w; ++s) // square w times + --pos; + mul_amm(r_tmp, r_cur, r_cur, m, mod_inv); // Square for this bit. + std::swap(r_cur, r_tmp); + if (!exp[pos]) + continue; + + const auto [v, lo] = window(pos); + for (auto b = lo; b != pos; ++b) // One more square for each remaining window bit. { 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); - } + mul_amm(r_tmp, r_cur, precomputed(v), m, mod_inv); + std::swap(r_cur, r_tmp); + pos = lo; } // Convert from Montgomery form: multiply by 1. Reuses precomputed(1) storage. diff --git a/test/unittests/precompiles_expmod_test.cpp b/test/unittests/precompiles_expmod_test.cpp index 6726ba77d5..b03fe58c6c 100644 --- a/test/unittests/precompiles_expmod_test.cpp +++ b/test/unittests/precompiles_expmod_test.cpp @@ -291,55 +291,70 @@ TEST_P(expmod, inputs) // 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. + // Sliding-window exponentiation in modexp_odd. One case per window width w=1..5. + // Each exponent is built as: top bit (1) | zero run of w+1 bits | one run of w bits + // | trailing zeros, so its windows exercise both the first (b^1) and last + // (b^(2^w-1)) precomputed odd powers, the zero run in between being wide enough to + // keep them in separate windows, traversed by squarings alone. Modulus is the + // secp256k1 field prime: odd, 4 words, so these also cover the mul_amm<4> + // specialization. + // exp_bits=6, w=1: plain binary square-and-multiply, no table. + {"03", "24", "fffffffffffffffffffffffffffffffffffffffffffffffffffffffefffffc2f", + "00000000000000000000000000000000000000000000000002153e468b91c6d1"}, + // exp_bits=10, w=2: windows hit b^1 and b^3. + {"03", "0230", "fffffffffffffffffffffffffffffffffffffffffffffffffffffffefffffc2f", + "e123f780b153ebd75b17a6e7a7133dba60d90a7dbc0f770f08af0055f8e2c7ed"}, + // exp_bits=30, w=3: windows hit b^1 and b^7. + {"03", "21c00000", "fffffffffffffffffffffffffffffffffffffffffffffffffffffffefffffc2f", + "6cc581d10c7d071216edf63238959949056d7cddf5a90711a7c7cdec6b3e861f"}, + // exp_bits=100, w=4: windows hit b^1 and b^15. + {"03", "083c0000000000000000000000", + "fffffffffffffffffffffffffffffffffffffffffffffffffffffffefffffc2f", + "7ff2c68783b688439f7c43de4cbfe265f8875ec726564a442c2cbd1244f6d99e"}, + // exp_bits=254 (mainnet-typical size), w=5: windows hit b^1 and b^31. + {"03", "207c000000000000000000000000000000000000000000000000000000000000", + "fffffffffffffffffffffffffffffffffffffffffffffffffffffffefffffc2f", + "40ea9ce0f6a2c94a7bec98114179d8e1a21287312a25c1fdd7bf46e3d723984a"}, + // Same exponent as the w=5 case above, with a 5-word modulus: the cases above only + // ever run through the mul_amm<4> specialization, this covers the generic + // std::dynamic_extent instantiation at w > 1. + {"03", "207c000000000000000000000000000000000000000000000000000000000000", + "fffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffefffffc2f", + "aa50260a96f69a722fc965bbfec20c21195eda68068b20e9899976f80ed8f6d4f6816bec10fc4ee6"}, + + // Random exponents straddling the width thresholds of the fixed-window + // implementation this replaced, which no longer coincide with the bands above: + // exp_bits 16..18 (w=2), 48..51 (w=3), 144..148 (w=4). Same modulus as above, + // except for the last case, which repeats exp_bits=148 with a 5-word modulus. {"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"},