Fast power in Python

Why this matters

Raising a number to a power is a basic step in many programs. Compound interest, hashing, checking whether a number is prime, and counting ways to arrange things all use it. A naive program multiplies the base by itself once for every unit of the exponent. That works for small exponents, and it becomes very slow for large ones.

Take the power 2 to the power 1,000,000. The naive method does 1,000,000 multiplications. The method in this lesson does about 27. The idea is simple once you see it: squaring a number doubles its exponent in one step. So you can reach a large exponent in a few steps, provided you keep track of the odd steps.

The recursion pattern from the first lesson of this module is used here. The base case is n equal to 0, and every other call makes the problem smaller by halving n. Each call does a small amount of work, so the number of calls is the number of halvings, which is about the number of binary digits of n.

The idea: squaring doubles the exponent

Start with a small case. You want 3 to the power 8. The naive method multiplies eight 3s together, which takes 7 multiplications between them.

Look at the exponent 8 instead. It is 2 times 4, and 4 is 2 times 2. So 3 to the power 8 is the same as the square of 3 to the power 4, and 3 to the power 4 is the square of 3 to the power 2. Work it from the bottom:

1. 3 to the power 2 is 3 times 3, which is 9. One multiplication.
2. 3 to the power 4 is 9 times 9, which is 81. One multiplication.
3. 3 to the power 8 is 81 times 81, which is 6,561. One multiplication.

That is 3 multiplications instead of 7. The trick is that each square is a power you need anyway.

The second case is an odd exponent, because squaring alone cannot reach it. For 3 to the power 5, first find 3 to the power 2, which is 9. Then square it to get 3 to the power 4, which is 81. The exponent 5 is 4 plus 1, so multiply 81 by one extra 3. The result is 243.

So there are two rules:

• If n is even, x to the power n is half times half, where half is x to the power n // 2.
• If n is odd, x to the power n is half times half times x.

Each recursive call halves n, so the number of calls is about the number of times you can divide n by 2 before it reaches 0. That is the base-2 logarithm of n, written log n. For n equal to 1,000,000 that is about 20 halvings.

Example 1: the recursive function

The function below follows the two rules. It also counts multiplications, so you can compare it with the naive loop. The naive loop is written first, so the comparison is clear. The counter mults is a global variable, and each function changes it with global mults.

The base case returns 1 when n is 0, because any number to the power 0 is 1. The expression n // 2 uses floor division, which drops the remainder. So 13 divided by 2 is 6, and the remainder 1 tells the function that 13 is odd. The check n % 2 == 0 looks at the remainder to decide which rule to use.

Notice that fast is called once, not twice. The result of the call is stored in half, and then it is used twice. That is the reason the function does roughly log n work instead of n.

The counter adds one for each squaring. It adds one more when the exponent is odd, for the extra multiplication by x. Press Run. The output is:

2^10 naive: 1024
2^10 fast: 1024
3^1 fast=3 naive=3 same=True mults fast=2 naive=1
3^2 fast=9 naive=9 same=True mults fast=3 naive=2
3^3 fast=27 naive=27 same=True mults fast=4 naive=3
3^5 fast=243 naive=243 same=True mults fast=5 naive=5
3^8 fast=6561 naive=6561 same=True mults fast=5 naive=8
3^13 fast=1594323 naive=1594323 same=True mults fast=7 naive=13
3^100 fast=515377520732011331036461129765621272702107522001 naive=515377520732011331036461129765621272702107522001 same=True mults fast=10 naive=100
n=1000 fast mults: 16
n=1000000 fast mults: 27

The small cases show the pattern. For 3^1 the fast function does 2 multiplications and the naive loop does 1. The fast function is slightly worse there, because on its last level it squares the value 1 without needing to, and then multiplies by 3. Its advantage starts to show from about 3^8 onward, where it does 5 multiplications and the loop does 8. At 3^13 it does 7 and the loop does 13.

For n equal to 1,000,000, the fast function does 27 multiplications. The naive loop would do 1,000,000. The count 27 is 20 halvings plus 7 extra multiplications for the odd steps. The number 1,000,000 has 7 ones in its binary form, and each one causes an extra multiplication. This is the reason the work is O(log n): the number of calls grows with the number of binary digits, not with the size of n.

The line for 3 to the power 100 shows an exact answer with 48 digits. Python integers have no size limit, so the power is never wrong from overflow. The cost of a very large power is the memory and time for the big numbers, which the Common mistakes section returns to.

Example 2: the iterative method and the bits of n

The recursive function works from the top: it halves n and waits for the answer to come back. There is also a loop version that works from the bottom. It reads the binary digits of n one at a time, from the lowest digit up.

Binary is the base-2 way of writing a number. Each digit is either 0 or 1, and each position stands for a power of 2. The number 13 is 1101 in binary, because 8 plus 4 plus 1 is 13. The positions stand for 1, 2, 4, and 8. So 13 is made from the positions with a 1 in them: 1, 4 and 8. Python writes a binary number with the prefix 0b, so bin(13) gives 0b1101.

That gives a second way to see fast power. Each binary digit of n that is 1 says that one power of x must be multiplied into the result. Those powers are x to the power 1, x to the power 2, x to the power 4, x to the power 8, and so on. The loop makes each of these powers by squaring the one before it.

The output is:

iterative(3, 13) = 1594323
matches 2 ** 1000: True

The expression e & 1 gives the lowest binary digit of e. The & operator compares the binary digits of two numbers, and with 1 it keeps only the lowest digit. If that digit is 1, the current power of base is multiplied into result. Then base is squared, so it moves to the next power. Then e >> 1 shifts all the binary digits one place to the right, which drops the lowest digit. The loop ends when e is 0.

Here is a trace for 3 to the power 13. The exponent 13 is 1101 in binary.

step  e (before)  lowest bit  result (after)  base (after squaring)
1     13          1           3               9
2     6           0           3               81
3     3           1           243             6561
4     1           1           1594323         43046721

At step 1, the lowest bit is 1, so result becomes 3. The base squares to 9, which is 3 to the power 2. At step 2, the lowest bit is 0, so result does not change. The base squares again to 81, which is 3 to the power 4. At step 3, the bit is 1, so result becomes 3 times 81, which is 243. That is 3 to the power 5, from the powers 1 and 4. The base squares to 6,561, which is 3 to the power 8. At step 4, the bit is 1 again, so result becomes 243 times 6,561, which is 1,594,323. That is 3 to the power 13.

The loop runs once for each binary digit of n. The number 13 has 4 binary digits, so it runs 4 times. The number 1,000 has 10 binary digits, so it runs 10 times. The int.bit_length() method returns the number of binary digits directly, so (1000).bit_length() is 10.

Example 3: large results and the modulo

Many programs do not want the exact value of a huge power. They want the answer with the remainder after dividing by a fixed number. Competitive programming problems often ask for the answer modulo 1,000,000,007 (the remainder after dividing by that number), which is a prime number. Using the remainder keeps every number small, so the program runs faster and uses less memory.

The remainder after dividing by a number is written % in Python. The function below takes the remainder after each multiplication. This keeps every intermediate value below the modulus. Python also has a built-in version that takes the modulus as a third argument, pow(x, n, MOD). The program checks the function against that built-in.

The output is:

2^10 mod p = 1024, pow = 1024, same = True
7^1000 mod p = 224787023, pow = 224787023, same = True
123456789^987654321 mod p = 652541198, pow = 652541198, same = True
10^1000000000 mod p = 142857001, pow = 142857001, same = True

Every line says same = True, so the function and the built-in pow give the same remainder in every case.

The reason every step reduces the value is size. Each multiplication of two large numbers makes a larger number, and a number with a million digits takes much longer to multiply than one with ten digits. Reducing after each step keeps the numbers small. A product of two numbers below 1,000,000,007 is about 10 to the power 18, which is small work for Python.

The expression x % MOD is needed because x may be bigger than the modulus. The base case returns 1 before any multiplication happens.

Python's % has one rule that differs from many other languages. The result of % takes the sign of the divisor, not the sign of the number. So (-3) % 5 is 2, which is the remainder a mathematician expects. Some other languages give -3 for the same expression. This matters when a negative base is reduced by a modulus.

Common mistakes

Calling the half-power twice

A tempting version calls the function twice instead of saving the result:

The output is:

n=8 two-calls mults=23 value=6561
n=16 two-calls mults=47 value=43046721
n=20 two-calls mults=51 value=3486784401

It gives the right answer, but it is slow. Each call makes two calls of half the size, so the total work doubles at each level. That is about n work, not log n. The version that saves half in a variable does 5, 6 and 7 multiplications for the same three inputs. For n equal to 16, the version with two calls does 47 multiplications, and the saving version does 6. The gap grows with n. The fix is to store the result of the call in a variable, as the main function does with half.

Forgetting the odd step

If the function squares half for an odd exponent and leaves out the extra multiplication by x, the answer is wrong. This version has the bug:

The output is:

2^5 wrong: 1
2^4 wrong: 1

The correct answers are 32 and 16. The bug shows up at the level where n is 1. There half is 1, and squaring 1 stays 1, so every later result is 1 as well. Test the function on a small odd exponent such as 1 or 3 to find this kind of bug early.

Missing base case

The base case must stop the recursion. If n equal to 0 is not handled, the function keeps calling itself with the same value. This version has no base case:

def no_base(x: int, n: int) -> int:
    half = no_base(x, n // 2)
    return half * half if n % 2 == 0 else half * half * x
 
 
print(no_base(2, 10))

The program stops with an error. The first lines of its output are these, and the middle lines repeat until the limit is reached:

Traceback (most recent call last):
  File "no_base.py", line 6, in <module>
    print(no_base(2, 10))
          ~~~~~~~^^^^^^^
  File "no_base.py", line 2, in no_base
...
  [Previous line repeated 996 more times]
RecursionError: maximum recursion depth exceeded

The value of n goes 10, 5, 2, 1, 0, and then stays at 0 forever, because 0 // 2 is still 0. Each call adds a frame, and Python stops the chain of calls at its limit of 1000 frames. The fix is the base case if n == 0: return 1 at the top of the function.

A negative exponent

The functions in this lesson handle n of 0 or more. For a negative n, the answer is 1 divided by x to the power of -n. The base case is never reached, because halving a negative number moves it towards -1, not towards 0. So a negative exponent fails the same way as a missing base case:

def fast(x: int, n: int) -> int:
    if n == 0:
        return 1
    half = fast(x, n // 2)
    return half * half if n % 2 == 0 else half * half * x
 
 
print(fast(2, -1))

The output ends with the same RecursionError as the missing base case. For negative exponents, use x ** n, which returns a float (2 ** -1 is 0.5), or check the sign at the top of the function and raise an error.

Overflow is not a problem in Python, but size is

Python does not wrap a large result around, so 3 ** 100 is always correct. A float does have a limit. A float power that is too big raises an error instead of wrapping, and a float product that is too big becomes inf. The expression 2.0 ** 1100 gives this message:

OverflowError: (34, 'Numerical result out of range')

Use integers for exact powers. Use the modulus version when the remainder is enough, because it keeps the numbers small.

Practice

Exercise 1

How many times does the while loop in iterative run for n equal to 1000?

Show the answer

It runs 10 times. The number 1000 is 1111101000 in binary, and that has 10 digits. The loop runs once for each digit. (1000).bit_length() returns 10 as well.

Exercise 2

Use the recursive function fast to find 5 to the power 6. Write down the return value of each call that matters.

Show the answer

fast(5, 1) returns 5, because half is fast(5, 0), which is 1, its square is 1, and 1 times 5 is 5. fast(5, 3) returns 125, because half is 5, its square is 25, and 25 times 5 is 125. fast(5, 6) returns 15,625, because half is 125 and its square is 15,625. So 5 to the power 6 is 15,625.

Exercise 3

The exponent 37 is 100101 in binary. Which powers of x are multiplied together to get x to the power 37?

Show the answer

The powers are x to the power 1, x to the power 4 and x to the power 32. Those are the positions with a 1 in 100101, counted from the right: positions 0, 2 and 5. Their exponents add up to 1 + 4 + 32, which is 37.

Exercise 4

Change MOD in the mod_pow program from 1,000,000,007 to 1000, and compute 2 to the power 10 with mod_pow. What is the output? Does the function need MOD to be a prime number?

Show the answer

The output is 24, because 1024 divided by 1000 leaves the remainder 24. The function does not need MOD to be prime. It only multiplies numbers and takes remainders, and those steps work for any positive modulus.

Exercise 5

Explain in your own words why the version that calls the half-power twice does about n multiplications, while the version that stores half does about log n.

Show the answer

The version with two calls makes two calls of size n / 2. Each of those makes two calls of size n / 4, and so on. The total number of calls doubles at each level, and the number of levels is about log n, so the total is about n. The version that stores half makes one call per level, so the total number of calls is about log n.

Summary

• Squaring the base doubles the exponent in one multiplication, so the number of steps grows with the number of binary digits of n, not with n.
• A recursive function with the two rules, even and odd, does about log n multiplications, plus one more for each odd step.
• The loop over the binary digits of n gives the same result from the lowest digit upward.
• Take the remainder after each multiplication when only the answer modulo a number is needed, and check the result against the built-in pow(x, n, MOD).
• Store the result of a recursive call in a variable rather than calling the function twice, and make sure the base case stops the recursion.