crypto: Use fixed-window exponentiation in modexp - #1618
Open
AskAlexSharov wants to merge 11 commits into
Open
Conversation
AskAlexSharov
force-pushed
the
alex/modexp-windowing
branch
2 times, most recently
from
August 1, 2026 13:00
7959fbf to
ef07630
Compare
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.
AskAlexSharov
force-pushed
the
alex/modexp-windowing
branch
from
August 1, 2026 13:13
ef07630 to
006ba7d
Compare
Sahil-4555
pushed a commit
to Sahil-4555/erigon
that referenced
this pull request
Aug 7, 2026
…192, 2^256) (erigontech#22940) See also: ipsilon/evmone#1618 ## Summary Add a fixed-width `uint256` square-and-multiply for MODEXP and route inputs there when the modulus is in `[2^192, 2^256)` and the base fits in 256 bits, instead of to the `evmone` default. This is a **pure addition** — all existing branches are unchanged from `main`. `uint256` avoids arbitrary-precision bookkeeping and the cgo boundary for this case, and is allocation-free. ## Routing ``` mod in {0, 1} -> trivial result base == 1 -> 1 modLen > 32, exp <= 1 byte -> math/big (unchanged from main) 2^192 <= mod < 2^256, baseLen <= 32 -> uint256 (this PR — new) otherwise -> evmone (default; unchanged from main) ``` Only that one class moves. Gas and EIP-7823 limits are untouched. The two routing bounds are not tuned constants, they are where the implementation stops paying: - **`mod >= 2^192`.** `uint256.Reciprocal` returns nothing usable when `m[3] == 0`, and `MulModWithReciprocal` then falls back to a full `udivrem` on every multiply. Since the operand is a byte field, the test is on the modulus *value*, not on `modLen` — a 32-byte field holding a 128-bit modulus stays on `evmone`. - **`baseLen <= 32`.** A wider base has to be folded in before exponentiation, which costs more than the whole `evmone` call. It is also a correctness precondition: `uint256.SetBytes` silently truncates above 32 bytes. ## Benchmarks — time (median ns/op) The tables compare backends directly. Which branch erigon actually uses per row: | row | erigon backend | |---|---| | 256-bit / * | **uint256** (this PR) | | 2048-bit / exp3 | math/big (case unchanged from main) | | 2048-bit / 65537, full | evmone (default) | `pbig` (geth's patched `math/big` fork) and GMP (`mpz_powm`, reth's optional backend, via a tight reused-handle binding) are shown for reference only; neither is used by erigon. **x86_64** | modulus / exponent | Go big | pbig | evmone | GMP | uint256 | |---|--:|--:|--:|--:|--:| | 256-bit / 65537 | 422 | 1,563 | 785 | 709 | 485 | | 256-bit / 64-bit exp | 6,264 | 4,557 | 4,961 | 1,715 | 3,197 | | 256-bit / 256-bit exp | 15,113 | 14,811 | 19,599 | 5,533 | 12,844 | | 2048-bit / exp3 | 3,733 | 3,711 | 5,775 | 4,797 | n/a | | 2048-bit / 65537 | 24,053 | 24,279 | 28,396 | 16,432 | n/a | | 2048-bit / full exp | 2,002,185 | 2,004,973 | 6,138,396 | 1,766,554 | n/a | **arm64** | modulus / exponent | Go big | pbig | evmone | GMP | uint256 | |---|--:|--:|--:|--:|--:| | 256-bit / 65537 | 325 | 1,362 | 468 | 712 | 386 | | 256-bit / 64-bit exp | 6,669 | 4,146 | 2,815 | 2,002 | 2,323 | | 256-bit / 256-bit exp | 13,632 | 13,398 | 11,341 | 7,114 | 9,831 | | 2048-bit / exp3 | 3,096 | 3,102 | 3,502 | 3,820 | n/a | | 2048-bit / 65537 | 19,463 | 20,462 | 15,444 | 12,250 | n/a | | 2048-bit / full exp | 1,692,423 | 1,673,792 | 3,215,552 | 1,289,299 | n/a | For the rows this PR touches, `uint256` is faster than the previous (`evmone`) path on both machines: e.g. x86 `785 -> 485` (65537), `19,599 -> 12,844` (256-bit exp); arm64 `468 -> 386`, `11,341 -> 9,831`. The 2048-bit rows are shown for completeness and are unaffected by this PR. ## Benchmarks — allocations, allocs/op (bytes/op), architecture-independent | modulus / exponent | Go big | pbig | evmone | GMP | uint256 | |---|--:|--:|--:|--:|--:| | 256-bit / 65537 | 8 (352) | 6 (424) | 1 (32) | 2 (16) | 0 (0) | | 256-bit / 64-bit exp | 12 (768) | 6 (424) | 1 (32) | 2 (56) | 0 (0) | | 256-bit / 256-bit exp | 23 (2,097) | 6 (480) | 1 (32) | 2 (56) | 0 (0) | | 2048-bit / exp3 | 8 (3,498) | 8 (3,498) | 1 (256) | 2 (296) | n/a | | 2048-bit / 65537 | 11 (5,227) | 9 (3,114) | 1 (256) | 2 (296) | n/a | | 2048-bit / full exp | 26 (12,941) | 9 (3,444) | 1 (256) | 2 (296) | n/a | In `Run` all paths write into one shared `result` buffer, so these are each algorithm's own temporary allocations. The `uint256` path is allocation-free. (`evmone`'s single Go alloc and GMP's two are the output slices; their working temporaries live in C, outside Go's allocator.) ## Review fix: narrowing the routing (51eb334) The first revision routed on `modLen <= 32` alone and folded an oversized base in byte by byte. @taratorio found that both moved input classes that `evmone` handles better. End-to-end `Run` benchmarks, Apple M4 Max, `-benchtime=300ms -count=6`: | case | before | after | |---|--:|--:| | **control** base 32 B, mod 256-bit, exp 65537 | 386 ns | 392 ns | | **control** base 32 B, mod 256-bit, exp 8 B | 1743 ns | 1753 ns | | **control** base 32 B, mod 2^192+237, exp 8 B | 1765 ns | 1775 ns | | base 32 B, mod 2^192−237, exp 8 B | 4651 ns | 2089 ns | | base 32 B, mod 128-bit, exp 8 B | 3731 ns | 1246 ns | | base 32 B, mod 128-bit, exp 65537 | 732 ns | 328 ns | | base 32 B, mod 64-bit, exp 8 B | 2094 ns | 599 ns | | base 64 B, mod 256-bit, exp 1 | 1313 ns | 191 ns | | base 128 B, mod 256-bit, exp 1 | 2553 ns | 272 ns | | base 1024 B, mod 256-bit, exp 1 | 19,923 ns | 1344 ns | | base 1024 B, mod 256-bit, exp 65537 | 20,413 ns | 1768 ns | | base 1024 B, mod 256-bit, exp 8 B | 21,797 ns | 3594 ns | The controls are the rows that must stay on the `uint256` path, and they do; the ~5 ns is the routing predicate itself. All other rows are back on `evmone`, i.e. at `main`'s numbers. Note the third control (`2^192+237`) versus the fourth row (`2^192−237`): adjacent moduli on opposite sides of the reciprocal boundary. Folding the base 32 bytes at a time (`b = b*(2^256 mod m) + chunk`) instead of guarding was measured too — around 1.6–2 µs for a 1024-byte base against `evmone`'s 1.3 µs, so it would not have recovered the class. ## Correctness - `TestModexpU256Applicable` pins the routing boundaries, including `2^192−1` vs `2^192`, moduli padded with leading zero bytes, and bases of 33 and 1024 bytes. - `TestModexpU256` cross-checks fixed vectors against `math/big` (odd/even/power-of-two moduli, base `0`/`1`/`>mod`, exponent `0`/small/full-width). - `TestModexpU256Random` fuzzes 20,000 random inputs against `math/big`, drawing moduli with randomly zeroed leading bytes so the boundary is hit from both sides. - Existing `TestPrecompiledModExp*` vectors exercise this path and pass. - The reduction is Barrett-style (`uint256.MulModWithReciprocal`), valid for odd and even moduli alike. ## Scope / follow-ups - Only the compute backend changes for the routed class; gas, EIP-7823 limits, result encoding, and the existing `math/big`/`evmone` branches are untouched. - A 4-bit windowed `uint256` gives a further ~1.5x for large exponents but regresses tiny ones, so it's left out to keep this minimal and strictly non-regressing for small exponents. - Moduli below `2^192` could be handled without division by reducing on a narrower fixed width (e.g. 192-bit operands with a 384/192 reduce), but `uint256` exposes no such primitive, so they stay on `evmone`.
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].
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).
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.
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.
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.
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.
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.
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.
There was a problem hiding this comment.
Pull request overview
This PR updates the modexp_odd implementation in evmone_precompiles to use fixed-window exponentiation (with a width chosen by exponent bit-length) to reduce the number of Montgomery multiplications for large exponents, and adds targeted regression/differential-style vectors to exercise the windowing thresholds and edge-window alignments.
Changes:
- Implement fixed-window exponentiation in
modexp_odd, including precomputingb^1 .. b^(2^w - 1)in Montgomery form and processingwexponent bits per multiply. - Increase/adjust scratch sizing and stack-buffer capacity accounting to accommodate the precomputed power table.
- Add unit test vectors covering each window width (
w=1..4) and each possible “leading partial window” width at the threshold boundaries.
Reviewed changes
Copilot reviewed 2 out of 2 changed files in this pull request and generated no comments.
| File | Description |
|---|---|
lib/evmone_precompiles/modexp.cpp |
Replaces binary square-and-multiply in modexp_odd with a fixed-window method; updates scratch/stack capacity calculations accordingly. |
test/unittests/precompiles_expmod_test.cpp |
Adds vectors designed to exercise window sizing thresholds and window-value coverage for the new exponentiation loop. |
💡 Add a code-review agent skill or configure MCP servers for context-aware, tailored reviews. Learn more in the docs.
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%).
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
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Summary
modexp_odduses binary square-and-multiply — one Montgomery multiply per set exponent bit. For large exponents that is roughly twice the multiplies a windowed method needs.This precomputes a small table of base powers (
b^1 .. b^(2^w - 1)in Montgomery form) and consumeswexponent bits per multiply. The window width scales with the exponent size so the table cost stays amortized even for a sparse exponent:w = 1(plain binary, no table) for exponents ≤ 16 bits;w = 2up to 48 bits,w = 3up to 144 bits,w = 4above.With
w = 1the loop is byte-for-byte the previous binary square-and-multiply.Benchmarks
evmone-precompiles-bench --benchmark_filter='modexp<expmod_execute_evmone>', AMD EPYC 4344P, gcc 15.2.0, Release:Small exponents (≤ 16 bits) are unchanged; none regress.
Cost
The power table adds
MODEXP_TABLE_MAX * nwords to the stack scratch buffer (STACK_CAPACITYand themodexp_oddscratch requirement are updated). At the EIP-7823 limit (n = 128words) that is ~15 KB of additional stack in the single modexp frame.The window widths and thresholds are simple, conservative choices; happy to tune them or switch to a sliding window (odd-power table, ~half the entries) if preferred.
Correctness
Adds
expmod.windowing_vs_gmp: a differential test comparing evmone against GMP across many exponent bit-lengths (crossing the window-width thresholds) and bit patterns (all window values), for odd and even moduli. Existingexpmodvectors andlarge_inputscontinue to pass.