Fast Exponentiation in Python: Computing Huge Powers in log n Steps

Ask Python for 2 ** 10**18 % 1000000007 and it will never answer. Not because the modulus is hard, but because ** is evaluated first, and 2 ** 10**18 is a number with 10^18 binary digits. Writing it down needs about 125 petabytes of memory. The machine will thrash and die long before the % runs.

Ask Python for pow(2, 10**18, 1000000007) and it answers instantly. Same value, 83 multiplications, every intermediate number under ten digits. The gap between those two lines is one of the largest practical speedups in all of elementary computing, and it comes from a single identity you already know: a^n for even n equals (a^(n/2))^2.

Halve the exponent, square the answer. Repeat. An exponent of 10^18 reaches zero in 60 halvings, not 10^18 decrements. That is fast exponentiation — also called binary exponentiation, exponentiation by squaring, or square-and-multiply — and it is the arithmetic that every RSA handshake and every Diffie-Hellman key exchange on the internet is built from.

The idea

Start with the obvious way to compute a^n: multiply by a, n times.

def naive_power(base: int, exponent: int) -> int:
    """Multiply base by itself, one multiplication per unit of the exponent."""
    result = 1
    for _ in range(exponent):
        result *= base
    return result


print(naive_power(3, 13))
print(naive_power(2, 10))

1594323
1024

That costs n multiplications. The cost is proportional to the value of the exponent, which is exponential in the number of digits you typed. A 19-digit exponent means a quintillion multiplications.

Now the identity. For any a and any even n:

a^n = (a^(n/2))^2

Read it right to left and it is a bargain. To get a^n, you do not need a^n. You need a^(n/2), and then one squaring. Half the problem, plus one multiplication.

Odd exponents cannot be halved cleanly, so peel one factor off first:

a^n = a × a^(n-1) when n is odd

Subtracting 1 from an odd number always leaves an even number, so an odd step is immediately followed by a halving. The two rules together, plus the base case a^0 = 1, define the whole algorithm:

• n is 0: the answer is 1.
• n is even: compute a^(n/2), square it.
• n is odd: compute a^(n-1), multiply by a.

Every two steps at worst, the exponent halves. Halving 10^18 takes 60 steps. That is the entire performance argument.

Watching it work

Take 3^13. The naive loop would do 12 multiplications. Here is what the halving rules do instead, working downwards from the top:

• 3^13 — 13 is odd, so 3^13 = 3 × 3^12. Need 3^12.
• 3^12 — even, so 3^12 = (3^6)^2. Need 3^6.
• 3^6 — even, so 3^6 = (3^3)^2. Need 3^3.
• 3^3 — odd, so 3^3 = 3 × 3^2. Need 3^2.
• 3^2 — even, so 3^2 = (3^1)^2. Need 3^1.
• 3^1 — odd, so 3^1 = 3 × 3^0. Need 3^0.
• 3^0 — the base case, 1.

Now unwind, filling in real numbers on the way back up: 3^0 = 1, 3^1 = 3, 3^2 = 9, 3^3 = 27, 3^6 = 729, 3^12 = 531441, 3^13 = 1594323.

Six lines, and only five of them are genuine multiplications — 3 × 3^0 is a multiplication by 1, which costs nothing in principle. Twelve multiplications became five.

The same thing, seen in binary

There is a second way to look at this that explains the name square-and-multiply, and it is the one the iterative code implements.

Write the exponent in binary. 13 is 1101, which says 13 = 8 + 4 + 1. Exponents add when powers multiply, so:

3^13 = 3^8 × 3^4 × 3^1

You can build 3^1, 3^2, 3^4, 3^8 by starting at 3 and squaring repeatedly — each squaring doubles the exponent, so four squarings reach 3^16. Then keep only the ones whose bit is set and multiply them together.

6561 × 81 × 3 = 1594323. Two multiplications to combine them, three squarings to build them: five, matching the recursion exactly. That is not a coincidence — the recursion's even steps are the squarings and its odd steps are the kept bits, read in the opposite order.

The number of squarings is the number of bits in n minus one. The number of combining multiplications is the number of 1 bits minus one. Both are at most log2(n), which is where the bound comes from.

The code

The recursive version reads exactly like the three rules.

def fast_power(base: int, exponent: int) -> int:
    """base ** exponent, computed with the halving identity.

    An even exponent squares the half power. An odd exponent peels off one
    factor of base, which always leaves an even exponent behind.
    """
    if exponent == 0:
        return 1
    if exponent % 2 == 0:
        half = fast_power(base, exponent // 2)
        return half * half
    return base * fast_power(base, exponent - 1)


print(fast_power(3, 13))
print(fast_power(2, 10))
print(fast_power(5, 0))

1594323
1024
1

Add a print and you can watch the walkthrough happen for real:

def traced_power(base: int, exponent: int) -> int:
    """The same recursion, printing each power as the calls unwind."""
    if exponent == 0:
        print(f"{base}^0 = 1")
        return 1
    if exponent % 2 == 0:
        half = traced_power(base, exponent // 2)
        print(f"{base}^{exponent} = ({base}^{exponent // 2})^2 = {half * half}")
        return half * half
    smaller = traced_power(base, exponent - 1)
    print(f"{base}^{exponent} = {base} * {base}^{exponent - 1} = {base * smaller}")
    return base * smaller


traced_power(3, 13)

3^0 = 1
3^1 = 3 * 3^0 = 3
3^2 = (3^1)^2 = 9
3^3 = 3 * 3^2 = 27
3^6 = (3^3)^2 = 729
3^12 = (3^6)^2 = 531441
3^13 = 3 * 3^12 = 1594323

The iterative version walks the exponent's bits from the bottom up. It keeps two values: result, the answer accumulated so far, and current, which holds base raised to the next power of two.

def fast_power_iterative(base: int, exponent: int, trace: bool = False) -> int:
    """Square-and-multiply, driven by the bits of the exponent from the bottom up."""
    result = 1
    current = base       # after step k this holds base ** (2 ** k)
    remaining = exponent

    while remaining > 0:
        bit = remaining % 2
        if bit == 1:
            result *= current
        if trace:
            print(f"remaining={remaining:<3} bit={bit}  current={current:<5} result={result}")
        remaining //= 2
        if remaining > 0:
            # The last squaring would never be used, so skip it.
            current *= current

    return result


print(fast_power_iterative(3, 13))
fast_power_iterative(3, 13, trace=True)

1594323
remaining=13  bit=1  current=3     result=3
remaining=6   bit=0  current=9     result=3
remaining=3   bit=1  current=81    result=243
remaining=1   bit=1  current=6561  result=1594323

How the code maps to the idea

remaining % 2 is the current low bit of the exponent, and remaining //= 2 is the halving — a right shift by one place. So the loop runs once per bit of n, which for 10^18 is 60 times.

current is the repeated-squaring chain. It starts at base and squares every iteration, so it passes through base^1, base^2, base^4, base^8 and so on. Iteration k always has current equal to base ** (2 ** k), which is exactly the value the binary decomposition wants when bit k is set.

result *= current is the "keep this one" decision. It fires only when the bit is 1, which is why the loop's multiplication count depends on how many 1 bits n has.

The order inside the loop matters. You must fold current into result before squaring current, because the bit you just examined refers to the current value, not the next one. Swapping those two lines silently produces base ** (2 * n) instead.

The final squaring is skipped by the if remaining > 0 guard. Nothing reads current after the last iteration, and that skipped squaring is the most expensive one, on the biggest numbers.

The recursion's depth is logarithmic, not linear. Each odd step is followed by an even step, so fast_power recurses at most 2 × log2(n) deep — about 120 frames for a 10^18 exponent, nowhere near Python's default limit of 1000.

Counting the multiplications directly makes the growth impossible to argue with:

def multiplication_count(exponent: int) -> int:
    """Multiplications square-and-multiply performs, without performing them."""
    count = 0
    remaining = exponent
    while remaining > 0:
        if remaining % 2 == 1:
            count += 1          # fold base ** (2 ** k) into the result
        remaining //= 2
        if remaining > 0:
            count += 1          # square to reach the next power of two
    return count


print(f"{'n':>19}  {'bits in n':>9}  {'naive':>19}  {'fast':>4}")
for exponent in (13, 100, 1000, 10 ** 6, 10 ** 18):
    print(f"{exponent:>19}  {exponent.bit_length():>9}  {exponent - 1:>19}  "
          f"{multiplication_count(exponent):>4}")

                  n  bits in n                naive  fast
                 13          4                   12     6
                100          7                   99     9
               1000         10                  999    15
            1000000         20               999999    26
1000000000000000000         60   999999999999999999    83

A quintillion multiplications become 83, and fast never exceeds twice bits in n.

Modular exponentiation

Everything above still produces enormous numbers: 7^1000 looks harmless and has 846 digits. In cryptography you never want a^n itself — you want a^n mod m, and there the numbers can be kept small at every step.

The rule that makes this work is that modular arithmetic survives multiplication:

(x × y) mod m = ((x mod m) × (y mod m)) mod m

So you can reduce after every single multiplication instead of at the end. No intermediate value ever exceeds m^2, which for a 2048-bit modulus means numbers of at most 4096 bits rather than numbers of n × 2048 bits.

def fast_power_mod(base: int, exponent: int, modulus: int) -> int:
    """base ** exponent % modulus, reducing after every multiplication."""
    if modulus == 1:
        return 0                      # everything is congruent to 0 mod 1
    result = 1
    current = base % modulus
    remaining = exponent

    while remaining > 0:
        if remaining % 2 == 1:
            result = (result * current) % modulus
        remaining //= 2
        if remaining > 0:
            current = (current * current) % modulus

    return result


print(fast_power_mod(7, 1000, 13), pow(7, 1000, 13))
print(len(str(7 ** 1000)))
print(fast_power_mod(2, 10 ** 18, 1_000_000_007), pow(2, 10 ** 18, 1_000_000_007))

9 9
846
719476260 719476260

Three details earn their place. base % modulus up front handles a base larger than the modulus. The modulus == 1 guard exists because every integer is congruent to 0 modulo 1, so result = 1 would otherwise be returned unchanged. And the reduction happens on both the accumulator and the squaring chain — miss either and the numbers grow without limit.

Use the built-in

Python has this in the language. pow(base, exponent) is the same as base ** exponent, and the three-argument pow(base, exponent, modulus) is modular exponentiation implemented in C, using a windowed variant of the same square-and-multiply algorithm. Since Python 3.8, a negative exponent with a modulus gives the modular inverse: pow(a, -1, m) is the value x with a × x ≡ 1 (mod m), and it raises ValueError when no inverse exists, which happens exactly when a and m share a factor.

print(pow(3, 13))
print(pow(3, 13, 100))

private_exponent = pow(17, -1, 3120)
print(private_exponent, (17 * private_exponent) % 3120)

try:
    pow(6, -1, 9)
except ValueError as error:
    print(f"ValueError: {error}")

print(pow(7, 11, 13), pow(7, -1, 13))

1594323
23
2753 1
ValueError: base is not invertible for the given modulus
2 2

That 2753 is not a random number: it is the RSA private exponent from the textbook key with primes 61 and 53, where the modulus is 3233 and e is 17. The last line shows the other route to an inverse — by Fermat's little theorem, a^(p-2) ≡ a^(-1) (mod p) for prime p, so pow(7, 11, 13) and pow(7, -1, 13) agree.

Complexity

Multiplications: O(log n). The loop runs once per bit of n, and the number of bits is floor(log2(n)) + 1. Each iteration does at most two multiplications: one squaring, plus one fold into the result when the bit is set. So the total is at most 2 × log2(n) + 1 and at least log2(n). The exact count is (bits - 1) + popcount(n), and the table above shows it: 10^18 has 60 bits and 24 of them set, giving 59 + 24 = 83.

The honest caveat: a multiplication is not O(1). Big O on multiplication counts is the right measure only when each multiplication has fixed cost. For pow(a, n, m) that holds — every operand stays below m, so each multiplication costs the same, and the total is O(log n) machine-word multiplications for fixed m.

Without a modulus it does not hold at all. a^n has about n × log2(a) bits, so the final squaring alone operates on numbers half that size, and CPython multiplies large integers with Karatsuba's algorithm at roughly O(d^1.585) for d-digit operands. The cost is then dominated by the last one or two multiplications, and it is proportional to the size of the answer — which no algorithm can beat, because you have to write the answer down.

Space: O(1) for the iterative version — three integers, regardless of n. The recursive version adds O(log n) stack frames. With a modulus, every one of those integers stays below m, so the space is genuinely constant in n.

Matrix exponentiation, and Fibonacci in O(log n)

Nothing in the algorithm cares that the values are numbers. It needs multiplication, associativity and an identity element — nothing else. Matrices have all three.

Linear recurrences can be written as matrix powers. The Fibonacci step F(n+1) = F(n) + F(n-1) is exactly what the matrix [[1, 1], [1, 0]] does to the vector (F(n), F(n-1)). Raise that matrix to the nth power and its top-right entry is F(n).

Matrix = list[list[int]]


def matrix_multiply(left: Matrix, right: Matrix) -> Matrix:
    """Multiply two 2x2 integer matrices."""
    return [
        [left[0][0] * right[0][0] + left[0][1] * right[1][0],
         left[0][0] * right[0][1] + left[0][1] * right[1][1]],
        [left[1][0] * right[0][0] + left[1][1] * right[1][0],
         left[1][0] * right[0][1] + left[1][1] * right[1][1]],
    ]


def matrix_power(matrix: Matrix, exponent: int) -> Matrix:
    """The same square-and-multiply loop, with matrices instead of numbers."""
    result = [[1, 0], [0, 1]]        # the identity matrix plays the role of 1
    current = matrix
    remaining = exponent

    while remaining > 0:
        if remaining % 2 == 1:
            result = matrix_multiply(result, current)
        remaining //= 2
        if remaining > 0:
            current = matrix_multiply(current, current)

    return result


def fibonacci(index: int) -> int:
    """F(index), using O(log index) matrix multiplications."""
    return matrix_power([[1, 1], [1, 0]], index)[0][1]


print([fibonacci(n) for n in range(11)])
print(fibonacci(90))
print(len(str(fibonacci(1_000_000))))

[0, 1, 1, 2, 3, 5, 8, 13, 21, 34, 55]
2880067194370816120
208988

The loop body is identical to fast_power_iterative with * replaced by matrix_multiply and 1 replaced by the identity matrix. The last line computes the millionth Fibonacci number — 208,988 digits long — with 39 matrix multiplications. The usual iterative Fibonacci would need a million additions on numbers of that size.

The same construction handles any recurrence of the form T(n) = c1×T(n-1) + ... + ck×T(n-k): build the k-by-k companion matrix and raise it to the nth power in O(k^3 log n) operations.

When to use it, and when not to

Use it whenever the exponent is large and the modulus is present. That is the entire domain: cryptography, hashing, counting problems modulo a prime like 10^9 + 7, and linear recurrences via matrices. Under a modulus there is no reason ever to write the naive loop.

Do not write it yourself in Python. pow(base, exponent, modulus) does the same thing in C, and beats a Python-level loop by a wide margin. The value in the code above is knowing what pow costs and why, not shipping it.

Do not use it for cryptographic code you deploy. The textbook loop branches on the bits of the exponent, so its running time and power draw leak the secret exponent — this is the basis of Paul Kocher's 1996 timing attacks on Diffie-Hellman and RSA. Real implementations use constant-time variants such as the Montgomery ladder; OpenSSL's is BN_mod_exp_mont_consttime. Use a vetted library.

Do not bother for small fixed exponents. x * x * x beats any general routine for a cube.

Do not use it to build an exact huge power you then throw away. If you only need the last few digits or the size of the answer, use a modulus or math.log instead. 2 ** 10**9 is a valid Python expression and a bad idea.

Where it shows up in the real world

RSA. Both encryption and decryption are single modular exponentiations: the ciphertext is m^e mod n and the plaintext is c^d mod n. The common public exponent 65537 is chosen precisely because it is 2^16 + 1 — binary 10000000000000001, two bits set — so encryption costs 16 squarings and one multiply. The private exponent is a full-length 2048-bit number, which is why signing is far slower than verifying.

Diffie-Hellman and its elliptic-curve variant. Every TLS handshake computes g^a mod p with a secret exponent, and elliptic-curve Diffie-Hellman runs the identical algorithm with point addition in place of multiplication, where it is called double-and-add.

Primality testing. The Miller-Rabin and Fermat tests are built on computing a^(n-1) mod n for candidate n. Every large prime used in cryptography was found by running modular exponentiation a few hundred times.

Rolling hashes. Rabin-Karp and Rabin fingerprints need base^(window - 1) mod m to drop the outgoing character. With a window of 10,000 that is 27 multiplications instead of 10,000.

Modular inverses. RSA key generation derives the private exponent d as the inverse of e modulo (p-1)(q-1), which is one pow(e, -1, phi) call.

Linear recurrences. Matrix exponentiation is the standard way to answer "the nth term of this recurrence, modulo p, for n up to 10^18".

Common mistakes

Computing the power first, then the modulus. (a ** b) % m is the single most common bug here. It is mathematically correct and computationally hopeless — Python builds the entire a^b before reducing. Always pow(a, b, m).

Forgetting to reduce the squaring chain. Reducing result but not current keeps the answer correct and lets current grow to the full unreduced size, which throws away the whole point.

Squaring before multiplying into the result. Inside the loop, result *= current must come before current *= current. Reversed, every set bit contributes one power too many.

Returning 1 for a modulus of 1. pow(x, y, 1) is 0, because every integer is congruent to 0 modulo 1. A hand-written version that starts result = 1 and never multiplies returns 1 instead. Guard it.

Negative exponents. fast_power(2, -3) never reaches exponent == 0 with the halving rules as written, because -3 is odd and -4 // 2 is -2, then -1, then -2 again. Either raise for negative exponents or convert them to 1 / fast_power(base, -exponent).

Using math.pow. It returns a float and loses exactness above 2^53. math.pow(3, 40) is not an integer answer. Use ** or pow for integers.

Practice

1. Rewrite fast_power so it raises a ValueError on a negative exponent, and confirm the error fires.
2. Write the top-down version that walks bin(exponent) from the leading bit, squaring the result each step and multiplying by the base when the bit is 1, then check it agrees with fast_power_iterative.
3. Compute the last ten digits of 2^(10^18) using a modulus of 10^10, and explain why the answer is even.
4. Implement a modular inverse using Fermat's little theorem for a prime modulus, and compare it against pow(a, -1, m) for every a from 1 to 100 with modulus 101.
5. Use matrix_power to compute the nth term of T(n) = 2×T(n-1) + 3×T(n-2) with T(0) = 0 and T(1) = 1, modulo 10^9 + 7, for n equal to one million.

Summary

Fast exponentiation is one identity applied repeatedly: halve the exponent and square, peel a factor off when the exponent is odd. That turns n multiplications into (bits - 1) + popcount(n), which is never more than 2 × log2(n) + 1. Reduce modulo m after every multiplication and the numbers stay small too, which is what makes public-key cryptography possible at all. In Python you write pow(base, exponent, modulus) — but knowing what those three arguments cost is the difference between code that returns in a microsecond and code that never returns.

Keep reading

• The Euclidean Algorithm — the other halving-style number-theory algorithm, and the extended version that also computes modular inverses.
• Sieve of Eratosthenes — where the primes used as moduli come from.
• Rabin-Karp String Matching — rolling hashes, which need a modular power of the base.
• Big O Notation — why log n and n are separated by more than most people expect.