Fast exponentiation
In short: Fast exponentiation computes aⁿ with O(log n) multiplications by squaring and reading the binary digits of n; taking the remainder after each step gives modular exponentiation.
Problem
The key of an HTTPS connection
When the browser opens an HTTPS page, it and the server agree on a shared secret key without sending it over the network. In the classic version (Diffie–Hellman) each side computes a modular power: gᵃ mod p, where a is a secret number tens of digits long. The same computation sits in RSA signatures and in the hash of a text, which multiplies by powers of a base.
The naive way is a loop that multiplies n times:
def power_loop(a, n, m): result = 1 for _ in range(n): result = result * a % m return resultWith a 77-digit exponent (for example 2²⁵⁵ − 19) that is about 6 · 10⁷⁶ multiplications. At a billion multiplications per second the loop would run for 10⁶⁷ seconds, and the universe is about 4 · 10¹⁷ seconds old. Yet the server makes thousands of such connections per second.
The idea
The square of the half
3¹³ does not need 13 multiplications. 3¹³ = 3⁶ · 3⁶ · 3, and 3⁶ = 3³ · 3³. Every step halves the exponent:
- aⁿ = (aⁿᐟ²)², if n is even;
- aⁿ = (aⁿᐟ²)² · a, if n is odd (n / 2 is integer division);
- a⁰ = 1: this is the base case.
def power(a, n, m): if n == 0: return 1 half = power(a, n // 2, m) square = half * half % m if n % 2 == 1: return square * a % m return square
print(power(3, 13, 1_000_000_007))print(power(5, 1_000_000_000, 1_000_000_007))1594323142848001The remainder mod m after every multiplication keeps the numbers small: (x · y) mod m is the same as ((x mod m) · (y mod m)) mod m. Here m = 1,000,000,007, a prime and a common bound in hashes and programming contests.
What does it cost? Every call divides n by 2, just like binary search, so the depth is about log₂ n: for n = 10⁹ that is 30 frames. The call is not a tail call (squaring comes after it), but from the first lesson you know that 30 frames are nothing for the stack. The work is at most two multiplications per level: O(log n) instead of O(n).
The idea
Without recursion: the binary digits of n
A loop can do the same. 13 in binary is 1101: 13 = 8 + 4 + 1, so 3¹³ = 3⁸ · 3⁴ · 3¹. The loop reads the digits of n from right to left (n % 2 is the last one, and n /= 2 removes it) and keeps a squared: 3¹, 3², 3⁴, 3⁸… On a digit 1 the current a goes into the result.
def power(a, n, m): result = 1 a %= m while n > 0: if n % 2 == 1: result = result * a % m a = a * a % m n //= 2 return result
print(power(3, 13, 1_000_000_007))print(power(5, 1_000_000_000, 1_000_000_007))1594323142848001The multiplications are as many as the digits of n (one squaring per digit) plus the ones in it (one multiplication of the result each). This is the version libraries use: no recursion, and no dependence on the size of the stack.
See it
The squares go into the result
def power(a, n, m): result = 1 a %= m while n > 0: if n % 2 == 1: result = result * a % m a = a * a % m n //= 2 return result
print(power(3, 13, 1_000_000_007)) # 1594323Step by step, in words
- result starts at 1. n = 13, which is 1101 in binary: each binary digit decides whether a goes into the result.
- n = 13 is odd: result = 1 · 3 = 3. Then a becomes 3² = 9 and n becomes 6.
- n = 6 is even: result stays 3. Only a is squared: 9 · 9 = 81 = 3⁴, and n = 3.
- n = 3: result = 3 · 81 = 243 = 3⁵. a = 6561 = 3⁸, n = 1.
- n = 1: result = 243 · 6561 = 1594323 = 3¹³, because 13 = 8 + 4 + 1. The loop still squares a once more: 3¹⁶ is a multiplication nobody uses.
- n = 0 and the loop stops: 7 multiplications instead of 13. For n = 10⁹ it is 43 instead of a billion.
See it
Counting the multiplications
Both versions with a counter, for exponents from 10 to 10,000,000, and the fast one for a billion too:
multiplications = 0
def power_loop(a, n, m): global multiplications result = 1 for _ in range(n): multiplications += 1 result = result * a % m return result
def power(a, n, m): global multiplications result = 1 a %= m while n > 0: if n % 2 == 1: multiplications += 1 result = result * a % m multiplications += 1 a = a * a % m n //= 2 return result
m = 1_000_000_007for n in [10, 1_000, 100_000, 10_000_000]: multiplications = 0 power_loop(5, n, m) slow = multiplications multiplications = 0 power(5, n, m) print(f"n = {n}: loop {slow}, fast {multiplications} multiplications")multiplications = 0power(5, 1_000_000_000, m)print(f"n = 1000000000: fast {multiplications} multiplications")n = 10: loop 10, fast 6 multiplicationsn = 1000: loop 1000, fast 16 multiplicationsn = 100000: loop 100000, fast 23 multiplicationsn = 10000000: loop 10000000, fast 32 multiplicationsn = 1000000000: fast 43 multiplicationsRun them here as well. Bench counts by the same rule (one multiplication mod m is one operation), so the numbers are the ones in the output above. We do not run a billion with the naive loop: Bench estimates it from the curve.
- Loop: n multiplications
O(n) - Squaring
O(log n)
| Algorithm | n = 10 | n = 1,000 | n = 100,000 | n = 10,000,000 | n = 1,000,000,000 |
|---|---|---|---|---|---|
| Loop: n multiplications | 10 | 1,000 | 100,000 | 10,000,000 | ≈ 1,000,000,000 |
| Squaring | 6 | 16 | 23 | 32 | 43 |
def power(a, n, m): result = 1 a %= m while n > 0: if n % 2 == 1: result = result * a % m a = a * a % m n //= 2 return resultThe naive loop grows like n, the fast one like the number of digits of n: an n ten times bigger adds about 3.3 binary digits, that is, 3 to 7 multiplications.
Try it
15 or 16
Which exponent needs more multiplications?
multiplications = 0
def power(a, n, m): global multiplications result = 1 a %= m while n > 0: if n % 2 == 1: multiplications += 1 result = result * a % m multiplications += 1 a = a * a % m n //= 2 return result
for n in [15, 16]: multiplications = 0 power(3, n, 1_000_000_007) print(f"{n}: {multiplications} multiplications")15 = 1111: four digits and four ones, 8 multiplications. 16 = 10000: five digits and one one, 6 multiplications. A bigger exponent does not mean more work: what counts are the digits and the ones.
The idea
When the product does not fit
The remainder keeps the numbers below m, but the product of two such numbers is up to m². If m² does not fit in the type, the computation is wrong, and silently so. Here is the same power with a bigger modulus:
m = 10_000_000_019print("power:", power(5, 1_000_000_000, m))print("pow: ", pow(5, 1_000_000_000, m))power: 475867028pow: 475867028int grows as needed, so the same code gives the correct answer with a big modulus too. The price is that every multiplication of a big number is slower. The built-in pow(a, n, m) does the same squaring, but in C.The idea
The built-in, and a 77-digit exponent
In real code a modular power is not written by hand:
secret = 2**255 - 19p = 2**127 - 1print("Digits in the exponent: " + str(len(str(secret))))print("Multiplications with squaring: " + str(secret.bit_length() + secret.bit_count()))print("5^secret mod p =", pow(5, secret, p))Digits in the exponent: 77Multiplications with squaring: 5085^secret mod p = 97285997395055909576908201857106207893pow(a, n, m) with three arguments does squaring and takes the remainder along the way. a ** n % m is correct, but first builds the whole aⁿ: with a 77-digit exponent it would never finish. 508 multiplications of 39-digit numbers instead of 6 · 10⁷⁶: that is why the key is agreed in milliseconds. Libraries need even fewer, because they look at several digits at a time, but the idea is the same.
Try it
Find the bug
This version returns the correct result. But for n = 1,000,000 it does 2,048,575 multiplications instead of 27. Which line is to blame?
Pick the line with the bug.
The half is computed twice. Every level doubles the calls: 2, 4, 8… and at the bottom there are about n. The depth is still log₂ n, so the stack notices nothing, but the work is back to O(n). The fix is to put the half in a variable and then half * half.
Try it
Fill in
When does the current a go into the result?
def power(a, n, m):
result = 1
a %= m
while n > 0:
if :
result = result * a % m
a = a * a % m
n //= 2
return resultWhen the last binary digit of n is 1, that is, when n is odd.
Exercise
Interest to the cent
Solve this in your editor.
What the deposit grows to. You put €1000 at 1% interest per month, which is added to the amount. Read the number of months and print the amount rounded to the cent. Compute it exactly, with integers: €1000 is 100,000 cents, and after n months the amount is 100,000 · 101ⁿ / 100ⁿ cents. Write your own fast exponentiation with int, without a modulus, and round half a cent up.
120After 120 months: €3300.39For 120 months the multiplications are few even with the naive loop. Here, though, the numbers grow (101¹²⁰ has 241 digits), so every next multiplication costs more, and saving multiplications matters even more.
More tasks
- Daily interest for 30 years. The same with 0.01% per day for 10,950 days. Count the multiplications of both versions.
- The last digits. Print the last 6 digits of 2¹⁰⁰⁰⁰⁰⁰⁰⁰⁰ (modulus 1,000,000), without building the whole number.
More precisely
More about powers
More precisely, if you are curious
Fibonacci numbers in O(log n). The same squaring works with a 2 × 2 matrix instead of a number: the n-th Fibonacci number comes out of the n-th power of [[1, 1], [1, 0]], with about 2 log₂ n matrix multiplications.
Elliptic curves. Today TLS most often agrees on the key with an elliptic curve (X25519). There, instead of multiplication there is "adding points", and instead of squaring there is "doubling", but the algorithm is the same: double-and-add over the binary digits of the secret number.
Constant time. The loop in the lesson multiplies the result only on a digit 1. Someone measuring the time can guess how many ones the secret number has. That is why cryptographic libraries do the same work for every digit (the Montgomery ladder), and why cryptography is not written by hand, but taken from a library.
- Python:
pow(a, -1, m)(since 3.8) gives the inverse of a mod m.a ** nwith integers is exact, but grows without limit.
Test