Fast Exponentiation
In short: write the exponent in binary and square-and-multiply computes base**exp % mod in O(log exp) time and O(1) extra space — one squaring of the base per bit, and one
extra multiply into the accumulator wherever that bit is 1. That turns a cost that grows
with the exponent's value into one that grows with its digit count, which is the
difference between instant and 30,000 years for an exponent like 10¹⁸. The same loop
works for anything associative you can multiply: numbers, matrices, permutations.
Later in the course Quarry gets a text index, and its core is a rolling hash. The
hash is a number computed from a window of characters, updated in O(1) as the window
slides one character to the right. To drop the character leaving on the left, the hash
needs B**(k-1) % p, where B is the hash base, k the window length and p a large
prime. Checksums, id scramblers and the modular inverses in the next lesson all need the
same thing: a big power, reduced modulo something, computed quickly.
Python has it built in, as pow(base, exp, mod). This lesson is about how it works,
because the idea behind it isn't really about numbers. It's about reading the
exponent's bits, and it applies to anything you can multiply: matrices, permutations,
functions. The interview corner computes the 10¹⁸th Fibonacci number with it.
Multiplying e times is too slow
The obvious loop multiplies e times:
For e = 45 that's fine. For e = 10**18, a perfectly ordinary exponent in hashing
and cryptography, it would take 30,000 years at a billion multiplications per second.
The loop's cost grows with the value of the exponent. We want it to grow with the
number of digits in the exponent.
Read the exponent in binary
Write the exponent in binary. 45 is 101101, which is 32 + 8 + 4 + 1, so
b**45 = b**32 · b**8 · b**4 · b**1
The powers b, b**2, b**4, b**8, b**16, b**32 each come from the one before by a single
squaring. So walk the exponent's bits from the lowest. Keep squaring the base, and
whenever the current bit is 1, multiply the base into an accumulator. That's the
square-and-multiply algorithm, also called binary exponentiation. Step through
7**45 % 1000:
Each bit of the exponent either multiplies the current base into the answer or skips it, and the base squares every step.
Bit 0 is 1, so acc = 1 × 7 = 7. Then square the base: 7² = 49. 6 cells. At index 5.
In code, e & 1 reads the lowest bit and e >>= 1 moves to the next, both from the
first lesson:
The loop runs once per bit of e, so it takes e.bit_length() passes, about 60 for
10**18:
e, steps = 1000, 0
while e:
e >>= 1
steps += 1
print(steps)Why taking % m at every step is allowed
(x · y) % m == ((x % m) · (y % m)) % m. Remainders multiply the same way the numbers
do, so you can reduce after every multiplication and never hold anything bigger than
m². Without that, 7**(10**18) has about 8 × 10¹⁷ digits, more memory than any
machine has, even in Python, whose ints never overflow. In C or Java the product would
overflow the 64-bit word first, which is lesson 1's problem again. There you keep
m below 2³², so m² fits in 64 bits, or use a 128-bit intermediate.
print(pow(3, 0, 7), pow(5, 3, 1))Here is the same loop with the order of two lines swapped. Squaring before using the base means every bit uses a power one step too high:
It prints 49 576, which are 3¹⁰ and 2²⁰ mod 1000. Bit 0 is multiplied by b² instead of b, so every power comes out doubled.
Trace a small case and watch acc, b and e move together:
The loop never stores the powers it squares, it just overwrites b. Keep them in a
list and you get a table of b, b², b⁴, … mod m, and the answer multiplies
together exactly the entries whose exponent bit is set. Run your own version below. The
done cells are the ones used. Try an exponent with many 1 bits (2**10 - 1), then
one with a single 1 bit (2**10):
Your table of repeated squares mod m. The answer is the product of exactly the cells whose exponent bit is 1.
How few multiplications does this take for a big exponent? For an exponent with x
decimal digits it needs at most about 2 · x · log₂10 of them, one squaring and at most
one multiply per bit. Quarry's hash kernel can afford 120 modular multiplications per
call:
| digits in the exponent | worst-case multiplications |
|---|---|
| 5 | 33.2 |
| 10 | 66.4 |
| 18 | 119.6 |
| 19 | 126.2 |
So 18 digits fit: exponents up to about 10¹⁸, which needs about 60 bits. The naive loop would need 10¹⁸ multiplications. That jump from linear in the value to linear in the digit count is the whole point of reading the exponent in binary.
It works for anything you can multiply
The loop only needs a multiplication that is associative: (x·y)·z = x·(y·z). So it
works just as well for matrices, permutations and function composition. Two things
change: the "1" you start from becomes an identity (the identity matrix, or the
permutation that moves nothing), and * becomes that structure's multiplication.
Interview corner — the 10¹⁸th Fibonacci number
Return F(n) mod 1_000_000_007 for n up to 10¹⁸. A loop is O(n), which is out of the
question. The trick is that one step of the recurrence is a matrix multiplication:
[F(k+1), F(k)] = [[1, 1], [1, 0]] × [F(k), F(k-1)]. So n steps are one matrix raised
to the nth power, and square-and-multiply does that in about 60 matrix products:
It's the same loop as mod_pow, with a different multiplication. Whenever something
must be applied a huge number of times and the applications compose, reach for the
exponent's bits.
Quick practice
def mod_pow(b, e, m):
acc = 1 % m
b %= m
while e:
if e & 1:
acc = acc * b % m
b = b * b % m
e >>= 1
return acc
print(mod_pow(2, 10, 1000))It prints 1 1. Everything mod 1 is 0, so mod_pow(5, 0, 1) should be 0, but with e = 0 the loop never runs to correct it.
🎯 Your turn
Write mod_pow(base, exp, mod) with square-and-multiply, without calling pow or
using **. It should give base ** exp % mod for exp >= 0 and mod >= 1, using
O(log exp) multiplications.
mod_pow(7, 45, 1000)→807mod_pow(5, 0, 1)→0
Hint — keep acc, starting at 1 % mod. While exp is non-zero: if exp & 1,
multiply base into acc, then square base and shift exp right. Reduce mod mod
after every multiplication.
Then press ▶ Run, try the Live App chips, and hit ✓ Check. ✅