Fast Modular Exponentiation Using Binary Exponentiation in Python

When computing \(a^k \pmod{p}\) for large numbers, naively multiplying \(a\) by itself \(k\) times before taking the remainder is extremely inefficient, as intermediate values become astronomically large.

Binary exponentiation (also known as exponentiation by squaring) is an algorithm that performs this computation efficiently. It is an essential technique for calculating modular exponentiation of large numbers in applications such as RSA encryption.

The key insight is to express the exponent \(k\) in binary, reducing the computational complexity to \(O(\log k)\) .

Example

Consider computing \(5^{21} \pmod{p}\) .

First, express the exponent \(21\) in binary: \(21 = 16 + 4 + 1 = 1 \cdot 2^4 + 0 \cdot 2^3 + 1 \cdot 2^2 + 0 \cdot 2^1 + 1 \cdot 2^0\) So the binary representation of \(21\) is \((10101)_2\) .

Using this, \(5^{21}\) can be decomposed as:

By sequentially computing powers of the form \(a^{2^i}\) (simply by squaring the previous result), we only need to multiply together the terms where the binary digit is 1. This dramatically reduces the number of multiplications.

Implementation

Binary exponentiation comes in two flavors: reading the exponent’s binary representation from the least significant bit (Right-to-Left) or from the most significant bit (Left-to-Right). Both achieve the same \(O(\log k)\) complexity — rather than stating this as a given, we derive it from the binary expansion of the exponent.

Why \(O(\log k)\) multiplications suffice

Expanding the exponent \(k\) in binary, with bit length \(L = \lfloor \log_2 k \rfloor + 1\) , gives

By the law of exponents \(a^{x+y} = a^x a^y\) ,

The factors \(a^{2^0}, a^{2^1}, \dots, a^{2^{L-1}}\) on the right-hand side can all be obtained starting from \(a^{2^0}=a\) by squaring the previous value:

which takes only \(L-1\) squarings. Multiplying together only the factors whose bit \(b_i = 1\) then yields \(a^k\) , so the total number of multiplications is at most “\(L-1\) squarings + \(L\) or fewer products,” i.e. \(O(L) = O(\log k)\) . This bound does not depend on which direction the bits are scanned — the underlying structure (“prepare the powers \(a^{2^i}\) , multiply in only where the bit is set”) is direction-agnostic. Below, we turn this structure into two concrete algorithms — scanning from the low end and from the high end — and prove correctness for each via an explicit loop invariant.

Right-to-Left (scanning from the least significant bit)

The Right-to-Left method reads the exponent’s bits starting from the LSB: whenever a bit is 1, it multiplies res by the current base (which holds \(a^{2^i}\) at that point), then squares base and shifts exp right by one bit.

Loop invariant: immediately before the \(i\) -th iteration of the while loop (\(i = 0, 1, 2, \dots\) ):

Proof (induction): at \(i=0\) , res\(=1=a^0\) , base\(=a=a^{2^0}\) , exp\(=k\) , so the invariant holds trivially. Assuming it holds at \(i\) , the iteration reads the least significant bit \(b_i = \texttt{exp} \bmod 2\) of exp, and multiplies res *= base only when \(b_i=1\) , so the updated res becomes

(using the identity \(k \bmod 2^{i+1} = (k \bmod 2^i) + b_i 2^i\) from the binary expansion). The loop then squares base (\(a^{2^i} \to a^{2^{i+1}}\) ) and right-shifts exp (\(\lfloor k/2^i \rfloor \to \lfloor k/2^{i+1} \rfloor\) ), so the invariant holds at \(i+1\) as well. The loop terminates when exp\(=0\) , i.e. at \(i=L\) ; since \(k < 2^L\) , \(k \bmod 2^L = k\) , so res\(=a^k\) . \(\blacksquare\)

Python Program

def power(base, exp, mod):
    """
    Efficiently calculates (base^exp) % mod using binary exponentiation.

    :param base: The base number.
    :param exp: The exponent.
    :param mod: The modulus.
    :return: The result of (base^exp) % mod.
    """
    res = 1
    base %= mod
    while exp > 0:
        # If the least significant bit of exp is 1, multiply res by base
        if exp % 2 == 1:
            res = (res * base) % mod

        # Square the base and right-shift the exponent by 1 bit
        base = (base * base) % mod
        exp //= 2

    return res

# --- Examples ---
# Compute 5^21 mod 99
k = 21
g = 5
p = 99
result = power(g, k, p)
print(f"{g}^{k} mod {p} = {result}") # -> 5^21 mod 99 = 20

# Example with very large numbers
k = 12345678901234567890
g = 987654321987654321
p = 1000000007
result = power(g, k, p)
print(f"Result for large numbers: {result}")

How the Program Works

1. res = 1: Initializes the result variable to 1.
2. while exp > 0:: Loops until the exponent exp becomes 0.
3. if exp % 2 == 1:: Checks whether the least significant bit of exp is 1.
4. res = (res * base) % mod: If the bit is 1, multiplies the current base into the result res.
5. base = (base * base) % mod: Squares base, progressively computing \(a, a^2, a^4, a^8, \dots\) .
6. exp //= 2: Integer-divides exp by 2, shifting to the next bit.
7. When the loop ends, res holds the final result.

Since the remainder is taken at each step, intermediate values never exceed the modulus, preventing overflow and enabling fast computation even with extremely large numbers (see “Edge Cases and Practical Notes” below for details).

Left-to-Right (scanning from the most significant bit)

The Left-to-Right method reads the exponent’s bits starting from the MSB: it squares res on every step, and multiplies in base whenever the current bit is 1. This is exactly Horner’s method, the same scheme used to evaluate a binary number as a decimal integer.

Let \(p_j = \lfloor k / 2^{L-j} \rfloor\) denote the value formed by the top \(j\) bits read so far (\(j = 0, \dots, L\) , with \(p_0 = 0\) and \(p_L = k\) ). The structure of binary notation gives the recurrence

Loop invariant: immediately before the \(j\) -th iteration (\(j=0,1,\dots,L-1\) ), res\(\equiv a^{p_j} \pmod m\) .

Proof: at \(j=0\) , res\(=1=a^{p_0}=a^0\) , so the invariant holds. Assuming it holds at \(j\) , the iteration first squares res (\(a^{p_j} \to a^{2p_j}\) ), then multiplies in base (\(=a\) ) if bit \(b_{L-1-j}\) is 1, so the updated res is

so the invariant holds at \(j+1\) . When the loop finishes at \(j=L\) , res\(=a^{p_L}=a^k\) . \(\blacksquare\)

def power_left_to_right(base, exp, mod):
    """
    Left-to-Right binary exponentiation, scanning from the most
    significant bit of the exponent.
    """
    if exp == 0:
        return 1 % mod
    base %= mod
    res = 1
    for bit in bin(exp)[2:]:  # binary string without the "0b" prefix, MSB first
        res = (res * res) % mod
        if bit == "1":
            res = (res * base) % mod
    return res

Both Right-to-Left and Left-to-Right perform exactly \(L-1\) squarings, and a multiplication for every bit that is set — that is, a number of multiplications equal to the exponent’s Hamming weight \(w(k) = \sum_i b_i\) . The total multiplication count therefore falls in the range \(L-1 \le (\text{multiplications}) \le 2L-1\) regardless of scan direction, confirming \(\Theta(\log k)\) either way. In practice, Right-to-Left is simpler to implement because squaring base and right-shifting exp are independent update steps, while Left-to-Right pairs naturally with implementations that treat the exponent as a bit string or bit array — including the Montgomery ladder discussed next.

Side-Channel Attacks and Constant-Time Exponentiation: the Montgomery Ladder

Both methods above share a weakness: the branch if bit == 1: multiply is a data-dependent branch. When the exponent is public (e.g. RSA’s encryption exponent \(e\) ), this is harmless. But when the exponent is secret — RSA’s private exponent \(d\) (see the RSA article ) or a secret scalar in elliptic-curve scalar multiplication (see the Elliptic Curve Cryptography article ’s double-and-add) — this branch can leak the secret through observable physical side effects such as execution time, power consumption, or electromagnetic emission (side channels).

The mechanics of timing attacks

For a fixed exponent bit length \(L\) , Right-to-Left always performs \(L-1\) squarings, but the number of multiplications equals the Hamming weight \(w(k)\) . If an attacker can measure execution time across many operations that reuse the same secret exponent (e.g. repeated decryption requests with a fixed RSA private key \(d\) ), the statistics reveal “more multiplications happened here” — i.e. higher Hamming weight — and with a more refined measurement setup, an attacker can recover the secret bit by bit. This is exactly the class of attack Paul Kocher demonstrated in 1996 against naive RSA/DH/DSS implementations, after which constant-time exponentiation became the standard defense in RSA/ECC implementations. The “Verification” section below measures this effect directly, showing that a difference in multiplication count translates into a measurable difference in wall-clock time.

Deriving the Montgomery ladder

The fix is not “eliminate the branch,” but rather “make the number and kind of operations performed independent of the secret bit’s value.” The Montgomery ladder achieves this by always maintaining two registers \(R_0, R_1\) and performing exactly one multiplication and one squaring on every single bit — the only thing that changes with the bit’s value is which register receives which result.

Using the same prefix \(p_j = \lfloor k/2^{L-j} \rfloor\) as in Left-to-Right, impose the loop invariant

i.e. “\(R_1\) is always exactly one power ahead of \(R_0\) .” When the next bit \(b_{L-1-j}\) advances the prefix to \(p_{j+1} = 2p_j + b_{L-1-j}\) :

• If \(b_{L-1-j}=0\) : \(p_{j+1}=2p_j\) . We need \(R_0' = a^{2p_j} = R_0^2\) and \(R_1' = a^{2p_j+1} = a^{p_j}\cdot a^{p_j+1} = R_0 \cdot R_1\) .
• If \(b_{L-1-j}=1\) : \(p_{j+1}=2p_j+1\) . We need \(R_0' = a^{2p_j+1} = R_0 \cdot R_1\) and \(R_1' = a^{2p_j+2} = (a^{p_j+1})^2 = R_1^2\) .

In both cases the work is identical: “one product \(R_0 \cdot R_1\) ” and “one squaring of one of the registers” — the operation count and type are exactly the same; only the destination of the product and the destination of the square differ. The initial state at \(j=0\) (\(p_0=0\) ) is \(R_0=a^0=1\) , \(R_1=a^1=a\) , and after processing all \(L\) bits (\(j=L\) , \(p_L=k\) ), \(R_0 = a^k\) is returned.

def power_montgomery_ladder(base, exp, mod):
    """
    Montgomery-ladder exponentiation: the operation count is independent
    of the exponent's bit pattern. Every bit does exactly one
    multiplication and one squaring; only the assignment target changes.
    """
    if mod == 1:
        return 0
    base %= mod
    if exp == 0:
        return 1 % mod

    bitlen = exp.bit_length()
    R0, R1 = 1, base
    for i in range(bitlen - 1, -1, -1):
        bit = (exp >> i) & 1
        if bit == 0:
            R1 = (R0 * R1) % mod
            R0 = (R0 * R0) % mod
        else:
            R0 = (R0 * R1) % mod
            R1 = (R1 * R1) % mod
    return R0

One important caveat: the if bit == 0: ... else: ... branch itself is still a Python-level branch, and can leak information through CPU branch prediction or cache timing on real hardware. Truly constant-time implementations (such as the Montgomery ladder over Curve25519 used by X25519) replace the branch with a bitwise conditional swap (cswap) — using a bitmask such as mask = -bit to swap without branching, e.g. R0, R1 = R0 ^ (mask & (R0 ^ R1)), R1 ^ (mask & (R0 ^ R1)). What this article’s Montgomery ladder solves is specifically the problem that the number and kind of operations depend on the secret bit’s value (i.e. the fact that naive binary exponentiation performs a number of multiplications proportional to the Hamming weight) — which is the dominant contributor to the timing difference measured in the next section. In production RSA/ECC implementations of secret-key operations, use libraries hardened down to the cswap level, such as cryptography or OpenSSL.

Edge Cases and Practical Notes

• Exponent equal to zero: \(a^0 = 1\) (by convention \(0^0 = 1\) as well; indeed Python’s pow(0, 0, m) returns 1 % m). Right-to-Left never executes the loop body when exp=0, so res=1 is returned correctly without any special-casing. Left-to-Right and the Montgomery ladder handle this case with an explicit early return, which is required for the latter since bin(0) / bit_length() produce no bits to iterate over.

• Negative exponents: computing \(a^{-n} \bmod m\) requires the modular inverse \(a^{-1}\) of \(a\) with respect to \(m\) (satisfying \(a \cdot a^{-1} \equiv 1 \pmod m\) , which exists only when \(\gcd(a,m)=1\) ), found via the extended Euclidean algorithm, and then \(a^{-n} \equiv (a^{-1})^n \pmod m\) — reducing the problem back to ordinary binary exponentiation. This is exactly the idea covered in the modular inverse section of the RSA article . Python’s three-argument pow supports negative exponents directly: pow(3, -1, 11) returns 4 (verified; \(3 \times 4 = 12 \equiv 1 \pmod{11}\) ), computed via the extended Euclidean algorithm. When \(\gcd(a,m) \ne 1\) no inverse exists — for instance pow(2, -1, 4) raises ValueError: base is not invertible for the given modulus (also verified).

• Modulus equal to one: \(a^k \bmod 1\) is always \(0\) . The Montgomery ladder implementation special-cases mod == 1 explicitly, though even without it base %= mod would force base to 0, giving the same result. Handling it explicitly makes the boundary condition self-documenting.

• Avoiding overflow: naively computing \(a^k\) first and only reducing modulo \(m\) afterward is infeasible, since \(a^k\) has roughly \(k \log_{10} a\) digits — for RSA-2048-scale \(a, k\) (hundreds to 2048 bits), this is astronomically large. As in this article’s code, reducing modulo \(m\) at every step keeps intermediate values within \([0, m)\) , so the computation runs in \(O(\log k)\) rounds of “fixed-size multiplication + reduction” instead. Forgetting to reduce modulo \(m\) in an RSA or DH implementation does not change correctness, but degrades running time and memory usage by orders of magnitude.

Verification: Correctness and Timing Measurements

Correctness

We verified that the Right-to-Left, Left-to-Right, and Montgomery-ladder implementations always agree with Python’s built-in pow(base, exp, mod). Testing 500 random cases with bit lengths in \(\{8, 16, 64, 256, 512\}\) , plus 8 edge cases (\(k=0\) , \(0^0\) , \(a=0\) , \(k=1\) , \(m=1\) , and a 20-digit exponent), all 508 test cases produced identical output across all four implementations (pow, Right-to-Left, Left-to-Right, Montgomery ladder), with zero mismatches. (Negative-exponent / modular-inverse behavior was checked separately, as described above.)

import random

random.seed(42)
cases = []
for _ in range(500):
    bits = random.choice([8, 16, 64, 256, 512])
    mod = random.randrange(1, 2**bits)
    base = random.randrange(0, 2**bits)
    exp = random.randrange(0, 2**bits)
    cases.append((base, exp, mod))
cases += [(5, 0, 99), (0, 0, 99), (0, 5, 99), (5, 1, 99), (5, 21, 1),
          (12345678901234567890, 1, 1000000007), (2, 3, 5),
          (987654321987654321, 12345678901234567890, 1000000007)]

mismatches = [
    (b, e, m) for b, e, m in cases
    if not (pow(b, e, m) == power(b, e, m)
            == power_left_to_right(b, e, m)
            == power_montgomery_ladder(b, e, m))
]
print(len(cases), len(mismatches))  # -> 508 0

Timing: dependence on the exponent’s bit pattern

We measured how execution time for Right-to-Left (ordinary binary exponentiation) and the Montgomery ladder responds to the exponent’s bit pattern. Using an RSA-2048-scale 2048-bit modulus, we constructed three exponents that all share the same 2048-bit length but very different Hamming weights:

• all_ones: every bit set (Hamming weight 2048, \(2^{2048}-1\) )
• sparse: only the most- and least-significant bits set (Hamming weight 2)
• alternating: alternating 1s and 0s (Hamming weight 1024)

Each condition was run for 150 trials, and — crucially — the three conditions were interleaved round-robin rather than measured in three separate blocks, so that thermal drift or OS jitter would affect every condition equally rather than biasing whichever condition happened to run first or last (running all trials for one condition before moving to the next introduces a systematic slow-down bias from CPU throttling over time).

Measuring the spread between conditions as (max median − min median) / min median: Right-to-Left shows a 121.6% spread (sparse at 27.22 ms vs. all_ones at 60.31 ms — more than doubling in time purely because the exponent has 1024 more set bits), while the Montgomery ladder shows only a 0.1% spread (all three conditions cluster around 59.7–59.8 ms), i.e. execution time is essentially independent of the bit pattern. This matches the theory from the previous section quantitatively: Right-to-Left’s multiplication count scales with Hamming weight, while the Montgomery ladder always performs a fixed number of operations. The measurement confirms directly, in real numbers, the premise underlying timing attacks — that a secret exponent with higher Hamming weight takes measurably longer to process.

import time, random, statistics

def make_exponents(bitlen):
    all_ones = (1 << bitlen) - 1
    sparse = (1 << (bitlen - 1)) | 1
    alternating = int("10" * (bitlen // 2), 2) | (1 << (bitlen - 1))
    return {"all_ones": all_ones, "sparse": sparse, "alternating": alternating}

def time_interleaved(func, base, exps, mod, repeats):
    samples = {label: [] for label in exps}
    for _ in range(repeats):
        for label, exp in exps.items():          # round-robin
            t0 = time.perf_counter()
            func(base, exp, mod)
            samples[label].append(time.perf_counter() - t0)
    return {l: (statistics.median(v), statistics.stdev(v)) for l, v in samples.items()}

random.seed(7)
bitlen = 2048
mod = random.getrandbits(bitlen) | 1 | (1 << (bitlen - 1))
base = random.getrandbits(bitlen) % mod
exps = make_exponents(bitlen)

r2l = time_interleaved(power, base, exps, mod, repeats=150)
mont = time_interleaved(power_montgomery_ladder, base, exps, mod, repeats=150)

Application in RSA Encryption

One of the most important applications of binary exponentiation is RSA encryption . RSA requires modular exponentiation of large integers for both encryption and decryption.

Here, \(e\) and \(d\) are integers hundreds to thousands of bits long, making computation infeasible without binary exponentiation. In particular, the private exponent \(d\) used during decryption is genuinely secret, making it a target for the side-channel attacks described above. Production RSA implementations (OpenSSL, the cryptography library, etc.) avoid naive binary exponentiation, instead combining the Montgomery ladder or Montgomery multiplication (a related technique that speeds up the modular reduction itself, avoiding division) with constant-time conditional selection so that \(d\) ’s bit pattern never leaks through timing.

Computational Complexity Comparison

Python’s built-in pow(base, exp, mod) function uses binary exponentiation internally and can be used directly in RSA implementations. That said, CPython’s pow does not guarantee constant-time behavior in every case, so production systems handling secret keys should prefer libraries such as cryptography that delegate to hardened implementations in OpenSSL.

Related Articles

• RSA Encryption: Theory and Python Implementation - Explains how modular exponentiation serves as the core operation in RSA encryption and decryption.
• Elliptic Curve ElGamal Encryption: Principles and Simplified Python Implementation - A public-key cryptosystem that applies modular exponentiation to scalar multiplication on elliptic curves.
• Diffie-Hellman Key Exchange Protocol: Theory and Python Implementation - Covers the Diffie-Hellman key exchange protocol based on modular exponentiation.
• Security Certification Comparison: CISSP vs Japan’s Registered Information Security Specialist - Compares security certifications whose exam scope includes cryptographic technologies such as RSA.
• Cryptography Roadmap: Symmetric, Public-Key, Hash, and Key Exchange Hub — A meta-hub that maps RSA, DH, and ElGamal, all of which use binary exponentiation as their core compute kernel.
• RSA Encryption: Theory and Python Implementation — Canonical application: \(c = m^e \bmod n\) is infeasible without binary exponentiation.
• Diffie-Hellman Key Exchange: Theory and Python Implementation — Fast \(g^a \bmod p\) via binary exponentiation is what makes DH practical, and is also the baseline assumed in DH’s discrete-log security analysis. Also covers X25519’s constant-time Montgomery ladder.
• Elliptic Curve Cryptography (ECC) in Python — The square-and-multiply structure of binary exponentiation maps directly onto double-and-add scalar multiplication on elliptic curves.