Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
106 changes: 61 additions & 45 deletions lib/evmone_precompiles/modexp.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -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).
Expand Down Expand Up @@ -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;
}

Expand All @@ -408,7 +420,7 @@ void modexp_odd(std::span<uint64_t> result, std::span<const uint64_t> 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].
Expand All @@ -433,46 +445,50 @@ void modexp_odd(std::span<uint64_t> result, std::span<const uint64_t> base, Expo
auto r_tmp = std::span<uint64_t, N>{u.first(n)};
const auto m = std::span<const uint64_t, N>{mod};

// 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)};
// base_mont^v, for odd v.
const auto precomputed = [table, n](size_t v) noexcept {
return std::span<uint64_t, N>{table.subspan((v / 2) * 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;
// Fill the precomputed table (precomputed(1) is already set).
if (table_size > 1)
{
mul_amm<N>(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<N>(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<size_t>(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<N>(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<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);
}
mul_amm<N>(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.
Expand Down
55 changes: 35 additions & 20 deletions test/unittests/precompiles_expmod_test.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Comment thread
chfast marked this conversation as resolved.
// 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"},
Expand Down