Bernstein-Yang's Inverse
Presentation: https://youtu.be/t0_iII-nKNw
Introduction
GCD-based algorithms are widely used as an efficient method for computing modular inverses. However, most of them involve input-dependent branching, which makes them vulnerable to side-channel attacks.
For example, an implementation of Stein’s binary GCD algorithm in OpenSSL was shown to be susceptible to an RSA key extraction attack. This kind of attack demonstrates that the execution path of conditional branches can leak secret information.
A simple countermeasure is to use Fermat’s Little Theorem, which computes the inverse as . This method can be implemented in constant time, but it is generally less efficient.
This raises the question: Is there a way to compute the modular inverse that is both constant-time secure and faster than the Fermat-based approach?
Background
Limitations of Existing Constant-Time GCD Algorithms
The Binary GCD algorithm (also known as Stein’s Algorithm) involves branching based on the parity of the input integers or comparisons between them, which causes the execution path to vary. For example:
If is even → right shift
If is even → right shift
If →
Branching based on whether or is even is relatively safe since it only inspects the least significant bit. However, branching on comparisons like often requires examining multiple bits (or limbs), which makes it difficult to ensure constant-time execution.
Protocol Explanation
The original paper also covers the case of computing inverses of polynomials, but for simplicity, we will restrict ourselves to the case of integers. The explanation closely follows the implementation described in https://github.com/bitcoin-core/secp256k1/blob/master/doc/safegcd_implementation.md, while adhering to the original paper's formulation of the divstep
update as much as possible.
GCD with divstep
divstep
When is an odd integer and is any integer, the divstep
function is defined as follows:
This operation is designed to execute in constant-time, relying only on fixed-bit values like and the least significant bit of to determine behavior, ensuring that only a single bit needs to be checked.
Stages
The divstep
operation can be viewed as a composition of two functions:
Conditional Swap
Condition: and is odd
Operation:
TODO(chokobole): explain the reason why helps to reduce to 0.
Elimination
Always performed
Purpose: to reduce the size of
Operation:
Since is always odd, is guaranteed, allowing simplification:
This division by 2 is always valid because:
If is odd, then is even (odd - odd = even)
If is even, the result is trivially divisible by 2
GCD Preservation
This logic is rooted in the classical Euclidean algorithm:
Thus, divstep
preserves the GCD while continually reducing the size of . Eventually, reaches 0, and the GCD remains in .
Simplified divstep
divstep
The behavior of divstep
can be more concretely classified as:
As previously discussed, this operation ultimately returns the GCD.
Code
def gcd(a, b):
"""Compute the GCD of an odd integer a and another integer b."""
assert a & 1 # require a to be odd
delta = 1 # additional state variable
while b != 0:
assert a & 1 # a will be odd in every iteration
if delta > 0 and b & 1:
delta, a, b = 1 - delta, b, (a - b) // 2
elif b & 1:
delta, a, b = 1 + delta, a, (b - a) // 2
else:
delta, a, b = 1 + delta, a, b // 2
return abs(a)
Modular Inverse with divstep
divstep
The idea of using divstep
to compute the modular inverse becomes most clear when we assume that the modulus p
is an odd prime. In this section, we explain the procedure for computing the inverse of an integer , i.e., .
Initial Setup
We initialize the variables as , , , and , and maintain the following invariants throughout:
At each iteration, we apply a divstep
transformation to the pair . Since is a prime and , we know that . This guarantees that repeated applications of divstep
will eventually reduce to , and to either or .
Formally, after a sufficient number of iterations, we reach:
At this point, because the algorithm maintains the invariant , we conclude:
Since , multiplying by gives the correct modular inverse of .
The divstep
function branches into three cases depending on and the parity of . We analyze how the invariants are preserved in each case.
Case 1: and is odd
We use the rule for Case 1. Accordingly, the values of and are updated as:
Since is guaranteed to be odd, the subtraction is always even, making the division by 2 valid.
We apply the same update rule to and as follows:
Case 2: and is odd
We use the rule for Case 2 Accordingly, the values of and are updated as:
Since is guaranteed to be odd, the subtraction is always even, making the division by 2 valid.
We apply the same update rule to and as follows:
Case 3: is even
We use the rule for Case 3 Accordingly, the values of and are updated as:
We apply the same update rule to and as follows:
Code
def div2(p, a):
"""Helper routine to compute a/2 mod p (where p is odd)."""
assert p & 1
if a & 1: # If a is odd, make it even by adding p.
a += p
# a must be even now, so clean division by 2 is possible.
return a // 2
def modinv(p, a):
"""Compute the inverse of a mod p (given that it exists, and p is odd)."""
assert p & 1
delta, u, v = 1, p, a
x1, x2 = 0, 1
while v != 0:
# Note that while division by two for u and v is only ever done on even inputs,
# this is not true for x1 and x2, so we need the div2 helper function.
if delta > 0 and v & 1:
delta, u, v = 1 - delta, v, (u - v) // 2
x1, x2 = x2, div2(p, x1 - x2)
elif v & 1:
delta, u, v = 1 + delta, u, (v - u) // 2
x1, x2 = x1, div2(p, x2 - x1)
else:
delta, u, v = 1 + delta, u, v // 2
x1, x2 = x1, div2(p, x2)
# Verify that the invariants x1=u/a mod p, x2=v/a mod p are maintained.
assert u % p == (a * x1) % p
assert v % p == (a * x2) % p
assert u == 1 or u == -1 # |u| is the GCD, it must be 1
# Because of invariant x1 = u/a mod p, 1/a = x1/u mod p. Since |u| = 1, x1/u = x1 * u.
return (x1 * u) % p
Batching Multiple Divsteps
Motivation
The divstep
operation involves division by 2 for the variables . Instead of performing division by 2 on and at every iteration, we can batch divstep
s together and apply them to and all at once. This reduces individual divisions into a single division.
Transition Matrix
A single divstep
can be represented using a matrix multiplication. The one-step update is expressed as:
Here, is or depending on and that make sure the expression inside the parentheses becomes divisible by 2.
The transition matrix used to update is the same one applied to . The key idea is to compute the transition matrix from and over steps, and then apply it to using matrix-vector multiplication as follows:
Here, and are correction terms chosen so that the expression inside the parentheses becomes divisible by . This will be explained in here.
The transition matrix is defined as follows:
Upper Bound on Iteration Count: Theorem 11.2
Let , which can usually be considered the bit length of the modulus. Then, the number of N-step matrix applications is bounded by:
Example: For the BN254 curve where the prime has 254 bits:
So if we set to 62, it has been proven that a modular inverse can be computed with at most 12 iterations of the loop.
Code
transition_matrix
: Computes the matrix afterN
divsteps, this also computes
def transition_matrix(delta, u, v):
"""Compute delta and transition matrix t after N divsteps (multiplied by 2^N)."""
m00, m01, m10, m11 = 1, 0, 0, 1 # start with identity matrix
for _ in range(N):
if delta > 0 and v & 1:
delta, u, v, m00, m01, m10, m11 = 1 - delta, v, (u - v) // 2, 2*m10, 2*m11, m00 - m10, m01 - m11
elif v & 1:
delta, u, v, m00, m01, m10, m11 = 1 + delta, u, (v - u) // 2, 2*m00, 2*m01, m10 - m00, m11 - m01
else:
delta, u, v, m00, m01, m10, m11 = 1 + delta, u, v // 2, 2*m00, 2*m01, m10, m11
return delta, u, v, (m00, m01, m10, m11)
div2n
&update_x1x2
: Applies the matrix to to compute
def div2n(p, p_inv, x):
"""Compute x/2^N mod p, given p_inv = 1/p mod 2^N."""
assert (p * p_inv) % 2**N == 1
m = (x * p_inv) % 2**N
x -= m * p
assert x % 2**N == 0
return (x >> N) % p
def update_x1x2(x1, x2, t, p, p_inv):
"""Multiply matrix t/2^N with [x1, x2], modulo p."""
m00, m01, m10, m11 = t
x1n, x2n = m00*x1 + m01*x2, m10*x1 + m11*x2
return div2n(p, p_inv, x1n), div2n(p, p_inv, x2n)
Why div2n
makes divisible by :
div2n
makes divisible by :We can express an integer as:
where is the remainder mod . To make divisible by , we need to eliminate the part.
To do this, we first compute:
Then, we subtract from :
which is now cleanly divisible by .
Additionally, since , the subtraction step preserves the modular equivalence:
modinv
: Combines all components to compute the modular inverse
def modinv(p, p_inv, x):
"""Compute the modular inverse of x mod p, given p_inv=1/p mod 2^N."""
assert p & 1
delta, u, v, x1, x2 = 1, p, x, 0, 1
while v != 0:
delta, u, v, t = transition_matrix(delta, u, v)
x1, x2 = update_x1x2(x1, x2, t, p, p_inv)
assert u == 1 or u == -1 # |u| must be 1
return (x1 * u) % p
Avoiding Modulus Operations
The original modinv
function already involved a modulus operation, but with the introduction of batched divstep
, expensive modular operations now occur at every loop iteration. In this section, we aim to eliminate those costly modulus computations.
Removing Modulus Operation in div2n
div2n
Original Purpose:
The div2n
function is used to reduce results into the range .
Optimization Insight:
Instead of restricting outputs to modulo , we can eliminate the mod
operation if we ensure that all intermediate values remain within a wider range, such as .
✅ Version 1: Allowing Range Expansion
def update_x1x2_optimized_ver1(x1, x2, t, p, p_inv):
"""Multiply matrix t/2^N with [x1, x2], modulo p, given p_inv = 1/p mod 2^N."""
m00, m01, m10, m11 = t
x1n, x2n = m00*x1 + m01*x2, m10*x1 + m11*x2
# Cancel out bottom N bits of x1n and x2n.
mx1n = ((x1n * p_inv) % 2**N)
mx2n = ((x2n * p_inv) % 2**N)
x1n -= mx1n * p
x2n -= mx2n * p
return x1n >> N, x2n >> N
🔍 Range Analysis of Version 1
Each transition matrix acts on the vector in one of the following ways:
After iterations, we get:
Since initially:
In update_x1x2_optimized_ver1
, we subtract a multiple of before division:
With , we obtain:
After division by :
Hence, the result stays within .
⚠️ Drawback
After each divstep
, the range expands:
This increasing range leads to:
Unnecessary arithmetic overhead
Larger intermediate values
Potential memory inefficiencies
✅ Version 2: Pre-Clamping Strategy
To prevent range blow-up, we can add to negative values to bring the range back to :
def update_x1x2_optimized_ver2(x1, x2, t, p, p_inv):
"""Multiply matrix t/2^N with [x1, x2], modulo p, given p_inv = 1/p mod 2^N."""
m00, m01, m10, m11 = t
# x1, x2 in (-2*p, p)
if x1 < 0:
x1 += p
if x2 < 0:
x2 += p
# x1, x2 in (-p, p)
x1n, x2n = m00*x1 + m01*x2, m10*x1 + m11*x2
mx1n = -((p_inv * x1n) % 2**N)
mx2n = -((p_inv * x2n) % 2**N)
x1n += mx1n * p
x2n += mx2n * p
return x1n >> N, x2n >> N
Avoiding Modulus Operation in modinv
modinv
After the final iteration, the results lie in . We want to map into the range using the following normalization:
def normalize(sign, v, p):
"""Compute sign*v mod p, where v is in range (-2*p, p); output in [0, p)."""
assert sign == 1 or sign == -1
if v < 0:
v += p
if sign == -1:
v = -v
if v < 0:
v += p
return v
if v < 0: v += p
if sign == -1: v = -v
if v < 0: v += p
Final Implementation
def modinv(p, p_inv, x):
"""Compute the modular inverse of x mod p, given p_inv=1/p mod 2^N."""
assert p & 1
delta, u, v, x1, x2 = 1, p, x, 0, 1
while v != 0:
delta, u, v, t = transition_matrix(delta, u, v)
x1, x2 = update_x1x2_optimized_ver2(x1, x2, t, p, p_inv)
assert u == 1 or u == -1
return normalize(u, x1, p)
Conclusion
Traditional GCD-based modular inverse algorithms are fast but suffer from a major drawback: they are vulnerable to side-channel attacks due to conditional branches and comparison operations. This vulnerability can be critical in environments where sensitive key material is involved, such as RSA or ECC.
This article introduced a robust alternative: the Bernstein–Yang divstep
-based modular inverse algorithm. This approach offers several compelling advantages:
✅ Constant-time execution: Designed to ensure execution time does not vary based on input values
✅ High performance: Faster than Fermat’s method and comparable to (or better than) classical GCD-based methods
✅ Optimized computation: Improves efficiency through batch processing via matrix multiplication and by avoiding modulus operations ↪ For additional optimizations, refer to Bitcoin Core’s safegcd implementation guide
Ultimately, the divstep
-based algorithm strikes an excellent balance between performance and security. It is especially well-suited for environments that demand resistance to timing attacks, such as cryptographic protocol implementations or hardware acceleration contexts.
References
Last updated