Fast Exponentiation: Binary Exponentiation in O(log n)

Computing x to the power n by multiplying x by itself n - 1 times is fine for n = 13 and hopeless for n = 10^18. Fast exponentiation, also called binary exponentiation or exponentiation by squaring, gets there in about 2 log n multiplications instead. The idea fits in one sentence, write the exponent in binary, and it shows up far beyond interview questions: modular arithmetic, cryptography, and computing the n-th term of a recurrence without a loop of length n.

The problem it solves

Three versions of the same question come up:

• Pow(x, n): compute x^n for a float x and an integer n, which may be negative.
• Modular power: compute base^exp mod m when exp is huge, as in "return the answer modulo 10^9 + 7".
• Fast recurrences: the n-th Fibonacci number or any linear recurrence, by raising a small matrix to the n-th power.

The naive loop is O(n). With n = 10^18, that is a quintillion multiplications. Squaring brings it down to about 60 squarings and at most 60 extra multiplications, because 10^18 has 60 binary digits.

The intuition

Write the exponent in binary. 13 is 1101, which means 13 = 8 + 4 + 1, so 3^13 = 3^8 · 3^4 · 3^1.

Now build a ladder of powers where each rung is the square of the one below: 3^1, 3^2, 3^4, 3^8. Each rung costs one multiplication and doubles the exponent, so reaching 3^8 takes three multiplications rather than seven. Then keep only the rungs whose binary digit is 1 and multiply them together.

The loop does both jobs in one pass from the lowest bit up. If the current bit of exp is 1, multiply the current rung into the result. Then square the rung and shift exp right by one. When exp reaches zero, every set bit has contributed its rung.

For modular powers, reduce by the modulus after every multiplication. Remainders survive multiplication, so (a · b) mod m equals ((a mod m) · (b mod m)) mod m, which keeps every intermediate number below m without changing the answer. That is why this algorithm sits in the Big-O cheat sheet at O(log n) time and O(1) space.

Watch it run

The animation computes 3^13. It opens with the cost of the plain way: twelve multiplications, and for a crypto-sized exponent that never finishes. It writes the exponent in binary, 13 = 8 + 4 + 1, three powers of two rather than thirteen ones. Then it builds the ladder from the bottom. The first rung is the base itself, 3^1 = 3. Squaring it doubles the exponent: 3^2 = 9, one multiplication, one rung. Squaring again gives 3^4 = 81, and once more gives 3^8 = 6561, so three multiplications have reached the eighth power. Now it keeps the rungs whose bit is on. The 8 place is on, so 6561 goes into the answer. The 4 place is on too: multiply by 81, and the running answer is 3^12. The 2 place is off, so 9 is skipped. The 1 place is on, so there is one last multiply by 3. The closing frame counts it: three squarings plus two multiplies, five instead of twelve, and the count grows with log n.

The code

The lesson's modular power. The last line gives the last three digits of 2^1000 without ever building the 302-digit number:

def power(base, exp, mod):
    result = 1
    while exp:                          # one pass per binary digit
        if exp & 1:                     # this digit is on -> keep the rung
            result = result * base % mod
        base = base * base % mod        # climb: b, b², b⁴, b⁸ …
        exp >>= 1
    return result

print(power(3, 13, 10 ** 9 + 7))  # 1594323
print(power(3, 13, 1000))         # 323
print(power(2, 1000, 1000))       # 376  -> last three digits of 2**1000

Counting the multiplications, as the animation does: one squaring per rung above the first, plus one multiply for every set bit after the first. Compare with n - 1 for the naive loop:

def count_mults(exp):
    squarings = exp.bit_length() - 1        # climb to the top rung
    keeps = bin(exp).count("1") - 1         # multiply in every other set bit
    return squarings + keeps

for n in (13, 1000, 2 ** 20, 10 ** 18):
    print(n, n - 1, count_mults(n))
# 13 12 5
# 1000 999 14
# 1048576 1048575 20
# 1000000000000000000 999999999999999999 82

The float version handles a negative exponent by inverting the base first, and the recursive version squares the half power. The key is one recursive call, not two:

def my_pow(x, n):
    if n < 0:
        x, n = 1 / x, -n                # x^-n is (1/x)^n
    result = 1.0
    while n:
        if n & 1:
            result *= x
        x *= x
        n >>= 1
    return result

print(my_pow(2.0, 10), my_pow(2.0, -2), my_pow(0.5, 0))   # 1024.0 0.25 1.0

def power_rec(base, exp, mod):
    if exp == 0:
        return 1 % mod
    half = power_rec(base, exp // 2, mod)   # one recursive call, not two
    half = half * half % mod
    return half * base % mod if exp % 2 else half

print(power_rec(3, 13, 1000))     # 323

The ladder only needs an associative multiplication, so it works on matrices too. Raising [[1, 1], [1, 0]] to the n-th power gives Fibonacci numbers in O(log n) matrix products:

def mat_mult(a, b, mod):
    return [[(a[0][0] * b[0][0] + a[0][1] * b[1][0]) % mod, (a[0][0] * b[0][1] + a[0][1] * b[1][1]) % mod],
            [(a[1][0] * b[0][0] + a[1][1] * b[1][0]) % mod, (a[1][0] * b[0][1] + a[1][1] * b[1][1]) % mod]]

def fib(n, mod=10 ** 9 + 7):
    result, m = [[1, 0], [0, 1]], [[1, 1], [1, 0]]
    while n:                            # the same ladder, with matrices as the "numbers"
        if n & 1:
            result = mat_mult(result, m, mod)
        m = mat_mult(m, m, mod)
        n >>= 1
    return result[0][1]

print(fib(10), fib(90, 10 ** 20))  # 55 2880067194370816120

Everything checked against brute force on 3,000 random cases: modular powers against a loop of exp multiplications and Python's built-in three-argument pow, Fibonacci against the plain loop, and the float version against **:

import random

def slow_power(base, exp, mod):         # brute force: exp multiplications
    result = 1 % mod
    for _ in range(exp):
        result = result * base % mod
    return result

def slow_fib(n):
    a, b = 0, 1
    for _ in range(n):
        a, b = b, a + b
    return a

random.seed(22)
ok = True
for _ in range(3000):
    base, exp, mod = random.randint(0, 10 ** 6), random.randint(0, 300), random.randint(1, 10 ** 9)
    want = slow_power(base, exp, mod)
    ok &= power(base, exp, mod) % mod == want == power_rec(base, exp, mod)
    ok &= want == pow(base, exp, mod)
    n = random.randint(0, 200)
    ok &= fib(n) == slow_fib(n) % (10 ** 9 + 7)
    x, k = random.uniform(0.5, 1.5), random.randint(-40, 40)
    ok &= abs(my_pow(x, k) - x ** k) <= 1e-9 * max(1.0, abs(x ** k))
print(ok)                               # True

The complexity

• Time: O(log n) multiplications: one squaring per bit of the exponent, plus at most one extra multiply per bit.
• Space: O(1) for the loop, O(log n) stack for the recursive version.
• With big numbers: without a modulus the result itself has about n · log x bits, so each multiplication gets slower as the numbers grow. The modulus is what keeps every step cheap.

Where it goes wrong

• Two recursive calls. Writing power(b, n // 2) * power(b, n // 2) recomputes the half and turns O(log n) back into O(n).
• Reducing only at the end. In languages with fixed-size integers the intermediate products overflow long before the final % mod. Reduce after every multiplication, and make sure the product of two values below m fits in your integer type.
• A modulus of 1. Anything mod 1 is 0, but the lesson's loop returns 1 when exp is 0. Starting from 1 % mod fixes it, which is why the check above reduces the result once more.
• Negating the minimum integer. In languages with 32-bit integers, -n overflows when n is the smallest negative value. Widen the type before flipping the sign.

When it shows up in interviews

"Implement pow(x, n)" is a standard medium, and "compute this modulo 10^9 + 7" appears in counting problems where the answer is astronomically large. It pairs naturally with the sieve of Eratosthenes in number theory rounds, and with modular inverses: when m is prime, the inverse of a is a^(m-2) mod m. The bit-by-bit loop is the same reading of binary digits used in counting set bits and bit masks.

How to say it in an interview

"I read the exponent in binary. I keep a running base that I square each step, so it goes through x, x^2, x^4, x^8, and whenever the current bit of n is 1 I multiply that base into the result. Then I shift n right. That is one squaring per bit, so O(log n) multiplications and O(1) space. For a negative n I invert x first. For a modular answer I reduce after every multiplication, which keeps the numbers small without changing the result."