Porting libsecp256k1's Modular Inversion to Pasta

Field inversion is one of the more expensive primitive operations in pasta_curves, the implementation of the Pasta fields used by Ironwood, Orchard, and halo2. Until recently, it computed an inverse using Fermat's little theorem:

This is simple and reliable, not to mention constant time, but expensive. For a roughly 255-bit field, exponentiation requires hundreds of dependent field squarings and multiplications. Pasta field elements are also represented in Montgomery form, so those operations include Montgomery reductions.

For settings where variable-time inversion is acceptable, we replaced this with a port of the modular inversion implementation from Bitcoin Core's libsecp256k1: Peter Dettman's signed-62 variable-time modular inversion, based on the Bernstein-Yang safegcd algorithm. We then specialized it for the Pasta fields and their Montgomery representation.

On our current Apple aarch64 baseline with the assembly backend, Fp::invert went from about 3.44 μs to 756 ns, a speedup of about 4.6×.

Inversion does not require exponentiation

Fermat's little theorem is one way to compute an inverse, but inversion is also an extended GCD problem. For a nonzero value modulo a prime ,

so an extended GCD computation can recover a coefficient corresponding to .

The Bernstein-Yang safegcd algorithm organizes this computation around small operations called divsteps. Conceptually, it maintains two integers and coefficient values . A sequence of divsteps transforms them with a small matrix

and after 62 divsteps, the integer state can be updated as

The corresponding coefficient state is updated similarly modulo . The decisions made by the next 62 divsteps depend only on the low 62 bits of and , so full-width arithmetic is not needed to discover what those 62 steps will do.

libsecp256k1 therefore computes the entire 2×2 transition matrix using ordinary machine-word arithmetic. Only after the matrix is known does it apply the result to the full-width values.

For a 255-bit integer, represented as five signed radix- limbs, that replaces 62 rounds of full-width updates with one. The signed representation is useful because the GCD computation naturally produces negative intermediate values. Using 62 bits per limb leaves enough room in a 64-bit word for the sign and the intermediate growth that occurs during these updates. This is the main reason the implementation is fast.

What the 62-bit kernel is doing

The variable-time kernel in libsecp256k1 does somewhat more than execute 62 elementary divsteps literally. When has trailing zero bits, it removes several of them at once. When is odd, small modular-inverse formulas choose a multiple of that cancels several low bits of . Depending on the state, the implementation can cancel up to four or six bits at a time.

During this process it updates only the small matrix coefficients . At the end of the batch, the matrix has determinant

and applying it to the full values produces numerators exactly divisible by .

The coefficient update has one additional problem. It needs to compute something of the form

while remaining equivalent modulo . Because is odd, it has an inverse modulo . If

then can be chosen from the low limb so that

and the division is an exact shift. This structure is taken directly from libsecp256k1. It uses ordinary integer multiplication and addition, but it does not perform a sequence of field multiplications and Montgomery reductions like the Fermat approach.

Staying in Montgomery form

The first specialization handles the representation of field elements. Like many high-performance prime-field implementations, pasta_curves stores a field element in Montgomery form, as the integer

where modulo the field modulus. Inverting the stored integer directly gives , which is off by two factors of from the stored form of the inverse, . Converting on either side of the inversion would supply the correction, but would add field multiplications and Montgomery reductions around it.

The standard alternative is to fold the correction into the coefficient seeds, since the extended GCD's coefficient updates are linear in and and scaling the initial coefficients scales the result. The Pasta implementation maintains the scaled invariants

beginning with

and the divstep machinery runs unchanged. At termination,

and the first invariant gives

the stored form of the inverse, with no conversion and no Montgomery reduction anywhere in the inversion.

The Pasta moduli are sparse in the signed-62 radix

The second specialization comes from a convenient property shared by the Pallas and Vesta field moduli. The signed-62 implementation represents the modulus as five radix- limbs, and both Pasta moduli have the form

In other words,

This matters during the coefficient updates: each update adds a correction so that its low 62 bits vanish, and with an arbitrary five-limb modulus, applying that correction involves multiplication by each modulus limb. For Pasta:

• is a general multiplication.
• is a general multiplication.
• is a shift.
• disappears.
• is another shift.

Only two of the five modulus limbs require real multiplications. The production kernels specialize directly for this shape, and compile-time assertions verify that both field moduli still have exactly the expected radix- representation. If the field constants were ever changed incompatibly, the code would stop compiling rather than silently using an invalid optimization.

Specializing the ends of the computation

There are also useful facts about the first and last 62-divstep batches. At the beginning,

so the first full-width update has much more structure than a generic iteration. We have a separate first-batch kernel that uses the sparse form of when updating , and skips products against the known-zero coefficient .

At the other end, the implementation updates before updating . Once an update produces

the GCD computation is finished and only the coefficient row that will become the inverse is still needed. The terminal batch therefore calls update_d_only_62 instead of computing both and .

This is a small reordering relative to upstream. The reference implementation in the test suite computes the state in the original order and verifies that the specialized driver produces the same result.

The variable-time driver also reduces the active limb count as and get smaller, so later full-width updates frequently operate on fewer than five limbs. None of these changes alters the underlying safegcd algorithm; they remove work that the particular Pasta invocation does not need.

Performance

We measured the change against the current Fermat implementation on Apple aarch64 with the aarch64-asm backend.

That is approximately a 4.6× reduction in inversion time. The comparison uses the current baseline: other recent arithmetic work had already made the previous Fermat implementation substantially faster, so older measurements of this change showed a larger relative improvement.

Against the current assembly field multiplication, an inversion is now roughly 43 multiplication-equivalents. That matters beyond direct calls to invert: inversions appear in affine point arithmetic, batch normalization, batch inversion, and other higher-level cryptographic operations (see our batch affine ladder work), and making the primitive cheaper changes which coordinate systems and batching strategies are profitable.