Modular Exponentiation, in depth: square-and-multiply traced by hand, overflow, exponent reduction, windows, Montgomery multiplication, constant time and RSA-CRT

The two facts that make it fast

Two facts make the problem tractable. First, reduction commutes with multiplication: (x * y) mod m = ((x mod m) * (y mod m)) mod m, so you may reduce after every step and no intermediate value ever exceeds (m-1)^2. Second, any exponent is a sum of powers of two, so b^e is a product of the values b^1, b^2, b^4, b^8, ... selected by the 1 bits of e, and each of those is the square of the previous one.

For e = 117 = 64 + 32 + 16 + 4 + 1, that means b^117 = b^64 * b^32 * b^16 * b^4 * b^1. Six squarings produce every power up to b^64, and four multiplications combine the five selected ones. The cost is proportional to the number of bits in the exponent, not its value: O(log e) multiplications instead of e - 1.

Two ways to walk the bits

There are two ways to walk the bits. Right-to-left keeps a running square of the base and multiplies it into the result when the current low bit is 1. Left-to-right reads from the top bit down, squares the accumulator at each step, and multiplies by the base when the bit is 1.

def modpow_rl(b, e, m):
    # Right-to-left binary exponentiation.
    if m == 1:
        return 0
    r, x = 1, b % m
    while e > 0:
        if e & 1:
            r = r * x % m
        x = x * x % m
        e >>= 1
    return r

def modpow_lr(b, e, m):
    # Left-to-right binary exponentiation.
    if m == 1:
        return 0
    r, b = 1, b % m
    for bit in bin(e)[2:]:
        r = r * r % m
        if bit == "1":
            r = r * b % m
    return r

The left-to-right form always multiplies by the same base, which is what makes the window methods below possible; the right-to-left form needs no knowledge of the exponent's length and suits streaming bits. Both handle e = 0 by returning 1, and the explicit m == 1 guard returns 0, which matches Python's built-in pow(x, 0, 1).

Worked example: 7^117 mod 23

Take 7^117 mod 23 with the left-to-right form. The exponent is 1110101 in binary. The diagram shows each step: square the accumulator, then multiply by 7 if the bit is 1.

The right-to-left run reaches the same answer by a different path. The running square x takes the values 7, 3, 9, 12, 6, 13, 8 (that is, 7^1, 7^2, 7^4, ... 7^64 mod 23), and the result multiplies in the ones at bits 0, 2, 4, 5 and 6, reading 7, 7, 17, 17, 10, 15, 5 after each bit. Both give 5.

A second check uses Fermat's little theorem, covered next: 23 is prime, so 7^22 = 1 mod 23, and 117 mod 22 = 7. Then 7^117 = 7^7 mod 23, and 7^7 = 823543 = 23 * 35806 + 5. Writing such cross-checks into tests catches most implementation bugs.

Counting the cost

For a k-bit exponent with w one bits, the left-to-right method starting from r = b at the top bit costs k - 1 squarings and w - 1 multiplications. A random exponent has about k/2 one bits, so the average is roughly 1.5k modular multiplications; for a 2048-bit RSA exponent that is around 3,000.

The cost of each multiplication depends on the size of m. When m fits a machine word it is a constant; for multi-precision numbers of n words, schoolbook multiplication costs about n^2 word operations plus a similar reduction. A full exponentiation with a modulus and exponent of k bits therefore costs on the order of k^3 bit operations with schoolbook arithmetic. That cubic shape is why RSA private-key operations are expensive and why the CRT trick below works.

Overflow in fixed-width integers

In languages with fixed-width integers, the product r * x needs twice the bits of m. With unsigned 64-bit arithmetic, the plain expression is only safe while m is below 2^32. The classic competitive-programming modulus 10^9 + 7 is safe; a 62-bit modulus used in Miller-Rabin is not. Widen the product instead:

#include <stdint.h>
typedef unsigned __int128 u128;            /* GCC and Clang */

static uint64_t mulmod(uint64_t a, uint64_t b, uint64_t m) {
    return (uint64_t)((u128)a * b % m);
}

uint64_t powmod(uint64_t b, uint64_t e, uint64_t m) {
    if (m == 1) return 0;
    uint64_t r = 1;
    b %= m;
    while (e) {
        if (e & 1) r = mulmod(r, b, m);
        b = mulmod(b, b, m);
        e >>= 1;
    }
    return r;
}

In Java, use Math.multiplyHigh or BigInteger.modPow; in JavaScript, Number loses exactness above 2^53, so switch to BigInt once m exceeds about 2^26. The silent version of this bug returns plausible wrong answers, which is why the deterministic Miller-Rabin implementation in the Miller-Rabin article uses a widened multiply.

Reducing the exponent, and when you may not

Large exponents can often be shrunk first. If p is prime and does not divide b, Fermat's little theorem gives b^(p-1) = 1 mod p, so e can be reduced modulo p - 1. Euler's theorem generalises this to b^phi(m) = 1 mod m, but only when gcd(b, m) = 1. Forgetting the condition gives wrong answers: with m = 12 and phi(12) = 4, the true value of 2^4 mod 12 is 4, while reducing the exponent to 0 gives 1.

When the base and modulus share factors, a weaker rule still holds: for e at least log2(m), b^e = b^((e mod phi(m)) + phi(m)) mod m. It is what makes power towers like a^(b^c) mod m computable by recursing on phi. Background on phi, inverses and the CRT is in number theory for programmers.

The same theorem gives modular inverses: for prime p, b^(p-2) mod p is the inverse of b. Combinatorics code computes nCr mod p this way from factorial tables; Lucas's theorem extends it to n larger than p. Python 3.8 and later also accept pow(b, -1, m) for inverses modulo any m coprime to b.

Going faster: windows and Montgomery multiplication

Sliding windows. Left-to-right processing can consume several bits at a time. Precompute the odd powers b^1, b^3, ..., b^(2^w - 1); then, scanning the exponent, square once per bit and multiply once per window of up to w bits that starts and ends with a 1. Squarings stay at about k, but multiplications drop from about k/2 to roughly k/(w+1) plus the 2^(w-1) table entries. For exponents of a few thousand bits, windows of around five bits are a common choice; the best value depends on the exponent size.

Montgomery multiplication. The expensive part of each step is the division in mod m. For an odd modulus, Montgomery's method picks R = 2^k greater than m and keeps numbers in the form aR mod m. Multiplying two such numbers and applying REDC, which divides by R using only multiplications, masks and shifts, yields another number in the same form. Conversion in and out costs a few operations, so it pays off over the thousands of multiplications in one exponentiation.

def montgomery_setup(m, k):
    R = 1 << k                      # R > m, gcd(R, m) = 1 because m is odd
    m_neg_inv = (-pow(m, -1, R)) % R
    return R, m_neg_inv

def redc(T, m, k, m_neg_inv):
    # Return T * R^-1 mod m for 0 <= T < m*R, with no division by m.
    mask = (1 << k) - 1
    u = ((T & mask) * m_neg_inv) & mask
    t = (T + u * m) >> k
    return t - m if t >= m else t

def modpow_mont(b, e, m, k):
    R, mi = montgomery_setup(m, k)
    x = b * R % m                   # into Montgomery form
    r = R % m                       # 1 in Montgomery form
    for bit in bin(e)[2:]:
        r = redc(r * r, m, k, mi)
        if bit == "1":
            r = redc(r * x, m, k, mi)
    return redc(r, m, k, mi)        # back out of Montgomery form

Real libraries implement REDC word by word on multi-precision limbs; the Python above shows the arithmetic, not the speed.

Secret exponents and constant time

When the exponent is a secret, such as an RSA private exponent, square-and-multiply leaks it. The multiply happens only for 1 bits, so execution time, power draw and cache behaviour correlate with the key. Paul Kocher's 1996 timing-attack paper showed this is practical. Variable-time code is therefore wrong for secrets even when it is mathematically correct.

The Montgomery ladder performs the same operations for every bit: one multiplication and one squaring, with the roles of two registers swapped by the bit value.

def ladder(b, e, m, bits):
    r0, r1 = 1, b % m
    for i in reversed(range(bits)):     # fixed length, independent of the key's value
        bit = (e >> i) & 1
        r0, r1 = cswap(r0, r1, bit)
        r1 = r0 * r1 % m
        r0 = r0 * r0 % m
        r0, r1 = cswap(r0, r1, bit)
    return r0

The ladder alone is not enough: cswap must be a branch-free masked swap, the multiplication itself must be constant-time, and table lookups in windowed versions must not index memory by secret bits. Python integers, pow and Java's BigInteger.modPow are not designed for this, so for real keys use a vetted cryptographic library and treat hand-rolled code as teaching material.

RSA and the Chinese remainder theorem

RSA decryption computes c^d mod n with n = p * q. Knowing the factors, you can instead compute two exponentiations modulo the half-size primes and combine them with the Chinese remainder theorem. Under the schoolbook cost model, an exponentiation with k-bit numbers costs about k^3, so two of size k/2 cost 2 * (k/2)^3 = k^3 / 4: roughly four times less work, before measuring any particular library.

A toy key shows the steps. With p = 61, q = 53, n = 3233, e = 17 and d = 2753, the message 65 encrypts to c = 2790. The CRT parameters are dp = d mod 60 = 53, dq = d mod 52 = 49 and qinv = 38. Then m1 = 2790^53 mod 61 = 4, m2 = 2790^49 mod 53 = 12, h = 38 * (4 - 12) mod 61 = 1, and m = 12 + 1 * 53 = 65.

CRT has a sharp edge. If a hardware or software fault corrupts one of the two halves of a signature, the faulty signature reveals a factor of n through a single gcd, as work by Boneh, DeMillo and Lipton and by Lenstra showed. Implementations verify the signature with the public exponent before releasing it.

Failure modes

What to do next

1. Implement both binary forms and test them against your language's built-in on random inputs, including e = 0, m = 1, negative bases and moduli near your word size.
2. Audit every hand-written modular multiply for overflow and widen it where m can exceed 2^32.
3. Wherever you reduce exponents, write down why the base and modulus are coprime, or switch to the extended rule.
4. Read the Miller-Rabin and Legendre and Jacobi symbol articles; both are built on this primitive.
5. For anything involving secret exponents, find where your system calls a cryptographic library and confirm no hand-rolled or general-purpose big-integer exponentiation touches key material.
6. If performance matters, profile before tuning: confirm exponentiation dominates, then try a windowed, Montgomery-based library routine rather than writing your own.